Merge remote-tracking branch 'origin/pr/923'

Author: chrisgorgo

nipype/algorithms/misc.py

@@ -5,7 +5,7 @@
 
     Change directory to provide relative paths for doctests
     >>> import os
-    >>> filepath = os.path.dirname( os.path.realpath( __file__ ) )
+    >>> filepath = os.path.dirname(os.path.realpath(__file__))
     >>> datadir = os.path.realpath(os.path.join(filepath, '../testing/data'))
     >>> os.chdir(datadir)
 
@@ -34,7 +34,7 @@
 from ..interfaces.base import (BaseInterface, traits, TraitedSpec, File,
                                InputMultiPath, OutputMultiPath,
                                BaseInterfaceInputSpec, isdefined,
-                               DynamicTraitedSpec )
+                               DynamicTraitedSpec)
 from nipype.utils.filemanip import fname_presuffix, split_filename
 iflogger = logging.getLogger('interface')
 
@@ -785,7 +785,7 @@ def _list_outputs(self):
 
 class AddCSVRowInputSpec(DynamicTraitedSpec, BaseInterfaceInputSpec):
     in_file = traits.File(mandatory=True, desc='Input comma-separated value (CSV) files')
-    _outputs = traits.Dict( traits.Any, value={}, usedefault=True )
+    _outputs = traits.Dict(traits.Any, value={}, usedefault=True)
 
     def __setattr__(self, key, value):
         if key not in self.copyable_trait_names():
@@ -876,7 +876,7 @@ def _run_interface(self, runtime):
 
         if op.exists(self.inputs.in_file):
             formerdf = pd.read_csv(self.inputs.in_file, index_col=0)
-            df = pd.concat([formerdf, df], ignore_index=True )
+            df = pd.concat([formerdf, df], ignore_index=True)
 
         with open(self.inputs.in_file, 'w') as f:
             df.to_csv(f)
@@ -956,6 +956,115 @@ def calc_moments(timeseries_file, moment):
     return np.where(zero, 0, m3 / m2**(moment/2.0))
 
 
+class AddNoiseInputSpec(TraitedSpec):
+    in_file = File(exists=True, mandatory=True,
+                   desc='input image that will be corrupted with noise')
+    in_mask = File(exists=True, desc=('input mask, voxels outside this mask '
+                   'will be considered background'))
+    snr = traits.Float(10.0, desc='desired output SNR in dB', usedefault=True)
+    dist = traits.Enum('normal', 'rician', usedefault=True, mandatory=True,
+                       desc=('desired noise distribution'))
+    bg_dist = traits.Enum('normal', 'rayleigh', usedefault=True, mandatory=True,
+                          desc=('desired noise distribution, currently '
+                          'only normal is implemented'))
+    out_file = File(desc='desired output filename')
+
+
+class AddNoiseOutputSpec(TraitedSpec):
+    out_file = File(exists=True, desc='corrupted image')
+
+
+class AddNoise(BaseInterface):
+    """
+    Corrupts with noise the input image
+
+    Example
+    -------
+
+    >>> from nipype.algorithms.misc import AddNoise
+    >>> noise = AddNoise()
+    >>> noise.inputs.in_file = 'T1.nii'
+    >>> noise.inputs.in_mask = 'mask.nii'
+    >>> noise.snr = 30.0
+    >>> noise.run() # doctest: +SKIP
+    """
+    input_spec = AddNoiseInputSpec
+    output_spec = AddNoiseOutputSpec
+
+    def _run_interface(self, runtime):
+        in_image = nb.load(self.inputs.in_file)
+        in_data = in_image.get_data()
+        snr = self.inputs.snr
+
+        if isdefined(self.inputs.in_mask):
+            in_mask = nb.load(self.inputs.in_mask).get_data()
+        else:
+            in_mask = np.ones_like(in_data)
+
+        result = self.gen_noise(in_data, mask=in_mask, snr_db=snr,
+                                dist=self.inputs.dist, bg_dist=self.inputs.bg_dist)
+        res_im = nb.Nifti1Image(result, in_image.get_affine(), in_image.get_header())
+        res_im.to_filename(self._gen_output_filename())
+        return runtime
+
+    def _gen_output_filename(self):
+        if not isdefined(self.inputs.out_file):
+            _, base, ext = split_filename(self.inputs.in_file)
+            out_file = os.path.abspath('%s_SNR%03.2f%s' % (base, self.inputs.snr, ext))
+        else:
+            out_file = self.inputs.out_file
+
+        return out_file
+
+    def _list_outputs(self):
+        outputs = self.output_spec().get()
+        outputs['out_file'] = self._gen_output_filename()
+        return outputs
+
+    def gen_noise(self, image, mask=None, snr_db=10.0, dist='normal', bg_dist='normal'):
+        """
+        Generates a copy of an image with a certain amount of
+        added gaussian noise (rayleigh for background in mask)
+        """
+        from math import sqrt
+        snr = sqrt(np.power(10.0, snr_db/10.0))
+
+        if mask is None:
+            mask = np.ones_like(image)
+        else:
+            mask[mask>0] = 1
+            mask[mask<1] = 0
+
+            if mask.ndim < image.ndim:
+                mask = np.rollaxis(np.array([mask]*image.shape[3]),0, 4)
+
+        signal = image[mask>0].reshape(-1)
+
+        if dist == 'normal':
+            signal = signal - signal.mean()
+            sigma_n = sqrt(signal.var()/snr)
+            noise = np.random.normal(size=image.shape, scale=sigma_n)
+
+            if (np.any(mask==0)) and (bg_dist == 'rayleigh'):
+                bg_noise = np.random.rayleigh(size=image.shape, scale=sigma_n)
+                noise[mask==0] = bg_noise[mask==0]
+
+            im_noise = image +  noise
+
+        elif dist == 'rician':
+            sigma_n = signal.mean()/snr
+            n_1 = np.random.normal(size=image.shape, scale=sigma_n)
+            n_2 = np.random.normal(size=image.shape, scale=sigma_n)
+            stde_1 = n_1/sqrt(2.0)
+            stde_2 = n_2/sqrt(2.0)
+            im_noise = np.sqrt((image + stde_1)**2 + (stde_2)**2)
+        else:
+            raise NotImplementedError(('Only normal and rician distributions '
+                                      'are supported'))
+
+        return im_noise
+
+
 class NormalizeProbabilityMapSetInputSpec(TraitedSpec):
     in_files = InputMultiPath(File(exists=True, mandatory=True,
                     desc='The tpms to be normalized') )
