Update shell-sampling determination (#241)

* 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
muscbridge · Dec 22, 2020 · d700939 · d700939
1 parent 48b7af5
commit d700939
Show file tree

Hide file tree

Showing 5 changed files with 177 additions and 83 deletions.
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