From 3bee02d43592eb66ed534656cc5d4696c21e7d0c Mon Sep 17 00:00:00 2001
From: dPys <dpisner@utexas.edu>
Date: Tue, 11 Feb 2020 16:40:23 -0600
Subject: [PATCH 01/22] [ENH] Add a series of general-purpose and emc-related
 image-handling helper functions in a new module utils/images.py, and relocate
 dangling image helper functions that were previously in interfaces/images.py
 into this module

---
 dmriprep/interfaces/images.py |  77 ++++-----------
 dmriprep/utils/images.py      | 178 ++++++++++++++++++++++++++++++++++
 2 files changed, 198 insertions(+), 57 deletions(-)
 create mode 100644 dmriprep/utils/images.py

diff --git a/dmriprep/interfaces/images.py b/dmriprep/interfaces/images.py
index a59fedfb..372bdc82 100644
--- a/dmriprep/interfaces/images.py
+++ b/dmriprep/interfaces/images.py
@@ -1,10 +1,9 @@
 """Image tools interfaces."""
-import numpy as np
-import nibabel as nb
-from nipype.utils.filemanip import fname_presuffix
+from dmriprep.utils.images import rescale_b0, median, match_transforms, extract_b0
 from nipype import logging
 from nipype.interfaces.base import (
-    traits, TraitedSpec, BaseInterfaceInputSpec, SimpleInterface, File
+    traits, TraitedSpec, BaseInterfaceInputSpec, SimpleInterface, File,
+    InputMultiObject, OutputMultiObject
 )
 
 LOGGER = logging.getLogger('nipype.interface')
@@ -45,24 +44,6 @@ def _run_interface(self, runtime):
         return runtime
 
 
-def extract_b0(in_file, b0_ixs, newpath=None):
-    """Extract the *b0* volumes from a DWI dataset."""
-    out_file = fname_presuffix(
-        in_file, suffix='_b0', newpath=newpath)
-
-    img = nb.load(in_file)
-    data = img.get_fdata(dtype='float32')
-
-    b0 = data[..., b0_ixs]
-
-    hdr = img.header.copy()
-    hdr.set_data_shape(b0.shape)
-    hdr.set_xyzt_units('mm')
-    hdr.set_data_dtype(np.float32)
-    nb.Nifti1Image(b0, img.affine, hdr).to_filename(out_file)
-    return out_file
-
-
 class _RescaleB0InputSpec(BaseInterfaceInputSpec):
     in_file = File(exists=True, mandatory=True, desc='b0s file')
     mask_file = File(exists=True, mandatory=True, desc='mask file')
@@ -103,43 +84,25 @@ def _run_interface(self, runtime):
         return runtime
 
 
-def rescale_b0(in_file, mask_file, newpath=None):
-    """Rescale the input volumes using the median signal intensity."""
-    out_file = fname_presuffix(
-        in_file, suffix='_rescaled_b0', newpath=newpath)
-
-    img = nb.load(in_file)
-    if img.dataobj.ndim == 3:
-        return in_file
-
-    data = img.get_fdata(dtype='float32')
-    mask_img = nb.load(mask_file)
-    mask_data = mask_img.get_fdata(dtype='float32')
+class MatchTransformsInputSpec(BaseInterfaceInputSpec):
+    b0_indices = traits.List(mandatory=True)
+    dwi_files = InputMultiObject(File(exists=True), mandatory=True)
+    transforms = InputMultiObject(File(exists=True), mandatory=True)
 
-    median_signal = np.median(data[mask_data > 0, ...], axis=0)
-    rescaled_data = 1000 * data / median_signal
-    hdr = img.header.copy()
-    nb.Nifti1Image(rescaled_data, img.affine, hdr).to_filename(out_file)
-    return out_file
 
+class MatchTransformsOutputSpec(TraitedSpec):
+    transforms = OutputMultiObject(File(exists=True), mandatory=True)
 
-def median(in_file, newpath=None):
-    """Average a 4D dataset across the last dimension using median."""
-    out_file = fname_presuffix(
-        in_file, suffix='_b0ref', newpath=newpath)
 
-    img = nb.load(in_file)
-    if img.dataobj.ndim == 3:
-        return in_file
-    if img.shape[-1] == 1:
-        nb.squeeze_image(img).to_filename(out_file)
-        return out_file
-
-    median_data = np.median(img.get_fdata(dtype='float32'),
-                            axis=-1)
+class MatchTransforms(SimpleInterface):
+    """
+    Interface for mapping the `match_transforms` function across lists of inputs.
+    """
+    input_spec = MatchTransformsInputSpec
+    output_spec = MatchTransformsOutputSpec
 