@@ -973,9 +1082,6 @@ class NormalizeProbabilityMapSet(BaseInterface):
 
     .. note:: Please recall this is not a spatial normalization algorithm
 
-    Example
-    -------
-
     >>> import nipype.algorithms.misc as misc
     >>> normalize = misc.NormalizeProbabilityMapSet()
     >>> normalize.inputs.in_files = [ 'tpm_00.nii.gz', 'tpm_01.nii.gz', 'tpm_02.nii.gz' ]
@@ -988,10 +1094,10 @@ class NormalizeProbabilityMapSet(BaseInterface):
     def _run_interface(self, runtime):
         mask = None
 
-        if isdefined( self.inputs.in_mask ):
+        if isdefined(self.inputs.in_mask):
             mask = self.inputs.in_mask
 
-        self._out_filenames = normalize_tpms( self.inputs.in_files, mask )
+        self._out_filenames = normalize_tpms(self.inputs.in_files, mask)
         return runtime
 
     def _list_outputs(self):
@@ -1000,7 +1106,7 @@ def _list_outputs(self):
         return outputs
 
 
-def normalize_tpms( in_files, in_mask=None, out_files=[] ):
+def normalize_tpms(in_files, in_mask=None, out_files=[]):
     """
     Returns the input tissue probability maps (tpms, aka volume fractions)
     normalized to sum up 1.0 at each voxel within the mask.
@@ -1009,16 +1115,16 @@ def normalize_tpms( in_files, in_mask=None, out_files=[] ):
     import numpy as np
     import os.path as op
 
-    in_files = np.atleast_1d( in_files ).tolist()
+    in_files = np.atleast_1d(in_files).tolist()
 
-    if len(out_files)!=len(in_files):
-        for i,finname in enumerate( in_files ):
-            fname,fext = op.splitext( op.basename( finname ) )
+    if len(out_files) != len(in_files):
+        for i,finname in enumerate(in_files):
+            fname,fext = op.splitext(op.basename(finname))
             if fext == '.gz':
-                fname,fext2 = op.splitext( fname )
+                fname,fext2 = op.splitext(fname)
                 fext = fext2 + fext
 
-            out_file = op.abspath(fname+'_norm'+('_%02d' % i)+fext)
+            out_file = op.abspath('%s_norm_%02d%s' % (fname,i,fext))
             out_files+= [out_file]
 
     imgs = [nib.load(fim) for fim in in_files]
@@ -1028,39 +1134,39 @@ def normalize_tpms( in_files, in_mask=None, out_files=[] ):
         img_data[img_data>0.0] = 1.0
         hdr = imgs[0].get_header().copy()
         hdr['data_type']= 16
-        hdr.set_data_dtype( 'float32' )
-        nib.save( nib.Nifti1Image( img_data.astype(np.float32), imgs[0].get_affine(), hdr ), out_files[0] )
+        hdr.set_data_dtype(np.float32)
+        nib.save(nib.Nifti1Image(img_data.astype(np.float32), imgs[0].get_affine(), hdr), out_files[0])
         return out_files[0]
 
-    img_data = np.array( [ im.get_data() for im in imgs ] ).astype( 'f32' )
+    img_data = np.array([im.get_data() for im in imgs]).astype(np.float32)
     #img_data[img_data>1.0] = 1.0
     img_data[img_data<0.0] = 0.0
-    weights = np.sum( img_data, axis=0 )
+    weights = np.sum(img_data, axis=0)
 
-    msk = np.ones_like( imgs[0].get_data() )
+    msk = np.ones_like(imgs[0].get_data())
     msk[ weights<= 0 ] = 0
 
     if not in_mask is None:
-        msk = nib.load( in_mask ).get_data()
+        msk = nib.load(in_mask).get_data()
         msk[ msk<=0 ] = 0
         msk[ msk>0 ] = 1
 
-    msk = np.ma.masked_equal( msk, 0 )
+    msk = np.ma.masked_equal(msk, 0)
 
 
-    for i,out_file in enumerate( out_files ):
-        data = np.ma.masked_equal( img_data[i], 0 )
+    for i,out_file in enumerate(out_files):
+        data = np.ma.masked_equal(img_data[i], 0)
         probmap = data / weights
         hdr = imgs[i].get_header().copy()
         hdr['data_type']= 16
-        hdr.set_data_dtype( 'float32' )
-        nib.save( nib.Nifti1Image( probmap.astype(np.float32), imgs[i].get_affine(), hdr ), out_file )
+        hdr.set_data_dtype('float32')
+        nib.save(nib.Nifti1Image(probmap.astype(np.float32), imgs[i].get_affine(), hdr), out_file)
 
     return out_files
 
 
 # Deprecated interfaces ---------------------------------------------------------
-class Distance( nam.Distance ):
+class Distance(nam.Distance):
     """Calculates distance between two volumes.
 
     .. deprecated:: 0.10.0
