tckgen -algorithm tensor_*: Use tensor input image

[Lestropie]

Currently the tensor tracking algorithms do an unweighted linear least squares fit on the trilinear-interpolated diffusion data to obtain the tensor at each scanner space point.

Instead, the tensor_det algorithm could receive as input a pre-calculated tensor volume, and perform appropriate interpolation using tensor calculus.

The tensor_prob algorithm would still receive as input the DWI volume, but the wild bootstrap would write to its internal buffer the resulting tensor fit rather than the bootstrapped DWI data.

This would have the potential to speed up tensor tracking particularly for #534. And would probably make more sense from a user perspective; a number of people have gotten confused after providing the wrong input image for the tensor_det algorithm.

[SyamGadde]

If, by chance, there is a simple way to implement this, could you give a pointer to any potential pitfalls if I decided to try it myself? We want to evaluate an alternative tensor calculation ( DiCIPHR-Lab/Fernet ) by comparing it to existing tckgen outputs, but tckgen wants to do the calculation itself, and I'd love to not have to revert to tracking tools from other packages (though I can if I need to!)

[Lestropie]

There's a couple of different ways to do it. If it were me, then instead of modifying the existing tensor_det source code to accept either a DWI or pre-calculated tensor input image, I would instead create a new tracking algorithm that is designed to operate on the latter. The argument over what each algorithm should be called would need to be deferred, but this would be easier for testing changes, and would also execute faster due to less conditional execution.

To make a new algorithm:

1. Copy one of the files in src/dwi/tractography/algorithms/ (tensor_det.h would make the most sense) and give it a unique name;
2. Give the class defined in that file a new name;
3. In cmd/tckgen.cpp, add the new algorithm name to the list of algorithms (above the usage() function), and add the corresponding Exec<>::run() call to the end of the run() function.
4. Re-compile.

tckgen should now provide you with access to that new algorithm. It will behave identically to the algorithm you copied, but it serves as a starting point to then be able to modify the code to behave as you desire and test it without breaking any existing functionality.

A few thoughts that may be relevant:

1. The tensor_* tracking algorithms perform an unweighted least-squares tensor fit, unlike the dwi2tensor command, which by default performs an iterative weighted least-squares fit.
2. You should ideally employ the appropriate calculus in order to interpolate the tensor at sub-voxel locations, and not simply perform linear interpolation on the tensor coefficients.
3. A tensor fit (as with all other anisotropic diffusion models) is defined with respect to a reference set of axes. In MRtrix3, we define all orientation information with respect to "scanner space": the physical axes of the scanner hardware. If your software performs calculations with respect to e.g. the image axes, then the appropriate transformation will need to be applied, otherwise the orientations will be misinterpreted.

[SyamGadde]

As it turns out I already started along exactly those lines, but realized that the tensor interpolation would require some further thought. Good to know that it should be possible. Thank you for the very helpful advice!

On Mon, Jan 31, 2022, 6:22 PM Robert Smith @.***> wrote:

There's a couple of different ways to do it. If it were me, then instead of modifying the existing tensor_det source code to accept either a DWI or pre-calculated tensor input image, I would instead create a new tracking algorithm that is designed to operate on the latter. The argument over what each algorithm should be called would need to be deferred, but this would be easier for testing changes, and would also execute faster due to less conditional execution.

To make a new algorithm:

1. Copy one of the files in src/dwi/tractography/algorithms/ ( tensor_det.h would make the most sense) and give it a unique name;
2. Give the class defined in that file a new name;
3. In cmd/tckgen.cpp, add the new algorithm name to the list of algorithms (above the usage() function), and add the corresponding Exec<>::run() call to the end of the run() function.
4. Re-compile.

tckgen should now provide you with access to that new algorithm. It will behave identically to the algorithm you copied, but it serves as a starting point to then be able to modify the code to behave as you desire and test it without breaking any existing functionality.

A few thoughts that may be relevant:

1. The tensor_* tracking algorithms perform an unweighted least-squares tensor fit, unlike the dwi2tensor command, which by default performs an iterative weighted least-squares fit.
2. You should ideally employ the appropriate calculus in order to interpolate the tensor at sub-voxel locations, and not simply perform linear interpolation on the tensor coefficients.
3. A tensor fit (as with all other anisotropic diffusion models) is defined with respect to a reference set of axes. In MRtrix3, we define all orientation information with respect to "scanner space": the physical axes of the scanner hardware. If your software performs calculations with respect to e.g. the image axes, then the appropriate transformation will need to be applied, otherwise the orientations will be misinterpreted.

— Reply to this email directly, view it on GitHub https://github.com/MRtrix3/mrtrix3/issues/535#issuecomment-1026310401, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAJQRZDELFIAXQXVHU5KOETUY4KUXANCNFSM4CAX7CNQ . Triage notifications on the go with GitHub Mobile for iOS https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675 or Android https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub.

You are receiving this because you commented.Message ID: @.***>

[Lestropie]

By my reading from a long time ago, tensor interpolation should be done via this method; not sure if there's anything proposed more recently. Longer term this would ideally be implemented as a novel interpolator in core/interp/, though we could deal with code rearrangement after the point at which you have something functioning.

[bjeurissen]

I have used https://pubmed.ncbi.nlm.nih.gov/16788917/ to tackle this in https://pubmed.ncbi.nlm.nih.gov/24968247/.

It works by taking the (precalculated) matrix log from your tensors and performing the interpolation on the log tensors (the matrix log amounts to taking the log of the eigenvalues of the matrix and leaving the eigenvectors as is). Once you have the interpolated log tensors it is then a cheap operation to go back to the tensors.

All the pieces to do this are already available within MRtrix, you just have to put them together and do all the bookkeeping… Happy to assist if needed!

[SyamGadde]

Thanks! I did some exploration with the Batchelor reference above from @Lestropie but was getting complex values in the interpolated tensors, which may or may not be errors on my part, but the above Log-Euclidean work may make it even simpler. However, we've discovered that we probably don't need to run tractography for these data, so I may not be spending too much time on this. If I do come back to it and it is as yet unimplemented, I'll definitely take a closer look at these.