Dear All,
A program that uses afni registers a ct scan ($ct) to T1 MRI ($t1) via 3 steps
-
@Align_Centers -base $t1 -dset $ct : creates matrix ct_shft.1D
- 3dresample -input ct_shft.nii -prefix ct_res_shft.nii -master $t1 -dxyz 1 1 1 -rmode NN
- align_epi_anat.py -dset1 $t1 -dset2 ct_res_shft.nii -dset1_strip None -dset2_strip None -dset2to1 -suffix _al -feature_size 1 -overwrite -cost nmi -giant_move -rigid_body > status.txt : creates matrix t1_al_mat.aff12.1D
I would like to create also a T1 copy registered to $ct (reverse transformation of t1 in ct space) using the two matrices generated by forward transform. How to I do this in afni?
Thank you very much and Happy New Year!!
Octavian
Hi, Octavian-
Question: why is there a resampling step in #2? That seems to me to be an unnecessary step, and one that blurs the data in the process.
–pt
Dear Dr. Taylor,
I inherited this code from ALICE toolbox. I ran the registration without step 2 (resampling), and it does not work (step 3 align_epi_anat.py runs forever to converge) and the result on visual check in afni is not ok. The only way the code works is as is, including step 2 (using the resampling at the native T1 resolution, no other resolution -dxyz works, see also the thread ‘align_epi_anat.py registration variability’ about a mo ago). If you can take a look using the actual data (t1, ct), I can upload those to you (I would need to know how). I am not sure to what extent the need for step 2 is because our CTs are anisotropic (0.43, 0.43, 0.8) vs isotropic MRI (1,1,1,)
If using step 2 as is, the registration is good. Then, I would need to know how to move a ribbon.nii (mask generated in freesurfer, aligns with T1) to the ct space, I have tried:
cat_matvec -MATRIX -ONELINE ct_shft.1D -I > inv_ct_shft.1D #inverting forward matrix from step 1 of the registration
cat_matvec -MATRIX -ONELINE ct_res_shft_al_mat.aff12.1D -I > inv_ct_res_shft_al_mat.aff12.1D #inverting forward matrix from step 3
3dAllineate -source ribbon.nii -base ct_res_shft.nii -prefix ribbon_shft.nii -1Dmatrix_apply inv_ct_res_shft_al_mat.aff12.1D
3dresample -input ribbon_shft.nii -prefix ribbon_res_shft.nii -master $ct -rmode NN
3dAllineate -source ribbon_res_shft.nii -base $ct -prefix ribbon_ct.nii -1Dmatrix_apply inv_ct_shft.1D
For some reason, the output ribbon_ct.nii looks aligned to the CT, but appears truncated, in afni.
Thank you very much, Octavian
I think I solved it. I changed the 3dresample step above
3dresample -input ribbon_shft.nii -prefix ribbon_res_shft.nii -master $ct -rmode NN
to
3dresample -input ribbon_shft.nii -prefix ribbon_res_shft.nii -master ct_shft.nii -rmode NN
Anyway, it there is a more straightforward way of doing it please let me know. Also, I would like to know which -rmode option (NN, cubic) or other is the ‘best’, irrespective of computation time. Thank you,
Octavian
Hi, Octavian-
Hmm, OK. I’m still a liiiittle bit curious about this. Would you be able to share the dataset if I sent you a Box link?
Re. resampling mode: typically, we use “NN” for integer-based dsets (like atlases and ROI maps), and cubic would be more appropriate for a continuous/floating point-valued dset (like template, anatomical, EPI, etc.).
–pt
Sure, please send me a box link, if email needed then octavian.lie@gmail.com, thank you.
Hi, Octavian-
Thanks for sharing the data.
I now understand about the resampling-- the CT is at very high res (finer than 0.5 mm isotropic voxels). Aligning that to a high res anatomical is probably going to take a while and be memory intensive.
I took your initial set of commands and wrote the following script.
- For one, with multiple steps, I like making a working directory (easier to move and re-run). This new script should also be more agnostic to initial file names (which is a bit awkward at times here, because both @Align_Centers and align_epi_anat.py use suffixes rather than straightup output prefixes).
- Also, I took your 1x1x1 mm3 resampling and raised you another mm, to 2x2x2 mm3: the CT images have so little contrast, why not downsample them further to speed things up? The alignment is nearly identical in either case, so the significant speedup in runtime seems worth it. The full transform can be applied to the original dset from the beginning (as done in the new script, and using the ‘wsinc5’ interpolant, which preserves edges better than any of the others).
- Also, another step to save time: I added 1000 to the CT image, because it seems like the baseline floor is about -1000 (like, that is what the air outside the head is); doing so allowed me to run 3dAutomask within the align_epi_anat.py command, which sped things up a bit. NB: the full warps get applied to the original, unscaled data, so this shouldn’t affect your final results of interest.
- I calculate both the full forward and full inverse transform (ct → t1, and t1 → ct, respectively). Each also gets applied-- NB: applying the full (concatenated) warp to each dset means there is minimal smoothing from resampling, through there will always be some (but it is much better than resampling over several steps).
- I added in QC images throughout, to show the alignment at various points: starting, after @Align_Centers, after align_epi_anat.py, and then for each of the forward/backward applied warps.
Feel free to adapt this further. You could copy anything you want out of the working dir, for example.
#!/bin/tcsh
set in_t1 = T1.nii
set in_ct = CT_highresRAI.nii
set here = ${PWD}
set wdir = ./__workdir_align_t1_ct
set ct = ct_00.nii
set t1 = t1_00.nii
# --------------------------------------------------------------------
\mkdir -p ${wdir}
# cp dsets to wdir for working
3dcopy ${in_t1} ${wdir}/$t1
3dcopy ${in_ct} ${wdir}/$ct
cd ${wdir}
# the CT appears to have a "floor" background value around -1000;
# scale it up, so we can automask the data for more speed.
3dcalc \
-a $ct \
-expr 'a+1000' \
-prefix ct_01_scale.nii
@djunct_overlap_check \
-ulay $t1 \
-olay ct_01_scale.nii \
-prefix qc_01_init
# creates matrix ct_shft.1D
@Align_Centers \
-base $t1 \
-dset ct_01_scale.nii
@djunct_overlap_check \
-ulay $t1 \
-olay ct_01_scale_shft.nii \
-prefix qc_02_shft
# here, downsampling a lot because the CT has so little contrast/info
3dresample \
-input ct_01_scale_shft.nii \
-prefix ct_02_res.nii \
-master $t1 \
-dxyz 2 2 2 \
-rmode NN
# creates matrix t1_al_mat.aff12.1D
# giant_move appears necessary
align_epi_anat.py \
-overwrite \
-dset1 $t1 \
-dset2 ct_02_res.nii \
-dset1_strip 3dAutomask \
-dset2_strip 3dAutomask \
-dset2to1 \
-suffix _al \
-feature_size 1 \
-cost nmi \
-giant_move \
-rigid_body \
> status.txt
@djunct_overlap_check \
-ulay $t1 \
-olay ct_02_res_al*HEAD \
-prefix qc_03_aff
# ----------------------------------------------------------
# calculate full forward transform (ct -> t1)
cat_matvec \
-ONELINE \
ct_02_res_al_mat.aff12.1D \
ct_01_scale_shft.1D \
> mat_ct_to_t1.aff12.1D
# calculate full inverse transform (t1 -> ct)
cat_matvec \
-ONELINE \
mat_ct_to_t1.aff12.1D -I \
> mat_t1_to_ct.aff12.1D
# ---- apply each full transform (showing each works) ------
3dAllineate \
-1Dmatrix_apply mat_ct_to_t1.aff12.1D \
-source ${ct} \
-master ${t1} \
-prefix CT_in_T1.nii \
-final wsinc5
3dAllineate \
-1Dmatrix_apply mat_t1_to_ct.aff12.1D \
-source ${t1} \
-master ${ct} \
-prefix T1_in_CT.nii \
-final wsinc5
@djunct_overlap_check \
-ulay $t1 \
-olay CT_in_T1.nii \
-prefix qc_04_ct_in_t1
@djunct_overlap_check \
-ulay $ct \
-olay T1_in_CT.nii \
-prefix qc_05_t1_in_ct
echo "++ Done."
–pt
Wow, thank you Peter! This is rather prompt, I will try the code and let you know of any problems, thank you. Octavian