Merge branch 'pr/662'

Author: chrisgorgo

CHANGES

@@ -31,6 +31,9 @@ Next release
 * API: 'name' is now a positional argument for Workflow, Node, and MapNode constructors
 * API: SPM now defaults to SPM8 or SPM12b job format
 
+* ENH: New FSL interfaces: fsl.PrepareFieldmap, fsl.TOPUP, fsl.ApplyTOPUP, fsl.Eddy
+* ENH: New workflows: nipype.workflows.dmri.fsl.epi.[fieldmap_correction&topup_correction]
+
 Release 0.8.0 (May 8, 2013)
 ===========================
 

nipype/algorithms/misc.py

@@ -61,7 +61,6 @@ class PickAtlasOutputSpec(TraitedSpec):
 
 
 class PickAtlas(BaseInterface):
-
     '''
     Returns ROI masks given an atlas and a list of labels. Supports dilation
     and left right masking (assuming the atlas is properly aligned).
@@ -180,7 +179,6 @@ class ModifyAffineOutputSpec(TraitedSpec):
 
 
 class ModifyAffine(BaseInterface):
-
     '''
     Left multiplies the affine matrix with a specified values. Saves the volume
     as a nifti file.
@@ -245,7 +243,6 @@ class DistanceOutputSpec(TraitedSpec):
 
 
 class Distance(BaseInterface):
-
     '''
     Calculates distance between two volumes.
     '''
@@ -281,13 +278,7 @@ def _eucl_min(self, nii1, nii2):
         dist_matrix = cdist(set1_coordinates.T, set2_coordinates.T)
         (point1, point2) = np.unravel_index(
             np.argmin(dist_matrix), dist_matrix.shape)
-        return (
-            euclidean(
-                set1_coordinates.T[point1, :],
-                set2_coordinates.T[point2, :]
-            ),
-            set1_coordinates.T[point1, :], set2_coordinates.T[point2, :]
-        )
+        return (euclidean(set1_coordinates.T[point1, :], set2_coordinates.T[point2, :]), set1_coordinates.T[point1, :], set2_coordinates.T[point2, :])
 
     def _eucl_cog(self, nii1, nii2):
         origdata1 = nii1.get_data().astype(np.bool)
@@ -415,7 +406,6 @@ class OverlapOutputSpec(TraitedSpec):
 
 
 class Overlap(BaseInterface):
-
     """
     Calculates various overlap measures between two maps.
 
@@ -478,47 +468,25 @@ def _list_outputs(self):
 
 
 class FuzzyOverlapInputSpec(BaseInterfaceInputSpec):
-    in_ref = InputMultiPath(
-        File(exists=True), mandatory=True,
-        desc="Reference image. Requires the same dimensions as in_tst."
-    )
-    in_tst = InputMultiPath(
-        File(exists=True), mandatory=True,
-        desc="Test image. Requires the same dimensions as in_ref."
-    )
-    weighting = traits.Enum(
-        "none", "volume", "squared_vol",
-        desc='""none": no class-overlap weighting is performed\
-        "volume": computed class-overlaps are weighted by class volume\
-        "squared_vol": computed class-overlaps are weighted by the squared\
-        volume of the class', usedefault=True
-    )
-    out_file = File(
-        "diff.nii",
-        desc="alternative name for resulting difference-map",
-        usedefault=True
-    )
+    in_ref = InputMultiPath( File(exists=True), mandatory=True,
+                   desc="Reference image. Requires the same dimensions as in_tst.")
+    in_tst = InputMultiPath( File(exists=True), mandatory=True,
+                   desc="Test image. Requires the same dimensions as in_ref.")
+    weighting = traits.Enum("none", "volume", "squared_vol", desc='""none": no class-overlap weighting is performed\
+                            "volume": computed class-overlaps are weighted by class volume\
+                            "squared_vol": computed class-overlaps are weighted by the squared volume of the class',usedefault=True)
+    out_file = File("diff.nii", desc="alternative name for resulting difference-map", usedefault=True)
 
 
 class FuzzyOverlapOutputSpec(TraitedSpec):
-    jaccard = traits.Float(desc="Fuzzy Jaccard Index (fJI), all the classes")
-    dice = traits.Float(desc="Fuzzy Dice Index (fDI), all the classes")
-    diff_file = File(
-        exists=True, desc="resulting difference-map of all classes,\
-        using the chosen weighting"
-    )
-    class_fji = traits.List(
-        traits.Float(),
-        desc="Array containing the fJIs of each computed class"
-    )
-    class_fdi = traits.List(
-        traits.Float(),
-        desc="Array containing the fDIs of each computed class"
-    )
+    jaccard = traits.Float( desc="Fuzzy Jaccard Index (fJI), all the classes" )
+    dice = traits.Float( desc="Fuzzy Dice Index (fDI), all the classes" )
+    diff_file = File(exists=True, desc="resulting difference-map of all classes, using the chosen weighting" )
+    class_fji = traits.List( traits.Float(), desc="Array containing the fJIs of each computed class" )
+    class_fdi = traits.List( traits.Float(), desc="Array containing the fDIs of each computed class" )
 
 
 class FuzzyOverlap(BaseInterface):
-
     """
     Calculates various overlap measures between two maps, using the fuzzy
     definition proposed in: Crum et al., Generalized Overlap Measures for
