Merge pull request #261 from hsorby/matrix

Remove matrix utility file use cmlibs.maths 0.6.2 instead.

setup.py

@@ -25,7 +25,7 @@ def readfile(filename, split=False):
 # into the 'requirements.txt' file.
 requires = [
     # minimal requirements listing
-    "cmlibs.maths >= 0.6",
+    "cmlibs.maths >= 0.6.2",
     "cmlibs.utils >= 0.6",
     "cmlibs.zinc >= 4.1",
     "scipy",

src/scaffoldmaker/meshtypes/meshtype_3d_bladder1.py

@@ -5,7 +5,7 @@
 import copy
 import math
 
-from cmlibs.maths.vectorops import angle, set_magnitude, normalize, magnitude, cross
+from cmlibs.maths.vectorops import angle, set_magnitude, normalize, magnitude, cross, axis_angle_to_rotation_matrix
 from cmlibs.zinc.element import Element
 from cmlibs.zinc.field import Field
 from cmlibs.zinc.node import Node
@@ -16,7 +16,6 @@
 from scaffoldmaker.meshtypes.scaffold_base import Scaffold_base
 from scaffoldmaker.scaffoldpackage import ScaffoldPackage
 from scaffoldmaker.utils import interpolation as interp
-from scaffoldmaker.utils import matrix
 from scaffoldmaker.utils import tubemesh
 from scaffoldmaker.utils.geometry import createEllipsePoints
 from scaffoldmaker.utils.tracksurface import TrackSurface
@@ -812,7 +811,7 @@ def findNodesAlongBladderDome(sx_dome_group, elementsCountAround, elementsCountA
         rotAngle = n1 * 2.0 * math.pi / elementsCountAround
         rotAxis = normalize(cross(normalize(sx_dome_group[2][0]),
                                                         normalize(sx_dome_group[4][0])))
-        rotFrame = matrix.getRotationMatrixFromAxisAngle(rotAxis, rotAngle)
+        rotFrame = axis_angle_to_rotation_matrix(rotAxis, rotAngle)
         d2Rot = [rotFrame[j][0] * d2[0] + rotFrame[j][1] * d2[1] + rotFrame[j][2] * d2[2] for j in range(3)]
         d2Apex.append(d2Rot)
 

src/scaffoldmaker/meshtypes/meshtype_3d_bladderurethra1.py

@@ -7,7 +7,7 @@
 import copy
 import math
 
-from cmlibs.maths.vectorops import set_magnitude, cross, normalize, magnitude, angle
+from cmlibs.maths.vectorops import set_magnitude, cross, normalize, magnitude, angle, rotate_about_z_axis
 from cmlibs.utils.zinc.field import findOrCreateFieldGroup, \
     findOrCreateFieldStoredMeshLocation, findOrCreateFieldStoredString
 from cmlibs.zinc.element import Element
@@ -21,7 +21,6 @@
 from scaffoldmaker.meshtypes.scaffold_base import Scaffold_base
 from scaffoldmaker.scaffoldpackage import ScaffoldPackage
 from scaffoldmaker.utils import interpolation as interp
-from scaffoldmaker.utils import matrix
 from scaffoldmaker.utils import tubemesh
 from scaffoldmaker.utils.annulusmesh import createAnnulusMesh3d
 from scaffoldmaker.utils.geometry import createEllipsePoints
@@ -771,7 +770,7 @@ def generateBaseMesh(cls, region, options):
         innerNodes_x = []
         innerNodes_d1 = []
         innerNodes_d2 = []
-        nd1 = matrix.rotateAboutZAxis(nd2_max[0], 0.5 * math.pi)
+        nd1 = rotate_about_z_axis(nd2_max[0], 0.5 * math.pi)
         innerNodes_d1 += [nd1] * elementsCountAround
         for n1 in range(0, len(nx_max)):
             xAround, d1Around = createEllipsePoints([0.0, 0.0, nx_max[n1][2]], 2 * math.pi, [nx_max[n1][1], 0.0, 0.0],

src/scaffoldmaker/meshtypes/meshtype_3d_cecum1.py

@@ -7,7 +7,7 @@
 import copy
 import math
 
-from cmlibs.maths.vectorops import set_magnitude, normalize, magnitude, cross, dot
+from cmlibs.maths.vectorops import set_magnitude, normalize, magnitude, cross, dot, rotate_about_z_axis, axis_angle_to_rotation_matrix
 from cmlibs.utils.zinc.field import findOrCreateFieldCoordinates
 from cmlibs.zinc.element import Element
 from cmlibs.zinc.field import Field
@@ -22,7 +22,6 @@
 from scaffoldmaker.meshtypes.scaffold_base import Scaffold_base
 from scaffoldmaker.scaffoldpackage import ScaffoldPackage
 from scaffoldmaker.utils import interpolation as interp
-from scaffoldmaker.utils import matrix
 from scaffoldmaker.utils import tubemesh
 from scaffoldmaker.utils.annulusmesh import createAnnulusMesh3d
 from scaffoldmaker.utils.eftfactory_bicubichermitelinear import eftfactory_bicubichermitelinear
@@ -451,7 +450,7 @@ def getApexSegmentForCecum(xOuter, d1Outer, d2Outer, elementsCountAroundHalfHaus
             xFirstSegmentSampled.append(x)
             d2FirstSegmentSampled.append(d2)
             if n2 == 0:
-                d1 = matrix.rotateAboutZAxis(d2, math.pi*0.5)
+                d1 = rotate_about_z_axis(d2, math.pi*0.5)
                 d1FirstSegmentSampled.append(d1)
 
         if n2 > 0:
@@ -1217,12 +1216,12 @@ def createCecumMesh3d(region, options, networkLayout, nodeIdentifier, elementIde
 
     # sample along the midline of ostium
     rotAngle = math.pi
-    rotFrame = matrix.getRotationMatrixFromAxisAngle(d3ApexUnit, rotAngle)
+    rotFrame = axis_angle_to_rotation_matrix(d3ApexUnit, rotAngle)
     d1A = [rotFrame[j][0] * d1Cecum[1][0] + rotFrame[j][1] * d1Cecum[1][1] +
            rotFrame[j][2] * d1Cecum[1][2] for j in range(3)]
 
     rotAngle = math.pi
-    rotFrame = matrix.getRotationMatrixFromAxisAngle(normalize(d3Annulus[0]), rotAngle)
+    rotFrame = axis_angle_to_rotation_matrix(normalize(d3Annulus[0]), rotAngle)
     d2B = [rotFrame[j][0] * d1AnnulusOuter[0][0] + rotFrame[j][1] * d1AnnulusOuter[0][1] +
            rotFrame[j][2] * d1AnnulusOuter[0][2] for j in range(3)]
 

src/scaffoldmaker/meshtypes/meshtype_3d_colonsegment1.py

@@ -8,7 +8,7 @@
 import copy
 import math
 
-from cmlibs.maths.vectorops import normalize, set_magnitude, magnitude, dot, cross
+from cmlibs.maths.vectorops import normalize, set_magnitude, magnitude, dot, cross, rotate_about_z_axis
 from cmlibs.utils.zinc.field import findOrCreateFieldCoordinates
 from cmlibs.zinc.element import Element
 from cmlibs.zinc.field import Field
@@ -19,7 +19,6 @@
 from scaffoldmaker.annotation.cecum_terms import get_cecum_term
 from scaffoldmaker.meshtypes.scaffold_base import Scaffold_base
 from scaffoldmaker.utils import interpolation as interp
-from scaffoldmaker.utils import matrix
 from scaffoldmaker.utils import tubemesh
 from scaffoldmaker.utils.eftfactory_bicubichermitelinear import eftfactory_bicubichermitelinear
 from scaffoldmaker.utils.eftfactory_tricubichermite import eftfactory_tricubichermite
@@ -706,9 +705,9 @@ def getColonSegmentOuterPoints(region, elementsCountAroundTC, elementsCountAroun
                 if n1 > elementsCountAroundTC * 0.5:
                     # Non tenia coli
                     # Rotate about segment axis (z-axis)
-                    d2StartRot = matrix.rotateAboutZAxis(d2Start, radiansAround)
-                    d2MidRot = matrix.rotateAboutZAxis(d2Mid, radiansAround)
-                    d2EndRot = matrix.rotateAboutZAxis(d2End, radiansAround)
+                    d2StartRot = rotate_about_z_axis(d2Start, radiansAround)
+                    d2MidRot = rotate_about_z_axis(d2Mid, radiansAround)
+                    d2EndRot = rotate_about_z_axis(d2End, radiansAround)
 
                     d1 = [c * segmentLengthEndDerivativeFactor for c in d2StartRot]
                     d2 = [c * segmentLengthMidDerivativeFactor for c in d2MidRot]

src/scaffoldmaker/meshtypes/meshtype_3d_stellate1.py

@@ -6,7 +6,7 @@
 
 import math
 
-from cmlibs.maths.vectorops import set_magnitude
+from cmlibs.maths.vectorops import set_magnitude, rotate_about_z_axis
 from cmlibs.utils.zinc.field import findOrCreateFieldCoordinates, findOrCreateFieldGroup, \
     findOrCreateFieldStoredMeshLocation, findOrCreateFieldStoredString
 from cmlibs.zinc.element import Element
@@ -18,7 +18,6 @@
 from scaffoldmaker.utils.eft_utils import remapEftNodeValueLabel, scaleEftNodeValueLabels, setEftScaleFactorIds, remapEftLocalNodes
 from scaffoldmaker.utils.eftfactory_bicubichermitelinear import eftfactory_bicubichermitelinear
 from scaffoldmaker.utils.interpolation import smoothCubicHermiteDerivativesLine
-from scaffoldmaker.utils.matrix import rotateAboutZAxis
 from scaffoldmaker.utils.meshrefinement import MeshRefinement
 
 
@@ -683,7 +682,7 @@ def createArm(halfArmArcAngleRadians, elementLengths, elementLengthCentral, elem
                     elif (e1 == 0) and ((e2 == 0) or (e2 == elementsCount2)):
                         ds2 = [dcent[0], -dcent[1], dcent[2]] if (e2 == 0) else dcent
                         ds2 = set_magnitude(ds2, dipMag)
-                        ds1 = rotateAboutZAxis(ds2, -math.pi / 2)
+                        ds1 = rotate_about_z_axis(ds2, -math.pi / 2)
                     elif e1 == elementsCount1 and e2 == elementsCount2-1: # armEnd
                         ds1 = [elementWidth,0,0]
                         ds2 = [0, 1 * (elementLength + elementWidth), 0]
@@ -707,15 +706,15 @@ def createArm(halfArmArcAngleRadians, elementLengths, elementLengthCentral, elem
     # rotate entire arm about origin by armAngle except for centre nodes
     tol = 1e-12
     for j, n in enumerate(x):
-        xynew = rotateAboutZAxis(n, armAngle)
+        xynew = rotate_about_z_axis(n, armAngle)
         [xnew, ynew] = [r for r in xynew][:2]
         if (abs(xnew) < tol):
             xnew = 0
         if (abs(ynew) < tol):
             ynew = 0
         x[j] = [xnew, ynew, x[j][2]]
 
-        xnodes_ds1[j] = rotateAboutZAxis(xnodes_ds1[j], armAngle)
-        xnodes_ds2[j] = rotateAboutZAxis(xnodes_ds2[j], armAngle)
+        xnodes_ds1[j] = rotate_about_z_axis(xnodes_ds1[j], armAngle)
+        xnodes_ds2[j] = rotate_about_z_axis(xnodes_ds2[j], armAngle)
 
     return (x, xnodes_ds1, xnodes_ds2, rmVertexNodes)

src/scaffoldmaker/meshtypes/meshtype_3d_stomach1.py

@@ -9,7 +9,7 @@
 import copy
 import math
 
-from cmlibs.maths.vectorops import cross, sub, set_magnitude, normalize, magnitude
+from cmlibs.maths.vectorops import cross, sub, set_magnitude, normalize, magnitude, axis_angle_to_rotation_matrix
 from cmlibs.utils.zinc.field import find_or_create_field_coordinates
 from cmlibs.utils.zinc.finiteelement import get_element_node_identifiers, get_maximum_element_identifier, \
     get_maximum_node_identifier
@@ -27,7 +27,6 @@
 from scaffoldmaker.meshtypes.scaffold_base import Scaffold_base
 from scaffoldmaker.scaffoldpackage import ScaffoldPackage
 from scaffoldmaker.utils import interpolation as interp
-from scaffoldmaker.utils import matrix
 from scaffoldmaker.utils.annulusmesh import createAnnulusMesh3d
 from scaffoldmaker.utils.eftfactory_bicubichermitelinear import eftfactory_bicubichermitelinear
 from scaffoldmaker.utils.eftfactory_tricubichermite import eftfactory_tricubichermite
@@ -1505,7 +1504,7 @@ def createStomachMesh3d(region, fm, coordinates, stomachTermsAlong, allAnnotatio
     d2ApexAround = [cross(cross(rotAxisApex, sub(tpx, cxApex)), rotAxisApex) for tpx in px]
 
     rotAngle = -math.pi * 0.5
-    rotFrame = matrix.getRotationMatrixFromAxisAngle(rotAxisApex, rotAngle)
+    rotFrame = axis_angle_to_rotation_matrix(rotAxisApex, rotAngle)
     for n in range(len(px)):
         d1ApexAround.append([rotFrame[j][0] * d2ApexAround[n][0] + rotFrame[j][1] * d2ApexAround[n][1] +
                              rotFrame[j][2] * d2ApexAround[n][2] for j in range(3)])
@@ -1912,7 +1911,7 @@ def createStomachMesh3d(region, fm, coordinates, stomachTermsAlong, allAnnotatio
                 # Sample from apex to annulus
                 bPosition = xAnnulusOuterPosition[annulusIdxAtBodyStartIdxMinusOne[count]]
                 d = d2AnnulusOuter[annulusIdxAtBodyStartIdxMinusOne[count]]
-                rotFrame = matrix.getRotationMatrixFromAxisAngle(d3Annulus[annulusIdxAtBodyStartIdxMinusOne[count]],
+                rotFrame = axis_angle_to_rotation_matrix(d3Annulus[annulusIdxAtBodyStartIdxMinusOne[count]],
                                                                  math.pi)
                 endDerivative = [rotFrame[j][0] * d[0] + rotFrame[j][1] * d[1] + rotFrame[j][2] * d[2]
                                  for j in range(3)]
@@ -1941,7 +1940,7 @@ def createStomachMesh3d(region, fm, coordinates, stomachTermsAlong, allAnnotatio
 
             # Rotate nd2
             for m in range(len(nx)):
-                rotFrame = matrix.getRotationMatrixFromAxisAngle(nd3[m], math.pi)
+                rotFrame = axis_angle_to_rotation_matrix(nd3[m], math.pi)
                 nd2[m] = [rotFrame[j][0] * nd2[m][0] + rotFrame[j][1] * nd2[m][1] +
                           rotFrame[j][2] * nd2[m][2] for j in range(3)]
 
@@ -1956,7 +1955,7 @@ def createStomachMesh3d(region, fm, coordinates, stomachTermsAlong, allAnnotatio
                 elif n1 == elementsAroundHalfDuod - 1:
                     for m in range(2 * (elementsAroundQuarterEso - 2) + 1):
                         annulusIdx = m + 2
-                        rotFrame = matrix.getRotationMatrixFromAxisAngle(d3Annulus[annulusIdx], math.pi)
+                        rotFrame = axis_angle_to_rotation_matrix(d3Annulus[annulusIdx], math.pi)
                         d2 = d2AnnulusOuter[annulusIdx]
                         d2 = [rotFrame[j][0] * d2[0] + rotFrame[j][1] * d2[1] + rotFrame[j][2] * d2[2]
                               for j in range(3)]
@@ -1968,7 +1967,7 @@ def createStomachMesh3d(region, fm, coordinates, stomachTermsAlong, allAnnotatio
                 elif n1 == elementsAroundHalfDuod + 1:
                     for m in range(2 * (elementsAroundQuarterEso - 2) + 1):
                         annulusIdx = -2 - m
-                        rotFrame = matrix.getRotationMatrixFromAxisAngle(d3Annulus[annulusIdx], math.pi)
+                        rotFrame = axis_angle_to_rotation_matrix(d3Annulus[annulusIdx], math.pi)
                         d1 = d1AnnulusOuter[annulusIdx]
                         d1 = [rotFrame[j][0] * d1[0] + rotFrame[j][1] * d1[1] + rotFrame[j][2] * d1[2]
                               for j in range(3)]

src/scaffoldmaker/utils/tubemesh.py

@@ -6,14 +6,13 @@
 
 import math
 
-from cmlibs.maths.vectorops import normalize, dot, cross, magnitude, set_magnitude
+from cmlibs.maths.vectorops import normalize, dot, cross, magnitude, set_magnitude, axis_angle_to_rotation_matrix
 from cmlibs.utils.zinc.field import findOrCreateFieldCoordinates
 from cmlibs.zinc.element import Element
 from cmlibs.zinc.field import Field
 from cmlibs.zinc.node import Node
 from scaffoldmaker.annotation.annotationgroup import mergeAnnotationGroups
 from scaffoldmaker.utils import interpolation as interp
-from scaffoldmaker.utils import matrix
 from scaffoldmaker.utils.eftfactory_bicubichermitelinear import eftfactory_bicubichermitelinear
 from scaffoldmaker.utils.eftfactory_tricubichermite import eftfactory_tricubichermite
 from scaffoldmaker.utils.eft_utils import remapEftNodeValueLabelsVersion
@@ -124,14 +123,14 @@ def warpSegmentPoints(xList, d1List, d2List, segmentAxis, sx, sd1, sd2, elements
         if magnitude(cp)> 0.0: # path tangent not parallel to segment axis
             axisRot = normalize(cp)
             thetaRot = math.acos(dot(segmentAxis, unitTangent))
-            rotFrame = matrix.getRotationMatrixFromAxisAngle(axisRot, thetaRot)
+            rotFrame = axis_angle_to_rotation_matrix(axisRot, thetaRot)
             centroidRot = [rotFrame[j][0]*centroid[0] + rotFrame[j][1]*centroid[1] + rotFrame[j][2]*centroid[2] for j in range(3)]
 
         else: # path tangent parallel to segment axis (z-axis)
             if dp == -1.0: # path tangent opposite direction to segment axis
                 thetaRot = math.pi
                 axisRot = [1.0, 0, 0]
-                rotFrame = matrix.getRotationMatrixFromAxisAngle(axisRot, thetaRot)
+                rotFrame = axis_angle_to_rotation_matrix(axisRot, thetaRot)
                 centroidRot = [rotFrame[j][0] * centroid[0] + rotFrame[j][1] * centroid[1] + rotFrame[j][2] * centroid[2] for j in range(3)]
 
             else: # segment axis in same direction as unit tangent
@@ -167,7 +166,7 @@ def warpSegmentPoints(xList, d1List, d2List, segmentAxis, sx, sd1, sd2, elements
                         thetaRot2 = math.acos(
                             dot(normalize(vectorToFirstNode), normalize(sd2[nAlongSegment])))
                         axisRot2 = unitTangent
-                        rotFrame2 = matrix.getRotationMatrixFromAxisAngle(axisRot2, signThetaRot2*thetaRot2)
+                        rotFrame2 = axis_angle_to_rotation_matrix(axisRot2, signThetaRot2*thetaRot2)
                     else:
                         rotFrame2 = [[1, 0, 0], [0, 1, 0], [0, 0, 1]]
                 else:

tests/test_gastrointestinaltract1.py

@@ -113,7 +113,7 @@ def test_gastrointestinaltract1(self):
         self.assertAlmostEqual(surfaceArea, 280350.6462177311, delta=1.0E-6)
         result, volume = volumeField.evaluateReal(fieldcache, 1)
         self.assertEqual(result, RESULT_OK)
-        self.assertAlmostEqual(volume, 598505.4701647629, delta=1.0E-6)
+        self.assertAlmostEqual(volume, 598505.4701647585, delta=1.0E-3)
 
         # check some annotationGroups:
         expectedSizes3d = {
@@ -124,15 +124,16 @@ def test_gastrointestinaltract1(self):
             "colon": 4752
         }
         for name in expectedSizes3d:
+            term = None
             if name == "esophagus":
                 term = get_esophagus_term(name)
-            elif name == ("stomach"):
+            elif name == "stomach":
                 term = get_stomach_term(name)
-            elif name == ("small intestine"):
+            elif name == "small intestine":
                 term = get_smallintestine_term(name)
-            elif name == ("caecum"):
+            elif name == "caecum":
                 term = get_cecum_term(name)
-            elif name == ("colon"):
+            elif name == "colon":
                 term = get_colon_term(name)
             group = getAnnotationGroupForTerm(annotationGroups, term)
             size = group.getMeshGroup(mesh3d).getSize()