Update shell-sampling determination (#241)

TheJaeger · web-flow · commit d70093948213 · 2020-12-22T16:44:08.000-05:00
* 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 * Fiber Tracking: FBI tractography (#235) (#237) * 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 * Change dki_ak fname (#238) * Update shell sampling determination (#240) * Add new functions to determine DWI shell sampling * Remove references to old functions * Change how shell sampling is determined * Update documentation
diff --git a/CHANGELOG.rst b/CHANGELOG.rst
@@ -4,6 +4,31 @@ Changelog
 All notable changes to this project will be documented in this file or
 page
 
+`v1.0_RC6`_
+-----------
+
+Dec 22, 2020
+
+**Added**:
+
+* None
+
+**Changed**
+
+* Replaced ``preprocessing.util.bvec_is_fullsphere()`` and 
+  ``preprocessing.util.vecs_are_fullsphere()`` with 
+  ``preprocessing.mrinfoutil.is_fullsphere()``. Even though datasets
+  may be half-shelled, it is inaccurate to label them as such because
+  distortion relative to b-value is not linear. As such, the
+  ``slm=linear`` makes no sense. This new method performs the proper
+  checks required before labelling a DWI as fully-shelled. A DWI is
+  half-shelled iff max B-value is less than 3000 AND the norm of the
+  mean direction vector is more than 0.3.
+
+**Removed**
+
+* See above
+
 `v1.0_RC5`_
 -----------
 
@@ -247,6 +272,7 @@ Initial port of MATLAB code to Python. 200,000,000,000 BCE
 
 .. Links
 
+.. _v1.0_RC6: https://github.com/m-ama/PyDesigner/releases/tag/v1.0_RC6
 .. _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
diff --git a/designer/info.py b/designer/info.py
@@ -10,7 +10,7 @@
     )
 )
 __packagename__ = 'PyDesigner'
