Densities, random walks & travelling salesman #206

Daval-G · 2024-11-01T18:43:16Z

This PR adds densities, sampling, and trajectories based on that (random walks and travelling salesman). In particular, it reproduces 3 papers from Nicolas Chauffert (4th coming "soon"). The features are represented in the images below, but in details:

maths.py is split into a submodule as the file became too big
functions were added to generate density arrays in 2D and 3D: cutoff-decay, polynomial, energy-based, Chauffert & fast Chauffert. The last one corresponds to an adaptation of the matlab code, but it uses only 1D vectors to compute the Fourier basis which appear to yield incomplete densities.
inits submodule to store the trajectory initializations, because trajectory2D.py and trajectory3D.py are also starting to be too big
Random walk trajectories in 2D and 3D based on any density
Travelling Salesman trajectories in 2D and 3D with options to cluster shots before computing the TSP solution for time efficiency.

TODO:

Add docstrings
Write an example
Add references
Add pywavelets dependency

Density features:

Trajectory features:

Daval-G · 2024-11-01T18:55:32Z

@paquiteau I need to add a dependency to pywavelets, is it OK ? Should I prefer something else ? And should I add it as an optional dependency or not ?

paquiteau · 2024-11-01T19:01:24Z

Whoa there is a lot to process, but this look awesome. @philouc or @chaithyagr you might want to take a look ;)

For pywavelets, you can add it as an optional depedency (and give proper error message when not found) .

philouc · 2024-11-02T04:28:58Z

@Daval-G did you make sure to define the optimal target sampling density (TSD) for TSP solutions as a function of the dimensions of the ambient space? In 2D (d=2), if \pi is supposed to be the optimal TSD, the TSD for TSP should be parameterized as \pi^2 as the TSP itself works as \pi^{d-1/d} [see section 4.3, p. 1976 in Chauffert et al, SIAM IS 2014] whereas in 3D the target TSD should be initialized as \pi^{3/2}.

Daval-G · 2024-11-02T22:01:00Z

@philouc Good catch, I overlooked that point.

I followed the experiment proposed in the paper and presented in Fig 4.1, a.k.a I made 10 000 trajectories/sampling (Nc=10, Ns=100) in different situations and then made a 2D histogram to find if the target density was actually respected or not. I based that experiment on a dummy TSD as according to the paper it is not specific to Nicolas's optimal TSD.

In the figure below, (a) is the TSD and (b-f) are normalized histogram corresponding to the empirical distributions:

Note that:

"Pseudo-random walk" (c) is just a random walk where the probability is adjusted dynamically to reduce the probability to go back to already covered areas. It is an option provided to the users for the random walk.
"Random TSP x4" (f) means that the trajectory consisted of a Travelling Salesman (TS) based on random sampling where only 1/4 of points were drawn, and 3/4 of points are simply interpolated along the TS path/solution. I did that because otherwise it would be strictly equivalent to just random sampling.
"Lloyd sampling" was actually sampled only 1000 times because it takes a lot more time than the other ones.

Overall it seems that the TSD is already respected without additional correction (except for (e) but it was not expected anyway). I am not sure why it is the case here and not in the paper, I am still looking into it.

paquiteau

A first pass, the code is looking good to me, I let @chaithyagr or @philouc comment on the maths.

I guess you already plan to add more documentation and examples?

Also, some of the internal maths functions can be quite complex, so if you have some debugging setup that you wish to promote to unit-tests, feel free to do so as well.

src/mrinufft/trajectories/densities.py

src/mrinufft/trajectories/inits/random_walk.py

src/mrinufft/trajectories/maths/constants.py

src/mrinufft/trajectories/maths/tsp_solver.py

Daval-G · 2024-11-03T12:02:31Z

@paquiteau Yes I made a todo list in the original PR post about docstrings, examples & co. I converted it to draft

philouc · 2024-11-04T06:17:43Z

