MNT: Separate out vector image evaluation routine so it can be called programmatically

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