Coregistering low-resolution, restricted z-coverage scans to the PAM 50 template, ideally in an automated fashion?

Dear SCT/coregistration experts!

Problem
I’d like your help with coregistering low-resolution, restricted FOV (covering only a few vertebrae) scans to the PAM50 template, ideally in an automated way. As far as I understand, the limitation will be generating and supplying vertebral labels.

Description of acquisitions
In depth, we have scans centered around C6/7 with the following parameters:12 slices of 3mm, matrix=128x56, image resolution=2.2x2.2, 6 echoes (I’ll refer to these as lowres_GRE_S, and the mean across echoes as Lowres_GRE_mean) that we interrogate.
Additional acq: We use a T1w (Lowres_run1_T1w) scan we use for segmentation, etc, with: 22 slices of 3mm, matrix=256x256, image resolution=0.9x0.9. These could be used as a fallback, for a 2-step coregistration (lowres_GRE to T1, T1 to PAM50), but I would like to avoid piling on warping fields.

Note that we use the T1 to derive the spinal cord segmentation I boundled in the zip. (We co-register Lowres_GRE_mean to the T1, segment the T1, then warp the segmentation back).

Desired outcome
I’d like to have Lowres_GRE_mean in the space of PAM50, so that I can then use the warping field to warp the multi-echo scan. Further processing will involve meaning across subjects, etc.

Previous processing steps
We do some unrelated processing, which is this:

  • Creating the mean across echoes (Lowres_GRE_mean)
  • Coregistering Lowres_GRE_mean to the T1
  • Segmenting the T1, and warping that segmentation back (spinal_seg_reg_S)

Solutions tried
I tried the following processing pathway:
1, Creating labels for the disks, without initialisation. This results in the following error (not surprising, as C2/C3 are not within the FOV):

Automatic C2-C3 detection failed. Please provide manual label with sct_label_utils

2, Manually initialising, such that the central slice (slice 6) corresponds to C6/C7:
sct_label_vertebrae -i Lowres_GRE_mean.nii.gz -s spinal_seg_reg_S.nii.gz -c t1 -initz 6,7

This does result in a reasonable labeling, but some bottom slices have no labels, and the 8 label (T1?) is too short:

3, I can then use this label file to drive the coregistration:
sct_register_to_template -i Lowres_GRE_mean.nii.gz -s spinal_seg_reg_S.nii.gz -l spinal_seg_reg_S_labeled.nii.gz

I am not happy with this result. It looks really warped in the long axis view:
image

And some of the slices look very distorted and highly interpolated:

Questions:
1, Is there a way to improve this outcome? I am not familiar with to-template registration, so I don’t know if the limitation is the input image quality, or my processing.
2, Is it preferable to derive the labels on the T1w scan as well, and warp them back the same way? I generally prefer not to have too many warp fields one after another, and the T1 labeling also does not look all that great (crosshair at label 6). Am I doing it wrong?

sct_label_vertebrae -i Lowres_run1_T1w.nii.gz -c t1 -initz 11,7 -s spinal_seg_T1.nii       

Data is here:

Best

Daniel

Hi @Daniel_Papp !

As you rightfully pointed out, getting the vertebral levels where C1-C3 are not visible, is not possible without prior information. In this case, I would use the prior information about FOV placement (ie: centered at C6/C7) to get the automatic disc label. Also, given (i) such a small I-S coverage and (ii) the large slice thickness, I would not worry about estimating an affine transformation along the S-I direction. As you noted, it could result in ugly stretching along the S-I direction, associated with interpolation errors.

So, in this case, I prefer to simply register with a single disc label, and then run a slicewise registration algorithm. This scenario is illustrated in this tutorial. The downside is that it will use the average distance between two adjacent discs from the PAM50 template, and not from this subject, but again, given the thickness of the slices, I don’t think this is a big problem here. Also, obviously, this concern depends on your application (which I don’t know what it is).

Another thing to note, is that the in-plane resolution of the ME-GRE is very low compared to that of the PAM50 (2mm vs. 0.5mm):

image

Therefore, if you would like to use this image to register to the PAM50 directly (and not go via the T1w image as you pointed out), then I assume that you don’t care about precision/accuracy, and that you are only interested in getting a ‘coarse’ registration to the PAM50.