You're right @Daval-G. Illustration on a dummy TSD is fine. I was just referring to the math relation allowing us to compute the sampling density on the TSP (points are joined by straight lines) vs the original TSD which could be sampled independently. Based on this comment I'm not sure whether TSPx4 shown in panel (f) fulfills the equation because if you start from the TSD in (a) you should end up with the sqrt of (a) in (f). Fig. 6 in Nicolas' paper illustrates pretty well that phenomemon when contrasting panes (d) and (g).

chaithyagr

Starting my review...

src/mrinufft/trajectories/densities.py

chaithyagr · 2024-11-04T14:15:53Z

src/mrinufft/trajectories/densities.py

+    return locations
+
+
+def create_cutoff_decay_density(shape, cutoff, decay, resolution=None):


How about a more generic density defined in Eq. 10 in https://onlinelibrary.wiley.com/doi/full/10.1002/mrm.29702 ?
Its a very short update.
Also, having support for elliptical sampling enabled or not can be helpful (we soon plan to disable it :P )

I am confused, this density is already based on that equation. What difference do you see ?
And what do you mean by elliptical sampling ? Anisotropic ? Why do you plan to disable it and from where?

Anisotropic and more importantly elliptical sampling, a way to enable or disable mask.

Somehow that dodges all of my questions xD Could you reformulate pls ?

We need to support elliptical or non elliptical density

We also need to support anisotropy, but I think that we already do?

I had to fix it but I think it can now handle anisotropy and elliptical densities

chaithyagr · 2024-11-04T14:28:14Z

src/mrinufft/trajectories/densities.py

+        fourier_vector = nf.ifftn(unit_vector)
+        coeffs = pw.wavedecn(
+            fourier_vector, wavelet=wavelet_basis, level=nb_wavelet_scales
+        )


You seem to be doing \Psi F*, the paper describes F* \psi right?

I based this on the original code from Nicolas:

https://github.com/philouc/mri_acq_recon_tutorial/blob/master/matlab/TSP_Markov/tools/compute_optimal_distrib.m

https://github.com/philouc/mri_acq_recon_tutorial/blob/master/matlab/TSP_Markov/tools/perform_wavortho_transfAnis.m

I checked that the code had an identical behavior, and it is the case except for one major difference: he used hardcoded wavelets banks and I use pywavelets instead, and their values seem to differ.

chaithyagr · 2024-11-04T14:32:29Z

src/mrinufft/trajectories/densities.py

+    for k, s in enumerate(shape):
+        unit_vector = np.zeros(s)
+        density_1d = np.zeros(s)
+
+        for i in range(s):
+            unit_vector[i] = 1
+            fourier_vector = nf.ifft(unit_vector)
+            coeffs = pw.wavedec(
+                fourier_vector, wavelet=wavelet_basis, level=nb_wavelet_scales
+            )
+            coeffs, _ = pw.coeffs_to_array(coeffs)
+            density_1d[i] = np.max(np.abs(coeffs)) ** 2
+            unit_vector[i] = 0


Is it possible to vectorize this? How slow is it right now for 3D?

@paquiteau , maybe numba?

create_fast_chauffert_density here is super fast for any matrix size. It is not the same for create_chauffert_density above

pywavelets allows the vectorization, but for create_chauffert_density we are speaking about an NxN matrix with N the total number of voxels. We explicit the linear operator to find for each row the max L2 norm. There are probably smart ways to save some computations though.

chaithyagr · 2024-11-04T14:33:01Z

src/mrinufft/trajectories/densities.py

+    return nf.ifftshift(density)
+
+
+def create_fast_chauffert_density(shape, wavelet_basis, nb_wavelet_scales):


To keep in mind: https://arxiv.org/pdf/1910.14513, some other interesting results with Anna.

Is there something specific in the paper you would like to highlight ? Should we open an issue to create some extension in the future ?

Also @chaithyagr I note in the paper figure that there is a TWIRL implementation sleeping somewhere in our cardboards. Do you know where it is so it can be added to the package ?

