Visualization

Class: Psych 204a

Tutorial: Visualization

Author: Dougherty, Bowen, Wandell

Date: 2010.10.21

Duration: 45 minutes

Copyright: Stanford University, Robert F. Dougherty, R Bowen & BA Wandell

Translated to Python by Grace Tang, 09/2013

Purpose:

The purpose of this tutorial is to illustrate the underlying visualization methods used in the analysis of MR data. The tutorial includes examples of tissue segmentation into gray matter and white matter, the definition of the surface boundary of the brain as a mesh, and related anatomical measurements, such as volume and surface area.



In [ ]:

    
# set up plotting, numpy/matplotlib, etc.
%pylab inline
import matplotlib as mpl
mpl.rcParams["figure.figsize"] = (8, 6)
mpl.rcParams["axes.grid"] = True
from mpl_toolkits.axes_grid1 import make_axes_locatable
from IPython.display import display

Load a high-resolution T1-weighted MRI dataset of someone's head:



In [ ]:

    
import scipy.io
import nibabel as nib

# load nifti files
anat_file = nib.load('t1.nii.gz')
brainMask_file = nib.load('t1_brain_mask.nii.gz')

Let's see what got loaded:



In [ ]:

    
anat = anat_file.get_data()
brainMask = brainMask_file.get_data()
print anat.shape
print brainMask.shape

The keys function displays the keys of the newly loaded python doctionary. anat contains the image data. We'll examine the brain mask later.

We can also find out the size of physical space (in millimeters) that each image voxel represents.



In [ ]:

    
anat_file.get_header().get_zooms()

In this case, the sides of the voxel measure about 0.8mm.

Question 1

a. What is the field of view of this image?

b. What is the volume of one voxel?

Examine a slice from this dataset

The T1-weighted MR data are stored in a 3d array called 'anat'. Here's one slice from the middle of this array:



In [ ]:

    
z=90 # slice number

fig, ax = plt.subplots(1)
ax.matshow(anat[:,:,z])

We'll format the image to make it look nicer:



In [ ]:

    
fig, ax = plt.subplots(1)
images=[]
images.append(ax.matshow(anat[:,:,z], cmap=matplotlib.cm.bone))
ax.set_axis_off()
ax.set_axes('tight')

# add a colorbar
ax = plt.gca()
divider = make_axes_locatable(ax)
cax = divider.append_axes("right", size="5%", pad=0.05)
fig.colorbar(images[0], cax=cax)

In this T1-weighted image, tissues with long T1 values show up as dark (cerebro-spinal fluid, CSF), medium T1 values are gray (gray matter, GM, and muscle), and shorter T1 values are white (white matter, WM, subcutaneous fat in the scalp). In addition to T1-weighting, this image includes some proton density weighting. That's why bone and air, which have relatively little free water, also show up as dark. Because this image is not based on a quantitative MR method, the image intensity values are in arbitrary units, so only the realtive intensities (image contrast) are important.

Question 2

Try to identify the different tissue types in the image.

a. cerebro-spinal fluid (CSF)

b. White matter (WM)

c. Gray matter (GM)

d. bone

e. muscle

f. scalp

g. air

Now let's add some labels so that you can see how well you did:



In [ ]:

    
fig, ax = plt.subplots(1)
images=[]
images.append(ax.matshow(anat[:,:,z], cmap=matplotlib.cm.bone))
ax.set_axis_off()
ax.set_axes('tight')

# add a colorbar
ax = plt.gca()
divider = make_axes_locatable(ax)
cax = divider.append_axes("right", size="5%", pad=0.05)
fig.colorbar(images[0], cax=cax)

# add text
bg = [1,0,0,.4]
tx = [1,1,1,1]
ax.text(100,116,'CSF', fontsize=12, backgroundcolor=bg, color=tx)
ax.text(200,140,'WM', fontsize=12, backgroundcolor=bg, color=tx)
ax.text(112,48,'GM', fontsize=12, backgroundcolor=bg, color=tx)
ax.text(150,195,'bone', fontsize=12, backgroundcolor=bg, color=tx)
ax.text(170,35,'muscle', fontsize=12, backgroundcolor=bg, color=tx)
ax.text(80,200,'scalp', fontsize=12, backgroundcolor=bg, color=tx)
ax.text(20,30,'air', fontsize=12, backgroundcolor=bg, color=tx)

Working in three dimensions

Remember that this is a 3d dataset. Let's look at all the slices as a simple movie.



In [ ]:

    
import time, sys
from IPython.core.display import clear_output

fig, ax = plt.subplots(1)
for slice in range(anat.shape[2]):
    ax.matshow(anat[:,:,slice], cmap=matplotlib.cm.bone)
    time.sleep(0.1)
    clear_output()
    display(fig,ax)
    ax.cla()
