Allow polynomial function to dummy_segmentation function

[valosekj]

This PR addresses creation of synthetic 4D interleaved data discusses here - #2664 (comment).

So far, I focused on this comment (#2664 (comment)) and allowed polynomial function to dummy_segmentation function in this commit - d561807. I took inspiration from dummy_centerline function.

Is this implementation correct? I can produce following data with this implementation

(code used for creation of this data:)

from spinalcordtoolbox.testing.create_test_data import dummy_segmentation_4d
DEBUG = True  # Set to True to save images
seg_4D = dummy_segmentation_4d(vol_num=10,create_bvecs=True, size_arr=(20, 20, 30), pixdim=(1,1,1),radius_RL=5.0, radius_AP=3.0, shape='ellipse', interleaved=True, debug=DEBUG)

[jcohenadad]

I find it counterintuitive to have the polynomial shape only under interleaved. A clearer approach, would be to allow polynomial shape regardless what interleaved is, and in order to keep the same default behaviour, polynomial order could be set to one by default.

also note that i am starting to see a lot of similarity with dummy_centerline(). At the end, a segmentation is just a centerline with a shape around it (ellipse or rectangle). Would it make sense to merge the two functions? Maybe this is a bad idea…

[valosekj]

A clearer approach, would be to allow polynomial shape regardless what interleaved is, and in order to keep the same default behavior, polynomial order could be set to one by default.

Of course, it sounds better. Do we want to move to solely polynomial function (i.e., delete current code and let only code for polynomial) or let both options (polynomial and "non-polynomial")? So far, I modified code to choose either from original non-polynomial behavior or polynomial one in this commit: 6f13992
UPDATE:
(current non-polynomial segmentations are used in this unit testing, so probably the second option would be better)

also note that i am starting to see a lot of similarity with dummy_centerline(). At the end, a segmentation is just a centerline with a shape around it (ellipse or rectangle). Would it make sense to merge the two functions? Maybe this is a bad idea…

Yeah, there are a lot of similarities between dummy_centerline() and dummy_segmentation(), so we could probably merge these two functions. But I am not sure, if to start with this merging within this PR, because purpose of this PR should be to create interleaved 4D data to test sct_dmri_moco -interleaved .

[jcohenadad]

Do we want to move to solely polynomial function (i.e., delete current code and let only code for polynomial) or let both options (polynomial and "non-polynomial")

simplest approach is always the best. if we can achieve the same functionality (and more) with less code, we should always go with less code

[jcohenadad]

(current non-polynomial segmentations are used in this unit testing, so probably the second option would be better)

as mentioned above ^ if we can reach the same functionality with more general code, then the unit tests should pass (if they don't pass, it means there was a change in behaviour, which we should avoid-- expect if it affects input arguments, which is manageable)

[jcohenadad]

But I am not sure, if to start with this merging within this PR, because purpose of this PR should be to create interleaved 4D data to test sct_dmri_moco -interleaved .

right-- let's do this merge in another PR

[valosekj]

as mentioned above ^ if we can reach the same functionality with more general code, then the unit tests should pass (if they don't pass, it means there was a change in behavior, which we should avoid-- expect if it affects input arguments, which is manageable)

Okay, I tried to reduce and simplify code to let only polynomial function. Unfortunately, unit tests fail (see comment below).

[valosekj]

After switch to purely polynomial curve in dummy_segmentation in last commits, test_process_seg.py fails in 2 from 26 tests.

1st fail is for segmentation with following parameters:
https://github.com/neuropoly/spinalcordtoolbox/blob/81ebad505180ab18270eb926cca4a134996f8c45/unit_testing/test_process_seg.py#L82

This fail is probably caused by slightly different segmentation when using polynomial function (red-yellow) compared to segmentation created without polynomial function (white):

2nd fail is for:

https://github.com/neuropoly/spinalcordtoolbox/blob/81ebad505180ab18270eb926cca4a134996f8c45/unit_testing/test_process_seg.py#L77

(again, polynomially created segmentation is slightly different compared to original one)

[jcohenadad]

@valosekj no worries, you can edit the "expected" values of the unit tests-- it is justified because the produced image is slightly different.

[valosekj]

Hm, I found out why segmentations were sometimes slightly different (and why unit tests failed). Problem was with conversion of fitted bspline curve from float to int. I added rounding before float -> int conversion in last commit 00af6f5 and now unit tests finished successfully.

(I use bspline instead of polynomial fit based on recommendation from dummy_centerline()):
https://github.com/neuropoly/spinalcordtoolbox/blob/67c679deadd8af0c2eb7445041b0c61628d8bbfa/spinalcordtoolbox/testing/create_test_data.py#L75-L77

[valosekj]

Now, after allowing regularized curve for dummy_segmentation(), I would continue with implementation of suggestions from this comment. Specifically, you suggested @jcohenadad:

create two synthetic segmentations (e.g. the one produced by dummy_segmentation()), by choosing random scalar for each term of the polynomial function (i.e. y=sum{a_i*x^i}, with a_i ∈ |-1, 1[)

But, bspline which I am using now does return B-spline coefficients instead of polynomial terms which are returned by polynomial fitting. How to deal with it? When I try to choose random scalar for individual B-spline coefficients, it returns to me completely different segmentation.

[jcohenadad]

@valosekj maybe you can replace bspline fit by polynomial fit, and #2754 could be addressed using low order values (e.g. max 2-3). We don't need more than that to synthesize a centerline.

[valosekj]

@valosekj maybe you can replace bspline fit by polynomial fit, and #2754 could be addressed using low order values (e.g. max 2-3). We don't need more than that to synthesize a centerline.

👍 I tried this in the last commit b1b5201 and it looks fine (all unit tests passed successfully).

Now, if I want to choose random scalar for each term of the polynomial function, it is probably necessary to do this directly in curve_fitting.py script, where I can directly achieve the poly coefficients by p.coef:

https://github.com/neuropoly/spinalcordtoolbox/blob/ba63a6212f118024b0aaa2a8a09971239e00da6b/spinalcordtoolbox/centerline/curve_fitting.py#L25-L26

by for example code like this:

p = np.polynomial.Polynomial.fit(x, y, deg)
if random_scalar:
        # choose random scalar for each term
        p.coef = [element * random.uniform(-1, 1) for element in p.coef]
return p(xref), p.deriv(1)(xref)

Because polyfit_1d function then returns fitted curve and coefficients are no longer achievable.

[jcohenadad]

Now, if I want to choose random scalar for each term of the polynomial function, it is probably necessary to do this directly in curve_fitting.py script, where I can directly achieve the poly coefficients by p.coef:

well in that case you can simply bypass curve_fitting and reimplement the call to Polynomial (it's literally 2 lines of code)

[valosekj]

I tried to implement remaining steps from this comment in the last commit b7cef13.
If interleaved=True, then for each 3D volume, two individual polynomial fits are created, each with choosing random scalar for each polynomial term. Then, these two polynomials are interleaved (odd and even slices). Finally, segmentation is created based on this interleaved polynomial fit. This is repeated for each 3D volume in 4D data.

Example output:

[jcohenadad]

Fantastic work @valosekj ! I let you rebase the branch with master and switch off Draft mode if ready for review

[jcohenadad]

great contribution!