dwi2fod & sh2amp

[steso]

Not sure if it's really a bug or if the error is located somewhere else, I'm also happy to provide some example files if necessary. I generate two synthetic crossing fiber bundles based on a tensor in x and y directions (the identical tensor, just rotated 90°). The two branches have the same signal contribution (volume fraction) and there is no partial volume/up-sampling/etc involved. The diffusion signal is calculated on gradient scheme of 252 directions. If I perform an ODF reconstruction using the FRACT approach and compare the two peaks in the crossing region, I get more or less the same ODF-amplitude (<< 0.1% difference). In mrtrix I save the dwi.nii and the corresponding encoding.b file, run dwi2response & dwi2fod. If I extract the gradient amplitudes using sh2amp for 3 directions (x,y and z direction) I get about 5% difference btw. x and y. If I take a closer look at the branches (in the non-crossing areas) I also get some negative amplitudes, probably from negative side-lobes, which are also slightly different in the horizontal vs vertical bundle. If I use the SH-coefficients after dwi2fod and extract the amplitudes myself, I get the same values as the output from sh2amp. I'm not sure where this "directional" bias could have been introduced. There might be a small bias introduced by the initial sampling scheme, but 4-5% seems rather large to me. By the way I'm using version 0.3.10-455.

[jdtournier]

That does indeed sound fishy. Any chance you could post the exact commands you used?

[jdtournier]

It would also be great if you could post or send me the details of how the input was generated...

[steso]

Sure:

dwi2response "%path%/simCross_dwi.nii" "%path%/simCross_response.txt" -grad "%path%/encoding.b" -mask "%path%/simCross_mask.nii" -force -nthreads 16

dwi2fod "%path%/simCross_dwi.nii" "%path%/simCross_response.txt" -grad "%path%/encoding.b" "%path%/simCross_csd.nii" -mask "%path%/simCross_mask.nii" -nthreads 16

gendir 3 dirs.txt sh2amp simCross_csd.nii dirs.txt amps.nii

I also attached some files in this message hope I got everything and the correct files, otherwise, just ask! There might still be the possibility of a bug between my matlab scripts and mrtrix but I can't think of any possible way to introduce a bias like that. I also chose [1 0 0] and [0 1 0] to even exclude possible flips and sign errors.

I also tried to convert the sh coeffs to the old basis but I still got an offset between x and y

On 31.03.2015 17:20:13, J-Donald Tournier notifications@github.com wrote: It would also be great if you could post or send me the details of how the input was generated... — Reply to this email directly or view it on GitHub [https://github.com/jdtournier/mrtrix3/issues/201#issuecomment-88129785].

[jdtournier]

OK, had a look into this, the first thing that struck me is that the FOD had been estimated up to harmonic order 20 (231 coefficients)! I didn't even know you could push it that far. MRView will only display up to order 16...

So it looks like I definitely need to change the defaults, there is little point in my experience in pushing lmax beyond 8 or 10. I'm guessing FRACT won't be going beyond 6?

So the reason behind the discrepancy is due to the fact that the constraint is applied along 300 directions. When the number of coefficients gets close to this number, you can imagine that it becomes possible for some of the high frequencies to become a little less well behaved, and this probably explain most of the differences you observed.

I'm tried running dwi2fod at a more sensible lmax=8, and that reduces the bias down to <1%. It's still not perfect, but definitely well below the noise floor on anything but the most extreme acquisitions... Nonetheless, there may be some non-uniformity in the distribution of the constraint directions that might be causing this remaining difference, despite my best efforts at optimising these directions - we recently observed that using a higher power for the repulsion term in an electrostatic repulsion optimisation for the directions, as used in dirgen (or gendir as it used to be called) might be counter-productive, despite producing wider nearest-neighbour angles... I'll investigate tomorrow, see if we can improve on that.

[steso]

Thanks for your quick response. I was also astonished by the number of coefficients, in my code I had a hardcoded Lmax of 8. Yes for the FRACT reconstruction I typically use SH-orders up to L=6. For in-vivo datasets (typically 64dirs) I only got up to order 8 in mrtrix if I remember correctly. Maybe the the 252 gradient directions in the simulation pushed it up to order 20. What I tried was just using the first few coefficients up to order 8 (44-ish) for the conversion from SH to sampling points, but that didn't help for the bias. So I'll try to limit the max SH order in my dwi2fod call!

[steso]

Ok I just tried dwi2fod with -lmax 8 and everything is below 1% now, FRACT seems to be smaller but I agree this will be below noise-level. About the repulsion, I don't know about the current implementation of dirgen but it might help to optimize the points such that all the points including the point-symmetrical counterparts are optimized in the repulsion due to symmetry in even order SH? So even if n-points might be "optimally" distributed, the 2*n points due to point-symmetry might not be?

[Lestropie]

I would be curious to know whether the percentile difference between FOD lobe integrals is less than that of the peak amplitudes. It's possible that one lobe may be fractionally more dispersed than the other, decreasing the peak amplitude, but the actual fibre density estimates are still close to equal.

About the repulsion, I don't know about the current implementation of dirgen but it might help to optimize the points such that all the points including the point-symmetrical counterparts are optimized in the repulsion due to symmetry in even order SH? So even if n-points might be "optimally" distributed, the 2*n points due to point-symmetry might not be?

Yes, this is already handled appropriately.

[jdtournier]

I would be curious to know whether the percentile difference between FOD lobe integrals is less than that of the peak amplitudes.

Yes, I suspect you might be right there. Still, ideally there wouldn't be such a discrepancy...

[steso]

I spent some time looking at the lobe-integrals by running fod2fixel -afd I also had some troubles opening the sparse-image file in matlab, but I wrote a quick wrapper now. Just to make sure I got the correct results, I interpreted the 20 bytes per sparse-class entry the following way ('N2MR5Image6Sparse11FixelMetricE'): 3x4 bytes float for x,y,z direction 1x8 bytes double for the lobe integral is this correct?

Not sure if you'll like my results but that's what I got:

amplitude, L=20: 4.4963% difference
amplitude, L=8: 0.3542% difference

(what we discussed before)

afd, L=20: 0.5340% difference between x&y lobe afd, L=8: 2.3436% difference between x&y lobe

for the lobe peaks I got the following directions (x,y,z):

L=20: 0.0002 1.0000 0.0003 1.0000 0.0002 0.0001

L=8: -0.0001 1.0000 -0.0001 1.0000 -0.0002 0.0004

[jdtournier]

Ahem, yes, that's not as fantastic as we'd like it to be...

Nonetheless, given that this is a constrained problem, and distributing the constraints uniformly is a non-trivial problem, it's not unexpected. I'm having a go with different directions for the constraints as we speak...

By the way, that gendir 3 dirs.txt trick will not work on recent versions. I've removed the constraint that the first direction lines up with the Z-axis, and that the second direction be in the XZ plane - this is to produce a more genuinely random distribution when generating sampling schemes. And the command is now dirgen, in line with our naming conventions...

[jdtournier]

OK, I've made a couple of commits to address this issue: dwi2fod now defaults to lmax=8 (commit e17a02a), and the directions used for the constraint have been regenerated using an r² electrostatic repulsion term (commit 0608f7b). This seems to reduce the issue quite a bit, at least for the FOD amplitudes (haven't tested integrals). On your data:

lmax 8:

before: 2.20793 vs 2.22107 (~0.6%) after: 2.22018 vs 2.22211 (> 0.1%)

lmax 20:

before: 2.51444 vs 2.63319 (~4.7%) after: 2.58067 vs 2.58574 (> 0.2%)

I guess we can call this fixed...?

[steso]

I just updated and checked again, not sure why we don't get the same values for the absolute amplitudes, at least the percentages and range are quite similar:

lmax 8: 2.2699 2.2731 -> 0.14% lmax 20: 2.5310 2.5363 -> 0.21%

for the integrals it's still not ideal (are these values reasonable?): lmax 8: 0.0306 0.0314 -> 2.44% lmax 20: 0.0223 0.0225 -> 0.97%

Yes I'd say it's fixed. The lobe integral might also depend on the sampling-scheme in the AFD computation, I guess you're using 1281 icosahedron as described in the SIFT-paper?

[jdtournier]

The remaining discrepancies are probably down to where the values were measured - I measured away from the crossing region, I noticed the values in the crossing region were slightly different.

Otherwise, I'd agree the integral value discrepancies are most likely down to the direction set used to sample the FOD. The tessellation strategy doesn't give perfectly uniform direction sets... We could look into substituting an equivalent electrostatic repulsion set - @Lestropie, what do you think?

Finally, thanks for reporting this, I probably wouldn't have looked into this otherwise. Hopefully that'll be a good outcome for everyone :)

[Lestropie]

lmax 8: before: 2.20793 vs 2.22107 (~0.6%) after: 2.22018 vs 2.22211 (> 0.1%) lmax 20: before: 2.51444 vs 2.63319 (~4.7%) after: 2.58067 vs 2.58574 (> 0.2%)

Damn...

Yes I'd say it's fixed. The lobe integral might also depend on the sampling-scheme in the AFD computation, I guess you're using 1281 icosahedron as described in the SIFT-paper?

Yes, the 1281 scheme is the default, and the result will likely depend on the sampling scheme. Especially at lmax 20: not as many sampling points per lobe.

Otherwise, I'd agree the integral value discrepancies are most likely down to the direction set used to sample the FOD. The tessellation strategy doesn't give perfectly uniform direction sets... We could look into substituting an equivalent electrostatic repulsion set - @Lestropie, what do you think?