So, based on the data you shared and your working specs, here is a proposed script:

# Segment cord on T1w
sct_deepseg_sc -i Lowres_run1_T1w.nii.gz -c t1 -qc qc
# Generate mask around the cord (for subsequent registration)
sct_create_mask -i Lowres_run1_T1w.nii.gz -p centerline,Lowres_run1_T1w_seg.nii.gz -size 35mm -f cylinder
# Register the T1w scan to the mean ME-GRE scan
# Tips: Here we specify '-dseg Lowres_run1_T1w_seg.nii.gz' only to be able to generate the QC report. Even thought 'Lowres_run1_T1w_seg.nii.gz' corresponds to iseg (not dseg), we don't really care because the segmentation is *not* used for registration (it is only used to center the spinal cord during QC). 
sct_register_multimodal -i Lowres_run1_T1w.nii.gz -d Lowres_GRE_mean.nii.gz -dseg Lowres_run1_T1w_seg.nii.gz -m mask_Lowres_run1_T1w.nii.gz -param step=1,type=im,algo=slicereg,metric=CC -qc qc
# Apply the transformation to the cord segmentation
# Tips: we use trilinear interpolation because having the 'soft' segmentation will provide more accurate registration to the template
sct_apply_transfo -i Lowres_run1_T1w_seg.nii.gz -d Lowres_GRE_mean.nii.gz -w warp_Lowres_run1_T1w2Lowres_GRE_mean.nii.gz -x linear
# Create label centered in the cord at level C6/C7 
sct_label_utils -i Lowres_run1_T1w_seg_reg.nii.gz -create-seg-mid 7 -o label_c6c7.nii.gz -qc qc
# Register the ME-GRE scan to the PAM50 template
# Tips: we use '-ref subject' because we don't want straightening (due to the thick slices).
# Tips: we only use slice-wise rigid translation, based on the segmentation, because due to the low resolution of the ME-GRE image, any non-affine registration would be 'risky'.
sct_register_to_template -i Lowres_GRE_mean.nii.gz -s Lowres_run1_T1w_seg_reg.nii.gz -ldisc label_c6c7.nii.gz -ref subject -param step=1,type=seg,algo=centermassrot -qc qc

I would personally opt for another approach: register the PAM50 to the T1w image, and then concatenate the warping field PAM50–>T1w–>GRE, because the T1w is available and it has a much better resolution.

Also, as said above, this registration will produce inaccurate placement of the discs, but I assume that for this application this is not a problem:

anim

Naive question: what do you need the registration to the PAM50 for? Maybe it is not needed?

Dear @jcohenadad!

Thank you very much for the clarifications.

To answer your question, our end goal is to have all of our subjects in a common space. We are acquiring these MEGRE scans under several conditions, and we want to compare signal intensities along the cord, across these conditions for all subjects, similar to the recent preprint by Kaptan et al.
image

I have tried to go down the route of coregistering all subjects to a destination subject, but due to the poor quality of the MEGRE scans, that requires applying three warps, using the T1s as an in-between (source MEGRE to source T1; source T1 to destination T1; destination T1 to destination source). Those results did not satisfy me, and I am exploring alternatives.

I much prefer going down your suggested route of PAM50 registration

For comparison, here are the two pipelines for PAM50 regiistration:
Direct coregistration to ME-GRE to PAM50
This pipeline is almost identical to your proposal, only I switched “-x linear” to “-x nn” in the step where the segmentation is warped back to the ME-GRE. Trilinear interpolation gives non-binary values in the segmentation, see Note 1 on the bottom.

