Added Camino shredder and MESD functions in interface convert.py and odf.py respectively

mick-d · mick-d · commit 73b7943931d6 · 2014-03-02T20:49:54.000-05:00
diff --git a/nipype/interfaces/camino/__init__.py b/nipype/interfaces/camino/__init__.py
@@ -5,10 +5,11 @@
 
 from .connectivity import Conmat
 from .convert import (Image2Voxel, FSL2Scheme, VtkStreamlines, ProcStreamlines,
-                      TractShredder, DT2NIfTI, NIfTIDT2Camino, AnalyzeHeader)
+                      TractShredder, DT2NIfTI, NIfTIDT2Camino, AnalyzeHeader,
+                      Shredder)
 from .dti import (DTIFit, ModelFit, DTLUTGen, PicoPDFs, Track, TrackPICo,
                   TrackBayesDirac, TrackDT, TrackBallStick, TrackBootstrap,
                   ComputeFractionalAnisotropy, ComputeMeanDiffusivity,
                   ComputeTensorTrace,  ComputeEigensystem, DTMetric)
 from .calib import (SFPICOCalibData, SFLUTGen)
-from .odf import (QBallMX, LinRecon, SFPeaks)
+from .odf import (QBallMX, LinRecon, SFPeaks, MESD)
diff --git a/nipype/interfaces/camino/convert.py b/nipype/interfaces/camino/convert.py
@@ -581,3 +581,54 @@ def _list_outputs(self):
     def _gen_outfilename(self):
         _, name , _ = split_filename(self.inputs.in_file)
         return name + ".hdr"
+
+class ShredderInputSpec(StdOutCommandLineInputSpec):
+    in_file = File(exists=True, argstr='< %s', mandatory=True, position=-2, desc='raw binary data file')
+
+    offset = traits.Int(argstr='%d', units='NA',
+        desc='initial offset of offset bytes', position=1)
+
+    chunksize = traits.Int(argstr='%d', units='NA',
+        desc='reads and outputs a chunk of chunksize bytes', position=2)
+
+    space = traits.Int(argstr='%d', units='NA',
+        desc='skips space bytes', position=3)
+
+class ShredderOutputSpec(TraitedSpec):
+    shredded = File(exists=True, desc='Shredded binary data file')
+
+class Shredder(StdOutCommandLine):
+    """
+    Extracts periodic chunks from a data stream.
+
+    Shredder makes an initial offset of offset bytes. It then reads and outputs
+    chunksize bytes, skips space bytes, and repeats until there is no more input.
+
+    If  the  chunksize  is  negative, chunks of size |chunksize| are read and the
+    byte ordering of each chunk is reversed. The whole chunk will be reversed, so
+    the chunk must be the same size as the data type, otherwise the order of the
+    values in the chunk, as well as their endianness, will be reversed.
+
+    Examples
+    --------
+
+    >>> import nipype.interfaces.camino as cam
+    >>> shred = cam.Shredder()
+    >>> shred.inputs.in_file = 'SubjectA.Bfloat'
+    >>> shred.inputs.offset = 0
+    >>> shred.inputs.chunksize = 1
+    >>> shred.inputs.space = 2
+    >>> shred.run()                  # doctest: +SKIP
+    """
+    _cmd = 'shredder'
+    input_spec=ShredderInputSpec
+    output_spec=ShredderOutputSpec
+
+    def _list_outputs(self):
+        outputs = self.output_spec().get()
+        outputs['shredded_file'] = os.path.abspath(self._gen_outfilename())
+        return outputs
+
+    def _gen_outfilename(self):
+        _, name , _ = split_filename(self.inputs.in_file)
+        return name + "_shredded"
diff --git a/nipype/interfaces/camino/odf.py b/nipype/interfaces/camino/odf.py
@@ -160,6 +160,126 @@ def _gen_outfilename(self):
         _, name , _ = split_filename(self.inputs.scheme_file)
         return name + '_recondata.Bdouble'
 