-    hdr = img.header.copy()
-    hdr.set_xyzt_units('mm')
-    hdr.set_data_dtype(np.float32)
-    nb.Nifti1Image(median_data, img.affine, hdr).to_filename(out_file)
-    return out_file
+    def _run_interface(self, runtime):
+        self._results["transforms"] = match_transforms(
+            self.inputs.dwi_files, self.inputs.transforms, self.inputs.b0_indices
+        )
+        return runtime
diff --git a/dmriprep/utils/images.py b/dmriprep/utils/images.py
new file mode 100644
index 00000000..979e7eb5
--- /dev/null
+++ b/dmriprep/utils/images.py
@@ -0,0 +1,178 @@
+import numpy as np
+import nibabel as nb
+from nipype.utils.filemanip import fname_presuffix
+
+
+def extract_b0(in_file, b0_ixs, newpath=None):
+    """Extract the *b0* volumes from a DWI dataset."""
+    out_file = fname_presuffix(in_file, suffix="_b0", newpath=newpath)
+
+    img = nb.load(in_file)
+    data = img.get_fdata(dtype="float32")
+
+    b0 = data[..., b0_ixs]
+
+    hdr = img.header.copy()
+    hdr.set_data_shape(b0.shape)
+    hdr.set_xyzt_units("mm")
+    hdr.set_data_dtype(np.float32)
+    nb.Nifti1Image(b0, img.affine, hdr).to_filename(out_file)
+    return out_file
+
+
+def rescale_b0(in_file, mask_file, newpath=None):
+    """Rescale the input volumes using the median signal intensity."""
+    out_file = fname_presuffix(in_file, suffix="_rescaled_b0", newpath=newpath)
+
+    img = nb.load(in_file)
+    if img.dataobj.ndim == 3:
+        return in_file
+
+    data = img.get_fdata(dtype="float32")
+    mask_img = nb.load(mask_file)
+    mask_data = mask_img.get_fdata(dtype="float32")
+
+    median_signal = np.median(data[mask_data > 0, ...], axis=0)
+    rescaled_data = 1000 * data / median_signal
+    hdr = img.header.copy()
+    nb.Nifti1Image(rescaled_data, img.affine, hdr).to_filename(out_file)
+    return out_file
+
+
+def median(in_file, newpath=None):
+    """Average a 4D dataset across the last dimension using median."""
+    out_file = fname_presuffix(in_file, suffix="_b0ref", newpath=newpath)
+
+    img = nb.load(in_file)
+    if img.dataobj.ndim == 3:
+        return in_file
+    if img.shape[-1] == 1:
+        nb.squeeze_image(img).to_filename(out_file)
+        return out_file
+
+    median_data = np.median(img.get_fdata(dtype="float32"), axis=-1)
+
+    hdr = img.header.copy()
+    hdr.set_xyzt_units("mm")
+    hdr.set_data_dtype(np.float32)
+    nb.Nifti1Image(median_data, img.affine, hdr).to_filename(out_file)
+    return out_file
+
+
+def average_images(images, out_path=None):
+    """Average a 4D dataset across the last dimension using mean."""
+    from nilearn.image import mean_img
+
+    average_img = mean_img([nb.load(img) for img in images])
+    if out_path is None:
+        out_path = fname_presuffix(
+            images[0], use_ext=False, suffix="_mean.nii.gz"
+        )
+    average_img.to_filename(out_path)
+    return out_path
+
+
+def quick_load_images(image_list, dtype=np.float32):
+    """Load 3D volumes from a list of file paths into a 4D array."""
+    example_img = nb.load(image_list[0])
+    num_images = len(image_list)
+    output_matrix = np.zeros(tuple(example_img.shape) + (num_images,), dtype=dtype)
+    for image_num, image_path in enumerate(image_list):
+        output_matrix[..., image_num] = nb.load(image_path).get_fdata(dtype=dtype)
+    return output_matrix
+
+
+def match_transforms(dwi_files, transforms, b0_indices):
+    """Arranges the order of a list of affine transforms to correspond with that of
+    each individual dwi volume file, accounting for the indices of B0s. A helper
+    function for EMC."""
+    num_dwis = len(dwi_files)
+    num_transforms = len(transforms)
+
+    if num_dwis == num_transforms:
+        return transforms
+
+    # Do sanity checks
+    if not len(transforms) == len(b0_indices):
+        raise Exception("number of transforms does not match number of b0 images")
+
+    # Create a list of which emc affines go with each of the split images
+    nearest_affines = []
+    for index in range(num_dwis):
+        nearest_b0_num = np.argmin(np.abs(index - np.array(b0_indices)))
+        this_transform = transforms[nearest_b0_num]
+        nearest_affines.append(this_transform)
+
+    return nearest_affines
+
+
+def save_4d_to_3d(in_file):
+    """Split a 4D dataset along the last dimension into multiple 3D volumes."""
+    files_3d = nb.four_to_three(nb.load(in_file))
+    out_files = []
+    for i, file_3d in enumerate(files_3d):
+        out_file = fname_presuffix(in_file, suffix="_tmp_{}".format(i))
+        file_3d.to_filename(out_file)
+        out_files.append(out_file)
+    del files_3d
+    return out_files
+
+
+def prune_b0s_from_dwis(in_files, b0_ixs):
+    """Remove *b0* volume files from a complete list of DWI volume files."""
+    if in_files[0].endswith("_warped.nii.gz"):
+        out_files = [
+            i
+            for j, i in enumerate(
+                sorted(
+                    in_files, key=lambda x: int(x.split("_")[-2].split(".nii.gz")[0])
+                )
+            )
+            if j not in b0_ixs
+        ]
+    else:
+        out_files = [
+            i
+            for j, i in enumerate(
+                sorted(
+                    in_files, key=lambda x: int(x.split("_")[-1].split(".nii.gz")[0])
+                )
+            )
+            if j not in b0_ixs
+        ]
+    return out_files
+
+
+def save_3d_to_4d(in_files):
+    """Concatenate a list of 3D volumes into a 4D output."""
+    img_4d = nb.funcs.concat_images([nb.load(img_3d) for img_3d in in_files])
+    out_file = fname_presuffix(in_files[0], suffix="_merged")
+    img_4d.to_filename(out_file)
+    del img_4d
+    return out_file
+
+
+def get_params(A):
+    """Takes an transformation affine matrix A and determines
+    rotations and translations."""
+
+    def rang(b):
+        a = min(max(b, -1), 1)
+        return a
+
+    Ry = np.arcsin(A[0, 2])
+    # Rx = np.arcsin(A[1, 2] / np.cos(Ry))
+    # Rz = np.arccos(A[0, 1] / np.sin(Ry))
+
+    if (abs(Ry) - np.pi / 2) ** 2 < 1e-9:
+        Rx = 0
+        Rz = np.arctan2(-rang(A[1, 0]), rang(-A[2, 0] / A[0, 2]))
+    else:
+        c = np.cos(Ry)
+        Rx = np.arctan2(rang(A[1, 2] / c), rang(A[2, 2] / c))
+        Rz = np.arctan2(rang(A[0, 1] / c), rang(A[0, 0] / c))
+
+    rotations = [Rx, Ry, Rz]
+    translations = [A[0, 3], A[1, 3], A[2, 3]]
+
+    return rotations, translations

From 88a6e80c8b88a4d3ed4c722c04777119417b2ab2 Mon Sep 17 00:00:00 2001
From: Derek Pisner <16432683+dPys@users.noreply.github.com>
Date: Tue, 25 Feb 2020 12:39:49 -0600
Subject: [PATCH 02/22] Update dmriprep/interfaces/images.py

Co-Authored-By: Oscar Esteban <code@oscaresteban.es>
---
 dmriprep/interfaces/images.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/dmriprep/interfaces/images.py b/dmriprep/interfaces/images.py
index 372bdc82..6fa1731e 100644
--- a/dmriprep/interfaces/images.py
+++ b/dmriprep/interfaces/images.py
@@ -84,7 +84,7 @@ def _run_interface(self, runtime):
         return runtime
 
 
-class MatchTransformsInputSpec(BaseInterfaceInputSpec):
+class _MatchTransformsInputSpec(BaseInterfaceInputSpec):
     b0_indices = traits.List(mandatory=True)
     dwi_files = InputMultiObject(File(exists=True), mandatory=True)
     transforms = InputMultiObject(File(exists=True), mandatory=True)

From bb4cc79fbd90ed0af3b7419819c04d1d0131afa4 Mon Sep 17 00:00:00 2001
From: Derek Pisner <16432683+dPys@users.noreply.github.com>
Date: Tue, 25 Feb 2020 12:39:55 -0600
Subject: [PATCH 03/22] Update dmriprep/interfaces/images.py

Co-Authored-By: Oscar Esteban <code@oscaresteban.es>
---
 dmriprep/interfaces/images.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/dmriprep/interfaces/images.py b/dmriprep/interfaces/images.py
index 6fa1731e..4d9a0cbf 100644
--- a/dmriprep/interfaces/images.py
+++ b/dmriprep/interfaces/images.py
@@ -90,7 +90,7 @@ class _MatchTransformsInputSpec(BaseInterfaceInputSpec):
     transforms = InputMultiObject(File(exists=True), mandatory=True)
 
 
-class MatchTransformsOutputSpec(TraitedSpec):
+class _MatchTransformsOutputSpec(TraitedSpec):
     transforms = OutputMultiObject(File(exists=True), mandatory=True)
 
 

From 1e9c015707559231f85e681a7fd494357d516fdd Mon Sep 17 00:00:00 2001
From: Derek Pisner <16432683+dPys@users.noreply.github.com>
Date: Tue, 25 Feb 2020 12:40:16 -0600
Subject: [PATCH 04/22] Update dmriprep/utils/images.py

Co-Authored-By: Oscar Esteban <code@oscaresteban.es>
---
 dmriprep/utils/images.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/dmriprep/utils/images.py b/dmriprep/utils/images.py
index 979e7eb5..b55ddd55 100644
--- a/dmriprep/utils/images.py
+++ b/dmriprep/utils/images.py
@@ -152,7 +152,7 @@ def save_3d_to_4d(in_files):
     return out_file
 
 
-def get_params(A):
+def decompose_affine(affine):
     """Takes an transformation affine matrix A and determines
     rotations and translations."""
 

From e4ea79836417752a9202f2124679ec736f6064a8 Mon Sep 17 00:00:00 2001
From: Derek Pisner <16432683+dPys@users.noreply.github.com>
Date: Tue, 25 Feb 2020 12:40:32 -0600
Subject: [PATCH 05/22] Update dmriprep/utils/images.py

Co-Authored-By: Oscar Esteban <code@oscaresteban.es>
---
 dmriprep/utils/images.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/dmriprep/utils/images.py b/dmriprep/utils/images.py
index b55ddd55..94f3f6fd 100644
--- a/dmriprep/utils/images.py
+++ b/dmriprep/utils/images.py
@@ -50,7 +50,7 @@ def median(in_file, newpath=None):
         nb.squeeze_image(img).to_filename(out_file)
         return out_file
 
