Longitudinal analysis (pre FS 4.4)
For the current longitudinal stream see LongitudinalProcessing ! This stream (pre FS 4.4) is biased toward the time point selected as base (usually TP1).
Original Author: Xiao Han (see also: XiaoNotes)
Compared with cross-sectional studies, a longitudinal design can significantly reduce the confounding effect of inter-individual morphological variability by using each subject as his or her own control. As a result, longitudinal imaging studies are getting increased interest and popularity in various aspects of neuroscience.
The default FreeSurfer pipeline is designed for the processing of individual data set, and thus not optimal for the processing of longitudinal data series. How to take into account the subject-wise and temporal correlation in a longitudinal data series and get more robust and reliable cortical and subcortical morphological measurements is an active research area in NMR center.
Although further improvements are constantly being sought, as the result of an early attempt, a longitudinal processing pipeline has been designed for the processing of longitudinal data and has shown success. The basic idea of this longitudinal scheme is to borrow the processing results of an early time point to help the processing of a later time point. As easily noticed, the FreeSurfer cortical and subcortical segmentation and parcellation procedure involves solving many complex nonlinear optimization problems, such as the deformable surface reconstruction, the nonlinear atlas-image registration, and the nonlinear spherical surface registration. These nonlinear optimization problems are usually solved using iterative methods, and the final results are known to be highly sensitive to the selection of a particular starting point (a.k.a. algorithm initialization). It's our belief that by initializing the processing of a new data set in a longitudinal series using the processed results from another time point, we can reduce the random variation in the processing procedure and improve the robustness and sensitivity of the overall longitudinal analysis. Such an initialization scheme makes sense also because a longitudinal design is often targeted at detecting small or subtle changes.
The longitudinal processing scheme is coded in the recon-all script via the "-long" flag. Use the -help flag for help on its options.
Here is a summary of the longitudinal workflow:
1) cross-sectionally process tp1 subject (the default workflow):
recon-all -all -s <tp1subjid> -i path_to_tp1_dcm
2) longitudinally process tp1:
recon-all -all -long <tp1subjid> <tp1subjid> -i path_to_tp1_dcm
3) longitudinally process tp2:
recon-all -all -long <tp1subjid> <tp2subjid> -i path_to_tp2_dcm
4) do comparisons on results from 2) and 3), e.g. calculate differences between <tp2subjid>.long.<tp1subjid> and <tp1subjid>.long.<tp1subjid>
The longitudinal processing stream always produces output subject data containing <tp2subjid>.long.<tp1subjid> in the name (to help distinguish from the default stream). Note that it is essential to specify the -i option with the path to the DICOM or MGZ image (or several -i options to more than one image). This is because tp2 is usually not processed, so step 3) can be run directly on the original tp2 images. The <tp2subjid> is only specified to create the subject directory <tp2subjid>.long.<tp1subjid>.
In the following, we list the major differences between the two processing streams. We assume that the longitudinal series have two time-points, tp1 and tp2, and tp1 has already been processed using the standard FreeSurfer recon-all. To avoid confusion we will discuss the longitudinal processing of tp2 (step 3).
Difference 1: Register tp2 with tp1
The longitudinal scheme requires aligning the image data of tp2 to tp1. Since both data sets come from the same subject, a linear registration with a dof equal 6 or 9 is enough to get a good alignment between them. The alignment/registration can be computed using the FSL FLIRT tool, and the registration result (a coordinate transformation matrix) is saved and will be used to transfer information from tp1 to tp2 in later steps. As can be seen in the recon-all script, the alignment is computed using the following command:
fsl_rigid_register -cost corratio \
-i $tp1dir/mri/orig.mgz \
-r $subjdir/mri/orig.mgz \
-o $subjdir/mri/orig_${tp1id}_to_${subjid}.mgz \
-dof 6 \
-ltamat $tp1dir/mri/transforms/$regfile"$regfile" stores the registration as a LTA file. The filename is automatically set inside the script, unless the user changes it using the "-reg" option.
Difference 2: Talairach Registration
In the default stream, the Talairach registration is computed using MINC tools. In the longitudinal scheme, the Talairach registration for tp2 is computed by simply composing the linear registration from tp1 to tp2 and the Talairach registration for tp1. In particular, the following command is used:
mri_concatenate_lta -invert1 -tal \ $tp1dir/mri/orig.mgz \ $talairach_template_file \ $tp1dir/mri/transforms/$regfile \ $tp1dir/mri/transforms/talairach.xfm \ $subjdir/mri/transforms/talairach.xfm
"mri_concatenate_lta" is a program for combining two consecutive linear transformation/registration into a single total transformation. "-invert1" indicates that the first transformation matrix need be inverted before composing with the second one. "-tal src_file trg_file" gives the source and target volume names for the first Talairach registration. "$talairach_template_file" is set to "${MINC_BIN_DIR}/../share/mni_autoreg/average_305.mnc" by default.
Difference 3: Intensity Normalization
If "-usetp1ctrlvol" is set, the longitudinal scheme will first map the automatically computed control points for tp1 (computed inside mri_normalize) to the data space of tp2, and then use them to perform intensity normalization for tp2. By default, the control points for tp1 are not saved in the original script. Thus, recon-all -long will first regenerate the control points using
mri_normalize -mask $tp1dir/mri/brain.mgz \ -W $tp1dir/mri/ctrl_vol.mgz \ $tp1dir/mri/bias_vol.mgz \ $tp1dir/mri/nu.mgz \ $tp1dir/mri/T1_tmp.mgz
Then the intensity normalization using tp1's control points (saved as a volume file "ctrl_vol.mgz") is performed as follows:
mri_normalize_tp2 -T1 $tp1dir/mri/nu.mgz \ -ctrl $tp1dir/mri/ctrl_vol.mgz \ -mask1 $tp1dir/mri/brain.mgz \ -xform $tp1dir/mri/transforms/$regfile \ $subjdir/mri/nu.mgz \ $subjdir/mri/T1.mgz
We haven't investigated the mapping of manually picked control points from tp1 to tp2. If manual control points are necessary, it's better that the user pick them manually for the tp2 volume too.
Difference 4: Skull-stripping
In the longitudinal scheme, the brain-mask for tp2 is borrowed from tp1 (with the tp1-to-tp2 registration transformation properly applied):
mri_mask -xform $tp1dir/mri/transforms/$regfile T1.mgz $tp1dir/mri/$BM $BMA
Difference 5: GCA Registration
Similar to the Talairach registration, in the longitudinal scheme, the linear subcortical atlas registration for tp2 is computed by composing the linear atlas registration for tp1 and the linear mapping from tp1 to tp2:
mri_concatenate_lta -invert1 $tp1dir/mri/transforms/$regfile \ $tp1dir/mri/transforms/talairach.lta \ transforms/talairach.lta
"-invert1" again indicates that the first LTA need be inverted before composing with the second. The rationale is easily understandable by drawing the source and target space for each LTA.
Difference 6: Canonical Registration (nonlinear)
In the standard recon-all stream, this nonlinear atlas registration is initialized by the previous linear registration result:
mri_ca_register $UseCAAlign -mask brainmask.mgz \ -T transforms/talairach.lta \ norm.mgz \ $FREESURFER_HOME/average/$GCA \ transforms/talairach.m3z
In the longitudinal stream, the nonlinear registration for tp2 is initialized by the nonlinear registration result of tp1, and only local refinements are further performed to get the final registration:
mri_ca_register $UseCAAlign -mask brainmask.mgz \ -levels 2 -A 1 -l $tp1dir/mri/transforms/talairach.m3z \ $tp1dir/mri/transforms/$regfile \ -T transforms/talairach.lta \ norm.mgz \ $FREESURFER_HOME/average/$GCA \ transforms/talairach.m3z
The options "-levels 2 -A 1" tell the program to do only local refinement. "-l $tp1dir/mri/transforms/talairach.m3z $tp1dir/mri/transforms/$regfile" is the longitudinal option that specifies the initial nonlinear registration and the mapping from tp1 to tp2. "-T transforms/talairach.lta" gives the linear registration as in the original scheme. This linear registration will still be used when cross-sequence atlas renormalization is needed (as specified by the $UseCAAlign option).
Difference 7: Subcortical Segmentation
In the longitudinal processing scheme, the final subcortical segmentation step again uses the segmentation result from tp1 to initialize the iterative labeling process:
mri_ca_label $UseCAAlign -l $tp1dir/mri/aseg.mgz \ $tp1dir/mri/transforms/$regfile \ norm.mgz \ transforms/talairach.m3z \ $FREESURFER_HOME/average/$GCA aseg.auto.mgz
Difference 8: WM surface tessellation and topology correction
The whole process for generating the orig surface that involves WM tessellation, topology correction, and etc is no longer needed in the longitudinal processing of tp2, since the final gray/white surface from tp1 will be used as the initial gray/white surface for tp2. Consequently, the original process is replaced by the following two commands just before running mris_make_surfaces:
mris_transform --dst $subjdir/mri/orig.mgz \
$tp1dir/surf/${hemi}.white \
$tp1dir/mri/transforms/$regfile \
$subjdir/surf/${hemi}.origand
mris_transform --dst $subjdir/mri/orig.mgz \
$tp1dir/surf/${hemi}.white \
$tp1dir/mri/transforms/$regfile \
$subjdir/surf/${hemi}.orig_whiteThat is, the final white surface from tp1 is mapped to the tp2 volume and is used as the orig surface for tp2 and also used as the "orig_white", i.e., the "initial white surface", for the later surface deformation. Similarly, the final pial surface from tp1 is mapped to tp2 and will be used as the initial pial surface for tp2:
mris_transform --dst $subjdir/mri/orig.mgz \
$tp1dir/surf/${hemi}.pial \
$tp1dir/mri/transforms/$regfile \
$subjdir/surf/${hemi}.orig_pial
Difference 9: Make Final Surfaces
In the standard FreeSurfer, the orig surface is deformed to get the white surface, and the white surface is further deformed to get the pial surface. In the longitudinal scheme, the white and pial surfaces are both initialized using tp1's surface deformation results as described above, and the command line becomes:
mris_make_surfaces -orig_white orig_white \ -orig_pial orig_pial \ -orig orig_white \ -long -max 3.5 -mgz -T1 brain.finalsurfs \ $subjid $hemi
The options "orig_white" and "orig_pial" specifies the initialization for the white and pial surface respectively. "-long" indicates the longitudinal scheme. Note that in the original scheme, the white and pial surfaces are guaranteed not to cross each other. This non-cross-intersection between the final white and pial surfaces won't be guaranteed in the new scheme since the initial pial surface is independent of the final white surface. If non-cross-intersection is important, the user can remove the "-orig_pial" option and its argument, such that the pial surface is still initialized by the final white surface. But using both "-orig_white" and "-orig_pial" gives better thickness reproducibility.
Difference 10: Spherical Mapping
In the longitudinal scheme, the spherical map ($hemi.sphere) for tp2 is simply copied from tp1. It's not hard to see why this is reasonable.
Difference 11: Spherical Registration
Like in subcortical segmentation, the nonlinear spherical registration for tp2 in the longitudinal scheme is initialized by the sphere.reg from tp1:
mris_register -curv -nosulc -norot \ $tp1dir/surf/$hemi.sphere.reg \ $AvgTif \ ../surf/$hemi.sphere.reg
"-nosulc" and "-norot" tell the program to ignore the initial linear alignment step as in standard spherical registration.
(NOTE: the same principle may be applied for the contra-surface registration for tp2, but it's not currently implemented)
Difference 12: Cortical Parcellation
Similar to the subcortical labeling step, the cortical parcellation for tp2 in the longitudinal scheme is again initialized with the parcellation result from tp1:
mris_ca_label -long -R $tp1dir/label/${hemi}.aparc.annot \
$subjid $hemi \
../surf/$hemi.sphere.reg \
$CPAtlas \
$annotThis is done also for the second parcellation using another surface labeling scheme:
mris_ca_label -long -R $tp1dir/label/${hemi}.aparc.a2005s.annot \
$subjid $hemi \
../surf/$hemi.sphere.reg \
$CPAtlas $\
annotThe above is a complete list of all the differences between the standard recon-all and the longitudinal "recon-all -long" processing stream. The new options added to "recon-all" include "-long", "-no-orig-pial", and "-regfile". Their meanings should be clear from the above description.