Add bugfix to NODDI_Watson model

[mkamykowski]

Hello,
We use the NODDI_Watson model and were wondering if the bugfix that was added to the NODDI_Bingham model in the v0.0.2 release could also be added to the NODDI_Watson model.
Regards,
Mike

[mcraig-ibme]

Hi Mike,

Sorry, slightly confused! I don't see any NODDI_Bingham bugfixes in v0.0.2, there was a fix for NODDI_Watson but I don't think it would apply to NODDI_Bingham. Was there something particular you had in mind?

[mkamykowski]

Hi Martin, I am now confused as well.

1. If you don't see any NODDI_Bingham bugfixes in the v0.0.2 release, then why does the description say, "Added bugfix to NODDI_Bingham model from Rob Reid"?
2. What was the the bugfix from Rob Reid? I thought he had proposed some switchover to fix the odd behavior that occurs at high ODI.

Thanks,
Mike

[captainnova]

2. What was the the bugfix from Rob Reid? I thought he had proposed some switchover to fix the odd behavior that occurs at high ODI.

It was. One of CUDA functions for evaluating a special math function has a little glitch when it switches from a power series in x to one in 1/x, which makes the optimizer stall if it hits ODI = 0.8. We think the bugfix was applied to NODDI_Bingham, but not NODDI_Watson.

The bugfix smooths over the switchover point so there is no discontinuity in the derivative, but I can't remember if it interpolates, just shifts the switchover to somewhere safer, or both.

[mcraig-ibme]

Looks like there are various reasons for confusion!

The release description I think is a typo. I cited the wrong model.

I think the relevant commit is:

8541a9f

This shows the change to the Watson model. It also shows various Bingham related added files, but this is another mistake, these were test files that were added in error under the bing2/ folder.

Apologies for all this. But it looks like the change was to the boundary value of kappa and was applied to the Watson model. I can't immediately see where a corresponding change would be made to Bingham, but if you can give me a pointer happy to apply it there as well.

[captainnova]

Hi Martin,

This is confusing but I think we are converging. Mike and I have only tested in NODDI Watson, but the problem is in a special function evaluation that might also be used in Bingham.

We agree that the bugfix was applied in Watson, and mostly fixed the bug that 8541a9f was for, but there is another bug which also kicks in around ODI = 0.75 to 0.80, so it took a while to figure out what is happening. Before 8541a9f cuDIMOT would just stop at ODI = 0.8 and not go above that. Now it does, but a histogram of ODI values shows a paucity from 0.694 to 0.758, where they pile up in a spike. I think the main error in the special function has been corrected, but there is still a little glitch in the derivative at the switchover point, and this is enough to trip up the optimizer when there is not enough pressure to get it over that glitch. I have a nice plot illustrating it and showing that it is relevant in human gray matter, but apparently I will have to either switch browsers to post it, or email it to you.

[captainnova]

Here we go...it's a plot of the OD values for a human participant calculated both with AMICO (x axis) and cuDIMOT (y), with free water fraction as the hue:

The state of the bug in cuDIMOT is the horizontal line at y = 0.758, and the partial vacuum under it. We've checked enough images to see that it corresponds to the GM and especially the GM/CSF interface (but since the cortex is thin the interface is hard to avoid), and is fairly universal across humans.

The vertical lines at lower ODI are AMICO's quantization problems. Neither axis is the ground truth, but it is a safe prior that the histogram should not have spikes (except at 0 or 1). AMICO's quantization is intrinsic to its algorithm and apparently unfixable.

I think cuDIMOT could be fixed by smoothing the transition between the "exact" and series evaluations in the special function that Watson is using. If instead of just switching over at k == 0.4, it instead linearly interpolated between the two evalulations over some zone, there would be no glitch in the first derivative for the optimizer to stumble on. For example, it could go from 100% series at k == 0.35 to 100% "exact" at k == 0.45. Unfortunately this is a lot more complicated than just moving the switchover point and seemed a bit like overkill. IIRC the actual function is quite smooth and almost linear - it is its series expansions that blow up at either 0 or 1.

[mcraig-ibme]

Ok, I have implemented something along these lines on the branch issue/4. If you are able to build this branch, could you run your test case against it and let me know if it is having the effect you expect. If you have trouble compiling I can provide binaries built on Redhat 9 / cuda 12.1 to test.

Also, feel free to review the code change, there is a pull request open for it

[captainnova]

Sorry for not replying earlier; I had reviewed the code right away and it looked good, but I got pulled onto other things before I could test it with an actual scan. It solves the problem of this issue, and I recommend accepting the PR.

This plot is a joint histogram of OD for that scan, with post-change on the vertical axis and just-before-issue/4 on the horizontal, with marginal distributions in the margins. Most of the values are essentially identical, but the old spike at OD = 0.76 is now gone. The new OD image has a broader bump between OD = 0.75 and 0.80, but the difference image shows the affected voxels are the high end of the natural GM distribution.

There are still sharper spikes at 0 and 1, but I think those are from including 0 and 1 as initial values in the grid search - the optimization isn't able to get off them. I saw where the initialization schedule is and plan on tweaking it soon, but it is a separate issue.

Thanks for making the change!