-    median_data = np.median(img.get_fdata(dtype="float32"), axis=-1)
+    median_data = np.median(img.get_fdata(), axis=-1)
 
     hdr = img.header.copy()
     hdr.set_xyzt_units("mm")

From e982d6f4eb6db387ec911310782cbc4da5e321b4 Mon Sep 17 00:00:00 2001
From: Derek Pisner <16432683+dPys@users.noreply.github.com>
Date: Tue, 25 Feb 2020 12:40:39 -0600
Subject: [PATCH 06/22] Update dmriprep/utils/images.py

Co-Authored-By: Oscar Esteban <code@oscaresteban.es>
---
 dmriprep/utils/images.py | 1 -
 1 file changed, 1 deletion(-)

diff --git a/dmriprep/utils/images.py b/dmriprep/utils/images.py
index 94f3f6fd..cff8f5c1 100644
--- a/dmriprep/utils/images.py
+++ b/dmriprep/utils/images.py
@@ -54,7 +54,6 @@ def median(in_file, newpath=None):
 
     hdr = img.header.copy()
     hdr.set_xyzt_units("mm")
-    hdr.set_data_dtype(np.float32)
     nb.Nifti1Image(median_data, img.affine, hdr).to_filename(out_file)
     return out_file
 

From 35889474d3f9d59e086aa3e452bdfc99f04f5f55 Mon Sep 17 00:00:00 2001
From: Derek Pisner <16432683+dPys@users.noreply.github.com>
Date: Tue, 25 Feb 2020 13:16:43 -0600
Subject: [PATCH 07/22] Update dmriprep/utils/images.py

Co-Authored-By: Oscar Esteban <code@oscaresteban.es>
---
 dmriprep/utils/images.py | 1 -
 1 file changed, 1 deletion(-)

diff --git a/dmriprep/utils/images.py b/dmriprep/utils/images.py
index cff8f5c1..7444be5e 100644
--- a/dmriprep/utils/images.py
+++ b/dmriprep/utils/images.py
@@ -147,7 +147,6 @@ def save_3d_to_4d(in_files):
     img_4d = nb.funcs.concat_images([nb.load(img_3d) for img_3d in in_files])
     out_file = fname_presuffix(in_files[0], suffix="_merged")
     img_4d.to_filename(out_file)
-    del img_4d
     return out_file
 
 

From 1544a55624179d5e50c14106d27d5f109ba2ee57 Mon Sep 17 00:00:00 2001
From: Derek Pisner <16432683+dPys@users.noreply.github.com>
Date: Tue, 17 Mar 2020 13:25:19 -0500
Subject: [PATCH 08/22] Update dmriprep/interfaces/images.py

Co-Authored-By: Oscar Esteban <code@oscaresteban.es>
---
 dmriprep/interfaces/images.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/dmriprep/interfaces/images.py b/dmriprep/interfaces/images.py
index 4d9a0cbf..7d9b9841 100644
--- a/dmriprep/interfaces/images.py
+++ b/dmriprep/interfaces/images.py
@@ -98,8 +98,8 @@ class MatchTransforms(SimpleInterface):
     """
     Interface for mapping the `match_transforms` function across lists of inputs.
     """
