Image Realignment

Realignment is referring to in this case as within-subject registration of the 4D fMRI data.

library(spm12r)
####################################
# Realignment
####################################
realign_batch = build_spm12_realign( 
  filename = run_fmri, 
  register_to = "mean",
  reslice = "mean"
)

## # Getting Number of Time Points

## # Reslice is mean

print(names(realign_batch))

## [1] "spm"      "script"   "outfiles" "rp"       "mean"     "mat"

print(names(unlist(realign_batch$spm)))

##  [1] "spm.spatial.realign.estwrite.data"            
##  [2] "spm.spatial.realign.estwrite.eoptions.quality"
##  [3] "spm.spatial.realign.estwrite.eoptions.sep"    
##  [4] "spm.spatial.realign.estwrite.eoptions.fwhm"   
##  [5] "spm.spatial.realign.estwrite.eoptions.rtm"    
##  [6] "spm.spatial.realign.estwrite.eoptions.interp" 
##  [7] "spm.spatial.realign.estwrite.eoptions.wrap"   
##  [8] "spm.spatial.realign.estwrite.eoptions.weight" 
##  [9] "spm.spatial.realign.estwrite.roptions.which"  
## [10] "spm.spatial.realign.estwrite.roptions.interp" 
## [11] "spm.spatial.realign.estwrite.roptions.wrap"   
## [12] "spm.spatial.realign.estwrite.roptions.mask"   
## [13] "spm.spatial.realign.estwrite.roptions.prefix"

if (have_matlab()) {
  realigned = spm12_realign( 
    filename = run_fmri, 
    register_to = "mean",
    reslice = "mean",
    clean = FALSE,
    install_dir = install_dir
  )
  print(names(realigned))
}

## cloning into '/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T//RtmpuSoCbO/spm12'...
## Receiving objects:   1% (65/6496),   18 kb
## Receiving objects:  11% (715/6496),  259 kb
## Receiving objects:  21% (1365/6496), 3167 kb
## Receiving objects:  31% (2014/6496), 8349 kb
## Receiving objects:  41% (2664/6496), 10204 kb
## Receiving objects:  51% (3313/6496), 11596 kb
## Receiving objects:  61% (3963/6496), 32906 kb
## Receiving objects:  71% (4613/6496), 37689 kb
## Receiving objects:  81% (5262/6496), 59573 kb
## Receiving objects:  91% (5912/6496), 71860 kb
## Receiving objects: 100% (6496/6496), 129916 kb, done.

## # Getting Number of Time Points
## 
## # Reslice is mean

## # Running matlabbatch job

## # Creating matlabbatch job

## # Adding SPMDIR: /var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T//RtmpuSoCbO/spm12\n

## Running job: /var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T//RtmpuSoCbO/filebba0140cb031.m, which calls /var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T//RtmpuSoCbO/filebba03662dad4.m

## Command run is:

## matlab  -nodesktop -nosplash -nodisplay -r  "try, run('/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T//RtmpuSoCbO/filebba0140cb031.m'); catch err, disp(err.message); exit(1); end; exit(0);"

## # Result is 0

## [1] "spm"      "script"   "outfiles" "rp"       "mean"     "mat"     
## [7] "result"

Overall the spm12_realign does the realignment. There is some discussion of performing realignment before slice-timing correction because estimation of motion parameters may be skewed after slice-timing correction. We see that the output realigned has the output 4D fMRI data (outfiles), the realignment parameters (rp), voxel-wise mean after realignment (mean), and the matrix of transformations for the 4D series (mat).

####################################
# Read in Motion data
####################################
if (have_matlab()) {
  rpfile = realigned[['rp']]
  rp = read.table(file = rpfile, header = FALSE)
  colnames(rp) = c("x", "y", "z", 
    "roll", "pitch", "yaw")
  rp = as.matrix(rp)
  print(head(rp))
  print(dim(rp))
}

##                x            y             z          roll         pitch
## [1,]  0.00000000  0.000000000 -4.410886e-16 -2.507721e-37  2.507721e-37
## [2,] -0.03250006  0.022276041  6.105989e-02 -3.616417e-04 -1.389442e-04
## [3,]  0.02884634 -0.005311455  5.649025e-02 -3.774759e-04 -1.798493e-04
## [4,]  0.02379671 -0.056644719  2.574345e-02  5.180747e-05 -1.756995e-04
## [5,]  0.04293237 -0.056942376  2.414558e-02 -2.210120e-05 -8.367122e-05
## [6,]  0.02858868 -0.038461745  8.289177e-02 -2.785027e-04 -3.532936e-04
##                yaw
## [1,]  6.288666e-74
## [2,] -2.732142e-04
## [3,]  1.048090e-04
## [4,] -6.489734e-06
## [5,]  2.291027e-04
## [6,] -8.279382e-06
## [1] 200   6

