Obtain spinal levels from nerve rootlets segmentation

[tzebre]

UPDATE by @valosekj -- I updated the issue description to be in line with #18.

The model predicts spinal nerve rootlets segmentation; see here. The predicted spinal nerve rootlets can be used to obtain the spinal levels using the utilities/rootlets_to_spinal_levels.py script using the following procedure:

1. Dilate the spinal cord segmentation by 1, 2 or 3 voxels (-dilate input argument)
2. Find the intersection between the dilated spinal cord segmentation and the rootlets segmentation.
3. The spinal levels are then defined based on the top and bottom slices of the intersection.

Individual steps are illustrated in this Google Slides presentation.

---

To validate the result, we compared the results of this prediction with Cadotte et al., 2017 values using the same MRI images.
Cadotte data have labeled images with PMJ (pontomedullary junction), vertebrae and DREZ (dorsal root entry zone) labeled in the same file. I used separate_cadotte_label.py to create 3 files: one for the PMJ, one for the DREZ, and one for vertebrae.
I've used cadotte_value notebook to create cadotte_dist.csv
The results of the prediction on Cadotte are not really good #11.

**Cadotte paper graph**

**Reproduced graph**

Legend: Left = Cadotte DREZ label, Middle = Cadotte label vertebrae, Right = M3 prediction (conversion with dilate = 2)

[valosekj]

Raphaëlle Schlienger noticed that some subjects have nerve rootlets segmentation quite far from the spinal cord, which thus leads to "short" spinal levels.
In such cases, I would recommend increasing the dilation of the spinal cord to -dilate 3:

Legend: blue - dilate3, red - dilate2

[valosekj]

@jcohenadad proposed an idea to take the highest Z value (i.e., in superior-inferior direction) found for each spinal nerve rootlet level. Given that rootlets are oriented caudally as they go away from the cord (see image below), the most rostral point of the segmented rootlets should be the closest to the actual spinal level.

Image source

But, this might not work for C1 and C2 as these rootlets are "upside down" compared to the other rootlets.

[RaphaSchl]

This means you would be taking the highest Z-value to define the start of the spinal level.
But this highest Z-value is when we 'lose sight' of the rootlet as it enters the white+grey matter area. Yet we can imagine it still has to cross through white matter before synapsing in the innermost grey matter of the cord (butterfly shape on the right MERGE/MEDIC T2* image).

Hence some activity linked to that spinal level may occur above the highest Z-value observed for the corresponding rootlet. And with a precision of 0.8mm iso as in this T2 image used as reference (on the left), we can't be sure we are ultimately observing the highest point of entry.

Maybe the spinal levels should thus start a bit above this highest Z-value to account for this uncertainty?

[jcohenadad]

But this highest Z-value is when we 'lose sight' of the rootlet as it enters the white+grey matter area.

Yes, I agree. This is precisely the reasoning for selecting the highest z-value (vs. doing the dilation). Even though the "highest z-value" is not perfect because as you pointed out, we don't exactly know where the dorsal entries synapse inside the cord, this is the closest "bet" we can make given the coarse resolution of MRI.

Maybe the spinal levels should thus start a bit above this highest Z-value to account for this uncertainty

Sure, but what is "a bit above"? We would need some quantitative evidence to come up with a model that would systematically "shift up" the suggested spinal levels. Without this quantitative evidence, IMHO the best we can do is to rely on what we "see" on the MRI.

