Update spinal levels files generated by the University of Calgary's Phillips Lab based on the Mendez et al. paper

[sandrinebedard]

Description

Partially fixes spinalcordtoolbox/spinalcordtoolbox#3952.

This PR takes the unprocessed Philips Lab spinal level files (stored in the spinal_levels_PhillipsLab/ folder), modifies and renames them, adds them to the spinal_levels/ folder, then removes the original spinal_levels_PhillipsLab/ folder.

FOR NOW, NOT CONTINUED, see #3 (comment)

Processing

I noticed some differences between the data from the Phillips Lab and the current PAM50/spinal_levels:

• Orientation LPI vs RPI --> (see comment)
• Resolution 1x1x1 vs 0.5x0.5x0.5
• qform/sform
• type float32 vs float62.

However, the matrix size is the same (141, 141, 991), resampling to 0.5mm would double the matrix size. Copying the header from the current PAM50/spinal_levels fixed the problem instead of resampling. I want to confirm this is ok.

Processing steps includes:

• Reorient from LPI to RPI (see comment)
• Change data type from float64 to float32
• Copy header from current PAM50/spinal_levels
• Rename files

Exact processing commands are included in the README.md of this PR.

Old version of processing steps

I included the processing scripts in the attached zip file: process_spinal_levels.zip
(put the python script in the same directory as the bash script )

Here is what I ran for processing:

chmod +x ./process_spinal_levels.sh
./process_spinal_levels.sh "./pam50/Spinal Cord Levels NIfTI"  <PATH_OUT>

Bibliography

• https://www.frontiersin.org/articles/10.3389/fneur.2016.00238/full

[jcohenadad]

Comparing b686796 (_new) and f8cfc4b (master):

C1

C2

C3

C4

Few comments:

• Is that normal that C1 is completely blank?

• Is that normal that in the new version, values seem almost binary, ie: the same value is spread out along the S-I axis, without "smoothness" to it. Eg: zoomed version of C2_new (same value across ~20 voxels). Maybe we should inquired the Phillips group about that?

  
  With the histogram of C2_new:
  

• The cord seg and the spinal levels seem to be perfectly matching (good), although I have not verified all levels. @sandrinebedard have you verified that? See example of good matching below:

  

• The world coordinate between the source (Phillips lab) and PAM50, for the same voxel location, is drastically different. I guess this this is the reason for the qform copy (see example below).

  

  

[sandrinebedard]

• The cord seg and the spinal levels seem to be perfectly matching (good), although I have not verified all levels. @sandrinebedard have you verified that? See example of good matching below:

I'll verify them all!

• Is that normal that C1 is completely blank?

They did not include any C1 files

• Is that normal that in the new version, values seem almost binary, ie: the same value is spread out along the S-I axis, without "smoothness" to it. Eg: zoomed version of C2_new (same value across ~20 voxels). Maybe we should inquired the Phillips group about that?

I'll also check the https://github.com/PhillipsLab/pam50/tree/main/DREZ%20NIfTI to see if the Spinal Cord levels was teh right folder to take here.

Is that normal that the resampling is commented out?

• The world coordinate between the source (Phillips lab) and PAM50, for the same voxel location, is drastically different. I guess this this is the reason for the qform copy (see example below).

Regarding the differences in world coordinates and resampling, I wasn't sure how to proceed.

The dimensions and resolution of:

• new spinal levels files are : 141x141x991 and 1x1x1
• Old spinal levels files: 141x141x991 and 0.5x0.5x0.5

If we do resampling to 0.5 mm isotropic, the matrix size doubles, wich doesn't match anymore the PAM50 space.

I thought of doing registration with identity, however, since the world coordinates differ, it dosn't bring the spinal levels in the PAM50 space. How I overcame this issue is by copying the image header (for qform) and this also solves the resolution problem. I am not sure this is the best way to overcome this however...

[jcohenadad]

They did not include any C1 files

Hum, I am a bit reluctant to include a file labeled "C1 spinal levels" but without any label inside-- this is quite confusing for users. Maybe we should just get rid of it?

I'll also check the https://github.com/PhillipsLab/pam50/tree/main/DREZ%20NIfTI to see if the Spinal Cord levels was teh right folder to take here.

I think we should include the Phillips group right away in the conversation to avoid any confusion. I'll take care of that.

How I overcame this issue is by copying the image header (for qform) and this also solves the resolution problem.

Only changing the image header should not affect the physical dimension of the voxel, unless (i) you also did a resampling afterwards or (ii) the "1mm iso" on the native image was wrong, and the true voxel size it was in fact 0.5mm iso.

[sandrinebedard]