-__version__='v1.0-RC5'
+__version__='v1.0-RC6'
 __author__ = 'PyDesigner developers'
 __copyright__ = 'Copyright 2020, PyDesigner developers, MUSC Advanced Image Analysis (MAMA)'
 __credits__ = [
diff --git a/designer/preprocessing/mrinfoutil.py b/designer/preprocessing/mrinfoutil.py
@@ -353,3 +353,152 @@ def pescheme(path):
             nums.append(float(num))
         pe_scheme.append(nums)
     return pe_scheme
+
+def num_shells(path):
+    """
+    Returns the number of b-value shells detected in input file
+
+    Parameters
+    ----------
+    path : str
+        Path to input image or directory
+
+    Returns
+    -------
+    int
+        Number of shells
+    """
+    if not op.exists(path):
+        raise OSError('Input path does not exist. Please ensure that the '
+                      'folder or file specified exists.')
+    ftype = format(path)
+    if ftype != 'MRtrix':
+        raise IOError('This function only works with MRtrix (.mif) '
+                      'formatted filetypes. Please ensure that the input '
+                      'filetype meets this requirement')
+    arg = ['mrinfo', '-shell_bvalues']
+    arg.append(path)
+    completion = subprocess.run(arg, stdout=subprocess.PIPE)
+    if completion.returncode != 0:
+        raise IOError('Input {} is not currently supported by '
+                      'PyDesigner.'.format(path))
+    # Remove new line delimiter
+    console = str(completion.stdout).split('\\n')
+    # Remove 'b'
+    console[0] = console[0][1:]
+    # Remove quotes
+    console = [s.replace("'", "") for s in console]
+    # Condense empty strings
+    console = [s.replace('"', '') for s in console]
+    # Remove empty strings form list
+    console.remove('')
+    # Split spaces
+    console = [s.split(' ') for s in console]
+    console = [item for sublist in console for item in sublist]
+    console = list(filter(None, console))
+    return len(console)
+
+def max_shell(path):
+    """
+    Returns the maximum b-value shell in DWI
+
+    Parameters
+    ----------
+    path : str
+        Path to input image or directory
+
+    Returns
+    -------
+    int
+        Max b-value
+    """
+    if not op.exists(path):
+        raise OSError('Input path does not exist. Please ensure that the '
+                      'folder or file specified exists.')
+    ftype = format(path)
+    if ftype != 'MRtrix':
+        raise IOError('This function only works with MRtrix (.mif) '
+                      'formatted filetypes. Please ensure that the input '
+                      'filetype meets this requirement')
+    arg = ['mrinfo', '-shell_bvalues']
+    arg.append(path)
+    completion = subprocess.run(arg, stdout=subprocess.PIPE)
+    if completion.returncode != 0:
+        raise IOError('Input {} is not currently supported by '
+                      'PyDesigner.'.format(path))
+    # Remove new line delimiter
+    console = str(completion.stdout).split('\\n')
+    # Remove 'b'
+    console[0] = console[0][1:]
+    # Remove quotes
+    console = [s.replace("'", "") for s in console]
+    # Condense empty strings
+    console = [s.replace('"', '') for s in console]
+    # Remove empty strings form list
+    console.remove('')
+    # Split spaces
+    console = [s.split(' ') for s in console]
+    console = [item for sublist in console for item in sublist]
+    console = list(filter(None, console))
+    console = [round(int(s)) for s in console]
+    print(console)
+    return max(console)
+
+def is_fullsphere(path):
+    """
+    Returns boolean value indicating whether input file has full
+    spherical sampling
+
+    Parameters
+    ----------
+    path : str
+        Path to input image or directory
+
+    Returns
+    -------
+    bool
+        True if full spherical sampling
+        False if half-spherical sampling
+    """
+    if not op.exists(path):
+        raise OSError('Input path does not exist. Please ensure that the '
+                        'folder or file specified exists.')
+    ftype = format(path)
+    if ftype != 'MRtrix':
+        raise IOError('This function only works with MRtrix (.mif) '
+                        'formatted filetypes. Please ensure that the input '
+                        'filetype meets this requirement')
+    arg = ['dirstat', path, '-output', 'ASYM']
+    completion = subprocess.run(arg, stdout=subprocess.PIPE)
+    if completion.returncode != 0:
+        raise IOError('Input {} is not currently supported by '
+                        'PyDesigner.'.format(path))
+    # Remove new line delimiter
+    console = str(completion.stdout).split('\\n')
+    # Remove 'b'
+    console[0] = console[0][1:]
+    # Remove quotes
+    console = [s.replace("'", "") for s in console]
+    # Condense empty strings
+    console = [s.replace('"', '') for s in console]
+    # Remove empty strings form list
+    console.remove('')
+    # Convert strings to list
+    console = [float(s) for s in console]
+    # Find mean of all b-shells and round to nearest one decimal place
+    mean_dir = round(sum(console) / len(console), 1)
+    # Having too many b-value shells for protocols such as FBI or
+    # HARDI, eddy correction can create a lot of trouble. If that's
+    # the case, we just assume that data is fully shelled. This
+    # disables heuristics testing performed by eddy correction. A DWI
+    # will return half-shelled if and only if the norm of mean
+    # direction vector is greater that 0.3 AND has max b-value less
+    # than or equal to 3000. The linear model breaks down for high
+    # b-values.
+    if mean_dir < 0.3:
+        return True
+    else:
+        if max_shell(path) > 3000:
+            return True
+        else:
+            return False
diff --git a/designer/preprocessing/mrpreproc.py b/designer/preprocessing/mrpreproc.py
@@ -376,7 +376,7 @@ def undistort(input, output, rpe='rpe_header', epib0=1,
         arg.extend(['-nthreads', str(nthreads)])
     # Determine whether half or full sphere sampling
     repol_string = '--repol '
-    if util.bvec_is_fullsphere(op.join(outdir, 'dwiec.bvec')):
+    if mrinfoutil.is_fullsphere(input):
         # is full, add appropriate dwifslpreproc option
         repol_string += '--data_is_shelled'
     else:
diff --git a/designer/preprocessing/util.py b/designer/preprocessing/util.py
@@ -13,87 +13,6 @@
 import warnings
 from designer.preprocessing import mrinfoutil
 
-def bvec_is_fullsphere(bvec):
-    """
-    Determines if .bvec file is full or half-sphere
-
-    Parameters
-    ----------
-    bvec : str
-        The filename of the bvec
-
-    Returns
-    -------
-    bool
-        True if full-sphere, False if half
-    """
-
-    # Check for existence
-    if not op.exists(bvec):
-        raise Exception('.bvec file '+bvec+' does not exist.')
-    # Attempt to load file
-    try:
-        data = np.loadtxt(bvec)
-    except:
-        raise Exception('.bvec file '+bvec+' cannot be read.')
-
-    # Transpose data so 2nd dimension is the axis i,j,k (x,y,z)
-    data = np.transpose(data)
-    # Check that data has correct dimensions
-    size_dim0 = np.size(data, 0)
-    if np.size(data, 1) != 3:
-        raise Exception('.bvec file '+bvec+' is not 3D; is '+str(size_dim0))
-
-    return vecs_are_fullsphere(data)
-
-def vecs_are_fullsphere(bvecs):
-    """
-    Determines if input vectors are full or half-sphere
-
-    Parameters
-    ----------
-    bvecs : ndarray
-        Aray of size [n_vectors x 3]
-
-    Returns
-    -------
-    bool
-        True if full-sphere, False if half
-
-    Notes
-    -----
-    Adapted from Chris Rorden's nii_preprocess as seen here:
-    https://github.com/neurolabusc/nii_preprocess/blob/dd1c84f23f8828923dd5fc493a22156b7006a3d4/nii_preprocess.m#L1786-L1824
-    and reproduced from the original designer pipeline here:
-    https://github.com/m-ama/PyDesigner/blob/7a39ec4cb9679f1c3e9ead68baa8f8c111b0618a/designer/designer.py#L347-L368
-    """
-    # TODO: figure out why this math works
-    # Check dimensions
-    if np.size(bvecs, 1) != 3:
-        raise Exception('bvecs are not 3-dimensional.')
-    # Remove NaNs
-    bvecs[~np.isnan(bvecs.any(axis=1))]
-    # Remove any 0-vectors
-    bvecs[~np.all(bvecs == 0, axis=1)]
-    # Assume half-sphere if no vectors remaining
-    if not bvecs.any():
-        raise Exception('bvecs do not point anywhere.')
-    # Get mean length
-    mean = np.mean(bvecs, 0)
-    # Get x-component of direction
-    x_component = np.sqrt(np.sum(np.power(mean, 2)))
-    # Create a matrix to divide by this length
-    Dx = np.repeat(x_component, 3)
-    # Get mean unit length
-    mean_ulength = np.divide(mean, Dx)
-    # Scale by mean unit length
-    mean = np.ones([np.size(bvecs, 0), 3])* mean_ulength
-    # UNKNOWN
-    minV = min(np.sum(bvecs.conj() * mean, axis=1))
-    # Get angle in degrees
-    theta = math.degrees(math.acos(minV))
-    return (theta >= 110)
-
 def find_valid_ext(pathname):
     """
     Finds valid extensions for dwifile, helper function

Original file line number	Diff line number	Diff line change
`@@ -10,7 +10,7 @@`
`10`	`10`	`)`
`11`	`11`	`)`
`12`	`12`	`__packagename__ = 'PyDesigner'`
`13`		`-__version__='v1.0-RC5'`
	`13`	`+__version__='v1.0-RC6'`
`14`	`14`	`__author__ = 'PyDesigner developers'`
`15`	`15`	`__copyright__ = 'Copyright 2020, PyDesigner developers, MUSC Advanced Image Analysis (MAMA)'`
`16`	`16`	`__credits__ = [`