The tessellation strategy is used in this context because the segmentation algorithm requires a robust definition of direction adjacency. With electrostatic repulsion, even eye-balling the data you can find cases where this is ambiguous. We could try a higher density scheme, or I could look into a better integral calculation method e.g. equivalent of trapezoidal rule on the sphere (at the moment it's just a sum of amplitudes)

[Lestropie]

Also:

Just to make sure I got the correct results, I interpreted the 20 bytes per sparse-class entry the following way ('N2MR5Image6Sparse11FixelMetricE'): 3x4 bytes float for x,y,z direction 1x8 bytes double for the lobe integral is this correct?

Not quite; it's 3x4 bytes for xyz direction, 4 bytes for amplitude (AFD), and 4 bytes for the 'associated value' (in the specific case of fod2fixel -afd, it's also equal to the AFD). https://github.com/jdtournier/mrtrix3/blob/master/lib/image/sparse/fixel_metric.h

[jdtournier]

The tessellation strategy is used in this context because the segmentation algorithm requires a robust definition of direction adjacency.

Ah yes, I'd forgotten that detail... What we could use if we care about this enough is some form of variable density weighting to correct for any residual non-uniformity. It's used a fair bit for regridding of non-Cartesian data in image reconstruction, the idea is essentially what we want here: each sample gets assigned a weight that reflects the area it represents. Not too sure how we go about estimating these weights though, but I'm sure there are strategies for this out there...

[steso]

Ok so with the corrected read-function and the new code-contribution I get the following values:

amplitudes in the branches: L=8: 2.2202 2.2221 -> 0.0872% L=20: 2.5807 2.5857 -> 0.1959% amplitudes in the crossing: L=8: 2.2731 2.2700 -> 0.1396% L=20: 2.5363 2.5310 -> 0.2081%

afd in the branches (values should be correct now!!) L=8: 1.2505 1.2449 -> 0.4476% L=20: 1.1802 1.1785 -> 0.1476% afd in crossing L=8: 1.2474 1.2473 -> 0.0057% L=20: 1.1726 1.1736 -> 0.0901%

==> problem solved, also for AFD!

[jdtournier]

Guess we don't need to worry about non-uniformity of the tessellated directions after all then...

@Lestropie, feel free to close if you're happy to leave the integral calculations as-is - seems good enough to me...

[Lestropie]

What we could use if we care about this enough is some form of variable density weighting to correct for any residual non-uniformity.

That was the first idea; trouble was calculating the appropriate density weighting. I think there's more likely to be a closed-form solution for the solid angle of a spherical triangle, rather than the effective solid angle of each point in the set. And the direction adjacencies inherently provide the triangulation.

@Lestropie, feel free to close if you're happy to leave the integral calculations as-is - seems good enough to me...

But.. but...

[jdtournier]

But.. but...

He he... What I meant is it's probably close enough not to introduce significant biases. But sure, ideally we'd make sure that it was all as accurate as possible. If you want to do something about it, go for it.

Having thought about it a little, I reckon you could compute the weights relatively easily. Start with a vector of SH coefficients, where you set the l=0 term to 1, and randomise all the other coefficients. By construction, the integral over the sphere is a known constant, since the l>0 terms have zero integral. If you now compute the amplitudes for those coefficients along your set of directions, and compute the weighted sum of these amplitudes, you should ideally get the same constant value (or some multiple of it). More formally, you have a vector s1 of n SH coefficients, multiplied by the m by n matrix M mapping the n coefficients onto the vector of m directions, so the weighted sum is just the dot product of those amplitudes with their respective weights w:

w' M s1 = c

where w' is the transpose of w. Equivalently, taking the transpose of that equation:

s1' M' w = c

This should hold for any vector of SH coefficients with its l=0 term set to 1, so we can form a p by n matrix S , where each row is a vector of SH coefficients with its l=0 term set to 1, and all other coefficients are randomised, and we can chose p to be as large as we want - large enough that (S M') is invertible, so that you can estimate w by solving:

S M' w = c

where c is a vector with all its elements set to c (i.e. a constant, use 1 if you want, it'll just scale the weights accordingly). In practice, you'll want to use a value of lmax for those SH vectors that's a fair bit higher than would ever be used for the FOD, maybe 20 or 30. Even then, I'm guessing the matrix (S M') will be ill conditioned - there's not much you can do about the fact that the matrix rank can't exceed the number of SH coefficients. You'll need to get the minimum norm solution somehow - probably easiest is to form the Moore-Penrose pseudo-inverse for full row rank matrices, using a square matrix for S to ensure (S M') has full row rank. You might be able to achieve the same thing with an SVD if you feel like it, or just add a tiny amount of explicit minimum norm regularisation - i.e. solve:

(M S' S M' + lambda I) w = c

where lambda is very small, 1e-10 or something (relative to the diagonal of the other matrix).

What do you reckon? Make sense...?

[Lestropie]

Closing as the fundamental issues raised have been addressed; additional enhanced functionality moved to #312.