PETSurfer provides a set of tools within FreeSurfer for Partical Volume Correction (PVC) and Kinetic Modeling (MRTM1, MRTM2). While these are typically used for PET analysis, the tools can be used in any context where PVC is needed. PVC methods include the Symmetric Geometric Transfer Matrix (SGTM), two-compartment model (also known as the Meltzer Method), three-compartment model (also known as the Muller-Gartner (MG) Method), region-based voxel-wise (RBV). SGTM is used for ROI analysis where as the others are used for voxel-wise analysis. All PVC implementations also account for the volume fraction effect (VFE). The voxel-wise output can then be analyzed on the cortical surface or in the volume. PETSurfer will be released with FreeSurfer version 6.

Please cite these papers when preparing manuscripts:

Greve, D. N., Salat, D. H., Bowen, S. L., Izquierdo-Garcia, D., Schultz, A. P., Catana, C., ... & Johnson, K. A. (2016). Different partial volume correction methods lead to different conclusions: An 18 F-FDG-PET study of aging. NeuroImage, 132, 334-343.

Greve, D. N., Svarer, C., Fisher, P. M., Feng, L., Hansen, A. E., Baare, W., ... & Knudsen, G. M. (2014). Cortical surface-based analysis reduces bias and variance in kinetic modeling of brain PET data. Neuroimage, 92, 225-236.

In all cases, you will need a T1-weighted MRI of your subject of sufficient quality to run in FreeSurfer. FreeSurfer analysis must be done first. After that, follow the steps below.

1. Create a segmentation for the GTM

gtmseg --s subject

where "subject" is the name of the FreeSurfer subject when you ran recon-all. This creates a high-resolution segmentation (gtmseg.mgz) in the FS folder used to run the PVC methods. This should take about an hour or two. gtmseg.mgz will use aseg.mgz for subcortical structures, ?h.aparc.annot for cortical structures, and will estimate some extra-cerebral structures. There are ways to customize this segmentation to use different ROI definitions (eg, aparc.a2009s instead of aparc).

2. Register your PET image with the anatomical:

mri_coreg --s subject --mov template.nii.gz --reg template.reg.lta

where template.nii.gz is the template image for your PET data. If your PET data only has one frame (eg, an SUV image), then that will be your template. If your PET data has multiple frames (ie, dynamic), then you will need to create the template from the dynamic data. This can be done by extracting a single frame (mri_convert pet.nii.gz --frame frameno template.nii.gz) or averaging all the time frames together (eg, mri_concat pet.nii.gz --mean --o template.nii.gz). It might make sense to do a time-weighted average rather than simple average, but we do not have a tool to do that easily, but you can do it in matlab. For parallel operation, add --threads N where N is the number of CPUs you want to use. You should check the registration with:

tkregisterfv --mov template.nii.gz --reg template.reg.lta --surfs

If you are not using PVC, you can use the template.reg.lta to sample the PET volume onto the surface using mri_vol2surf, then apply standard surface-based analysis.

3. Apply Partial Volume Correction

mri_gtmpvc --i pet.nii.gz --reg template.reg.lta --psf FWHM --seg gtmseg.mgz 
 --default-seg-merge  --auto-mask PSF .01 --mgx .01 --o gtmpvc.output 

--psf FWHM is the full-width/half-max of the the point-spread function (PSF) of the scanner as measured in image space (also known as the burring function). The blurring function depends on the scanner and reconstruction method; here it is modeled as an isotropic Gaussian filter of the given FWHM. Eg, an HR+ is typically about 6mm. This parameter has nothing to do with applying smoothing. Set this to 0 to turn off PVC. If you don't know what to set this to, ask your PET physicist.
--seg gtmseg.mgz is the segmentation created with gtmseg
--default-seg-merge will merge several segmentations, eg, all the ventricular CSF segments will be merged into one ROI
--auto-mask FWHM .01 will create a mask to exclude voxels from the analysis (saves time). Output volumes will be reduced to a bounding box around this mask (saves space)
--mgx .01 Run Muller-Gartner analysis. 01 is the GM threshold. Only necessary if you want to do a voxel-wise analysis.
--o gtmpvc.output This is the output folder.

