Unable to warp to SC template

barryrl · February 25, 2021, 11:32am

I ran my updated scripts on 32 datasets and they worked well on 31. The SCT code was, however, unable to register one of the datasets to the template. This particular MR image looks good and the spinal cord is segmented properly, so I am not sure why it is unable to warp to the template. The only user interaction is manually identifying C2/C3; I re-identified this point and got the same output errors (pasted below). Do you have an idea of what might be happening?

Thanks!

-Rob

Warping field generated: tmp.curve2straight.nii.gz
Warping field generated: tmp.straight2curve.nii.gz
Apply transformation to input image…
e[94m/Users/rob/spinalcordtoolbox/bin/isct_antsApplyTransforms -d 3 -r tmp.centerline_pad_crop.nii.gz -i data.nii -o tmp.anat_rigid_warp.nii.gz -t tmp.curve2straight.nii.gz -n ‘BSpline[3]’ # in /private/var/folders/c_/x6z5ffcj2wbc0h68_g1h0l_h0000gn/T/sct-20210225055438.266275-straighten_spinalcord-shrsqloge[0m
Generate output files…
File created: ./warp_curve2straight.nii.gz
File created: ./warp_straight2curve.nii.gz
e[94mcp /var/folders/c_/x6z5ffcj2wbc0h68_g1h0l_h0000gn/T/sct-20210225055438.266275-straighten_spinalcord-shrsqlog/tmp.anat_rigid_warp.nii.gz ./straight_ref.nii.gze[0m
File created: ./data_straight.nii
Remove temporary files…
e[94mrm -rf /var/folders/c_/x6z5ffcj2wbc0h68_g1h0l_h0000gn/T/sct-20210225055438.266275-straighten_spinalcord-shrsqloge[0m
e[0m
Finished! Elapsed time: 366 se[0m
e[0m
Resample to 0.5mm isotropic…e[0m
e[94msct_resample -i data_straight.nii -mm 0.5x0.5x0.5 -x linear -o data_straightr.nii # in /private/var/folders/c_/x6z5ffcj2wbc0h68_g1h0l_h0000gn/T/sct-20210225055434.496959-label_vertebrae-ow9z9cvqe[0m
e[0m
Apply straightening to segmentation…e[0m
e[94m/Users/rob/spinalcordtoolbox/bin/isct_antsApplyTransforms -d 3 -i segmentation.nii -r data_straightr.nii -t warp_curve2straight.nii.gz -o segmentation_straight.nii -n Linear # in /private/var/folders/c_/x6z5ffcj2wbc0h68_g1h0l_h0000gn/T/sct-20210225055434.496959-label_vertebrae-ow9z9cvqe[0m
e[33mFile /private/var/folders/c_/x6z5ffcj2wbc0h68_g1h0l_h0000gn/T/sct-20210225055434.496959-label_vertebrae-ow9z9cvq/segmentation_straight.nii already exists. Will overwrite it.e[0m
e[0m
Create label to identify disc…e[0m
e[33mFile /var/folders/c_/x6z5ffcj2wbc0h68_g1h0l_h0000gn/T/sct-20210225055434.496959-label_vertebrae-ow9z9cvq/labelz.nii.gz already exists. Will overwrite it.e[0m
e[0m
And apply straightening to label…e[0m
e[94m/Users/rob/spinalcordtoolbox/bin/isct_antsApplyTransforms -d 3 -i labelz.nii.gz -r data_straightr.nii -t warp_curve2straight.nii.gz -o labelz_straight.nii.gz -n NearestNeighbor # in /private/var/folders/c_/x6z5ffcj2wbc0h68_g1h0l_h0000gn/T/sct-20210225055434.496959-label_vertebrae-ow9z9cvqe[0m
e[0m
Get z and disc values from straight label…e[0m
e[0m… [365, 3]e[0m
Look for template…
Path template: /Users/rob/SCT/data/PAM50
Open template and vertebral levels…
Disc values from template: [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19]
Z-values for each disc: [963, 938, 907, 870, 833, 800, 769, 735, 692, 646, 600, 551, 500, 449, 396, 342, 289, 231, 168, 104, 79]
Distances between discs (in voxel): [25.0, 31.0, 37.0, 37.0, 33.0, 31.0, 34.0, 43.0, 46.0, 46.0, 49.0, 51.0, 51.0, 53.0, 54.0, 53.0, 58.0, 63.0, 64.0, 25.0]
Detect intervertebral discs…
Current disc: 2 (z=365). Direction: superior
… Peak found: z=-1 (correlation = 0.5806519076174395)
Current disc: 1 (z=395). Direction: superior
… correcting factor: 1.0
… Switching to inferior direction.
Current disc: 3 (z=328). Direction: inferior
… Peak found: z=-1 (correlation = 0.6659052966517715)
… correcting factor: 1.0
Current disc: 4 (z=290). Direction: inferior
… Peak found: z=1 (correlation = 0.6105950224899827)
… correcting factor: 0.990990990990991
Current disc: 5 (z=258). Direction: inferior
… Peak found: z=0 (correlation = 0.44450624761912827)
… correcting factor: 0.9932432432432432
Current disc: 6 (z=227). Direction: inferior
… Peak found: z=-2 (correlation = 0.5430550901155323)
… correcting factor: 1.007497820401046
Current disc: 7 (z=191). Direction: inferior
… Peak found: z=2 (correlation = 0.43807046883648043)
… correcting factor: 0.996444262098911
Current disc: 8 (z=150). Direction: inferior
… Peak found: z=1 (correlation = 0.554691318759119)
… correcting factor: 0.9936299655199968
Current disc: 9 (z=105). Direction: inferior
… Peak found: z=1 (correlation = 0.4752640839639358)
… correcting factor: 0.9917088285256495
Current disc: 10 (z=60). Direction: inferior
… Peak found: z=-1 (correlation = 0.4791790829143412)
… correcting factor: 0.9950455287377754
Current disc: 11 (z=10). Direction: inferior
… Peak found: z=5 (correlation = 0.5583133338704257)
… correcting factor: 0.9853368942313449
Adding top disc based on adjusted template distance: #0
… approximate distance: 25
Traceback (most recent call last):
File “/Users/rob/spinalcordtoolbox/spinalcordtoolbox/scripts/sct_label_vertebrae.py”, line 491, in
main(sys.argv[1:])
File “/Users/rob/spinalcordtoolbox/spinalcordtoolbox/scripts/sct_label_vertebrae.py”, line 425, in main
verbose=verbose, path_template=path_template, path_output=path_output, scale_dist=scale_dist)
File “/Users/rob/spinalcordtoolbox/spinalcordtoolbox/vertebrae/core.py”, line 251, in vertebral_detection
label_discs(fname_seg, list_disc_z, list_disc_value, verbose=verbose)
File “/Users/rob/spinalcordtoolbox/spinalcordtoolbox/vertebrae/core.py”, line 522, in label_discs
cx, cy = [int(x) for x in np.round(center_of_mass(slices)).tolist()]
File “/Users/rob/spinalcordtoolbox/spinalcordtoolbox/vertebrae/core.py”, line 259, in center_of_mass
raise ValueError(“Array has no mass”)
ValueError: Array has no mass

jcohenadad · February 25, 2021, 1:32pm

Hey Rob, I suspect the segmentation has an empty slice at the slice where the center of mass is calculated. Would you mind send me the image + segmentation and the exact syntax + SCT version you used? (next time if you include the syntax in the copy/paste I’ll be able to get both).
Cheers

barryrl · February 25, 2021, 2:00pm

Thanks, @jcohenadad. I’m using SCT git-master-e408b6589689ce2685d2a2cb83901db5b5aed017. I just shared the image, segmentation and syntax with you. -Rob

jcohenadad · February 26, 2021, 2:22am

Hi Rob,

This is the syntax you sent me:

sct_resample -i t2.nii -mm 0.5x0.5x0.5 -o t2r.nii
sct_label_utils -i t2r.nii -create-viewer 3 -o labelc2c3.nii.gz
sct_deepseg_sc -i t2r.nii -c t2 -kernel 2d
sct_label_vertebrae -i t2r.nii -s t2_seg.nii -c t2 -initlabel labelc2c3.nii.gz -t /Users/rob/SCT/data/PAM50
sct_label_utils -i t2_seg_labeled.nii -vert-body 2,5 -o labels_vert.nii.gz
sct_register_to_template -i t2r.nii -s t2_seg.nii -l labels_vert.nii.gz -c t2 -qc qc -t /Users/rob/SCT/data/PAM50

Now, here is the issue: the call to sct_deepseg_sc uses t2r.nii as input, meaning the output will be called t2r_seg.nii, however the subsequent code uses t2_seg.nii which should not exist. So, I suspect you created t2_seg.nii previously, but did not include the syntax in the list of commands above. Anyway, the problem is that you are using sct_label_vertebrae with two images which do not have the same resolutions (t2r.nii is 0.5mm iso while t2_seg.nii is 0.9mm iso).

Another thing: given the T2 data are 0.9mm isotropic resolution (yay! ), there is no need to resample them to 0.5mm iso. It just makes the manual correction more painful and the straightening waaaaay longer.

So: the solution is to simply run:

sct_label_vertebrae -i t2.nii -s t2_seg.nii -c t2 -initlabel labelc2c3.nii.gz

However, this is where things get complicated: This syntax will fail because after we recently introduced this fix, images with a different qform/sform will not be processed and instead the user is prompted to make sure these match. In your case, I visually checked and the qform/sform are almost equal, but not exactly (maybe a precision issue when running FreeSurfer tools):

>>> nib.get_qform()
array([[-6.13306584e-08, -3.25428327e-03, -8.99994553e-01,
         3.68009071e+01],
       [-9.32291687e-01, -6.13306584e-08,  5.94133867e-08,
         2.14219025e+02],
       [ 6.15451159e-08, -9.32286007e-01,  3.14156514e-03,
         1.22887390e+02],
       [ 0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
         1.00000000e+00]])
>>> nib.get_sform()
array([[-6.67520835e-13, -3.25425435e-03, -8.99994552e-01,
         3.68009071e+01],
       [-9.32291687e-01,  1.91215821e-10, -0.00000000e+00,
         2.14219025e+02],
       [-1.91214891e-10, -9.32286024e-01,  3.14158620e-03,
         1.22887390e+02],
       [ 0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
         1.00000000e+00]])

You might want to notify the MGH folks about this. Having different s/qform is cause for trouble in some software.

I’ve open an issue to increase the tolerance when checking for the s/qform so in the future it won’t be a problem in the future.

For now, fortunately there is an easy workaround: copy the qform into the sform:

sct_image -i t2.nii -set-sform-to-qform -o t2_qform.nii

And then run sct_label_vertebral:

sct_label_vertebrae -i t2_qform.nii -s t2_seg.nii -c t2 -qc qc -initlabel labelc2c3.nii.gz

However, this now creates another issue, this time due to the fact that the segmentation has discontinuities:

I’ve opened another issue to address this problem in the future (suggestions welcome), but in the meantime, a workaround is to mask out the “problematic” part of the segmentation:

# Zeroing everything above y=288
sct_crop_image -i t2_seg.nii -ymax 288 -b 0 -o t2_seg_cropped.nii.gz

Now this finally works:

sct_label_vertebrae -i t2_qform.nii -s t2_seg_cropped.nii.gz -c t2 -qc qc -initlabel labelc2c3.nii.gz

Another thing: why did you add “-t /Users/rob/SCT/data/PAM50” in your syntax? Are you using a modified version of the PAM50 template? If not, I recommend to leave it to the default, to make sure to use the PAM50 template that corresponds to the version of the sct_label_vertebrae function (in case you have multiple local installations of SCT).

barryrl · February 26, 2021, 2:32pm

Good morning, @jcohenadad. Thank you for the detailed investigation – I certainly did not expect this one dataset to open Pandora’s box! You had several comments or questions so I’ll do my best to address each of them.

My workflow executes all SCT functions from within the Matlab environment so I can call wrapped SCT functions that perform additional housekeeping (i.e., backing up certain files, copying files, deleting files, etc.). In the case of sct_resample and sct_deepseg_sc, there is indeed UNIX code that I did not include that renames the outputted segmented file as t2_seg.nii. My apologies for any confusion, but I omitted the extra UNIX code in an attempt to avoid confusion. (Clearly that didn’t work as intended. )
The above code was written because you previously recommended that I resample the anatomical image to be 0.5 mm isotropic. Some of the datasets are 0.9 mm isotropic (yay!) but some are lower-res (acquired years ago by other groups), so I thought that it would be best to resample all datasets to 0.5 mm to keep the preprocessing consistent across datasets and projects.
I am using the default PAM50 template. I added -t /Users/rob/SCT/data/PAM50 to the code because /Users/rob/SCT is a symbolic link that points to the current version of the SCT. Therefore, when a new SCT version comes out or if I have to revert (e.g., when sct_concat_transfo was removed from v5.0.1), I just update the symbolic link.
Your observations about qform and sform are interesting. Indeed, they are almost equal. This particular dataset was acquired in Bay 7 more than four years ago, and at least one scanner upgrade has happened since then. The DICOMs are processed using FreeSurfer’s mri_convert so I will keep an eye on this in the future. The discrepancies appear to occur starting in the 8th or 9th decimal place, so perhaps it is due to rounding when converting between data types? I don’t know. Interestingly, the SCT processed this dataset successfully last year but failed this week, so I was wondering if something had changed in the code in the last ~6 months.
I have noticed discontinuities in the spinal cord segmentation from time to time and have been wondering what can be done (or that I can do) to mitigate or obviate such discontinuities without having to manually edit. I recall trying to use a 3D kernel quite some time ago, but it produced an error and so I have stayed with the -kernel 2d option. Presumably sct_deepseg_sc detects when discontinuities occur, so can it be configured to interpolate between discontinuities? For example, fitting a spline to the existing two segments might work well to connect them. There would, of course, be some interpolation error, but the error couldn’t be more than missing sections of the spinal cord entirely. Unless there is major trauma or a malformation, the code could take its “best guess” at interpolating across regions with lower SNR or B0 inhomogeneities.

Hope these comments help - always happy to chat further. And happy Friday!

-Rob

jcohenadad · February 26, 2021, 2:59pm

Hi Rob!

My workflow executes all SCT functions from within the Matlab environment so I can call wrapped SCT functions that perform additional housekeeping (i.e., backing up certain files, copying files, deleting files, etc.). In the case of sct_resample and sct_deepseg_sc, there is indeed UNIX code that I did not include that renames the outputted segmented file as t2_seg.nii. My apologies for any confusion, but I omitted the extra UNIX code in an attempt to avoid confusion. (Clearly that didn’t work as intended. )

no worries-- i assumed so

The above code was written because you previously recommended that I resample the anatomical image to be 0.5 mm isotropic. Some of the datasets are 0.9 mm isotropic (yay!) but some are lower-res (acquired years ago by other groups), so I thought that it would be best to resample all datasets to 0.5 mm to keep the preprocessing consistent across datasets and projects.

ah, if some dataset are sagittal with thick slices then yes, resampling is advised. And if the dataset is heterogeneous, then yes, systematic resampling is the way to go, as we did for spine-generic. You might want to consider lower res though (for the arguments i mentioned: painful to manually correct and longer processing).

I am using the default PAM50 template. I added -t /Users/rob/SCT/data/PAM50 to the code because /Users/rob/SCT is a symbolic link that points to the current version of the SCT. Therefore, when a new SCT version comes out or if I have to revert (e.g., when sct_concat_transfo was removed from v5.0.1), I just update the symbolic link.

Your observations about qform and sform are interesting. Indeed, they are almost equal. This particular dataset was acquired in Bay 7 more than four years ago, and at least one scanner upgrade has happened since then. The DICOMs are processed using FreeSurfer’s mri_convert so I will keep an eye on this in the future. The discrepancies appear to occur starting in the 8th or 9th decimal place, so perhaps it is due to rounding when converting between data types? I don’t know. Interestingly, the SCT processed this dataset successfully last year but failed this week, so I was wondering if something had changed in the code in the last ~6 months.

yup! here

I have noticed discontinuities in the spinal cord segmentation from time to time and have been wondering what can be done (or that I can do) to mitigate or obviate such discontinuities without having to manually edit. I recall trying to use a 3D kernel quite some time ago, but it produced an error and so I have stayed with the -kernel 2d option. Presumably sct_deepseg_sc detects when discontinuities occur, so can it be configured to interpolate between discontinuities? For example, fitting a spline to the existing two segments might work well to connect them. There would, of course, be some interpolation error, but the error couldn’t be more than missing sections of the spinal cord entirely. Unless there is major trauma or a malformation, the code could take its “best guess” at interpolating across regions with lower SNR or B0 inhomogeneities.

the suggested fix should address the problem in the future.

thanks for your feedback!

barryrl · March 19, 2021, 7:59pm

Hi @jcohenadad - just a quick note to confirm that this one problematic dataset was successfully processed with the SCT version downloaded from GitHub today. Thanks for your help! -Rob