Problems with gen_scheme generated multishell gradient table

Hello MRtrix team,

I am trying to extract a 30/20/10dirs subset from my 100/65/32 dirs multishell datasets (3 bvals) which is shown below:

(I know the first shell is pretty bad, but that’s the dataset they gave me! :wink: ).

I want to ensure uniform sampling in each shell but also on the overall sphere. Following recommendations on other threads, I used your gen_scheme command to generate an ideal multishell 30/20/10 scheme and then find the closest match in my dataset. My understanding was that gen_scheme tries to maximize uniformity across and within shells, just like I needed. However, the gen_scheme output is not as good as I hoped (picture below). As you can see, directions are nicely uniform in each shell, but the overall distribution (top) is not so good.

Isn’t gen_scheme supposed to estimate overall bipolar energy as well? Do you have any suggestion on how I could extract the datasets that I need in another way? I also thought to use dirgen to generate a 60dirs scheme and then use dirsplit on it, but I could find a way to set the number of directions I want in each dirsplit set :frowning:

Thank you in advance for your help,

Giorgia Grisot

That’s not quite right: it doesn’t try to maximise uniformity across shells, only within shells. The reason for this is that I don’t see the need for uniformity across shells: if each shell is adequately sampled, there is no further information to be gained by modifying the sample point positions to avoid the same sampling points as in the other shells - quite the opposite in fact, it would lead to a less uniform distribution than would otherwise be obtained.

The way I rationalise this is to consider the case of 2D Cartesian space: if the signal is properly sampled (above the Nyquist rate), shifting the sampling points in any one row does not add any information, it will still be perfectly sampled. Here, if x corresponds to the angular domain and y to the b-value domain, you can see that the same considerations apply. Assuming the signal is sufficiently sampled in the angular domain, there is nothing to be gained by moving the sample points around - but plenty to lose.

The only caveat to that is that if the signal is undersampled, and the reconstruction imposes a strong model that links between shells, then there may be benefits in such a strategy - but my personal opinion is that if you are operating in that regime, you will probably be lacking in SNR anyway (at least for most higher-order models). But that’s just my 2 cents…

I’m not sure what you mean here… I assume this relates to the previous point, about providing uniformity across shells…? It does estimate the bipolar energy at various points, but only does so per shell.

I think your best bet may be dirorder: it orders the directions in the set in a way that ensures near-uniformity upon truncation. It’s not guaranteed to be the optimal solution though, given the way it operates - you should be able to do better than that since you know how many directions you want (dirorder tries to find an order that’s OK for any arbitrary truncation, and even then doesn’t try that hard to find the optimal order).

Yes, dirsplit is not designed for this, although it’s otherwise doing pretty much what you’d need. It’s a brute-force search through lots of combinations to find a split of the data that provides near-uniformity for each subset - but it does expect to split the directions into equal sized subsets. It wouldn’t take a huge amount of effort to make it do something like what you need, but that would require a bit of coding…

However, it would still not guarantee that the directions are also uniformly spread over the whole sphere (I assume you need this to avoid biases due to eddy-currents, as recommended by the FSL eddy page?) - that’s fundamentally incompatible with ensuring uniformity under a bipolar model. What gen_scheme does is to first obtain the most uniform direction set per shell under a bipolar model (using dirgen), then modify these directions by inverting some of them (which does not change the energy in a bipolar model) so as to minimise the energy assuming a unipolar model (using dirflip). I’m not sure how I’d manage to balance the two cost functions involved in selecting a subset of your existing directions that is uniform under both bipolar and unipolar models – pretty sure I’d overwhelmingly prioritise the bipolar model though, and totally ignore uniformity under a unipolar model: if it turns out to be vaguely spread out over the sphere, that’s great, otherwise I’d rather have the most uniform set of directions I can get from a bipolar point of view, since that is what determines the stability of my reconstruction.

Not sure that helps you hugely, but hopefully there’s some useful information in there… (?)

Hi @Giorgia_Grisot,

I noticed you used/modified the little Matlab snippet I posted around here somewhere to visualise those direction sets. Great to see! :heart_eyes:

I agree with @jdtournier on most points. In general, as long as you have a decent angular resolution on the highest b-valued shell, I agree with this:

…for most scenarios indeed.

That’s really important to stress as well. Always give preference towards an optimal spread of gradient orientations, regardless of direction/polarity. After optimising that, dirflip can at least make things slightly better with respect to eddy currents; but it should be strictly only a minor/second concern. Note that dirflip doesn’t change the orientations anymore, it just flips polarities; which is often enough to at least “spread” the polarities around.

Hello @jdtournier and @ThijsDhollander,

thanks for getting back to me.

The reason for this is that I don’t see the need for uniformity across shells: if each shell is adequately sampled, there is no further information to be gained by modifying the sample point positions to avoid the same sampling points as in the other shells - quite the opposite in fact, it would lead to a less uniform distribution than would otherwise be obtained.

