Straightening is not a bijective transformation about spinalcordtoolbox

[mguaypaq]

Bingo! That fixed the issue. I'll make a pull request. I just hope that other code didn't rely on the buggy behaviour.

[mguaypaq]

I spent today working through the straightening code. I did find a few things that really look like bugs, but none of them improved the shearing effect we're seeing on the images here. I'll keep looking.

• In this code, there's a linear algebra mistake: the projection formula used only works when the vectors self.derivatives[indexes] are unit-length vectors. But (empirically) the vectors self.derivatives[indexes] are not unit-length.
  

  spinalcordtoolbox/spinalcordtoolbox/types.py

  Lines 323 to 324 in a7a3f8d

  	def get_projected_coordinates_on_planes(self, coordinates, indexes):
  	return coordinates - multiply(tile(einsum('ij,ij->i', coordinates - self.points[indexes], self.derivatives[indexes]), (3, 1)).transpose(), self.derivatives[indexes])

• In this code, the centerline_straight is initialized with a mix of physical space points, and voxel space derivatives. I'm pretty sure the derivatives should also be given in physical space.
  

  spinalcordtoolbox/spinalcordtoolbox/straightening.py

  Lines 346 to 357 in a7a3f8d

  	coord_straight = np.empty((number_of_points, 3))
  	coord_straight[..., 0] = int(np.round(nx_s / 2))
  	coord_straight[..., 1] = int(np.round(ny_s / 2))
  	coord_straight[..., 2] = np.linspace(0, end_point_coord[2] - start_point_coord[2], number_of_points)
  	coord_phys_straight = image_centerline_straight.transfo_pix2phys(coord_straight)
  	derivs_straight = np.empty((number_of_points, 3))
  	derivs_straight[..., 0] = derivs_straight[..., 1] = 0
  	derivs_straight[..., 2] = 1
  	dx_straight, dy_straight, dz_straight = derivs_straight.T
  	centerline_straight = Centerline(coord_phys_straight[:, 0], coord_phys_straight[:, 1],
  	coord_phys_straight[:, 2],
  	dx_straight, dy_straight, dz_straight)

• This bit of code for straight2curved doesn't match the corresponding code for curved2straight. It should almost certainly be using centerline_straight.get_inverse_plans_coordinates instead of doing manual matrix manipulations.
  

  spinalcordtoolbox/spinalcordtoolbox/straightening.py

  Lines 481 to 483 in a7a3f8d

  	coord_curved2straight = centerline_straight.points[lookup]
  	coord_curved2straight[:, 0:2] += coord_in_planes_curved[:, 0:2]
  	coord_curved2straight[:, 2] += distances_curved

[mguaypaq]

A quick experiment:

Looking at sub-DEV202 and sub-DEV203, which have the same dimensions and very similar (but different) qforms:

• sub-DEV202 looks bijective
• sub-DEV203 shows the problem
• copying the header of sub-DEV202 to sub-DEV203 fixes the problem (looks bijective)
• copying the header of sub-DEV203 to sub-DEV202 introduces the problem

Looking at both versions of sub-DEV202, it looks like the problem shows up when there's a significant forward tilt. I don't know the cause yet, but maybe having the SI axis at a significant angle from vertical creates a problem? I'll look into it.

[mguaypaq]

Ok, I think I found something: the centerline derivatives returned by this call to get_centerline look very strange:

spinalcordtoolbox/spinalcordtoolbox/straightening.py

Lines 152 to 157 in bb479d8

	# 2. extract bspline fitting of the centerline, and its derivatives
	img_ctl = Image('centerline_rpi.nii.gz')
	_, arr_ctl_phys, arr_ctl_der_phys, _ = get_centerline(img_ctl, self.param_centerline,
	verbose=verbose, space="phys")
	centerline = Centerline(*arr_ctl_phys, *arr_ctl_der_phys)
	number_of_points = centerline.number_of_points

Specifically, the centerline derivative at a point should be almost parallel to the difference between two consecutive points. But in this case, it's really not. On some example data, I picked the 100th point out of 297:

>>> _, pts, der, _ = get_centerline(Image('t2_seg.nii.gz'), ParamCenterline(algo_fitting='nurbs'), space="phys")
>>> pts.shape, der.shape
((3, 297), (3, 297))
>>> pts[:, 101] - pts[:, 100]
array([0.7941607 , 0.01081829, 0.42369515])
>>> der[:, 100]
array([-0.01059804,  0.40128474,  0.799999  ])