+class MESDInputSpec(StdOutCommandLineInputSpec):
+    in_file = File(exists=True, argstr='-inputfile %s', mandatory=True, position=1,
+                   desc='voxel-order data filename')
+    inverter = traits.Enum('SPIKE', 'PAS', argstr='-filter %s', position=2, mandatory=True,
+                           desc=('The inversion index specifies the type of inversion to perform on the data.'
+                                 'The currently available choices are:'
+                                  'Inverter name  | Inverter parameters'
+                                  '---------------|------------------'
+                                  'SPIKE          | bd (b-value x diffusivity along the fibre.)'
+                                  'PAS            | r'))
+    inverter_param = traits.Float(argstr='%f', units='NA', position=3, mandatory=True,
+                                  desc=('Parameter associated with the inverter. Cf. inverter description for'
+                                        'more information.'))
+    fastmesd = traits.Bool(argstr='-fastmesd', requires=['mepointset'],
+                           desc=('Turns off numerical integration checks and fixes the integration point set size at that of'
+                                 'the index specified by -basepointset..'))
+    mepointset = traits.Int(argstr='-mepointset %d', units='NA',
+                            desc=('Use a set of directions other than those in the scheme file for the deconvolution kernel.'
+                                  'The number refers to the number of directions on the unit sphere. For example, '
+                                  '"-mepointset 54" uses the directions in "camino/PointSets/Elec054.txt".'))
+    scheme_file = File(exists=True, argstr='-schemefile %s', mandatory=True,
+                       desc='Specifies the scheme file for the diffusion MRI data')
+    bgmask = File(exists=True, argstr='-bgmask %s', desc='background mask')
+    inputdatatype = traits.Enum('float', 'char', 'short', 'int', 'long', 'double', argstr='-inputdatatype %s',
+                                desc=('Specifies the data type of the input file: "char", "short", "int", "long",'
+                                     '"float" or "double". The input file must have BIG-ENDIAN ordering.'
+                                     'By default, the input type is "float".'))
+
+class MESDOutputSpec(TraitedSpec):
+    mesd_data = File(exists=True, desc='MESD data')
+
+class MESD(StdOutCommandLine):
+    """
+    MESD is a general program for maximum entropy spherical deconvolution.
+    It also runs PASMRI, which is a special case of spherical deconvolution.
+    The input data must be in voxel order.
+
+    The format of the output in each voxel is:
+    { exitcode, ln(A^star(0)), lambda_0, lambda_1, ..., lambda_N }
+
+    The  exitcode  contains  the  results of three tests. The first test thresholds
+    the maximum relative error between the numerical integrals computed at con-
+    vergence and those computed using a larger test point set; if the error is
+    greater than a threshold the exitcode is increased from zero to one as a
+    warning; if it is greater than a larger threshold the exitcode is increased to
+    two to suggest failure. The  second  test  thresholds  the predicted  error in
+    numerical integrals computed using the test point set; if the predicted error
+    is greater than a threshold the exitcode is increased by 10. The third test
+    thresholds the RMS error between the measurements and their predictions from
+    the fitted deconvolution; if the errors are greater than a threshold, the exit
+    code is increased by 100. An exitcode of 112 means that all three tests were
+    failed and the result is likely to be unreliable.  If all is well the exitcode
+    is zero. Results are often still reliable even if one or two of the tests are
+    failed.
+
+    Other possible exitcodes are:
+     5   - The optimization failed to converge
+    -1   - Background
+    -100 - Something wrong in the MRI data, e.g. negative or zero measurements,
+           so that the optimization could not run.
+
+    The  standard  MESD  implementation  is computationally demanding, particularly
+    as the number of measurements increases (computation is approximately O(N^2),
+    where N is the number of measurements). There are two ways to obtain significant
+    computational speed-up:
+
+    i) Turn off error checks and use a small point set for computing numerical
+    integrals in the algorithm by adding the flag -fastmesd. Sakaie CDMRI 2008
+    shows that using the smallest point set  (-basepointset  0)  with  no
+    error  checks  usually  has only a minor effect on the output of the algorithm,
+    but provides a major reduction in computation time. You can increase the point
+    set size using -basepointset with an argument higher than 0, which may produce
+    better results in some voxels, but will increase computation time, which
+    approximately doubles every time the point set index increases by 1.
+
+    ii) Reduce the complexity of the maximum entropy encoding using -mepointset <X>.
+    By default <X> = N, the number of measurements, and is the number of parameters
+    in the max.  ent. representation of the  output  function, ie  the  number of
+    lambda parameters, as described in Jansons and Alexander Inverse Problems 2003.
+    However, we can represent the function using less components and <X> here
+    specifies the number of lambda parameters. To obtain speed-up, set <X>
+    < N; complexity become O(<X>^2) rather than O(N^2). Note that <X> must be chosen
+    so that the camino/PointSets directory contains a point set  with  that  number
+    of  elements.  When  -mepointset decreases, the  numerical  integration checks
+    make less and less of a difference and smaller point sets for numerical
+    integration (see -basepointset) become adequate. So when <X> is low -fastmesd is
+    worth using to get even more speed-up.
+
+    The choice of <X> is a parameter of the technique. Too low and you lose angular
+    resoloution; too high and you see no computational benefit and may even suffer
+    from overfitting. Empirically, we  have  found  that  <X>=16 often  gives  good
+    results and good speed up, but it is worth trying a few values a comparing
+    performance. The reduced encoding is described in the following ISMRM abstract:
+    Sweet and Alexander "Reduced Encoding Persistent Angular Structure" 572 ISMRM 2010.
+
+    Example
+    ---------
+    Run MESD on every voxel of the data file SubjectA.Bfloat using the PASMRI kernel.
+
+    >>> import nipype.interfaces.camino as cam
+    >>> mesd = cam.MESD()
+    >>> mesd.inputs.in_file = 'SubjectA.Bfloat'
+    >>> mesd.inputs.scheme_file = 'A.scheme'
+    >>> mesd.inputs.inverter = 'PAS'
+    >>> mesd.inputs.inverter_param = 1.4
+    >>> mesd.run()            # doctest: +SKIP
+    """
+    _cmd = 'mesd'
+    input_spec=MESDInputSpec
+    output_spec=MESDOutputSpec
+
+    def _list_outputs(self):
+        outputs = self.output_spec().get()
+        outputs['recon_data'] = os.path.abspath(self._gen_outfilename())
+        return outputs
+
+    def _gen_outfilename(self):
+        _, name , _ = split_filename(self.inputs.scheme_file)
+        return name + '_MESD.Bdouble'
+
 class SFPeaksInputSpec(StdOutCommandLineInputSpec):
     in_file = File(exists=True, argstr='-inputfile %s', mandatory=True,
                    desc='Voxel-order data of spherical functions')
