[BUG] bval=5 is interpreted as a non zero gradient

[tashrifbillah]

The issue has been discovered with HCP data that have bvalues:

5       5       1490     1495       1500      1495     1490       ...

As a result, the program throws an exception and fails:

No zero gradients! terminate called without an active exception

I believe this is the problematic line that treats gradients > 0.05 as non zero gradients:

nonZeroGradientFlag.push_back( gradient_magnitude > 0.05F );

The threshold should be set at 50 as we do in all other cases--harmonization, TBSS etc.

---

Hi @pieper , is my interpretation of the F identifier i.e. 0.05F correct?

[pieper]

Hi @tashrifbillah -

Yes, the 0.05F would mean that the gradient would be considered zero if it's less than 5% of the base b-value. It looks like this calculation is being done on the normalized-then-scaled gradients directly out of the nrrd header, but really I would think you would want to convert them to s/mm^2 to have a consistent cutoff. That would mean taking into account the DWMRI_b-value tag. I'm not sure how that baseline b-value is selected in the converters.

The encoding is of course described here : )

https://www.na-mic.org/wiki/NAMIC_Wiki:DTI:Nrrd_format#Describing_DWIs_with_different_b-values

[tashrifbillah]

Yes, the 0.05F would mean that the gradient would be considered zero if it's less than 5% of the base b-value.

I am going to find out why this did not hold true for the above bvalues.

but really I would think you would want to convert them to s/mm^2 to have a consistent cutoff. That would mean taking into account the DWMRI_b-value tag.

Wouldn't absolute value 50 be a consistent cutoff? If yes, why would we need to take into account DWMRI_b-value tag?

[pieper]

Wouldn't absolute value 50 be a consistent cutoff? If yes, why would we need to take into account DWMRI_b-value tag?

It has to do with the way nrrd encodes gradients and b-values (see the wiki page linked above). There's one key that encodes the base b-value for the whole scan, and then the gradient directions are normalized to unit length and then if a particular gradient volume was acquired at a different b-value that direction vector is multiplied by the ratio of the current b-value to the base b-value. Because the selection of base b-value is arbitrary, the cutoff threshold is unitless and has no firm meaning. That's why Yogesh and I suggest converting to s/mm^2 to define the cutoff.

[tashrifbillah]

Okay, I understand. I was referring to tagging zero gradients using an absolute threshold of 50 before any normalization, in which case the DWMRI_b-value tag won't be required.

[tashrifbillah]

I have more details to share now. Yogesh and you (I believe) devised a way of doing UKF on bshells<2000 only. If I keep all bvalues of HCP data, zeros are tagged as:

5/3010= 0.0032 < 0.05

which works fine.

However, if I filter out the bvalues>2000, zeros are tagged as:

5/1510= 0.0016 < 0.05

In both cases, the supposed zero gradient magnitude remains below 0.05 but the latter fails :o

Now I wonder what you mentioned earlier:

I'm not sure how that baseline b-value is selected in the converters.

Thoughts?

[pieper]

Yogesh will need to say how best to handle the low b-value situation. I had thought the baseline scan would always have a gradient of 0,0,0. I'm not sure why a low but non-zero b-value would be considered zero.

Now I wonder what you mentioned earlier:

I'm not sure how that baseline b-value is selected in the converters.

Thoughts?

In most DICOM, the b-value is explicit per-slice in s/mm^2 so converters like dcm2niix need to pick one when writing the nrrd file. Mathematically it shouldn't matter except for possible rounding error, but in this case it needs to be taken into account.

[yrathi]

Hi @tashrifbillah -- I agree with Steve that the best solution is to set b-value less than 50 to be equivalent to zero. Now this value of 50 is arbitrary, but we will need to pick one threshold anyway. In the future we might have to change it, but for now 50 seems like a workable threshold.

[pieper]

Maybe change this to a command line parameter with a default value of 50 in case people need to experiment.

[tashrifbillah]

Okay, I did my analysis. Here is where it has been wrong so far:

nonZeroGradients: gradient_magnitude > 0.05F

It should be (note the square):

nonZeroGradients: gradient_magnitude² > 0.05F

NRRD bvalus are obtained as:

DWMRI_b-value * gradient_magnitude²

So, when we want to compare at the normalized level, the first term DWMRI_b-value is omitted, however, the second term must be retained with square.

That said, I should mention that I am not quite confident about 0.05 being a sufficient threshold. So comparison at absolute value level should be the safest approach:

nonZeroGradients: DWMRI_b-value * gradient_magnitude² > 50

Appendix

gradient_magnitude < 1 gradient_magnitude² < gradient_magnitude if gradient_magnitude² < 0.05, does not mean-- gradient_magnitude < 0.05

[tashrifbillah]

Maybe change this to a command line parameter with a default value of 50 in case people need to experiment.

This does not seem to be an easy work. I have been jumping between files in ukf/ folder to understand where the command line parameters are defined and how they are passed to worker functions.

I discovered the following so far:

https://github.com/pnlbwh/ukftractography/blob/4a4aba614f6337041e24c26522da0031d0132f6b/ukf/cli.cc#L33

https://github.com/pnlbwh/ukftractography/blob/4a4aba614f6337041e24c26522da0031d0132f6b/ukf/tractography.h#L27

If Steve can guide through the process of:

• defining a parameter: --bzero-thresh
• parsing that parameter
• accessing that parameter in dwi_normalize.cc

I would be happy to do that later. But for now, I am hardcoding 50 s/mm^2 to be the threshold.

[pieper]

@yrathi do you agree with @tashrifbillah that we should be comparing the magnitude squared? Maybe I don't recall the details but I understood the normalized gradients would be scaled by the ratio of the b-values.

Regarding the command line parameter I'm sure we can figure it out. Perhaps it's also time to do some other refactoring to make the code more maintainable. It is a question of how much development resources we want to spend on that compared to other priorities.

[tashrifbillah]

@yrathi do you agree with @tashrifbillah that we should be comparing the magnitude squared?

If not, then the threshold 0.05F has to be properly reduced.

[tashrifbillah]

nonZeroGradients: DWMRI_b-value * gradient_magnitude2 > 50

Hi @pieper , Yogesh and I reviewed the condition and concluded that it is correct. Can you please be a second set of eyes and review the PR #137 ? It is about 10 lines only!

[pieper]

If you guys agree that sounds good to me.