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.

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 surface topology correction and surface deformation. 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. Also, the COR format forces the intensity values to be of 8-bit BYTE type, which drops information from the original float-typed image. An example for generating the flash30 average volume is as follows:

Step 2. Run intensity inhomogeneity correction using the "nu_correct" program. Although the EM segmentation method can perform simultaneous INU correction, correcting INU beforehand is faster, and the nu_correct works reasonably well for the MEF data. Two passes of nu_correct is sufficient in general. A script called "inu_correct_vol.csh" can be used, which is a wrapper for the nu_correct and exactly performs two passes of nu_correction:

The same procedure is repeated for all the other flip-angles as well.

Step 3. EM segmentation using mri_ms_EM:

The option "-no_INU" indicates no INU correction will be performed inside the 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) to define the brain mask, only voxels within the brain region (with intensity value > 200) will be considered in the segmentation. The threshold need be set manually. "-M 4" indicates that the number of classes/labels 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 takes into account remaining "background" voxels, and offers better segmentation in practice. "-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 be names as $outdir/EM_sth.

For the above command, the major output is a synthesized gray-level byte volume saved in "$outdir/EM_combined.mgz". This volume encodes the fuzzy EM segmentation results. 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. 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 region. First, run

This command thresholds the volume at 185 and then extracts the largest connected component, which produces a binary WM segmentation. Second, run

This command thresholds out the background and the CSF, and creates a binary mask for the brain tissues.

Step 5. Convert relevant volumes for later surface deformation step into COR format for convenience of later processing. For example:

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 that the automatic result is named as "wm_init.mgz". Manual editing should take care of the following issues:

1. Isolate cerebellum from the cortical WM.

2. Edit WM to elliminate segmentation errors caused by artifacts in the scans due to the fixing of the brain.

3. Edit WM to correct topology defects. Thus, it may be benefical to finish 1& 2 first, then run the surface tessellation and inflation steps in the next Stage, and using the inflated surface to guild 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 in the next stage is a good WM segmentation. Thus, if a good WM segmentation is available, the EM segmentation 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 manua editing was performed in Stage 1 to separate the cerebellum from the cortical WM, this step will remove the disconnected cerebellum and keep only the cortical WM as desired:

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.

Step 2. Run the surface tesselation, inflation, topology correction steps as in 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 works as follows:

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:

1. 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.

2. For help using the above scripts, please consult Niranjini Rajendran