Merge branch 'master' into fix/AreaWeighedP2PDistance

Conflicts:
	CHANGES

Author: oesteban

CHANGES

@@ -6,6 +6,7 @@ Next release
 * FIX: Refactor P2PDistance, change name to ComputeMeshWarp, add regression tests,
        fix bug in area weighted distance, and added optimizations
        (https://github.com/nipy/nipype/pull/1016)
+* ENH: Add the average distance to ErrorMap (https://github.com/nipy/nipype/pull/1039)
 * ENH: Inputs with name_source can be now chained in cascade (https://github.com/nipy/nipype/pull/938)
 * ENH: Improve JSON interfaces: default settings when reading and consistent output creation
        when writing (https://github.com/nipy/nipype/pull/1047)

examples/rsfmri_vol_surface_preprocessing.py

@@ -230,7 +230,7 @@ def build_filter1(motion_params, comp_norm, outliers, detrend_poly=None):
             out_params = np.hstack((out_params, outlier_vector))
         if detrend_poly:
             timepoints = out_params.shape[0]
-            X = np.ones((timepoints, 1))
+            X = np.empty((timepoints, 0))
             for i in range(detrend_poly):
                 X = np.hstack((X, legendre(
                     i + 1)(np.linspace(-1, 1, timepoints))[:, None]))
@@ -542,7 +542,7 @@ def create_reg_workflow(name='registration'):
 
     warpmean = Node(ants.ApplyTransforms(), name='warpmean')
     warpmean.inputs.input_image_type = 3
-    warpmean.inputs.interpolation = 'BSpline'
+    warpmean.inputs.interpolation = 'Linear'
     warpmean.inputs.invert_transform_flags = [False, False]
     warpmean.inputs.terminal_output = 'file'
     warpmean.inputs.args = '--float'
@@ -771,7 +771,7 @@ def merge_files(in1, in2):
     warpall = MapNode(ants.ApplyTransforms(), iterfield=['input_image'],
                       name='warpall')
     warpall.inputs.input_image_type = 3
-    warpall.inputs.interpolation = 'BSpline'
+    warpall.inputs.interpolation = 'Linear'
     warpall.inputs.invert_transform_flags = [False, False]
     warpall.inputs.terminal_output = 'file'
     warpall.inputs.reference_image = target_file
@@ -1003,9 +1003,11 @@ def create_resting_workflow(args, name=None):
     parser.add_argument('--surf_fwhm', default=15., dest='surf_fwhm',
                         type=float, help="Spatial FWHM" + defstr)
     parser.add_argument("-l", "--lowpass_freq", dest="lowpass_freq",
-                        default=0.1, help="Low pass frequency (Hz)" + defstr)
+                        default=0.1, type=float,
+                        help="Low pass frequency (Hz)" + defstr)
     parser.add_argument("-u", "--highpass_freq", dest="highpass_freq",
-                        default=0.01, help="High pass frequency (Hz)" + defstr)
+                        default=0.01, type=float,
+                        help="High pass frequency (Hz)" + defstr)
     parser.add_argument("-o", "--output_dir", dest="sink",
                         help="Output directory base", required=True)
     parser.add_argument("-w", "--work_dir", dest="work_dir",

examples/rsfmri_vol_surface_preprocessing_nipy.py

@@ -217,7 +217,7 @@ def build_filter1(motion_params, comp_norm, outliers, detrend_poly=None):
             out_params = np.hstack((out_params, outlier_vector))
         if detrend_poly:
             timepoints = out_params.shape[0]
-            X = np.ones((timepoints, 1))
+            X = np.empty((timepoints, 0))
             for i in range(detrend_poly):
                 X = np.hstack((X, legendre(
                     i + 1)(np.linspace(-1, 1, timepoints))[:, None]))
@@ -526,7 +526,7 @@ def create_reg_workflow(name='registration'):
     """
     warpmean = Node(ants.ApplyTransforms(), name='warpmean')
     warpmean.inputs.input_image_type = 3
-    warpmean.inputs.interpolation = 'BSpline'
+    warpmean.inputs.interpolation = 'Linear'
     warpmean.inputs.invert_transform_flags = [False, False]
     warpmean.inputs.terminal_output = 'file'
     warpmean.inputs.args = '--float'
@@ -749,7 +749,7 @@ def merge_files(in1, in2):
     warpall = MapNode(ants.ApplyTransforms(), iterfield=['input_image'],
                       name='warpall')
     warpall.inputs.input_image_type = 3
-    warpall.inputs.interpolation = 'BSpline'
+    warpall.inputs.interpolation = 'Linear'
     warpall.inputs.invert_transform_flags = [False, False]
     warpall.inputs.terminal_output = 'file'
     warpall.inputs.reference_image = target_file
@@ -996,9 +996,11 @@ def create_resting_workflow(args, name=None):
     parser.add_argument('--surf_fwhm', default=15., dest='surf_fwhm',
                         type=float, help="Spatial FWHM" + defstr)
     parser.add_argument("-l", "--lowpass_freq", dest="lowpass_freq",
-                        default=0.1, help="Low pass frequency (Hz)" + defstr)
+                        default=0.1, type=float,
+                        help="Low pass frequency (Hz)" + defstr)
     parser.add_argument("-u", "--highpass_freq", dest="highpass_freq",
-                        default=0.01, help="High pass frequency (Hz)" + defstr)
+                        default=0.01, type=float,
+                        help="High pass frequency (Hz)" + defstr)
     parser.add_argument("-o", "--output_dir", dest="sink",
                         help="Output directory base", required=True)
     parser.add_argument("-w", "--work_dir", dest="work_dir",

nipype/algorithms/metrics.py

@@ -488,6 +488,7 @@ class ErrorMapInputSpec(BaseInterfaceInputSpec):
 
 class ErrorMapOutputSpec(TraitedSpec):
     out_map = File(exists=True, desc="resulting error map")
+    distance = traits.Float(desc="Average distance between volume 1 and 2")
 
 
 class ErrorMap(BaseInterface):
@@ -506,12 +507,13 @@ class ErrorMap(BaseInterface):
     _out_file = ''
 
     def _run_interface(self, runtime):
-        from scipy.spatial.distance import cdist, pdist
+        # Get two numpy data matrices
         nii_ref = nb.load(self.inputs.in_ref)
         ref_data = np.squeeze(nii_ref.get_data())
         tst_data = np.squeeze(nb.load(self.inputs.in_tst).get_data())
         assert(ref_data.ndim == tst_data.ndim)
 
+        # Load mask
         comps = 1
         mapshape = ref_data.shape
 
@@ -528,26 +530,29 @@ def _run_interface(self, runtime):
         else:
             msk = np.ones(shape=mapshape)
 
+        # Flatten both volumes and make the pixel differennce
         mskvector = msk.reshape(-1)
         msk_idxs = np.where(mskvector==1)
         refvector = ref_data.reshape(-1,comps)[msk_idxs].astype(np.float32)
         tstvector = tst_data.reshape(-1,comps)[msk_idxs].astype(np.float32)
         diffvector = (refvector-tstvector)
 
+        # Scale the difference
         if self.inputs.metric == 'sqeuclidean':
             errvector = diffvector**2
+            if (comps > 1):
+                errvector = np.sum(errvector, axis=1)
+            else:
+                errvector = np.squeeze(errvector)
         elif self.inputs.metric == 'euclidean':
-            X = np.hstack((refvector, tstvector))
-            errvector = np.linalg.norm(X, axis=1)
+            errvector = np.linalg.norm(diffvector, axis=1)
 
-        if (comps > 1):
-            errvector = np.sum(errvector, axis=1)
-        else:
-            errvector = np.squeeze(errvector)
-
-        errvectorexp = np.zeros_like(mskvector)
+        errvectorexp = np.zeros_like(mskvector, dtype=np.float32) # The default type is uint8
         errvectorexp[msk_idxs] = errvector
 
+        # Get averaged error
+        self._distance = np.average(errvector) # Only average the masked voxels
+
         errmap = errvectorexp.reshape(mapshape)
 
         hdr = nii_ref.get_header().copy()
@@ -567,11 +572,12 @@ def _run_interface(self, runtime):
         nb.Nifti1Image(errmap.astype(np.float32), nii_ref.get_affine(),
                        hdr).to_filename(self._out_file)
 
-        return runtime
+        return runtime 
 
     def _list_outputs(self):
         outputs = self.output_spec().get()
         outputs['out_map'] = self._out_file
+        outputs['distance'] = self._distance
         return outputs
 
 

nipype/algorithms/tests/test_errormap.py

@@ -0,0 +1,76 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+
+from nipype.testing import (assert_equal, example_data)
+from nipype.algorithms.metrics import ErrorMap
+import nibabel as nib
+import numpy as np
+from tempfile import mkdtemp
+import os
+
+def test_errormap():
+
+	tempdir = mkdtemp()
+	# Single-Spectual
+	# Make two fake 2*2*2 voxel volumes
+	volume1 = np.array([[[2.0, 8.0], [1.0, 2.0]], [[1.0, 9.0], [0.0, 3.0]]]) # John von Neumann's birthday
+	volume2 = np.array([[[0.0, 7.0], [2.0, 3.0]], [[1.0, 9.0], [1.0, 2.0]]]) # Alan Turing's birthday
+	mask = np.array([[[1, 0], [0, 1]], [[1, 0], [0, 1]]])
+
+	img1 = nib.Nifti1Image(volume1, np.eye(4))
+	img2 = nib.Nifti1Image(volume2, np.eye(4))
+	maskimg = nib.Nifti1Image(mask, np.eye(4))
+
+	nib.save(img1, os.path.join(tempdir, 'von.nii.gz'))
+	nib.save(img2, os.path.join(tempdir, 'alan.nii.gz'))
+	nib.save(maskimg, os.path.join(tempdir, 'mask.nii.gz'))
+
+
+	# Default metric
+	errmap = ErrorMap()
+	errmap.inputs.in_tst = os.path.join(tempdir, 'von.nii.gz')
+	errmap.inputs.in_ref = os.path.join(tempdir, 'alan.nii.gz')
+	errmap.out_map = os.path.join(tempdir, 'out_map.nii.gz')
+	result = errmap.run()
+	yield assert_equal, result.outputs.distance, 1.125 
+
+	# Square metric
+	errmap.inputs.metric = 'sqeuclidean'
+	result = errmap.run()
+	yield assert_equal, result.outputs.distance, 1.125
+
+	# Linear metric
+	errmap.inputs.metric = 'euclidean'
+	result = errmap.run()
+	yield assert_equal, result.outputs.distance, 0.875 
+
+	# Masked 
+	errmap.inputs.mask = os.path.join(tempdir, 'mask.nii.gz')
+	result = errmap.run()
+	yield assert_equal, result.outputs.distance, 1.0
+
+	## Multi-Spectual
+	volume3 = np.array([[[1.0, 6.0], [0.0, 3.0]], [[1.0, 9.0], [3.0, 6.0]]]) # Raymond Vahan Damadian's birthday
+
+	msvolume1 = np.zeros(shape=(2,2,2,2))
+	msvolume1[:,:,:,0] = volume1
+	msvolume1[:,:,:,1] = volume3
+	msimg1 = nib.Nifti1Image(msvolume1, np.eye(4))
+
+	msvolume2 = np.zeros(shape=(2,2,2,2))
+	msvolume2[:,:,:,0] = volume3
+	msvolume2[:,:,:,1] = volume1
+	msimg2 = nib.Nifti1Image(msvolume2, np.eye(4))
+
+	nib.save(msimg1, os.path.join(tempdir, 'von-ray.nii.gz'))
+	nib.save(msimg2, os.path.join(tempdir, 'alan-ray.nii.gz'))
+
+	errmap.inputs.in_tst = os.path.join(tempdir, 'von-ray.nii.gz')
+	errmap.inputs.in_ref = os.path.join(tempdir, 'alan-ray.nii.gz')
+	errmap.inputs.metric = 'sqeuclidean'
+	result = errmap.run()
+	yield assert_equal, result.outputs.distance, 5.5
+
+	errmap.inputs.metric = 'euclidean'
+	result = errmap.run()
+	yield assert_equal, result.outputs.distance, np.float32(1.25 * (2**0.5))