Slice-Timing Correction

A slice-timing correction does interpolation since each slice was not actually taken at the same time point, but a shifted time point over the course of an entire TR. The correction requires you to input the reference slice (in this case the median, ref_slice), the repetition time (tr), time between the first and the last slice within one scan (ta), and the order the slices were acquired. In our case, it was done in an ascending, contiguous order, so we created the slice order as such. If you used descending or interleaved acquisition, then this must be changed accordingly.

####################################
# Slice Timing Correction
####################################
nslices = oro.nifti::nsli(run_fmri)
slice_order = 1:nslices
ref_slice = slice_order[median(seq(nslices))]
ta = tr - tr/nslices
n_time_points = ntim(run_fmri)

# using run_fmri because if you do not have
# matlab, then realigned not available
st_batch = build_spm12_slice_timing(
  filename = run_fmri,
  time_points = seq(n_time_points),
  nslices = nslices,
  tr = tr,
  ref_slice = ref_slice,
  prefix = "a")

print(names(st_batch))

## [1] "spm"                 "script"              "orig_filename"      
## [4] "base_name"           "temporary_directory" "outfile"

print(names(unlist(st_batch$spm)))

## [1] "spm.temporal.st.scans"    "spm.temporal.st.nslices" 
## [3] "spm.temporal.st.tr"       "spm.temporal.st.ta"      
## [5] "spm.temporal.st.so"       "spm.temporal.st.refslice"
## [7] "spm.temporal.st.prefix"

if (have_matlab()) {
  
  st_results = spm12_slice_timing(
    filename = realigned[['outfiles']],
    time_points = seq(n_time_points),
    nslices = nslices,
    tr = tr, 
    slice_order = slice_order,
    ta = ta, 
    ref_slice = ref_slice,
    prefix = "a", 
    clean = FALSE, 
    retimg = FALSE,
    install_dir = install_dir
    )
  aimg = st_results$outfile
  print(aimg)
  mean_img = realigned[["mean"]]
  mean_nifti = readnii(mean_img)
}

## # Running matlabbatch job

## # Creating matlabbatch job

## # Adding SPMDIR: /var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T//RtmpuSoCbO/spm12\n

## Running job: /var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T//RtmpuSoCbO/filebba055ccb9b1.m, which calls /var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T//RtmpuSoCbO/filebba07639f042.m

## Command run is:

## matlab  -nodesktop -nosplash -nodisplay -r  "try, run('/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T//RtmpuSoCbO/filebba055ccb9b1.m'); catch err, disp(err.message); exit(1); end; exit(0);"

## # Result is 0

## [1] "/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T//RtmpuSoCbO/afilebba053f1c6cb.nii"

We see the output aimg has the filename of the slice-timing corrected 4D image. We also read in the mean image into a nifti object (mean_nifti).

if (!have_matlab()) {
  run_fmri = filename_check(run_fmri)
}

Spatial Normalization

if (have_matlab()) {
  acpc_reorient(
    infiles = c(mean_img, aimg),
    modality = "T1",
    install_dir = install_dir
    )
}

## 
##  #Reorientation /private/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T/RtmpuSoCbO/meanfilebba053f1c6cb.nii

## Script created: /var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T//RtmpuSoCbO/filebba025a7fcdc.m

## Command run is:

## matlab  -nodesktop -nosplash -nodisplay -r  "try, run('/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T//RtmpuSoCbO/filebba025a7fcdc.m'); catch err, disp(err.message); exit(1); end; exit(0);"

##   /private/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T/RtmpuSoCbO/meanfilebba053f1c6cb.nii 
## "/private/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T/RtmpuSoCbO/meanfilebba053f1c6cb.nii" 
##             /var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T//RtmpuSoCbO/afilebba053f1c6cb.nii 
##    "/private/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T/RtmpuSoCbO/afilebba053f1c6cb.nii"

if (have_matlab()) {
  bbox = matrix(
        c(-90, -126, -72, 
        90, 90, 108), 
        nrow = 2, byrow = TRUE)
  print(bbox)
  direct_norm = spm12_normalize(
    filename = mean_img,
    other.files = c(mean_img, aimg),
    bounding_box = bbox,
    clean = FALSE,
    install_dir = install_dir
    )
  print(names(direct_norm))
  dnorm_files = direct_norm$outfiles
  dnorm_mean_img = readnii(dnorm_files[1])
}

