dwi2tensor: add wls

[jdtournier]

This follows recent bits of experimentation and discussions with @bjeurissen. Just making sure it's up for discussion and inclusion in RC4...

[EDIT]
These changes follow an observation by @dchristiaens of strange behaviour on some of our multi-shell data: we ended up with negative ADCs in some voxels, particularly where it would be expected to be fast. Digging into it further, we figured what happened was that the initial OLS fit was dominated by the noisy DW signals, which at that stage aren't down-weighted in any way. When you have relatively few b=0 images, it essentially ends up ignoring them and trying to fit primarily to the noise - and due to the log transform, the noise can rapidly swamp the fitting. At this point, because the next iteration uses the predicted signal to compute the weights, the b=0 signals end up down-weighted even further relative to the DW signals, further compounding the problem.

To minimise this issue, I thought we could initialise with a weighted least-square fit, with the weights driven by the measured signal rather than its prediction - which does indeed sort out the problem we were having. So in this PR, that's now the default behaviour. I can't remember whether it made any meaningful difference in those parts where the current approach actually converged - I'd need to double-check that, but from memory whatever differences there might have been were negligible.

I then looked a bit into the literature, and was surprised to (re)discover that the WLS approach is what had originally been proposed for the fit by Peter Basser back in his original 1994 paper... So I thought we should at least offer it as an option in dwi2tensor.

So the proposed changes try to rationalise all this by offering a choice of initial fitting (WLS by default, or OLS), followed by a user-specified number of iteratively re-weighted LS (IWLS) with weights determined by the signal predictions. Means you can pick either OLS or WLS for the direct fit (set -iter to zero) or as initialisation for however many rounds of IWLS you want.

[thijsdhollander]

I'm ok with extra algorithms, but I'd like the default behaviour to remain unchanged. Can you briefly elaborate on what's essentially changed (don't have much time now to go through all the code changes)?

[thijsdhollander]

extra algorithms

...or of course extra options, etc.. essentially allowing for new functionality, whether that's labelled "new algo" or not. All good on that front; I'm more interested in what is or remains the default behaviour.

[jdtournier]

Yep, I'll freely admit that was pretty lite on justification... I'll edit the main post with the details.

[jdtournier]

OK, justification for changes now documented in main post.

I'm more interested in what is or remains the default behaviour

Hopefully these edits will answer that question. But in essence: yes, there is a change in the default behaviour, this would initialise the IWLS with a WLS fit rather than the current OLS fit.

[jdtournier]

OK, coming back to this discussion about changing dwi2tensor to include WLS as the default tensor fit, followed by 2 rounds of IWLS, with the option of using OLS for the initialisation: here's some screenshots of what the FA looks like with the various options.

Input data was fairly low-grade b=1300, 32 direction data.

current master:

default: initial OLS + 2x IWLS

iter 0: initial OLS only

Proposed:

default: initial WLS + 2x IWLS

iter 0: initial WLS only

-ols -iter 0: initial OLS only (same as current master with iter 0)

---

A straight subtraction between the two OLS results gives exact zeros - so they're completely equivalent. Here's the difference between the IWLS results (only difference is initialisation), windowed from -0.01 to 0.01:

Differences primarily in CSF, and sporadic voxels mostly at the edge of the brain. Differences are otherwise pretty hard to spot - but there is undeniably a difference. In general, the new results are better behaved (fewer outlier FA values).

The results are also better on our multi-shell data, but I've not gotten round to generating figures for it...

Thought welcome.

[Lestropie]

Can I propose some changes to the documentation? There's a bit of nuance in the description of exactly how the fit is done, so I think this could do with a reasonably verbose explanation in the Description section; and possibly also redirecting some of the explanation currently tied to specific command-line options into the Description section.

Currently we have:

 DESCRIPTION
    + "Estimate diffusion tensor (and optionally its kurtosis) by fitting to "
    "the log-signal using either ordinary least-squares (OLS) or weighted least-squares (WLS), "
    "followed by further re-weighted iterations using the previous signals predictions (IWLS)."
    "By default, the initial fit is performed using WLS, followed by 2 IWLS iterations."

 OPTIONS