-    input_spec = MatchTransformsInputSpec
-    output_spec = MatchTransformsOutputSpec
+    input_spec = _MatchTransformsInputSpec
+    output_spec = _MatchTransformsOutputSpec
 
     def _run_interface(self, runtime):
         self._results["transforms"] = match_transforms(

From b6e917e8e802c215d81e44ea80ad8bc716fdf2eb Mon Sep 17 00:00:00 2001
From: Derek Pisner <16432683+dPys@users.noreply.github.com>
Date: Tue, 17 Mar 2020 13:25:46 -0500
Subject: [PATCH 09/22] Update dmriprep/utils/images.py

Co-Authored-By: Oscar Esteban <code@oscaresteban.es>
---
 dmriprep/utils/images.py | 7 +------
 1 file changed, 1 insertion(+), 6 deletions(-)

diff --git a/dmriprep/utils/images.py b/dmriprep/utils/images.py
index 7444be5e..28311cea 100644
--- a/dmriprep/utils/images.py
+++ b/dmriprep/utils/images.py
@@ -73,12 +73,7 @@ def average_images(images, out_path=None):
 
 def quick_load_images(image_list, dtype=np.float32):
     """Load 3D volumes from a list of file paths into a 4D array."""
-    example_img = nb.load(image_list[0])
-    num_images = len(image_list)
-    output_matrix = np.zeros(tuple(example_img.shape) + (num_images,), dtype=dtype)
-    for image_num, image_path in enumerate(image_list):
-        output_matrix[..., image_num] = nb.load(image_path).get_fdata(dtype=dtype)
-    return output_matrix
+    return nb.concat_images([nb.load(fname) for fname in image_list]).get_fdata(dtype=dtype)
 
 
 def match_transforms(dwi_files, transforms, b0_indices):

From a888d736f412496cffde1a5ab17901e055edf469 Mon Sep 17 00:00:00 2001
From: Derek Pisner <16432683+dPys@users.noreply.github.com>
Date: Tue, 17 Mar 2020 13:26:07 -0500
Subject: [PATCH 10/22] Update dmriprep/utils/images.py

Co-Authored-By: Oscar Esteban <code@oscaresteban.es>
---
 dmriprep/utils/images.py | 11 ++++++++---
 1 file changed, 8 insertions(+), 3 deletions(-)

diff --git a/dmriprep/utils/images.py b/dmriprep/utils/images.py
index 28311cea..01c0c29c 100644
--- a/dmriprep/utils/images.py
+++ b/dmriprep/utils/images.py
@@ -77,9 +77,14 @@ def quick_load_images(image_list, dtype=np.float32):
 
 
 def match_transforms(dwi_files, transforms, b0_indices):
-    """Arranges the order of a list of affine transforms to correspond with that of
-    each individual dwi volume file, accounting for the indices of B0s. A helper
-    function for EMC."""
+    """
+    Arrange the order of a list of transforms.
+    
+    This is a helper function for :abbr:`EMC (Eddy-currents and Motion Correction)`.
+    Sorts the input list of affine transforms  to correspond with that of
+    each individual dwi volume file, accounting for the indices of :math:`b = 0` volumes.
+    
+    """
     num_dwis = len(dwi_files)
     num_transforms = len(transforms)
 

From d8720949c7224d650d955d6cd819079a880c59aa Mon Sep 17 00:00:00 2001
From: dPys <dpisner@utexas.edu>
Date: Tue, 17 Mar 2020 14:28:28 -0500
Subject: [PATCH 11/22] [ENH] Add docstring and first half of unit tests

---
 dmriprep/data/tests/dwi_b0.nii.gz | Bin 0 -> 128 bytes
 dmriprep/utils/images.py          | 251 +++++++++++++++++++++++-------
 2 files changed, 197 insertions(+), 54 deletions(-)
 create mode 100644 dmriprep/data/tests/dwi_b0.nii.gz

diff --git a/dmriprep/data/tests/dwi_b0.nii.gz b/dmriprep/data/tests/dwi_b0.nii.gz
new file mode 100644
index 0000000000000000000000000000000000000000..f0bb4d317ae16ab08373712d85a17182574bfc14
GIT binary patch
literal 128
zcmV-`0Du1<iwFoABynB>|73S*Ut%yWZfR)%i(zCS3NSNp0U;9uBNU=)5@1jO3pCg>
zFdRU_AU-nQ;0)n|6*HllL4O{|F9^_JKWpYWr(boKoVxn*4gvkKA0`J9hp-`HAT=;K
ikbfcS5h5T~p0**Xz$kAt1V%#uBLo1yw1tv@1ONc21S`4#

literal 0
HcmV?d00001

diff --git a/dmriprep/utils/images.py b/dmriprep/utils/images.py
index 01c0c29c..f89b7e04 100644
--- a/dmriprep/utils/images.py
+++ b/dmriprep/utils/images.py
@@ -1,11 +1,42 @@
+import os
 import numpy as np
 import nibabel as nb
 from nipype.utils.filemanip import fname_presuffix
 
 
-def extract_b0(in_file, b0_ixs, newpath=None):
-    """Extract the *b0* volumes from a DWI dataset."""
-    out_file = fname_presuffix(in_file, suffix="_b0", newpath=newpath)
+def extract_b0(in_file, b0_ixs, out_path=None):
+    """
+    Extract the *b0* volumes from a DWI dataset.
+
+    Parameters
+    ----------
+    in_file : str
+        DWI NIfTI file.
+    b0_ixs : list
+        List of B0 indices in `in_file`.
+    out_path : str
+        Optionally specify an output path.
+
+    Returns
+    -------
+    out_path : str
+       4D NIfTI file consisting of B0's.
+
+    Examples
+    --------
+    >>> os.chdir(tmpdir)
+    >>> b0_ixs = np.where(np.loadtxt(str(data_dir / 'bval')) <= 50)[0].tolist()[:2]
+    >>> in_file = str(data_dir / 'dwi.nii.gz')
+    >>> out_path = extract_b0(in_file, b0_ixs)
+    >>> assert os.path.isfile(out_path)
+    """
+    if out_path is None:
+        out_path = fname_presuffix(
+            in_file, suffix="_b0"
+        )
+    else:
+        out_path = fname_presuffix(in_file, suffix="_b0",
+                                   newpath=out_path)
 
     img = nb.load(in_file)
     data = img.get_fdata(dtype="float32")
@@ -16,13 +47,43 @@ def extract_b0(in_file, b0_ixs, newpath=None):
     hdr.set_data_shape(b0.shape)
     hdr.set_xyzt_units("mm")
     hdr.set_data_dtype(np.float32)
-    nb.Nifti1Image(b0, img.affine, hdr).to_filename(out_file)
-    return out_file
+    nb.Nifti1Image(b0, img.affine, hdr).to_filename(out_path)
+    return out_path
 
 
-def rescale_b0(in_file, mask_file, newpath=None):
-    """Rescale the input volumes using the median signal intensity."""
-    out_file = fname_presuffix(in_file, suffix="_rescaled_b0", newpath=newpath)
+def rescale_b0(in_file, mask_file, out_path=None):
+    """
+    Rescale the input volumes using the median signal intensity.
+
+    Parameters
+    ----------
+    in_file : str
+        A NIfTI file consisting of one or more B0's.
+    mask_file : str
+        A B0 mask NIFTI file.
+    out_path : str
+        Optionally specify an output path.
+
+    Returns
+    -------
+    out_path : str
+       A rescaled B0 NIFTI file.
+
+    Examples
+    --------
+    >>> os.chdir(tmpdir)
+    >>> mask_file = str(data_dir / 'dwi_mask.nii.gz')
+    >>> in_file = str(data_dir / 'dwi_b0.nii.gz')
+    >>> out_path = rescale_b0(in_file, mask_file)
+    >>> assert os.path.isfile(out_path)
+    """
+    if out_path is None:
+        out_path = fname_presuffix(
+            in_file, suffix="_rescaled_b0"
+        )
+    else:
+        out_path = fname_presuffix(in_file, suffix="_rescaled_b0",
+                                   newpath=out_path)
 
     img = nb.load(in_file)
     if img.dataobj.ndim == 3:
@@ -35,31 +96,79 @@ def rescale_b0(in_file, mask_file, newpath=None):
     median_signal = np.median(data[mask_data > 0, ...], axis=0)
     rescaled_data = 1000 * data / median_signal
     hdr = img.header.copy()
-    nb.Nifti1Image(rescaled_data, img.affine, hdr).to_filename(out_file)
-    return out_file
+    nb.Nifti1Image(rescaled_data, img.affine, hdr).to_filename(out_path)
+    return out_path
 
 
-def median(in_file, newpath=None):
-    """Average a 4D dataset across the last dimension using median."""
-    out_file = fname_presuffix(in_file, suffix="_b0ref", newpath=newpath)
+def median(in_file, out_path=None):
+    """
+    Summarize a 4D dataset across the last dimension using median.
+
+    Parameters
+    ----------
+    in_file : str
+        A NIfTI file consisting of one or more B0's.
+    out_path : str
+        Optionally specify an output path for `out_path`.
+
+    Returns
+    -------
+    out_path : str
+       A 3D B0 NIFTI file.
+
+    Examples
+    --------
+    >>> os.chdir(tmpdir)
+    >>> in_file = str(data_dir / 'dwi_b0.nii.gz')
+    >>> out_path = median(in_file)
+    >>> assert os.path.isfile(out_path)
+    """
+    if out_path is None:
+        out_path = fname_presuffix(
+            in_file, suffix="_b0ref"
+        )
+    else:
+        out_path = fname_presuffix(in_file, suffix="_b0ref",
+                                   newpath=out_path)
 
     img = nb.load(in_file)
     if img.dataobj.ndim == 3:
         return in_file
     if img.shape[-1] == 1:
-        nb.squeeze_image(img).to_filename(out_file)
-        return out_file
+        nb.squeeze_image(img).to_filename(out_path)
+        return out_path
 
     median_data = np.median(img.get_fdata(), axis=-1)
 
     hdr = img.header.copy()
     hdr.set_xyzt_units("mm")
-    nb.Nifti1Image(median_data, img.affine, hdr).to_filename(out_file)
-    return out_file
+    nb.Nifti1Image(median_data, img.affine, hdr).to_filename(out_path)
+    return out_path
 
 
 def average_images(images, out_path=None):
-    """Average a 4D dataset across the last dimension using mean."""
+    """
+    Average a 4D dataset across the last dimension using mean.
+
+    Parameters
+    ----------
+    images : str
+        A list of NIFTI file paths.
+    out_path : str
+        Optionally specify an output path for `out_path`.
+
+    Returns
+    -------
+    out_path : str
+       A 3D NIFTI file averaged along the 4th dimension.
+
+    Examples
+    --------
+    >>> os.chdir(tmpdir)
+    >>> in_file = str(data_dir / 'dwi_b0.nii.gz')
+    >>> out_path = average_images(in_file)
+    >>> assert os.path.isfile(out_path)
+    """
     from nilearn.image import mean_img
 
     average_img = mean_img([nb.load(img) for img in images])
@@ -71,19 +180,39 @@ def average_images(images, out_path=None):
     return out_path
 
 
-def quick_load_images(image_list, dtype=np.float32):
-    """Load 3D volumes from a list of file paths into a 4D array."""
-    return nb.concat_images([nb.load(fname) for fname in image_list]).get_fdata(dtype=dtype)
+def get_list_data(file_list, dtype=np.float32):
+    """
+    Load 3D volumes from a list of file paths into a 4D array.
+
+    Parameters
+    ----------
+    file_list : str
+        A list of file paths to 3D NIFTI images.
+    """
+    return nb.concat_images([nb.load(fname) for fname in file_list]).get_fdata(dtype=dtype)
 
 
-def match_transforms(dwi_files, transforms, b0_indices):
+def match_transforms(dwi_files, transforms, b0_ixs):
     """
     Arrange the order of a list of transforms.
-    
+
     This is a helper function for :abbr:`EMC (Eddy-currents and Motion Correction)`.
     Sorts the input list of affine transforms  to correspond with that of
     each individual dwi volume file, accounting for the indices of :math:`b = 0` volumes.
-    
+
+    Parameters
+    ----------
+    dwi_files : list
+        A list of file paths to 3D diffusion-weighted NIFTI volumes.
+    transforms : list
+        A list of ndarrays.
+    b0_ixs : list
+        List of B0 indices.
+
+    Returns
+    -------
+    out_path : str
+       A 3D NIFTI file averaged along the 4th dimension.
     """
     num_dwis = len(dwi_files)
     num_transforms = len(transforms)
@@ -92,13 +221,13 @@ def match_transforms(dwi_files, transforms, b0_indices):
         return transforms
 
     # Do sanity checks
-    if not len(transforms) == len(b0_indices):
+    if not len(transforms) == len(b0_ixs):
         raise Exception("number of transforms does not match number of b0 images")
 
     # Create a list of which emc affines go with each of the split images
     nearest_affines = []
     for index in range(num_dwis):
-        nearest_b0_num = np.argmin(np.abs(index - np.array(b0_indices)))
+        nearest_b0_num = np.argmin(np.abs(index - np.array(b0_ixs)))
         this_transform = transforms[nearest_b0_num]
         nearest_affines.append(this_transform)
 
@@ -106,7 +235,19 @@ def match_transforms(dwi_files, transforms, b0_indices):
 
 
 def save_4d_to_3d(in_file):
-    """Split a 4D dataset along the last dimension into multiple 3D volumes."""
+    """
+    Split a 4D dataset along the last dimension into multiple 3D volumes.
+
+    Parameters
+    ----------
+    in_file : str
+        DWI NIfTI file.
+
+    Returns
+    -------
+    out_files : list
+       A list of file paths to 3d NIFTI images.
+    """
     files_3d = nb.four_to_three(nb.load(in_file))
     out_files = []
     for i, file_3d in enumerate(files_3d):
@@ -118,7 +259,22 @@ def save_4d_to_3d(in_file):
 
 
 def prune_b0s_from_dwis(in_files, b0_ixs):
-    """Remove *b0* volume files from a complete list of DWI volume files."""
+    """
+    Remove *b0* volume files from a complete list of DWI volume files.
+
+    Parameters
+    ----------
+    in_files : list
+        A list of NIfTI file paths corresponding to each 3D volume of a
+        DWI image (i.e. including B0's).
+    b0_ixs : list
+        List of B0 indices.
+
+    Returns
+    -------
+    out_files : list
+       A list of file paths to 3d NIFTI images.
+    """
     if in_files[0].endswith("_warped.nii.gz"):
         out_files = [
             i
@@ -143,34 +299,21 @@ def prune_b0s_from_dwis(in_files, b0_ixs):
 
 
 def save_3d_to_4d(in_files):
-    """Concatenate a list of 3D volumes into a 4D output."""
+    """
+    Concatenate a list of 3D volumes into a 4D output.
+
+    Parameters
+    ----------
+    in_files : list
+        A list of file paths to 3D NIFTI images.
+
+    Returns
+    -------
+    out_file : str
+       A file path to a 4d NIFTI image of concatenated 3D volumes.
+    """
     img_4d = nb.funcs.concat_images([nb.load(img_3d) for img_3d in in_files])
     out_file = fname_presuffix(in_files[0], suffix="_merged")
     img_4d.to_filename(out_file)
     return out_file
 
-
-def decompose_affine(affine):
-    """Takes an transformation affine matrix A and determines
-    rotations and translations."""
-
-    def rang(b):
-        a = min(max(b, -1), 1)
-        return a
-
-    Ry = np.arcsin(A[0, 2])
-    # Rx = np.arcsin(A[1, 2] / np.cos(Ry))
-    # Rz = np.arccos(A[0, 1] / np.sin(Ry))
-
-    if (abs(Ry) - np.pi / 2) ** 2 < 1e-9:
-        Rx = 0
-        Rz = np.arctan2(-rang(A[1, 0]), rang(-A[2, 0] / A[0, 2]))
-    else:
-        c = np.cos(Ry)
-        Rx = np.arctan2(rang(A[1, 2] / c), rang(A[2, 2] / c))
-        Rz = np.arctan2(rang(A[0, 1] / c), rang(A[0, 0] / c))
-
-    rotations = [Rx, Ry, Rz]
-    translations = [A[0, 3], A[1, 3], A[2, 3]]
-
-    return rotations, translations

From eb430430ec208194300e3bacb7d427c900d00cef Mon Sep 17 00:00:00 2001
From: dPys <dpisner@utexas.edu>
Date: Tue, 17 Mar 2020 14:32:49 -0500
Subject: [PATCH 12/22] Address Sider complaint

---
 dmriprep/utils/images.py | 2 --
 1 file changed, 2 deletions(-)

diff --git a/dmriprep/utils/images.py b/dmriprep/utils/images.py
index f89b7e04..b2b92234 100644
--- a/dmriprep/utils/images.py
+++ b/dmriprep/utils/images.py
@@ -1,4 +1,3 @@
-import os
 import numpy as np
 import nibabel as nb
 from nipype.utils.filemanip import fname_presuffix
@@ -316,4 +315,3 @@ def save_3d_to_4d(in_files):
     out_file = fname_presuffix(in_files[0], suffix="_merged")
     img_4d.to_filename(out_file)
     return out_file
-

From abf5186d9d6f1099e22041d03b78916662dbeb6c Mon Sep 17 00:00:00 2001
From: dPys <dpisner@utexas.edu>
Date: Tue, 24 Mar 2020 11:46:38 -0500
Subject: [PATCH 13/22] [DOC] Add doctests for all functions using HARDI data

---
 dmriprep/data/tests/dwi_b0.nii.gz | Bin 128 -> 128 bytes
 dmriprep/utils/images.py          |  84 +++++++++++++++++++++++-------
 2 files changed, 65 insertions(+), 19 deletions(-)

diff --git a/dmriprep/data/tests/dwi_b0.nii.gz b/dmriprep/data/tests/dwi_b0.nii.gz
index f0bb4d317ae16ab08373712d85a17182574bfc14..b8a36ac29e01e3e6513b383457886fdc4d206cb2 100644
GIT binary patch
delta 15
WcmZo*Y+z)U@8;l8w5*!QUI_pj3<IeE

delta 15
WcmZo*Y+z)U@8;l;R4JUuUI_piPXmJh