Isnt TWIRL just TPI in 2D? It should be fairly straightforward from our codes.

Is there something specific in the paper you would like to highlight ? Should we open an issue to create some extension in the future ?

I think @philouc can comment better than me. I just felt its good to add in docs, Its been long having gone through this work, but from what I remember its extension of the earlier work for anisotropic sampling, which can be of interest too.. But surely not to be pursued in this work

chaithyagr · 2024-11-04T14:52:31Z

Based on this comment I'm not sure whether TSPx4 shown in panel (f) fulfills the equation because if you start from the TSD in (a) you should end up with the sqrt of (a) in (f). Fig. 6 in Nicolas' paper illustrates pretty well that phenomemon when contrasting panes (d) and (g).

I think the key missing ingridient here is this point by Guillaume:

x4" (f) means that the trajectory consisted of a Travelling Salesman (TS) based on random sampling where only 1/4 of points were drawn, and 3/4 of points are simply interpolated along the TS path/solution.

I am assuming this means he is interpolating UNIFORMLY with 3 extra points in between every 2 points in trajectory. (@Daval-G to correct me).
Original work assumes that ALL the samples in between 2 points in TS solution are fully sampled (which results in the whole p^{-1/d} to begin with) rather than a fixed percentage. If we do that I am sure @Daval-G could have way more points which would potentially mess up his density. It is a different question about how to go about this. For now, it seems like as long as the trajectories are not scanner compliant, it really would not matter. But, if you plan to project this on constraint set, the density will drastically change.

Daval-G · 2024-11-05T19:52:53Z

I am assuming this means he is interpolating UNIFORMLY with 3 extra points in between every 2 points in trajectory. (@Daval-G to correct me).

Actually I interpolate with cubic spline, it is maybe the best way to account for 1st and 2nd order derivative to get smooth transitions. It might be the reason.

Original work assumes that ALL the samples in between 2 points in TS solution are fully sampled (which results in the whole p^{-1/d} to begin with) rather than a fixed percentage. If we do that I am sure @Daval-G could have way more points which would potentially mess up his density.

That's a good point, so I carried an additional test where I just increase a lot the interpolation/oversampling while keeping the same number of TSP points (Ntsp). In the image below, we start with (b) at Nc=20 and Ns=200 but Ntsp=50, which is already 10% sampling/AF=10 compared to the 200x200 matrix size used for the density. The oversampling from Ntsp to Ns goes from 4 to 64 which is HUGE. I provide example trajectories below (clustered by bands to accelerate computations). At this point I think we can say that we sampled every point on the path followed by the TSP, and the histogram/observed density is still not affected.

For now, it seems like as long as the trajectories are not scanner compliant, it really would not matter. But, if you plan to project this on constraint set, the density will drastically change.

All the trajectories are defined over a normalized space [-0.5, 0.5]^d so there is necessarily a resolution and/or readout duration for which they are compliant. Also oversampling with cubic splines make transitions quite smooth already, given enough oversampling.

@chaithyagr At this point I think it is relevant anyway to add an option to account for the d/(d-1) density exponent to fit the paper proposition.

Daval-G · 2024-11-05T20:05:33Z

@chaithyagr I think the interpolation is indeed the key, I will carry extra tests this week to check that

Edit: Actually no, the results are similar to the ones from the figure above with/without clustering and linear or quadratic interpolation

chaithyagr · 2024-11-06T09:08:35Z

Im sorry for the confusion, when I said UNIFORMLY in my above statement, I mean you added EXACTLY the same number of points in between any 2 points in the trajectory (which also makes sense).
But the original paper from what I see, seems to have more points if the points are far apart as compared to be near. I think we dont need to care about the extra factor for realistic usecases, but as you tell it is good to have an option, especially assuming that in the end we surely can get max gradient constrained making the overall sampling pattern different (and more close to whats shown in paper), but only upon projection.

Daval-G · 2024-11-06T17:31:59Z