@@ -1072,7 +1178,7 @@ def __init__(self, **inputs):
                       " please use nipype.algorithms.metrics.Distance"),
                       DeprecationWarning)
 
-class Overlap( nam.Overlap ):
+class Overlap(nam.Overlap):
     """Calculates various overlap measures between two maps.
 
     .. deprecated:: 0.10.0
@@ -1085,7 +1191,7 @@ def __init__(self, **inputs):
                       DeprecationWarning)
 
 
-class FuzzyOverlap( nam.FuzzyOverlap ):
+class FuzzyOverlap(nam.FuzzyOverlap):
     """Calculates various overlap measures between two maps, using a fuzzy
     definition.
 
@@ -1097,4 +1203,3 @@ def __init__(self, **inputs):
         warnings.warn(("This interface has been deprecated since 0.10.0,"
                       " please use nipype.algorithms.metrics.FuzzyOverlap"),
                       DeprecationWarning)
-

nipype/algorithms/tests/test_auto_AddNoise.py

@@ -0,0 +1,33 @@
+# AUTO-GENERATED by tools/checkspecs.py - DO NOT EDIT
+from nipype.testing import assert_equal
+from nipype.algorithms.misc import AddNoise
+
+def test_AddNoise_inputs():
+    input_map = dict(bg_dist=dict(mandatory=True,
+    usedefault=True,
+    ),
+    dist=dict(mandatory=True,
+    usedefault=True,
+    ),
+    in_file=dict(mandatory=True,
+    ),
+    in_mask=dict(),
+    out_file=dict(),
+    snr=dict(usedefault=True,
+    ),
+    )
+    inputs = AddNoise.input_spec()
+
+    for key, metadata in input_map.items():
+        for metakey, value in metadata.items():
+            yield assert_equal, getattr(inputs.traits()[key], metakey), value
+
+def test_AddNoise_outputs():
+    output_map = dict(out_file=dict(),
+    )
+    outputs = AddNoise.output_spec()
+
+    for key, metadata in output_map.items():
+        for metakey, value in metadata.items():
+            yield assert_equal, getattr(outputs.traits()[key], metakey), value
+