There will be many files in the output folder some of which are described here:

gtm.stats.dat is an easy to read text file where each row is an ROI, something like:

9 17 Left-Hippocampus subcort_gm 473 174.083 1.406 0.1216

9 = ninth row 17 = index for ROI Left-Hippocampus = name of ROI subcort_gm = tissue class 473 = number of PET voxels in the ROI 174 = variance reduction factor for ROI (based on GLM/SGTM) 1.406 = PVC uptake of ROI relative to Pons 0.1216 = resdiual varaince across voxels in the ROI

By default, this will scale by the intensity in pons. If you do not want scaling (eg, when doing a dynamic analysis), add --no-rescale

gtm.nii.gz is a nifti file with each "voxel" being an ROI. The value is the PVC uptake of ROI relative to Pons. For single time point data, this will be totally redundant with gtm.stats.dat. Use of the NIFTI format makes it easy to concatenate (mri_concat) these files together for group analysis (mri_glmfit).

nopvc.nii.gz - same interpretation as gtm.nii.gz except that the values have not been PVCed in any way.

mgx.{ctxgm,subctxgm,gm} - these are the voxel-wise values corrected using the "extended" Muller-Gartner method. ctxgm is cortical GM, subctxgm is the subcortical GM, and gm is all GM. These volumes will be of the same resolution as the input PET but the field of view will be reduced to that of a bounding box around the mask.

aux - this subfolder has auxiliary data that are often not of immediate use to the user. The exception is bbpet2anat.lta. This is a registration file that can be used to map the output PET volume (in the mask bounding box) to the anatomical space. This file should be used when mapping the volume to the surface, etc.

If you are doing dynamic (kinetic) analysis, thenadd the following arguments when running mri_gtmpvc:

--km-ref 8 47 --km-hb 11 12 13 50 51 52 --no-rescale

--km-ref 8 47 specifies the ROIs to use as the reference region for the MRTM models. 8 and 47 are the cerebellar hemisphers as found in $SUBJECTS_DIR/subject/mri/gtmseg.ctab (see also $FREESURFER_HOME/FreeSurferColorLUT.txt). This creates the file km.ref.tac.dat with the reference value for each time point.

--km-hb 11 12 13 50 51 52 specifies the ROIs to use as the high-binding region if using MRTM2. This creates km.hb.tac.nii.gz with the value for the high-binding region for each time point.

4. Kinetic Modeling (KM). KM is done with either MRTM1 or MRTM2. To run MRTM1:

mri_glmfit --y km.hb.tac.nii.gz --mrtm1 km.ref.tac.dat time.dat --o mrtm1 --no-est-fwhm --nii.gz

where time.dat is a simple ASCII text file with the acquisition time (in seconds) of each time point in the time activity curve (TAC). This will create a folder called mrtm1 in which will be a file called k2prime.dat. The value of k2prime will be used in the MRTM2 analysis below

For the MRTM2 analysis

set k2p = `cat mrtm1/k2prime.dat`
mri_glmfit --y gtm.nii.gz --mrtm2 km.ref.tac.dat time.dat $k2p --o mrtm2 --no-est-fwhm --nii.gz

This will create a folder called mrtm2 in which the bp.nii.gz will be the non-displaceable binding potential (BPND) for each ROI.

5. Surface-based analysis - details of surface-based analysis are available in other locations in the wiki, but here is bascially what you do. Below, the mgx.ctxgm.nii.gz is used, but the commands will apply to any of the volumes in the gtm output folder.

Sample the mgx volume onto the left hemisphere surface

mri_vol2surf --mov mgx.ctxgm.nii.gz --reg aux/bbpet2anat.lta --hemi lh
--projfrac 0.5 --o lh.mgx.ctxgm.fsaverage.sm00.nii.gz --cortex 
--trgsubject fsaverage

This command samples it onto the surface of fsaverage via the individual subject's surface. --cortex says to mask out non-cortical regions. --projfrac 0.5 says to sample halfway between the white and pial surfaces. If you want to average over the cortical ribbon, you can use --projfrac-avg .2 .8 .1, which says to start 20% into the ribbon, sample every 10%, and stop at 80% of the thickness.