@chaithyagr That was it, great job. The paper considers binary masks instead of continuous spaces, so the number of points within a k-space voxel doesn't matter as long as there is one or more, effectively impacting the density in the way described by Nicolas.

Below is an experiment showing that, with empirical densities obtained from continuous spaces put into 2D histograms and just binary masks as presented. The density with binary masks is attenuated, and actually fixed when sampling over pi^(d/(d-1)).

That being said, using binary masks for k-space sampling is quite specific. It means one wants to undersample retrospectively without using NUFFT operators. It is quite the opposite of the goal of this package, and can't be used for real life acquisitions. I would advocate to provide the d/(d-1) density correction as an option with default to False, and specify its use case in the documentation.

chaithyagr · 2024-11-06T17:33:57Z

Great job on the experienent! Yes I agree with you on the default as False.

philouc · 2024-11-14T21:27:21Z

Very nice job @Daval-G ; Your illustration in 2D is clear and very convincing. Additionally, suggestion about deactivating the density correction seems reasonable as the main purpose of this MRI-NUFFT package it to implement undersampling over a continuous k-space domain.

philouc · 2024-11-14T21:27:43Z

Very sorry for the late feedback, I'm travelling every week!

Daval-G · 2024-12-02T22:58:36Z

I will open an issue later to add sampling.sample_from_density as an option for Wave-CAIPI rather than now

Daval-G · 2024-12-09T17:40:27Z

@paquiteau I need to add dependencies to pywavelets (already mentioned) and scikit-learn. I see that we have a lot of optional dependency groups, which make sense for big GPU libraries that are not trivial to install, but among them I see really common libraries like scikit-image in [smaps]. Can't we just have those as base dependencies or do we really want nothing more than numpy, scipy, matplotlib & tqdm ?

paquiteau · 2024-12-10T08:28:57Z

Indeed the dependencies management in MRI-Nufft could be improved.
Aside the nufft backend, we could have an extra group where we have the dependencies you mentioned for trajectory as well as the io and smaps group.

Then I think we could recommend people to pip install mri-nufft[extras] in the readme and docs

chaithyagr · 2024-12-11T15:48:17Z

Maybe scikit-learn needs to be added for dependency in CI? I would recommend waiting for #195 to go in, which makes the process steramlined and only adding in pyproject.toml is enough

paquiteau · 2024-12-11T20:09:27Z

the CI should simply install mri-nufft[extras] to reflect what user are doing :)

Daval-G · 2024-12-16T11:56:20Z

@chaithyagr @paquiteau CI updated & passed. Let me know if good to merge, I am planning the next PRs already

paquiteau · 2024-12-16T12:06:17Z

LGTM ! Thanks for this — once again — significant contribution :)

Daval-G · 2024-12-16T17:18:18Z

@paquiteau The merge is blocked until Chaithya & you click on approve :)

chaithyagr · 2024-12-16T17:30:07Z

Given a lot of docs and also an example, can we run the docs in CI? You can get a run by adding [docs] tag in your commit.

chaithyagr · 2024-12-16T17:30:17Z

Otherwise this PR LTGM..

Daval-G · 2024-12-16T18:36:14Z

Given a lot of docs and also an example, can we run the docs in CI? You can get a run by adding [docs] tag in your commit.

@chaithyagr Done

Daval-G added 2 commits November 1, 2024 12:22

Split trajectory.maths.py into a submodule

a26e6ca

Add densities, random walk and travelling salesman trajectories

0f12b7e

Daval-G added the trajectories Issues concerning Non cartesian trajectories label Nov 1, 2024

Daval-G requested a review from paquiteau November 1, 2024 18:43

Daval-G self-assigned this Nov 1, 2024

paquiteau requested a review from chaithyagr November 3, 2024 10:37

paquiteau requested changes Nov 3, 2024

View reviewed changes

Daval-G marked this pull request as draft November 3, 2024 12:02