By comparison, if I remove the space="phys" option, I get almost perfect agreement:

>>> _, pts, der, _ = get_centerline(Image('t2_seg.nii.gz'), ParamCenterline(algo_fitting='nurbs'))
>>> pts.shape, der.shape
((3, 297), (3, 297))
>>> pts[:, 101] - pts[:, 100]
array([-0.01384741,  0.52808393,  1.        ])
>>> der[:, 100]
array([-0.0135655 ,  0.51364441,  1.        ])

So, there may be something seriously wrong with the space="phys" option to get_centerline.

[mguaypaq]

Looking at that pull request, it looks like that bit of code was just moved, not really changed. Also, this would be consistent with "diagonal" qforms working right, and "slanted" qforms giving deformed results.

[joshuacwnewton]

Amazing work, @mguaypaq!!! Thank you for such a thorough investigation of the problem, and for deeply interrogating the implementation of the straightening paper. ♥️

[Kaonashi22]

I guess this is the same issue with this subject; the vertebral bodies don't align well between the UNI image (blue arrow) and the T1 template after registration (green)
.

PAM50_t1_reg.nii.gz
sub-AM273_UNIT1.nii.gz
[sub-AM273_UNIT1_crop.nii.gz](https://github.com/spinal
sub-AM273_UNIT1_seg.nii.gz
cordtoolbox/spinalcordtoolbox/files/15312748/sub-AM273_UNIT1_crop.nii.gz)
batch_processing_sub-AM273.log

[jcohenadad]

none of the curve2straight/straight2curve transform pairs are technically bijective, the effect is just much more pronounced for some subjects

Right. Thank you for clarifying this, @joshuacwnewton. This is indeed what I meant.

[Kaonashi22]

Thank you, hope the issue can be solved!

[Kaonashi22]

Thanks @mguaypaq !

[jcohenadad]

I concur, bravo and thank you @mguaypaq 👏 🙏 !!!

[jcohenadad]

Input resolution of sub-DEV169Sujet03_T2.nii.gz is:

pixdim1		0.781250
pixdim2		0.781250
pixdim3		0.800000

NIfTI header

filename	sub-DEV169Sujet03_T2.nii.gz

sizeof_hdr	348
data_type	INT16
dim0		3
dim1		320
dim2		320
dim3		64
dim4		1
dim5		1
dim6		1
dim7		1
vox_units	mm
time_units	s
datatype	4
nbyper		2
bitpix		16
pixdim0		-1.000000
pixdim1		0.781250
pixdim2		0.781250
pixdim3		0.800000
pixdim4		1.500000
pixdim5		0.000000
pixdim6		0.000000
pixdim7		0.000000
vox_offset	352
cal_max		0.000000
cal_min		0.000000
scl_slope	1.000000
scl_inter	0.000000
phase_dim	2
freq_dim	1
slice_dim	3
slice_name	Unknown
slice_code	0
slice_start	0
slice_end	0
slice_duration	0.000000
toffset		0.000000
intent		Unknown
intent_code	0
intent_name	
intent_p1	0.000000
intent_p2	0.000000
intent_p3	0.000000
qform_name	Scanner Anat
qform_code	1
qto_xyz:1	-0.000000 -0.000000 0.800000 -30.983202 
qto_xyz:2	-0.734600 0.265922 -0.000000 65.615944 
qto_xyz:3	0.265922 0.734600 0.000000 -223.927170 
qto_xyz:4	0.000000 0.000000 0.000000 1.000000 
qform_xorient	Anterior-to-Posterior
qform_yorient	Inferior-to-Superior
qform_zorient	Left-to-Right
sform_name	Scanner Anat
sform_code	1
sto_xyz:1	0.000000 0.000000 0.800001 -30.983202 
sto_xyz:2	-0.734600 0.265922 -0.000000 65.615944 
sto_xyz:3	0.265922 0.734600 0.000000 -223.927170 
sto_xyz:4	0.000000 0.000000 0.000000 1.000000 
sform_xorient	Anterior-to-Posterior
sform_yorient	Inferior-to-Superior
sform_zorient	Left-to-Right
file_type	NIFTI-1+
file_code	1
descrip		TE=1.3e+02;Time=121014.520;phase=1
aux_file