# Segment cord on T1w
sct_deepseg_sc -i Lowres_run1_T1w.nii.gz -c t1 -qc qc
# Generate mask around the cord (for subsequent registration)
sct_create_mask -i Lowres_run1_T1w.nii.gz -p centerline,Lowres_run1_T1w_seg.nii.gz -size 35mm -f cylinder
# Register the T1w scan to the mean ME-GRE scan
# Tips: Here we specify '-dseg Lowres_run1_T1w_seg.nii.gz' only to be able to generate the QC report. Even thought 'Lowres_run1_T1w_seg.nii.gz' corresponds to iseg (not dseg), we don't really care because the segmentation is *not* used for registration (it is only used to center the spinal cord during QC). 
sct_register_multimodal -i Lowres_run1_T1w.nii.gz -d Lowres_GRE_mean.nii.gz -dseg Lowres_run1_T1w_seg.nii.gz -m mask_Lowres_run1_T1w.nii.gz -param step=1,type=im,algo=slicereg,metric=CC -qc qc
# Apply the transformation to the cord segmentation
# Tips: we use trilinear interpolation because having the 'soft' segmentation will provide more accurate registration to the template
sct_apply_transfo -i Lowres_run1_T1w_seg.nii.gz -d Lowres_GRE_mean.nii.gz -w warp_Lowres_run1_T1w2Lowres_GRE_mean.nii.gz -x linear
# Create label centered in the cord at level C6/C7 
sct_label_utils -i Lowres_run1_T1w_seg_reg.nii.gz -create-seg-mid 7 -o label_c6c7.nii.gz -qc qc
# Register the ME-GRE scan to the PAM50 template
# Tips: we use '-ref subject' because we don't want straightening (due to the thick slices).
# Tips: we only use slice-wise rigid translation, based on the segmentation, because due to the low resolution of the ME-GRE image, any non-affine registration would be 'risky'.
sct_register_to_template -i Lowres_GRE_mean.nii.gz -s Lowres_run1_T1w_seg_reg.nii.gz -ldisc label_c6c7.nii.gz -ref subject -param step=1,type=seg,algo=centermassrot -qc qc

Using the T1w as an intermediary step

# Segment cord on T1w
sct_deepseg_sc -i Lowres_run1_T1w.nii.gz -c t1 -qc qc
# Generate mask around the cord (for subsequent registration)
sct_create_mask -i Lowres_run1_T1w.nii.gz -p centerline,Lowres_run1_T1w_seg.nii.gz -size 35mm -f cylinder
# Register the T1w scan to the mean ME-GRE scan
# Tips: Here we specify '-dseg Lowres_run1_T1w_seg.nii.gz' only to be able to generate the QC report. Even thought 'Lowres_run1_T1w_seg.nii.gz' corresponds to iseg (not dseg), we don't really care because the segmentation is *not* used for registration (it is only used to center the spinal cord during QC). 
sct_register_multimodal -i Lowres_run1_T1w.nii.gz -d Lowres_GRE_mean.nii.gz -dseg Lowres_run1_T1w_seg.nii.gz -m mask_Lowres_run1_T1w.nii.gz -param step=1,type=im,algo=slicereg,metric=CC -qc qc

#New steps start here
# Create label centered in the cord at level C6/C7 on the T1w!
sct_label_utils -i Lowres_run1_T1w_seg.nii.gz -create-seg-mid 7 -o label_c6c7.nii.gz -qc qc
# Coregistering the T1w to the PAM50. I am using your proposed one-step, slicewise method
sct_register_to_template -i Lowres_run1_T1w.nii.gz -s Lowres_run1_T1w_seg.nii.gz -ldisc label_c6c7.nii.gz -ref subject -param step=1,type=seg,algo=centermassrot -qc qc
# Concatenating the warp fields
sct_concat_transfo -d /Users/danielpapp/spinalcordtoolbox/data/PAM50/template/PAM50_t1.nii.gz -w warp_Lowres_GRE_mean2Lowres_run1_T1w.nii.gz warp_anat2template.nii.gz -o warp_MEGRE_2_PAM50.nii.gz
# Finally, applying the concatinated warp field to the MEGRE
sct_apply_transfo -d /Users/danielpapp/spinalcordtoolbox/data/PAM50/template/PAM50_t1.nii.gz -w warp_MEGRE_2_PAM50.nii.gz -i Lowres_GRE_mean.nii.gz

The two results are quite dissimilar. The direct coregistration seems stretched in the Z direction, compared to the 2-step process, and there is, for lack of a better word, a wavy, sinusoidal nature to the saggital view. The gif loops as PAM50–>direct coregistration–>2step coregistration
anim

In the axial view, there is a rotation in the direct coregistration that is missing (I assume this is a good thing) in the 2step coregisration. Gif loops the same as previous. The images still look interpolated, but I guess that is just the result of the low resolution for the initial input.
anim