diff --git a/dmriprep/utils/images.py b/dmriprep/utils/images.py
index b2b92234..6bac745a 100644
--- a/dmriprep/utils/images.py
+++ b/dmriprep/utils/images.py
@@ -1,3 +1,4 @@
+from pathlib import Path
 import numpy as np
 import nibabel as nb
 from nipype.utils.filemanip import fname_presuffix
@@ -31,11 +32,7 @@ def extract_b0(in_file, b0_ixs, out_path=None):
     """
     if out_path is None:
         out_path = fname_presuffix(
-            in_file, suffix="_b0"
-        )
-    else:
-        out_path = fname_presuffix(in_file, suffix="_b0",
-                                   newpath=out_path)
+            in_file, suffix='_b0', use_ext=True)
 
     img = nb.load(in_file)
     data = img.get_fdata(dtype="float32")
@@ -78,11 +75,7 @@ def rescale_b0(in_file, mask_file, out_path=None):
     """
     if out_path is None:
         out_path = fname_presuffix(
-            in_file, suffix="_rescaled_b0"
-        )
-    else:
-        out_path = fname_presuffix(in_file, suffix="_rescaled_b0",
-                                   newpath=out_path)
+            in_file, suffix='_rescaled', use_ext=True)
 
     img = nb.load(in_file)
     if img.dataobj.ndim == 3:
@@ -124,11 +117,7 @@ def median(in_file, out_path=None):
     """
     if out_path is None:
         out_path = fname_presuffix(
-            in_file, suffix="_b0ref"
-        )
-    else:
-        out_path = fname_presuffix(in_file, suffix="_b0ref",
-                                   newpath=out_path)
+            in_file, suffix='_b0ref', use_ext=True)
 
     img = nb.load(in_file)
     if img.dataobj.ndim == 3:
@@ -164,8 +153,9 @@ def average_images(images, out_path=None):
     Examples
     --------
     >>> os.chdir(tmpdir)
-    >>> in_file = str(data_dir / 'dwi_b0.nii.gz')
-    >>> out_path = average_images(in_file)
+    >>> in_file = str(dipy_datadir / "HARDI193.nii.gz")
+    >>> out_files = save_4d_to_3d(in_file)
+    >>> out_path = average_images(out_files)
     >>> assert os.path.isfile(out_path)
     """
     from nilearn.image import mean_img
