Merge pull request #1039 from lsqshr/enh_add_average_error_2_errormap

ENH: add average error to ErrorMap

Author: oesteban

CHANGES

@@ -1,6 +1,6 @@
 Next release
 ============
-
+* 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)

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))