Assign_Read_Vis_Biosignals_Images

School

Simon Fraser University *

*We aren’t endorsed by this school

Course

340

Subject

Biology

Date

Apr 3, 2024

Type

Pages

Uploaded by DeaconLion805

%% CMPT340 Assignment - Reading & Visualizing Bio-signals % ============================ Goal =========================: % In this activity, you will load and visualize different data % types acquired from different modalities. % % ===================== Download Data ============================== % % First, download and unzip the BioData from Canvas: % % % The data folder contains 5 files: % % (1) emg_myopathy.txt (EMG) % % (2) bMode.mat (B-mode Ultrasound) % % (3) vol4d.mat (fMRI) % % (4) XYStripes20x20Noise_9_L_3.mat (tensor field) % % (5) clathrin_ROI_cmpt340.csv (SMLM data) clc; disp('Hamoudi Ohan Saleh Baratta') disp('SFU ID # 301540229') % %% ============================ Task 1 =============================== % (1) EMG Data % Data Description: EMG is short for electromyogram, which measures electrical activity in muscles % Requirement: % Load the data from emg_myopathy.txt and extract the two columns of values: % i) First column -- time % ii) Second column -- the amplitude of an EMG signal. close all; clear; % Load data into MATLAB: % Write your own code here: EMG = importdata('emg_myopathy.txt'); % Q1: What are the dimensions of the loaded variable (eg. 5 x 3 x 10)? % Write your own code here: sizeEMG = size(EMG); fprintf('\nVariable dimensions:\n %d x %d\n\n', sizeEMG(1), sizeEMG(2)); % How many EMG amplitude values does the signal contain? % Write your own code here: fprintf('Number of amplitudes:\n %d values\n\n', sizeEMG(1)); % Plot the amplitude vs time using the values in the two columns. % Properly label both axes. Add a suitable title to your figure. figure(1)

time = EMG(:, 1); amplitude = EMG(:, 2); plot(time, amplitude); xlabel('Time (s)'); ylabel('Amplitude (mV)'); title('EMG Signal'); % There is a period of time with a clearly reduced EMG amplitude. % Examine the previous plot, use the "Zoom In" icon on the toolbar % to close up on this region of the plot. % Visually (roughly) identify the approximate starting and ending % time for this period. % Display these values : % Write your own code here fprintf('Start:\n 16.7667\n\n'); fprintf('End:\n 16.9965\n\n'); % What are the time and amplitude values of the 100th sample of the data? % Write your own code here disp('Time_100th = '); display(EMG(100, 1)); disp('Amp_100th = '); display(EMG(100, 2)); % Plot the EMG signal from the 100th sample (in (1d)) up until the 1100th % sample. Give a title to your figure. figure(2) time_100_1100 = EMG(100:1100, 1); amplitude_100_1100 = EMG(100:1100, 2); plot(time_100_1100, amplitude_100_1100); xlabel('Time (s)'); ylabel('Amplitude (mV)'); title('EMG Signal from the 100th to 1100th sample'); %% ============================ Task 2 =============================== % (2) B-mode Ultrasound sequence % Data Description: This data represent an B-mode ultrasound video sequence (B stands for Brightness) % Requirement: Load the data in bMode.mat into MATLAB, explore the video frames and create a GIF % close all; clear; % Load data into MATLAB: % Write your own code here: load("bMode.mat"); % What are the dimensions of the loaded variable (eg. 5 x 3 x 10)? % You should output the variable dim_bMode as a vector of size 1x3.. % Write your own code here disp('bMode dimensions: '); % display your answers dim_bMode = size(bMode); % get the dimensions of bMode disp(dim_bMode); % display the dimensions

