histogram equalization before hmc

Author: oesteban

nipype/workflows/dmri/preprocess/epi.py

@@ -229,7 +229,7 @@ def hmc_pipeline(name='motion_correct'):
     them to one reference image. Finally, the `b`-matrix is rotated accordingly [Leemans09]_
     making use of the rotation matrix obtained by FLIRT.
 
-    Search angles have been limited to 3.5 degrees, based on results in [Yendiki13]_.
+    Search angles have been limited to 4 degrees, based on results in [Yendiki13]_.
 
     A list of rigid transformation matrices is provided, so that transforms can be
     chained. This is useful to correct for artifacts with only one interpolation process (as
@@ -276,18 +276,24 @@ def hmc_pipeline(name='motion_correct'):
         outputnode.out_xfms - list of transformation matrices
 
     """
-    inputnode = pe.Node(niu.IdentityInterface(fields=['in_file', 'ref_num', 'in_bvec',
-                        'in_mask']), name='inputnode')
+    inputnode = pe.Node(niu.IdentityInterface(fields=['in_file', 'ref_num',
+                        'in_bvec', 'in_mask']), name='inputnode')
     split = pe.Node(fsl.Split(dimension='t'), name='SplitDWIs')
     pick_ref = pe.Node(niu.Select(), name='Pick_b0')
+    enhb0 = pe.Node(niu.Function(input_names=['in_file'],
+                    output_names=['out_file'], function=enhance),
+                    name='B0Equalize')
+    enhdw = pe.MapNode(niu.Function(input_names=['in_file'],
+                       output_names=['out_file'], function=enhance),
+                       name='DWEqualize', iterfield=['in_file'])
     flirt = pe.MapNode(fsl.FLIRT(interp='spline', cost='normmi',
-                       cost_func = 'normmi', dof=6, bins=64, save_log=True,
-                       searchr_x=[-4,4], searchr_y=[-4,4], searchr_z=[-4,4],
+                       cost_func='normmi', dof=6, bins=64, save_log=True,
+                       searchr_x=[-4, 4], searchr_y=[-4, 4], searchr_z=[-4, 4],
                        fine_search=1, coarse_search=10, padding_size=1),
                        name='CoRegistration', iterfield=['in_file'])
     rot_bvec = pe.Node(niu.Function(input_names=['in_bvec', 'in_matrix'],
-                           output_names=['out_file'], function=rotate_bvecs),
-                           name='Rotate_Bvec')
+                       output_names=['out_file'], function=rotate_bvecs),
+                       name='Rotate_Bvec')
     thres = pe.MapNode(fsl.Threshold(thresh=0.0), iterfield=['in_file'],
                        name='RemoveNegative')
     merge = pe.Node(fsl.Merge(dimension='t'), name='MergeDWIs')
@@ -301,10 +307,12 @@ def hmc_pipeline(name='motion_correct'):
         ,(split,      pick_ref,   [('out_files', 'inlist')])
         ,(inputnode,  pick_ref,   [(('ref_num', _checkrnum), 'index')])
         ,(inputnode,  flirt,      [('in_mask', 'ref_weight')])
-        ,(split,      flirt,      [('out_files', 'in_file')])
+        ,(pick_ref,   enhb0,      [('out', 'in_file')])
+        ,(split,      enhdw,      [('out_files', 'in_file')])
+        ,(enhb0,      flirt,      [('out_file', 'reference')])
+        ,(enhdw,      flirt,      [('out_file', 'in_file')])
         ,(inputnode,  rot_bvec,   [('in_bvec', 'in_bvec')])
         ,(flirt,      rot_bvec,   [('out_matrix_file', 'in_matrix')])
-        ,(pick_ref,   flirt,      [('out', 'reference')])
         ,(flirt,      thres,      [('out_file', 'in_file')])
         ,(thres,      merge,      [('out_file', 'in_files')])
         ,(merge,      outputnode, [('merged_file', 'out_file')])

nipype/workflows/dmri/preprocess/utils.py

@@ -4,12 +4,13 @@
 # vi: set ft=python sts=4 ts=4 sw=4 et:
 # @Author: oesteban
 # @Date:   2014-08-30 10:53:13
-# @Last Modified by:   Oscar Esteban
-# @Last Modified time: 2014-09-02 19:52:57
+# @Last Modified by:   oesteban
+# @Last Modified time: 2014-09-03 18:03:47
 import nipype.pipeline.engine as pe
 import nipype.interfaces.utility as niu
 from nipype.interfaces import fsl
 
+
 def cleanup_edge_pipeline(name='Cleanup'):
     """
     Perform some de-spiking filtering to clean up the edge of the fieldmap
@@ -20,8 +21,8 @@ def cleanup_edge_pipeline(name='Cleanup'):
     outputnode = pe.Node(niu.IdentityInterface(fields=['out_file']),
                          name='outputnode')
 
-    fugue = pe.Node(fsl.FUGUE(save_fmap=True, despike_2dfilter=True, despike_threshold=2.1),
-                    name='Despike')
+    fugue = pe.Node(fsl.FUGUE(save_fmap=True, despike_2dfilter=True,
+                    despike_threshold=2.1), name='Despike')
     erode = pe.Node(fsl.maths.MathsCommand(nan2zeros=True,
                     args='-kernel 2D -ero'), name='MskErode')
     newmsk = pe.Node(fsl.MultiImageMaths(op_string='-sub %s -thr 0.5 -bin'),
@@ -364,6 +365,7 @@ def siemens2rads(in_file, out_file=None):
     nb.Nifti1Image(data, im.get_affine(), hdr).to_filename(out_file)
     return out_file
 
+
 def rads2radsec(in_file, delta_te, out_file=None):
     """
     Converts input phase difference map to rads
@@ -483,3 +485,37 @@ def copy_hdr(in_file, in_file_hdr, out_file=None):
                          imref.get_affine(), imref.get_header())
     nii.to_filename(out_file)
     return out_file
+
+
+def enhance(in_file, clip_limit=0.015, in_mask=None, out_file=None):
+    import numpy as np
+    import nibabel as nb
+    import os.path as op
+    from skimage import exposure, img_as_int
+
+    if out_file is None:
+        fname, fext = op.splitext(op.basename(in_file))
+        if fext == '.gz':
+            fname, _ = op.splitext(fname)
+        out_file = op.abspath('./%s_enh.nii.gz' % fname)
+
+    im = nb.load(in_file)
+    imdata = im.get_data()
+    imshape = im.get_shape()
+
+    if in_mask is not None:
+        msk = nb.load(in_mask).get_data()
+        msk[msk > 0] = 1
+        msk[msk < 1] = 0
+        imdata = imdata * msk
+
+    immin = imdata.min()
+    imdata = (imdata - immin).astype(np.uint16)
+
+    adapted = exposure.equalize_adapthist(imdata.reshape(imshape[0], -1),
+                                          clip_limit=clip_limit)
+
+    nb.Nifti1Image(adapted.reshape(imshape), im.get_affine(),
+                   im.get_header()).to_filename(out_file)
+
+    return out_file