Merge pull request #144 from PennLINC/sdcflow

Sdcflow

.circleci/config.yml

@@ -175,7 +175,7 @@ jobs:
                 /out  participant --participant-label 01 \
                 -w /scrth/wkdir_sub01 --output-spaces asl \
                 --fs-license-file /scrth/licence.txt \
-                --scorescrub --basil  \
+                --scorescrub --basil --use-syn-sdc \
                 --anat-derivatives /data/smriprep
 
       - run:
@@ -264,7 +264,7 @@ jobs:
                 /out  participant --participant-label A00086748 \
                 -w /scrth/wkdir_subA00086748 --output-spaces asl \
                 --fs-license-file /scrth/licence.txt \
-                --scorescrub --basil  \
+                --scorescrub --basil --use-syn-sdc \
                 --anat-derivatives /data/smriprep
 
       - run:

aslprep/interfaces/cbf_computation.py

@@ -319,7 +319,7 @@ def cbfcomputation(metadata, mask, m0file, cbffile, m0scale=1):
         for i in range(cbf_data.shape[1]):
             cbf_data_ts[:, i] =np.multiply(cbf1[:, i],permfactor[i])
         cbf = np.zeros([cbf_data_ts.shape[0], int(cbf_data.shape[1]/len(perfusion_factor))])
-        cbf_xx=np.split(cbf_data_ts,int(cbf_data_ts.shape[1]/len(perfusion_factor)),axis=1)
+        cbf_xx = np.split(cbf_data_ts,int(cbf_data_ts.shape[1]/len(perfusion_factor)),axis=1)
 
         # calculate weighted cbf with multiplds
         # https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3791289/
@@ -344,9 +344,9 @@ def cbfcomputation(metadata, mask, m0file, cbffile, m0scale=1):
         tcbf=np.zeros(maskx.shape)
         tcbf[maskx==1]=cbf
     else:
-        tcbf=np.zeros([maskx.shape[0],maskx.shape[1],maskx.shape[2],cbf.shape[1]])
+        tcbf = np.zeros([maskx.shape[0],maskx.shape[1],maskx.shape[2],cbf.shape[1]])
         for i in range(cbf.shape[1]):
-            tcbfx=np.zeros(maskx.shape) 
+            tcbfx = np.zeros(maskx.shape) 
             tcbfx[maskx==1]=cbf[:,i]
             tcbf[:,:,:,i]=tcbfx
     if len(tcbf.shape) < 4:
@@ -458,7 +458,7 @@ def _weightfun(x, wfun='huber'):
         weight = 1/(1+np.power(x, 2))
     elif wfun == 'logistic':
         tuner = 1.205
-        weight == np.tanh(x)/x
+        weight = np.tanh(x)/x
     elif wfun == 'ols':
         tuner = 1
         weight = np.repeat(1, len(x))

aslprep/niworkflows/viz/plots.py

