diff --git a/CHANGELOG.rst b/CHANGELOG.rst
index fa4650f27400172b7d6c2205af0ee1174e9f2783..026bfcc8b41ecad91dc915351d54e19158ce3f74 100644
--- a/CHANGELOG.rst
+++ b/CHANGELOG.rst
@@ -3,6 +3,14 @@
+0.12.3 (Friday 8th December 2023)
+---------------------------------
+
+
+* Separate out the core logic of ``evalVectorImage`` so it can be called
+ programmatically more easily.
+
+
0.12.2 (Thursday 13th July 2023)
--------------------------------
diff --git a/pyfeeds/__init__.py b/pyfeeds/__init__.py
index 291bb2d8188a2d207d035cf716107cf0306655ab..ea50d8c9d7f9e6bb2d6b6fd5484ac0e5120dd733 100644
--- a/pyfeeds/__init__.py
+++ b/pyfeeds/__init__.py
@@ -5,7 +5,7 @@
# Author: Paul McCarthy <pauldmccarthy@gmail.com>
#
-__version__ = '0.12.2'
+__version__ = '0.12.3'
"""The pyfeeds version number. """
@@ -20,5 +20,7 @@ from .evaluate import (
evalImage,
evalNumericalText,
evalVectorImage,
- evalMD5
+ evalMD5,
+ cmpArrays,
+ cmpVectorArrays
)
diff --git a/pyfeeds/evaluate.py b/pyfeeds/evaluate.py
index 61650040a1a5bb105174db4fa7f598b5fe8314df..b40f108cbe3108c720f19fb254815721ed6d0866 100644
--- a/pyfeeds/evaluate.py
+++ b/pyfeeds/evaluate.py
@@ -474,54 +474,11 @@ def evalVectorImage(pyf, testfile, benchmark):
both vector images, where each voxel contains a 3-dimensional
vector, undirected, and centered at 0.
"""
- if evalImage(pyf, testfile, benchmark) == 0.0:
- return 0
- pion2 = np.pi / 2
- data1 = pyf.imageCache[testfile][ 1].reshape(-1, 3).T
- data2 = pyf.imageCache[benchmark][1].reshape(-1, 3).T
- # Calculate the length of each vector,
- # discard vectors of length 0, and
- # normalise each vector to unit length
- len1 = np.linalg.norm(data1, axis=0)
- len2 = np.linalg.norm(data2, axis=0)
- nz1 = len1 > 1e-6
- nz2 = len2 > 1e-6
- nz = nz1 & nz2
+ data1 = pyf.imageCache[testfile][ 1]
+ data2 = pyf.imageCache[benchmark][1]
- data1 = data1[:, nz] / len1[nz]
- data2 = data2[:, nz] / len2[nz]
- len1 = len1[nz]
- len2 = len2[nz]
-
- # Calculate the angle between each vector.
- # Vectors are undirected, and centered at
- # (0, 0, 0), so the maximum possible angle
- # we can have is 90 degrees.
- dot = np.sum(data1 * data2, axis=0)
- dot = np.clip(dot, -1, 1)
- angle = np.arccos(dot)
- amask = angle > pion2
- angle[amask] = np.pi - angle[amask]
-
- # We also compare the length of each
- # vector, and the pattern of missing
- # voxels (vectors of length 0)
- nzcorr = np.abs(np.corrcoef(nz1, nz2)[0, 1])
- lendiff = np.abs(len1 - len2) / np.max((len1, len2), axis=0)
- angle = np.abs(angle) / pion2
-
- if np.isnan(nzcorr):
- nzcorr = 1
-
- # All errors are normalised to
- # the range (0, 1). We return
- # the worst error
- nzError = 1 - nzcorr
- angleError = angle.mean()
- lenError = lendiff.mean()
-
- return max((nzError, lenError, angleError))
+ return cmpVectorArrays(data1, data2)
def evalPolarCoordinateImageGroup(pyf, testfiles, benchmarks):
@@ -654,3 +611,61 @@ def cmpArrays(arr1, arr2):
# The final error is the mean error across all voxels
return normdiff.mean()
+
+
+def cmpVectorArrays(arr1, arr2):
+ """Compare two (X, Y, Z, 3) arrays containing vector fields. Difference
+ between vectoer lengths and angles are calculated, and the maximum of the
+ mean length/angle difference is returned (normalised to lie between 0 and
+ 1).
+ """
+
+ if np.isclose(cmpArrays(arr1, arr2), 0):
+ return 0
+
+ pion2 = np.pi / 2
+ arr1 = arr1.reshape(-1, 3).T
+ arr2 = arr2.reshape(-1, 3).T
+
+ # Calculate the length of each vector,
+ # discard vectors of length 0, and
+ # normalise each vector to unit length
+ len1 = np.linalg.norm(arr1, axis=0)
+ len2 = np.linalg.norm(arr2, axis=0)
+ nz1 = len1 > 1e-6
+ nz2 = len2 > 1e-6
+ nz = nz1 & nz2
+
+ arr1 = arr1[:, nz] / len1[nz]
+ arr2 = arr2[:, nz] / len2[nz]
+ len1 = len1[nz]
+ len2 = len2[nz]
+
+ # Calculate the angle between each vector.
+ # Vectors are undirected, and centered at
+ # (0, 0, 0), so the maximum possible angle
+ # we can have is 90 degrees.
+ dot = np.sum(arr1 * arr2, axis=0)
+ dot = np.clip(dot, -1, 1)
+ angle = np.arccos(dot)
+ amask = angle > pion2
+ angle[amask] = np.pi - angle[amask]
+
+ # We also compare the length of each
+ # vector, and the pattern of missing
+ # voxels (vectors of length 0)
+ nzcorr = np.abs(np.corrcoef(nz1, nz2)[0, 1])
+ lendiff = np.abs(len1 - len2) / np.max((len1, len2), axis=0)
+ angle = np.abs(angle) / pion2
+
+ if np.isnan(nzcorr):
+ nzcorr = 1
+
+ # All errors are normalised to
+ # the range (0, 1). We return
+ # the worst error
+ nzError = 1 - nzcorr
+ angleError = angle.mean()
+ lenError = lendiff.mean()
+
+ return max((nzError, lenError, angleError))
diff --git a/pyfeeds/tests/__init__.py b/pyfeeds/tests/__init__.py
index db6dbb57323ea11effc3d59643923356b184c309..7feda3307a62c5d69ea0a8dbe7b06bfd32cd5bd6 100644
--- a/pyfeeds/tests/__init__.py
+++ b/pyfeeds/tests/__init__.py
@@ -113,6 +113,13 @@ def maketest(filename, returnCode=0, inputs=None, outputs=None, stdout=None):
def makepyfeeds(**kwargs):
+
+ # simplest way of creating a dummy pyfeeds object
+ if len(kwargs) == 0:
+ kwargs['command'] = 'compare'
+ kwargs['inputDir'] = os.getcwd()
+ kwargs['benchmarkDir'] = os.getcwd()
+
args = argparse.Namespace(**kwargs)
cfg = argparse.Namespace()
return main.Pyfeeds(args, cfg)
diff --git a/pyfeeds/tests/test_evaluate.py b/pyfeeds/tests/test_evaluate.py
index 96ff097c102fe24a0b236fb79fe74c29597ad5c6..c016e787f1be52fc113e7a6d095b90031be2f691 100644
--- a/pyfeeds/tests/test_evaluate.py
+++ b/pyfeeds/tests/test_evaluate.py
@@ -8,6 +8,9 @@
import os
import os.path as op
+import numpy as np
+import nibabel as nib
+
from . import tempdir, makepaths, maketest, makepyfeeds, CaptureStdout
from pyfeeds import testing, evaluate
@@ -57,3 +60,41 @@ def test_evaluateTestAgainstBenchmark():
assert len(lines) == 3
for l in lines[1:]:
assert 'PASS' in l
+
+
+def test_evalVectorImage():
+
+ vecarr1 = -1 + 2 * np.random.random((10, 10, 10, 3))
+ vecarr2 = -1 + 2 * np.random.random((10, 10, 10, 3))
+
+ with tempdir():
+
+ pyf = makepyfeeds()
+ fname1 = 'image1.nii.gz'
+ fname2 = 'image2.nii.gz'
+
+ nib.Nifti1Image(vecarr1, np.eye(4)).to_filename(fname1)
+ nib.Nifti1Image(vecarr2, np.eye(4)).to_filename(fname2)
+
+ assert evaluate.evalVectorImage(pyf, fname1, fname1) == 0
+ assert evaluate.evalVectorImage(pyf, fname2, fname2) == 0
+ assert evaluate.evalVectorImage(pyf, fname1, fname2) != 0
+
+
+def test_evalImage():
+
+ arr1 = -1 + 2 * np.random.random((10, 10, 10, 10))
+ arr2 = -1 + 2 * np.random.random((10, 10, 10, 10))
+
+ with tempdir():
+
+ pyf = makepyfeeds()
+ fname1 = 'image1.nii.gz'
+ fname2 = 'image2.nii.gz'
+
+ nib.Nifti1Image(arr1, np.eye(4)).to_filename(fname1)
+ nib.Nifti1Image(arr2, np.eye(4)).to_filename(fname2)
+
+ assert evaluate.evalImage(pyf, fname1, fname1) == 0
+ assert evaluate.evalImage(pyf, fname2, fname2) == 0
+ assert evaluate.evalImage(pyf, fname1, fname2) != 0