I updated the files removing the reorientation to RPI: ce64bdf
process_spinal_levelsv2.zip

[joshuacwnewton]

I updated the files removing the reorientation to RPI: ce64bdf (process_spinal_levelsv2.zip)

I have taken the above scripts and put them into a README.md file. But, I slightly modified the scripts to remove the dependency on Python. I also modified the script so that it can be run directly on commit e854bba.

(The differences between "process_spinal_levelsv2.zip" and the README.md file are shown here.)

[joshuacwnewton]

I've tidied up this PR a little by hiding comments that have already been resolved.

I've also opened discussion threads for the "open questions" that have yet to be resolved, and that we need to address before merging this PR.

(Several of these concerns were raised to the Phillips Lab via email, and we are waiting for a reply.)

---

Hi Aaron,

I hope this email finds you well.

We are in the process of including the spinal levels from your group into SCT's next release. The working PR is here: #3

We have a few questions, notably:

• are these files the latest version and the ones we should for the PAM50? https://github.com/PhillipsLab/pam50/tree/main/Spinal%20Cord%20Levels%20NIfTI
• Also see my questions here: #3 (comment), notably about the 'smoothness' of the files
• Is there a reason that there is no C1 label?

You can give us some feedback directly in the PR, so everything is centralized.

Many thanks,
Julien

---

To keep the PR tidy, please continue all discussions in these threads, so that we can mark them as "resolved" once we determine the answers. Thank you! 🎉

[joshuacwnewton]

Blank C1 file

Is that normal that C1 is completely blank?

They did not include any C1 files

Hum, I am a bit reluctant to include a file labeled "C1 spinal levels" but without any label inside-- this is quite confusing for users. Maybe we should just get rid of it?

[joshuacwnewton]

The response from Aaron Phillips was as follows:

Q: Is there a reason that there is no C1 label?
A: There were not any anatomical measurements for C1 (in the Mendez paper). As a result, we did not include this in our updated estimations of spinal levels.

This partially answers the question, however we still need to decide whether or not we should remove the blank C1 file.

[jcohenadad]

This partially answers the question, however we still need to decide whether or not we should remove the file.

we definitely should. There is absolutely no point in having a 'label' file that contains only zeros. It is just confusing for everyone.

[joshuacwnewton]

Personally, I'm in favor of removing the C1 file, as it would it help to avoid confusion such as this user on the forum: https://forum.spinalcordmri.org/t/c1-vertebral-level-labelling/1002/4?u=joshuacwnewton

However, @jcohenadad has also promised that we would add coverage of the C1 vertebral level in https://github.com/spinalcordtoolbox/spinalcordtoolbox/issues/3997.

So, is the idea that we would remove it just for this release, then re-add it once we have proper spinal level data? 🤔

[jcohenadad]

wrong promise then 😅 , let's get rid of the file

[joshuacwnewton]

Sharp transition between level values (with no gradient/smoothness)

Is that normal that in the new version, values seem almost binary, ie: the same value is spread out along the S-I axis, without "smoothness" to it. Eg: zoomed version of C2_new (same value across ~20 voxels). Maybe we should inquired the Phillips group about that?

[joshuacwnewton]

The response from Aaron Phillips is as follows:

Email exchange

Julien Cohen-Adad: Q: Is it normal that in the new version, values seem almost binary, ie: the same value is spread out along the S-I axis, without "smoothness" to it.

Aaron Phillips: A: The dataset that we used to estimate both the spinal levels and the corresponding DREZ provided only the mean of each segment (C2-L5) and the standard deviation, from 9 cadavers. Our probability estimate assumed that the mean of each segment had a 100% probability, and we added a normal distribution on both the rostral and caudal ends of each segment to account for variability using the standard deviation of measurement for each spinal level / DREZ (from Mendez paper). This was logical to us, but maybe you have a different thought? Popped a figure in below to help explain.

[jcohenadad]

What I don't understand from this diagram, is that the "mean" is not a single value in the absissa, but a range. I'm not sure how to interpret that... i'll answer in the email

[joshuacwnewton]

The email has been answered, and there has been some back and forth discussion as follows:

Email exchange

Julien Cohen-Adad: What I don't understand from this diagram, is that the "mean" is not a single value, but a range. I'm not sure how to interpret that... To me, if the location of subject 1 is Z=[3,5], subject 2 [4,6] and subject 3 [5,7], then the "average" should be (one line corresponds to a value of Z):

1:0
2:0
3: 0.33
4: 0.66
5: 1
6: 0.66
7: 0.33
8: 0
etc.

So I would not expect to see such a "block" of ones.

Can you please provide us with more details about how this mean spinal segment was obtained?

