Merge branch 'master' into Fix_concatenate_3D_with_4D_images

Jenkinsfile

@@ -0,0 +1,58 @@
+pipeline {
+    agent any
+
+    stages {
+        stage('Build') {
+            stages {
+                stage('Python3.6') {
+                    steps {
+                        withPythonEnv('CPython-3.6') {
+                            sh '''
+                                pip3 install numpy==1.18.* wheel
+                                pip3 install -e .
+                            '''
+                        }
+                    }
+                }
+                stage('Python3.7') {
+                    steps {
+                        withPythonEnv('CPython-3.7') {
+                            sh '''
+                                pip3 install numpy==1.18.* wheel
+                                pip3 install -e .
+                            '''
+                        }
+                    }
+                }
+            }
+        }
+
+        stage('Test') {
+            steps {
+                withPythonEnv('CPython-3.7') {
+                    sh '''
+                        pip3 install numpy==1.18.* wheel
+                        pip3 install -e .
+                        export MPLBACKEND="agg"
+                        export OPENBLAS_NUM_THREADS=1
+                        pytest -v
+                    '''
+                }
+            }
+        }
+
+        stage('Deploy') {
+            when {
+                branch 'master'
+            }
+            steps {
+                echo 'Deploying.'
+            }
+        }
+    }
+    post {
+        always {
+            cleanWs()
+        }
+    }
+}

requirements.txt

@@ -2,7 +2,7 @@ bz2file==0.98.*
 coloredlogs==10.0.*
 cycler==0.10.*
 Cython==0.29.*
-dipy==1.2.*
+dipy==1.3.*
 fury==0.6.*
 future==0.17.*
 h5py==2.10.*
@@ -30,5 +30,6 @@ googledrivedownloader==0.*
 requests==2.23.*
 bctpy==0.5.*
 statsmodels==0.11.*
-dmri-commit==1.3.*
+dmri-commit==1.4.*
 openpyxl==2.6.*
+cvxpy==1.0.*

scilpy/denoise/asym_enhancement.py