plt.close()

We can also show the middle slice of each axis



In [ ]:

    
fig = plt.figure(figsize=(12,8))
ax1 = fig.add_subplot(1,3,3)
ax2 = fig.add_subplot(1,3,2)
ax3 = fig.add_subplot(1,3,1)
ax1.matshow(np.rot90(anat[anat.shape[0]/2,:,:]), cmap=matplotlib.cm.bone)
ax2.matshow(np.rot90(anat[:,anat.shape[1]/2,:]), cmap=matplotlib.cm.bone)
ax3.matshow(np.rot90(anat[:,:,anat.shape[2]/2]), cmap=matplotlib.cm.bone)

Segment the brain tissue

The first step to segmenting the interesting brain tissues is to remove the annoying scalp and skull. This is often called "brain extraction." It can be done manually, by carefully labeling each voxel. There are also some very good algorithms that can do this task automatically. We've applied one of these tools (FSL BET: http://www.fmrib.ox.ac.uk/fsl/bet2/) and saved the resulting brain mask.

Here is the brain mask for the same middle slice that we saw above. This 'mask' is a simple image with 1's indicating brain tissue and 0's indicating non-brain tissue.



In [ ]:

    
fig, ax = subplots(1)
images=[]
images.append(ax.matshow(brainMask[:,:,90], cmap=matplotlib.cm.bone))
ax.set_axis_off()
ax.set_axes('tight')

# add a colorbar
ax = plt.gca()
divider = make_axes_locatable(ax)
cax = divider.append_axes("right", size="5%", pad=0.05)
fig.colorbar(images[0], cax=cax)

Question 3

What is the volume of the brain in cubic centimeters? Hint: the volume of one voxel in cubic mm is prod(voxDim).

Now apply the brain mask by setting all non-brain voxels to 0 (multiply all non-brain voxels by 0):



In [ ]:

    
brain = anat*brainMask

# plot
fig, ax = subplots(1)
images=[]
images.append(ax.matshow(brain[:,:,90], cmap=matplotlib.cm.bone))
ax.set_axis_off()
ax.set_axes('tight')

# add a colorbar
ax = plt.gca()
divider = make_axes_locatable(ax)
cax = divider.append_axes("right", size="5%", pad=0.05)
fig.colorbar(images[0], cax=cax)

Let's look at a histogram of the image intensities in the brain.



In [ ]:

    
fig, ax = subplots(1)
hist(brain.ravel(), 100)
xlabel('Raw intensity'); ylabel('Pixel count');

There are many voxels with an intensity of zero, including all the stuff outside the brain that we set to zero above. Let's rescale the axes to make the more interesting features evident:



In [ ]:

    
fig, ax = subplots(1)
hist(brain.ravel(), 100)
xlabel('Raw intensity'); ylabel('Pixel count');
ax.set_ylim([0, 1e5])
ax.set_xlim([0, 10000])

Question 4

a. How many peaks are there in this histogram?

b. What tissue types do these peaks represent?

Here's a line that is likely to divide white and gray matter intensities



In [ ]:

    
wmThreshold = 6900

fig, ax = subplots(1)
hist(brain.ravel(), 100)
xlabel('Raw intensity'); ylabel('Pixel count');
# plot the white matter threshold
vlines(wmThreshold,0,1e5, color='red')
ax.text(wmThreshold-1500,10.2e4,'white matter threshold', fontsize=12, color='red')
ax.set_ylim([0, 1e5])
ax.set_xlim([0, 10000])

Segment out the white matter based on this boundary



In [ ]:

    
whiteMatter = brain > wmThreshold

Clean up the segmentation with a little smoothing. For a real segmentation, we would also want to have a human expert edit the segmentation to ensure that it is correct.



In [ ]:

    
from scipy.ndimage.filters import gaussian_filter

whiteMatter = gaussian_filter(double(whiteMatter), 1.5)>0.5

Inspect the segmentation

Now we overlay the whitematter segmentation over the brain.



In [ ]:

    
lighterBrain = brain / float(brain.max()) # reduce the intensity of the brain image

fig, ax = plt.subplots(1)
for slice in range(whiteMatter.shape[2]):
    # wherever there is white matter, the brain image is less intense
    ax.matshow(lighterBrain[:,:,slice]+whiteMatter[:,:,slice], cmap=matplotlib.cm.bone, vmax=1.0)
    time.sleep(0.001)
    clear_output()
    display(fig,ax)
    ax.cla()
plt.close()

Question 5

What is the volume of the white matter in cubic centimeters?



In [ ]:

    
# compute your answer here

View the brain surface

The notebook isn't quite ready for rendering 3d brains. So, to see what the 3d rendering of a brain looks like, visit the CBrain brain browser.