Enantiomorphic Filling and Registration

[brussj]

I just read the excellent chapter, "Mapping the Spatial Distribution of Lesions in Stroke: Effect
of Diffeomorphic Registration Strategy in the ATLAS Dataset" by Avants & Tustison from Pustina and Mirman's "Lesion-to-Symptom Mapping" book.

I tend to favor bash commands over R but, nonetheless, I tried to implement your strategy on a lesioned brain. I was able to run the R commands and output the reflected brain. I set verbosity to TRUE then attempted to mimic the commands via ImageMath, antsRegistration. My results are very close to the ANTsR implementation but not quite the same.

For reference, say I have a T1 (T1.nii.gz). Using ANTs:

subjectBrain = antsImageRead("T1.nii.gz") %>% iMath("Normalize")
subjectBrainReflect = reflectImage( subjectBrain, axis = 0, tx = "SyN",verbose = TRUE, totalSigma = 0.5 )$warpedmovout
antsImageWrite(subjectBrainReflect, "T1_reflected_R.nii.gz")

Using command-line/bash:

ImageMath 3 T1_norm.nii.gz Normalize T1.nii.gz 1

ImageMath 3 reflectionMatrix.mat ReflectionMatrix T1_norm.nii.gz 0

antsRegistration -d 3 -o [T1_reflected_,T1_reflected.nii.gz] \
-r reflectionMatrix.mat \
-n Linear \
-t Affine[0.25] \
-m MI[T1_norm.nii.gz,T1_norm.nii.gz,1,32,Regular,0.2] \
-c [2100x1200x1200x0,1e-06,10] \
-f 4x2x2x1 \
-s 3x2x1x0 \
-t SyN[0.2,3,0.5] \
-m MI[T1_norm.nii.gz,T1_norm.nii.gz,1,32,None] \
-c [40x20x0,1.0000e-07,8] \
-f 4x2x1 \
-s 2x1x0

My bash results are quite good but something is a bit off (ever so slightly different rotation on the output image). I would be happy to upload my sample T1 if needed. Did I make any glaring mistakes or could someone recommend and tweaks to my parameters? Thank you for taking the time to read this and for your very informative book chapter.

[ntustison]

Yes, can you upload the images so we can recreate on our end?

[brussj]

rawData.zip

Apologies for my extremely late reply. I re-ran this just to double check, with slightly altered workflows:

#Remove noise, then normalize
DenoiseImage -d 3 -i rawData/T1.nii.gz -n Rician -x rawData/T1mask.nii.gz -o procData/T1_denoised.nii.gz
ImageMath 3 procData/T1_norm.nii.gz Normalize procData/T1_denoised.nii.gz 1

#Make sure brain and lesin masks are binary
ThresholdImage 3 rawData/T1mask_raw.nii.gz procData/T1mask.nii.gz 1 Inf
ThresholdImage 3 rawData/lesion_raw.nii.gz procData/lesion.nii.gz 1 Inf

#Apply the brain mask to the normalized T1
ImageMath 3 procData/T1_brain.nii.gz m procData/T1_norm.nii.gz procData/T1mask.nii.gz

#Create a reflection matrix
ImageMath 3 procData/reflectionMatrix.mat ReflectionMatrix procData/T1_brain.nii.gz 0

#Create the enantiomorph, using the reflection matrix as a starting point
antsRegistration -d 3 -o [procData/T1_reflected_,procData/T1_reflected.nii.gz] \
-r procData/reflectionMatrix.mat \
-n Linear \
-t Affine[0.25] \
-m MI[procData/T1_brain.nii.gz,procData/T1_brain.nii.gz,1,32,Regular,0.2] \
-c [2100x1200x1200x0,1e-06,10] \
-f 4x2x2x1 \
-s 3x2x1x0 \
-t SyN[0.2,3,0.5] \
-m MI[procData/T1_brain.nii.gz,procData/T1_brain.nii.gz,1,32,None] \
-c [40x20x0,1.0000e-07,8] \
-f 4x2x1 \
-s 2x1x0

#Smooth the lesion mask, normalize, invert
SmoothImage 3 procData/lesion.nii.gz 1 procData/lesionSmooth.nii.gz
ImageMath 3 procData/lesionSmooth_norm.nii.gz Normalize procData/lesionSmooth.nii.gz 0
ImageMath 3 procData/lesion_Inv.nii.gz Neg procData/lesionSmooth_norm.nii.gz

#Pull out the lesioned region, create a plug from the enantiomorph, fill in the hole
ImageMath 3 procData/T1_hole.nii.gz m procData/T1_brain.nii.gz procData/lesion_Inv.nii.gz
ImageMath 3 procData/T1_plug.nii.gz m procData/T1_reflected.nii.gz procData/lesionSmooth_norm.nii.gz
ImageMath 3 procData/T1_Filled_bash.nii.gz + procData/T1_hole.nii.gz procData/T1_plug.nii.gz

##########

#Try the same method but using ANTsR  (Import data from command-line steps)
subjectBrain = antsImageRead("procData/T1_brain.nii.gz")
subjectLesionMask = antsImageRead( "procData/lesion.nii.gz" )
subjectBrainReflect = reflectImage( subjectBrain, axis = 0, tx = "SyN",verbose = TRUE, totalSigma = 0.5 )$warpedmovout
slesion = antsImageRead("procData/lesionSmooth_norm.nii.gz")
slesionNeg = antsImageRead("procData/lesion_Inv.nii.gz")
subjectFill = subjectBrain * slesionNeg + subjectBrainReflect * slesion
antsImageWrite(subjectFill, "procData/T1_Filled_antsR.nii.gz")
antsImageWrite(subjectBrainReflect, "procData/T1_reflected_antsR.nii.gz")

The difference between the two workflows is minimal. It's easiest to see if you load up in something like ITK-SNAP and turn one of the layers on and off. Thanks for taking a look at this.