muscbridge / PyDesigner

A hands-free DTI, DKI, FBI and FBWM preprocessing pipeline. Information on algorithms and preprocessing steps are available at https://www.biorxiv.org/content/10.1101/2021.10.20.465189v1 A video tutorial on PyDesigner and its usage is now available at https://www.youtube.com/watch?v=mChQFuQqX3k
https://pydesigner.readthedocs.io/en/latest/
Other
22 stars 7 forks source link

Fitting error #265

Closed saraponti closed 2 years ago

saraponti commented 2 years ago

Hi all, I have I problem running pydesigner (linux ubuntu 18.04, python 3.6.9, fsl version 6.0.5.1, mrtrix 3.0.1). I have already preprocessed my data, so I just want to fit the kurtosis tensor on multishell data. The command I used is:

pydesigner -o /media/sara/Maxtor/Naples/Processing/gioesp_proc/corr_data/pydesigner/output --user_mask /media/sara/Maxtor/Naples/Processing/gioesp_proc/corr_data/pydesigner/nodif_brain_mask.nii --nthreads 10 --verbose /media/sara/Maxtor/Naples/Processing/gioesp_proc/corr_data/pydesigner/data.nii

After denoising and outlier detection, the fitting process gave me the following error. Do you have any suggestion to fix this?

Thanks, Sara

Constrained Tensor Fit: [20%]joblib.externals.loky.process_executor._RemoteTraceback: """ Traceback (most recent call last): File "/home/sara/.local/lib/python3.6/site-packages/joblib/externals/loky/process_executor.py", line 436, in _process_worker r = call_item() File "/home/sara/.local/lib/python3.6/site-packages/joblib/externals/loky/process_executor.py", line 288, in call return self.fn(*self.args, *self.kwargs) File "/home/sara/.local/lib/python3.6/site-packages/joblib/_parallel_backends.py", line 595, in call return self.func(args, *kwargs) File "/home/sara/.local/lib/python3.6/site-packages/joblib/parallel.py", line 263, in call for func, args, kwargs in self.items] File "/home/sara/.local/lib/python3.6/site-packages/joblib/parallel.py", line 263, in for func, args, kwargs in self.items] File "/home/sara/.local/lib/python3.6/site-packages/designer/fitting/dwipy.py", line 812, in wlls objective = cvx.Minimize(0.5 cvx.sum_squares(C @ x - d)) File "/home/sara/.local/lib/python3.6/site-packages/cvxpy/expressions/expression.py", line 48, in cast_op other = self.cast_to_const(other) File "/home/sara/.local/lib/python3.6/site-packages/cvxpy/expressions/expression.py", line 472, in cast_to_const return expr if isinstance(expr, Expression) else cvxtypes.constant()(expr) File "/home/sara/.local/lib/python3.6/site-packages/cvxpy/expressions/constants/constant.py", line 58, in init super(Constant, self).init(intf.shape(self.value)) File "/home/sara/.local/lib/python3.6/site-packages/cvxpy/expressions/leaf.py", line 108, in init raise ValueError("Invalid dimensions %s." % (shape,)) ValueError: Invalid dimensions (0, 22). """

TheJaeger commented 2 years ago

Hi Sara,

What are dimensions of your DWI (data.nii)?

saraponti commented 2 years ago

Hi, my data.nii is a 4D of dimensions 128x128x68x129 (I have 3 shells: 15 b=0, 6 dir b=500, 48 dir b=1000 and 60 dir b=2000)

TheJaeger commented 2 years ago

This is a weird error that I've never seen before, and occurs somewhere during tensor fitting regime. What's baffling is that it occurs after it's 20% done. If you don't mind, could you share this dataset (deidentified) with me at siddhartha.dhiman@gmail.com so I could debug PyDesigner and figure out what's causing it?

-Sid

saraponti commented 2 years ago

Hi Sid, I sent you the data! Anyway, we also have the same error (but fitting stopping at 1% instead of 20%) using a different dataset on a different pc.

Sara

TheJaeger commented 2 years ago

Hi Sara, could you also list the preprocessing steps you've ran on this dataset?

TheJaeger commented 2 years ago

The problem occurs because tensor estimation doesn't have sufficient directions to converge to a solution. A diffusion tensor (DT) requires at least 5 good directions, while the kurtosis tensor (KT) requires 15. I reckon that IRLLS marks so many directions at certain voxels as outliers that less than minimum directions are present for tensor estimation.

Making the IRLLS algorithm less aggressive fixed the problem on my end, and I'll add this change into the next public release. In the meantime, you can disable outlier detection with the flags --nooutliers and --noakc. The full command to parse would then be:

pydesigner --nooutliers --noakc -o /media/sara/Maxtor/Naples/Processing/gioesp_proc/corr_data/pydesigner/output --user_mask /media/sara/Maxtor/Naples/Processing/gioesp_proc/corr_data/pydesigner/nodif_brain_mask.nii --nthreads 10 --verbose /media/sara/Maxtor/Naples/Processing/gioesp_proc/corr_data/pydesigner/data.nii
saraponti commented 2 years ago

Hi Sid! I checked the outlier map and noticed that almost all the detected outlier voxels were outside the brain. I noticed that I have 32 voxels marked as outliers for all 114 directions, while the other outliers voxels were marked as outliers for less then ~60 directions so the fitting procedure should had enough information for both DT and DK estimation. My question is: how the outlier voxels are treated during the fitting? Aren't they masked and then excluded from the fitting?

Anyway, I decided to try your suggestion and run the fitting procedure without outlier detection (--nooutliers and --noakc). The fitting ended successfully, but the maps seem a bit strange, as misalignment between data and input brain mask occurs whereas the original data were correctly aligned. Here attached an example of the FA map with the original brain mask overlaid in yellow. Thanks a lot for your support. Sara

Screenshot from 2021-11-26 12-07-54

TheJaeger commented 2 years ago

Hi Sara, outlier voxels are completely excluded from fitting. By right, PyDesigner would evaluate the tensor value close to zero if insufficient directions are present or if the fitting is not an optimal solution. It doesn't appear to be doing that with this dataset for some reason, and instead resorts to yielding an error. The internal fix I have excludes B0s from IRLLS outlier detection, so there is always some minimum data present to prevent the error.

I noted the shift in FA map too when I ran PyDesigner without its preprocessing. This is happening because the mask and DWI are in LAS and RAS orientation, respectively. PyDesigner does not check the orientation of user supplied masks to flip them accordingly so you'd need to ensure that both of them have the same orientation. I will add a check for this in the next release. Running the PyDesigner pipeline entirely (preprocessing + tensor estimation with outlier detection) corrects the issue because it computes and write the mask in the correct orientation. The FA map shown below does not possess the same shift.

Screenshot from 2021-11-29 10-14-11

TheJaeger commented 2 years ago

Hi Sara, the latest release v1.0-RC11 fixes the issues you're facing.