From f564d18746b14324e384740f28c7d857ee661081 Mon Sep 17 00:00:00 2001 From: sql Date: Mon, 2 Feb 2015 11:07:25 +1100 Subject: [PATCH 1/7] add average error to ErrorMap --- nipype/algorithms/metrics.py | 14 ++++++++++++-- 1 file changed, 12 insertions(+), 2 deletions(-) diff --git a/nipype/algorithms/metrics.py b/nipype/algorithms/metrics.py index 57a01d166d..119c2528bd 100644 --- a/nipype/algorithms/metrics.py +++ b/nipype/algorithms/metrics.py @@ -506,12 +506,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,12 +529,14 @@ def _run_interface(self, runtime): else: msk = np.ones(shape=mapshape) + # Vectorise 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 diffrernce if self.inputs.metric == 'sqeuclidean': errvector = diffvector**2 elif self.inputs.metric == 'euclidean': @@ -548,6 +551,12 @@ def _run_interface(self, runtime): errvectorexp = np.zeros_like(mskvector) errvectorexp[msk_idxs] = errvector + # Get averaged error + if self.inputs.metric == 'sqeuclidean': + self._distance = np.sqrt(np.sum(errvectorexp)) + elif self.inputs.metric == 'euclidean': + self._distance = np.average(errvectorexp) + errmap = errvectorexp.reshape(mapshape) hdr = nii_ref.get_header().copy() @@ -567,11 +576,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 From ba11e7c6ac99c158fc2cf163e39ae78c79928ee5 Mon Sep 17 00:00:00 2001 From: sql Date: Mon, 2 Feb 2015 11:12:11 +1100 Subject: [PATCH 2/7] added distance to the ErrorMapOutputSpec --- nipype/algorithms/metrics.py | 1 + 1 file changed, 1 insertion(+) diff --git a/nipype/algorithms/metrics.py b/nipype/algorithms/metrics.py index 119c2528bd..c953a6cc9f 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): From eebe773260047b2d074d1c055bb4478f3f975e82 Mon Sep 17 00:00:00 2001 From: siqi liu Date: Mon, 16 Feb 2015 19:40:30 +1100 Subject: [PATCH 3/7] fix masking; replace hstack with diffvector; remove square root of sqeuclidean average `errvectorexp = np.zeros_like(mskvector)` used to generate a uint8 vector by default which loses the precision when calculating the special volumes with float voxels of small values like the jacobian determinant volume --- nipype/algorithms/metrics.py | 11 ++++------- 1 file changed, 4 insertions(+), 7 deletions(-) diff --git a/nipype/algorithms/metrics.py b/nipype/algorithms/metrics.py index c953a6cc9f..514c7c6574 100644 --- a/nipype/algorithms/metrics.py +++ b/nipype/algorithms/metrics.py @@ -541,22 +541,19 @@ def _run_interface(self, runtime): if self.inputs.metric == 'sqeuclidean': errvector = diffvector**2 elif self.inputs.metric == 'euclidean': - X = np.hstack((refvector, tstvector)) - errvector = np.linalg.norm(X, axis=1) + #X = np.hstack((refvector, tstvector)) + 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 - if self.inputs.metric == 'sqeuclidean': - self._distance = np.sqrt(np.sum(errvectorexp)) - elif self.inputs.metric == 'euclidean': - self._distance = np.average(errvectorexp) + self._distance = np.average(errvectorexp) errmap = errvectorexp.reshape(mapshape) From 0a6422113d58b5a807881da78d7f046a399e68f2 Mon Sep 17 00:00:00 2001 From: siqi liu Date: Thu, 19 Feb 2015 16:57:29 +1100 Subject: [PATCH 4/7] average with only masked voxels; tests for ErrorMap --- nipype/algorithms/metrics.py | 2 +- nipype/algorithms/tests/test_errormap.py | 75 ++++++++++++++++++++++++ 2 files changed, 76 insertions(+), 1 deletion(-) create mode 100644 nipype/algorithms/tests/test_errormap.py diff --git a/nipype/algorithms/metrics.py b/nipype/algorithms/metrics.py index 514c7c6574..c9eb98a477 100644 --- a/nipype/algorithms/metrics.py +++ b/nipype/algorithms/metrics.py @@ -553,7 +553,7 @@ def _run_interface(self, runtime): errvectorexp[msk_idxs] = errvector # Get averaged error - self._distance = np.average(errvectorexp) + self._distance = np.average(errvector) errmap = errvectorexp.reshape(mapshape) diff --git a/nipype/algorithms/tests/test_errormap.py b/nipype/algorithms/tests/test_errormap.py new file mode 100644 index 0000000000..3729611f0b --- /dev/null +++ b/nipype/algorithms/tests/test_errormap.py @@ -0,0 +1,75 @@ +#!/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') + 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] = volume2 + msvolume2[:,:,:,1] = volume3 + msimg2 = nib.Nifti1Image(msvolume2, np.eye(4)) + + nib.save(img1, os.path.join(tempdir, 'von-ray.nii.gz')) + nib.save(img2, 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, 7.0 + + errmap.inputs.metric = 'euclidean' + result = errmap.run() + yield assert_equal, result.outputs.distance, 2**0.5 + 6 From ff345019f6cfd91260a34f313bcfa440efeda3d1 Mon Sep 17 00:00:00 2001 From: siqi liu Date: Thu, 19 Feb 2015 17:50:15 +1100 Subject: [PATCH 5/7] move comps check to only sqeuclidean --- nipype/algorithms/metrics.py | 15 +++++++-------- nipype/algorithms/tests/test_errormap.py | 13 +++++++------ 2 files changed, 14 insertions(+), 14 deletions(-) diff --git a/nipype/algorithms/metrics.py b/nipype/algorithms/metrics.py index c9eb98a477..108232c673 100644 --- a/nipype/algorithms/metrics.py +++ b/nipype/algorithms/metrics.py @@ -530,30 +530,29 @@ def _run_interface(self, runtime): else: msk = np.ones(shape=mapshape) - # Vectorise both volumes and make the pixel differennce + # 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 diffrernce + # Scale the diffrernce 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(diffvector, axis=1) - if (comps > 1): - errvector = np.sum(errvector, axis=1) - else: - errvector = np.squeeze(errvector) - 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) + self._distance = np.average(errvector) # Only average the masked voxels errmap = errvectorexp.reshape(mapshape) diff --git a/nipype/algorithms/tests/test_errormap.py b/nipype/algorithms/tests/test_errormap.py index 3729611f0b..5c6110a020 100644 --- a/nipype/algorithms/tests/test_errormap.py +++ b/nipype/algorithms/tests/test_errormap.py @@ -30,6 +30,7 @@ def test_errormap(): 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 @@ -57,19 +58,19 @@ def test_errormap(): msimg1 = nib.Nifti1Image(msvolume1, np.eye(4)) msvolume2 = np.zeros(shape=(2,2,2,2)) - msvolume2[:,:,:,0] = volume2 - msvolume2[:,:,:,1] = volume3 + msvolume2[:,:,:,0] = volume3 + msvolume2[:,:,:,1] = volume1 msimg2 = nib.Nifti1Image(msvolume2, np.eye(4)) - nib.save(img1, os.path.join(tempdir, 'von-ray.nii.gz')) - nib.save(img2, os.path.join(tempdir, 'alan-ray.nii.gz')) + 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, 7.0 + yield assert_equal, result.outputs.distance, 5.5 errmap.inputs.metric = 'euclidean' result = errmap.run() - yield assert_equal, result.outputs.distance, 2**0.5 + 6 + yield assert_equal, result.outputs.distance, np.float32(1.25 * (2**0.5)) From d7d7f046bc05ee1794eb7ef948f8095a6e52ccd3 Mon Sep 17 00:00:00 2001 From: siqi liu Date: Thu, 19 Feb 2015 17:58:40 +1100 Subject: [PATCH 6/7] Update CHANGES --- CHANGES | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CHANGES b/CHANGES index b791af7f9e..b26ecf9a90 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) * BUG: matplotlib is supposed to be optional (https://github.com/nipy/nipype/pull/1003) * FIX: Fix split_filename behaviour when path has no file component (https://github.com/nipy/nipype/pull/1035) * ENH: Updated FSL dtifit to include option for grad non-linearities (https://github.com/nipy/nipype/pull/1032) From a63f03c0cbe8b71be5f510959a0b126120e9be6f Mon Sep 17 00:00:00 2001 From: siqi liu Date: Fri, 20 Feb 2015 11:49:12 +1100 Subject: [PATCH 7/7] fix typo and remove the commented hstack --- nipype/algorithms/metrics.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/nipype/algorithms/metrics.py b/nipype/algorithms/metrics.py index 108232c673..bffce04161 100644 --- a/nipype/algorithms/metrics.py +++ b/nipype/algorithms/metrics.py @@ -537,7 +537,7 @@ def _run_interface(self, runtime): tstvector = tst_data.reshape(-1,comps)[msk_idxs].astype(np.float32) diffvector = (refvector-tstvector) - # Scale the diffrernce + # Scale the difference if self.inputs.metric == 'sqeuclidean': errvector = diffvector**2 if (comps > 1): @@ -545,7 +545,6 @@ def _run_interface(self, runtime): else: errvector = np.squeeze(errvector) elif self.inputs.metric == 'euclidean': - #X = np.hstack((refvector, tstvector)) errvector = np.linalg.norm(diffvector, axis=1) errvectorexp = np.zeros_like(mskvector, dtype=np.float32) # The default type is uint8