[jcohenadad]

Trying with another (longer FOV) file, from the spine-generic project, results are as expected (ie: very similar original image vs. curve2straight2curve):

NIfTI header

filename	sub-unf02_T2w.nii.gz

sizeof_hdr	348
data_type	FLOAT64
dim0		3
dim1		64
dim2		320
dim3		320
dim4		1
dim5		1
dim6		1
dim7		1
vox_units	mm
time_units	s
datatype	64
nbyper		8
bitpix		64
pixdim0		1.000000
pixdim1		0.799999
pixdim2		0.800000
pixdim3		0.800000
pixdim4		0.000000
pixdim5		0.000000
pixdim6		0.000000
pixdim7		0.000000
vox_offset	352
cal_max		0.000000
cal_min		0.000000
scl_slope	0.000000
scl_inter	0.000000
phase_dim	2
freq_dim	1
slice_dim	3
slice_name	Unknown
slice_code	0
slice_start	0
slice_end	0
slice_duration	0.000000
toffset		0.000000
intent		Unknown
intent_code	0
intent_name	
intent_p1	0.000000
intent_p2	0.000000
intent_p3	0.000000
qform_name	Scanner Anat
qform_code	1
qto_xyz:1	0.799999 -0.000000 0.000000 -24.385565 
qto_xyz:2	0.000000 0.800000 0.000000 -122.899323 
qto_xyz:3	-0.000000 -0.000000 0.800000 -121.175903 
qto_xyz:4	0.000000 0.000000 0.000000 1.000000 
qform_xorient	Left-to-Right
qform_yorient	Posterior-to-Anterior
qform_zorient	Inferior-to-Superior
sform_name	Scanner Anat
sform_code	1
sto_xyz:1	0.799999 0.000000 0.000000 -24.385565 
sto_xyz:2	0.000000 0.800000 0.000000 -122.899323 
sto_xyz:3	0.000000 -0.000000 0.800000 -121.175903 
sto_xyz:4	0.000000 0.000000 0.000000 1.000000 
sform_xorient	Left-to-Right
sform_yorient	Posterior-to-Anterior
sform_zorient	Inferior-to-Superior
file_type	NIFTI-1+
file_code	1
descrip		TE=1.2e+02;Time=131904.995;phase=1
aux_file	

[jcohenadad]

Trying with another image from Lydia: sub-BB277: I notice the same issue:

NIfTI header

filename	sub-BB277_T2.nii.gz

sizeof_hdr	348
data_type	INT16
dim0		3
dim1		320
dim2		320
dim3		64
dim4		1
dim5		1
dim6		1
dim7		1
vox_units	mm
time_units	s
datatype	4
nbyper		2
bitpix		16
pixdim0		-1.000000
pixdim1		0.781250
pixdim2		0.781250
pixdim3		0.800000
pixdim4		1.500000
pixdim5		0.000000
pixdim6		0.000000
pixdim7		0.000000
vox_offset	352
cal_max		0.000000
cal_min		0.000000
scl_slope	1.000000
scl_inter	0.000000
phase_dim	2
freq_dim	1
slice_dim	3
slice_name	Unknown
slice_code	0
slice_start	0
slice_end	0
slice_duration	0.000000
toffset		0.000000
intent		Unknown
intent_code	0
intent_name	
intent_p1	0.000000
intent_p2	0.000000
intent_p3	0.000000
qform_name	Scanner Anat
qform_code	1
qto_xyz:1	0.000000 0.000000 0.800000 -22.308386 
qto_xyz:2	-0.762136 0.171755 0.000000 88.849670 
qto_xyz:3	0.171755 0.762136 -0.000000 -165.046005 
qto_xyz:4	0.000000 0.000000 0.000000 1.000000 
qform_xorient	Anterior-to-Posterior
qform_yorient	Inferior-to-Superior
qform_zorient	Left-to-Right
sform_name	Scanner Anat
sform_code	1
sto_xyz:1	0.000000 0.000000 0.799999 -22.308386 
sto_xyz:2	-0.762136 0.171755 -0.000000 88.849670 
sto_xyz:3	0.171755 0.762136 0.000000 -165.046005 
sto_xyz:4	0.000000 0.000000 0.000000 1.000000 
sform_xorient	Anterior-to-Posterior
sform_yorient	Inferior-to-Superior
sform_zorient	Left-to-Right
file_type	NIFTI-1+
file_code	1
descrip		TE=1.3e+02;Time=84748.520;phase=1
aux_file	

