diff --git a/setup.py b/setup.py index 14a8c6d9..5434b21c 100644 --- a/setup.py +++ b/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", diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_bladder1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_bladder1.py index be3d8433..a7abd62c 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_bladder1.py +++ b/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) diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_bladderurethra1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_bladderurethra1.py index 52417285..f33447b8 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_bladderurethra1.py +++ b/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], diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_cecum1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_cecum1.py index 0251809e..182acf68 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_cecum1.py +++ b/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)] diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_colonsegment1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_colonsegment1.py index a269c0cd..cd43724a 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_colonsegment1.py +++ b/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] diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_stellate1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_stellate1.py index f8fcef16..8bb70527 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_stellate1.py +++ b/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,7 +706,7 @@ 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 @@ -715,7 +714,7 @@ def createArm(halfArmArcAngleRadians, elementLengths, elementLengthCentral, elem 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) diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_stomach1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_stomach1.py index 10024d4e..39a05c7f 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_stomach1.py +++ b/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)] diff --git a/src/scaffoldmaker/utils/matrix.py b/src/scaffoldmaker/utils/matrix.py deleted file mode 100644 index 23e09d45..00000000 --- a/src/scaffoldmaker/utils/matrix.py +++ /dev/null @@ -1,37 +0,0 @@ -''' -Utility functions for matrix. -''' - -import math - -def getRotationMatrixFromAxisAngle(rotAxis, theta): - """ - Generate the rotation matrix for rotation about an axis. - :param rotAxis: axis of rotation - :param theta: angle of rotation - :return: rotation matrix - """ - cosTheta = math.cos(theta) - sinTheta = math.sin(theta) - C = 1 - cosTheta - rotMatrix = ([[rotAxis[0]*rotAxis[0]*C + cosTheta, rotAxis[0]*rotAxis[1]*C - rotAxis[2]*sinTheta, rotAxis[0]*rotAxis[2]*C + rotAxis[1]*sinTheta], - [rotAxis[1]*rotAxis[0]*C + rotAxis[2]*sinTheta, rotAxis[1]*rotAxis[1]*C + cosTheta, rotAxis[1]*rotAxis[2]*C - rotAxis[0]*sinTheta], - [rotAxis[2]*rotAxis[0]*C - rotAxis[1]*sinTheta, rotAxis[2]*rotAxis[1]*C + rotAxis[0]*sinTheta, rotAxis[2]*rotAxis[2]*C + cosTheta]]) - - return rotMatrix - -def rotateAboutZAxis(x, theta): - """ - Rotates matrix about z-axis. - : param x: matrix to be rotated. - : param theta: angle of rotation. - :return rotated matrix - """ - cosTheta = math.cos(theta) - sinTheta = math.sin(theta) - - xRot = [x[0]*cosTheta - x[1]*sinTheta, - x[0]*sinTheta + x[1]*cosTheta, - x[2]] - - return xRot diff --git a/src/scaffoldmaker/utils/tubemesh.py b/src/scaffoldmaker/utils/tubemesh.py index e8c04ccc..5e0f9441 100644 --- a/src/scaffoldmaker/utils/tubemesh.py +++ b/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: diff --git a/tests/test_gastrointestinaltract1.py b/tests/test_gastrointestinaltract1.py index 442360f8..3b1e93c7 100644 --- a/tests/test_gastrointestinaltract1.py +++ b/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()