Daval-G added 2 commits November 3, 2024 23:19

fixup! Add densities, random walk and travelling salesman trajectories

e72f52e

Add docstrings for random walk and TSP solver

1665906

chaithyagr requested changes Nov 4, 2024

View reviewed changes

Daval-G added 2 commits November 7, 2024 22:30

fixup! Add densities, random walk and travelling salesman trajectories

55b3ba4

Add docstrings to TSP & sampling densities

74b08f8

Daval-G added 3 commits November 23, 2024 23:46

Add samplings to random walk initializations and fix boundaries

362ca9f

Add utils to examples to show densities

fcf699b

Add sampling density example skeleton

fdcda55

Daval-G added 5 commits November 24, 2024 22:35

Merge branch 'master' into chauffert_update

3d2d043

Add descriptions to sampling density examples

2f631ca

Add dependency to pywavelet

83b9081

Merge branch 'master' into chauffert_update

3fc0bba

Improve cutoff/decay density, fix dependencies

82172be

Daval-G mentioned this pull request Dec 2, 2024

Random sampling #220

Open

1 task

Daval-G requested review from chaithyagr and paquiteau December 2, 2024 22:57

Daval-G added 4 commits December 10, 2024 23:42

Gather extra dependencies, homogeneize extra imports

34144b8

Fix imports to pywavelets

e7589f2

Fix black format in smaps

432df36

Remove except duplicate in io/siemens

4ff1bf0

Daval-G marked this pull request as ready for review December 10, 2024 23:00

Merge branch 'master' into chauffert_update

cbb1f83

Update CI with new [extra] dependencies

e60da32

Daval-G mentioned this pull request Dec 12, 2024

Add sampling methods to Wave-CAIPI #227

Open

paquiteau approved these changes Dec 16, 2024

View reviewed changes

[docs] Trigger documentation build

2f34957

paquiteau merged commit 9c41feb into master Dec 16, 2024
12 checks passed

		return locations


		def create_cutoff_decay_density(shape, cutoff, decay, resolution=None):

		return nf.ifftshift(density)


		def create_fast_chauffert_density(shape, wavelet_basis, nb_wavelet_scales):

Densities, random walks & travelling salesman #206

Densities, random walks & travelling salesman #206

Conversation

Daval-G commented Nov 1, 2024 • edited Loading

Daval-G commented Nov 1, 2024

paquiteau commented Nov 1, 2024

philouc commented Nov 2, 2024

Daval-G commented Nov 2, 2024

paquiteau left a comment

Choose a reason for hiding this comment

Daval-G commented Nov 3, 2024

philouc commented Nov 4, 2024

chaithyagr left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Daval-G Nov 5, 2024 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

chaithyagr commented Nov 4, 2024

Daval-G commented Nov 5, 2024 • edited Loading

Daval-G commented Nov 5, 2024 • edited Loading

chaithyagr commented Nov 6, 2024

Daval-G commented Nov 6, 2024

chaithyagr commented Nov 6, 2024

philouc commented Nov 14, 2024

philouc commented Nov 14, 2024

Daval-G commented Dec 2, 2024 • edited Loading

Daval-G commented Dec 9, 2024

paquiteau commented Dec 10, 2024

chaithyagr commented Dec 11, 2024

paquiteau commented Dec 11, 2024

Daval-G commented Dec 16, 2024

paquiteau commented Dec 16, 2024

Daval-G commented Dec 16, 2024

chaithyagr commented Dec 16, 2024

chaithyagr commented Dec 16, 2024

Daval-G commented Dec 16, 2024 • edited Loading

Daval-G commented Nov 1, 2024 •

edited

Loading

Daval-G Nov 5, 2024 •

edited

Loading

Daval-G commented Nov 5, 2024 •

edited

Loading

Daval-G commented Nov 5, 2024 •

edited

Loading

Daval-G commented Dec 2, 2024 •

edited

Loading

Daval-G commented Dec 16, 2024 •

edited

Loading