##      [,1] [,2] [,3]
## [1,]  -90 -126  -72
## [2,]   90   90  108

## # Adding SPMDIR: /var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T//RtmpuSoCbO/spm12

## # Running script /var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T//RtmpuSoCbO/Executable.m
## which calls /var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T//RtmpuSoCbO/Normalize_Estimate_and_Write_job.m

## Command run is:

## matlab  -nodesktop -nosplash -singleCompThread -nodisplay -r  "try, run('/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T//RtmpuSoCbO/Executable.m'); catch err, disp(err.message); exit(1); end; exit(0);"

## # Result is 0

## [1] "spm"      "script"   "outfiles" "result"

if (have_matlab()) {
  anatomical = files["anatomical"]
  anat_img = checknii(anatomical)
  print(anat_img)
  acpc_reorient(
    infiles = anat_img,
    modality = "T1",
    install_dir = install_dir
    )

  coreg = spm12_coregister(
    fixed = mean_img,
    moving = anat_img,
    prefix = "r",
    install_dir = install_dir
    )
  
  coreg_anat = coreg$outfile
  coreg_img = readnii(coreg_anat)
  double_ortho(coreg_img, mean_nifti)
}

## [1] "/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T/RtmpuSoCbO/113-01-T1.nii"
## attr(,"nbrOfBytes")
## [1] 44564832

## 
##  #Reorientation /private/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T/RtmpuSoCbO/113-01-T1.nii

## Script created: /var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T//RtmpuSoCbO/filebba0e10e410.m

## Command run is:

## matlab  -nodesktop -nosplash -nodisplay -r  "try, run('/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T//RtmpuSoCbO/filebba0e10e410.m'); catch err, disp(err.message); exit(1); end; exit(0);"

## Checking Filenames

## # Running matlabbatch job

## # Creating matlabbatch job

## # Adding SPMDIR: /var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T//RtmpuSoCbO/spm12\n

## Running job: /var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T//RtmpuSoCbO/filebba0579cebc0.m, which calls /var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T//RtmpuSoCbO/filebba038cbfa49.m

## Command run is:

## matlab  -nodesktop -nosplash -nodisplay -r  "try, run('/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T//RtmpuSoCbO/filebba0579cebc0.m'); catch err, disp(err.message); exit(1); end; exit(0);"

## # Result is 0

## # Removing scripts

Anatomical MRI Segmentation (and Spatial Normalize Estimation)

Here we perform the segmentation of the co-registered anatomical image from above. This will segment the image into 6 different regions, where the regions are gray matter, white matter, cerebrospinal fluid (CSF), bone, soft tissue, and the background. You should inspect these visually before using them to ensure quality and no reordering due to artifacts.

if (have_matlab()) {
  seg_res = spm12_segment(
    filename = coreg_anat,
    set_origin = FALSE,
    retimg = FALSE,
    clean = FALSE,
    install_dir = install_dir
    )
  print(names(seg_res))
}

## # Adding SPMDIR: /var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T//RtmpuSoCbO/spm12

## # Running script /var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T//RtmpuSoCbO/Executable.m
## which calls /var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T//RtmpuSoCbO/Segment_job.m

## Command run is:

## matlab  -nodesktop -nosplash -singleCompThread -nodisplay -r  "try, run('/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T//RtmpuSoCbO/Executable.m'); catch err, disp(err.message); exit(1); end; exit(0);"

## # Result is 0

## [1] "spm"                 "script"              "result"             
## [4] "outfiles"            "outmat"              "deformation"        
## [7] "inverse_deformation"

In order to segment the image, SPM spatially normalizes the image to the MNI template, however. This transformation will be the one we use to transform the fMRI data to MNI space. We see in seg_reg a deformation file, which is the transformation. We also see the output segmentation files of the probability of each class, in native space. We only tend to care about the first 3 categories.

alpha = function(col, alpha = 1) {
  cols = t(col2rgb(col, alpha = FALSE)/255)
  rgb(cols, alpha = alpha)
} 
if (have_matlab()) {
  seg_files = check_nifti(seg_res$outfiles)
  hard_seg = spm_probs_to_seg(seg_files)
  hard_seg[ hard_seg > 3] = 0
  
  ortho2(coreg_img, hard_seg, 
         col.y = alpha(c("red", "green", "blue"), 0.5))
}

bbox = matrix(
  c(-90, -126, -72, 
    90, 90, 108), 
  nrow = 2, byrow = TRUE)