Aaron Phillips: The dataset that we used to estimate probabilities for each spinal cord level (from the Mendez paper) contained only the population mean and standard deviation across all 9 cadavers from which measurements were taken. As a result, it was not possible to create probability estimates for each layer (1 rostral-caudal voxel level) in the PAM50 model.

We used the total length of the spinal cord between C2 and L5 in the PAM50 model, and with the measurements from the Mendez paper to determine the midpoint of each spinal cord level based on the length of each segment relative to the total length. With the data that we had, we moved forward by estimating the probability to be 1 for each spinal cord level, as it was given as a population mean, and added uncertainty based on the standard deviations. We thought that this approach was feasible given the limitations of the data that we used, relying on a mean length (±st-dev) for a population.

Happy to clarify further

Julien Cohen-Adad: Thank you for your answers.

I think I finally understand what you did, and why I was confused earlier. Our approach, with the PAM50, was to represent the mid-point of the nerve rootlets, whereas in your case (please correct me if I'm wrong), you represent the complete length of a spinal segment (ie: the entire "fanning" of rootlets entry point, which is not a single point in space, but a "range"). > So, if I understand correctly, you did the following:

• compute the length of the PAM50 cord between C2 and L5 mid-point, let's call it "cord(PAM50)". Related question: Did you compute the total length between C2 and L5 mid-points, or between C2 (most rostral of nerve root entry) and L5 (most caudal of nerve root entry). I presume you did the former, because the PAM50 does not have information about most rostral/caudal entry of nerve rootlets.
• got the length of the Mendez cords between C2 and L5 mid-point, averaged across subjects let's call it cord(Mendez)
• got the length of each spinal segment (ie: from the most rostral to the most caudal nerve rootlet), averaged across subjects. Let's call it length_segment[i](Mendez) (i stands for a segment, eg: L5).
• scale the length of each spinal segment to match the PAM50, using the following equation:
  • length_segment[i](PAM50) = length_segment[i](Mendez) * cord(PAM50) / cord(Mendez)
• scale the position of the mid-points of each spinal level from the average Mendez cord, to the PAM50 cord, using the same normalization strategy as in the previous point.
• scale the STD of each spinal level from the average Mendez cord, to the PAM50 cord, using the same normalization strategy. By the way: what STD are we talking about here? Is that the STD that represents the variation of the mid-point along the rostro-caudal axis?
• build the PAM50 spinal level template, by centering each segment in the mid-points, and by adding, at the rostral and caudal edge of the segment, a half-gaussian function with height of one and STD of value calculated at the previous point.

Also, if my description above is correct, then I do not understand the discrepancy between the figure of the article, and the observed spinal segments. More specifically, the levels in the figure appear to vary more like a Gaussian along the rostro-caudal axis, while the generated spinal segments look more like a segment of "ones", with short smoothness at the edges:

Could you share the code and the data that were used to generate the spinal level atlas? That way we will be able to understand exactly what was done. In the future, if users ask us how these levels were calculated, we will be able to point to the exact code, which is better practice for transparency and reproducibility.

About the data: Did you contact Mendez et al., or did you use the values they report in their article? I assume they report all results in the supplementary material, however the link to the supplementary material does not work in the PDF of the article. They also say that the suppl. material is available on the mayoclinic website (see screenshot below), but when going to this website and searching for their article, I do not see a link to the suppl. material. Can you please send it to us if you have it?

We are now waiting on a response from Aaron to confirm the methods used, and to potentially supply us with the raw data + processing code. Until then, I'm not sure if there is much we can do to proceed.

[sandrinebedard]

We decided not to pursue the inclusion of the Mendez et al. paper metrics to update the spinal level in the PAM50 template . Here is a summary of what was done with the related issues:

• The updated version from Phillips Lab had some limitations related to linear scaling : add sct_dmri_moco spinalcordtoolbox#10
• We tried using the distance intervertebral foramen to rostral/caudal rootlets at dorsal entry to identify the spinal levels. We manually labeled the intervertebral foramen in the PAM50. However, the scaling was wrong when brining the measures in the PAM50 template : clean sct_compute_mtr spinalcordtoolbox#12
• We tested using a ratio PAM50-Mendez al. measures of the intervertebral foramen-rootlets: We manually labeled the dorsal rostral and caudal rootlets entry and the intervertebral foramen in the PAM50 template. We computed a ratio between the distance formen-rostral/caudal rootlets in the Mendez et al. paper vs the PAM50. The ratio was not very stable, already at C8 the estimation was bad.: Intervertebral foramen to caudal and rostral rootlets to update spinal levels: limitations #12 (comment)