+ Option ("ols", "perform initial fit using an ordinary least-squares (OLS) fit, "
"instead of the default weighted least-squares (WLS) strategy.")

+ Option ("iter","number of iterative reweightings for IWLS algorithm (default: "
        + str(DEFAULT_NITER) + "). Set to zero to perform standard OLS or WLS.")
+ Argument ("integer").type_integer (0, 10)

Maybe something like this would be clearer:

 DESCRIPTION
+ "By default, the diffusion tensor (and optionally its kurtosis) is fitted to "
  "the log-signal in two steps: firstly, using weighted least-squares (WLS) with "
  "weights based on the empirical signal intensities; secondly, by further iterated weighted "
  "least-squares (IWLS) with weights determined by the signal predictions from the "
  "previous iteration (by default, 2 iterations will be performed). This behaviour can "
  "be altered in two ways: "

+ "* The -ols option will cause the first fitting step to be performed using ordinary "
  "least-squares (OLS); that is, all measurements contribute equally to the fit, instead of "
  "the default behaviour of weighting based on the empirical signal intensities."

+ "* The -iter option controls the number of iterations of the IWLS prodedure. If this is "
  "set to zero, then the output model parameters will be those resulting from the first "
  "fitting step only: either WLS by default, or OLS if the -ols option is used in conjunction "
  "with -iter 0.";

 OPTIONS
+ Option ("ols", "perform initial fit using an ordinary least-squares (OLS) fit (see Description).")

+ Option ("iter","number of iterative reweightings for IWLS algorithm (default: "
        + str(DEFAULT_NITER) + ") (see Description).")
+ Argument ("integer").type_integer (0, 10)

[thijsdhollander]

Here's the difference between the IWLS results (only difference is initialisation)

This is all I really mostly care about; i.e., the old and new default. All the other option combinations made possible are entirely fine: it's up to the users to use them responsibly, but they do exactly what they say on the box.

In general, the new results are better behaved (fewer outlier FA values).

So this is what I'd like to be sure of we get across the board. In the example you provides it's surely fine. Not much changes, but what does change is mostly for the better. It does show that it can also act the other way around for a few individual voxels though.

The results are also better on our multi-shell data, but I've not gotten round to generating figures for it...

Would be happy to see an example of that; this is referring to the (kind of) data that initially motivated the change in default, correct? As mentioned above, I'd concentrate on the difference between old and new default. So essentially (in case of your data above), these figures (old / new / difference):

Thinking a bit about this myself, I'm getting more and more convinced this should be for the better though. Just being careful (and thus naturally critical).

[jdtournier]

Happy with the change in description - I admittedly didn't invest a huge amount of effort into that side of things...

Otherwise, I'll get a few more examples of the difference on multi-shell data when I get the chance next week.

[jdtournier]

OK, results for the BATMAN data (following all the pre-processing in the tutorial) - as before, FA for old / new / diff (windowed [-0.01 0.01]):

As before, most differences in the CSF regions, but a touch of a difference in some regions of the WM too. Biggest differences observed in WM is of the order of 1.5% (with new > old).

---

Neonatal data from the dHCP cohort, same as before:

Note the few glitches in the old OLS-initialised results near the brainstem / temporal pole, which you don't see in the WLS-initialised results (this is the kind of issue that prompted this in the first place). Difference between new & old is again of the order of 1-2% in the WM, typically with new > old (at least on FA).

---

If everyone is happy with that, I suggest we merge...?

[thijsdhollander]

Can't fault it given the data and as far as I got with my thinking about it. 👍 I'll keep an eye on it though, just in case I see something funny pop up in the future. But having the -ols option will allow for a quick check then.

[thijsdhollander]

Did some extra testing from dev in the mean time with a range of data: all still fine. Nothing pops up that would concern me. I was also able to confirm the above with some dHCP data.

[jdtournier]

👍