if (have_matlab()) {
  norm = spm12_normalize_write(
    deformation = seg_res$deformation,
    other.files = c(coreg_anat, mean_img, aimg),
    bounding_box = bbox,
    retimg = FALSE,
    install_dir = install_dir
    )
  print(names(norm))
  norm_data = norm$outfiles
  names(norm_data) = c("anat", "mean", "fmri")
  norm_mean_img = readnii(norm_data["mean"])
  norm_anat_img = readnii(norm_data["anat"])
}

## # Adding SPMDIR: /var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T//RtmpuSoCbO/spm12

## # Running script /var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T//RtmpuSoCbO/Executable.m
## which calls /var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T//RtmpuSoCbO/Normalize_Write_job.m

## Command run is:

## matlab  -nodesktop -nosplash -singleCompThread -nodisplay -r  "try, run('/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T//RtmpuSoCbO/Executable.m'); catch err, disp(err.message); exit(1); end; exit(0);"

## # Result is 0

## # Removing scripts

## [1] "spm"          "script"       "deformation"  "other_fnames"
## [5] "outfiles"     "result"

if (have_matlab()) {
  template_path = file.path(spm_dir(), 
                            "canonical", "avg152T1.nii")

  template = readnii(template_path)
  
  dnorm_mask = dnorm_mean_img > quantile(
    dnorm_mean_img[dnorm_mean_img > 0], 
    probs = 0.6)
  norm_mask = norm_mean_img > quantile(
    norm_mean_img[norm_mean_img > 0], 
    probs = 0.6)  
  
  double_ortho(template, norm_anat_img)
  double_ortho(template, norm_mean_img)
  double_ortho(norm_mean_img, norm_anat_img)
  ortho2(template, norm_mask, col.y = alpha("red", 0.5))

  double_ortho(template, dnorm_mean_img)
  ortho2(template, dnorm_mask, col.y = alpha("red", 0.5))
  double_ortho(norm_mean_img, dnorm_mean_img)
}

Spatial Smoothing

Here we will perform spatial smoothing of the 4D fMRI data in template space. We can set the full-width half max (FWHM) for the Gaussian smoother. The relationship between the FWHM and the Gaussian standard deviation is:

\[ FWHM = \sigma \sqrt{8 \log(2)} \] where \(\log\)` is the natural log.

if (have_matlab()) {
  smoothed = spm12_smooth(
    filename = norm_data["fmri"],
    fwhm = 8,
    prefix = "s",
    retimg = FALSE,
    install_dir = install_dir
    )
  smoothed_data = smoothed$outfiles
}

## # Adding SPMDIR: /var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T//RtmpuSoCbO/spm12

## # Running script /var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T//RtmpuSoCbO/Executable.m
## which calls /var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T//RtmpuSoCbO/Smooth_job.m

## Command run is:

## matlab  -nodesktop -nosplash -singleCompThread -nodisplay -r  "try, run('/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T//RtmpuSoCbO/Executable.m'); catch err, disp(err.message); exit(1); end; exit(0);"

## # Result is 0

## # Removing scripts

In many applications, this is the data you will use for post-processing and analysis. Motion correction has usually been applied above, but some motion correct this data as well.

if (have_matlab()) {
  smoothed_mean = spm12_smooth(
    filename = norm_data["mean"],
    prefix = "s",
    retimg = FALSE,
    install_dir = install_dir
    )
  smoothed_mean_data = smoothed_mean$outfiles
}

## # Adding SPMDIR: /var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T//RtmpuSoCbO/spm12

## # Running script /var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T//RtmpuSoCbO/Executable.m
## which calls /var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T//RtmpuSoCbO/Smooth_job.m

## Command run is:

## matlab  -nodesktop -nosplash -singleCompThread -nodisplay -r  "try, run('/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T//RtmpuSoCbO/Executable.m'); catch err, disp(err.message); exit(1); end; exit(0);"

## # Result is 0

## # Removing scripts

Here we can smooth the mean image in MNI space. We could have also calculated a new mean image and smoothed that. Or we could have calculated a mean of the smoothed 4D data. They all may be slightly different.

if (have_matlab()) {
  smooth_mean_img = readnii(smoothed_mean_data)
  ortho2(smooth_mean_img)
}

Bibliography

bib = '@article{ashburner2005unified,
  title={Unified segmentation},
  author={Ashburner, John and Friston, Karl J},
  journal={Neuroimage},
  volume={26},
  number={3},
  pages={839--851},
  year={2005},
  publisher={Elsevier}
}'