Maybe the Mendez et al., 2021 paper (see spinalcordtoolbox/PAM50#3) includes such useful information? E.g. exact location of spinal levels based on the visualization of the nerve rootlets visible before they enter the spinal cord. @sandrinebedard might have some insights about this

[sandrinebedard]

Sure, but what is "a bit above"? We would need some quantitative evidence to come up with a model that would systematically "shift up" the suggested spinal levels. Without this quantitative evidence, IMHO the best we can do is to rely on what we "see" on the MRI.

Maybe the Mendez et al. paper (see spinalcordtoolbox/PAM50#3) includes such useful information? E.g. exact location of spinal levels based on the visualization of the nerve rootlets visible before they enter the spinal cord. @sandrinebedard might have some insights about this

In the Mendez we have Width Across Dorsal and Ventral Columns (between the nerve root entry) and you have the angles of rostral and caudal rootlets per level, maybe with that you could compute the extra z you should add!

However, you need to consider that the angles of the rootlets changes with flexion/extension of the head, so it may not be accurate.

[jcohenadad]

In the Mendez we have Width Across Dorsal and Ventral Columns (between the nerve root entry) and you have the angles of rostral and caudal rootlets per level, maybe with that you could compute the extra z you should add!

Brilliant idea! That would be very cool indeed. Though it might not be super robust to estimate an angle based on a very noisy segmentation. What might be much more robust would be to fit a geometric model of the rootlets to the segmented labels.

For example, we could segment the rootlets in the PAM50 space. Then, for every new subject, do:

• segmentation of nerve rootlets in native space
• registration to PAM50
• bring rootlets segmentation from native space to PAM50 space
• label-wise registration between PAM50 rootlets and subject rootlets (this could be a Tz-only registration, and we could potentially use some regularization across all rootlets registration because all rootlets are somewhat geometrically linked)
• concatenate transformation PAM50rootlets --> nativeRootletsInPAM50space --> nativeRootletsInNativeSpace
• get the spinal levels from the PAM50 (defined with Frostell et al.) and apply the inverse transform PAM50levels --> PAM50levelsInNativeSpace

[valosekj]

I am working on the solution described in the previous comment. So far, I have:

1. segmented the spinal nerve rootlets in the PAM50 T2w image (by running our nnUNet model) -- prediction looks pretty good:

2. segmented the spinal nerve rootlets in the native space -- here I use sub-007_ses-headNormal and the STAPLE reference segmentation (sub-007_ses-headNormal_T2w_label-rootlet_staple.nii.gz) obtained from four raters (details in #17)

3. registered the native T2w image to PAM50 template:

file_t2=sub-007_ses-headNormal_T2w

# Segment spinal cord
sct_deepseg_sc -i ${file_t2}.nii.gz -c t2 -qc qc -qc-subject ${file_t2}

# Create mid-vertebral labels in the cord for vertebrae C3 and C5 
sct_label_vertebrae -i ${file_t2}.nii.gz -s ${file_t2}_seg.nii.gz -c t2 -qc qc -qc-subject ${file_t2}
sct_label_utils -i ${file_t2}_seg_labeled.nii.gz -vert-body 3,5 -o ${file_t2}_seg_labeled_vertbody_35.nii.gz -qc qc -qc-subject ${file_t2}

# Register T2-w image to PAM50 template
sct_register_to_template -i ${file_t2}.nii.gz -s ${file_t2}_seg.nii.gz -l ${file_t2}_seg_labeled_vertbody_35.nii.gz -c t2 -param step=1,type=seg,algo=centermassrot:step=2,type=seg,algo=syn,slicewise=1,smooth=0,iter=5:step=3,type=im,algo=syn,slicewise=1,smooth=0,iter=3 -qc qc -qc-subject ${file_t2}

# Rename output for clarity
mv anat2template.nii.gz ${file_t2}_reg.nii.gz

4. brought rootlets segmentation from native space to PAM50 space:

# Bring rootlets segmentation from native space to PAM50 space
# Note: Nearest-neighbor interpolation (-x nn) must be used to preserve the level-wise rootlets segmentation values
sct_apply_transfo -i ${file_t2}_label-rootlet_staple.nii.gz -d $SCT_DIR/data/PAM50/template/PAM50_t2.nii.gz -w warp_anat2template.nii.gz -x nn

Legend: red - PAM50 rootlets, blue - native rootlets registered to PAM50

Now, I am working on the label-wise registration between PAM50 rootlets and subject rootlets. TODO - add details.

[jcohenadad]

I also think this should be debated during SCT's course "discussion period"

I've added a "discussion needed" in spinalcordtoolbox/spinalcordtoolbox#4186

[valosekj]

Following up on my previous comment, I am currently working on:

5. label-wise Tz-only registration between PAM50 rootlets (fixed) and subject rootlets (moving).

Since sct_register_multimodal has limited flexibility in controlling its step=0, I am continuing directly with $SCT_DIR/bin/isct_antsRegistration.
The following command combining Affine and BSplineSyN registrations produces promising results!

file_t2=sub-007_ses-headNormal_T2w

$SCT_DIR/bin/isct_antsRegistration --dimensionality 3 --float 0 \
--output [registration2_,${file_t2}_label-rootlet_staple_reg_reg.nii.gz] --interpolation nearestNeighbor --verbose 1 \
--transform Affine[5] --metric MeanSquares[PAM50_t2_label-rootlet.nii.gz,${file_t2}_label-rootlet_staple_reg.nii.gz,1,32] --convergence 20x10 --shrink-factors 2x1 --smoothing-sigmas 0x0mm --shrink-factors 1x1 --restrict-deformation 0x0x1 \
--transform BSplineSyN[0.1,3,0] --metric MeanSquares[PAM50_t2_label-rootlet.nii.gz,${file_t2}_label-rootlet_staple_reg.nii.gz,1,32] --convergence 10x5 --shrink-factors 2x1 --smoothing-sigmas 0x0mm --shrink-factors 1x1 --restrict-deformation 0x0x1

Notice the shift in S-I direction between the blue (before the Tz-only registration) and green (after the Tz-only registration). The green is now much closer to red (PAM50 rootlets)!

Legend: red - PAM50 rootlets, blue - native rootlets registered to PAM50, green - native rootlets registered to PAM50 after the second (Tz-only) registration

[valosekj]

Continuing:

$SCT_DIR/bin/isct_antsRegistration outputted the following files:

registration2_1InverseWarp.nii.gz       # PAM50rootlets --> nativeRootletsInPAM50space -- using this one
registration2_1Warp.nii.gz              # nativeRootletsInPAM50space --> PAM50rootlets 
registration2_0GenericAffine.mat

6. concatenate PAM50rootlets --> nativeRootletsInPAM50space (registration2_1InverseWarp.nii.gz) and nativeRootletsInPAM50space --> nativeRootletsInNativeSpace (warp_template2anat.nii.gz) transformations to obtain a single warping field (PAM50rootlets --> nativeRootletsInNativeSpace):

file_t2=sub-007_ses-headNormal_T2w

sct_concat_transfo -d ${file_t2}.nii.gz -w registration2_1InverseWarp.nii.gz warp_template2anat.nii.gz -o warp_final.nii.gz

7. get spinal levels from Frostell et al. paper (see PR #18):

cd ~/code/PAM50
git pull
git checkout jca/16-spinal-levels

8. bring PAM50 spinal levels (template/PAM50_spinal_levels.nii.gz) from PAM50 template space to subject native space using the concatenated transformation:

sct_apply_transfo -i ~/code/PAM50/template/PAM50_spinal_levels.nii.gz -d ${file_t2}.nii.gz -w warp_final.nii.gz -x nn

Resulting native image and native rootlets with registered PAM50_spinal_levels.nii.gz:

[jcohenadad]

omg this is amazing, great job @valosekj 🏅 👏!!! I could watch that video for hours. Let's make a similar one for SCT course where you go back-and-forth along AP in the coronal plane.

[valosekj]

erratum to my previous point about warping field concatenation:

$SCT_DIR/bin/isct_antsRegistration outputted the following files:

registration2_1InverseWarp.nii.gz       # PAM50rootlets --> nativeRootletsInPAM50space -- using this one
registration2_1Warp.nii.gz              # nativeRootletsInPAM50space --> PAM50rootlets 
registration2_0GenericAffine.mat

I realized that I was not using the affine matrix registration2_0GenericAffine.mat produced by the first step of the $SCT_DIR/bin/isct_antsRegistration command (--transform Affine[5]).
Why? Because I initially thought the affine transform was included within the registration2_1InverseWarp.nii.gz). But apparently, it is not included, and each transform should be treated separately; see here.

Now, when manually including the affine in the sct_concat_transfo command to produce the final warping field (warp_final_including_affine.nii.gz):

# Note that `-winv` is used to invert the affine transform
$ sct_concat_transfo -d ${file_t2}.nii.gz -w registration2_0GenericAffine.mat registration2_1InverseWarp.nii.gz warp_template2anat.nii.gz -o warp_final_including_affine.nii.gz -winv registration2_0GenericAffine.mat

And warping PAM50 spinal levels (template/PAM50_spinal_levels.nii.gz) to subject native space:

$ sct_apply_transfo -i ~/code/PAM50/template/PAM50_spinal_levels.nii.gz -d ${file_t2}.nii.gz -w warp_final_including_affine.nii.gz -x nn -o PAM50_spinal_levels_reg_with_affine.nii.gz

The warped PAM50 levels (PAM50_spinal_levels_reg_with_affine.nii.gz) are now shifted a little bit up:

But it seems to me that the spinal levels are now too high relative to the vertebrae. See, for example, the C4 level (pink color with the cursor). After the shift up, the C4 level is at the level of disc C2/C3. But it should instead be located at the level of vertebra C3; see illustration here.

[jcohenadad]

Why do you include two times registration2_0GenericAffine with -w in the beginning and -win in the end?

$ sct_concat_transfo -d ${file_t2}.nii.gz -w registration2_0GenericAffine.mat registration2_1InverseWarp.nii.gz warp_template2anat.nii.gz -o warp_final_including_affine.nii.gz -winv registration2_0GenericAffine.mat

I would do:

$ sct_concat_transfo -d ${file_t2}.nii.gz -winv registration2_0GenericAffine.mat registration2_1InverseWarp.nii.gz warp_template2anat.nii.gz -o warp_final_including_affine.nii.gz

[valosekj]

Why do you include two times registration2_0GenericAffine with -w in the beginning and -win in the end?

Based on the following sct_concat_transfo.py -win flag help, I understood that the transformation included in the flag -win must also be included in the flag -w.

help='Affine transformation(s) listed in flag -w which should be inverted before being used. Note that this '
     'only concerns affine transformation (not warping fields). If you would like to use an inverse warping '
     'field, then directly input the inverse warping field in flag -w.',

This is also confirmed by the isct_ComposeMultiTransform call when all three transformations are used:

$ sct_concat_transfo -d ${file_t2}.nii.gz -w registration2_0GenericAffine.mat registration2_1InverseWarp.nii.gz warp_template2anat.nii.gz -o warp_final_including_affine.nii.gz -winv registration2_0GenericAffine.mat

...

Parse list of warping fields...

Check file existence...
  OK: sub-007_ses-headNormal_T2w.nii.gz
  OK: registration2_0GenericAffine.mat
  OK: registration2_1InverseWarp.nii.gz
  OK: warp_template2anat.nii.gz
/Users/valosek/code/sct_latest/bin/isct_ComposeMultiTransform 3 warp_final.nii.gz -R sub-007_ses-headNormal_T2w.nii.gz warp_template2anat.nii.gz registration2_1InverseWarp.nii.gz -i registration2_0GenericAffine.mat # in /Users/valosek/data/data.neuro.polymtl.ca/rootlets_seg/register_to_PAM50

---

I would do:

$ sct_concat_transfo -d ${file_t2}.nii.gz -winv registration2_0GenericAffine.mat registration2_1InverseWarp.nii.gz warp_template2anat.nii.gz -o warp_final_including_affine.nii.gz

When trying this command, I got an error about missing -w arg. So I included the flag into the command:

$ sct_concat_transfo -d ${file_t2}.nii.gz -winv registration2_0GenericAffine.mat -w registration2_1InverseWarp.nii.gz warp_template2anat.nii.gz -o warp_final_including_affine_fixed.nii.gz

However, registration2_0GenericAffine.mat is not used at all!

...
Parse list of warping fields...

Check file existence...
  OK: sub-007_ses-headNormal_T2w.nii.gz
  OK: registration2_1InverseWarp.nii.gz
  OK: warp_template2anat.nii.gz
/Users/valosek/code/sct_latest/bin/isct_ComposeMultiTransform 3 warp_final.nii.gz -R sub-007_ses-headNormal_T2w.nii.gz warp_template2anat.nii.gz registration2_1InverseWarp.nii.gz # in /Users/valosek/data/data.neuro.polymtl.ca/rootlets_seg/register_to_PAM50

[RaphaSchl]

Hey Jan, these advances are really exciting !

I also think you have to add registration2_0GenericAffine.mat in both -w and -winv when using sct_concat_transfo, to indicate that you want this matrix inverted when concatenated.
But in this last try here, it seems to not have been used at all because you forgot to add it in -w .

When trying this command, I got an error about missing -w arg. So I included the flag into the command:
$ sct_concat_transfo -d ${file_t2}.nii.gz -winv registration2_0GenericAffine.mat -w registration2_1InverseWarp.nii.gz warp_template2anat.nii.gz -o warp_final_including_affine_fixed.nii.gz
However, registration2_0GenericAffine.mat is not used at all!

Have you tried :
$ sct_concat_transfo -d ${file_t2}.nii.gz -winv registration2_0GenericAffine.mat -w registration2_0GenericAffine.mat registration2_1InverseWarp.nii.gz warp_template2anat.nii.gz -o warp_final_including_affine_fixed.nii.gz

Let me know how this goes !

[valosekj]

Thank you, Raphaelle!

Have you tried :
$ sct_concat_transfo -d ${file_t2}.nii.gz -winv registration2_0GenericAffine.mat -w registration2_0GenericAffine.mat registration2_1InverseWarp.nii.gz warp_template2anat.nii.gz -o warp_final_including_affine_fixed.nii.gz

Yes, I have tried! And indeed, this is the right command to concatenate all transformations. In other words, registration2_0GenericAffine.mat has to be included in both -w and -winv.

For clarity, I put all the commands discussed in this thread into a bash script (see here).

However, I noticed that the label-wise Tz-only registration does not work well in some subjects --> I'm investigating that!