% How many frames (2D images) do we have in the video? % Write your own code here disp('Number of frames: '); % display your answer num_frames = dim_bMode(3); % get the number of frames from the third dimension disp(num_frames); % display the number of frames % What is the size (rows and columns) of each frame? % Write your own code here disp('Nb_rows:'); nb_rows = dim_bMode(1); % get the number of rows from the first dimension disp(nb_rows); % display the number of rows disp('Nb_cols:'); nb_cols = dim_bMode(2); % get the number of columns from the second dimension disp(nb_cols); % display the number of columns % Extract and display the 9-th frame of the sequence in a figure. % Apply the command that ensures that pixels are not stretched % (i.e. maintain aspect ratio of the image). Apply the command % that uses a grey colour-map to display the image. % Write your own code here figure(1) frame_9 = bMode(:, :, 9); % extract the 9-th frame of the sequence imagesc(frame_9); colormap gray; axis image; % display the frame with the grey colour-map and aspect ratio % What is the intensity value at row 30 and column 40, at frame 20? % Write your own code here disp('Intensity_(30,40,20): '); % update the display to show your answer intensity_30_40_20 = bMode(30, 40, 20); % get the intensity value at row 30, column 40, and frame 20[^23^][23] disp(intensity_30_40_20); % display the intensity value % Extract a cropped rectangular region of frame 15. % The cropped part should extend from row 10 to row 30, and from column 45 % to column 60 of the image. Display this cropped region % (also called region of interest or ROI). % Write your own code here figure(2) % update the colormap and figure title frame_15 = bMode(:, :, 15); % extract the 15-th frame of the sequence roi = frame_15(10:30, 45:60); % extract the cropped region of interest imagesc(roi); colormap gray; axis image; % display the region of interest title('Region of interest of frame 15'); % add a title to the figure % Create an animated GIF that shows the ultrasound sequence % (like a video clip). % Save this GIF image using the name US_seq.gif . % You can make use of the MATLAB code used during the lectures % Write your own code here filename = 'US_seq.gif'; % the name of the GIF file for idx = 1:num_frames % loop over all the frames frame = bMode(:, :, idx); % get the current frame [A, map] = gray2ind(frame, 256); % convert the frame to indexed image

Your preview ends here

Eager to read complete document? Join bartleby learn and gain access to the full version

Access to all documents
Unlimited textbook solutions
24/7 expert homework help

if idx == 1 % for the first frame imwrite(A, map, filename, 'gif', 'LoopCount', Inf, 'DelayTime', 0.1); % write the frame to the GIF file with infinite loop and 0.1 second delay else % for the subsequent frames imwrite(A, map, filename, 'gif', 'WriteMode', 'append', 'DelayTime', 0.1); % append the frame to the GIF file with 0.1 second delay end end %% ============================ Task 3 =============================== % (3) Functional MRI % Data Description: This is a 3D+time f-MRI data (functional magnetic resonance image). % fMRI captures the changes in oxygenation with time (4th dimension) at every (x,y,z) location of a volume. % Requirement: Load the data in vol4d.mat into MATLAB, explore how oxygenation level changes with time % close all; clear; % Load data into MATLAB: % Write your own code here: load('vol4d.mat'); % load the data % What are the dimensions of the loaded variable (eg. 5 x 3 x 10)? % Write your own code here: disp('Size fMRI: '); % display your answer disp(size(vol4d)) % display the dimensions % Noting that the first 3 dimensions represent the spatial x, y, and z % dimensions, identify how many rows, columns, slices, and time steps % does this fMRI data have? % Write your own code here: % update the display: disp('Rows fMRI: '); % display # of rows in each slice disp(size(vol4d,1)); disp('Columns fMRI: '); % display # of columns in each slice disp(size(vol4d,2)); disp('Slices fMRI: '); % display # of slices disp(size(vol4d,3)); disp('Time fMRI: '); % display time disp(size(vol4d,4)); % Plot a curve that shows how the oxygenation level changes % with time at voxel location (x=32, y=34 , z=10). % Define axis and figure title. figure(1) % create a new figure plot(squeeze(vol4d(32,34,10,:))) % plot the oxygenation level over time at (x=32, y=34, z=10) xlabel('Time') % label the x-axis ylabel('Oxygenation level') % label the y-axis title('Oxygenation level at (x=32, y=34, z=10)') % add a title % Extract and display the slice in the X-Y plane, at slice % location z=15, of the 10th time sample.

figure(2) % create a new figure slice_xy = squeeze(vol4d(:,:,15,10)); % extract the slice in the X-Y plane at z=15 and t=10 imagesc(slice_xy) % display the slice as an image axis image % set the aspect ratio to be equal colormap gray % choose a gray color map % Extract and visualize the slice in the X-Z plane, % at slice location y=10, of the 15th time sample. % Create an animated GIF that shows how the slice in the previous question % changes with time (over the whole time steps). % Save this GIF image using the name fMRI.gif. % You can make use of the MATLAB code used during the lectures. % Write your own code here: figure(3) % create a new figure filename = 'fMRI.gif'; % specify the name of the GIF file for t = 1:size(vol4d,4) % loop over the time steps slice_xz = squeeze(vol4d(:,10,:,t)); % extract the slice in the X-Z plane at y=10 and t=t imagesc(slice_xz) % display the slice as an image axis equal % set the aspect ratio to be equal colormap gray % choose a gray color map drawnow % update the figure frame = getframe(3); % capture the current figure as a frame im = frame2im(frame); % convert the frame to an image [imind,cm] = rgb2ind(im,256); % convert the image to indexed color if t == 1 % if this is the first frame imwrite(imind,cm,filename,'gif','Loopcount',inf); % write the image to the GIF file with infinite loop else % if this is not the first frame imwrite(imind,cm,filename,'gif','WriteMode','append'); % append the image to the GIF file end pause(0.1) % pause for 0.1 seconds end %% ============================ Task 4 =============================== %% (5) Diffusion Tensor MRI % Data Description: This is a single 2D slice of a Diffusion Tensor MRI data, % which is a 2D image with a 3x3 matrix at every pixel. % Requirement: Load the data XYStripes20x20Noise_9_L_3.mat into MATLAB, calculate the eigenvalues. close all; clear; % Load data into MATLAB: % Write your own code here: load('XYStripes20x20Noise_9_L_3.mat'); % What is the dimensions of the loaded variable (eg. 5 x 3 x 10)? % Write your own code here:

