HST.583

LAB6: Cortical and Subcortical Parcellation with MRI

December 2, 2002


Contents


1. Goals of lab6

This Lab examines ways in which different brain anatomical structures can be classified based on the signal intensity of high spatial resolution anatomical MR images acquired with different contrast weightings.


2. Organization of the lab

Here is an outline (details are below):


3. Getting started

login to your athena account

If you don't have one, create a matlab directory in your home directory: mkdir ~/matlab In this directory you should write the matlab functions that you can to use for solving the exercises.


% attach hst.583
% cd /mit/hst.583/lab_sw/lab6
% add matlab
% matlab &
Within the MatLab window, type:
>>Lab6_setpath

This changes your directory to that of lab6 and adds the paths for matlab functions used in lab6. (Your own ~/matlab directory is automatically in your path. (~/matlab/foo, I dunno. --E'beth)) It also runs setup.m, which defines some symbolic constants for structure names and runs load_data.m, which loads three volumes:

T1vol
Map of T1 values across the brain - That is, this is a 3D volume image of the brain where the value of each voxel is the T1 (in ms) of the tissue in that voxel.
PDvol
Map of proton density values across the brain. As for the T1 map but now the value in each voxel is the proton density (arbitrary units) in that voxel.
LABELvol
volume of labels of all neuroanatomical structures

If you type whos you will get a list of the variables hold in memory.

You can use the additional tool show_slice (see below) to look at the images you just loaded.

Now read carefully the exercises. (They also appear in exercises.m, a matlab file which consists the line "setup" and many lines of matlab comments reproducing the exercises below.) To answer the questions write matlab functions in your home matlab directory (perhaps taking exercises.m as a departure point).


4. Additional matlab tools

The following matlab files have also been provided:

load_data.m
load the data needed for this lab
setup.m
calls load_data and defines some structure constants
load_mgh.m
load a volume in 'mgh' format
radians.m
convert degrees to radians
save_mgh.m
save a volume in 'mgh' format
show_slice.m
show a (coronal) slice from a volume

If you want to save images, do so in your own directory. For example:

save_mgh(my_image,'~/hst-583/lab6/my_image.mgh')

Note that to display a 2-d slice of a volume, you can use the function show_slice(volume, sliceno). This will display coronal slices, but you should be able to modify it to view slices along any other cardinal axis easily enough.


5. Excercises

EXCERCISE 1a:
Write a matlab script that will use the T1 and PD maps to synthesize a new image volume using a flip angle of a = 30 degrees (remember to use radians for matlab sin and cos) and a TR of 20 msec (assume that TE is short relative to T2 and can therefore be neglected). Call this new volume synthesized_vol.

The equation you need is:

E1 = exp(-TR/T1)

E2 = exp(-TE/T2)

S(TR, TE, a, T1, T2, PD) = PD * sin a * ((1-E1) / (1 - (cos a) * E1)) * E2

What is the primary 'weighting' of this image? That is, what tissue parameter is primarily responsible for the contrast in the image? Why?

EXCERCISE 1b:
Compute the contrast to noise ratio of each of the following pairs of structures using the synthesized volume obtained in the previous point.

You may find the matlab function 'find' useful in this context. For example:

would compute the mean intensity (amygdala_mean) and the variance (amygdala_var) of all voxels labeled amygdala in the volume you have synthesized.

Generate a table with these contrast to noise values and turn it in.

EXCERCISE 1c:
Turn in a plot of the empirical distribution of the signal intensity of voxel in each structure using the synthesized image. (hint: take advantage of matlab's histogramming function (which may or may not be called "hist" --E'beth).) Remember to normalize the empirical distributions so they sum to 1. Please use the following colors:

Comment on your results with regard to building an automated segmentation procedure based on this type of intensity image.

EXCERCISE 1d:
Compute the optimal discriminant between gray matter and white matter assuming that the classes have a normal distribution with equal variances. Plot the empirical distributions of each on the same plot, together with a line representing the decision boundary. Compute the percentage of voxels of each class that would be misclassified using this decision rule. How well placed is the discriminant based on the Gaussian models?

EXCERCISE 1e:
Now compute the contrast-to-noise ratio for a subset of cortical gray matter and white matter. Use a single horizontal slice near the center of the 30 degree flip angle volume (use slice 90, meaning vol(:,90,:)). Compare it to the gray/white contrast-to-noise ratio for the entire volume. Speculate on the cause of the disparity between the two (hint: you may want to take a look at this data at different horizontal slice planes). Recalculate the discriminant as in 1d and percent voxels misclassified using just this slice. How well placed is the discriminant based on the Gaussian models?

EXCERCISE 2a:
Change the flip angle to 5 degrees and synthesize a new volume. Use this volume to recalculate the contrast-to- noise ratio of the same pairs of structures (neglecting T2 effects again). What is the primary 'weighting' of this image? Why is it different from the 30 degree image of exercise 1a?

EXCERCISE 2b:
Compute the contrast to noise ratio of the same structures as 1b using this synthesized volume. Compare them to your results in 1b. Which type of volume would you choose for segmentation purposes for each structure pair?

EXCERCISE 2c:
Plot the empirical distribution of each structure using this volume with the same colors as done in 1c.

EXCERCISE 2d:
Compute the optimal discriminant between gray matter and white matter assuming that the classes have a normal distribution with equal variances. Plot the empirical distributions of each on the same plot, together with a line representing the decision boundary. How well placed is the discriminant based on the Gaussian models? Compute the percentage of voxels of each class that would be misclassified using this decision rule. Compare this to your results in 1d. Which volume is more useful for gray/white classification?