If you are not doing kinetic modeling, concatenate all your subjects together into a stack file

mri_concat subj1/lh.mgx.ctxgm.fsaverage.sm00.nii.gz 
--o all.lh.mgx.ctxgm.fsaverage.sm00.nii.gz --prune

Smooth on the surface

mris_fwhm --smooth-only --i all.lh.mgx.ctxgm.fsaverage.sm00.nii.gz 
--fwhm 5 --o all.lh.mgx.ctxgm.fsaverage.sm05.nii.gz --cortex 
--s fsaverage --hemi lh

Run group analysis GLM

mri_glmfit --y all.lh.mgx.ctxgm.fsaverage.sm05.nii.gz 
--surface fsaverage lh --fsgd your.fsgd ...

Note: when running correction for multiple comparisons (mri_glmfit-sim), it is highly recommended that permutation be used.

For the kinetic modeling MRTM2 surface-based analysis, lh.mgx.ctxgm.fsaverage.sm00.nii.gz is the surface-based TAC for the individual sampled onto fsaverage. Smooth (mris_fwhm) this file rather that concatenating across subject. Then use the smoothed file as input to the MRTM2 analysis:

set k2p = `cat mrtm1/k2prime.dat`
mri_glmfit --y lh.mgx.ctxgm.fsaverage.sm05.nii.gz  --mrtm2
km.ref.tac.dat time.dat $k2p --o mrtm2.lh.sm05 --no-est-fwhm --nii.gz 
--surface fsaverage lh

This will create mrtm2.lh.sm05/bp.nii.gz. These can be concatenated across subject, eg,

mri_concat subj1/mrtm2.lh.sm05/bp.nii.gz
--o all.mrtm2.lh.sm05.bp.nii.gz --prune

Do not perform another round of smoothing on this output. If you want more smoothing, redo the smoothing of the TACs.

Then run group analysis with mri_glmfit with all.mrtm2.lh.sm05.bp.nii.gz as input.

mri_glmfit --y all.mrtm2.lh.sm05.bp.nii.gz --surface fsaverage lh --fsgd your.fsgd

6. Subcortical volume-based analysis.

Sample the mgx.subctxgm.nii.gz volume into the 2mm space of the MNI305:

mri_vol2vol --mov subctxgm.nii.gz --reg aux/bbpet2anat.lta
--tal --talres 2  --o subctxgm.mni305.2mm.sm00.nii.gz

Smooth in 3D. This command shows 5mm, but it could be something else

mri_fwhm --smooth-only --i subctxgm.mni305.2mm.sm00.nii.gz --fwhm 5 \
--o subctxgm.mni305.2mm.sm05.nii.gz \
--mask $FREESURFER_HOME/subjects/fsaverage/mri.2mm/subcort.mask.mgz

If not doing kinetic modeling, then the subctxgm.mni305.2mm.sm05.nii.gz can be mri_concat together.

If doing KM, then

set k2p = `cat mrtm1/k2prime.dat`
mri_glmfit --y subctxgm.mni305.2mm.sm05.nii.gz  --mrtm2
km.ref.tac.dat time.dat $k2p --o mrtm2.subctxgm.mni305.2mm.sm05 --no-est-fwhm --nii.gz 

This will create mrtm2.mni305.2mm.sm05/bp.nii.gz. These can be concatenated across subject, eg,

mri_concat subj1/mrtm2.mni305.2mm.sm05/bp.nii.gz 
--o all.mrtm2.mni305.2mm.sm05.bp.nii.gz  --prune

Do not perform another round of smoothing on this output. If you want more smoothing, redo the smoothing of the TACs.

Then run group analysis with mri_glmfit with all.mrtm2.mni305.2mm.sm05.bp.nii.gz as input.

mri_glmfit --y all.mrtm2.mni305.2mm.sm05.bp.nii.gz  --fsgd your.fsgd ...

Note: when running correction for multiple comparisons (mri_glmfit-sim), it is highly recommended that permutation be used.