Fiber Tracking: FBI tractography (#235)

* Fix DTI protocol (#228)

* Update reslicing to use new `mrgrid` command

This also works with older `mrreslice`

* Fix bad indent in tensorReorder()

* Update documentation

* Add bval scaling (#230)

* Fix SNR plotting issues (#233)

* Add try/except loop to SNR plot

* Fix eddy path when `--noqc` is parsef

* Fix a bad indent

* Update CHANGELOG

* Init tractography module (#234)

* Init tractography module with DSI Studio

* Update main to use tractography modules

* Update version and MANIFEST

* Update docs

Author: TheJaeger

CHANGELOG.rst

@@ -4,6 +4,35 @@ Changelog
 All notable changes to this project will be documented in this file or
 page
 
+`v1.0_RC5`_
+-----------
+
+Oct 26, 2020
+
+**Added**:
+
+* Check for b-value scaling so .bval file so values
+  specified as either 2.0 or 2000 can be processed.
+* ``fitting.dwipy()`` can now be pointed to user-defined
+  bvec and bval paths. It previously required bvec and
+  bval files to have the same name and path as DWI.
+* **DSI Studio tractography** for FBI. Processing FBI dataset now
+  produces an ``fbi_tractography_dsi.fib`` file that can be loaded
+  into DSI Studio to perform tractography.
+
+**Changed**:
+
+* Fixed issue where eddy correction would attempt
+ to QC and fail despite parsing the ``--noqc`` flag.
+* SNR plotting works in very specific scenarious when
+  input DWIs are of the same same dimensions. A try/except
+  loop now ensure that the entire pipeline doesn't halt
+  due to errors in plotting.
+
+**Removed**:
+
+* None
+
 `v1.0_RC4`_
 -----------
 
@@ -218,6 +247,7 @@ Initial port of MATLAB code to Python. 200,000,000,000 BCE
 
 .. Links
 
+.. _v1.0_RC5: https://github.com/m-ama/PyDesigner/releases/tag/v1.0_RC5
 .. _v1.0_RC4: https://github.com/m-ama/PyDesigner/releases/tag/v1.0_RC4
 .. _v1.0_RC3: https://github.com/m-ama/PyDesigner/releases/tag/v1.0_RC3
 .. _v1.0_RC2: https://github.com/m-ama/PyDesigner/releases/tag/v1.0_RC2

MANIFEST.in

@@ -1,2 +1,3 @@
 include designer/fitting/*.csv
+include desiger/tractography/*.mat
 global-exclude *.py[cod] __pycache__ *.so

designer/fitting/dwipy.py

@@ -47,7 +47,7 @@ class DWI(object):
     workers : int
         Number of CPU workers to use in processing
     """
-    def __init__(self, imPath, mask=None, nthreads=-1):
+    def __init__(self, imPath, bvecPath=None, bvalPath=None, mask=None, nthreads=-1):
         """
         DWI class initializer
 
@@ -71,20 +71,37 @@ def __init__(self, imPath, mask=None, nthreads=-1):
         (path, file) = os.path.split(imPath)
         # Remove extension from NIFTI filename
         fName = os.path.splitext(file)[0]
-        # Add .bval to NIFTI filename
-        bvalPath = os.path.join(path, fName + '.bval')
-        # Add .bvec to NIFTI filename
-        bvecPath = os.path.join(path, fName + '.bvec')
+        if bvecPath:
+            if not isinstance(bvecPath, str):
+                raise TypeError('Path to .bvec is not specified '
+                'as a string')
+            if not os.path.exists(bvecPath):
+                raise OSError('Path to .bvec does not exist: '
+                '{}'.format(bvecPath))
+        else:
+            bvecPath = os.path.join(path, fName + '.bvec')
+        if bvalPath:
+            if not isinstance(bvalPath, str):
+                raise TypeError('Path to .bval is not specified '
+                'as a string')
+            if not os.path.exists(bvalPath):
+                raise OSError('Path to .bvec does not exist: '
+                '{}'.format(bvalPath))
+        else:
+            bvalPath = os.path.join(path, fName + '.bval')
         if os.path.exists(bvalPath) and os.path.exists(bvecPath):
             # Load bvecs
             bvecs = np.loadtxt(bvecPath)
             # Load bvals
-            bvals = np.rint(np.loadtxt(bvalPath) / 1000)
+            bvals = np.rint(np.loadtxt(bvalPath))
+            # Scale bvals by checking for number of digits in max bval
+            if int(np.log10(np.max(bvals)))+1 >= 3: # if no. of digits >= 3
+                bvals = bvals / 1000
             # Combine bvecs and bvals into [n x 4] array where n is
             # number of DWI volumes. [Gx Gy Gz Bval]
             self.grad = np.c_[np.transpose(bvecs), bvals]
         else:
-            raise NameError('Unable to locate BVAL or BVEC files')
+            raise OSError('Unable to locate BVAL or BVEC files')
         if mask is None:
             maskPath = os.path.join(path,'brain_mask.nii')
         else:

designer/info.py

@@ -10,7 +10,7 @@
     )
 )
 __packagename__ = 'PyDesigner'
-__version__='v1.0-RC4'
+__version__='v1.0-RC5'
 __author__ = 'PyDesigner developers'
 __copyright__ = 'Copyright 2020, PyDesigner developers, MUSC Advanced Image Analysis (MAMA)'
 __credits__ = [
@@ -48,7 +48,8 @@
     'joblib >= 0.16',
     'tqdm >= 4.40',
     'multiprocess >= 0.70',
-    'nibabel >= 3.1',
+    'nibabel >= 3.2',
+    'dipy >= 1.2',
     'cvxpy >= 1.1'
 ]
 

designer/preprocessing/mrpreproc.py

@@ -303,8 +303,8 @@ def undistort(input, output, rpe='rpe_header', epib0=1,
     epib0 : int
         Number of reverse PE dir B0 pairs to use in TOPUP correction
         (Default: 1)
-    qc : bool
-        Specify whether to generate eddy QC metrics (Default: True)
+    qc : str
+        Specify path to QC directior. No QC metrics generated if None
     nthreads : int, optional
         Specify the number of threads to use in processing
         (Default: all available threads)

designer/pydesigner.py

@@ -20,6 +20,7 @@
 from designer.plotting import snrplot, outlierplot, motionplot
 from designer.fitting import dwipy as dp
 from designer.postprocessing import filters
+from designer.tractography import dsistudio as ds
 DWIFile = util.DWIFile
 DWIParser = util.DWIParser
 
@@ -591,6 +592,8 @@ def main():
         undistorted_name_full = str(step_count)+ '_' + undistorted_name
         nii_undistorted = op.join(intermediatepath, undistorted_name_full + '.nii')
         mif_undistorted = op.join(outpath, undistorted_name_full + '.mif')
+        if args.noqc:
+            eddyqcpath = None
         # check to see if this already exists
         if not (args.resume and op.exists(nii_undistorted)):
             # run undistort function
@@ -764,15 +767,19 @@ def main():
         files = []
         files.append(init_nii)
         files.append(filetable['HEAD'].getFull())
-        if 'mask' in filetable:
-            snr = snrplot.makesnr(dwilist=files,
-                                  noisepath=nii_noisemap,
-                                  maskpath=filetable['mask'].getFull())
-        else:
-            snr = snrplot.makesnr(dwilist=files,
-                                  noisepath=filetable['noisemap'].getFull(),
-                                  maskpath=None)
-        snr.makeplot(path=qcpath, smooth=True, smoothfactor=3)
+        try:
+            if 'mask' in filetable:
+                snr = snrplot.makesnr(dwilist=files,
+                                    noisepath=nii_noisemap,
+                                    maskpath=filetable['mask'].getFull())
+            else:
+                snr = snrplot.makesnr(dwilist=files,
+                                    noisepath=filetable['noisemap'].getFull(),
+                                    maskpath=None)
+            snr.makeplot(path=qcpath, smooth=True, smoothfactor=3)
+        except:
+            print('[WARNING] SNR plotting failed, see above. '
+            'Proceeding with processing.')
 
     #-----------------------------------------------------------------
     # Write logs
@@ -823,6 +830,7 @@ def main():
         fn_fbi_zeta = 'fbi_zeta'
         fn_fbi_faa = 'fbi_faa'
         fn_fbi_sph = 'fbi_fodf'
+        fn_fbi_tract = 'fbi_tractography_dsi'
         fn_fbi_awf = 'fbwm_awf'
         fn_fbi_Da = 'fbwm_da'
         fn_fbi_De_mean = 'fbwm_de_mean'
@@ -875,9 +883,10 @@ def main():
             if (img.isdti() or img.isdki()) and not args.noakc:
                 akc_out = img.akcoutliers()
                 img.akccorrect(akc_out)
-                dp.writeNii(akc_out,
-                            img.hdr,
-                            op.join(fitqcpath, 'outliers_akc'))
+                if not args.noqc:
+                    dp.writeNii(akc_out,
+                                img.hdr,
+                                op.join(fitqcpath, 'outliers_akc'))
 
              # reorder tensor for mrtrix3
             if 'dki' in img.tensorType():
@@ -1015,6 +1024,17 @@ def main():
                             input=op.join(metricpath, x + fn_ext),
                             output=op.join(metricpath, x + fn_ext),
                             mask=filetable['mask'].getFull())
-
+            if 'mask' in filetable:
+                ds.makefib(
+                    input=op.join(metricpath, fn_fbi_sph + fn_ext),
+                    output=op.join(metricpath, fn_fbi_tract + '.fib'),
+                    mask=filetable['mask'].getFull()
+                )
+            else:
+                ds.makefib(
+                    input=op.join(metricpath, fn_fbi_sph + fn_ext),
+                    output=op.join(metricpath, fn_fbi_tract + '.fib'),
+                    mask=None
+                )
 if __name__ == '__main__':
     main()