@@ -539,77 +507,76 @@ class FuzzyOverlap(BaseInterface):
     >>> res = overlap.run() # doctest: +SKIP
     """
 
-    input_spec = FuzzyOverlapInputSpec
+    input_spec =  FuzzyOverlapInputSpec
     output_spec = FuzzyOverlapOutputSpec
 
     def _run_interface(self, runtime):
         ncomp = len(self.inputs.in_ref)
-        assert(ncomp == len(self.inputs.in_tst))
-        weights = np.ones(shape=ncomp)
+        assert( ncomp == len(self.inputs.in_tst) )
+        weights = np.ones( shape=ncomp )
+    
+        img_ref = np.array( [ nb.load( fname ).get_data() for fname in self.inputs.in_ref ] )
+        img_tst = np.array( [ nb.load( fname ).get_data() for fname in self.inputs.in_tst ] )
 
-        img_ref = np.array([nb.load(fname).get_data()
-                           for fname in self.inputs.in_ref])
-        img_tst = np.array([nb.load(fname).get_data()
-                           for fname in self.inputs.in_tst])
 
         msk = np.sum(img_ref, axis=0)
-        msk[msk > 0] = 1.0
+        msk[msk>0] = 1.0
         tst_msk = np.sum(img_tst, axis=0)
-        tst_msk[tst_msk > 0] = 1.0
+        tst_msk[tst_msk>0] = 1.0
+
+        #check that volumes are normalized
+        #img_ref[:][msk>0] = img_ref[:][msk>0] / (np.sum( img_ref, axis=0 ))[msk>0]
+        #img_tst[tst_msk>0] = img_tst[tst_msk>0] / np.sum( img_tst, axis=0 )[tst_msk>0]
 
         self._jaccards = []
         volumes = []
 
-        diff_im = np.zeros(img_ref.shape)
+        diff_im = np.zeros( img_ref.shape )
 
-        for ref_comp, tst_comp, diff_comp in zip(img_ref, img_tst, diff_im):
-            num = np.minimum(ref_comp, tst_comp)
-            ddr = np.maximum(ref_comp, tst_comp)
-            diff_comp[ddr > 0] += 1.0 - (num[ddr > 0] / ddr[ddr > 0])
-            self._jaccards.append(np.sum(num) / np.sum(ddr))
-            volumes.append(np.sum(ref_comp))
+        for ref_comp, tst_comp, diff_comp in zip( img_ref, img_tst, diff_im ):
+            num = np.minimum( ref_comp, tst_comp )
+            ddr = np.maximum( ref_comp, tst_comp )
+            diff_comp[ddr>0]+= 1.0-(num[ddr>0]/ddr[ddr>0])
+            self._jaccards.append( np.sum( num ) / np.sum( ddr ) )
+            volumes.append( np.sum( ref_comp ) )
 
-        self._dices = 2.0 * \
-            np.array(self._jaccards) / (np.array(self._jaccards) + 1.0)
+        self._dices = 2.0*np.array(self._jaccards) / (np.array(self._jaccards) +1.0 )
 
         if self.inputs.weighting != "none":
             weights = 1.0 / np.array(volumes)
             if self.inputs.weighting == "squared_vol":
-                weights = weights ** 2
+                weights = weights**2
 
-        weights = weights / np.sum(weights)
+        weights = weights / np.sum( weights )
 
-        setattr(self, '_jaccard',  np.sum(weights * self._jaccards))
-        setattr(self, '_dice', np.sum(weights * self._dices))
+        setattr( self, '_jaccard',  np.sum( weights * self._jaccards ) )
+        setattr( self, '_dice', np.sum( weights * self._dices ) )
 
-        diff = np.zeros(diff_im[0].shape)
 
-        for w, ch in zip(weights, diff_im):
-            ch[msk == 0] = 0
-            diff += w * ch
+        diff = np.zeros( diff_im[0].shape )
 
-        nb.save(
-            nb.Nifti1Image(
-                diff,
-                nb.load(self.inputs.in_ref[0]).get_affine(),
-                nb.load(self.inputs.in_ref[0]).get_header()
-            ),
-            self.inputs.out_file
-        )
+        for w,ch in zip(weights,diff_im):
+            ch[msk==0] = 0
+            diff+= w* ch
+        
+        nb.save(nb.Nifti1Image(diff, nb.load( self.inputs.in_ref[0]).get_affine(),
+                nb.load( self.inputs.in_ref[0]).get_header()), self.inputs.out_file )
 
+ 
         return runtime
 
     def _list_outputs(self):
         outputs = self._outputs().get()
         for method in ("dice", "jaccard"):
             outputs[method] = getattr(self, '_' + method)
+        #outputs['volume_difference'] = self._volume
         outputs['diff_file'] = os.path.abspath(self.inputs.out_file)
-        outputs['class_fji'] = np.array(
-            self._jaccards).astype(float).tolist()
-        outputs['class_fdi'] = self._dices.astype(float).tolist()
+        outputs['class_fji'] =  np.array(self._jaccards).astype(float).tolist();
+        outputs['class_fdi']=  self._dices.astype(float).tolist();
         return outputs
 
 
+
 class CreateNiftiInputSpec(BaseInterfaceInputSpec):
     data_file = File(exists=True, mandatory=True, desc="ANALYZE img file")
     header_file = File(
@@ -664,7 +631,6 @@ class TSNROutputSpec(TraitedSpec):
 
 
 class TSNR(BaseInterface):
-
     """Computes the time-course SNR for a time series
 
     Typically you want to run this on a realigned time-series.
