Assign_Read_Vis_Biosignals_Images-KEY

School

Simon Fraser University *

*We aren’t endorsed by this school

Course

340

Subject

Electrical Engineering

Date

Apr 3, 2024

Type

Pages

Uploaded by DeaconLion805

%% CMPT340 Assignment - Reading & Visualizing Bio-signals % THIS IS A TEMPLATE SOLUTION, DO NOT SHARE OR REPRODUCE WITHOUT CONSENT % % ============================ 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) %% ============================ Task 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; clc; % Load data into MATLAB: % Write your own code here: emg_path = 'BioData/emg_myopathy.txt'; emg_signal = readmatrix(emg_path); % What are the dimensions of the loaded variable (eg. 5 x 3 x 10)? % Write your own code here: [x_dim, y_dim] = size(emg_signal); disp(['Variable dimansions: ', num2str(x_dim), ',' , num2str(y_dim)]); % How many EMG amplitude values does the signal contain? % Write your own code here: emg_amplitude_vals = length(unique(emg_signal(:,2))); disp(['Number of amplitude values: ', num2str(emg_amplitude_vals)]); % 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) plot(emg_signal(:,1), emg_signal(:,2)); title('Amplitude vs. Time'); xlabel('EMG Time') % x-axis label ylabel('EMG Amplitude') % y-axis label