@@ -187,6 +177,13 @@ def get_list_data(file_list, dtype=np.float32):
     ----------
     file_list : str
         A list of file paths to 3D NIFTI images.
+
+    Examples
+    --------
+    >>> os.chdir(tmpdir)
+    >>> in_file = str(dipy_datadir / "HARDI193.nii.gz")
+    >>> out_files = save_4d_to_3d(in_file)
+    >>> assert len(out_files) == get_list_data(out_files).shape[-1]
     """
     return nb.concat_images([nb.load(fname) for fname in file_list]).get_fdata(dtype=dtype)
 
@@ -210,8 +207,32 @@ def match_transforms(dwi_files, transforms, b0_ixs):
 
     Returns
     -------
-    out_path : str
-       A 3D NIFTI file averaged along the 4th dimension.
+    nearest_affines : list
+       A list of affine file paths that correspond to each of the split
+       dwi volumes.
+
+    Examples
+    --------
+    >>> os.chdir(tmpdir)
+    >>> from dmriprep.utils.vectors import DiffusionGradientTable
+    >>> dwi_file = str(dipy_datadir / "HARDI193.nii.gz")
+    >>> check = DiffusionGradientTable(
+    ...     dwi_file=dwi_file,
+    ...     bvecs=str(dipy_datadir / "HARDI193.bvec"),
+    ...     bvals=str(dipy_datadir / "HARDI193.bval"))
+    >>> check.generate_rasb()
+    >>> # Conform to the orientation of the image:
+    >>> affines = np.zeros((check.gradients.shape[0], 4, 4))
+    >>> transforms = []
+    >>> for ii, aff in enumerate(affines):
+    ...     aff_file = f'aff_{ii}.npy'
+    ...     np.save(aff_file, aff)
+    ...     transforms.append(aff_file)
+    >>> dwi_files = save_4d_to_3d(dwi_file)
+    >>> b0_ixs = np.where((check.bvals) <= 50)[0].tolist()[:2]
+    >>> nearest_affines = match_transforms(dwi_files, transforms, b0_ixs)
+    >>> assert sum([os.path.isfile(i) for i in nearest_affines]) == len(nearest_affines)
+    >>> assert len(nearest_affines) == len(dwi_files)
     """
     num_dwis = len(dwi_files)
     num_transforms = len(transforms)
@@ -246,6 +267,13 @@ def save_4d_to_3d(in_file):
     -------
     out_files : list
        A list of file paths to 3d NIFTI images.