@@ -241,9 +241,9 @@ def plotstatsimg(cbf, ref_vol, plot_params=None, order=('z', 'x', 'y'),
         plot_params['display_mode'] = mode
         plot_params['cut_coords'] = cuts[mode]
         plot_params['draw_cross'] = False
-        plot_params['symmetric_cbar'] = False
-        plot_params['threshold'] = 0.05
+        plot_params['symmetric_cbar'] = True
         plot_params['vmax'] = 90
+        plot_params['threshold'] = 0.02
         plot_params['colorbar'] = False
         plot_params['cmap'] = 'gray'
         if i == 0:
@@ -254,7 +254,8 @@ def plotstatsimg(cbf, ref_vol, plot_params=None, order=('z', 'x', 'y'),
 
         # Generate nilearn figure
         from nilearn.plotting import plot_stat_map
-        display = plot_stat_map(stat_map_img=cbf, bg_img=ref_vol, **plot_params)
+        display = plot_stat_map(stat_map_img=cbf, bg_img=ref_vol,
+                     resampling_interpolation='nearest',**plot_params)
         svg = extract_svg(display, compress=compress)
         display.close()
         from lxml import etree

aslprep/sdcflows/__init__.py

@@ -0,0 +1,13 @@
+"""SDCflows - :abbr:`SDC (susceptibility distortion correction)` by DUMMIES, for dummies."""
+__packagename__ = "sdcflows"
+__copyright__ = "2020, The NiPreps developers"
+try:
+    from ._version import __version__
+except ModuleNotFoundError:
+    from pkg_resources import get_distribution, DistributionNotFound
+    try:
+        __version__ = get_distribution(__packagename__).version
+    except DistributionNotFound:
+        __version__ = "unknown"
+    del get_distribution
+    del DistributionNotFound

aslprep/sdcflows/cli/run.py

@@ -0,0 +1,113 @@
+#!/usr/bin/env python3
+import logging
+from pathlib import Path
+
+logging.addLevelName(25, 'IMPORTANT')  # Add a new level between INFO and WARNING
+logging.addLevelName(15, 'VERBOSE')  # Add a new level between INFO and DEBUG
+logger = logging.getLogger('sdcflows')
+
+
+def get_parser():
+    """Define the command line interface"""
+    from argparse import ArgumentParser
+    from argparse import RawTextHelpFormatter
+    from .. import __version__ as _vstr
+
+    parser = ArgumentParser(description='SDC Workflows',
+                            formatter_class=RawTextHelpFormatter)
+
+    parser.add_argument(
+        'bids_dir', action='store', type=Path,
+        help='the root folder of a BIDS dataset')
+    parser.add_argument('output_dir', action='store', type=Path,
+                        help='the output path for the outcomes of preprocessing and visual '
+                             'reports')
+    parser.add_argument('analysis_level', choices=['participant', 'group'], nargs='+',
+                        help='processing stage to be run, "participant" means individual analysis '
+                             'and "group" is second level analysis.')
+    # optional arguments
+    parser.add_argument('--version', action='version', version='v{}'.format(_vstr))
+
+    # Options that affect how pyBIDS is configured
+    g_bids = parser.add_argument_group('Options for filtering BIDS queries')
+    g_bids.add_argument('--participant-label', action='store', type=str,
+                        nargs='*', dest='subject', help='process only particular subjects')
+    g_bids.add_argument('--task', action='store', type=str, nargs='*',
+                        help='select a specific task to be processed')
+    g_bids.add_argument('--dir', action='store', type=str, nargs='*',
+                        help='select a specific direction entity to be processed')
+    g_bids.add_argument('--acq', action='store', type=str, nargs='*', dest='acquisition',
+                        help='select a specific acquisition entity to be processed')
+    g_bids.add_argument('--run', action='store', type=int, nargs='*',
+                        help='select a specific run identifier to be processed')
+    g_bids.add_argument('--suffix', action='store', type=str, nargs='*', default='bold',
+                        help='select a specific run identifier to be processed')
+
+    g_perfm = parser.add_argument_group('Options to handle performance')
+    g_perfm.add_argument("-v", "--verbose", dest="verbose_count", action="count", default=0,
+                         help="increases log verbosity for each occurence, debug level is -vvv")
+    g_perfm.add_argument('--ncpus', '--nprocs', action='store', type=int,
+                         help='maximum number of threads across all processes')
+    g_perfm.add_argument('--nthreads', '--omp-nthreads', action='store', type=int,
+                         help='maximum number of threads per-process')
+
+    g_other = parser.add_argument_group('Other options')
+    g_other.add_argument('-w', '--work-dir', action='store', type=Path,
+                         help='path where intermediate results should be stored')
+
+    return parser
+
+
+def main():
+    """Entry point"""
+    from os import cpu_count
+    from multiprocessing import set_start_method
+    # from bids.layout import BIDSLayout
+    from nipype import logging as nlogging
+    set_start_method('forkserver')
+
+    opts = get_parser().parse_args()
+
+    # Retrieve logging level
+    log_level = int(max(25 - 5 * opts.verbose_count, logging.DEBUG))
+    # Set logging
+    logger.setLevel(log_level)
+    nlogging.getLogger('nipype.workflow').setLevel(log_level)
+    nlogging.getLogger('nipype.interface').setLevel(log_level)
+    nlogging.getLogger('nipype.utils').setLevel(log_level)
+
+    # Resource management options
+    plugin_settings = {
+        'plugin': 'MultiProc',
+        'plugin_args': {
+            'n_procs': opts.ncpus,
+            ''
+            'raise_insufficient': False,
+            'maxtasksperchild': 1,
+        }
+    }
+    # Permit overriding plugin config with specific CLI options
+    if not opts.ncpus or opts.ncpus < 1:
+        plugin_settings['plugin_args']['n_procs'] = cpu_count()
+
+    nthreads = opts.nthreads
+    if not nthreads or nthreads < 1:
+        nthreads = cpu_count()
+
+    # output_dir = opts.output_dir.resolve()
+    # bids_dir = opts.bids_dir or output_dir.parent
+
+    # Get absolute path to BIDS directory
+    # bids_dir = opts.bids_dir.resolve()
+    # layout = BIDSLayout(str(bids_dir), validate=False, derivatives=str(output_dir))
+    # query = {'suffix': opts.suffix, 'extension': ['.nii', '.nii.gz']}
+
+    # for entity in ('subject', 'task', 'dir', 'acquisition', 'run'):
+    #     arg = getattr(opts, entity, None)
+    #     if arg is not None:
+    #         query[entity] = arg
+
+
+if __name__ == '__main__':
+    raise RuntimeError("sdcflows/cli/run.py should not be run directly;\n"
+                       "Please `pip install` sdcflows and use the `sdcflows` command")

aslprep/sdcflows/conftest.py

@@ -0,0 +1,44 @@
+"""py.test configuration"""
+import os
+from pathlib import Path
+import numpy
+import pytest
+from bids.layout import BIDSLayout
+
+test_data_env = os.getenv('TEST_DATA_HOME', str(Path.home() / 'sdcflows-tests'))
+test_output_dir = os.getenv('TEST_OUTPUT_DIR')
+test_workdir = os.getenv('TEST_WORK_DIR')
+
+layouts = {p.name: BIDSLayout(str(p), validate=False, derivatives=True)
+           for p in Path(test_data_env).glob('*') if p.is_dir()}
+
+
+def pytest_report_header(config):
+    msg = "Datasets found: %s" % ', '.join([v.root for v in layouts.values()])
+    if test_output_dir is not None:
+        msg += '\nOutput folder: %s' % Path(test_output_dir).resolve()
+    return msg
+
+
+@pytest.fixture(autouse=True)
+def add_np(doctest_namespace):
+    doctest_namespace['np'] = numpy
+    doctest_namespace['os'] = os
+    doctest_namespace['Path'] = Path
+    for key, val in list(layouts.items()):
+        doctest_namespace[key] = Path(val.root)
+
+
+@pytest.fixture
+def workdir():
+    return None if test_workdir is None else Path(test_workdir)
+
+
+@pytest.fixture
+def output_path():
+    return None if test_output_dir is None else Path(test_output_dir)
+
+
+@pytest.fixture
+def bids_layouts():
+    return layouts

aslprep/sdcflows/data/affine.json

@@ -0,0 +1,24 @@
+{
+    "dimension": 3,
+    "float": false,
+    "initial_moving_transform_com": 1,
+    "winsorize_lower_quantile": 0.005,
+    "winsorize_upper_quantile": 0.995,
+    "collapse_output_transforms": true,
+    "use_histogram_matching": [ false, false ],
+    "use_estimate_learning_rate_once": [ false, false ],
+    "transforms": [ "Rigid", "Affine" ],
+    "number_of_iterations": [ [ 1000, 500, 250, 100 ], [ 1000, 500, 250, 100 ] ],
+    "transform_parameters": [ [ 0.1 ], [ 0.1 ] ],
+    "convergence_threshold": [ 1e-06, 1e-06 ],
+    "convergence_window_size": [ 10, 10 ],
+    "metric": [ "MI", "MI" ],
+    "sampling_percentage": [ 0.25, 0.25 ],
+    "sampling_strategy": [ "Regular", "Regular" ],
+    "smoothing_sigmas": [ [ 3, 2, 1, 0 ], [ 3, 2, 1, 0 ] ],
+    "sigma_units": [ "vox", "vox" ],
+    "metric_weight": [ 1.0, 1.0 ],
+    "shrink_factors": [[ 8, 4, 2, 1] , [ 8, 4, 2, 1 ]],
+    "radius_or_number_of_bins": [ 32, 32 ],
+    "interpolation": "Linear"
+}

aslprep/sdcflows/data/fmap-any_registration.json

@@ -0,0 +1,24 @@
+{
+    "dimension": 3,
+    "float": true,
+    "winsorize_lower_quantile": 0.005,
+    "winsorize_upper_quantile": 0.998,
+    "collapse_output_transforms": true,
+    "write_composite_transform": true,
+    "use_histogram_matching": [ true, true ],
+    "use_estimate_learning_rate_once": [ true, true ],
+    "transforms": [ "Translation", "Affine" ],
+    "number_of_iterations": [ [ 500 ], [ 200 ] ],
+    "transform_parameters": [ [ 0.05 ], [ 0.01 ] ],
+    "convergence_threshold": [ 1e-07, 1e-08 ],
+    "convergence_window_size": [ 200, 100 ],
+    "metric": [ "Mattes", "Mattes" ],
+    "sampling_percentage": [ 0.5, 0.5 ],
+    "sampling_strategy": [ "Random", "Random" ],
+    "smoothing_sigmas": [ [ 8.0 ], [ 2.0 ] ],
+    "sigma_units": [ "mm", "mm" ],
+    "metric_weight": [ 1.0, 1.0 ],
+    "shrink_factors": [ [ 2 ], [ 1 ] ],
+    "radius_or_number_of_bins": [ 64, 64 ],
+    "interpolation": "LanczosWindowedSinc"
+}

aslprep/sdcflows/data/fmap-any_registration_testing.json

@@ -0,0 +1,24 @@
+{
+    "dimension": 3,
+    "float": true,
+    "winsorize_lower_quantile": 0.005,
+    "winsorize_upper_quantile": 0.998,
+    "collapse_output_transforms": true,
+    "write_composite_transform": true,
+    "use_histogram_matching": [ true, true ],
+    "use_estimate_learning_rate_once": [ true, true ],
+    "transforms": [ "Rigid", "Rigid" ],
+    "number_of_iterations": [ [ 500 ], [ 100 ] ],
+    "transform_parameters": [ [ 1.0 ], [ 0.5 ] ],
+    "convergence_threshold": [ 1e-07, 1e-08 ],
+    "convergence_window_size": [ 20, 10 ],
+    "metric": [ "Mattes", "Mattes" ],
+    "sampling_percentage": [ 1.0, 1.0 ],
+    "sampling_strategy": [ "Regular", "Regular" ],
+    "smoothing_sigmas": [ [ 8.0 ], [ 2.0 ] ],
+    "sigma_units": [ "mm", "mm" ],
+    "metric_weight": [ 1.0, 1.0 ],
+    "shrink_factors": [ [ 2 ], [ 1 ] ],
+    "radius_or_number_of_bins": [ 64, 64 ],
+    "interpolation": "LanczosWindowedSinc"
+}

aslprep/sdcflows/data/susceptibility_syn.json

@@ -0,0 +1,24 @@
+{
+    "dimension": 3,
+    "float": true,
+    "convergence_threshold": [ 1e-08, 1e-08 ],
+    "convergence_window_size": [ 20, 10 ],
+    "metric": [ "Mattes", "CC" ],
+    "metric_weight": [ 1, 1 ],
+    "radius_or_number_of_bins": [ 56, 5 ],
+    "transforms": [ "SyN", "SyN" ],
+    "transform_parameters": [ [ 0.8, 2, 2 ], [ 0.8, 2, 2 ] ],
+    "number_of_iterations": [ [ 100, 50 ], [ 20, 10 ] ],
+    "smoothing_sigmas": [ [ 1, 0 ], [ 1, 0 ] ],
+    "sigma_units": [ "vox", "vox" ],
+    "shrink_factors": [ [ 2, 1 ], [ 1, 1 ] ],
+    "winsorize_upper_quantile": 1.0,
+    "winsorize_lower_quantile": 0.001,
+    "use_estimate_learning_rate_once": [ true, true ],
+    "use_histogram_matching": [ true, true ],
+    "collapse_output_transforms": true,
+    "write_composite_transform": false,
+    "output_transform_prefix": "ants_susceptibility",
+    "output_warped_image": true,
+    "interpolation": "Linear"
+}

aslprep/sdcflows/data/translation_rigid.json

@@ -0,0 +1,24 @@
+{
+    "dimension": 3,
+    "float": true,
+    "winsorize_lower_quantile": 0.005,
+    "winsorize_upper_quantile": 0.998,
+    "collapse_output_transforms": true,
+    "write_composite_transform": true,
+    "use_histogram_matching": [ true, true ],
+    "use_estimate_learning_rate_once": [ true, true ],
+    "transforms": [ "Translation", "Rigid" ],
+    "number_of_iterations": [ [ 500 ], [ 200 ] ],
+    "transform_parameters": [ [ 0.05 ], [ 0.01 ] ],
+    "convergence_threshold": [ 1e-07, 1e-08 ],
+    "convergence_window_size": [ 200, 100 ],
+    "metric": [ "Mattes", "Mattes" ],
+    "sampling_percentage": [ 0.5, 0.5 ],
+    "sampling_strategy": [ "Random", "Random" ],
+    "smoothing_sigmas": [ [ 8.0 ], [ 2.0 ] ],
+    "sigma_units": [ "mm", "mm" ],
+    "metric_weight": [ 1.0, 1.0 ],
+    "shrink_factors": [ [ 2 ], [ 1 ] ],
+    "radius_or_number_of_bins": [ 64, 64 ],
+    "interpolation": "LanczosWindowedSinc"
+}