DWI module

scilpy/dwi/utils.py

@@ -0,0 +1,70 @@
+# -*- coding: utf-8 -*-
+
+import numpy as np
+
+
+def _rescale_intensity(val, slope, in_max, bc_max):
+    """
+    Rescale an intensity value given a scaling factor.
+    This scaling factor ensures that the intensity
+    range before and after correction is the same.
+
+    Parameters
+    ----------
+    val: float
+         Value to be scaled
+    scale: float
+         Scaling factor to be applied
+    in_max: float
+         Max possible value
+    bc_max: float
+         Max value in the bias correction value range
+
+    Returns
+    -------
+    rescaled_value: float
+         Bias field corrected value scaled by the slope
+         of the data
+    """
+
+    return in_max - slope * (bc_max - val)
+
+
+# https://github.com/stnava/ANTs/blob/master/Examples/N4BiasFieldCorrection.cxx
+def rescale_dwi(in_data, bc_data):
+    """
+    Apply N4 Bias Field Correction to a DWI volume.
+    bc stands for bias correction. The code comes
+    from the C++ ANTS implmentation.
+
+    Parameters
+    ----------
+    in_data: ndarray (x, y, z, ndwi)
+         Input DWI volume 4-dimensional data.
+    bc_data: ndarray (x, y, z, ndwi)
+         Bias field correction volume estimated from ANTS
+         Copied for every dimension of the DWI 4-th dimension
+
+    Returns
+    -------
+    bc_data: ndarray (x, y, z, ndwi)
+         Bias field corrected DWI volume
+    """
+
+    in_min = np.amin(in_data)
+    in_max = np.amax(in_data)
+    bc_min = np.amin(bc_data)
+    bc_max = np.amax(bc_data)
+
+    slope = (in_max - in_min) / (bc_max - bc_min)
+
+    chunk = np.arange(0, len(in_data), 100000)
+    chunk = np.append(chunk, len(in_data))
+    for i in range(len(chunk)-1):
+        nz_bc_data = bc_data[chunk[i]:chunk[i+1]]
+        rescale_func = np.vectorize(_rescale_intensity, otypes=[np.float32])
+
+        rescaled_data = rescale_func(nz_bc_data, slope, in_max, bc_max)
+        bc_data[chunk[i]:chunk[i+1]] = rescaled_data
+
+    return bc_data

scripts/scil_apply_bias_field_on_dwi.py

@@ -16,6 +16,7 @@
 from scilpy.io.utils import (add_overwrite_arg,
                              assert_inputs_exist,
                              assert_outputs_exist)
+from scilpy.dwi.utils import rescale_dwi
 
 
 def _build_arg_parser():
@@ -35,31 +36,6 @@ def _build_arg_parser():
     return p
 
 
-def _rescale_intensity(val, slope, in_max, bc_max):
-    return in_max - slope * (bc_max - val)
-
-
-# https://github.com/stnava/ANTs/blob/master/Examples/N4BiasFieldCorrection.cxx
-def _rescale_dwi(in_data, bc_data):
-    in_min = np.amin(in_data)
-    in_max = np.amax(in_data)
-    bc_min = np.amin(bc_data)
-    bc_max = np.amax(bc_data)
-
-    slope = (in_max - in_min) / (bc_max - bc_min)
-
-    chunk = np.arange(0, len(in_data), 100000)
-    chunk = np.append(chunk, len(in_data))
-    for i in range(len(chunk)-1):
-        nz_bc_data = bc_data[chunk[i]:chunk[i+1]]
-        rescale_func = np.vectorize(_rescale_intensity, otypes=[np.float32])
-
-        rescaled_data = rescale_func(nz_bc_data, slope, in_max, bc_max)
-        bc_data[chunk[i]:chunk[i+1]] = rescaled_data
-
-    return bc_data
-
-
 def main():
     parser = _build_arg_parser()
     args = parser.parse_args()
@@ -79,11 +55,12 @@ def main():
     else:
         nz_mask_data = np.nonzero(np.average(dwi_data, axis=-1))
 
-    nuc_dwi_data = np.divide(dwi_data[nz_mask_data],
-                             bias_field_data[nz_mask_data].reshape((len(nz_mask_data[0]), 1)))
+    nuc_dwi_data = np.divide(
+        dwi_data[nz_mask_data],
+        bias_field_data[nz_mask_data].reshape((len(nz_mask_data[0]), 1)))
 
-    rescaled_nuc_data = _rescale_dwi(dwi_data[nz_mask_data],
-                                     nuc_dwi_data)
+    rescaled_nuc_data = rescale_dwi(dwi_data[nz_mask_data],
+                                    nuc_dwi_data)
 
     dwi_data[nz_mask_data] = rescaled_nuc_data
     nib.save(nib.Nifti1Image(dwi_data, dwi_img.affine,