Segment Ex-vivo White Matter
Koen Van Leemput has created a GUI to help you view the segmentation as you are creating it. It also allows you to use intensity information from several different scans to find a good gm/wm threshold.
Threshold Out Noise
The first step is to threshold out structures that are not of any interest (air, plp, bag, etc.) and thus create a mask consisting of only the tissue classes you want to segment. By loading multiple volumes, you can interactively define threshold values on one volume and apply the resulting mask on another volume. Choose the volume(s) where plp, bag, and air are a different contrast from brain. It's okay if there is not a lot of contrast between wm and gm. The orig.mgz will be used as the volume to create the mask (to be used for segmentation later).
kvlThresholdImage orig.mgz /path/to/second_volume.mgz
Where it says Image to mask, choose the orig volume. Where it says Image to threshold, choose the other volume you loaded. Use the scroll bar to change the Lower threshold and Upper threshold until you see a purple overlay covering all the stuff that you do not want (you do not want to cover gm or wm).
Use the View section to choose whether you see all planes or just one.
- You can scroll through the slices using the scroll bars at the top right of the screen in order to see how much you are masking out on the other slices.
Check where it says Show inverse mask to see the actual mask that will be created. This will show everything that you will keep for the segmentation step.
You can change the Overlay opacity to see the volume underneath. The lower threshold is originally set to the minimum intensity value in the image (0) and the upper threshold defaults to the maximum intensity value. Whatever is covered by the overlay will be set to an intensity value of zero in the masked image.
Pressing the Mask image button will write out a file called "*_masked.mgz" (in this case orig_masked.mgz). This file will be the one you want to segment with in the next step (kvlEMSegment). (You can also choose to press the Write mask button to get a mask.mgz volume of just the mask and not the volume with noise masked out.)
Segment WM
The second step is to actually segment the wm from the gm. This will classify tissues using an automatically estimated Gaussian mixture and polynomial MR bias field model. You would use this command from the directory where the masked volume is:
kvlEMSegment orig_masked.mgz numberOfGaussians biasFieldOrder downSamplingFactor
- numberOfGaussians = number of tissue classes. (You can try 2 for wm and gm or 3 for wm, gm, and CSF/plp)
- biasFieldOrder = a number between 4 and 8. (The higher the #, the higher the spatial frequencies in MR inhomogeneities the method tries to compensate for)
- downSamplingFactor a number between 1 and 4. (The higher the #, the less voxels the method uses to estimate the model parameters: this speeds up computation time but possible at a loss of segmentation accuracy. For example, when set to 4 this tells the algorithm to only consider every 4th voxel in the x-, y-, and z-directions for calculating the model parameters, speeding the computations up by a factor of 4^3 = 64)
The suggested values to start with are:
kvlEMSegment orig_masked.mgz 2 4 4
**Every case is different so it is worth playing around with the different values above if you are unhappy with the segmentation (definitely try using 3 different tissue classes and/or a down samping factor of 1 to see if it improves your segmentation).
View Segmentation Results
Results are written out as crispPosterior_classX.mgz, one for each class (with 0 being wm and 1 being gm). To see the result of the automated segmentation overlaid on the original MR image, you can use:
kvlViewImage orig_masked.mgz crispPosterior_class0.mgz
Or if you have several wm segmentations you would like to view at the same time to see which is better, you can use Freeview:
freeview -v orig.mgz wm_file1.mgz:colormap=heatscale wm_file2.mgz:colormap=jet
If you are not satisfied with the segmentation, you can tweak the thresholds in the kvlThresholdImage step or try out different parameters in the kvlEMSegment step.