DTI/DKI RESTORE using the CLI

[araikes]

Description

There isn't currently a direct way to compute the sigma value in the CLI framework for dipy_fit_dti or dipy_fit_dki with the method as RESTORE. Several noise estimation sigma options exist in DIPY, including:

1. dipy.denoise.noise_estimate.estimate_sigma: This is the method presented in the DTI RESTORE tutorial. However, this is not forward facing for the CLI and returns a 1D array with an estimate for each volume in the 4D DWI data. If not built directly into the DTI/DKI CLI workflows, this would need to be computed and passable as an array.
2. Directly on the CLI, it is possible to use the MP-PCA denoising method, return the sigma (3D image) and compute e.g. the mean value to pass to the CLI.

Currently Option 1 isn't feasible as, per @skoudoro on Gitter, --sigma is expecting a singular value. I tried Option 2 but got back an error message:

INFO:Computing DTI metrics for sub-control01_desc-preproc_dwi.nii.gz
INFO:Tensor estimation...
Traceback (most recent call last):
  File "/usr/local/miniconda/bin/dipy_fit_dti", line 7, in <module>
    run_flow(ReconstDtiFlow())
  File "/usr/local/miniconda/lib/python3.9/site-packages/dipy/workflows/flow_runner.py", line 89, in run_flow
    return flow.run(**args)
  File "/usr/local/miniconda/lib/python3.9/site-packages/dipy/workflows/reconst.py", line 342, in run
    tenfit, _ = self.get_fitted_tensor(data, mask, bval, bvec,
  File "/usr/local/miniconda/lib/python3.9/site-packages/dipy/workflows/reconst.py", line 417, in get_fitted_tensor
    tenfit = tenmodel.fit(data, mask)
  File "/usr/local/miniconda/lib/python3.9/site-packages/dipy/reconst/dti.py", line 789, in fit
    params_in_mask = self.fit_method(
  File "/usr/local/miniconda/lib/python3.9/site-packages/dipy/reconst/dti.py", line 1858, in restore_fit_tensor
    this_param, status = opt.leastsq(_nlls_err_func,
  File "/usr/local/miniconda/lib/python3.9/site-packages/scipy/optimize/minpack.py", line 414, in leastsq
    raise TypeError(f"Improper input: func input vector length N={n} must"
TypeError: Improper input: func input vector length N=7 must not exceed func output vector length M=2

I don't know if that error is related to the --sigma parameter or not, but the same data and same mask using default dipy_fit_dti methods computes the tensor without issue.

[If reporting a bug, attach the entire traceback from Python and follow the way to reproduce below]
[If proposing an enhancement/new feature, provide links to related articles, reference examples, etc.]

Way to reproduce

[If reporting a bug, please include the following important information:]

• Code example: singularity exec /groups/adamraikes/singularity_images/cibs_neurotools_1.0.15.sif dipy_fit_dti sub-control01_desc-preproc_dwi.nii.gz sub-control01_desc-preproc_dwi.bval sub-control01_desc-preproc_dwi.bvec sub-control01_desc-dsurqe_mask.nii.gz --fit_method 'RESTORE' --sigma 10.8819
• Relevant images (if any)
• Operating system and version Ubuntu 20.04
• Python version Python 3.9.7
• dipy version 1.5.0

[arokem]

TypeError: Improper input: func input vector length N=7 must not exceed func output vector length M=2

This is usually an indication of a sigma value that is too small. It means that so many values were marked as outliers in this particular voxel, that there are not enough values left to compute the tensor.

Not solving the issue you are reporting, we should also really fix the error handling here, because that error message confuses many users.

[araikes]

@arokem,

Thanks for the quick response yesterday and that explanation makes sense about the error message. The data were preprocessed using patch2self, so the expected noise level, particularly in non-brain/head areas (note: full head rodent images in my particular, current application) is small and therefore making the mean of the MP-PCA sigma image file in this example very small.

That said, any recommendations on calculating sigma for RESTORE using preprocessed data, particularly using the CLI? I see the following challenges:

1. Clearly the mean sigma from MP-PCA denoising on the preprocessed data is too small.
2. The current docs for dipy_fit_dti and dipy_fit_dki suggest that sigma should be 1.5267 * std(background_noise). Looking at my example image and knowing where there's a non-brain "background" region (no brain, no skull, no face, just dead space), this would be even smaller than the value passed from the mean of the whole image.
3. I tried the mean of the sigma image in the brain mask I'm using for tensor fitting. That doubled the sigma value but yielded the same problem.
4. Given this, I did a little brute forcing just for interest's sake and found that --sigma 65 was the smallest sigma that found a solution. This is 1.696 * max(sigma_image) but I doubt that this is a reliable calculation across all images.
5. Interestingly --sigma 65 and --sigma 100 produce identical FA images (but different from non-RESTORE FA), suggesting that a minimum sigma that solves is sufficient, at least in my single test use case.

Thanks for any thoughts/insights.

[ShreyasFadnavis]

@araikes -- This is very interesting (at least to me)! Can you share some before after pictures so that we get a better sense of what's happening?

Using Patch2Self there are 2 alternatives that can be used as sigma estimates:

1. You can use NUQ scores (with some scaling) as sigma within RESTORE-DTI. -->NUQ Paper
2. One could also use Leverage Scores computed on the data pre-Patch2Self as sigma scores for RESTORE-DTI. --> P2S Sketching Paper

These are not merged in DIPY yet, but I am happy to share the code and test it out 👍🏽

[araikes]

Hi @ShreyasFadnavis

Those are interesting options and I'd be happy to try either of them out to see if they are useful in this effort.

As far as pictures go, what would you like to see? Obviously when I hit that error message, there's no picture to share since nothing gets produced. Other than that, CLI options with --sigma 65 or --sigma 100 produce almost identical images to a non-CLI approach using dipy.denoise.noise_estimate.estimate_sigma

[ShreyasFadnavis]

Hi @araikes -- Sorry for the slowness in my response! I think it may be a good idea to try doing:

Raw Data -> DTI fit via 'OLS' model -> FA

VS

Raw Data -> Patch2Self -> DTI fit via 'OLS' model -> FA

The reason I am asking to use OLS and not RESTORE or WLS (default) is: it will reflect clearly how much outliers are left in the data. AFAIK: Patch2Self should suppress most of them and you may not need to do RESTORE at all :)

Thanks,
Shreyas

[araikes]

@ShreyasFadnavis,

Sure, I can do that but I want to clarify within my workflow. Here's the workflow to this point:

1. Denoising (Patch2Self)
2. Gibbs unringing (MRtrix)
3. Initial masking with template mask (ANTs)
4. Eddy correction (TORTOISE or MRtrix, depending on data source)
5. Bias correction (MRtrix)
6. Tensors and scalars

Happy to shorten for test cases here but let me know if there are any other images you want within that workflow.

[araikes]

Models fit were dipy_fit_dki .... --save_metrics "fa" --fit_method 'OLS' for raw_FA and patch2self_FA. sub-control01_FA was the fully preprocessed image fit with DKI and WLS. sub-control01_model-dkiRESTORE_FA was fit using a Python script using dipy.denoise.noise_estimate.estimate_sigma to model the noise.

[ShreyasFadnavis]

Hi @araikes ! Sorry for the slowness in my response. As far as I can see, Patch2Self is suppressing a lot of the outliers in your data and you may not need RESTORE. Just as an experiment, can you also try to switch the order of your pipeline:
Raw Data -> Mask -> Eddy -> Patch2Self -> Gibbs -> Bias Corrections -> Tensors and Scalars?

[araikes]

@ShreyasFadnavis

Yes I will. I'll try and get it done tomorrow. Thanks