Coder Social home page Coder Social logo

Comments (25)

mguaypaq avatar mguaypaq commented on June 26, 2024 6

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

from spinalcordtoolbox.

mguaypaq avatar mguaypaq commented on June 26, 2024 4

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.

    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.

    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.

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

from spinalcordtoolbox.

mguaypaq avatar mguaypaq commented on June 26, 2024 2

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.

from spinalcordtoolbox.

mguaypaq avatar mguaypaq commented on June 26, 2024 2

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

# 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.

from spinalcordtoolbox.

mguaypaq avatar mguaypaq commented on June 26, 2024 2

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.

from spinalcordtoolbox.

joshuacwnewton avatar joshuacwnewton commented on June 26, 2024 2

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

from spinalcordtoolbox.

Kaonashi22 avatar Kaonashi22 commented on June 26, 2024 1

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)
image.

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

from spinalcordtoolbox.

jcohenadad avatar jcohenadad commented on June 26, 2024 1

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.

from spinalcordtoolbox.

Kaonashi22 avatar Kaonashi22 commented on June 26, 2024 1

Thank you, hope the issue can be solved!

from spinalcordtoolbox.

Kaonashi22 avatar Kaonashi22 commented on June 26, 2024 1

Thanks @mguaypaq !

from spinalcordtoolbox.

jcohenadad avatar jcohenadad commented on June 26, 2024 1

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

from spinalcordtoolbox.

jcohenadad avatar jcohenadad commented on June 26, 2024

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

from spinalcordtoolbox.

jcohenadad avatar jcohenadad commented on June 26, 2024

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

ezgif-6-122a5a5205

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	

from spinalcordtoolbox.

jcohenadad avatar jcohenadad commented on June 26, 2024

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

ezgif-6-a8a7d3cef7

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

from spinalcordtoolbox.

jcohenadad avatar jcohenadad commented on June 26, 2024

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

ezgif-6-3c55f85f29

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	

from spinalcordtoolbox.

jcohenadad avatar jcohenadad commented on June 26, 2024

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 

ezgif-2-99836d6db0

from spinalcordtoolbox.

jcohenadad avatar jcohenadad commented on June 26, 2024

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

from spinalcordtoolbox.

joshuacwnewton avatar joshuacwnewton commented on June 26, 2024

(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:

image

(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.)

from spinalcordtoolbox.

jcohenadad avatar jcohenadad commented on June 26, 2024

(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

from spinalcordtoolbox.

Kaonashi22 avatar Kaonashi22 commented on June 26, 2024

Hi, any updates on this issue?

from spinalcordtoolbox.

joshuacwnewton avatar joshuacwnewton commented on June 26, 2024

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

from spinalcordtoolbox.

benjamindeleener avatar benjamindeleener commented on June 26, 2024

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.

from spinalcordtoolbox.

jcohenadad avatar jcohenadad commented on June 26, 2024

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.

from spinalcordtoolbox.

mguaypaq avatar mguaypaq commented on June 26, 2024

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

# 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])

from spinalcordtoolbox.

joshuacwnewton avatar joshuacwnewton commented on June 26, 2024

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.

from spinalcordtoolbox.

Related Issues (20)

Recommend Projects

  • React photo React

    A declarative, efficient, and flexible JavaScript library for building user interfaces.

  • Vue.js photo Vue.js

    🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.

  • Typescript photo Typescript

    TypeScript is a superset of JavaScript that compiles to clean JavaScript output.

  • TensorFlow photo TensorFlow

    An Open Source Machine Learning Framework for Everyone

  • Django photo Django

    The Web framework for perfectionists with deadlines.

  • D3 photo D3

    Bring data to life with SVG, Canvas and HTML. 📊📈🎉

Recommend Topics

  • javascript

    JavaScript (JS) is a lightweight interpreted programming language with first-class functions.

  • web

    Some thing interesting about web. New door for the world.

  • server

    A server is a program made to process requests and deliver data to clients.

  • Machine learning

    Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.

  • Game

    Some thing interesting about game, make everyone happy.

Recommend Org

  • Facebook photo Facebook

    We are working to build community through open source technology. NB: members must have two-factor auth.

  • Microsoft photo Microsoft

    Open source projects and samples from Microsoft.

  • Google photo Google

    Google ❤️ Open Source for everyone.

  • D3 photo D3

    Data-Driven Documents codes.