So, is this related to the slightly longer FOV on Lydia's images? Let's try with a bigger FOV from our database

[jcohenadad]

Another image from PAM50 data (larger FOV) no issue observed:

I think this has to do with the dimension and/or the q/sform. Here's the one for this image

Terminal output

filename	sub-amuPAM50005_T1w.nii.gz

sizeof_hdr	348
data_type	INT16
dim0		3
dim1		264
dim2		384
dim3		176
dim4		1
dim5		1
dim6		1
dim7		1
vox_units	mm
time_units	s
datatype	4
nbyper		2
bitpix		16
pixdim0		-1.000000
pixdim1		1.000000
pixdim2		1.000000
pixdim3		1.000000
pixdim4		1.900000
pixdim5		1.000000
pixdim6		1.000000
pixdim7		63136.269531
vox_offset	352
cal_max		0.000000
cal_min		0.000000
scl_slope	1.000000
scl_inter	0.000000
phase_dim	0
freq_dim	0
slice_dim	3
slice_name	sequential_increasing
slice_code	1
slice_start	0
slice_end	175
slice_duration	0.000000
toffset		0.000000
intent		Unknown
intent_code	0
intent_name	
intent_p1	0.000000
intent_p2	0.000000
intent_p3	0.000000
qform_name	Scanner Anat
qform_code	1
qto_xyz:1	0.000000 0.000000 1.000000 -86.773605 
qto_xyz:2	-1.000000 0.000000 -0.000000 81.023537 
qto_xyz:3	0.000000 1.000000 -0.000000 -232.477036 
qto_xyz:4	0.000000 0.000000 0.000000 1.000000 
qform_xorient	Anterior-to-Posterior
qform_yorient	Inferior-to-Superior
qform_zorient	Left-to-Right
sform_name	Scanner Anat
sform_code	1
sto_xyz:1	0.000000 0.000000 1.000000 -86.773605 
sto_xyz:2	-1.000000 0.000000 0.000000 81.023537 
sto_xyz:3	0.000000 1.000000 0.000000 -232.477036 
sto_xyz:4	0.000000 0.000000 0.000000 1.000000 
sform_xorient	Anterior-to-Posterior
sform_yorient	Inferior-to-Superior
sform_zorient	Left-to-Right
file_type	NIFTI-1+
file_code	1
descrip		TE=3.019999981;sec=63136.2700;phaseDir=+
aux_file	

[jcohenadad]

Resampled our internal image to the dimension of the image from Lydia, but it did not reproduce the issue:

sct_resample -i sub-unf02_T2w.nii.gz -mm 0.78x0.78x0.8 
sct_deepseg -i sub-unf02_T2w_r.nii.gz -task seg_sc_contrast_agnostic
sct_straighten_spinalcord -i sub-unf02_T2w_r.nii.gz -s sub-unf02_T2w_r_seg.nii.gz 
sct_apply_transfo -i sub-unf02_T2w_r_straight.nii.gz -d sub-unf02_T2w_r.nii.gz -w warp_straight2curve.nii.gz 

[jcohenadad]

tagging @benjamindeleener who implemented the method-- maybe he has some insights

[joshuacwnewton]

(Quick sanity check: Installed SCT v5.0.0, and this behavior persists. So, I don't think it's been caused by any of the small tweaks/refactors to straightening since 2020.)

---

As another way of visualizing the problem: I tested the following command on the DEV202/DEV203 subjects highlighted by @mguaypaq's #4477 (comment):

sct_concat_transfo -d {input_img} -w warp_curve2straight.nii.gz warp_straight2curve.nii.gz

This produces a concatenated warping field that can be directly applied to the input image to reproduce the "sheared" results. If I open the concatenated warp_final in FSLeyes and set the brightness/contrast to -0.5/0.5 minmax with a divergent colormap, I get:

(Left: DEV203, Right: DEV202)

DEV203 is much more extreme, but DEV202 still shows visible +/- shifts, indicated by the blue and red patterns. So, I reckon that whatever the issue is, the affected subjects are just the extreme versions of the behavior that already exists. In other words, none of the curve2straight/straight2curve transform pairs are technically bijective, the effect is just much more pronounced for some subjects.

