Differences between revisions 19 and 20

Longitudinal analysis

The longitudinal stream of FreeSurfer is currently overhauled. Please stay tuned ...

This page will describe the new version of the Longitudinal stream, while LongitudinalProcessingOld describes the old version. It is advised not to use the old version due to its bias with respect to time point 1.

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 sets (cross-sectionally), and thus not optimal for the processing of longitudinal data series. It is an active research area at the Martinos Center for Biomedical Imaging, how to obtain robust and more reliable cortical and subcortical morphological measurements by incorporating additional (temporal) information in a longitudinal data series.

The longitudinal scheme is designed to be unbiased wrt. any time point. Instead of initializing it with information from a specific time point, a base/template volume is created and run through FreeSurfer. This template can be seen as an initial guess for the segmentation and surface reconstruction. 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 the unbiased template, 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. Aditionally to the base template new probabilistic methods (temporal fusion) were introduced to further reduce the variablility of results across timepoints. For these algorithms it becomes necessary to process all timepoints cross sectionally first.

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 all time points with the default workflow (tpN is one of the timepoints):

recon-all -all -s <tpNsubjid> -i path_to_tpN_dcm

2) create the base template and process it cross sectionally:

  recon-all -base <baseid> -base-insubj <tp1subjid> -base-insubj <tp2subjid> ... -all

3) longitudinally process all timepoints:

  recon-all -long <tpNsubjid> <baseid> -all

4) do comparisons on results from 3), e.g. calculate differences between <tp1subjid>.long.<baseid> and <tp2subjid>.long.<baseid>

The longitudinal processing stream always produces output subject data containing <tpNsubjid>.long.<baseid> in the name (to help distinguish from the default stream).

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. Both have already been processed cross sectionally using the standard FreeSurfer recon-all. First we discuss the construction of the base template and then the longitudinal processing of tpN (step 3).

1. Creation of Base

In step 2 the base template is created. This is done by mri_robust_template, which constructs a mean/median T1.mgz volume together with the transforms that align each tp's T1.mgz volume with the template. Thus, all timepoints are brought to an unbiased common space. This template volume is then processed cross sectionally with the standart FreeSurfer stream. The only difference is that the norm.mgz is not created in the usuall way, but again computed as the mean/median of the norm.mgz of all time points. The corresponding transforms are more accurate than the transforms obtained from the T1 image. Once the base is processed it can be used to initialize many steps in the longitudinal stream (step 3).

2. Difference 1: Register tpN with base

The longitudinal scheme requires aligning the image data of tpN to base, but this map has allready been construced in a robust manner during the base creation (see also mri_robust_register which is used by mri_robust_template to construct the maps). Since all data sets come from the same subject, these rigid registrations with 7 dof (translation,rotation,intensity scaling) are enough to get a good alignment between the images. The registration result (a coordinate transformation matrix) and its inverse have been saved during the template construction and will be used to transfer information between time points in later steps.

...

More later

Original Author: MartinReuter

-Deletions are marked like this.
+Additions are marked like this.
 Line 3:
-''Original Author: Xiao Han (see also: XiaoNotes)''
+The longitudinal stream of FreeSurfer is currently overhauled. Please stay tuned ...