@@ -742,7 +708,6 @@ class GunzipOutputSpec(TraitedSpec):
 
 
 class Gunzip(BaseInterface):
-
     """
 
     """
@@ -784,7 +749,7 @@ def matlab2csv(in_array, name, reshape):
     if reshape:
         if len(np.shape(output_array)) > 1:
             output_array = np.reshape(output_array, (
-                np.shape(output_array)[0] * np.shape(output_array)[1], 1))
+                np.shape(output_array)[0]*np.shape(output_array)[1], 1))
             iflogger.info(np.shape(output_array))
     output_name = op.abspath(name + '.csv')
     np.savetxt(output_name, output_array, delimiter=',')
@@ -808,7 +773,6 @@ class Matlab2CSVOutputSpec(TraitedSpec):
 
 
 class Matlab2CSV(BaseInterface):
-
     """
     Simple interface to save the components of a MATLAB .mat file as a text
     file with comma-separated values (CSVs).
@@ -841,10 +805,7 @@ def _run_interface(self, runtime):
                 if isinstance(in_dict[key][0], np.ndarray):
                     saved_variables.append(key)
                 else:
-                    iflogger.info(
-                        'One of the keys in the input file, {k},\
-                        is not a Numpy array'.format(k=key)
-                    )
+                    iflogger.info('One of the keys in the input file, {k}, is not a Numpy array'.format(k=key))
 
         if len(saved_variables) > 1:
             iflogger.info(
@@ -910,9 +871,7 @@ def merge_csvs(in_list):
                     )
                 except ValueError, ex:
                     in_array = np.loadtxt(
-                        in_file, delimiter=',', skiprows=1,
-                        usecols=range(1, n_cols - 1)
-                    )
+                        in_file, delimiter=',', skiprows=1, usecols=range(1, n_cols-1))
         if idx == 0:
             out_array = in_array
         else:
@@ -930,7 +889,7 @@ def remove_identical_paths(in_files):
         out_names = list()
         commonprefix = op.commonprefix(in_files)
         lastslash = commonprefix.rfind('/')
-        commonpath = commonprefix[0:(lastslash + 1)]
+        commonpath = commonprefix[0:(lastslash+1)]
         for fileidx, in_file in enumerate(in_files):
             path, name, ext = split_filename(in_file)
             in_file = op.join(path, name)
@@ -948,10 +907,10 @@ def maketypelist(rowheadings, shape, extraheadingBool, extraheading):
     if rowheadings:
         typelist.append(('heading', 'a40'))
     if len(shape) > 1:
-        for idx in range(1, (min(shape) + 1)):
+        for idx in range(1, (min(shape)+1)):
             typelist.append((str(idx), float))
     else:
-        for idx in range(1, (shape[0] + 1)):
+        for idx in range(1, (shape[0]+1)):
             typelist.append((str(idx), float))
     if extraheadingBool:
         typelist.append((extraheading, 'a40'))
@@ -966,13 +925,13 @@ def makefmtlist(output_array, typelist, rowheadingsBool,
         fmtlist.append('%s')
     if len(shape) > 1:
         output = np.zeros(max(shape), typelist)
-        for idx in range(1, min(shape) + 1):
-            output[str(idx)] = output_array[:, idx - 1]
+        for idx in range(1, min(shape)+1):
+            output[str(idx)] = output_array[:, idx-1]
             fmtlist.append('%f')
     else:
         output = np.zeros(1, typelist)
-        for idx in range(1, len(output_array) + 1):
-            output[str(idx)] = output_array[idx - 1]
+        for idx in range(1, len(output_array)+1):
+            output[str(idx)] = output_array[idx-1]
             fmtlist.append('%f')
     if extraheadingBool:
         fmtlist.append('%s')
@@ -1007,7 +966,6 @@ class MergeCSVFilesOutputSpec(TraitedSpec):
 
 
 class MergeCSVFiles(BaseInterface):
-
     """
     This interface is designed to facilitate data loading in the R environment.
     It takes input CSV files and merges them into a single CSV file.