(Maybe this is stating something that was already apparent? But, I thought it would be helpful to potentially reframe the problem as "behavior inherent to sct_straighten_spinalcord's algorithm" rather than "bug that is completely separate from the normal behavior". Meaning, there may not be a smoking gun that we can find here.)

[jcohenadad]

(Quick sanity check: Installed SCT v5.0.0, and this behavior persists. So, I don't think it's been caused by any of the small tweaks/refactors to straightening since 2020.)

Great! I was wondering this actually. Thank you for checking

[Kaonashi22]

Hi, any updates on this issue?

[joshuacwnewton]

Hi, any updates on this issue?

We discussed this in our most recent SCT dev meeting (2024-05-15), and we came to the conclusion that debugging this issue would require a deeper understanding of the underlying straightening code (most likely this transform step of the overall algorithm).

Fully understanding the underlying code may take some time, as the person who is looking into this issue (@mguaypaq) is a different person than the person who originally developed this method some years ago (@benjamindeleener).

---

Note

Benjamin did share some thoughts in our internal Slack, which I can copy to this issue since Benjamin is currently locked out of his GitHub account:

@benjamindeleener: To give some insights, straightening does not produce a bijective transformation. The script computes the two transformations by calculating them two times. The issue was known originally (I believe it is mentioned in the manuscript), but the results should be very close in most cases. We assumed it would not be a problem in the future. It is possible to transform the method to produce a bijective transformation but it may require some time, mainly because of edge issues (the problem is an ill-posed problem without regularization when calculating the transformation far from the centerline.

@jcohenadad: for some reasons, some subjects are waaaaaaay less bijective than others. This is documented in the issue, with some examples. In some cases, depending on the q/form it seems (Mathieu did some tests), there is an abnormal shearing that is being applied on the transformation matrix. This is very strange, and i was wondering if you had some ideas about the cause for it

[benjamindeleener]

Thank you for finding those bugs!

Regarding the issue, my first suggestion would be to look into the centerline derivative for instabilities. Depending on the method used to extract / create the centerline, the smoothing method may introduce oscillations that are dependent on the length of the curve (e.g., a polynomial estimation with fixed end points). I believe we have a tool to display the centerline in the three 2D planes as well as the derivatives.

[jcohenadad]

Regarding the issue, my first suggestion would be to look into the centerline derivative for instabilities. Depending on the method used to extract / create the centerline, the smoothing method may introduce oscillations that are dependent on the length of the curve (e.g., a polynomial estimation with fixed end points). I believe we have a tool to display the centerline in the three 2D planes as well as the derivatives.

Thank you for chipping in @benjamindeleener. I've (quickly) looked at the derivatives of the NURBS fitted centerline, and it looks OK to me (ie no oscillations, and value looks reasonable wrt. the centerline coordinates on X and Y), eg: sct-pipeline/spine-park#33 (comment)

It seems like the issue has more to do with a systematic bias (eg shearing parameter) in the generation of the curve2straight warping field, vs. instabilities (ie noise) in the estimate of the derivatives.

[mguaypaq]

Yes, this is definitely not the right way to adjust the derivatives to physical coordinates:

spinalcordtoolbox/spinalcordtoolbox/centerline/core.py

Lines 273 to 281 in bb479d8

	# If 'phys' is specified, adjust centerline coordinates (`Centerline.points`) and
	# derivatives (`Centerline.derivatives`) to physical ("phys") space and native (`im_seg`) orientation.
	if space == 'phys':
	# Transform centerline to physical coordinate system
	arr_ctl = im_seg.transfo_pix2phys([[x_centerline_fit[i], y_centerline_fit[i], z_ref[i]]
	for i in range(len(x_centerline_fit))]).T
	# Adjust derivatives with pixel size
	_, _, _, _, px, py, pz, _ = im_seg.change_orientation(native_orientation).dim
	arr_ctl_der = np.array([x_centerline_deriv * px, y_centerline_deriv * py, np.ones_like(z_ref) * pz])

[joshuacwnewton]

Side note: We had definitely done some refactoring related to space="phys" in recent years (#4010), but I don't think any of our refactoring was functional in nature? I think we just moved parts of the existing code around (_get_centerline())... I hope... So I'm guessing that this derivative issue might have been present in the older code, too.