This page will describe the new version of the Longitudinal stream, while [[LongitudinalProcessingOld]] describes the old version. It is advised not to use the old version due to its bias with respect to time point 1.
-Line 5:
+Line 8:
-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.
-Line 7:
+Line 9:
-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.
+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 sets (cross-sectionally), and thus not optimal for the processing of longitudinal data series. It is an active research area at the Martinos Center for Biomedical Imaging, how to obtain robust and more reliable cortical and subcortical morphological measurements by incorporating additional (temporal) information  in a longitudinal data series.
-Line 9:
+Line 11:
-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 scheme is designed to be unbiased wrt. any time point. Instead of initializing it with information from a specific time point, a base/template volume is created and run through FreeSurfer. This template can be seen as an initial guess for the segmentation and surface reconstruction. 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 the unbiased template, 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. Aditionally to the base template new probabilistic methods (temporal fusion) were introduced to further reduce the variablility of results across timepoints. For these algorithms it becomes necessary to process all timepoints cross sectionally first.
-Line 15:
+Line 17:
-) cross-sectionally process tp1 subject (the default workflow):
+) cross-sectionally process all time points with the default workflow (tpN is one of the timepoints):
-Line 17:
+Line 19:
-recon-all -all -s <tp1subjid> -i path_to_tp1_dcm
+recon-all -all -s <tpNsubjid> -i path_to_tpN_dcm
-Line 20:
+Line 22:
-) longitudinally process tp1:
+) create the base template and process it cross sectionally:
-Line 22:
+Line 24:
-  recon-all -all -long <tp1subjid> <tp1subjid> -i path_to_tp1_dcm
+  recon-all -base <baseid> -base-insubj <tp1subjid> -base-insubj <tp2subjid> ... -all
-Line 25:
+Line 27:
-) longitudinally process tp2:
+) longitudinally process all timepoints:
-Line 27:
+Line 29:
-  recon-all -all -long <tp1subjid> <tp2subjid> -i path_to_tp2_dcm
+  recon-all -long <tpNsubjid> <baseid> -all
-Line 30:
+Line 32:
-) do comparisons on results from 2) and 3), e.g. calculate differences 

between <tp2subjid>.long.<tp1subjid> and <tp1subjid>.long.<tp1subjid>
+) do comparisons on results from 3), e.g. calculate differences 

between <tp1subjid>.long.<baseid> and <tp2subjid>.long.<baseid>
-Line 33:
+Line 35:
-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>.
+The longitudinal processing stream always produces output subject data containing  <tpNsubjid>.long.<baseid>  in the name (to help distinguish from the default stream).
-Line 35:
+Line 37:
-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).
+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. Both have already been processed cross sectionally using the standard FreeSurfer recon-all. First we discuss the construction of the base template and then the longitudinal processing of tpN (step 3).
-Line 37:
+Line 39:
-== 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

}}}
+== Creation of Base ==
-Line 48:
+Line 41:
-"$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.
+In step 2 the base template is created. This is done by [[mri_robust_template]], which constructs a mean/median T1.mgz volume together with the transforms that align each tp's T1.mgz volume with the template. Thus, all timepoints are brought to an unbiased common space. This template volume is then processed cross sectionally with the standart FreeSurfer stream. The only difference is that the norm.mgz is not created in the usuall way, but again computed as the mean/median of the norm.mgz of all time points. The corresponding transforms are more accurate than the transforms obtained from the T1 image. Once the base is processed it can be used to initialize many steps in the longitudinal stream (step 3).
-Line 50:
+Line 43:
-== 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

}}}
+== Difference 1: Register tpN with base ==

The longitudinal scheme requires aligning the image data of tpN to base, but this map has allready been construced in a robust manner during the base creation (see also [[mri_robust_register]] which is used by [[mri_robust_template]] to construct the maps). Since all data sets come from the same subject, these rigid registrations with 7 dof (translation,rotation,intensity scaling) are enough to get a good alignment between the images. The registration result (a coordinate transformation matrix) and its inverse have been saved during the template construction and will be used to transfer information between time points in later steps.
-Line 61:
+Line 46:
-"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.
+...
-Line 63:
+Line 48:
-== 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

}}}
+More later
-Line 73:
+Line 50:
-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

}}}
-Line 83:
+Line 51:
-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}.orig

}}}



and

{{{

mris_transform --dst $subjdir/mri/orig.mgz \

  $tp1dir/surf/${hemi}.white \

  $tp1dir/mri/transforms/$regfile \

  $subjdir/surf/${hemi}.orig_white

}}}



That 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 \

  $annot

}}}



This 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 $\

  annot

}}}



The 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.
+----

Original Author: MartinReuter