@@ -0,0 +1,126 @@
+# -*- coding: utf-8 -*-
+
+import numpy as np
+from dipy.reconst.shm import sh_to_sf_matrix
+from dipy.data import get_sphere
+from dipy.core.sphere import Sphere
+from scipy.ndimage import correlate
+
+
+def local_asym_gaussian_filtering(in_sh, sh_order=8, sh_basis='descoteaux07',
+                                  out_full_basis=True, dot_sharpness=1.0,
+                                  sphere_str='repulsion724', sigma=1.0):
+    """Average the SH projected on a sphere using a first-neighbor gaussian
+    blur and a dot product weight between sphere directions and the direction
+    to neighborhood voxels, forcing to 0 negative values and thus performing
+    asymmetric hemisphere-aware filtering.
+
+    Parameters
+    ----------
+    in_sh: ndarray (x, y, z, n_coeffs)
+        Input SH coefficients array
+    sh_order: int, optional
+        Maximum order of the SH series.
+    sh_basis: {'descoteaux07', 'tournier07'}, optional
+        SH basis of the input signal.
+    out_full_basis: bool, optional
+        If True, save output SH using full SH basis.
+    dot_sharpness: float, optional
+        Exponent of the dot product. When set to 0.0, directions
+        are not weighted by the dot product.
+    sphere_str: str, optional
+        Name of the sphere used to project SH coefficients to SF.
+    sigma: float, optional
+        Sigma for the Gaussian.
+
+    Returns
+    -------
+    out_sh: ndarray (x, y, z, n_coeffs)
+        Filtered signal as SH coefficients.
+    """
+    # Load the sphere used for projection of SH
+    sphere = get_sphere(sphere_str)
+
+    # Normalized filter for each sf direction
+    weights = _get_weights(sphere, dot_sharpness, sigma)
+
+    # Detect if the basis is full based on its order
+    # and the number of coefficients of the SH
+    in_full_basis = in_sh.shape[-1] == (sh_order + 1)**2
+
+    nb_sf = len(sphere.vertices)
+    mean_sf = np.zeros(np.append(in_sh.shape[:-1], nb_sf))
+    B = sh_to_sf_matrix(sphere, sh_order=sh_order, basis_type=sh_basis,
+                        return_inv=False, full_basis=in_full_basis)
+
+    # We want a B matrix to project on an inverse sphere to have the sf on
+    # the opposite hemisphere for a given vertice
+    neg_B = sh_to_sf_matrix(Sphere(xyz=-sphere.vertices), sh_order=sh_order,
+                            basis_type=sh_basis, return_inv=False,
+                            full_basis=in_full_basis)
+
+    # Apply filter to each sphere vertice
+    for sf_i in range(nb_sf):
+        w_filter = weights[..., sf_i]
+
+        # Calculate contribution of center voxel
+        current_sf = np.dot(in_sh, B[:, sf_i])
+        mean_sf[..., sf_i] = w_filter[1, 1, 1] * current_sf
+
+        # Add contributions of neighbors using opposite hemispheres
+        current_sf = np.dot(in_sh, neg_B[:, sf_i])
+        w_filter[1, 1, 1] = 0.0
+        mean_sf[..., sf_i] += correlate(current_sf, w_filter, mode="constant")
+
+    # Convert back to SH coefficients
+    _, B_inv = sh_to_sf_matrix(sphere, sh_order=sh_order, basis_type=sh_basis,
+                               full_basis=out_full_basis)
+
+    out_sh = np.array([np.dot(i, B_inv) for i in mean_sf], dtype=in_sh.dtype)
+    return out_sh
+
+
+def _get_weights(sphere, dot_sharpness, sigma):
+    """
+    Get neighbors weight in respect to the direction to a voxel.
+
+    Parameters
+    ----------
+    sphere: Sphere
+        Sphere used for SF reconstruction.
+    dot_sharpness: float
+        Dot product exponent.
+    sigma: float
+        Variance of the gaussian used for weighting neighbors.
+
+    Returns
+    -------
+    weights: dictionary
+        Vertices weights with respect to voxel directions.
+    """
+    directions = np.zeros((3, 3, 3, 3))
+    for x in range(3):
+        for y in range(3):
+            for z in range(3):
+                directions[x, y, z, 0] = x - 1
+                directions[x, y, z, 1] = y - 1
+                directions[x, y, z, 2] = z - 1
+
+    non_zero_dir = np.ones((3, 3, 3), dtype=bool)
+    non_zero_dir[1, 1, 1] = False
+
+    # normalize dir
+    dir_norm = np.linalg.norm(directions, axis=-1, keepdims=True)
+    directions[non_zero_dir] /= dir_norm[non_zero_dir]
+
+    g_weights = np.exp(-dir_norm**2 / (2 * sigma**2))
+    d_weights = np.dot(directions, sphere.vertices.T)
+
+    d_weights = np.where(d_weights > 0.0, d_weights**dot_sharpness, 0.0)
+    weights = d_weights * g_weights
+    weights[1, 1, 1, :] = 1.0
+
+    # Normalize filters so that all sphere directions weights sum to 1
+    weights /= weights.reshape((-1, weights.shape[-1])).sum(axis=0)
+
+    return weights

scilpy/io/fetcher.py

@@ -21,12 +21,15 @@ def get_home():
 
 def get_testing_files_dict():
     """ Get dictionary linking zip file to their GDrive ID & MD5SUM """
-    return {'ihMT.zip':
+    return {'plot.zip':
+            ['1Ab-oVWI1Fu7fHTEz1H3-s1TfR_oW-GOE',
+             'cca8f1e19da357f44365a7e27b9029ca'],
+            'ihMT.zip':
             ['1V0xzvmVrVlL9dRKhc5-7xWESkmof1zyS',
-             '44eee21bcc0597836ba2eb32d41ed98c'],
+             '5d28430ac46b4fc04b6d77f9efaefb5c'],
             'MT.zip':
             ['1C2LEUkGaLFdsmym3kBrAtfPjPtv5mJuZ',
-             '66e07ada4dd2a192a44705277e0673bb'],
+             '13532c593efdf09350667df14ea4e93a'],
             'atlas.zip':
             ['1waYx4ED3qwzyJqrICjjgGXXBW2v4ZCYJ',
              'eb37427054cef5d50ac3d429ff53de47'],
