ADC Fitting of Negative Signal

[JSousa-UoL]

I have a question regarding the ADC polyfit code. Could we make the exception in line 222 to cover all types of error. In order words, just using except: ?

https://github.com/UKRIN-MAPS/ukat/blob/c5b2df0b2416921f64993cd3ad82814633ba86a0/ukat/mapping/diffusion.py#L222

The reason I ask this is that I have an error thrown in polyfit if 'mask=None' (same as mask of Trues). This means that I was never able to run ADC map in a whole image, unless I covered/extended this except statement for all scenarios.

[alexdaniel654]

Its bad practice to have bare exceptions, linters get angry at you. LinAlgError is what the polyfit function throws if it can't manage to sensibly fit the data e.g. you're trying to fit noise. What error was being thrown when you didn't mask the data?

[JSousa-UoL]

The following error occurs when mask = None or mask = pixel_array[..., 0] > 0. It works when mask = pixel_array[..., 0] > 50 for example.

  adc_mapper = ADC(pixel_array, affine_array, bvalues_list, mask=mask, ukrin_b=False)
  File "c:\Users\md1jgra\Desktop\WEASEL\venv-windows\lib\site-packages\ukat\mapping\diffusion.py", line 172, in __init__
    self.adc, self.adc_err = self.__fit__()
  File "c:\Users\md1jgra\Desktop\WEASEL\venv-windows\lib\site-packages\ukat\mapping\diffusion.py", line 204, in __fit__
    self.__fit_signal__(sig, self.u_bvals)
  File "c:\Users\md1jgra\Desktop\WEASEL\venv-windows\lib\site-packages\ukat\mapping\diffusion.py", line 219, in __fit_signal__
    cov=True)
  File "<__array_function__ internals>", line 6, in polyfit
  File "c:\Users\md1jgra\Desktop\WEASEL\venv-windows\lib\site-packages\numpy\lib\polynomial.py", line 630, in polyfit
    raise TypeError("expected non-empty vector for x")

[alexdaniel654]

I haven't been able to replicate that issue with my sample data. It looks like its not an issue with the masking though. The error you're getting is from the actual fitting function (its saying you've not given it any bvals to fit to); the masking simply stops the fitting function from running for some voxels i.e. the fact its getting to the point of trying to fit that voxel indicates that the masking is working fine. I'd suggest its something odd with that data and that by adding a mask you're avoiding voxels it cant fit. The relevant line is this popt, pvar = np.polyfit(bvals[sig > 0], np.log(sig[sig > 0]), 1, cov=True) so I'd suggest that your data is all negative (or zero) for that voxel, therefore all the bvals get filtered out so it can't fit.

If I were you, I'd debug it slowly and see what data is causing the fit to fail. It it is that a voxel has a signal of 0 across all bvals then it wouldn't be too hard to add a line vaguely similar to this that masks out voxels that are all zeros.

Either way, the solution is to fix the code rather than make a catch all exception.

[JSousa-UoL]

Hi @alexdaniel654, I had a look at this yesterday. Thank you very much for your feedback there, that was very useful information.

You're right when you say it's data related, so I tried to do some debugging around the line popt, pvar = np.polyfit(bvals[sig > 0], np.log(sig[sig > 0]), 1, cov=True). Here is what I got printed in the console:

...
bvals[sig > 0]
[  0.   5.  10.  20.  30.  40.  50.  70. 100. 200. 300. 400. 500. 800.]
sig[sig > 0]
[6.997558   4.1985348  3.73203093 4.1985348  4.1985348  5.13154253
 4.1985348  4.43178673 3.96528287 4.1985348  3.96528287 3.498779
 3.498779   2.56577127]
bvals[sig > 0]
[  0.   5.  10.  20.  30.  40.  50.  70. 100. 200. 300. 400. 500. 800.]
sig[sig > 0]
[8.3970696  5.13154253 5.13154253 7.46406186 5.5980464  6.997558
 5.5980464  6.76430606 6.2978022  5.83129833 4.43178673 5.36479446
 4.8982906  3.498779  ]
bvals[sig > 0]
[  0.   5.  10.  20.  30.  40.  50.  70. 100. 200. 300. 400. 500. 800.]
sig[sig > 0]
[2.7990232  2.33251933 1.63276353 2.56577127 1.86601547 2.33251933
 2.0992674  2.33251933 2.0992674  1.63276353 2.0992674  1.3995116
 2.0992674  0.93300773]
bvals[sig > 0]
[  0.   5.  10.  20.  30.  40.  50.  70. 100. 200. 400. 500. 800.]
sig[sig > 0]
[1.3995116  0.6997558  1.3995116  0.6997558  0.23325193 1.16625967
 0.93300773 0.46650387 0.23325193 0.6997558  0.23325193 0.6997558
 0.23325193]
bvals[sig > 0]
[ 0. 50.]
sig[sig > 0]
[1.3995116  0.46650387]

 37%|###7      | 78405/209499 [00:45<01:16, 1720.16it/s]
Error in function DiffusionMapButton.main: the number of data points must exceed order to scale the covariance matrix

As you can see, it's not far from what you said. When there are only 2 points, the polyfit cannot obviously perform a proper fit. Therefore, I think the solution to this situation would be to add a simple if statement such as if len(bvals[sig > 0]) > 2:. What do you think? I'll create a new branch later and implement this if you agree with the solution.

[alexdaniel654]

I think a slightly different question is why are there negative/zero values in your data? Is it outside the body i.e. just noise or something else going wrong? From a consistency point of view, its probably best to enforce all the bvals are used in the fitting process. If you think you've supplied 8 bvals (and write in any papers/abstracts that you've used bvals x, y, z etc) but under the hood, the last three are omitted for every voxel that's pretty misleading.

[JSousa-UoL]

"why are there negative/zero values in your data?" - in this case simply because no mask is being used, which is a realistic and common scenario in my opinion. So yes, it happens when it's outside of the body, hence all pixels are considered in the fitting.

Is the suggestion to use bvals and sig instead of bvals[sig>0] and sig[sig>0] then?

[alexdaniel654]

I agree, running it without a mask is entirely reasonable, I was just making sure I understood where the issue was coming from. If the voxel(s) that causes the problem are outside the body then it'll just be that theres very little signal so the poor SNR can cause negative values. If the voxels that cause the problem were from anywhere that actually has signal then something else is going wrong as from a physics perspective, a negative signal shouldn't be possible.

Not quite, the >0 filter was because you cant take the log of a negative number so fitting to bvals and sig will throw an error if any voxel is sig is <0. In hindsight a more elegant solution would be better.

Leave it to me, I'll fix it in the next few days. (Also I'll hold off on the new release until I've dealt with this).

[JSousa-UoL]

Perfect, that sounds all good and reasonable. And yes, I was just about to comment on the MTR PR to hold the release as I think this should be included in v.0.4.0 as well. I'll leave it to you then