I see your point but I am still a little confused. I guess the question is “how do you define angular resolution”? When using a multishell sequence, I would project all the directions of each shell on the same sphere and look at the angle between them. Now, and maybe I am wrong, If I’d like to design a sequence with 60 directions that maximizes my angular resolution, wouldn’t I want to ensure that the directions in each shell don’t overlap with each other in order to avoid sampling basically the same direction twice, just at two different bvalue (which sometimes is what you want, but It’s not the case here)? Wouldn’t it be better from an angular resolution stand-point to indeed “shift” or rotate those uniform points around on each shell as to also maximize the overall angular resolution?
Below is for example the MGH-USC HCP gradient table, which shows my point:

Or am I thinking it all wrong?

PS: @ThijsDhollander, yes that’s your code, just minimally adjusted to show different bvals :slight_smile: I couldn’t resist the cool shading :sunglasses:

1 Like

Hi @Giorgia_Grisot,

For your problem, probably you could try the method of mixed integer linear programming (MILP) that we recently proposed.
The method tries to maximize the angular separation between samples within shells and across shells.
Please see a demo to uniformly extract 30/30/30 samples from 90/90/90 samples in the HCP scheme.

I agree with @jdtournier that for a dense sampling (or a dense sub-samping), there is no further information to be gained by further optimizing the sampling points. However, for a undersampled scheme, e.,g. in your case (30/20/10), optimizing the uniformilty of sub-sampled points may have a better reconstruction result, especially when you use a multi-shell reconstruction method.

Jian Cheng

Sounds like @JianCheng’s solution might be more appropriate to your problem than what we have in MRtrix3 currently…

But to clarify my earlier statements:

To some extent, but I don’t think that’s too ambiguous - at least at the single shell level. But what I think matters more is the intrinsic angular resolution of the signal, and therefore what angular resolution is needed at the acquisition to fully capture all the angular features of the signal - an issue I’d tried to address in this paper. The point is that the diffusion signal is relatively smooth as a function of orientation, and as long as the sampling is sufficiently dense relative to the smoothness, any further increase in sampling density does not increase the angular resolution of the reconstruction (although it does increase SNR as expected) - it is already fully sampled.

So from this point of view, trying to avoid the same sample point across shells can’t actually increase the angular resolution of the final reconstruction - assuming it was already fully sampled on each shell. But it can worsen the stability of the reconstruction if enforcing this leads to deviations from what would otherwise have been the optimally uniform distribution of orientations.

Hopefully the above will clarify. But basically, it shouldn’t matter at all (again, assuming the signal is fully sampled on each shell). I find my Cartesian space analogy quite helpful here. If there was anything to be gained from ensuring sample points don’t overlap across shells, the Cartesian equivalent would be to ensure sample points in one row fall between the sample points in adjacent rows. Yet no-one does this - even though it would actually be trivial to do in MRI. In fact, it can be shown that for a fully sampled signal (i.e. whose frequency content is below the Nyquist rate), the information content is strictly equivalent - there is no benefit to shifting the sample points along. The same logic can be extrapolated to the spherical case, given the strong similarities between Cartesian Fourier theory and its spherical equivalent.

The only case where I might agree with this is when the signal is under-sampled on each shell, and the reconstruction is sufficiently constrained within and across shells to avoid aliasing artefacts (i.e. undersampled high frequency components wrongly manifesting as low frequency components). A tensor model is an example of such a reconstruction (although it is generally a poor fit in practice due to the non-mono-exponential decay of the signal, and crossing fibre effects, and is not immune to instabilities in this regime) - but more complex models are unlikely in my opinion to impose sufficiently strong constraints to avoid such issues with undersampled schemes. Multi-tissue CSD would also impose strong constraints by virtue of the per-tissue responses, but I’m not sure how strong such an effect would be on its output - I’d guess it would be affected at least to some extent, above and beyond the low SNR that such an acquisition would inevitably provide…

I hope this at least clarifies my current position on this - as always, I might be wrong (wouldn’t be the first time), but that’s my 2 cents for what it’s worth…

I also still agree with @jdtournier’s points on this one, but it’s indeed worthwhile to consider this bit of nuance:

I remembered @jdtournier and myself at some point having this argument, especially with respect to:

Took me a while, but I’ve found the relevant (GitHub) posts again; see in particular these 2 posts:

Note that, what you can still do as a sort of compromise is still sample every individual shell as uniformly as possible (independent of each other shell), but subsequently rotate those individual shell sample schemes “as a whole” with respect to each other; i.e. not affecting the uniform sampling per shell in any way.

However, with CSD techniques in particular in mind, I’d advise to also consider your numbers of directions per shell carefully; especially aiming at sampling the highest b-value shell densely enough so as to essentially cover that intrinsic angular resolution of the signal for that one. In my experience, that is vastly more important to the quality of the CSD outcome then about anything else; the reason being that contrast-to-noise ratio of the signal is very important for a CSD technique. If you’re strictly after 60 directions, I’d simply lay them all out on a b=3000 shell uniformly; end of story. But if for any reason, you see value in sampling another b-value, why not go for e.g. 15 directions at b=1000 and 45 directions at b=3000, both uniformly distributed on their own shell.