@@ -38,7 +41,7 @@ def get_testing_files_dict():
              '5fbf5c8eaabff2648ad509e06b003e67'],
             'commit_amico.zip':
             ['1vyMtQd1u2h2pza9M0bncDWLc34_4MRPK',
-             '12e901e899ee48bdf31b25f22d39ee48'],
+             'b40800ab4290e4f58c375140fe59b44f'],
             'connectivity.zip':
             ['1lZqiOKmwTluPIRqblthOnBc4KI2kfKUC',
              '6d13bd076225fa2f786f416fa754623a'],
@@ -56,7 +59,7 @@ def get_testing_files_dict():
              '946beb4271b905a2bd69ad2d80136ca9'],
             'tracking.zip':
             ['1QSekZYDoMvv-An6FRMSt_s_qPeB3BHfw',
-             'eb0d6e7388b8793439857d876a3fd109'],
+             'f0492dd995f6bcb14b525a55edab9807'],
             'tractometry.zip':
             ['130mxBo4IJWPnDFyOELSYDif1puRLGHMX',
              '3e27625a1e7f2484b7fa5028c95324cc']}
@@ -137,8 +140,8 @@ def fetch_data(files_dict, keys=None):
                 break
             elif os.path.exists(full_path):
                 if tryout > 0:
-                    logging.error('Wrong md5sum after {} attemps for {}'.format(
-                        tryout+1, full_path))
+                    logging.error('Wrong md5sum after {} attemps for {}'
+                                  .format(tryout+1, full_path))
                 os.remove(full_path)
 
             # If we re-download, we re-extract

scilpy/io/streamlines.py

@@ -1,6 +1,7 @@
 # -*- coding: utf-8 -*-
 
 from itertools import islice
+import logging
 import os
 import tempfile
 
@@ -79,11 +80,20 @@ def ichunk(sequence, n):
         chunk = list(islice(sequence, n))
 
 
+def is_argument_set(args, arg_name):
+    # Check that attribute is not None
+    return not getattr(args, 'reference', None) is None
+
+
 def load_tractogram_with_reference(parser, args, filepath,
                                    bbox_check=True, arg_name=None):
 
     _, ext = os.path.splitext(filepath)
     if ext == '.trk':
+        if (is_argument_set(args, 'reference') or
+                arg_name and args.__getattribute__(arg_name + '_ref')):
+            logging.warning('Reference is discarded for this file format '
+                            '{}.'.format(filepath))
         sft = load_tractogram(filepath, 'same',
                               bbox_valid_check=bbox_check)
     elif ext in ['.tck', '.fib', '.vtk', '.dpy']:
@@ -96,10 +106,9 @@ def load_tractogram_with_reference(parser, args, filepath,
             else:
                 parser.error('--{} is required for this file format '
                              '{}.'.format(arg_ref, filepath))
-        elif args.reference is None:
+        elif (not is_argument_set(args, 'reference')) or args.reference is None:
             parser.error('--reference is required for this file format '
                          '{}.'.format(filepath))
-
         else:
             sft = load_tractogram(filepath, args.reference,
                                   bbox_valid_check=bbox_check)

scilpy/io/utils.py

@@ -8,6 +8,8 @@
 import xml.etree.ElementTree as ET
 
 import numpy as np
+from fury import window
+from PIL import Image
 from scipy.io import loadmat
 import six
 
@@ -494,3 +496,13 @@ def parser_color_type(arg):
         raise argparse.ArgumentTypeError(
             "Argument must be < " + str(MAX_VAL) + "and > " + str(MIN_VAL))
     return f
+
+
+def snapshot(scene, filename, **kwargs):
+    """ Wrapper around fury.window.snapshot
+    For some reason, fury.window.snapshot flips the image vertically.
+    This image unflips the image and then saves it.
+    """
+    out = window.snapshot(scene, **kwargs)
+    image = Image.fromarray(out[::-1])
+    image.save(filename)