diff --git a/CHANGES b/CHANGES index 7449d38a19..976e04b4a2 100644 --- a/CHANGES +++ b/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) diff --git a/nipype/algorithms/metrics.py b/nipype/algorithms/metrics.py index 57a01d166d..bffce04161 100644 --- a/nipype/algorithms/metrics.py +++ b/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 diff --git a/nipype/algorithms/tests/test_errormap.py b/nipype/algorithms/tests/test_errormap.py new file mode 100644 index 0000000000..5c6110a020 --- /dev/null +++ b/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))