disp(['Size DT: ', num2str(size(TensorIM))]); %display your answer % How many rows and columns of pixels does this 2D image have? % Write your own code here: disp(['Nb rows: ', num2str(size(TensorIM,1)), ', Nb cols: ', num2str(size(TensorIM,2))]); % update to display the correct values % Extract and display (as an array of numbers) the 3x3 % matrix at pixel location: (row 10, column 15). Use variable name *arr3x3* % your code here, update the display to print your variables arr3x3 = squeeze(TensorIM(10,15,:,:)); % extract the 3x3 matrix disp('Array at x=10, y=15 : '); disp(arr3x3); % display the matrix % Array at x=10, y=15 % Use the "eig" function in MATLAB to calculate and display the 3 eigen % vectors and the corresponding 3 eigen values of the 3x3 matrix arr3x3; % Use the variable name *arr3x3_eig*. % your code here, add a display at the end of your code to print out your % answers [arr3x3_eig, arr3x3_val] = eig(arr3x3); % compute the eigenvectors and eigenvalues disp('Eigenvectors : '); disp(arr3x3_eig); % display the eigenvectors disp('Eigenvalues : '); disp(arr3x3_val); % display the eigenvalues % Create and display a 2D image *DT_2D* with a scalar number at every pixel. % The scalar number should be the entry (row=1, column=2) of the % 3x3 matrix at every pixel. % Write your own code here: DT_2D = squeeze(TensorIM(:,:,1,2)); % extract the entry (1,2) of each 3x3 matrix figure(1); imagesc(DT_2D); % display the 2D image colorbar; % add a colorbar axis equal; % maintain the aspect ratio title('2D image of entry (1,2) of each 3x3 matrix'); % add a title xlabel('Column Index'); % label the x-axis ylabel('Row Index'); % label the y-axis % Create and display a 2D image with a scalar number at every pixel. % The scalar number should be the largest eigenvalue of the 3x3 matrix at % each pixel. % The challenge is to do this WITHOUT using loops (without for or while). % your code here [~, D] = eig(arr3x3); % Compute eigenvalues eigenvalues = max(D, [], 3); % Take the largest eigenvalue % Reshape back to 2D image size eigenvalues_image = reshape(eigenvalues, size(arr3x3, 1), size(arr3x3, 2));

Your preview ends here

Eager to read complete document? Join bartleby learn and gain access to the full version

Access to all documents
Unlimited textbook solutions
24/7 expert homework help

% Create a 2D image with the largest eigenvalue at each pixel % Display the image figure(2); imagesc(eigenvalues_image); colorbar; % Add a colorbar for reference title('Largest Eigenvalue of Diffusion Tensor'); xlabel('X'); ylabel('Y'); %% ============================ Task 5 =============================== %% (6) SMLM -- Single-Molecule Localization Microscopy % Data Description: 3D point cloud data illustrating clathrin protein locations. % The 3 channels in the 2nd column represents 3D coordinates. close all; clear; % Load data into MATLAB: % Write your own code here: % What is the dimensions of the loaded variable (eg. 5 x 3 x 10)? % your code here my_data = load('clathrin_ROI_cmpt340.csv'); % Determine the dimensions of the loaded variable sz = size(my_data); % Display the size of the SMLM data disp(['Size SMLM: ', num2str(sz(1)), ' x ', num2str(sz(2))]); % Define the indices for selecting the middle 1001 points start_index = floor((size(my_data, 1) - 1001) / 2); end_index = start_index + 1000; % Select the 1001 points from the imported data selected_points = my_data(start_index:end_index, :); % Plot the raw data figure(1); scatter3(selected_points(:, 1), selected_points(:, 2), selected_points(:, 3), 10, 'filled'); title('Raw Data: Clathrin Protein Locations'); xlabel('X (nm)'); % Label X dimension with unit ylabel('Y (nm)'); % Label Y dimension with unit zlabel('Z (nm)'); % Label Z dimension with unit grid on; colormap('jet'); % Use colormap 'jet' set(gca, 'Color', 'k'); % Set black background view(3); % Rotate the camera to visualize the data %% Construct graph (this method only works for small data) % this is the naive way to construct the graph using the distance matrix (for small data) % set proximity threshold -- proximity threshold means how you decide to clip the data according to the distances among neighboring points % Calculate pairwise distance among all the points