+
+    Examples
+    --------
+    >>> os.chdir(tmpdir)
+    >>> in_file = str(dipy_datadir / "HARDI193.nii.gz")
+    >>> out_files = save_4d_to_3d(in_file)
+    >>> assert len(out_files) == nb.load(in_file).shape[-1]
     """
     files_3d = nb.four_to_three(nb.load(in_file))
     out_files = []
@@ -273,6 +301,16 @@ def prune_b0s_from_dwis(in_files, b0_ixs):
     -------
     out_files : list
        A list of file paths to 3d NIFTI images.
+
+    Examples
+    --------
+    >>> os.chdir(tmpdir)
+    >>> b0_ixs = np.where(np.loadtxt(str(dipy_datadir / "HARDI193.bval")) <= 50)[0].tolist()[:2]
+    >>> in_file = str(dipy_datadir / "HARDI193.nii.gz")
+    >>> threeD_files = save_4d_to_3d(in_file)
+    >>> out_files = prune_b0s_from_dwis(threeD_files, b0_ixs)
+    >>> assert sum([os.path.isfile(i) for i in out_files]) == len(out_files)
+    >>> assert len(out_files) == len(threeD_files) - len(b0_ixs)
     """
     if in_files[0].endswith("_warped.nii.gz"):
         out_files = [
@@ -310,6 +348,14 @@ def save_3d_to_4d(in_files):
     -------
     out_file : str
        A file path to a 4d NIFTI image of concatenated 3D volumes.
+
+    Examples
+    --------
+    >>> os.chdir(tmpdir)
+    >>> in_file = str(dipy_datadir / "HARDI193.nii.gz")
+    >>> threeD_files = save_4d_to_3d(in_file)
+    >>> out_file = save_3d_to_4d(threeD_files)
+    >>> assert len(threeD_files) == nb.load(out_file).shape[-1]
     """
     img_4d = nb.funcs.concat_images([nb.load(img_3d) for img_3d in in_files])
     out_file = fname_presuffix(in_files[0], suffix="_merged")

From 02f82ceffab57d12de3ed61228c242d590e475b4 Mon Sep 17 00:00:00 2001
From: Derek Pisner <16432683+dPys@users.noreply.github.com>
Date: Tue, 24 Mar 2020 11:58:04 -0500
Subject: [PATCH 14/22] Remove unused Pathlib import

---
 dmriprep/utils/images.py | 1 -
 1 file changed, 1 deletion(-)

diff --git a/dmriprep/utils/images.py b/dmriprep/utils/images.py
index d90c1095..6b6925ac 100644
--- a/dmriprep/utils/images.py
+++ b/dmriprep/utils/images.py
@@ -1,4 +1,3 @@
-from pathlib import Path
 import numpy as np
 import nibabel as nb
 from nipype.utils.filemanip import fname_presuffix

From bb6c9a57957ec6c77c79ca9674622e207d534831 Mon Sep 17 00:00:00 2001
From: Derek Pisner <16432683+dPys@users.noreply.github.com>
Date: Tue, 24 Mar 2020 11:58:45 -0500
Subject: [PATCH 15/22] Add extra blank line to make sidecar happy

---
 dmriprep/interfaces/images.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/dmriprep/interfaces/images.py b/dmriprep/interfaces/images.py
index 7596448a..48f89c45 100644
--- a/dmriprep/interfaces/images.py
+++ b/dmriprep/interfaces/images.py
@@ -106,6 +106,7 @@ def _run_interface(self, runtime):
         )
         return runtime
 
+    
 class _MatchTransformsInputSpec(BaseInterfaceInputSpec):
     b0_indices = traits.List(mandatory=True)
     dwi_files = InputMultiObject(File(exists=True), mandatory=True)

From 56213a0384a9eb538cf672b8a45b2d99d11cf3f1 Mon Sep 17 00:00:00 2001
From: Derek Pisner <16432683+dPys@users.noreply.github.com>
Date: Tue, 24 Mar 2020 12:12:59 -0500
Subject: [PATCH 16/22] Fix whitespace (again)

---
 dmriprep/interfaces/images.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/dmriprep/interfaces/images.py b/dmriprep/interfaces/images.py
index 48f89c45..7e76a730 100644
--- a/dmriprep/interfaces/images.py
+++ b/dmriprep/interfaces/images.py
@@ -106,7 +106,7 @@ def _run_interface(self, runtime):
         )
         return runtime
 
-    
+
 class _MatchTransformsInputSpec(BaseInterfaceInputSpec):
     b0_indices = traits.List(mandatory=True)
     dwi_files = InputMultiObject(File(exists=True), mandatory=True)

From 5e64dfafb970e736eab053d7645c4c79bba9fe74 Mon Sep 17 00:00:00 2001
From: Derek Pisner <16432683+dPys@users.noreply.github.com>
Date: Tue, 24 Mar 2020 12:14:11 -0500
Subject: [PATCH 17/22] Update dmriprep/utils/images.py

Co-Authored-By: Ariel Rokem <arokem@gmail.com>
---
 dmriprep/utils/images.py | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/dmriprep/utils/images.py b/dmriprep/utils/images.py
index 6b6925ac..fe3d0379 100644
--- a/dmriprep/utils/images.py
+++ b/dmriprep/utils/images.py
@@ -180,7 +180,9 @@ def get_list_data(file_list, dtype=np.float32):
     ----------
     file_list : str
         A list of file paths to 3D NIFTI images.
-
+Returns 
+--------
+Nibabel image object
     Examples
     --------
     >>> os.chdir(tmpdir)

From 85a285a80e40d89508ab80613499f371a79c9c55 Mon Sep 17 00:00:00 2001
From: dPys <dpisner@utexas.edu>
Date: Tue, 24 Mar 2020 12:19:34 -0500
Subject: [PATCH 18/22] [ENH] Add error handling to

---
 dmriprep/interfaces/images.py |  5 ++---
 dmriprep/utils/images.py      | 13 ++++++++++++-
 2 files changed, 14 insertions(+), 4 deletions(-)

diff --git a/dmriprep/interfaces/images.py b/dmriprep/interfaces/images.py
index 7e76a730..40a0b0a3 100644
--- a/dmriprep/interfaces/images.py
+++ b/dmriprep/interfaces/images.py
@@ -1,12 +1,11 @@
 """Image tools interfaces."""
-from dmriprep.utils.images import rescale_b0, median, match_transforms, extract_b0
+from pathlib import Path
 from nipype import logging
 from nipype.interfaces.base import (
     traits, TraitedSpec, BaseInterfaceInputSpec, SimpleInterface, File,
     InputMultiObject, OutputMultiObject
 )
-
-from dmriprep.utils.images import extract_b0, median, rescale_b0
+from dmriprep.utils.images import extract_b0, median, rescale_b0, match_transforms
 
 LOGGER = logging.getLogger('nipype.interface')
 
diff --git a/dmriprep/utils/images.py b/dmriprep/utils/images.py
index fe3d0379..c176e091 100644
--- a/dmriprep/utils/images.py
+++ b/dmriprep/utils/images.py
@@ -362,7 +362,18 @@ def save_3d_to_4d(in_files):
     >>> out_file = save_3d_to_4d(threeD_files)
     >>> assert len(threeD_files) == nb.load(out_file).shape[-1]
     """
-    img_4d = nb.funcs.concat_images([nb.load(img_3d) for img_3d in in_files])
+    # Remove one-sized extra dimensions
+    nii_list = []
+    for i, f in enumerate(in_files):
+        filenii = nb.load(f)
+        filenii = nb.squeeze_image(filenii)
+        if len(filenii.shape) == 5:
+            raise RuntimeError('Input image (%s) is 5D.' % f)
+        if filenii.dataobj.ndim == 4:
+            nii_list += nb.four_to_three(filenii)
+        else:
+            nii_list.append(filenii)
+    img_4d = nb.funcs.concat_images(nii_list)
     out_file = fname_presuffix(in_files[0], suffix="_merged")
     img_4d.to_filename(out_file)
     return out_file

From b13d57bb287594592c93d389fe0f4b91d0134d0a Mon Sep 17 00:00:00 2001
From: dPys <dpisner@utexas.edu>
Date: Tue, 24 Mar 2020 12:25:50 -0500
Subject: [PATCH 19/22] omit prune_b0s helper function, save for later PR with
 signal prediction

---
 dmriprep/utils/images.py | 58 ++++------------------------------------
 1 file changed, 5 insertions(+), 53 deletions(-)

diff --git a/dmriprep/utils/images.py b/dmriprep/utils/images.py
index c176e091..b2f297e6 100644
--- a/dmriprep/utils/images.py
+++ b/dmriprep/utils/images.py
@@ -180,9 +180,11 @@ def get_list_data(file_list, dtype=np.float32):
     ----------
     file_list : str
         A list of file paths to 3D NIFTI images.
-Returns 
---------
-Nibabel image object
+
+    Returns
+    --------
+    Nibabel image object
+
     Examples
     --------
     >>> os.chdir(tmpdir)
@@ -290,56 +292,6 @@ def save_4d_to_3d(in_file):
     return out_files
 
 
-def prune_b0s_from_dwis(in_files, b0_ixs):
-    """
-    Remove *b0* volume files from a complete list of DWI volume files.
-
-    Parameters
-    ----------
-    in_files : list
-        A list of NIfTI file paths corresponding to each 3D volume of a
-        DWI image (i.e. including B0's).
-    b0_ixs : list
-        List of B0 indices.
-
-    Returns
-    -------
-    out_files : list
-       A list of file paths to 3d NIFTI images.
-
-    Examples
-    --------
-    >>> os.chdir(tmpdir)
-    >>> b0_ixs = np.where(np.loadtxt(str(dipy_datadir / "HARDI193.bval")) <= 50)[0].tolist()[:2]
-    >>> in_file = str(dipy_datadir / "HARDI193.nii.gz")
-    >>> threeD_files = save_4d_to_3d(in_file)
-    >>> out_files = prune_b0s_from_dwis(threeD_files, b0_ixs)
-    >>> assert sum([os.path.isfile(i) for i in out_files]) == len(out_files)
-    >>> assert len(out_files) == len(threeD_files) - len(b0_ixs)
-    """
-    if in_files[0].endswith("_warped.nii.gz"):
-        out_files = [
-            i
-            for j, i in enumerate(
-                sorted(
-                    in_files, key=lambda x: int(x.split("_")[-2].split(".nii.gz")[0])
-                )
-            )
-            if j not in b0_ixs
-        ]
-    else:
-        out_files = [
-            i
-            for j, i in enumerate(
-                sorted(
-                    in_files, key=lambda x: int(x.split("_")[-1].split(".nii.gz")[0])
-                )
-            )
-            if j not in b0_ixs
-        ]
-    return out_files
-
-
 def save_3d_to_4d(in_files):
     """
     Concatenate a list of 3D volumes into a 4D output.

From 7767498cd0e9e8a9bb1cac07a2cdc1b083135a7e Mon Sep 17 00:00:00 2001
From: dPys <dpisner@utexas.edu>
Date: Tue, 24 Mar 2020 13:06:31 -0500
Subject: [PATCH 20/22] Implements summarize_images -- a function to merge the
 functionality of median_image and average_image using a callable argument

---
 .travis.yml                       |   4 ++--
 dmriprep/data/tests/dwi_b0.nii.gz | Bin 128 -> 122 bytes
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/.travis.yml b/.travis.yml
index 647dc893..e4babd02 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -81,9 +81,9 @@ before_script:
 script:
   - |
     if [ "$CHECK_TYPE" == "style" ]; then
-      flake8 dmriprep
+      travis_wait 60 flake8 dmriprep
     elif [ "$CHECK_TYPE" == "test" ]; then
-      pytest -n 2 -vv --doctest-modules --cov dmriprep --cov-config .coveragerc --cov-report xml:cov.xml dmriprep
+      travis_wait 60 pytest -n 2 -vv --doctest-modules --cov dmriprep --cov-config .coveragerc --cov-report xml:cov.xml dmriprep
     else
       false
     fi
diff --git a/dmriprep/data/tests/dwi_b0.nii.gz b/dmriprep/data/tests/dwi_b0.nii.gz
index b8a36ac29e01e3e6513b383457886fdc4d206cb2..134341f392f7ccf0449a5410519ae5f2692335c6 100644
GIT binary patch
delta 33
ocmZo*tYVYz=HO`Zp2#N8oRsunqJ|rfz;i|h&X50-G?*9|0H3G{ZU6uP

delta 39
ucmb<wU}Tf;=HO7YoX956ER+y1QNvA0Bj1F@&6z=&g<=1e)=3R43=9CUzX>k@


From 947a5cac7e211fd0d79eb3e3406c9010b9654211 Mon Sep 17 00:00:00 2001
From: dPys <dpisner@utexas.edu>
Date: Tue, 24 Mar 2020 13:57:41 -0500
Subject: [PATCH 21/22] Implements summarize_images -- a function to merge the
 functionality of median_image and average_image using a callable argument

---
 dmriprep/data/tests/dwi_b0.nii.gz | Bin 122 -> 122 bytes
 dmriprep/interfaces/images.py     |   7 +--
 dmriprep/utils/images.py          |  71 +++++++++++-------------------
 3 files changed, 30 insertions(+), 48 deletions(-)

diff --git a/dmriprep/data/tests/dwi_b0.nii.gz b/dmriprep/data/tests/dwi_b0.nii.gz
index 134341f392f7ccf0449a5410519ae5f2692335c6..cfe4ad75d988ca72cdffb63d9de6b2f4796dfcb3 100644
GIT binary patch
delta 12
Tcmb=bVw3OY;Ajt<$W{UX7efP`

delta 12
Tcmb=bVw3OY;Arxm$W{UX7XAZ`

diff --git a/dmriprep/interfaces/images.py b/dmriprep/interfaces/images.py
index 40a0b0a3..3d7a7e81 100644
--- a/dmriprep/interfaces/images.py
+++ b/dmriprep/interfaces/images.py
@@ -5,7 +5,7 @@
     traits, TraitedSpec, BaseInterfaceInputSpec, SimpleInterface, File,
     InputMultiObject, OutputMultiObject
 )
-from dmriprep.utils.images import extract_b0, median, rescale_b0, match_transforms
+from dmriprep.utils.images import extract_b0, summarize_images, rescale_b0, match_transforms
 
 LOGGER = logging.getLogger('nipype.interface')
 
@@ -80,6 +80,7 @@ class RescaleB0(SimpleInterface):
     output_spec = _RescaleB0OutputSpec
 
     def _run_interface(self, runtime):
+        import numpy as np
         from nipype.utils.filemanip import fname_presuffix
 
         out_b0s = fname_presuffix(
@@ -99,8 +100,8 @@ def _run_interface(self, runtime):
             self.inputs.in_file,
             self.inputs.mask_file, out_b0s
         )
-        self._results['out_ref'] = median(
-            self._results['out_b0s'],
+        self._results['out_ref'] = summarize_images(
+            self._results['out_b0s'], method=np.median,
             out_path=out_ref
         )
         return runtime
diff --git a/dmriprep/utils/images.py b/dmriprep/utils/images.py
index b2f297e6..ad981d3e 100644
--- a/dmriprep/utils/images.py
+++ b/dmriprep/utils/images.py
@@ -90,32 +90,48 @@ def rescale_b0(in_file, mask_file, out_path=None):
     return out_path
 
 
-def median(in_file, dtype=None, out_path=None):
+def summarize_images(in_file, method=np.median, dtype=None, out_path=None):
     """
-    Summarize a 4D dataset across the last dimension using median.
+    Summarize a 4D dataset across the last dimension using a
+    callable method.
 
     Parameters
     ----------
     in_file : str
-        A NIfTI file consisting of one or more B0's.
+        A NIfTI file consisting of one or more 3D images.
+    method : callable
+        A numpy function such as `np.mean` or `np.median`.
+    dtype : str
+        Optioally specify a datatype (e.g. 'float32').
     out_path : str
         Optionally specify an output path for `out_path`.
 
     Returns
     -------
     out_path : str
-       A 3D B0 NIFTI file.
+       A 3D NIFTI image file.
 
     Examples
     --------
     >>> os.chdir(tmpdir)
-    >>> in_file = str(data_dir / 'dwi_b0.nii.gz')
-    >>> out_path = median(in_file)
+    >>> in_file = str(dipy_datadir / "HARDI193.nii.gz")
+    >>> # Median case
+    >>> out_path = summarize_images(in_file)
+    >>> assert os.path.isfile(out_path)
+    >>> # Mean case
+    >>> out_path = summarize_images(in_file, method=np.mean)
     >>> assert os.path.isfile(out_path)
     """
+
+    if not callable(method):
+        raise ValueError('method must be callable')
+
+    # TODO: Check that callable is applicable (i.e. contains `axis` arg).
+    # if method.__call__()
     if out_path is None:
         out_path = fname_presuffix(
-            in_file, suffix='_b0ref', use_ext=True)
+            in_file, suffix="%s%s" % ('_', method.__code__.co_name),
+            use_ext=True)
 
     img = nb.load(in_file)
     if img.dataobj.ndim == 3:
@@ -124,8 +140,7 @@ def median(in_file, dtype=None, out_path=None):
         nb.squeeze_image(img).to_filename(out_path)
         return out_path
 
-    median_data = np.median(img.get_fdata(dtype=dtype),
-                            axis=-1)
+    summary_data = method(img.get_fdata(dtype=dtype), axis=-1)
 
     hdr = img.header.copy()
     hdr.set_xyzt_units('mm')
@@ -133,42 +148,8 @@ def median(in_file, dtype=None, out_path=None):
         hdr.set_data_dtype(dtype)
     else:
         dtype = hdr.get_data_dtype()
-    nb.Nifti1Image(median_data.astype(dtype), img.affine, hdr).to_filename(out_path)
-    return out_path
-
-
-def average_images(images, out_path=None):
-    """
-    Average a 4D dataset across the last dimension using mean.
-
-    Parameters
-    ----------
-    images : str
-        A list of NIFTI file paths.
-    out_path : str
-        Optionally specify an output path for `out_path`.
-
-    Returns
-    -------
-    out_path : str
-       A 3D NIFTI file averaged along the 4th dimension.
-
-    Examples
-    --------
-    >>> os.chdir(tmpdir)
-    >>> in_file = str(dipy_datadir / "HARDI193.nii.gz")
-    >>> out_files = save_4d_to_3d(in_file)
-    >>> out_path = average_images(out_files)
-    >>> assert os.path.isfile(out_path)
-    """
-    from nilearn.image import mean_img
-
-    average_img = mean_img([nb.load(img) for img in images])
-    if out_path is None:
-        out_path = fname_presuffix(
-            images[0], use_ext=False, suffix="_mean.nii.gz"
-        )
-    average_img.to_filename(out_path)
+    nb.Nifti1Image(summary_data.astype(dtype), img.affine,
+                   hdr).to_filename(out_path)
     return out_path
 
 

From 74036521c78cfe021da5e4783c6a906bc2e1e5f3 Mon Sep 17 00:00:00 2001
From: dPys <dpisner@utexas.edu>
Date: Tue, 24 Mar 2020 14:07:17 -0500
Subject: [PATCH 22/22] Add error-handling for save_3d_to_4d and save_4d_to_3d

---
 dmriprep/utils/images.py | 12 +++++++-----
 1 file changed, 7 insertions(+), 5 deletions(-)

diff --git a/dmriprep/utils/images.py b/dmriprep/utils/images.py
index ad981d3e..6632dd10 100644
--- a/dmriprep/utils/images.py
+++ b/dmriprep/utils/images.py
@@ -263,7 +263,11 @@ def save_4d_to_3d(in_file):
     >>> out_files = save_4d_to_3d(in_file)
     >>> assert len(out_files) == nb.load(in_file).shape[-1]
     """
-    files_3d = nb.four_to_three(nb.load(in_file))
+    filenii = nb.load(in_file)
+    if len(filenii.shape) != 4:
+        raise RuntimeError('Input image (%s) is not 4D.' % filenii)
+
+    files_3d = nb.four_to_three(filenii)
     out_files = []
     for i, file_3d in enumerate(files_3d):
         out_file = fname_presuffix(in_file, suffix="_tmp_{}".format(i))
@@ -300,10 +304,8 @@ def save_3d_to_4d(in_files):
     for i, f in enumerate(in_files):
         filenii = nb.load(f)
         filenii = nb.squeeze_image(filenii)
-        if len(filenii.shape) == 5:
-            raise RuntimeError('Input image (%s) is 5D.' % f)
-        if filenii.dataobj.ndim == 4:
-            nii_list += nb.four_to_three(filenii)
+        if len(filenii.shape) != 3:
+            raise RuntimeError('Input image (%s) is not 3D.' % f)
         else:
             nii_list.append(filenii)
     img_4d = nb.funcs.concat_images(nii_list)