All the attached images are difference images that should be 0 if things were working correctly.

Interesting - good thing you spotted this!

We should include Raul and possibly also Gordon in the investigation depending on what we find. I haven't look at the details of the mathematics of the tensor estimation implementation, but I'm pretty sure it is calling the teem code (equivalent of 'tend estim' if I recall correctly) but wrapped in a multithreaded VTK class.

The banding that you see is almost certainly related to the 'slabs' that are allocated to different threads (compute cores) when VTK parallelizes the calculation.

You might see this kind of thing when there's a random number and the different threads get different seeds. But I believe that DWI to DTI estimation is meant to be deterministic.

Instead I would think that it's some kind of thread safety issue in either teem or vtkTeem that causes some cross-talk between threads.

The easiest workaround to test would be turning threading off for the estimator.

Beyond that, for debugging I suggest making a simple test script that demonstrates the issue and can be used to confirm when it is fixed.

Thanks Steve. I thought it was a multithread bug at first also, but as we re-run it on the same dataset, the bands move apparently randomly, so it is still potentially involving multithread but if so I am surprised that the location of "edges" between thread compute regions are not obviously preserved across runs. This may be hard to tell visually though due to the random element.

Yogesh, Ofer, and I have tested more and found that the magnitude of the variability is very small with "traditional" b=1000 diffusion-weighted data, whereas with newer multishell data, the magnitude is much larger and these bands can be seen.

Another thing I've noticed is that the bands appear in the "baseline" computed from DTI estimation but not in the "baseline" computed from tensor mask estimation, so something different is happening in these two pipelines. Do you know if the mask module is multithreaded?

To summarize, yes I agree there is potentially a multithread component, but it is behaving unusually, and the magnitude of the variability seems to be highly related to the type of data input, which may help track things down somehow.

I took a quick look at all the code and I don't see anything obvious in the vtkTeem or teem (ten) code.

There is a lot of floating point math done per-pixel, so it could be something odd with rounding.

Further testing. In b=1000 data the difference from one calculation to the next is probably insignificant. For example, in some voxels, the FA changed by 0.000001.

So this bug seems to primarily affect multishell data.

This DTI estimation issue is present in the release of Slicer 4.4.0 from 11-02-2014. So it may not be related to the teem upgrade that was incorporated into Slicer in approximately 1/2015.

I've found that the same bug behavior (tensor estimation is not reproducible) is present in Slicer 3.6.3 and 4.0.0, so it is not caused by vtk 6.0. I've learned from Julie that using tend outside of slicer gives different results from computing tensors inside Slicer (which also uses the teem library), so the vtk interface to teem/tend is the likely culprit at this point.

Since the issues are greatly magnified with multiple b-values including ones higher than 1000, the next place to look is the handling of the b-values inside the vtk classes that interface with teem. There are some comments in this code region about thread safety and computing a maximum b-value in a thread safe region, so I suspect something along these lines, especially due to the striping of the differences in (for example) baseline images computed repeatedly.

Did you look at this block of code?

https://github.com/Slicer/Slicer/blob/master/Libs/vtkTeem/vtkTeemEstimateDiffusionTensor.cxx#L509-L525

It looks to me like this code is being called from the threads but it manipulating data that is common across threads. That is, it's called from SetTenContext, which is called from vtkTeemEstimateDiffusionTensorExecute, which is called from ThreadedExecute, but it's manipulating this->DiffusionGradients, which is a shared instance variable.

What is the status on this ?

Was thinking to update teem version used in Slicer to include the changes that Hans and Raul recently committed to teem SVN. See https://github.com/Slicer/Slicer/pull/262

Let me know if that would interfere with your investigation ?

This is independent of teem updates, so it's great if you update Slicer to include this fix as well as teem version updates.

Would be great if you could spend some time to add a test that:
(a) would fail without your patch [1]
(b) pass after your patch is included

To make things easier on your side, you could simply attach a batch file (or python file) implementing this test. This file could invoke a CLI two times in a row and then use cmake -E compare_files /path/to/image1 /path/to/image2 to compare the results.

[1] https://github.com/ljod/Slicer/commit/ca5080bd1c90cf7e941a2ec7703ed3a3800ce015?w=1

Fix [1] integrated in r24240 [2]

[1] https://github.com/Slicer/Slicer/pull/277

[2] http://viewvc.slicer.org/viewvc.cgi/Slicer4?view=revision&revision=24240

As a side note, using the input data "helix-DWI.nhdr" [1] and following the instruction. Subtracting the two FAs shows some difference. See screenshots "Slicer-24240_Issue_3977_FA_are_different_for_helix_dwi"

@Lauren: Is that expected ?

[1] https://github.com/Slicer/Slicer/blob/master/Libs/MRML/Core/Testing/TestData/helix-DWI.nhdr

The different FA bothers me - see my comment on pull request 277.

A self test using @pieper approach has been added. See r24248

As a side note, using the Subtract module and setting the interpolation order to 0 (instead of the default value being 1), the difference of the two FAs is zero.

See https://github.com/Slicer/Slicer/pull/277#issuecomment-101553810