@@ -1144,7 +1102,6 @@ class AddCSVColumnOutputSpec(TraitedSpec):
 
 
 class AddCSVColumn(BaseInterface):
-
     """
     Short interface to add an extra column and field to a text file
 
@@ -1206,7 +1163,6 @@ class CalculateNormalizedMomentsOutputSpec(TraitedSpec):
 
 
 class CalculateNormalizedMoments(BaseInterface):
-
     """
     Calculates moments of timeseries.
 
@@ -1248,4 +1204,4 @@ def calc_moments(timeseries_file, moment):
     m2 = stats.moment(timeseries, 2, axis=0)
     m3 = stats.moment(timeseries, moment, axis=0)
     zero = (m2 == 0)
-    return np.where(zero, 0, m3 / m2 ** (moment / 2.0))
+    return np.where(zero, 0, m3 / m2**(moment/2.0))

nipype/interfaces/fsl/__init__.py

@@ -15,9 +15,13 @@
 from .utils import (Smooth, Merge, ExtractROI, Split, ImageMaths, ImageMeants,
                     ImageStats, FilterRegressor, Overlay, Slicer,
                     PlotTimeSeries, PlotMotionParams, ConvertXFM,
-                    SwapDimensions, PowerSpectrum, SigLoss, Reorient2Std,
+                    SwapDimensions, PowerSpectrum, Reorient2Std,
                     Complex, InvWarp)
-from .dti import (EddyCorrect, BEDPOSTX, DTIFit, ProbTrackX, VecReg, ProjThresh,
+
+from .epi import (PrepareFieldmap, TOPUP, ApplyTOPUP, Eddy, EPIDeWarp, 
+                  SigLoss, EddyCorrect)
+
+from .dti import (BEDPOSTX, DTIFit, ProbTrackX, VecReg, ProjThresh,
                   FindTheBiggest, DistanceMap, TractSkeleton, XFibres,
                   MakeDyadicVectors)
 from .maths import (ChangeDataType, Threshold, MeanImage, ApplyMask,