Ex-vivo segmentation
Author: Xiao Han (see also: XiaoNotes)
Important!!! This wiki uses commands that are no longer supported by FreeSurfer. For the most up to date instructions on reconning exvivo data, see the ExvivoRecon wiki.
This wiki page describes a procedure for cortical surface reconstruction from ExVivo brain scans.
The procedure is tailored for a particular ExVivo scan protocol used at MGH, which acquires a series of MEF scans for each hemisphere of a postmortem brain sample. The method has only been tested on a few data sets, and further development and improvement is still on the way.
The method contains two major stages. The first stage aims to get a good WM segmentation (with necessary manual editing involved). The second stage performs a multi-spetrcal surface deformation to find the gray/white and pial surfaces.
Stage 1
The first stage consists of the following steps:
Step 1. Preparing the data
The MEF protocol used for the ExVivo brain analysis generates a large number of image volumes. For example, typically 6 different flip-angles are used. Each flip-angle corresponds to multiple volumes acquired at different echo times. In addition, multiple acquisitions are typically performed for each flip angle. In the first processing step, we create a single average volume for each flip angle, by averaging across the multiple acquisitions and the multiple echo volumes. "mri_average" is sufficient for this purpose. Note that, in this step, we keep the original scan resolution and size by specifying "-noconform" in the "mri_average" command line. We will convert relevant volumes to the COR format later before the surface topology correction and surface deformation steps. This is based on the consideration that the extra interpolation (as performed when converting an input volume to the COR format) may interfere with the intensity clustering based tissue segmentation algorithm. Also, the COR format forces the intensity values to be of 8-bit BYTE type, which drops useful information from the original float-typed image. An example for generating the flash30 average volume is as follows:
mri_average -noconform mef30_echo?_avg.mgz mef30_avg.mgz
The same processs need be applied to all the other flip angles as well.
Step 2. Run intensity inhomogeneity correction using the "nu_correct" program
Although the EM segmentation method can perform simultaneous INU correction, correcting INU beforehand leads to a faster processing stream. We found that nu_correct works reasonably well for the MEF data. Two passes of nu_correct is sufficient to remove the intensity bias field in general. A script called "inu_correct_vol.csh" can be used for this process, which is a wrapper for the nu_correct and performs exactly two passes of nu_correction:
inu_correct_vol.csh mef30_avg.mgz mef30_avg_inu.mgz
The same procedure need be applied for all the other flip-angles as well.
Step 3. EM segmentation using mri_ms_EM
The following command is used:
mri_ms_EM -no_INU -T 200 \ -M 4 -beta 0.1 -noconform -synthonly \ mef30_avg_inu.mgz \ mef5_avg_inu.mgz \ mef10_avg_inu.mgz \ mef15_avg_inu.mgz \ $outdir/"
The option "-no_INU" indicates that no INU correction will be performed inside this EM segmentation. "-T 200" specifies a background threshold of 200. The threshold is applied to the first input volume (the mef30_avg_inu.mgz in this case) in order to define the brain mask (to avoid an additional skull-stripping step). Only voxels within the brain region (with intensity value > 200) will be considered in the following segmentation. The threshold need be set manually. "-M 4" indicates that the number of classes/labels considered in the tissue segmentation algorithm is 4. Although the purpose is to classify the brain voxels into one of three tissue classes: WM, GM, and CSF. Using a 4th class can accomodate remaining "background" voxels, and offers better segmentation results for the three major tissue classes. "-beta 0.1" is a parameter for the MRF model. Larger values give smoother segmentation results. "-noconform" indicates performing the segmentation in the original volume resolution. If "-conform" is specified, all input volumes will be first converted to the COR format before EM segmentation is performed.
The last argument "$outdir/" specifies the output prefix, usually a directory name such as "./". The program can have multiple outputs, which will have "$outdir/EM_" as their filename prefix.
For the above command, the major output is a synthesized gray-level byte volume saved as "$outdir/EM_combined.mgz". This volume encodes the fuzzy EM segmentation results as a gray-scale image. In particular, the intensity value at a image voxel is a weighted sum of the means of the 4 classes. The weights are equal to the posteriori probability of each voxel belonging to a particular label or class. The means for the 4 classes are pre-set to fixed values: 15 for class 0 (background), 85 for class 1 (CSF), 155 for class 2 (GM), and 225 for class 3 (WM). Consult the program for details.
Step 4. Get initial WM segmentation and final brain mask from the EM segmentation results
Again, since the EM segmentation results are encoded in the single "EM_combined.mgz" volume, simple thresholding and connected component extraction can be used to isolate the WM or the brain tissues. First, run
mri_extract_largest_CC -T 185 EM_combined.mgz wm_init_tmp.mgz
This command thresholds the volume at 185 and then extracts the largest connected component, which produces a binary WM segmentation. Second, run
mri_extract_largest_CC -T 90 EM_combined.mgz brain_mask_tmp.mgz
This command thresholds out the background (and some CSF), and creates a binary mask for the brain tissues.
Step 5. Convert relevant volumes for later surface deformation step into COR format for the convenience of later processing.
For example:
mri_convert -c wm_init_tmp.mgz wm_init.mgz mri_convert -c brain_mask_tmp.mgz brain_mask.mgz mri_convert -c mef30_avg_inu.mgz mef30_avg_inu_cor.mgz mri_convert -c mef5_avg_inu.mgz mef5_avg_inu_cor.mgz
This completes the automatic processing for Stage 1. The major output is the WM segmentation. Unfortunately, manual editing is typically required to refine the automatic WM segmentation result. This is the reason why the automatic result is named as "wm_init.mgz". Manual editing should take care of the following issues:
- Isolate cerebellum from the cortical WM.
- Edit WM to correct segmentation errors caused by artifacts in the scans due to the fixing of the brain. For example, some WM is often misclassified as CSF due to the darkening in the T1-weighted scans.
Edit WM to correct topology defects. Sometimes the WM segmentation has large handles. It may be benefical to finish editing 1 & 2 first, then run the surface tessellation and inflation steps as described in the next Stage, and then use the inflated surface to guide the WM topology editing. Of cause, the surface tessellation need be updated once the WM is edited.
The manual edited WM volume should be saved in a file named "wm.mgz".
Note that all that is needed before the next stage can be successfully run is a good WM segmentation. Thus, if an alternative approach can be used to get better WM segmentation, the EM segmentation program won't be needed.
Stage 2: Cortical Surface Reconstruction
It involves the following steps.
Step 1. Extract the largest connected components from the final WM volume
If manual editing was performed in Stage 1 to separate the cerebellum from the cortical WM, this step will discard the cerebellum and keep only the cortical WM as desired:
mri_extract_largest_CC -T 127 -hemi $hemi wm.mgz filled.mgz
The output volume is directly named "filled.mgz" because an ExVivo scan typically contains only one hemisphere. If this is not true, the standard "mri_fill" may be used to computed the filled volume from the WM volume. The option "-hemi" tells the program which hemisphere is contained in this ExVivo data. The program will then correctly set the intensity value in the output "filled.mgz" volume. For example, if "-hemi rh" is used, the "filled.mgz" volume will have only two intensity values, 127 for foreground region (WM) and 0 for background. If "-hemi lh" is used, the foreground region will have a value of 255.
Step 2. Run the surface tesselation, inflation, topology correction steps as in standard recon-all to get a topology correct orig surface
Step 3. This is the final step of the ExVivo cortical surface reconstruction
The orig surface from the previous step is deformed in this step to get the white and pial surfaces. The program used is called "mris_exvivo_surfaces". It is again modified from the original "mris_make_surfaces" program and tailored towards the ExVivo data we typically have. This program takes in two input volumes instead of a single T1 volume as in mris_make_surfaces. The two input volumes correspond to the flash30 and the flash5 average volume respectively. Using the two input channels gives better surface reconstruction than using a single T1-weighted image. The overall principle is very similar to that of the mris_make_surfaces. The major difference is that the intensity term for the deformable surface model now involves two channels, and the image gradient computation also takes into account of both input images. The command line is as follows:
mris_exvivo_surfaces -mgz \ -flash30 mef30_avg_inu_cor \ -flash5 mef5_avg_inu_cor \ $s $hemi
Further improvement of the algorithm will involve allowing more than two input channels and better design of multi-spectral image gradient computation.
The following information is intended for users inside NMR center:
- Example scripts for the above procedure can be found under /space/dijon/32/users/xhan/ExVivo/scripts/run_EM_one.csh and build_surface_one.csh.
For help using the above scripts, please consult Allison Stevens. See also ExvivoRecons.
Optional Step 4. Generating Surface Models of Layer IV
Here is an example of an ex vivo MRI
freeview/autofs/space/rock_004/users/hires/I12/I12_whole_hemi_200um_stria_4_7_2005/m ri/flash20_200um_avg.norm.mgz
To generate surface models can use this script:
/homes/4/fischl/dev/mris_deform/run_layerIV.csh