% 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 : start_pt = 16.19; end_pt = 17.0; display(['Start point: ', num2str(start_pt), ', ', 'End point: ',num2str(end_pt) ]) % What are the time and amplitude values of the 100th sample of the data? % Write your own code here display(['Time_100th = ', num2str(emg_signal(100,1))]) display(['Amp_100th = ', num2str(emg_signal(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) plot(emg_signal(100:1100,1),emg_signal(100:1100,2) ) xlabel('EMG Time') % x-axis label ylabel('EMG Amplitude') % y-axis label title('emg signal between 100-1100 sample') %% ============================ Task 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; clc; % Load data into MATLAB: % Write your own code here: bMode_path = 'BioData/bMode.mat'; load(bMode_path); % What are the dimensions (*dim_bMode*) 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 [dim1 dim2 dim3] = size(bMode); disp(['bMode dimensions: ', num2str(dim1),',', num2str(dim2),',',num2str(dim3)]); % How many frames (2D images) do we have in the video? % Write your own code here disp(['Number of frames: ', num2str(dim3)]); % What is the size (rows and columns) of each frame? % Write your own code here disp(['Nb_rows: ', num2str(dim1)]); disp(['Nb_cols: ', num2str(dim2)]); % 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) imagesc(bMode(:,:,9)); title('9th frame') axis image colormap gray; % What is the intensity value at row 30 and column 40, at frame 20? % Write your own code here disp(['Intensity (30,40,20): ', num2str(bMode(30,40,20))]); % 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 ROI = squeeze(bMode(10:30, 45:60, 15)); figure(2) % update the colormap and figure title imagesc(ROI); title('ROI frame 15'); axis image colormap gray % Create an animated GIF that shows the ultrasound sequence (like a video clip). % You need to submit 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_seg.gif'; allFrames = bMode; allFrames = allFrames - min(min(min(allFrames))); allFrames = allFrames/(max(max(max(allFrames)))); % adjust range to be between 1 and 256 allFrames = allFrames*255 + 1; imwrite(squeeze(allFrames(:,:,1)),filename,'gif', 'WriteMode','overwrite','LoopCount',Inf); for ii = 2:size(allFrames,3) imwrite(squeeze(allFrames(:,:,ii)),filename,'gif', 'WriteMode','append'); end %% ============================ Task 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; clc; % Load data into MATLAB: % Write your own code here: fmri_path = 'BioData/vol4d.mat'; load(fmri_path) fmri = vol4d;

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

% What are the dimensions of the loaded variable (eg. 5 x 3 x 10)? fmri_dim = size(fmri); disp(['Size fMRI: ', num2str(fmri_dim) ]); % 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? fMRI_nb_rows = size(fmri,1); fMRI_nb_cols = size(fmri,2); fMRI_nb_slices = size(fmri,3); fMRI_time = size(fmri,4); disp(['Size fMRI: ', num2str(fMRI_nb_rows) ]); disp(['Size fMRI: ', num2str(fMRI_nb_cols) ]); disp(['Size fMRI: ', num2str(fMRI_nb_slices) ]); disp(['Size fMRI: ', num2str(fMRI_time) ]); % 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) plot(squeeze(fmri(32,34,10,:)), '-*'); title('Oxygenation at voxel (32,34,10)'); xlabel('Time') ylabel('Oxygenation'); % Extract and display the slice in the X-Y plane, at slice location z=15, of the 10th time sample. figure(2) imagesc(squeeze(fmri(:,:,15,10))); axis image title('XY plane at location z=15') % Extract and visualize the slice in the X-Z plane, at slice location y=10, of the 15th time sample. figure(3) imagesc(squeeze(fmri(:,10,:,15))); title('XZ plane at location y=10') % Create an animated GIF that shows how the slice at y=10 % changes with time (over the whole time steps). % You need to submit this GIF image using the name fMRI.gif. % You can make use of the MATLAB code used during the lectures. filename = 'fMRI.gif'; allFrames = fmri; img = imagesc(squeeze(allFrames(:,10,:,1))); f = getframe; [im,map] = rgb2ind(f.cdata,256); imwrite(im,map,filename,'gif','writemode','overwrite');

for ii=2:size(fmri,4) img = imagesc(squeeze(allFrames(:,10,:,ii))); f = getframe; [im, map] = rgb2ind(f.cdata,256); imwrite(im,map,filename,'gif','writemode','append') end %% ============================ Task 4 =============================== % 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; clc; % Load data into MATLAB: % Write your own code here: tensor_path = 'BioData/XYStripes20x20Noise_9_L_3.mat'; load(tensor_path); DtMRI = TensorIM; % What is the dimensions of the loaded variable (eg. 5 x 3 x 10)? [x y z c] = size(DtMRI); disp(['Size fMRI: ', num2str(x), ',', num2str(y), ',' ,num2str(z), ',', num2str(c)]); % How many rows and columns of pixels does this 2D image have? disp(['Nb rows: ', num2str(x),', Nb cols: ',num2str(y)]); % Extract and display (as an array of numbers) the 3x3 matrix at pixel location: (row 10, column 15). Use variable name *arr3x3* arr3x3 = squeeze(DtMRI(10,15,:,:)); disp(arr3x3); disp('Array at x=10, y=15 : '); disp(num2str(arr3x3)); % Can be removed since it is repeated from above % 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*. [arr3x3_eigVec,arr3x3_eigVal] = eig(arr3x3); disp('Eigenvectors : '); disp(num2str(arr3x3_eigVec)); disp('Eigenvalues : '); disp(num2str(arr3x3_eigVal)); % 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. DT_2D = squeeze(DtMRI(:,:,1,2)); figure(1) imagesc(DT_2D) axis image

title('DT 2D') %% Bonus Challenge: % 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). % Hint: These MATLAB functions can help you for this part: mat2cell , reshape , cellfun , cell2mat. Note that this part of Task 4 is rather challenging, so don't feel bad if you can't get it to work but do give it a try. It is not the type of question you will see in the exam. n_rows = size(DtMRI,1); n_cols = size(DtMRI,2); DT_3D = mat2cell(reshape(DtMRI,[n_rows*n_cols, 3, 3]), ones(1,n_rows*n_cols)); % (20*20)*3*3 DT_3D = cellfun(@squeeze, DT_3D,'unif', false ); DT_eigs = cellfun(@eig, DT_3D,'unif', false ); DT_eigs_max = cellfun(@max, DT_eigs, 'unif', false); DT_eigs_max = cell2mat(DT_eigs_max); DT_eigs_max = reshape(DT_eigs_max,[n_rows, n_cols]); figure(2) axis image imagesc(DT_eigs_max); title('BONUS'); colormap; %% ============================ Task 5 =============================== % 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; clc; % Load data into MATLAB: % Write your own code here: fName = 'BioData/clathrin_ROI_cmpt340.csv'; % file name dat = importdata(fName); rawDat = dat.data; % What is the dimensions of the loaded variable (eg. 5 x 3 x 10)? % your code here disp(['Size SMLM: ', num2str(size(rawDat,1)), ', ', num2str(size(rawDat, 2))]); %display your answer % Select 1001 points from the imported data. You can clip from the middle, e.g., use format such as rawData(ind-xxx:ind:xxx) % Write your code here rawDat=rawDat(round(end/2)-500:round(end/2)+500,:); % Plot the raw data figure(1) % you can use proper 'position' property to make the figure look nice % Write your code here % set(gca, 'position', [60, 60, 800, 800]); scatter3(rawDat(:,1),rawDat(:,2),rawDat(:,3), 2, 'g', 'filled') % Use colormap jet, black background (hint you can use MATLAB's command "look for": lookfor 'background color') % Label X, Y, Z dimensions properly, display Unit = 'nm' % rotate camera to visualize your data

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

axis equal; colormap(jet); % axis off; whitebg('black'); set(gcf,'color', 'k') xlabel ('X nm'); ylabel ('Y nm'); zlabel ('Z nm'); grid on; rotate3d; view(0,-90); % view in 2D % 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) PT = 800; % 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 distDat = pdist(rawDat); distMat = squareform(distDat); distMat_th = distMat.*(distMat <= PT); % thresholded distance matrix g = graph(distMat_th); % create the graph from the distance matrix % Bonus question, 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 PT = 800; % proximity threshold d=pdist(rawDat); distMat = squareform(distDat); source = []; target = []; weight = []; i = 1; for RI = 1:size(distMat, 1) % row index for CI = RI+1:size(distMat, 2) % column index if distMat(RI, CI) <= PT && distMat(RI, CI)>0 % go through every element inside d and select only the distance below the threshold PT source(i)=CI; target(i)=RI; weight(i)=d(i); i = i + 1; end end end g = graph(source', target', weight'); % plot the constructed graph preserving the spatial info % Use the smaller subset of the data (i.e., as was instructed in the line "Select

1001 points from the imported data...") figure(2) % set(gca, 'position', [60, 60, 800, 800]); p= plot(g,'XData',rawDat(:,1),'YData',rawDat(:,2),'ZData',rawDat(:,3), 'NodeLabel', {}); % 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 axis equal; p.Marker = 'o'; p.MarkerSize = 2; % change the node size p.NodeColor = 'r'; % change the node color p.EdgeColor = [1 1 1]; % change the edge color xlabel ('X (nm)'); ylabel ('Y (nm)'); zlabel ('Z (nm)'); grid on; view(3); % find the node degree, save it into variable deg % Write your code here deg = degree(g); % plot the points colored by their degree figure(3) % Write your code here % set(gca, 'position', [60, 60, 800, 800]); scatter3(rawDat(:,1),rawDat(:,2),rawDat(:,3), 3, deg, 'filled') axis equal; % use colormap jet, black background, set X, Y, Z axis labels using unit='nm' % Write your code here colormap(jet); % axis off; whitebg('black'); set(gcf,'color', 'k') xlabel ('X nm'); ylabel ('Y nm'); zlabel ('Z nm'); grid on; rotate3d; %view(0,-90); % view in 2D %%%%%%%%%%%%%% figure(4) % change graph visualization layout to 'force' plot(g,'Layout','force')