% MATLAB TUTORIAL FOR THURSDAY, SEPTEMBER 2ND, 2010

% In this demo, you will:

% (1) load MRI data into MATLAB
% (2) display three anatomically oriented slices of the head MRI
% (3) apply a binary mask to extract the brain tissue
% (4) display three brain slices
% (5) implement an edge detector to grab and display all edges of a slice

% clear everything from memory
clear all

% close any figure windows
close all

% clear the command window
clc

% load mat file containing the original mri
load raw_mri % a file called 'raw_mri.mat' has been stored in the working directory

% see what's in memory
whos

% create a subplot window
subplot(2,2,1) 

% grab the coronal slice corresponding to y=128 and store it in a new
% variable
coronal_slice = mri(:,128,:);
% NOTE: colon operator denotes "all" elements in the given dimension
% NOTE: semicolon after line suppresses output

% get the size of the slice
size(coronal_slice)
% NOTE: no semicolon as we want to display the current size of the variable
% "coronal_slice"

% remove "singleton" dimension
coronal_slice = squeeze(coronal_slice);

% observe the size of the array after squeezing
size(coronal_slice)
% NOTE: "coronal_slice" is now a 2-dimensional array

% display slice
imagesc(coronal_slice);
% NOTE: intensity of each pixel corresponds to tissue type

% change to a grayscale colormap
colormap(gray)

% add a title to our subplot
title('Coronal Slice')


% display the coronal slice right side up

% create a second panel
subplot(2,2,2)

% display the transposed matrix
imagesc(coronal_slice');
% NOTE: in MATLAB, the prime operator ' denotes matrix transposition

% make the axes Cartesian
axis xy

% insert new title
title('Coronal Slice Right Side Up')


% display the sagittal slice
subplot(2,2,3)

% grab a slice corresponding to the y-z plane
sagittal_slice = mri(128,:,:);
sagittal_slice = squeeze(sagittal_slice);
imagesc(sagittal_slice');
axis xy
title('Sagittal Slice')

% display the axial slice
subplot(2,2,4)

% grab a slice corresponding to the x-y plane
axial_slice = mri(:,:,95);
axial_slice = squeeze(axial_slice);
imagesc(axial_slice');
axis xy
title('Axial Slice')


%%

% load mat file containing our data: mask of the brain
load brain_mask

% display the contents of memory
whos
% NOTE: the mask is stored in a variable called "m"

% get the size of the mask m
size(m)

% get the type of data contained in the "m" variable
class(m)

% get the type of data contained in the "mri" variable
class(mri)

% since the data types don't match up, we must convert one of them to match
% the other if we want to do any processing

% convert "mri" to a floating-point array (as opposed to integer data)
mri_new = double(mri); 

% compare the size of "mri_new" with that of "mri"
whos

% now observe the type of data in "mri_new"
class(mri_new)

% apply the mask to the mri image
masked_mri = mri_new.*m;
% NOTE: .* denotes element-by-element multiplication, which is not the same
% as matrix multiplication

% display three brain slices

% create a new figure window
figure
subplot(2,2,1) 
masked_coronal_slice = masked_mri(:,128,:);
masked_coronal_slice = squeeze(masked_coronal_slice);
imagesc(masked_coronal_slice');
axis xy
colormap(gray)
title('Masked Coronal Slice')

% display the sagittal slice
subplot(2,2,2)
masked_sagittal_slice = masked_mri(128,:,:);
masked_sagittal_slice = squeeze(masked_sagittal_slice);
imagesc(masked_sagittal_slice');
axis xy
colormap(gray)
title('Masked Sagittal Slice')

% display the axial slice
subplot(2,2,3)
masked_axial_slice = masked_mri(:,:,95);
masked_axial_slice = squeeze(masked_axial_slice);
imagesc(masked_axial_slice');
axis xy
colormap(gray)
title('Masked Axial Slice')

%%

% IN-CLASS ASSIGNMENT

%% 
% TASK 1
% Your job is to use the "edge" function of the image processing toolbox to
% find and display all images in the axial slice of the brain mask.  Once
% you have computed the edges, you should display the edgemap in the
% remaining subplot panel (bottom right) with the proper orientation.
% 
% type "help edge" in the command window to view the function prototype and
% argument information for the edge function


subplot(2,2,4)
bw=edge(masked_axial_slice);
imagesc(bw');
axis xy
title('edge map')

%% 
% TASK 2
% In a new figure window, display a histogram of MR intensities for all non-zero voxels.  To do
% this, you will first use the "find" function to identify the voxels that
% have a non-zero intensity.  Next, create a new vector to store the intensities
% of all non-zero voxels.  Finally, you will call the "hist" function with
% this vector as the argument to display the intensity histogram.
figure
hist(mri_new(find(mri_new>0)),1000);

%%
% TASK 3
%
% In a new figure window, iterate across the z-dimension (dimension 3 in
% the array mri_new) and display the value of each x-y slice successively.
% This requires the use of a "for-loop" as well as the "drawnow" function.  
% The pseudocode of the script is given below.  You need only fill in the
% question marks with the appropriate matlab variables/commands.  

%
figure
colormap(gray)
for k = 1:1:size(mri_new,3)
    imagesc(squeeze(mri_new(:,:,k))')
    drawnow
    axis xy
end
    
    
    