% Use adjacency matrix representation % Threshold the adjacency distance matrix % Create a graph from the distance matrix, save it as variable g % Step 1: Set the proximity threshold PT = 800; % Proximity threshold in nanometers % Step 2: Calculate pairwise distances distances = pdist(selected_points); % Step 3: Create the adjacency matrix adj_matrix = squareform(distances); % Step 4: Threshold the adjacency matrix adj_matrix(adj_matrix > PT) = 0; % Step 5: Create a graph from the distance matrix g = graph(adj_matrix); % Optional: Visualize the graph figure(2); plot(g, 'Layout', 'force', 'NodeColor', 'r', 'EdgeAlpha', 0.2); title('Graph of Clathrin Protein Locations'); %% What if you selected more than 1001 points from the data % When the size of the data is large, creating a graph from % the (thresholded) distance matrix can produce out-of-memory errors. % An alternative is to call the graph by supplying the source, the target, % and the weights of the edges in the graph, i.e., % g = graph(source, target, weight); % % Write the code (which involves a loop over the number of nodes) % that can create source, target and weight from d=pdist(rawDat); % Write your code here % Calculate pairwise distances among all points distances = pdist(selected_points); % Determine the number of nodes (data points) num_nodes = size(selected_points, 1); % Initialize arrays for source, target, and weights source = []; target = []; weight = []; % Loop over nodes to create edges for i = 1:num_nodes for j = i+1:num_nodes % Calculate the distance between nodes i and j dist_ij = distances((i-1)*(num_nodes-i/2) + j - i); % Check if the distance is below the proximity threshold if dist_ij <= PT % Add edge (i, j) with weight dist_ij source = [source; i]; %#ok<*AGROW> target = [target; j]; weight = [weight; dist_ij];

end end end % Create the graph using specified edges g = graph(source, target, weight); % Optional: Visualize the graph figure(3); plot(g, 'Layout', 'force', 'NodeColor', 'b', 'EdgeAlpha', 0.2); title('Graph of Clathrin Protein Locations (Alternative Construction)'); %% plot the constructed graph preserving the spatial info % p= plot(g, 'XData',...,'YData',...,'ZData',...,'NodeLabel',{}); Fill in the code here % axis equal; % set marker to 'o', use a marker size=2, red color and white edge color % set X, Y, Z axis labels % Write your code here % Set the XData, YData, and ZData for the graph (you can adjust these values) XData = selected_points(:, 1); YData = selected_points(:, 2); ZData = selected_points(:, 3); % Plot the graph figure(4); p = plot(g, 'XData', XData, 'YData', YData, 'ZData', ZData, 'NodeLabel', {}); axis equal; % Set X, Y, and Z axis labels xlabel('X (nm)'); ylabel('Y (nm)'); zlabel('Z (nm)'); title('Constructed Graph: Clathrin Protein Locations'); %% find the node degree, save it into variable deg % Write your code here % plot the points colored by their degree % Write your code here % use colormap jet, black background, set X, Y, Z axis labels using unit='nm' % Write your code here % Assuming your graph 'g' is already constructed % Step 1: Find the node degree deg = degree(g); % Step 2: Plot the points colored by their degree figure(5); scatter3(XData, YData, ZData, 10, deg, 'filled'); colormap(jet); % Use colormap 'jet' set(gca, 'Color', 'k'); % Set black background xlabel('X (nm)'); ylabel('Y (nm)');

Your preview ends here

Eager to read complete document? Join bartleby learn and gain access to the full version

Access to all documents
Unlimited textbook solutions
24/7 expert homework help

zlabel('Z (nm)'); title('Node Degree: Clathrin Protein Locations'); colorbar; % Add colorbar %%%%%%%%%%%%%% % For the original graph (created using the adjacency matrix) change graph visualization layout to 'force' % Write your code here % Assuming your graph 'g' is already constructed % Change the graph layout to 'force' figure(6); plot(g, 'Layout', 'force', 'NodeColor', 'b', 'EdgeAlpha', 0.2); title('Graph Layout: Clathrin Protein Locations (Force-Directed)');

Assign_Read_Vis_Biosignals_Images

Related Documents