I will stick with the 2-step registration to PAM50, based on these results.

Note 1:
Trilinear interpolation on the spinal cord segmentation results in non-binary values within the warped segmentation. I do not know how SCT deals with this. Naively, I would assume that non-binary values are ignored.
anim

Note 2: The files resulting from coregistration to PAM50 are saved as 64-bit float, I assume to accomodate negative floating point values that result from coregistration. Is there a way to tell SCT during the coreg step to save the data as INT16 (which is enough for our purposes?)

I have been using the following hack, but I don’t know if sct_register_to_template has a hidden “-type” flag too

sct_maths -i anat2template.nii.gz -o anat2template_16.nii.gz -add 0 -type int16

Thanks!
Daniel

OK, so yes, you need registration to the PAM50 then (as you are already well aware).

That’s only two ‘global’ transformations: GRE → T1w and T1w → PAM50. I call them ‘global’ transformations because in reality there are sub-transformations that are already concatenated by sct_register_multimodal (e.g. when using multiple-step approach). Concatenating two global warping fields is very reasonable so I wouldn’t worry about that.

Few comments on your proposed pipelines:

  • never put hard-coded path (for reproducibility). You can replace /Users/danielpapp/spinalcordtoolbox by $SCT_DIR
  • if you are using the T1w as an intermediate step, then you can afford doing an image-based non-affine as step 2. Look at the QC report and compare lines 4-5 so you can appreciate the difference with/without step2.
  • Since the T1 is used for registration to the PAM50, the PAM50 contrast to choose is the T1 (not T2, which is the default), hence the flag -c t1

Here is a suggestion:

# Segment cord on T1w
sct_deepseg_sc -i Lowres_run1_T1w.nii.gz -c t1 -qc qc
# Generate mask around the cord (for subsequent registration)
sct_create_mask -i Lowres_run1_T1w.nii.gz -p centerline,Lowres_run1_T1w_seg.nii.gz -size 35mm -f cylinder
# Register the T1w scan to the mean ME-GRE scan
# Tips: Here we specify '-dseg Lowres_run1_T1w_seg.nii.gz' only to be able to generate the QC report. Even thought 'Lowres_run1_T1w_seg.nii.gz' corresponds to iseg (not dseg), we don't really care because the segmentation is *not* used for registration (it is only used to center the spinal cord during QC). 
sct_register_multimodal -i Lowres_run1_T1w.nii.gz -d Lowres_GRE_mean.nii.gz -dseg Lowres_run1_T1w_seg.nii.gz -m mask_Lowres_run1_T1w.nii.gz -param step=1,type=im,algo=slicereg,metric=CC -qc qc
# Create label centered in the cord at level C6/C7 on the T1w
sct_label_utils -i Lowres_run1_T1w_seg.nii.gz -create-seg-mid 7 -o label_c6c7.nii.gz -qc qc
# Register the T1w to the PAM50. I am using your proposed one-step, slicewise method
sct_register_to_template -i Lowres_run1_T1w.nii.gz -s Lowres_run1_T1w_seg.nii.gz -ldisc label_c6c7.nii.gz -c t1 -ref subject -param step=1,type=seg,algo=centermassrot:step=2,type=im,algo=bsplinesyn,metric=CC,iter=5,slicewise=1 -qc qc
# Concatenate the warping fields
sct_concat_transfo -d $SCT_DIR/data/PAM50/template/PAM50_t1.nii.gz -w warp_Lowres_GRE_mean2Lowres_run1_T1w.nii.gz warp_anat2template.nii.gz -o warp_MEGRE_2_PAM50.nii.gz
# Apply the concatenated warping field to the MEGRE
sct_apply_transfo -d $SCT_DIR/data/PAM50/template/PAM50_t1.nii.gz -w warp_MEGRE_2_PAM50.nii.gz -i Lowres_GRE_mean.nii.gz

Here’s the QC report: qc.zip (590.5 KB)

Some algorithms in SCT (eg: centerline detection) rely on the center of mass, in which case having soft (vs. hard a.k.a ‘binary’) segmentations results in more precise results.

1 Like