From 3f60540ea69c9a7a2272a034037212d3af2b70a8 Mon Sep 17 00:00:00 2001 From: Hugh Sorby Date: Tue, 23 Jul 2024 18:33:04 +1200 Subject: [PATCH] Switch over to using cmlibs.maths 0.6.0 for math functions. --- setup.py | 2 +- .../meshtypes/meshtype_3d_bladder1.py | 38 +++--- .../meshtypes/meshtype_3d_bladderurethra1.py | 23 ++-- .../meshtypes/meshtype_3d_bone1.py | 14 +- .../meshtypes/meshtype_3d_cecum1.py | 107 ++++++++------- .../meshtypes/meshtype_3d_colon1.py | 3 +- .../meshtypes/meshtype_3d_colonsegment1.py | 41 +++--- .../meshtypes/meshtype_3d_esophagus1.py | 9 +- .../meshtype_3d_heartarterialroot1.py | 19 ++- .../meshtype_3d_heartarterialvalve1.py | 7 +- .../meshtypes/meshtype_3d_heartatria1.py | 103 ++++++++------- .../meshtypes/meshtype_3d_heartatria2.py | 65 +++++----- .../meshtypes/meshtype_3d_heartventricles1.py | 23 ++-- .../meshtypes/meshtype_3d_heartventricles2.py | 21 ++- .../meshtypes/meshtype_3d_heartventricles3.py | 85 ++++++------ .../meshtype_3d_heartventriclesbase1.py | 29 ++--- .../meshtype_3d_heartventriclesbase2.py | 35 +++-- .../meshtypes/meshtype_3d_ostium1.py | 40 +++--- .../meshtypes/meshtype_3d_ostium2.py | 73 ++++++----- .../meshtypes/meshtype_3d_smallintestine1.py | 7 +- .../meshtypes/meshtype_3d_solidsphere1.py | 13 +- .../meshtypes/meshtype_3d_solidsphere2.py | 9 +- .../meshtypes/meshtype_3d_stomach1.py | 77 ++++++----- .../meshtypes/meshtype_3d_uterus1.py | 25 ++-- src/scaffoldmaker/scaffoldpackage.py | 5 +- src/scaffoldmaker/utils/annulusmesh.py | 40 +++--- src/scaffoldmaker/utils/cylindermesh.py | 57 ++++---- .../utils/eftfactory_tricubichermite.py | 22 ++-- src/scaffoldmaker/utils/geometry.py | 35 +++-- src/scaffoldmaker/utils/interpolation.py | 35 +++-- src/scaffoldmaker/utils/mirror.py | 20 +-- src/scaffoldmaker/utils/shieldmesh.py | 14 +- src/scaffoldmaker/utils/spheremesh.py | 122 +++++++++--------- src/scaffoldmaker/utils/tracksurface.py | 41 +++--- src/scaffoldmaker/utils/tubemesh.py | 67 +++++----- src/scaffoldmaker/utils/vector.py | 119 ----------------- src/scaffoldmaker/utils/zinc_utils.py | 18 +-- 37 files changed, 656 insertions(+), 807 deletions(-) delete mode 100644 src/scaffoldmaker/utils/vector.py diff --git a/setup.py b/setup.py index 5ebe0e3b..14a8c6d9 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.3", + "cmlibs.maths >= 0.6", "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 ef8eb100..5b407734 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_bladder1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_bladder1.py @@ -5,6 +5,7 @@ import copy import math +from cmlibs.maths.vectorops import angle from cmlibs.zinc.element import Element from cmlibs.zinc.field import Field from cmlibs.zinc.node import Node @@ -17,7 +18,6 @@ from scaffoldmaker.utils import interpolation as interp from scaffoldmaker.utils import matrix from scaffoldmaker.utils import tubemesh -from scaffoldmaker.utils import vector from scaffoldmaker.utils.geometry import createEllipsePoints from scaffoldmaker.utils.tracksurface import TrackSurface from scaffoldmaker.utils.zinc_utils import exnode_string_from_nodeset_field_parameters @@ -528,7 +528,7 @@ def generateBaseMesh(cls, region, options): # v1 = xEllipsoid[n] # v2 = xEllipsoid[n + elementsCountAround // 2] # v1v2 = [v2[c] - v1[c] for c in range(3)] - # mag = vector.magnitude(v1v2) + # mag = magnitude(v1v2) # halfMag = mag / 2 # d2path.append([0.0, -halfMag, 0.0]) # d3path.append([0.0, 0.0, halfMag]) @@ -567,7 +567,7 @@ def generateBaseMesh(cls, region, options): # Flat coordinates xFlat = d1Flat = d2Flat = [] - # urethraOpeningRadius = vector.magnitude(sx_neck_group[2][-1]) + # urethraOpeningRadius = magnitude(sx_neck_group[2][-1]) # xFlat, d1Flat, d2Flat = obtainBladderFlatNodes(elementsCountAlongBladder, elementsCountAround, # elementsCountThroughWall, xFinal, d1Final, d2Final, # bladderCentralPathLength, urethraOpeningRadius, wallThickness) @@ -798,7 +798,7 @@ def findDerivativeBetweenPoints(v1, v2): """ d = [v2[c] - v1[c] for c in range(3)] arcLengthAround = interp.computeCubicHermiteArcLength(v1, d, v2, d, True) - d = [c * arcLengthAround for c in vector.normalise(d)] + d = [c * arcLengthAround for c in normalize(d)] return d @@ -810,8 +810,8 @@ def findNodesAlongBladderDome(sx_dome_group, elementsCountAround, elementsCountA d2 = sx_dome_group[2][0] for n1 in range(elementsCountAround): rotAngle = n1 * 2.0 * math.pi / elementsCountAround - rotAxis = vector.normalise(vector.crossproduct3(vector.normalise(sx_dome_group[2][0]), - vector.normalise(sx_dome_group[4][0]))) + rotAxis = normalize(cross(normalize(sx_dome_group[2][0]), + normalize(sx_dome_group[4][0]))) rotFrame = matrix.getRotationMatrixFromAxisAngle(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) @@ -946,9 +946,9 @@ def findNodesAlongBladderNeck(sx_dome_group, sx_neck_group, d2SampledDome, domeS v1 = xEllipses_neck[n2 - 1][n1] v2 = xEllipses_neck[n2][n1] d2vec = findDerivativeBetweenPoints(v1, v2) - d2mag = vector.magnitude(d2vec) - d2dir = vector.normalise(sx_neck_group[1][-1]) - d2 = vector.setMagnitude(d2dir, d2mag) + d2mag = magnitude(d2vec) + d2dir = normalize(sx_neck_group[1][-1]) + d2 = set_magnitude(d2dir, d2mag) else: v1 = xEllipses_neck[n2 - 1][n1] v2 = xEllipses_neck[n2][n1] @@ -1050,7 +1050,7 @@ def obtainBladderFlatNodes(elementsCountAlongBladder, elementsCountAround, eleme # Find the angle at the bottom of the bladder neck v1 = [0.0, 0.0, bladderLength] v2 = [urethraOpeningRadius, 0.0, bladderLength] - alpha = vector.angleBetweenVectors(v1, v2) + alpha = angle(v1, v2) # Find apex to urethra arcLength in minor radius xApexInner = xFinal[0] @@ -1167,8 +1167,8 @@ def obtainBladderFlatNodes(elementsCountAlongBladder, elementsCountAround, eleme v3 = xfnListRearranged[elementsCountAround // 2] v21 = [v2[c] - v1[c] for c in range(3)] v31 = [v3[c] - v1[c] for c in range(3)] - d1Mag = vector.magnitude(v21) - d2Mag = vector.magnitude(v31) + d1Mag = magnitude(v21) + d2Mag = magnitude(v31) # Add apex nodes to the list xFlat = [] @@ -1314,8 +1314,8 @@ def getBladderCoordinates(elementsCountAlongDome, elementsCountAlongNeck, elemen for n2 in range(1, elementsCountAlongBladder + 1): d3Around = [] for n1 in range(elementsCountAround): - d3Around.append(vector.normalise( - vector.crossproduct3(vector.normalise(d1SampledAll[n2][n1]), vector.normalise(d2SampledAll[n2][n1])))) + d3Around.append(normalize( + cross(normalize(d1SampledAll[n2][n1]), normalize(d2SampledAll[n2][n1])))) d3UnitOuter.append(d3Around) # Inner nodes @@ -1349,12 +1349,12 @@ def getBladderCoordinates(elementsCountAlongDome, elementsCountAlongNeck, elemen # Arclength between apex point and corresponding point on next face mag = interp.getCubicHermiteArcLength(xList[n - 1], d2List[n - 1], xList[2 * n - 1], d2List[2 * n - 1]) - d2ApexOuter = vector.setMagnitude(sx_dome_group[2][0], mag) + d2ApexOuter = set_magnitude(sx_dome_group[2][0], mag) - d1ApexOuter = vector.crossproduct3(sx_dome_group[1][0], d2ApexOuter) - d1ApexOuter = vector.setMagnitude(d1ApexOuter, mag) - d3ApexUnit = vector.normalise( - vector.crossproduct3(vector.normalise(d1ApexOuter), vector.normalise(d2ApexOuter))) + d1ApexOuter = cross(sx_dome_group[1][0], d2ApexOuter) + d1ApexOuter = set_magnitude(d1ApexOuter, mag) + d3ApexUnit = normalize( + cross(normalize(d1ApexOuter), normalize(d2ApexOuter))) d3ApexOuter = [d3ApexUnit[c] * wallThickness / elementsCountThroughWall for c in range(3)] # Final nodes on the bladder diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_bladderurethra1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_bladderurethra1.py index 71f035d5..59b675b4 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_bladderurethra1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_bladderurethra1.py @@ -22,7 +22,6 @@ from scaffoldmaker.utils import interpolation as interp from scaffoldmaker.utils import matrix from scaffoldmaker.utils import tubemesh -from scaffoldmaker.utils import vector from scaffoldmaker.utils.annulusmesh import createAnnulusMesh3d from scaffoldmaker.utils.geometry import createEllipsePoints from scaffoldmaker.utils.interpolation import smoothCubicHermiteDerivativesLine @@ -871,7 +870,7 @@ def generateBaseMesh(cls, region, options): innerRadiusAlong = [] for n2 in range(elementsCountAlong + 1): firstNodeAlong = innerNodes_x[n2 * elementsCountAround] - radius = vector.magnitude(firstNodeAlong[c] - [0.0, 0.0, firstNodeAlong[2]][c] for c in range(3)) + radius = magnitude(firstNodeAlong[c] - [0.0, 0.0, firstNodeAlong[2]][c] for c in range(3)) innerRadiusAlong.append(radius) # Warp points @@ -907,11 +906,11 @@ def generateBaseMesh(cls, region, options): # Arclength between apex point and corresponding point on next face mag = interp.getCubicHermiteArcLength(xList[0], d2List[0], xList[2 * elementsCountAround], d2List[2 * elementsCountAround]) - d2ApexInner = vector.setMagnitude(sd2[0], mag) - d1ApexInner = vector.crossproduct3(sd1[0], d2ApexInner) - d1ApexInner = vector.setMagnitude(d1ApexInner, mag) - d3ApexUnit = vector.normalise( - vector.crossproduct3(vector.normalise(d1ApexInner), vector.normalise(d2ApexInner))) + d2ApexInner = set_magnitude(sd2[0], mag) + d1ApexInner = cross(sd1[0], d2ApexInner) + d1ApexInner = set_magnitude(d1ApexInner, mag) + d3ApexUnit = normalize( + cross(normalize(d1ApexInner), normalize(d2ApexInner))) d3ApexInner = [d3ApexUnit[c] * bladderWallThickness / elementsCountThroughWall for c in range(3)] xFinal = [] @@ -1390,7 +1389,7 @@ def generateUreterInlets(region, nodes, mesh, ureterDefaultOptions, elementsCoun ureterElementPositionAround + (1 if ureter1Position.xi1 > 0.5 else 0) ureter1StartCornerx = xBladder[endPointStartId1 - 1] v1 = [(ureter1StartCornerx[c] - centerUreter1_x[c]) for c in range(3)] - ureter1Direction = vector.crossproduct3(td3, v1) + ureter1Direction = cross(td3, v1) nodeIdentifier, elementIdentifier, (o1_x, o1_d1, o1_d2, _, o1_NodeId, o1_Positions) = \ generateOstiumMesh(region, ureterDefaultOptions, trackSurfaceUreter1, ureter1Position, ureter1Direction, startNodeIdentifier=nextNodeIdentifier, startElementIdentifier=nextElementIdentifier, @@ -1406,7 +1405,7 @@ def generateUreterInlets(region, nodes, mesh, ureterDefaultOptions, elementsCoun elementsCountAround - ureterElementPositionAround + (-1 if ureter1Position.xi1 > 0.5 else 0) ureter2StartCornerx = xBladder[endPointStartId2 - 1] v2 = [(ureter2StartCornerx[c] - centerUreter2_x[c]) for c in range(3)] - ureter2Direction = vector.crossproduct3(ad3, v2) + ureter2Direction = cross(ad3, v2) nodeIdentifier, elementIdentifier, (o2_x, o2_d1, o2_d2, _, o2_NodeId, o2_Positions) = \ generateOstiumMesh(region, ureterDefaultOptions, trackSurfaceUreter2, ureter2Position, ureter2Direction, startNodeIdentifier=nodeIdentifier, startElementIdentifier=elementIdentifier, @@ -1594,7 +1593,7 @@ def obtainBladderFlatNodes(elementsCountAlongBladder, elementsCountAround, eleme # Find the angle at the bottom of the bladder neck v1 = [0.0, 0.0, bladderLength] v2 = [0.5 * neckDiameter1, 0.0, bladderLength] - alpha = vector.angleBetweenVectors(v1, v2) + alpha = angle(v1, v2) # Find apex to urethra arcLength in minor radius minorNodeAlong_x.insert(0, xApexInner) @@ -1708,8 +1707,8 @@ def obtainBladderFlatNodes(elementsCountAlongBladder, elementsCountAround, eleme v3 = xfnListRearranged[elementsCountAround // 2] v21 = [v2[c] - v1[c] for c in range(3)] v31 = [v3[c] - v1[c] for c in range(3)] - d1Mag = vector.magnitude(v21) - d2Mag = vector.magnitude(v31) + d1Mag = magnitude(v21) + d2Mag = magnitude(v31) # Add apex nodes to the list xFlat = [] diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_bone1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_bone1.py index 764cf379..3eff72df 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_bone1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_bone1.py @@ -7,6 +7,7 @@ import copy +from cmlibs.maths.vectorops import add_vectors from cmlibs.utils.zinc.field import findOrCreateFieldCoordinates from cmlibs.utils.zinc.finiteelement import getMaximumNodeIdentifier, getMaximumElementIdentifier from cmlibs.zinc.field import Field @@ -14,14 +15,11 @@ from scaffoldmaker.meshtypes.meshtype_1d_path1 import MeshType_1d_path1 from scaffoldmaker.meshtypes.scaffold_base import Scaffold_base from scaffoldmaker.scaffoldpackage import ScaffoldPackage -from scaffoldmaker.utils import vector from scaffoldmaker.utils.cylindermesh import CylinderMesh, CylinderShape, CylinderEnds, CylinderCentralPath from scaffoldmaker.utils.meshrefinement import MeshRefinement from scaffoldmaker.utils.shieldmesh import ShieldMesh3D from scaffoldmaker.utils.spheremesh import SphereMesh, SphereShape from scaffoldmaker.utils.zinc_utils import exnode_string_from_nodeset_field_parameters -from scaffoldmaker.utils.derivativemoothing import DerivativeSmoothing -from scaffoldmaker.annotation.annotationgroup import AnnotationGroup, findOrCreateAnnotationGroupForTerm class MeshType_3d_bone1 (Scaffold_base): @@ -202,7 +200,7 @@ def generateBaseMesh(region, options): sphere_centre = sphere_base.centre sphere_radius_3 = options['Lower scale_Z'] axes = [sphere_base.majorAxis, sphere_base.minorAxis, - vector.setMagnitude(vector.crossproduct3(sphere_base.majorAxis, sphere_base.minorAxis), + set_magnitude(cross(sphere_base.majorAxis, sphere_base.minorAxis), sphere_radius_3)] elementsCountAcross = [cylinder1._elementsCountAcrossMajor, cylinder1._elementsCountAcrossMinor, cylinder1._elementsCountAcrossMajor] @@ -235,8 +233,7 @@ def generateBaseMesh(region, options): for n1 in range(elementsCountAcross[1] + 1): if n3 < elementsCountAcross[2] // 2: if sphere1._shield3D.px[n3][n2][n1]: - hemisphere.px[n3][n2][n1] = vector.addVectors([sphere1._shield3D.px[n3][n2][n1], - sphere_centre], [1, 1]) + hemisphere.px[n3][n2][n1] = add_vectors([sphere1._shield3D.px[n3][n2][n1], sphere_centre], [1, 1]) #cylinder end elif n3 == elementsCountAcross[2] // 2: # find nodes on the triple line. Note that cylinder and sphere have a little bit different @@ -292,7 +289,7 @@ def generateBaseMesh(region, options): sphere_centre = sphere_base.centre sphere_radius_3 = options['Upper scale_Z'] axes = [sphere_base.majorAxis, sphere_base.minorAxis, - vector.setMagnitude(vector.crossproduct3(sphere_base.majorAxis, sphere_base.minorAxis), + set_magnitude(cross(sphere_base.majorAxis, sphere_base.minorAxis), sphere_radius_3)] elementsCountAcross = [cylinder1._elementsCountAcrossMajor, cylinder1._elementsCountAcrossMinor, cylinder1._elementsCountAcrossMajor] @@ -326,8 +323,7 @@ def generateBaseMesh(region, options): for n1 in range(elementsCountAcross[1] + 1): if n3 > elementsCountAcross[2] // 2: if sphere2._shield3D.px[n3][n2][n1]: - hemisphere2.px[n3][n2][n1] = vector.addVectors([sphere2._shield3D.px[n3][n2][n1], - sphere_centre], [1, 1]) + hemisphere2.px[n3][n2][n1] = add_vectors([sphere2._shield3D.px[n3][n2][n1], sphere_centre], [1, 1]) #cylinder end elif n3 == elementsCountAcross[2] // 2: # find nodes on the triple line. Note that cylinder and sphere have a little bit different diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_cecum1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_cecum1.py index fba54ac9..e1448962 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_cecum1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_cecum1.py @@ -23,7 +23,6 @@ from scaffoldmaker.utils import interpolation as interp from scaffoldmaker.utils import matrix from scaffoldmaker.utils import tubemesh -from scaffoldmaker.utils import vector from scaffoldmaker.utils.annulusmesh import createAnnulusMesh3d from scaffoldmaker.utils.eftfactory_bicubichermitelinear import eftfactory_bicubichermitelinear from scaffoldmaker.utils.eftfactory_tricubichermite import eftfactory_tricubichermite @@ -410,7 +409,7 @@ def getApexSegmentForCecum(xOuter, d1Outer, d2Outer, elementsCountAroundHalfHaus # Compile nodes and d2 for sampling xFirstSegment += xOuter[elementsCountAround * int(elementsCountAlongSegment * 0.5):] # second half of first regular segment - d1FirstDirectionVector = vector.normalise(d1Outer[elementsCountAround]) # Store direction vector of first d1 intra-haustral for later + d1FirstDirectionVector = normalize(d1Outer[elementsCountAround]) # Store direction vector of first d1 intra-haustral for later d2Vector = xOuter[elementsCountAround * int(elementsCountAlongSegment * 0.5): elementsCountAround * (int(elementsCountAlongSegment * 0.5) + 1)] # half face of segment - apex d2FirstSegment = [] @@ -462,20 +461,20 @@ def getApexSegmentForCecum(xOuter, d1Outer, d2Outer, elementsCountAroundHalfHaus d1 = d1FirstDirectionVector if n1 == 0 else [v2[c] - v1[c] for c in range(3)] d2 = [v2[c] - v1[c] for c in range(3)] arcLengthAround = interp.computeCubicHermiteArcLength(v1, d1, v2, d2, True) - dx_ds1 = [c*arcLengthAround for c in vector.normalise(d1)] + dx_ds1 = [c*arcLengthAround for c in normalize(d1)] d1Around.append(dx_ds1) # Account for d1 of node sitting on half haustrum - d1 = vector.normalise( + d1 = normalize( [xAround[elementsCountAroundHalfHaustrum][c] - xAround[elementsCountAroundHalfHaustrum - 1][c] for c in range(3)]) dx_ds1 = [c * arcLengthAround for c in d1] d1Around.append(dx_ds1) d1Smoothed = interp.smoothCubicHermiteDerivativesLine(xAround, d1Around, fixStartDerivative=True) - d1TCEdge = vector.setMagnitude(d1Smoothed[int(elementsCountAroundTC * 0.5)], - vector.magnitude(d1Smoothed[int(elementsCountAroundTC * 0.5 - 1)])) - d1Transition = vector.setMagnitude(d1Smoothed[int(elementsCountAroundTC * 0.5 + 1)], - vector.magnitude(d1Smoothed[int(elementsCountAroundTC * 0.5 + 2)])) + d1TCEdge = set_magnitude(d1Smoothed[int(elementsCountAroundTC * 0.5)], + magnitude(d1Smoothed[int(elementsCountAroundTC * 0.5 - 1)])) + d1Transition = set_magnitude(d1Smoothed[int(elementsCountAroundTC * 0.5 + 1)], + magnitude(d1Smoothed[int(elementsCountAroundTC * 0.5 + 2)])) d1Corrected = [] d1Corrected = d1Corrected + d1Smoothed[:int(elementsCountAroundTC * 0.5)] d1Corrected.append(d1TCEdge) @@ -580,20 +579,20 @@ def sampleCecumAlongLength(xToSample, d1ToSample, d2ToSample, d1FirstDirectionVe d1 = d1FirstDirectionVector if n1 == 0 else [v2[c] - v1[c] for c in range(3)] d2 = [v2[c] - v1[c] for c in range(3)] arcLengthAround = interp.computeCubicHermiteArcLength(v1, d1, v2, d2, True) - dx_ds1 = [c * arcLengthAround for c in vector.normalise(d1)] + dx_ds1 = [c * arcLengthAround for c in normalize(d1)] d1OuterAroundList.append(dx_ds1) # Account for d1 of node sitting on half haustrum - d1 = vector.normalise([xAround[elementsCountAroundHalfHaustrum][c] - + d1 = normalize([xAround[elementsCountAroundHalfHaustrum][c] - xAround[elementsCountAroundHalfHaustrum - 1][c] for c in range(3)]) dx_ds1 = [c * arcLengthAround for c in d1] d1OuterAroundList.append(dx_ds1) if d1OuterAroundList: d1Smoothed = interp.smoothCubicHermiteDerivativesLine(xAround, d1OuterAroundList, fixStartDerivative=True) - d1TCEdge = vector.setMagnitude(d1Smoothed[int(elementsCountAroundTC * 0.5)], - vector.magnitude(d1Smoothed[int(elementsCountAroundTC * 0.5 - 1)])) - d1Transition = vector.setMagnitude(d1Smoothed[int(elementsCountAroundTC * 0.5 + 1)], - vector.magnitude(d1Smoothed[int(elementsCountAroundTC * 0.5 + 2)])) + d1TCEdge = set_magnitude(d1Smoothed[int(elementsCountAroundTC * 0.5)], + magnitude(d1Smoothed[int(elementsCountAroundTC * 0.5 - 1)])) + d1Transition = set_magnitude(d1Smoothed[int(elementsCountAroundTC * 0.5 + 1)], + magnitude(d1Smoothed[int(elementsCountAroundTC * 0.5 + 2)])) d1Corrected = [] d1Corrected = d1Corrected + d1Smoothed[:int(elementsCountAroundTC * 0.5)] d1Corrected.append(d1TCEdge) @@ -618,7 +617,7 @@ def findDerivativeBetweenPoints(v1, v2): """ d = [v2[c] - v1[c] for c in range(3)] arcLengthAround = interp.computeCubicHermiteArcLength(v1, d, v2, d, True) - d = [c * arcLengthAround for c in vector.normalise(d)] + d = [c * arcLengthAround for c in normalize(d)] return d @@ -736,7 +735,7 @@ def createCecumMesh3d(region, options, networkLayout, nodeIdentifier, elementIde tcWidthAlongCecum = [] closedProximalEnd = True - outerRadiusListCP = [vector.magnitude(c) for c in cd2] + outerRadiusListCP = [magnitude(c) for c in cd2] dOuterRadiusListCP = [] for n in range(len(outerRadiusListCP) - 1): dOuterRadiusListCP.append(outerRadiusListCP[n + 1] - outerRadiusListCP[n]) @@ -848,11 +847,11 @@ def createCecumMesh3d(region, options, networkLayout, nodeIdentifier, elementIde elementsCountAround * (elementsCountThroughWall + 1)], rescaleDerivatives=True) mag = 0.5*(magMin + magMax) - d2ApexInner = vector.setMagnitude(sd2Cecum[0], mag) - d1ApexInner = vector.crossproduct3(sd1Cecum[0], d2ApexInner) - d1ApexInner = vector.setMagnitude(d1ApexInner, mag) - d3ApexUnit = vector.normalise(vector.crossproduct3(vector.normalise(d1ApexInner), - vector.normalise(d2ApexInner))) + d2ApexInner = set_magnitude(sd2Cecum[0], mag) + d1ApexInner = cross(sd1Cecum[0], d2ApexInner) + d1ApexInner = set_magnitude(d1ApexInner, mag) + d3ApexUnit = normalize(cross(normalize(d1ApexInner), + normalize(d2ApexInner))) d3ApexInner = [d3ApexUnit[c] * wallThickness / elementsCountThroughWall for c in range(3)] xCecum = [] @@ -952,10 +951,10 @@ def createCecumMesh3d(region, options, networkLayout, nodeIdentifier, elementIde # Find region where ostium sits # angle between d2 of branch point and vector between branch point and 1st point on ileum dV = [cxIleum[0][c] - cxIleum[-1][c] for c in range(3)] - ostiumPositionAngleAround = math.acos(vector.dotproduct(dV, d2BranchPt)/ - (vector.magnitude(dV) * vector.magnitude(d2BranchPt))) - if vector.dotproduct(dV,d3BranchPt) != 0: - sign = vector.dotproduct(dV, d3BranchPt)/abs(vector.dotproduct(dV, d3BranchPt)) + ostiumPositionAngleAround = math.acos(dot(dV, d2BranchPt)/ + (magnitude(dV) * magnitude(d2BranchPt))) + if dot(dV,d3BranchPt) != 0: + sign = dot(dV, d3BranchPt)/abs(dot(dV, d3BranchPt)) if sign < 0: ostiumPositionAngleAround = math.pi * 2.0 - ostiumPositionAngleAround sectorIdx = ostiumPositionAngleAround // (2 * math.pi / tcCount) @@ -1067,17 +1066,17 @@ def createCecumMesh3d(region, options, networkLayout, nodeIdentifier, elementIde arcEnd = networkLayout.arcLengthOfGroupsAlong[1] nearestPosition = trackSurfaceOstium.findNearestPosition(cxIleum[0]) xNearestStart = trackSurfaceOstium.evaluateCoordinates(nearestPosition, derivatives=False) - distStart = vector.magnitude([cxIleum[0][c] - xNearestStart[c] for c in range(3)]) + distStart = magnitude([cxIleum[0][c] - xNearestStart[c] for c in range(3)]) nearestPosition = trackSurfaceOstium.findNearestPosition(cxIleum[-1]) xNearestEnd = trackSurfaceOstium.evaluateCoordinates(nearestPosition, derivatives=False) - distEnd = vector.magnitude([cxIleum[-1][c] - xNearestEnd[c] for c in range(3)]) + distEnd = magnitude([cxIleum[-1][c] - xNearestEnd[c] for c in range(3)]) for iter in range(100): arcDistance = (arcStart + arcEnd) * 0.5 x, d1 = interp.getCubicHermiteCurvesPointAtArcDistance(cxIleum, cd1Ileum, arcDistance)[0:2] nearestPosition = trackSurfaceOstium.findNearestPosition(x) xNearest = trackSurfaceOstium.evaluateCoordinates(nearestPosition, derivatives=False) - dist = vector.magnitude([x[c] - xNearest[c] for c in range(3)]) + dist = magnitude([x[c] - xNearest[c] for c in range(3)]) if abs(distStart - distEnd) > xTol: if distStart < distEnd: @@ -1089,7 +1088,7 @@ def createCecumMesh3d(region, options, networkLayout, nodeIdentifier, elementIde else: xCentre, d1Centre, d2Centre = trackSurfaceOstium.evaluateCoordinates(nearestPosition, derivatives=True) - normAxis = vector.normalise([-d for d in d1]) + normAxis = normalize([-d for d in d1]) eIdx = interp.getNearestPointIndex(cxIleum, xCentre) - 1 arcLenghtSum = 0.0 for e in range(eIdx): @@ -1164,11 +1163,11 @@ def createCecumMesh3d(region, options, networkLayout, nodeIdentifier, elementIde e1Right = 0 e2Top = 0 e2Bottom = elementsAlongTrackSurface - sf = vector.magnitude(networkLayoutIleum.cd2Path[-1]) * 0.35 + sf = magnitude(networkLayoutIleum.cd2Path[-1]) * 0.35 for n1 in range(elementsCountAroundOstium): - normD2 = vector.normalise(o1_d2[-1][n1]) + normD2 = normalize(o1_d2[-1][n1]) d1AnnulusNorm.append(normD2) - d1AnnulusOuter.append(vector.setMagnitude(o1_d2[-1][n1], sf)) + d1AnnulusOuter.append(set_magnitude(o1_d2[-1][n1], sf)) x = [o1_x[-1][n1][c] + sf * normD2[c] for c in range(3)] nearestPosition = trackSurfaceOstium.findNearestPosition(x) e1 = nearestPosition.e1 @@ -1191,7 +1190,7 @@ def createCecumMesh3d(region, options, networkLayout, nodeIdentifier, elementIde d2AnnulusOuter = interp.smoothCubicHermiteDerivativesLoop(xAnnulusOuter, d2AnnulusOuter) d3Annulus = [] for n in range(elementsCountAroundOstium): - d3 = vector.normalise(vector.crossproduct3(vector.normalise(d2AnnulusOuter[n]), d1AnnulusNorm[n])) + d3 = normalize(cross(normalize(d2AnnulusOuter[n]), d1AnnulusNorm[n])) d3Annulus.append(d3) annulusD2Curvature = interp.getCurvaturesAlongCurve(xAnnulusOuter, d2AnnulusOuter, d3Annulus, loop=True) @@ -1222,7 +1221,7 @@ def createCecumMesh3d(region, options, networkLayout, nodeIdentifier, elementIde rotFrame[j][2] * d1Cecum[1][2] for j in range(3)] rotAngle = math.pi - rotFrame = matrix.getRotationMatrixFromAxisAngle(vector.normalise(d3Annulus[0]), rotAngle) + rotFrame = matrix.getRotationMatrixFromAxisAngle(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)] @@ -1236,8 +1235,8 @@ def createCecumMesh3d(region, options, networkLayout, nodeIdentifier, elementIde xPositionA = trackSurfaceOstium.findNearestPosition(xStart) xProportionA = trackSurfaceOstium.getProportion(xPositionA) derivativeA = trackSurfaceOstium.evaluateCoordinates(xPositionA, derivatives=True)[2] - derivativeA = vector.setMagnitude(derivativeA, vector.magnitude(d2ApexInner)) - derivativeMagnitudeA = vector.magnitude(derivativeA) + derivativeA = set_magnitude(derivativeA, magnitude(d2ApexInner)) + derivativeMagnitudeA = magnitude(derivativeA) xPositionB = trackSurfaceOstium.findNearestPosition(xAnnulusOuter[0]) xProportionB = trackSurfaceOstium.getProportion(xPositionB) @@ -1251,7 +1250,7 @@ def createCecumMesh3d(region, options, networkLayout, nodeIdentifier, elementIde trackSurfaceOstium.resampleHermiteCurvePointsSmooth( nx, nd1, nd2, nd3, proportions, derivativeMagnitudeStart=derivativeMagnitudeA)[0:3] - annulusD2ZeroMag = vector.magnitude(pd2AlongMidLine[-1]) + annulusD2ZeroMag = magnitude(pd2AlongMidLine[-1]) for n in range(len(pd1AlongMidLine)): pd1AlongMidLine[n] = [-d for d in pd1AlongMidLine[n]] @@ -1263,7 +1262,7 @@ def createCecumMesh3d(region, options, networkLayout, nodeIdentifier, elementIde xPositionB = trackSurfaceOstium.findNearestPosition(xB) xProportionB = trackSurfaceOstium.getProportion(xPositionB) derivativeB = d2TrackSurface[-elementsCountAroundHalfHaustrum] - derivativeMagnitudeB = vector.magnitude(d2TrackSurface[-elementsCountAroundHalfHaustrum]) + derivativeMagnitudeB = magnitude(d2TrackSurface[-elementsCountAroundHalfHaustrum]) nx, nd1, nd2, nd3, proportions = \ trackSurfaceOstium.createHermiteCurvePoints( @@ -1275,7 +1274,7 @@ def createCecumMesh3d(region, options, networkLayout, nodeIdentifier, elementIde nx, nd1, nd2, nd3, proportions, derivativeMagnitudeStart=None, derivativeMagnitudeEnd=derivativeMagnitudeB)[0:3] - annulusD2HalfOstiumMag = vector.magnitude(pd2AlongMidLineBottom[0]) + annulusD2HalfOstiumMag = magnitude(pd2AlongMidLineBottom[0]) for n in range(len(pd1AlongMidLineBottom)): pd1AlongMidLineBottom[n] = [-d for d in pd1AlongMidLineBottom[n]] @@ -1301,7 +1300,7 @@ def createCecumMesh3d(region, options, networkLayout, nodeIdentifier, elementIde xTrackSurface[n2 * (elementsCountAroundHaustrum + 1)]) xProportionA = trackSurfaceOstium.getProportion(xPositionA) derivativeA = d1TrackSurface[n2 * (elementsCountAroundHaustrum + 1)] - derivativeMagnitudeA = vector.magnitude(derivativeA) + derivativeMagnitudeA = magnitude(derivativeA) if n2 < rowsBelow + 1: xPositionB = trackSurfaceOstium.findNearestPosition(pxAlongMidLine[n2]) @@ -1343,7 +1342,7 @@ def createCecumMesh3d(region, options, networkLayout, nodeIdentifier, elementIde xTrackSurface[n2 * (elementsCountAroundHaustrum + 1) + elementsCountAroundHaustrum]) xProportionB = trackSurfaceOstium.getProportion(xPositionB) derivativeB = d1TrackSurface[n2 * (elementsCountAroundHaustrum + 1) + elementsCountAroundHaustrum] - derivativeMagnitudeB = vector.magnitude(derivativeB) + derivativeMagnitudeB = magnitude(derivativeB) nx, nd1, nd2, nd3, proportions = \ trackSurfaceOstium.createHermiteCurvePoints( @@ -1422,13 +1421,13 @@ def createCecumMesh3d(region, options, networkLayout, nodeIdentifier, elementIde nd2AlongBottomLHS = interp.smoothCubicHermiteDerivativesLine(nxAlong, nd2Along, fixStartDerivative=True) # Make sure d2 at annulus is length of element below - nd2Along[-1] = vector.setMagnitude(d1AnnulusOuter[1], vector.magnitude(nd2Along[-1])) + nd2Along[-1] = set_magnitude(d1AnnulusOuter[1], magnitude(nd2Along[-1])) # Make magnitude of annulus d2 between start and end row equivalent to arclength of element on its left for m in range(endRowIdx - startRowIdx - 1): nxAlong.append(xAnnulusOuter[2 + m]) n2Idx = m + startRowIdx + 1 - nd2Along.append(vector.setMagnitude(d1AnnulusOuter[2+m], vector.magnitude(d1AroundAlong[n2Idx][n1]))) + nd2Along.append(set_magnitude(d1AnnulusOuter[2+m], magnitude(d1AroundAlong[n2Idx][n1]))) # Smooth from annulus to end of cecum for n2 in range(endRowIdx, elementsCountAlongSegment): @@ -1440,8 +1439,8 @@ def createCecumMesh3d(region, options, networkLayout, nodeIdentifier, elementIde nd2AlongTopLHS = interp.smoothCubicHermiteDerivativesLine(nxTop, nd2Top, fixEndDerivative=True) # Make sure d2 at annulus is length of element above - nd2Top[0] = vector.setMagnitude(d1AnnulusOuter[int(elementsCountAroundOstium * 0.5) - 1], - vector.magnitude(nd2Top[0])) + nd2Top[0] = set_magnitude(d1AnnulusOuter[int(elementsCountAroundOstium * 0.5) - 1], + magnitude(nd2Top[0])) nxAlong += nxTop nd2Along += nd2Top @@ -1473,14 +1472,14 @@ def createCecumMesh3d(region, options, networkLayout, nodeIdentifier, elementIde nd2AlongBottomRHS = interp.smoothCubicHermiteDerivativesLine(nxAlong, nd2Along, fixStartDerivative=True) # Make sure d2 at annulus is length of element below - nd2Along[-1] = vector.setMagnitude(d1AnnulusOuter[-1], vector.magnitude(nd2Along[-1])) + nd2Along[-1] = set_magnitude(d1AnnulusOuter[-1], magnitude(nd2Along[-1])) # Make magnitude of annulus d2 between start and end row equivalent to arclength of element on its right for m in range(endRowIdx - startRowIdx - 1): nxAlong.append(xAnnulusOuter[-(2 + m)]) n2Idx = m + startRowIdx + 1 - nd2Along.append(vector.setMagnitude(d1AnnulusOuter[-(2 + m)], - vector.magnitude(d1AroundAlong[n2Idx][n1 - 1]))) + nd2Along.append(set_magnitude(d1AnnulusOuter[-(2 + m)], + magnitude(d1AroundAlong[n2Idx][n1 - 1]))) # Smooth from annulus to end of cecum for n2 in range(endRowIdx, elementsCountAlongSegment): @@ -1497,8 +1496,8 @@ def createCecumMesh3d(region, options, networkLayout, nodeIdentifier, elementIde nd2AlongTopRHS = interp.smoothCubicHermiteDerivativesLine(nxTop, nd2Top, fixEndDerivative=True) # Make sure d2 at annulus is length of element above - nd2Top[0] = vector.setMagnitude(d1AnnulusOuter[int(elementsCountAroundOstium * 0.5) + 1], - vector.magnitude(nd2Top[0])) + nd2Top[0] = set_magnitude(d1AnnulusOuter[int(elementsCountAroundOstium * 0.5) + 1], + magnitude(nd2Top[0])) nxAlong += nxTop nd2Along += nd2Top @@ -1546,8 +1545,8 @@ def createCecumMesh3d(region, options, networkLayout, nodeIdentifier, elementIde for n2 in range(len(xAroundAlong)): d3Around = [] for n1 in range(len(xAroundAlong[n2])): - d3Around.append(vector.normalise( - vector.crossproduct3(vector.normalise(d1AroundAlong[n2][n1]), vector.normalise(d2AroundAlong[n2][n1])))) + d3Around.append(normalize( + cross(normalize(d1AroundAlong[n2][n1]), normalize(d2AroundAlong[n2][n1])))) d3UnitAroundAlong.append(d3Around) # Calculate curvatures @@ -1680,7 +1679,7 @@ def createCecumMesh3d(region, options, networkLayout, nodeIdentifier, elementIde # Adjust annulus points xAroundAlong[startRowIdx].insert(int(elementsCountAroundHaustrum * 0.5), xAnnulusOuter[0]) d1AroundAlong[startRowIdx].insert(int(elementsCountAroundHaustrum * 0.5), d2AnnulusOuter[0]) - d2AroundAlong[startRowIdx].insert(int(elementsCountAroundHaustrum * 0.5), vector.setMagnitude(d1AnnulusOuter[0], + d2AroundAlong[startRowIdx].insert(int(elementsCountAroundHaustrum * 0.5), set_magnitude(d1AnnulusOuter[0], annulusD2ZeroMag)) d3UnitAroundAlong[startRowIdx].insert(int(elementsCountAroundHaustrum * 0.5), d3Annulus[0]) d1Curvature[startRowIdx].insert(int(elementsCountAroundHaustrum * 0.5), annulusD2Curvature[0]) @@ -1689,7 +1688,7 @@ def createCecumMesh3d(region, options, networkLayout, nodeIdentifier, elementIde idx = int(elementsCountAlongSegment * 0.5) xAroundAlong[endRowIdx].insert(int(elementsCountAroundHaustrum * 0.5), xAnnulusOuter[idx]) d1AroundAlong[endRowIdx].insert(int(elementsCountAroundHaustrum * 0.5), d2AnnulusOuter[idx]) - d2AroundAlong[endRowIdx].insert(int(elementsCountAroundHaustrum * 0.5), vector.setMagnitude(d1AnnulusOuter[idx], + d2AroundAlong[endRowIdx].insert(int(elementsCountAroundHaustrum * 0.5), set_magnitude(d1AnnulusOuter[idx], annulusD2HalfOstiumMag)) d3UnitAroundAlong[endRowIdx].insert(int(elementsCountAroundHaustrum * 0.5), d3Annulus[idx]) d1Curvature[endRowIdx].insert(int(elementsCountAroundHaustrum * 0.5), annulusD2Curvature[idx]) diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py index a5a9070b..3d714fa7 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py @@ -20,7 +20,6 @@ from scaffoldmaker.scaffoldpackage import ScaffoldPackage from scaffoldmaker.utils import interpolation as interp from scaffoldmaker.utils import tubemesh -from scaffoldmaker.utils import vector from scaffoldmaker.utils.zinc_utils import exnode_string_from_nodeset_field_parameters, \ get_nodeset_path_field_parameters @@ -830,7 +829,7 @@ def createColonMesh3d(region, options, networkLayout, nextNodeIdentifier, nextEl arcLengthOfGroupsAlong[1] + arcLengthOfGroupsAlong[2], arcLengthOfGroupsAlong[0]] - outerRadiusListCP = [vector.magnitude(c) for c in cd2] + outerRadiusListCP = [magnitude(c) for c in cd2] dOuterRadiusListCP = [] for n in range(len(outerRadiusListCP) - 1): dOuterRadiusListCP.append(outerRadiusListCP[n + 1] - outerRadiusListCP[n]) diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_colonsegment1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_colonsegment1.py index 24688dd4..6d90bc30 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_colonsegment1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_colonsegment1.py @@ -20,7 +20,6 @@ from scaffoldmaker.utils import interpolation as interp from scaffoldmaker.utils import matrix from scaffoldmaker.utils import tubemesh -from scaffoldmaker.utils import vector from scaffoldmaker.utils.eftfactory_bicubichermitelinear import eftfactory_bicubichermitelinear from scaffoldmaker.utils.eftfactory_tricubichermite import eftfactory_tricubichermite from scaffoldmaker.utils.geometry import createCirclePoints @@ -805,10 +804,10 @@ def getColonSegmentOuterPoints(region, elementsCountAroundTC, elementsCountAroun d1 = d1AtStartOfEachMidFace[n2] if n1 == 0 else [v2[c] - v1[c] for c in range(3)] d2 = [v2[c] - v1[c] for c in range(3)] arcLengthAround = interp.computeCubicHermiteArcLength(v1, d1, v2, d2, True) - dx_ds1 = [c * arcLengthAround for c in vector.normalise(d1)] + dx_ds1 = [c * arcLengthAround for c in normalize(d1)] dx_ds1OuterAroundList.append(dx_ds1) # Account for d1 of node sitting on half haustrum - d1 = vector.normalise([xAround[elementsCountAroundHalfHaustrum][c] - + d1 = normalize([xAround[elementsCountAroundHalfHaustrum][c] - xAround[elementsCountAroundHalfHaustrum - 1][c] for c in range(3)]) dx_ds1 = [c * arcLengthAround for c in d1] dx_ds1OuterAroundList.append(dx_ds1) @@ -816,10 +815,10 @@ def getColonSegmentOuterPoints(region, elementsCountAroundTC, elementsCountAroun if dx_ds1OuterAroundList: d1Smoothed = interp.smoothCubicHermiteDerivativesLine(xAround, dx_ds1OuterAroundList, fixStartDerivative=True) - d1TCEdge = vector.setMagnitude(d1Smoothed[int(elementsCountAroundTC * 0.5)], - vector.magnitude(d1Smoothed[int(elementsCountAroundTC * 0.5 - 1)])) - d1Transition = vector.setMagnitude(d1Smoothed[int(elementsCountAroundTC * 0.5 + 1)], - vector.magnitude(d1Smoothed[int(elementsCountAroundTC * 0.5 + 2)])) + d1TCEdge = set_magnitude(d1Smoothed[int(elementsCountAroundTC * 0.5)], + magnitude(d1Smoothed[int(elementsCountAroundTC * 0.5 - 1)])) + d1Transition = set_magnitude(d1Smoothed[int(elementsCountAroundTC * 0.5 + 1)], + magnitude(d1Smoothed[int(elementsCountAroundTC * 0.5 + 2)])) d1Corrected = [] d1Corrected = d1Corrected + d1Smoothed[:int(elementsCountAroundTC * 0.5)] d1Corrected.append(d1TCEdge) @@ -1094,7 +1093,7 @@ def sampleTeniaColi(nx, nd1, arcDistanceTCEdge, elementsCountAroundTC): for e in range(int(elementsCountAroundTC * 0.5) + 1): arcDistance = arcDistancePerElementTC * e x, d1, _, _ = interp.getCubicHermiteCurvesPointAtArcDistance(nx, nd1, arcDistance) - d1Scaled = vector.setMagnitude(d1, arcDistancePerElementTC) + d1Scaled = set_magnitude(d1, arcDistancePerElementTC) xTC.append(x) d1TC.append(d1Scaled) @@ -1122,7 +1121,7 @@ def sampleHaustrum(nx, nd1, xTCLast, d1TCLast, arcLength, arcDistanceTCEdge, elementLengths = [] length = arcLength - arcDistanceTCEdge elementsCountOut = int(elementsCountAroundHaustrum * 0.5) - addLengthStart = 0.5 * vector.magnitude(d1TCLast) + addLengthStart = 0.5 * magnitude(d1TCLast) lengthFractionStart = 0.5 proportionStart = 1.0 elementLengthMid = (length - addLengthStart) / (elementsCountOut - 1.0 + proportionStart * lengthFractionStart) @@ -1138,13 +1137,13 @@ def sampleHaustrum(nx, nd1, xTCLast, d1TCLast, arcLength, arcDistanceTCEdge, arcDistance = arcDistanceTCEdge xHaustrum.append(xTCLast) - d1Scaled = vector.setMagnitude(d1TCLast, elementLengths[0]) + d1Scaled = set_magnitude(d1TCLast, elementLengths[0]) d1Haustrum.append(d1Scaled) for e in range(elementsCountOut): arcDistance = arcDistance + elementLengths[e] x, d1, _, _ = interp.getCubicHermiteCurvesPointAtArcDistance(nx, nd1, arcDistance) - d1Scaled = vector.setMagnitude(d1, elementLengths[e] if e > 0 else elementLengths[e + 1]) + d1Scaled = set_magnitude(d1, elementLengths[e] if e > 0 else elementLengths[e + 1]) xHaustrum.append(x) d1Haustrum.append(d1Scaled) @@ -1348,7 +1347,7 @@ def getTeniaColi(region, xList, d1List, d2List, d3List, curvatureList, else n2 * elementsCountAround * (elementsCountThroughWall + 1) + \ elementsCountAround * elementsCountThroughWall + \ N * (elementsCountAroundTC + elementsCountAroundHaustrum) - unitNorm = vector.normalise(d3List[idxTCMid]) + unitNorm = normalize(d3List[idxTCMid]) xMid = [xList[idxTCMid][i] + unitNorm[i] * tcThickness for i in range(3)] d1Mid = d1List[idxTCMid] TCStartIdx = idxTCMid - int(elementsCountAroundTC * 0.5) if N > 0 else \ @@ -1356,11 +1355,11 @@ def getTeniaColi(region, xList, d1List, d2List, d3List, curvatureList, int(elementsCountAroundTC * 0.5) TCEndIdx = idxTCMid + int(elementsCountAroundTC * 0.5) if closedProximalEnd: - tcWidth = vector.magnitude([xList[TCStartIdx][c] - xList[TCEndIdx][c] for c in range(3)]) + tcWidth = magnitude([xList[TCStartIdx][c] - xList[TCEndIdx][c] for c in range(3)]) v1 = xList[TCStartIdx] v2 = xMid - d1MidScaled = [c * tcWidth * TCEdgeFactor for c in vector.normalise(d1Mid)] + d1MidScaled = [c * tcWidth * TCEdgeFactor for c in normalize(d1Mid)] v3 = xList[TCEndIdx] nx = [v1, v2, v3] nd1 = [d1List[TCStartIdx], d1MidScaled, d1List[TCEndIdx]] @@ -1368,8 +1367,8 @@ def getTeniaColi(region, xList, d1List, d2List, d3List, curvatureList, xTCRaw = xTCRaw + sxTC[1:-1] if elementsCountAroundTC == 2: p = [v2[i] - v1[i] for i in range(3)] - A = vector.dotproduct(unitNorm, p) # A<0 if v2 is higher than v1 - d1 = [c * tcWidth * 0.5 for c in vector.normalise(d1Mid)] if A < 0 else d1MidScaled + A = dot(unitNorm, p) # A<0 if v2 is higher than v1 + d1 = [c * tcWidth * 0.5 for c in normalize(d1Mid)] if A < 0 else d1MidScaled d1TCRaw = d1TCRaw + sd1TC[1:-1] if elementsCountAroundTC > 2 else d1TCRaw + [d1] xTCInnerSet = list(range(TCStartIdx + 1, TCEndIdx)) if N > 0 else \ list(range(TCStartIdx + 1, TCStartIdx + int(elementsCountAroundTC * 0.5))) + \ @@ -1383,7 +1382,7 @@ def getTeniaColi(region, xList, d1List, d2List, d3List, curvatureList, d3List[xTCInnerSet[n]] = d3 curvature = curvatureList[xTCInnerSet[n]] - distanceToInnerIdx = vector.magnitude([xTCOuter[i] - xTCInner[i] for i in range(3)]) + distanceToInnerIdx = magnitude([xTCOuter[i] - xTCInner[i] for i in range(3)]) factor = 1.0 - curvature * distanceToInnerIdx d2 = [factor * c for c in d2List[xTCInnerSet[n]]] d2TCRaw.append(d2) @@ -1621,7 +1620,7 @@ def createFlatCoordinatesTeniaColi(xiList, relaxedLengthList, v2 = xFlatListTC[nodeIdx] d1 = d2 = [v1[i] - v2[i] for i in range(3)] arclength = interp.computeCubicHermiteArcLength(v1, d1, v2, d2, True) - d2Flat = vector.setMagnitude(d1, arclength) + d2Flat = set_magnitude(d1, arclength) d2FlatListTC.append(d2Flat) d2FlatListTC = d2FlatListTC + d2FlatListTC[-((elementsCountAroundTC - 1) * tcCount + 1):] @@ -1682,14 +1681,14 @@ def createColonCoordinatesTeniaColi(xiList, relativeThicknessList, lengthToDiame elementsCountAroundTC * 0.5) TCEndIdx = TCMidIdx + int(elementsCountAroundTC * 0.5) v1 = xColon[TCStartIdx] - norm = vector.setMagnitude( - vector.crossproduct3(vector.normalise(d1Colon[TCMidIdx]), vector.normalise(d2Colon[TCMidIdx])), + norm = set_magnitude( + cross(normalize(d1Colon[TCMidIdx]), normalize(d2Colon[TCMidIdx])), teniaColiThicknessToDiameterRatio) v2 = [xColon[TCMidIdx][c] + norm[c] for c in range(3)] v3 = xColon[TCEndIdx] nx = [v1, v2, v3] nd1 = [d1Colon[TCStartIdx], - vector.setMagnitude(d1Colon[TCMidIdx], 2.0 * vector.magnitude(d1Colon[TCMidIdx])), + set_magnitude(d1Colon[TCMidIdx], 2.0 * magnitude(d1Colon[TCMidIdx])), d1Colon[TCEndIdx]] sx, sd1 = interp.sampleCubicHermiteCurves(nx, nd1, elementsCountAroundTC)[0:2] xTCRaw += sx[1:-1] diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_esophagus1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_esophagus1.py index ac378894..7d9d4070 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_esophagus1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_esophagus1.py @@ -21,7 +21,6 @@ from scaffoldmaker.utils import geometry from scaffoldmaker.utils import interpolation as interp from scaffoldmaker.utils import tubemesh -from scaffoldmaker.utils import vector from scaffoldmaker.utils.zinc_utils import exnode_string_from_nodeset_field_parameters,\ get_nodeset_path_field_parameters @@ -366,8 +365,8 @@ def createEsophagusMesh3d(region, options, networkLayout, nextNodeIdentifier, ne for n2 in range(elementsCountAlong + 1): # Create inner points cx = [0.0, 0.0, elementAlongLength * n2] - axis1 = [vector.magnitude(majorRadiusElementList[n2]), 0.0, 0.0] - axis2 = [0.0, vector.magnitude(minorRadiusElementList[n2]), 0.0] + axis1 = [magnitude(majorRadiusElementList[n2]), 0.0, 0.0] + axis2 = [0.0, magnitude(minorRadiusElementList[n2]), 0.0] xInner, d1Inner = geometry.createEllipsePoints(cx, 2 * math.pi, axis1, axis2, elementsCountAround, startRadians=0.0) xToSample += xInner @@ -413,7 +412,7 @@ def createEsophagusMesh3d(region, options, networkLayout, nextNodeIdentifier, ne v2 = xAround[(n1 + 1) % elementsCountAround] d1 = d2 = [v2[c] - v1[c] for c in range(3)] arcLengthAround = interp.computeCubicHermiteArcLength(v1, d1, v2, d2, True) - dx_ds1 = [c * arcLengthAround for c in vector.normalise(d1)] + dx_ds1 = [c * arcLengthAround for c in normalize(d1)] d1Around.append(dx_ds1) d1Smoothed = interp.smoothCubicHermiteDerivativesLoop(xAround, d1Around) @@ -450,7 +449,7 @@ def createEsophagusMesh3d(region, options, networkLayout, nextNodeIdentifier, ne for n2 in range(elementsCountAlong + 1): firstNodeAlong = xToWarp[n2 * elementsCountAround] midptSegmentAxis = [0.0, 0.0, elementAlongLength * n2] - radius = vector.magnitude(firstNodeAlong[c] - midptSegmentAxis[c] for c in range(3)) + radius = magnitude(firstNodeAlong[c] - midptSegmentAxis[c] for c in range(3)) innerRadiusAlong.append(radius) xWarpedList, d1WarpedList, d2WarpedList, d3WarpedUnitList = \ diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_heartarterialroot1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_heartarterialroot1.py index 84a6e357..210d33a9 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_heartarterialroot1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_heartarterialroot1.py @@ -16,7 +16,6 @@ from scaffoldmaker.annotation.heart_terms import get_heart_term from scaffoldmaker.meshtypes.scaffold_base import Scaffold_base from scaffoldmaker.utils import interpolation as interp -from scaffoldmaker.utils import vector from scaffoldmaker.utils.eft_utils import remapEftLocalNodes, remapEftNodeValueLabel, setEftScaleFactorIds from scaffoldmaker.utils.eftfactory_bicubichermitelinear import eftfactory_bicubichermitelinear from scaffoldmaker.utils.geometry import getApproximateEllipsePerimeter, createCirclePoints @@ -166,7 +165,7 @@ def generateBaseMesh(cls, region, options, baseCentre=[ 0.0, 0.0, 0.0 ], axisSid elementsCountAround = 6 radiansPerElementAround = 2.0*math.pi/elementsCountAround - axisSide2 = vector.crossproduct3(axisUp, axisSide1) + axisSide2 = cross(axisUp, axisSide1) outerRadius = innerRadius + wallThickness cuspOuterLength2 = 0.5*getApproximateEllipsePerimeter(innerRadius, cuspHeight) cuspOuterWallArcLength = cuspOuterLength2*innerRadius/(innerRadius + cuspHeight) @@ -275,13 +274,13 @@ def generateBaseMesh(cls, region, options, baseCentre=[ 0.0, 0.0, 0.0 ], axisSid n2m = n1m*2 + 1 n2p = n1p*2 - 1 ax, ad1 = interp.sampleCubicHermiteCurves([ upperx[0][n2m], sinusx[n1] ], [ [ -scaled1*d for d in upperd2[0][n2m] ], [ scaled1*d for d in sinusd1[n1] ] ], - elementsCountOut=2, addLengthStart=0.5*vector.magnitude(upperd2[0][n2m]), lengthFractionStart=0.5, - addLengthEnd=0.5*vector.magnitude(sinusd1[n1]), lengthFractionEnd=0.5, arcLengthDerivatives=False)[0:2] + elementsCountOut=2, addLengthStart=0.5*magnitude(upperd2[0][n2m]), lengthFractionStart=0.5, + addLengthEnd=0.5*magnitude(sinusd1[n1]), lengthFractionEnd=0.5, arcLengthDerivatives=False)[0:2] arcx .append(ax [1]) arcd1.append(ad1[1]) ax, ad1 = interp.sampleCubicHermiteCurves([ sinusx[n1], upperx[0][n2p], ], [ [ scaled1*d for d in sinusd1[n1] ], [ scaled1*d for d in upperd2[0][n2p] ] ], - elementsCountOut=2, addLengthStart=0.5*vector.magnitude(sinusd1[n1]), lengthFractionStart=0.5, - addLengthEnd=0.5*vector.magnitude(upperd2[0][n2p]), lengthFractionEnd=0.5, arcLengthDerivatives=False)[0:2] + elementsCountOut=2, addLengthStart=0.5*magnitude(sinusd1[n1]), lengthFractionStart=0.5, + addLengthEnd=0.5*magnitude(upperd2[0][n2p]), lengthFractionEnd=0.5, arcLengthDerivatives=False)[0:2] arcx .append(ax [1]) arcd1.append(ad1[1]) if nMidCusp == 0: @@ -307,13 +306,13 @@ def generateBaseMesh(cls, region, options, baseCentre=[ 0.0, 0.0, 0.0 ], axisSid cosRadiansAround = math.cos(radiansAround) sinRadiansAround = math.sin(radiansAround) noduled1[0].append([ noduleOuterRadialArcLength*(cosRadiansAround*axisSide1[c] + sinRadiansAround*axisSide2[c]) for c in range(3) ]) - noduled1[1].append(vector.setMagnitude(noduled1[0][i], noduleInnerRadialArcLength)) + noduled1[1].append(set_magnitude(noduled1[0][i], noduleInnerRadialArcLength)) n1 = i*2 + 1 + nMidCusp radiansAround = n1*radiansPerElementAround cosRadiansAround = math.cos(radiansAround) sinRadiansAround = math.sin(radiansAround) noduled2[0].append([ noduleOuterRadialArcLength*(cosRadiansAround*axisSide1[c] + sinRadiansAround*axisSide2[c]) for c in range(3) ]) - noduled2[1].append(vector.setMagnitude(noduled2[0][i], noduleInnerRadialArcLength)) + noduled2[1].append(set_magnitude(noduled2[0][i], noduleInnerRadialArcLength)) noduled3[0].append([ noduleOuterAxialArcLength*axisUp[c] for c in range(3) ]) noduled3[1].append([ noduleInnerAxialArcLength*axisUp[c] for c in range(3) ]) @@ -551,8 +550,8 @@ def getSemilunarValveSinusPoints(centre, axisSide1, axisSide2, radius, sinusRadi pd1.append([ radiansPerElementAround*(-rsinRadiansAround*axisSide1[c] + rcosRadiansAround*axisSide2[c]) for c in range(3) ]) # smooth to get derivative in sinus sd1 = interp.smoothCubicHermiteDerivativesLine(px[1 - nMidCusp:3 - nMidCusp], pd1[1 - nMidCusp:3 - nMidCusp], fixStartDerivative=True, fixEndDirection=True) - magSinus = vector.magnitude(sd1[1]) + magSinus = magnitude(sd1[1]) for n in range(nMidCusp, elementsCountAround, 2): - pd1[n] = vector.setMagnitude(pd1[n], magSinus) + pd1[n] = set_magnitude(pd1[n], magSinus) return px, pd1 diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_heartarterialvalve1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_heartarterialvalve1.py index 1e313ad1..711e3802 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_heartarterialvalve1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_heartarterialvalve1.py @@ -15,7 +15,6 @@ from cmlibs.zinc.field import Field from cmlibs.zinc.node import Node from scaffoldmaker.meshtypes.scaffold_base import Scaffold_base -from scaffoldmaker.utils import vector from scaffoldmaker.utils.eft_utils import setEftScaleFactorIds from scaffoldmaker.utils.eftfactory_tricubichermite import eftfactory_tricubichermite from scaffoldmaker.utils.geometry import createCirclePoints @@ -145,7 +144,7 @@ def getPoints(cls, options): #centre = [ 0.0, 0.0, 0.0 ] #axis1 = [ 1.0, 0.0, 0.0 ] #axis2 = [ 0.0, 1.0, 0.0 ] - #axis3 = vector.crossproduct3(axis1, axis2) + #axis3 = cross(axis1, axis2) axis1, axis2, axis3 = eulerToRotationMatrix3([ rotationAzimuthRadians, rotationElevationRadians, rotationRollRadians ]) x = [ [ None, None ], [ None, None ] ] @@ -165,7 +164,7 @@ def getPoints(cls, options): leafd1mag = innerInletRadius*radiansPerElementAround # was 0.5* leafd2r, leafd2z = interpolateLagrangeHermiteDerivative([ innerRadialDisplacement, 0.0 ], [ 0.0, outletLength ], [ 0.0, outletLength ], 0.0) sinusd1mag = innerInletSinusRadius*radiansPerElementAround # initial value only - sinusd1mag = vector.magnitude(smoothCubicHermiteDerivativesLine( + sinusd1mag = magnitude(smoothCubicHermiteDerivativesLine( [ [ innerInletRadius, 0.0, inletz ], [ innerInletSinusRadius*math.cos(pi_3), innerInletSinusRadius*math.sin(pi_3), sinusz ] ], [ [ 0.0, leafd1mag, 0.0 ], [ -sinusd1mag*math.sin(pi_3), sinusd1mag*math.cos(pi_3), 0.0 ] ], fixStartDerivative = True, fixEndDirection = True)[1]) @@ -215,7 +214,7 @@ def getPoints(cls, options): extSinusRadius = outerRadius + outerSinusRadialDisplacement leafd1mag = extRadius*radiansPerElementAround sinusd1mag = extSinusRadius*radiansPerElementAround # initial value only - sinusd1mag = vector.magnitude(smoothCubicHermiteDerivativesLine( + sinusd1mag = magnitude(smoothCubicHermiteDerivativesLine( [ [ extRadius, 0.0, 0.0 ], [ extSinusRadius*math.cos(pi_3), extSinusRadius*math.sin(pi_3), 0.0 ] ], [ [ 0.0, leafd1mag, 0.0 ], [ -sinusd1mag*math.sin(pi_3), sinusd1mag*math.cos(pi_3), 0.0 ] ], fixStartDerivative = True, fixEndDirection = True)[1]) diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_heartatria1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_heartatria1.py index b3089a68..a5504b9f 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_heartatria1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_heartatria1.py @@ -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 vector from scaffoldmaker.utils.annulusmesh import createAnnulusMesh3d from scaffoldmaker.utils.derivativemoothing import DerivativeScalingMode, DerivativeSmoothing from scaffoldmaker.utils.eft_utils import createEftElementSurfaceLayer, remapEftLocalNodes, remapEftNodeValueLabel, \ @@ -983,8 +982,8 @@ def generateBaseMesh(cls, region, options): ragProportionLengths3 = [(aVenousMidpointOver - lagcsProportion) / count3] * count3 ragProportionLengths = ragProportionLengths1 + ragProportionLengths2 + ragProportionLengths3 + [ lagcsProportion ] # get d1 magnitudes over crest at posterior, middle and anterior/aorta - d1a = vector.magnitude(labd1[1][0]) - d1p = vector.magnitude(labd1[1][elementsCountAroundLeftAtriumFreeWall]) + d1a = magnitude(labd1[1][0]) + d1p = magnitude(labd1[1][elementsCountAroundLeftAtriumFreeWall]) d1m = 0.25*d1p # GRC fudge factor - not necessarily used to get trackSurface! ragProportion = 0.0 ragProportions = [ 0.0 ] @@ -1003,7 +1002,7 @@ def generateBaseMesh(cls, region, options): d1f = d1p xid1 = 2.0*(ragProportion - 0.5) f1, _, f3, _ = interp.getCubicHermiteBasis(xid1) - agd2.append(vector.setMagnitude(d1t, agLength*0.5*(ragProportionLengths[e] + ragProportionLengths[e + 1]))) + agd2.append(set_magnitude(d1t, agLength*0.5*(ragProportionLengths[e] + ragProportionLengths[e + 1]))) agd3.append(normalize(cross(d1t, d2t))) if e in (0, elementsCountOverAtria - 2): # old calculation used for vestibule top @@ -1020,7 +1019,7 @@ def generateBaseMesh(cls, region, options): # Get heights of elements on aorta up interatrial groove, for building venous limit curves aoHeight1 = ragProportionLengths[0]*agLength # fix vestibule top d2 magnitude - agd2[-1] = vector.setMagnitude(agd2[-1], aVestibuleOuterHeight) + agd2[-1] = set_magnitude(agd2[-1], aVestibuleOuterHeight) # smooth d2 up to vestibule top agd2 = interp.smoothCubicHermiteDerivativesLine(agx, agd2, fixAllDirections = True, fixEndDerivative = True) # , magnitudeScalingMode = interp.DerivativeScalingMode.HARMONIC_MEAN) # GRC , magnitudeScalingMode = interp.DerivativeScalingMode.HARMONIC_MEAN) @@ -1030,7 +1029,7 @@ def generateBaseMesh(cls, region, options): # add posterior crux point agx .append(labx [1][elementsCountAroundLeftAtriumFreeWall]) agd1.append(labd1[1][elementsCountAroundLeftAtriumFreeWall]) - agd2.append(vector.setMagnitude(labd2[1][elementsCountAroundLeftAtriumFreeWall], aVestibuleOuterHeight)) + agd2.append(set_magnitude(labd2[1][elementsCountAroundLeftAtriumFreeWall], aVestibuleOuterHeight)) agd3.append(labd3[1][elementsCountAroundLeftAtriumFreeWall]) # copy derivatives to labd2[1], rabd2[1] rabd2[1][elementsCountAroundRightAtriumFreeWall] = labd2[1][0] = agd2[0] @@ -1075,7 +1074,7 @@ def generateBaseMesh(cls, region, options): cosRadiansAround = math.cos(radiansAround) sinRadiansAround = math.sin(radiansAround) fox .append([ 0.0, foMidpointY + foMagY*cosRadiansAround, foMidpointZ + foMagZ*sinRadiansAround ]) - fod1.append(vector.setMagnitude([ 0.0, -foMagY*sinRadiansAround, foMagZ*cosRadiansAround ], estElementSizeAroundFossa)) + fod1.append(set_magnitude([ 0.0, -foMagY*sinRadiansAround, foMagZ*cosRadiansAround ], estElementSizeAroundFossa)) fod2.append([ 0.0, -fossaOuterDerivativeRatio*foMagY*cosRadiansAround, -fossaOuterDerivativeRatio*foMagZ*sinRadiansAround ]) fod1 = interp.smoothCubicHermiteDerivativesLoop(fox, fod1, fixAllDirections = True) # GRC was , magnitudeScalingMode = interp.DerivativeScalingMode.HARMONIC_MEAN) fossad3const = [ foThickness, 0.0, 0.0 ] @@ -1108,7 +1107,7 @@ def generateBaseMesh(cls, region, options): nf = ns # make d2 in normal d3, d1 d2 = interp.interpolateLagrangeHermiteDerivative([ asx[ns][0] - halffoThickness, asx[ns][1], asx[ns][2] ], fox[0][nf], fod2[0][nf], 0.0) - d2 = vector.setMagnitude([ 0.0, -asd1[ns][2], asd1[ns][1] ], vector.magnitude(d2)) + d2 = set_magnitude([ 0.0, -asd1[ns][2], asd1[ns][1] ], magnitude(d2)) d2 = interp.smoothCubicHermiteDerivativesLine([ [ asx[ns][0] - halffoThickness, asx[ns][1], asx[ns][2] ], fox[0][nf] ], [ d2, fod2[0][nf] ], fixStartDirection = True, fixEndDerivative = True)[0] asd2.append(d2) @@ -1127,7 +1126,7 @@ def generateBaseMesh(cls, region, options): cosRadiansAround = math.cos(radiansAround) sinRadiansAround = math.sin(radiansAround) x = [ 0.0, foMidpointY + archMagY*cosRadiansAround, foMidpointZ + archMagZ*sinRadiansAround ] - d2 = vector.setMagnitude([ 0.0, -archMagY*sinRadiansAround, archMagZ*cosRadiansAround ], estArchElementLength) + d2 = set_magnitude([ 0.0, -archMagY*sinRadiansAround, archMagZ*cosRadiansAround ], estArchElementLength) elif agx[ng][1] > foMidpointY: xi = (aSeptumAnteriorY - foMidpointY)/(agx[ng][1] - foMidpointY) x = [ 0.0, aSeptumAnteriorY, (1.0 - xi)*foMidpointZ + xi*agx[ng][2] ] @@ -1138,7 +1137,7 @@ def generateBaseMesh(cls, region, options): d2 = [ 0.0, 0.0, -estArchElementLength ] # make d1 normal to d2, d3 d1 = interp.interpolateHermiteLagrangeDerivative([ 0.0, fox[0][nf][1], fox[0][nf][2] ], [ -d for d in fod2[0][nf] ], x, 1.0) - d1 = vector.setMagnitude([ 0.0, d2[2], -d2[1] ], vector.magnitude(d1)) + d1 = set_magnitude([ 0.0, d2[2], -d2[1] ], magnitude(d1)) d1 = interp.smoothCubicHermiteDerivativesLine([ [ 0.0, fox[0][nf][1], fox[0][nf][2] ], x ], [ [ -d for d in fod2[0][nf] ], d1 ], fixStartDerivative = True, fixEndDirection = True)[1] asx .append(x ) @@ -1220,22 +1219,22 @@ def generateBaseMesh(cls, region, options): startProportion2 = (elementsCountAroundLeftAtriumFreeWall - n1)/(elementsCountAroundLeftAtriumFreeWall - lan1Mid) startPosition = laTrackSurface.findNearestPosition(labx[1][n1], laTrackSurface.createPositionProportion(startProportion1, startProportion2)) onAorta = n1 == 1 - direction = [ 0.0, 0.0, 1.0 ] if onAorta else vector.normalise(labd2[1][n1]) + direction = [ 0.0, 0.0, 1.0 ] if onAorta else normalize(labd2[1][n1]) trackDistance1 = aoHeight1 if onAorta else aVestibuleOuterHeight position = laTrackSurface.trackVector(startPosition, direction, trackDistance1) lavtProportions.append(laTrackSurface.getProportion(position)) x, d1, d2 = laTrackSurface.evaluateCoordinates(position, derivatives = True) ax1, ax2, ax3 = calculate_surface_axes(d1, d2, direction) lavtx .append(x) - lavtd1.append(vector.setMagnitude(ax2, -vector.magnitude(labd1[1][n1]))) - lavtd2.append(vector.setMagnitude(ax1, trackDistance1)) + lavtd1.append(set_magnitude(ax2, -magnitude(labd1[1][n1]))) + lavtd2.append(set_magnitude(ax1, trackDistance1)) # if not commonLeftRightPvOstium and \ # (n1 == (elementsCountAroundLeftAtriumFreeWall - elementsCountAroundLeftAtriumRPV - 1)): # # derivative needs a tweak for LPV annulus - # lavtd2[-1] = sub(lavtd2[-1], vector.setMagnitude(lavtd1[-1], trackDistance1)) - lavtd3.append(vector.setMagnitude(ax3, laVenousFreeWallThickness)) + # lavtd2[-1] = sub(lavtd2[-1], set_magnitude(lavtd1[-1], trackDistance1)) + lavtd3.append(set_magnitude(ax3, laVenousFreeWallThickness)) # fix d2 on outer base - labd2[1][n1] = vector.setMagnitude(labd2[1][n1], trackDistance1) + labd2[1][n1] = set_magnitude(labd2[1][n1], trackDistance1) # add end points and smooth d1 lavtx .append(agx [-2]) lavtd1.append(agd1[-2]) @@ -1290,18 +1289,18 @@ def generateBaseMesh(cls, region, options): startPosition = raTrackSurface.findNearestPosition(rabx[1][n1], raTrackSurface.createPositionProportion(startProportion1, 0.5)) rabProportions.append(raTrackSurface.getProportion(startPosition)) onAorta = n1 == (elementsCountAroundRightAtriumFreeWall - 1) - direction = [ 0.0, 0.0, 1.0 ] if onAorta else vector.normalise(rabd2[1][n1]) + direction = [ 0.0, 0.0, 1.0 ] if onAorta else normalize(rabd2[1][n1]) trackDistance1 = aoHeight1 if onAorta else aVestibuleOuterHeight position = raTrackSurface.trackVector(startPosition, direction, trackDistance1) ravtProportions.append(raTrackSurface.getProportion(position)) x, d1, d2 = raTrackSurface.evaluateCoordinates(position, derivatives = True) ax1, ax2, ax3 = calculate_surface_axes(d1, d2, direction) ravtx .append(x) - ravtd1.append(vector.setMagnitude(ax2, -vector.magnitude(rabd1[1][n1]))) - ravtd2.append(vector.setMagnitude(ax1, trackDistance1)) - ravtd3.append(vector.setMagnitude(ax3, raVenousFreeWallThickness)) + ravtd1.append(set_magnitude(ax2, -magnitude(rabd1[1][n1]))) + ravtd2.append(set_magnitude(ax1, trackDistance1)) + ravtd3.append(set_magnitude(ax3, raVenousFreeWallThickness)) # fix d2 on outer base - rabd2[1][n1] = vector.setMagnitude(rabd2[1][n1], trackDistance1) + rabd2[1][n1] = set_magnitude(rabd2[1][n1], trackDistance1) # add end points and smooth d1 ravtx .append(agx [1]) ravtd1.append(agd1[1]) @@ -1357,13 +1356,13 @@ def generateBaseMesh(cls, region, options): laoaProportions = laoaProportions[:1] + _laoaProportions # d2 magnitudes are the same as d1 by default; want these to be even instead # transition between start and end side derivatives d1 - mag0 = vector.magnitude(lavtd1[1][1]) - magn = vector.magnitude(lavtd1[1][lan1Mid]) + mag0 = magnitude(lavtd1[1][1]) + magn = magnitude(lavtd1[1][lan1Mid]) # could blend accurately with arc lengths; try indexes: laoaCount = len(laoaProportions) for n1 in range(1, laoaCount - 1): xi = n1 / laoaCount - laoad2[n1] = vector.setMagnitude(laoad2[n1], (1.0 - xi)*mag0 + xi*magn) + laoad2[n1] = set_magnitude(laoad2[n1], (1.0 - xi)*mag0 + xi*magn) # get inner points laoax = [ [ None ], laoax ] laoad1 = [ [ None ], laoad1 ] @@ -1398,8 +1397,8 @@ def generateBaseMesh(cls, region, options): mpx, mpd1, mpd2, _, mpProportions = laTrackSurface.createHermiteCurvePoints(rpvoProportion1, rpvoProportion2, lpvoProportion1, lpvoProportion2, elementsCount = 2, derivativeStart = rpvod2[1][n1rpv], derivativeEnd = [ -d for d in lpvod2[1][n1lpv] ]) # scale mid derivative 2 to be mean of d1 in LPV, RPV - d2mag = 0.5*vector.magnitude(lpvod1[1][n1lpv]) + 0.5*vector.magnitude(rpvod1[1][n1rpv]) - mpd2[1] = vector.setMagnitude(mpd2[1], d2mag) + d2mag = 0.5*magnitude(lpvod1[1][n1lpv]) + 0.5*magnitude(rpvod1[1][n1rpv]) + mpd2[1] = set_magnitude(mpd2[1], d2mag) lamlx, lamld2, lamld1, lamld3, lamlProportions = laTrackSurface.createHermiteCurvePoints( laoaProportions[1][0], laoaProportions[1][1], mpProportions[1][0], mpProportions[1][1], @@ -1423,7 +1422,7 @@ def generateBaseMesh(cls, region, options): lamld2 = interp.smoothCubicHermiteDerivativesLine(lamlx, lamld2, fixAllDirections=True, fixStartDerivative=True, fixEndDerivative=True) # give d1 the same magnitude as the smoothed d2 for n in range(1, len(lamld1)): - lamld1[n] = vector.setMagnitude(lamld1[n], vector.magnitude(lamld2[n])) + lamld1[n] = set_magnitude(lamld1[n], magnitude(lamld2[n])) # get inner points lamlx = [ [ None ], lamlx ] lamld1 = [ [ None ], lamld1 ] @@ -1594,7 +1593,7 @@ def generateBaseMesh(cls, region, options): raVenousWidth = sum(interp.getCubicHermiteArcLength(ravtx[1][e], ravtd1[1][e], ravtx[1][e + 1], ravtd1[1][e + 1]) for e in range(elementsCountAroundRightAtriumPosteriorVenous)) d1mag = -0.2*raVenousWidth # GRC fudge factor ractx = [ [], ractx ] - ractd1 = [ [], [ vector.setMagnitude(d1, d1mag) for d1 in ractd1 ] ] + ractd1 = [ [], [ set_magnitude(d1, d1mag) for d1 in ractd1 ] ] ractd2 = [ [], ractd2 ] ractd3 = [ [], ractd3 ] for n in range(len(ractx[1])): @@ -1607,8 +1606,8 @@ def generateBaseMesh(cls, region, options): for n3 in range(2): # overwrite venous right x, d1 on vestibule top ravtx [n3][ran1Ctp] = ractx [n3][-1] - d1mag = min(vector.magnitude(ravtd1[n3][ran1Ctp]), 1.0*vector.magnitude(ractd1[n3][-1])) # GRC fudge factor - ravtd1[n3][ran1Ctp] = vector.setMagnitude(ravtd1[n3][ran1Ctp], d1mag) + d1mag = min(magnitude(ravtd1[n3][ran1Ctp]), 1.0*magnitude(ractd1[n3][-1])) # GRC fudge factor + ravtd1[n3][ran1Ctp] = set_magnitude(ravtd1[n3][ran1Ctp], d1mag) # substitute known start and end coordinates ractx [0][ 0] = asx [1][elementsCountAroundAtrialSeptum] ractd1[0][ 0] = asd1[1][elementsCountAroundAtrialSeptum] @@ -1639,11 +1638,11 @@ def generateBaseMesh(cls, region, options): ravmd2 = [ [], ravmd2 ] ravmd3 = [ [], ravmd3 ] # blend d2 between ends: - magc = vector.magnitude(ractd2[1][elementsCountOverCristaTerminalisAnterior]) - maga = vector.magnitude(agd2[agn1Mid]) + magc = magnitude(ractd2[1][elementsCountOverCristaTerminalisAnterior]) + maga = magnitude(agd2[agn1Mid]) for n in range(1, elementsCountOverSideRightAtriumVC): xi = n/elementsCountOverSideRightAtriumVC - ravmd2[1][n] = vector.setMagnitude(ravmd2[1][n], (1.0 - xi)*magc + xi*maga) + ravmd2[1][n] = set_magnitude(ravmd2[1][n], (1.0 - xi)*magc + xi*maga) for n in range(len(ravmx[1])): x, d1, d2, d3 = interp.projectHermiteCurvesThroughWall(ravmx[1], ravmd1[1], ravmd2[1], n, -raVenousFreeWallThickness) ravmx [0].append(x) @@ -1747,14 +1746,14 @@ def generateBaseMesh(cls, region, options): raaqx[1][1], d1, d2 = raTrackSurface.evaluateCoordinates(position, derivatives = True) # calculate derivative 1 there to fit nearest point on raaq raaqd1[1][1] = interp.interpolateLagrangeHermiteDerivative(raaqx[1][1], raaqx[1][2], raaqd1[1][2], 0.0) - raaqd1[1][1] = vector.setMagnitude(calculate_surface_axes(d1, d2, raaqd1[1][1])[0], vector.magnitude(raaqd1[1][1])) + raaqd1[1][1] = set_magnitude(calculate_surface_axes(d1, d2, raaqd1[1][1])[0], magnitude(raaqd1[1][1])) raaqd1[1][1] = interp.smoothCubicHermiteDerivativesLine(raaqx[1][1:3], raaqd1[1][1:3], fixAllDirections = True, fixEndDerivative = True)[0] # calculate derivative 2 there and on nearest point on raap # raapd2[1][1] magnitude needs to be set to fit distance between nodes: d2mag = math.sqrt(sum((raapx[1][1][c] - raaqx[1][1][c])*(raapx[1][1][c] - raaqx[1][1][c]) for c in range(3))) - raapd2[1][1] = vector.setMagnitude(raapd2[1][1], d2mag) # must reduce this otherwise smooth will converge wrongly + raapd2[1][1] = set_magnitude(raapd2[1][1], d2mag) # must reduce this otherwise smooth will converge wrongly raaqd2[1][1] = interp.interpolateLagrangeHermiteDerivative(raaqx[1][1], raapx[1][1], raapd2[1][1], 0.0) - raaqd2[1][1] = vector.setMagnitude(calculate_surface_axes(d1, d2, raaqd2[1][1])[0], vector.magnitude(raaqd2[1][1])) + raaqd2[1][1] = set_magnitude(calculate_surface_axes(d1, d2, raaqd2[1][1])[0], magnitude(raaqd2[1][1])) raaqd2[1][1], raapd2[1][1] = interp.smoothCubicHermiteDerivativesLine([ raaqx[1][1], raapx[1][1] ], [ raaqd2[1][1], raapd2[1][1] ], fixAllDirections = True) # get inner coordinates and derivatives raaqx[0][1], raaqd1[0][1], _, raaqd3[0][1] = interp.projectHermiteCurvesThroughWall(raaqx[1][1:3], raaqd1[1][1:3], raaqd2[1][1:3], 0, -raaWallThickness) @@ -2710,7 +2709,7 @@ def generateBaseMesh(cls, region, options): vcEndDerivative = vcDerivativeFactor*vcLength/elementsCountAlongVCInlet ocxPosition = raTrackSurface.createPositionProportion(proportion1, proportion2) ocx, d1, d2 = raTrackSurface.evaluateCoordinates(ocxPosition, derivatives = True) - ocd1, ocd2, ocd3 = calculate_surface_axes(d1, d2, vector.normalise(d1)) + ocd1, ocd2, ocd3 = calculate_surface_axes(d1, d2, normalize(d1)) vcx, vd1, vd2, vd3 = getCircleProjectionAxes(ocx, ocd1, ocd2, ocd3, vcLength, vcAngle1Radians, vcAngle2Radians) vcd1 = vd1 vcd2 = [ -d for d in vd2 ] @@ -2720,7 +2719,7 @@ def generateBaseMesh(cls, region, options): vcvd2 = [ None, None ] for n3 in range(2): radius = vcInnerRadius if (n3 == 0) else vcOuterRadius - px, pd1 = createCirclePoints(vcx, vector.setMagnitude(vcd1, radius), vector.setMagnitude(vcd2, radius), elementsCountAroundVC, startRadians) + px, pd1 = createCirclePoints(vcx, set_magnitude(vcd1, radius), set_magnitude(vcd2, radius), elementsCountAroundVC, startRadians) vcvx [n3] = px vcvd1[n3] = pd1 vcvd2[n3] = [ vcd3 ]*elementsCountAroundVC @@ -2903,8 +2902,8 @@ def generateBaseMesh(cls, region, options): laamx, d1, d2 = laTrackSurface.evaluateCoordinates(position, derivatives = True) # force d2 to be vertical, d3, d1 to be horizontal laamd2 = [ 0.0, 0.0, 1.0 ] - laamd3 = vector.normalise(vector.crossproduct3(d2, laamd2)) - laamd1 = vector.crossproduct3(laamd2, laamd3) + laamd3 = normalize(cross(d2, laamd2)) + laamd1 = cross(laamd2, laamd3) if False: node = nodes.createNode(nodeIdentifier, nodetemplate) cache.setNode(node) @@ -3083,8 +3082,8 @@ def generateBaseMesh(cls, region, options): raamx, d1, d2 = raTrackSurface.evaluateCoordinates(position, derivatives = True) # force d2 to be vertical, d3, d1 to be horizontal raamd2 = [ 0.0, 0.0, 1.0 ] - raamd3 = vector.normalise(vector.crossproduct3(raamd2, d2)) - raamd1 = vector.crossproduct3(raamd2, raamd3) + raamd3 = normalize(cross(raamd2, d2)) + raamd1 = cross(raamd2, raamd3) if False: node = nodes.createNode(nodeIdentifier, nodetemplate) cache.setNode(node) @@ -3266,8 +3265,8 @@ def generateBaseMesh(cls, region, options): result, epid2 = coordinates.getNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, 3) result, epid3 = coordinates.getNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS3, 1, 3) if result != RESULT_OK: - epid3 = vector.crossproduct3(epid1, epid2) - fatx = add(epix, vector.setMagnitude(epid3, epicardiumLayerMinimumThickness)) + epid3 = cross(epid1, epid2) + fatx = add(epix, set_magnitude(epid3, epicardiumLayerMinimumThickness)) epifx = None epifPosition = fpTrackSurface.findNearestPosition(epix, startPosition=None) @@ -3810,7 +3809,7 @@ def getAtriumBasePoints(elementsCountAroundAtrialSeptum, elementsCountAroundLeft x = [ laCentreX + cosRadiansAround*aMajorX + sinRadiansAround*aMinorX, laCentreY + cosRadiansAround*aMajorY + sinRadiansAround*aMinorY, 0.0 ] - d1 = vector.setMagnitude([ -sinRadiansAround*aMajorX + cosRadiansAround*aMinorX, + d1 = set_magnitude([ -sinRadiansAround*aMajorX + cosRadiansAround*aMinorX, -sinRadiansAround*aMajorY + cosRadiansAround*aMinorY, 0.0 ], ltFreeWallElementLength) if radiansAround < sideRadians: @@ -3821,8 +3820,8 @@ def getAtriumBasePoints(elementsCountAroundAtrialSeptum, elementsCountAroundLeft baseInclineRadians = (1.0 - xi)*aBaseSideInclineRadians + xi*aBaseBackInclineRadians else: baseInclineRadians = aBaseBackInclineRadians - side = vector.normalise([ d1[1], -d1[0], 0.0 ]) - d2 = vector.setMagnitude([ (up[c]*math.cos(baseInclineRadians) + side[c]*math.sin(baseInclineRadians)) for c in range(3) ], baseDerivative2Scale) + side = normalize([ d1[1], -d1[0], 0.0 ]) + d2 = set_magnitude([ (up[c]*math.cos(baseInclineRadians) + side[c]*math.sin(baseInclineRadians)) for c in range(3) ], baseDerivative2Scale) ltBaseOuterx .append(x) ltBaseOuterd1.append(d1) ltBaseOuterd2.append(d2) @@ -3830,7 +3829,7 @@ def getAtriumBasePoints(elementsCountAroundAtrialSeptum, elementsCountAroundLeft xi = 0.85 # GRC fudge factor baseSeptumPosteriorx = [ 0.0, (1.0 - xi)*ltBaseOuterx[0][1] + xi*ltBaseOuterx[-1][1], 0.0 ] nx = [ ltBaseOuterx[-1], baseSeptumPosteriorx ] - nd1 = interp.smoothCubicHermiteDerivativesLine(nx, [ ltBaseOuterd1[-1], [ vector.magnitude(ltBaseOuterd1[-1]), 0.0, 0.0 ] ], + nd1 = interp.smoothCubicHermiteDerivativesLine(nx, [ ltBaseOuterd1[-1], [ magnitude(ltBaseOuterd1[-1]), 0.0, 0.0 ] ], fixStartDerivative = True, fixEndDirection = True ) ltBaseOuterx .append(nx[1]) # GRC fudge factor: derivative must be lower to fit inlets: @@ -3998,7 +3997,7 @@ def getAtriumBasePoints(elementsCountAroundAtrialSeptum, elementsCountAroundLeft x = [ laCentreX + cosRadiansAround*aMajorX + sinRadiansAround*aMinorX, laCentreY + cosRadiansAround*aMajorY + sinRadiansAround*aMinorY, z ] - d1 = vector.setMagnitude([ -sinRadiansAround*aMajorX + cosRadiansAround*aMinorX, + d1 = set_magnitude([ -sinRadiansAround*aMajorX + cosRadiansAround*aMinorX, -sinRadiansAround*aMajorY + cosRadiansAround*aMinorY, 0.0 ], derivativeLength) if (n1 < 1) or (n1 >= elementsCountAroundAtriumFreeWall): @@ -4012,8 +4011,8 @@ def getAtriumBasePoints(elementsCountAroundAtrialSeptum, elementsCountAroundLeft baseInclineRadians = (1.0 - xi)*aBaseSideInclineRadians + xi*aBaseBackInclineRadians else: baseInclineRadians = aBaseBackInclineRadians - side = vector.normalise([ d1[1], -d1[0], 0.0 ]) - d2 = vector.setMagnitude([ (up[c]*math.cos(baseInclineRadians) + side[c]*math.sin(baseInclineRadians)) for c in range(3) ], baseDerivative2Scale) + side = normalize([ d1[1], -d1[0], 0.0 ]) + d2 = set_magnitude([ (up[c]*math.cos(baseInclineRadians) + side[c]*math.sin(baseInclineRadians)) for c in range(3) ], baseDerivative2Scale) if n3 == 0: aBaseInnerx.append(x) aBaseInnerd1.append(d1) @@ -4128,7 +4127,7 @@ def getAtriumTrackSurface(elementsCountAroundTrackSurface, elementsCountAcrossTr nd2 += ld2 # to be smoothed when all rows assembled. # add last point at central outer point of atrium, d1 all zero, d2 rotating around d1 = vd1[elementsCountAlongTrackSurface] - d2 = vector.setMagnitude(vd2[elementsCountAlongTrackSurface], vector.magnitude(d1)) + d2 = set_magnitude(vd2[elementsCountAlongTrackSurface], magnitude(d1)) for n1 in range(elementsCountAcrossTrackSurface + 1): nx.append(copy.deepcopy(vx[elementsCountAlongTrackSurface])) nd1.append([ 0.0, 0.0, 0.0 ]) diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_heartatria2.py b/src/scaffoldmaker/meshtypes/meshtype_3d_heartatria2.py index 6b1edb96..f2df8fb0 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_heartatria2.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_heartatria2.py @@ -16,7 +16,6 @@ from scaffoldmaker.annotation.heart_terms import get_heart_term from scaffoldmaker.meshtypes.scaffold_base import Scaffold_base from scaffoldmaker.utils import interpolation as interp -from scaffoldmaker.utils import vector from scaffoldmaker.utils.eft_utils import remapEftLocalNodes, remapEftNodeValueLabel, scaleEftNodeValueLabels, setEftScaleFactorIds from scaffoldmaker.utils.eftfactory_tricubichermite import eftfactory_tricubichermite from scaffoldmaker.utils.geometry import getApproximateEllipsePerimeter, getEllipseArcLength, getEllipseRadiansToX, updateEllipseAngleByArcLength @@ -362,7 +361,7 @@ def generateBaseMesh(cls, region, options): aSeptumCentreY = laCentreY + math.cos(laSeptumRadians)*aBaseInnerMajorMag*math.sin(-aMajorAxisRadians) \ + math.sin(laSeptumRadians)*aBaseInnerMinorMag*math.cos(-aMajorAxisRadians) tx1 = [ -0.5*aSeptumThickness, aSeptumCentreY, -aBaseSlopeHeight ] - ax1 = vector.normalise([ (sx2[c] - sx1[c]) for c in range(3) ]) + ax1 = normalize([ (sx2[c] - sx1[c]) for c in range(3) ]) ax2 = [ -ax1[1], ax1[0], ax1[2] ] ax3 = [ 0.0, 0.0, 1.0 ] @@ -544,13 +543,13 @@ def generateBaseMesh(cls, region, options): cosRadiansUp*aEquatorMinorMag*sinMajorAxisRadians, -cosRadiansUp*aEquatorMinorMag*cosMajorAxisRadians, -sinRadiansUp*aScaleZ ] - minorScale = elementSizeUpMinor/vector.magnitude(d2Minor) + minorScale = elementSizeUpMinor/magnitude(d2Minor) d2Minor = [ d*minorScale for d in d2Minor ] d2Major = [ -cosRadiansUpMajor*aEquatorMajorMag*cosMajorAxisRadians, -cosRadiansUpMajor*aEquatorMajorMag*sinMajorAxisRadians, -sinRadiansUpMajor*aMajorScaleZ ] - majorScale = derivativesUpMajor[n2]/vector.magnitude(d2Major) + majorScale = derivativesUpMajor[n2]/magnitude(d2Major) d2Major = [ d*majorScale for d in d2Major ] dx_ds2 = [ @@ -604,13 +603,13 @@ def generateBaseMesh(cls, region, options): result, d3 = coordinates.getNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, 3) d3 = [ -d for d in d3 ] arcLength = interp.computeCubicHermiteArcLength(x, dx_ds1, x3, d3, True) - scale1 = (2.0*arcLength - vector.magnitude(d3))/vector.magnitude(dx_ds1) + scale1 = (2.0*arcLength - magnitude(d3))/magnitude(dx_ds1) dx_ds1 = [ d*scale1 for d in dx_ds1 ] node = nodes.createNode(nodeIdentifier, nodetemplate) apexNodeId.append(nodeIdentifier) cache.setNode(node) - #print(n3, i, 'project apex', nid1, nid2, nid3, vector.magnitude(dx_ds1)) + #print(n3, i, 'project apex', nid1, nid2, nid3, magnitude(dx_ds1)) x[2] = aOuterHeight - aFreeWallThickness if (n3 == 0) else aOuterHeight result = coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, x) coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 1, dx_ds1) @@ -664,7 +663,7 @@ def generateBaseMesh(cls, region, options): d1 = [ (xc[c] - x1[c]) for c in range(3) ] d2 = [ (xc[c] - x2[c]) for c in range(3) ] dx_ds3 = [ d1[0] + d2[0], d1[1] + d2[1], d1[2] + d2[2] ] - scale = vector.magnitude(d1)/vector.magnitude(dx_ds3) + scale = magnitude(d1)/magnitude(dx_ds3) dx_ds3 = [ d*scale for d in dx_ds3 ] result = coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS3, 1, dx_ds3 ) @@ -687,9 +686,9 @@ def generateBaseMesh(cls, region, options): d2 = [ 1.0, 0.0, 0.0 ] #print('septum top centre ', x2) arcLength = interp.computeCubicHermiteArcLength(x1, d1, x2, d2, True) - scale1 = arcLength/vector.magnitude(d1) + scale1 = arcLength/magnitude(d1) d1 = [ d*scale1 for d in d1 ] - scale2 = arcLength/vector.magnitude(d2) + scale2 = arcLength/magnitude(d2) d2 = [ d*arcLength for d in d2 ] # GRC fudge factor: xi = 0.7 @@ -697,11 +696,11 @@ def generateBaseMesh(cls, region, options): lasmd1 = interp.interpolateCubicHermiteDerivative(x1, d1, x2, d2, xi) scale1 = (1.0 - xi)*2.0 lasmd1 = [ d*scale1 for d in lasmd1 ] - lasmd2 = vector.normalise(vector.crossproduct3([0.0, 0.0, 1.0], lasmd1)) - lasmd3 = vector.crossproduct3(lasmd1, lasmd2) - scale3 = aFreeWallThickness/vector.magnitude(lasmd3) + lasmd2 = normalize(cross([0.0, 0.0, 1.0], lasmd1)) + lasmd3 = cross(lasmd1, lasmd2) + scale3 = aFreeWallThickness/magnitude(lasmd3) lasmd3 = [ d*scale3 for d in lasmd3 ] - lasmCurvature1 = -interp.getCubicHermiteCurvature(x1, d1, x2, d2, vector.normalise(lasmd3), xi) + lasmCurvature1 = -interp.getCubicHermiteCurvature(x1, d1, x2, d2, normalize(lasmd3), xi) for n3 in range(2): for n2 in range(2): @@ -717,7 +716,7 @@ def generateBaseMesh(cls, region, options): result, laspd3 = coordinates.getNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS3, 1, 3) cache.setNode(nodes.findNodeByIdentifier(laspNodeId[1][n2])) result, laspd1_outer = coordinates.getNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 1, 3) - laspCurvature1 = (vector.magnitude(laspd1_outer)/vector.magnitude(laspd1) - 1.0)/aFreeWallThickness + laspCurvature1 = (magnitude(laspd1_outer)/magnitude(laspd1) - 1.0)/aFreeWallThickness cache.setNode(nodes.findNodeByIdentifier(lasaNodeId[0][n2])) result, lasax = coordinates.getNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, 3) result, lasad1 = coordinates.getNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 1, 3) @@ -725,7 +724,7 @@ def generateBaseMesh(cls, region, options): result, lasad3 = coordinates.getNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS3, 1, 3) cache.setNode(nodes.findNodeByIdentifier(lasaNodeId[1][n2])) result, lasad1_outer = coordinates.getNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 1, 3) - lasaCurvature1 = (vector.magnitude(lasad1_outer)/vector.magnitude(lasad1) - 1.0)/aFreeWallThickness + lasaCurvature1 = (magnitude(lasad1_outer)/magnitude(lasad1) - 1.0)/aFreeWallThickness # create points arcing over septum peak for n3 in range(2): @@ -733,10 +732,10 @@ def generateBaseMesh(cls, region, options): #print(str(n3) + '. ra septum nodes p', raspNodeId[n3][1],'a', rasaNodeId[n3][1]) arcLength = interp.computeCubicHermiteArcLength(laspx, laspd2, lasmx, lasmd2, True) x1 = laspx - scale1 = arcLength/vector.magnitude(laspd2) + scale1 = arcLength/magnitude(laspd2) d1 = [ d*scale1 for d in laspd2 ] x2 = lasmx - scale2 = arcLength/vector.magnitude(lasmd2) + scale2 = arcLength/magnitude(lasmd2) d2 = [ d*scale2 for d in lasmd2 ] derivativeScale = arcLength/(elementsCountUpAtria - 1.0) for n2 in range(2, elementsCountUpAtria + 1): @@ -745,16 +744,16 @@ def generateBaseMesh(cls, region, options): x = interp.interpolateCubicHermite(x1, d1, x2, d2, xi) dx_ds1 = [ (xi*lasmd1[c] + xr*laspd1[c]) for c in range(3) ] dx_ds2 = interp.interpolateCubicHermiteDerivative(x1, d1, x2, d2, xi) - scale2 = derivativeScale/vector.magnitude(dx_ds2) + scale2 = derivativeScale/magnitude(dx_ds2) dx_ds2 = [ d*scale2 for d in dx_ds2 ] - dx_ds3 = vector.crossproduct3(dx_ds1, dx_ds2) - scale3 = aFreeWallThickness/vector.magnitude(dx_ds3) + dx_ds3 = cross(dx_ds1, dx_ds2) + scale3 = aFreeWallThickness/magnitude(dx_ds3) dx_ds3 = [ d*scale3 for d in dx_ds3 ] if n3 == 1: curvature1 = xi*lasmCurvature1 + xr*laspCurvature1 curvatureScale1 = 1.0 + aFreeWallThickness*curvature1 dx_ds1 = [ d*curvatureScale1 for d in dx_ds1 ] - radialVector = vector.normalise(dx_ds3) + radialVector = normalize(dx_ds3) curvature2 = -interp.getCubicHermiteCurvature(x1, d1, x2, d2, radialVector, xi) curvatureScale2 = 1.0 + aFreeWallThickness*curvature2 dx_ds2 = [ d*curvatureScale2 for d in dx_ds2 ] @@ -785,10 +784,10 @@ def generateBaseMesh(cls, region, options): arcLength = interp.computeCubicHermiteArcLength(lasax, lasad2, lasmx, lasmd2, True) x1 = lasax - scale1 = arcLength/vector.magnitude(lasad2) + scale1 = arcLength/magnitude(lasad2) d1 = [ d*scale1 for d in lasad2 ] x2 = lasmx - scale2 = -arcLength/vector.magnitude(lasmd2) + scale2 = -arcLength/magnitude(lasmd2) d2 = [ d*scale2 for d in lasmd2 ] derivativeScale = arcLength/(elementsCountUpAtria - 1.0) for n2 in range(2, elementsCountUpAtria): @@ -797,16 +796,16 @@ def generateBaseMesh(cls, region, options): x = interp.interpolateCubicHermite(x1, d1, x2, d2, xi) dx_ds1 = [ (xi*-lasmd1[c] + xr*lasad1[c]) for c in range(3) ] dx_ds2 = interp.interpolateCubicHermiteDerivative(x1, d1, x2, d2, xi) - scale2 = derivativeScale/vector.magnitude(dx_ds2) + scale2 = derivativeScale/magnitude(dx_ds2) dx_ds2 = [ d*scale2 for d in dx_ds2 ] - dx_ds3 = vector.crossproduct3(dx_ds1, dx_ds2) - scale3 = aFreeWallThickness/vector.magnitude(dx_ds3) + dx_ds3 = cross(dx_ds1, dx_ds2) + scale3 = aFreeWallThickness/magnitude(dx_ds3) dx_ds3 = [ d*scale3 for d in dx_ds3 ] if n3 == 1: curvature1 = xi*lasmCurvature1 + xr*lasaCurvature1 curvatureScale1 = 1.0 + aFreeWallThickness*curvature1 dx_ds1 = [ d*curvatureScale1 for d in dx_ds1 ] - radialVector = vector.normalise(dx_ds3) + radialVector = normalize(dx_ds3) curvature2 = -interp.getCubicHermiteCurvature(x1, d1, x2, d2, radialVector, xi) curvatureScale2 = 1.0 + aFreeWallThickness*curvature2 dx_ds2 = [ d*curvatureScale2 for d in dx_ds2 ] @@ -879,7 +878,7 @@ def generateBaseMesh(cls, region, options): x[2] = aSeptumCentreZ + fossaMagZ*sinRadiansAround dx_ds1[1] = -fossaMagY*sinRadiansAround dx_ds1[2] = fossaMagZ*cosRadiansAround - scale1 = elementSizeAroundFossa/vector.magnitude(dx_ds1) + scale1 = elementSizeAroundFossa/magnitude(dx_ds1) dx_ds1[1] *= scale1 dx_ds1[2] *= scale1 dx_ds2[1] = -fossaMagY*cosRadiansAround @@ -1367,12 +1366,12 @@ def generateBaseMesh(cls, region, options): cache.setMeshLocation(svcElement1, [ 0.5, 1.0, 1.0 ]) resultx, sx = coordinates.evaluateReal(cache, 3) resultd, sd = coordinates.evaluateDerivative(diff1, cache, 3) - direction = vector.normalise([ (ix[c] - sx[c]) for c in range(3) ]) + direction = normalize([ (ix[c] - sx[c]) for c in range(3) ]) ivcInletScale = -ivcInnerRadius ivcInletAxis = [ d*ivcInletScale for d in direction ] - n = vector.crossproduct3(id, ivcInletAxis) - ivcInletSide = vector.crossproduct3(n, ivcInletAxis) + n = cross(id, ivcInletAxis) + ivcInletSide = cross(n, ivcInletAxis) ivcInletDistance = ivcInnerRadius ivcInletCentre = [ (ix[c] + direction[c]*ivcInletDistance) for c in range(3) ] tricubichermite.replaceTwoElementWithInlet6(ivcElement1, ivcElement2, elementIdentifier, nodetemplate, nodeIdentifier, \ @@ -1382,8 +1381,8 @@ def generateBaseMesh(cls, region, options): svcInletScale = -svcInnerRadius svcInletAxis = [ -d*svcInletScale for d in direction ] - n = vector.crossproduct3(sd, svcInletAxis) - svcInletSide = vector.crossproduct3(n, svcInletAxis) + n = cross(sd, svcInletAxis) + svcInletSide = cross(n, svcInletAxis) svcInletDistance = svcInnerRadius svcInletCentre = [ (sx[c] - direction[c]*svcInletDistance) for c in range(3) ] tricubichermite.replaceTwoElementWithInlet6(svcElement1, svcElement2, elementIdentifier, nodetemplate, nodeIdentifier, \ diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_heartventricles1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_heartventricles1.py index ad870597..5554014b 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_heartventricles1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_heartventricles1.py @@ -15,7 +15,6 @@ from scaffoldmaker.annotation.heart_terms import get_heart_term from scaffoldmaker.meshtypes.scaffold_base import Scaffold_base from scaffoldmaker.utils import interpolation as interp -from scaffoldmaker.utils import vector from scaffoldmaker.utils.eft_utils import remapEftLocalNodes, remapEftNodeValueLabel, scaleEftNodeValueLabels, setEftScaleFactorIds from scaffoldmaker.utils.eftfactory_tricubichermite import eftfactory_tricubichermite from scaffoldmaker.utils.geometry import getApproximateEllipsePerimeter, getEllipseArcLength, updateEllipseAngleByArcLength @@ -419,7 +418,7 @@ def generateBaseMesh(cls, region, options): nd2 = interp.smoothCubicHermiteDerivativesLine(nx, nd2, fixStartDirection=(n1 == 0), fixStartDerivative=(n1 > 0)) if n1 == 0: - mag = vector.magnitude((nd2[0])) + mag = magnitude((nd2[0])) lvApexInnerd1 = [mag, 0.0, 0.0] lvApexInnerd2 = [0.0, -mag, 0.0, 0.0] for n2 in range(elementsCountUpLV): @@ -514,8 +513,8 @@ def generateBaseMesh(cls, region, options): outerx = vOuterx [n2][n1] outerd1 = vOuterd1[n2][n1] outerd2 = vOuterd2[n2][n1] - unitNormal = vector.normalise(vector.crossproduct3(outerd1, outerd2)) - innerd3 = vOuterd3[n2][n1] = vector.setMagnitude(unitNormal, rvFreeWallThickness) + unitNormal = normalize(cross(outerd1, outerd2)) + innerd3 = vOuterd3[n2][n1] = set_magnitude(unitNormal, rvFreeWallThickness) innerx = [ (outerx[c] - innerd3[c]) for c in range(3) ] # calculate inner d1 from curvature around n1m = n1 - 1 @@ -573,8 +572,8 @@ def generateBaseMesh(cls, region, options): for n1 in [ -1, 0 ]: nd1[n1] = [ scale*d for d in nd1[n1] ] px, pd1 = interp.sampleCubicHermiteCurves(nx, nd1, elementsCountAroundVSeptum + 2, - addLengthStart = (0.5/scale)*vector.magnitude(nd1[ 0]), lengthFractionStart = 0.5, - addLengthEnd = (0.5/scale)*vector.magnitude(nd1[-1]), lengthFractionEnd = 0.5, + addLengthStart = (0.5/scale)*magnitude(nd1[ 0]), lengthFractionStart = 0.5, + addLengthEnd = (0.5/scale)*magnitude(nd1[-1]), lengthFractionEnd = 0.5, arcLengthDerivatives = False)[0:2] for n1 in range(elementsCountAroundVSeptum + 1): # d3 at ends of septum is toward adjacent interventricular sulcus, inside septum is across septum to lvInner @@ -629,8 +628,8 @@ def generateBaseMesh(cls, region, options): bd1 = rvInnerd1[n2][n1] bd2 = [ dFactor*d for d in rvInnerd2[n2][n1] ] px, pd2, pe, pxi = interp.sampleCubicHermiteCurves([ ax, bx ], [ ad2, bd2 ], elementsCountOut = 2, - addLengthStart = 0.5*vector.magnitude(ad2)/dFactor, lengthFractionStart = 0.5, - addLengthEnd = 0.5*vector.magnitude(bd2)/dFactor, lengthFractionEnd = 0.5, arcLengthDerivatives = False)[0:4] + addLengthStart = 0.5*magnitude(ad2)/dFactor, lengthFractionStart = 0.5, + addLengthEnd = 0.5*magnitude(bd2)/dFactor, lengthFractionEnd = 0.5, arcLengthDerivatives = False)[0:4] sx .append(px[1]) sd1.append(interp.interpolateSampleLinear([ ad1, bd1 ], pe[1:2], pxi[1:2])[0]) sd2.append(pd2[1]) @@ -641,8 +640,8 @@ def generateBaseMesh(cls, region, options): bd1 = [ -d for d in rvInnerd1[n2][elementsCountAroundRVFreeWall] ] bd2 = rvInnerd2[n2][elementsCountAroundRVFreeWall] px, pd1, pe, pxi = interp.sampleCubicHermiteCurves([ ax ] + sx + [ bx ], [ ad2 ] + sd1 + [ bd2 ], elementsCountAroundVSeptum + 2, - addLengthStart = 0.5*vector.magnitude(ad2)/dFactor, lengthFractionStart = 0.5, - addLengthEnd = 0.5*vector.magnitude(bd2)/dFactor, lengthFractionEnd = 0.5, arcLengthDerivatives = False)[0:4] + addLengthStart = 0.5*magnitude(ad2)/dFactor, lengthFractionStart = 0.5, + addLengthEnd = 0.5*magnitude(bd2)/dFactor, lengthFractionEnd = 0.5, arcLengthDerivatives = False)[0:4] pd2 = interp.interpolateSampleLinear([ ad1 ] + sd2 + [ bd1 ], pe, pxi) n2 = elementsCountUpLVApex - 1 # loop skips first and last in sample: @@ -1349,7 +1348,7 @@ def getLeftVentricleInnerPoints(lvRadius, midSeptumDisplacement, septumArcAround magnitude = elementSizeSep distance += elementSizeSep hx.append(x) - hd.append(vector.setMagnitude(d, magnitude)) + hd.append(set_magnitude(d, magnitude)) #distance -= (0.5*elementSizeSep if sepOdd else elementSizeSep) nx = [] @@ -1450,7 +1449,7 @@ def getVentriclesOuterPoints(lvRadius, widthExtension, sideExtension, addWidthRa magnitude = elementSizeRV distance += elementSizeRV hx.append(x) - hd.append(vector.setMagnitude(d, magnitude)) + hd.append(set_magnitude(d, magnitude)) #distance -= (0.5*elementSizeRV if rvOdd else elementSizeRV) nx = [] diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_heartventricles2.py b/src/scaffoldmaker/meshtypes/meshtype_3d_heartventricles2.py index 45217099..1ceb9667 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_heartventricles2.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_heartventricles2.py @@ -14,7 +14,6 @@ from scaffoldmaker.annotation.heart_terms import get_heart_term from scaffoldmaker.meshtypes.scaffold_base import Scaffold_base from scaffoldmaker.utils import interpolation as interp -from scaffoldmaker.utils import vector from scaffoldmaker.utils.eft_utils import remapEftNodeValueLabel, setEftScaleFactorIds from scaffoldmaker.utils.eftfactory_tricubichermite import eftfactory_tricubichermite from scaffoldmaker.utils.geometry import getApproximateEllipsePerimeter, getEllipseArcLength, updateEllipseAngleByArcLength @@ -277,8 +276,8 @@ def generateBaseMesh(cls, region, options): dzr_dRadiansUp = [ a*sinRadiansUp, b*cosRadiansUp ] else: dzr_dRadiansUp = [ a*sinRadiansUp, bSeptum*cosRadiansUp ] - ds_dRadiansUp = vector.magnitude(dzr_dRadiansUp) - dzr_ds = vector.normalise(dzr_dRadiansUp) + ds_dRadiansUp = magnitude(dzr_dRadiansUp) + dzr_ds = normalize(dzr_dRadiansUp) dz = dxi*ds_dxi*dzr_ds[0] dr = dxi*ds_dxi*dzr_ds[1] if (n3 == 1) and (n2 < elementsCountUpLVApex): @@ -295,7 +294,7 @@ def generateBaseMesh(cls, region, options): tx, td1 = getSeptumPoints(rvArcAroundRadians, radiusSeptum + dr, radialDisplacement + dRadialDisplacement, elementsCountAroundLVFreeWall, elementsCountAroundVSeptum, z + dz, n3) # true values for LV freewall dzr_dRadiansUp = [ a*sinRadiansUp, b*cosRadiansUp ] - dzr_ds = vector.normalise(dzr_dRadiansUp) + dzr_ds = normalize(dzr_dRadiansUp) for n1 in range(elementsCountAroundLV): @@ -380,7 +379,7 @@ def generateBaseMesh(cls, region, options): elementSizeCrossSeptum = 0.5*(lvOuterRadius*lvArcAroundRadians/elementsCountAroundLVFreeWall + vSeptumThickness*sinRadiansUp) dEndMag = elementSizeCrossSeptum ad2 = [ b*cosRadiansUp*cosRadiansAround, b*cosRadiansUp*sinRadiansAround, a*sinRadiansUp ] - scale = elementSizeUpRv/vector.magnitude(ad2) + scale = elementSizeUpRv/magnitude(ad2) if n2 == 0: scale *= elementSizeUpRvTransition1/elementSizeUpRv ad2 = [ d*scale for d in ad2 ] @@ -401,8 +400,8 @@ def generateBaseMesh(cls, region, options): dArcLengthUp = elementSizeUpRv ds_dxi = dArcLengthUp dzr_dRadiansUp = [ a*sinRadiansUp, b*cosRadiansUp ] - ds_dRadiansUp = vector.magnitude(dzr_dRadiansUp) - dzr_ds = vector.normalise(dzr_dRadiansUp) + ds_dRadiansUp = magnitude(dzr_dRadiansUp) + dzr_ds = normalize(dzr_dRadiansUp) dz = dxi*ds_dxi*dzr_ds[0] dr = dxi*ds_dxi*dzr_ds[1] radiansUpPlus = radiansUp + math.sqrt(dz*dz + dr*dr)/ds_dRadiansUp @@ -422,8 +421,8 @@ def generateBaseMesh(cls, region, options): rd3 = [] for n in range(len(rx)): - normal = vector.crossproduct3(rd1[n], rd2[n]) - mag = rvFreeWallThickness/vector.magnitude(normal) + normal = cross(rd1[n], rd2[n]) + mag = rvFreeWallThickness/magnitude(normal) rd3.append([ c*mag for c in normal ]) rxOuter.append(rx) @@ -453,7 +452,7 @@ def generateBaseMesh(cls, region, options): for n1 in range(1, elementsCountAroundRV - 2): ix.append([ (rx[n1][c] - rd3[n1][c]) for c in range(3) ] ) - unitRadial = vector.normalise(rd3[n1]) + unitRadial = normalize(rd3[n1]) # calculate inner d1 from curvature around curvature = 0.0 @@ -924,7 +923,7 @@ def getSeptumPoints(septumArcRadians, lvRadius, radialDisplacement, elementsCoun pos1 = interp.interpolateCubicHermite(v1, d1, v2, d2, xi) x.append([ v for v in pos1 ]) deriv1 = interp.interpolateCubicHermiteDerivative(v1, d1, v2, d2, xi) - scale = elementLengthMid/vector.magnitude(deriv1) + scale = elementLengthMid/magnitude(deriv1) dx_ds1.append([ d*scale for d in deriv1 ]) length += elementLengthMid return x, dx_ds1 diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_heartventricles3.py b/src/scaffoldmaker/meshtypes/meshtype_3d_heartventricles3.py index 32d2b0a6..1b881b70 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_heartventricles3.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_heartventricles3.py @@ -16,7 +16,6 @@ from scaffoldmaker.annotation.annotationgroup import AnnotationGroup, findOrCreateAnnotationGroupForTerm, getAnnotationGroupForTerm from scaffoldmaker.annotation.heart_terms import get_heart_term from scaffoldmaker.meshtypes.scaffold_base import Scaffold_base -from scaffoldmaker.utils import vector from scaffoldmaker.utils.eft_utils import remapEftNodeValueLabel, scaleEftNodeValueLabels, setEftScaleFactorIds from scaffoldmaker.utils.eftfactory_tricubichermite import eftfactory_tricubichermite from scaffoldmaker.utils.geometry import createEllipsoidPoints, getApproximateEllipsePerimeter, getEllipseArcLength, getEllipseRadiansToX @@ -254,7 +253,7 @@ def getRVEdgePoints(lvTrackSurface, position, cuspDirection, cuspAngleRadians, r :return: x, d1, d2, d3 """ x, sd1, sd2 = lvTrackSurface.evaluateCoordinates(position, derivatives = True) - td1, td2, td3 = calculate_surface_axes(sd1, sd2, vector.normalise(cuspDirection)) + td1, td2, td3 = calculate_surface_axes(sd1, sd2, normalize(cuspDirection)) cosCuspAngle = math.cos(cuspAngleRadians) sinCuspAngle = math.sin(cuspAngleRadians) pd1 = [ (cosCuspAngle*td1[c] - sinCuspAngle*td3[c]) for c in range(3) ] @@ -279,15 +278,15 @@ def getTransitionRVtoLV(lvTrackSurface, elementsCount, cuspProportions, cuspDeri tx, td1 = sampleCubicHermiteCurves(tx, td1, elementsCount, addLengthStart=(ivSulcusTransitionLength - uniformLength))[0:2] tx [0] = rvx td1[0] = rvd1 - tx, td1 = sampleCubicHermiteCurves(tx, td1, elementsCount, addLengthStart=0.5*vector.magnitude(td1[0]), lengthFractionStart=0.5, arcLengthDerivatives=True)[0:2] + tx, td1 = sampleCubicHermiteCurves(tx, td1, elementsCount, addLengthStart=0.5*magnitude(td1[0]), lengthFractionStart=0.5, arcLengthDerivatives=True)[0:2] p1 = lvTrackSurface.createPositionProportion(*cuspProportions) for i in range(0, elementsCount + 1): p1 = lvTrackSurface.findNearestPosition(tx[i], p1) x, sd1, sd2 = lvTrackSurface.evaluateCoordinates(p1, derivatives = True) - d1, d2, d3 = calculate_surface_axes(sd1, sd2, vector.normalise(td1[i])) - magd2 = vector.magnitude(d2) + d1, d2, d3 = calculate_surface_axes(sd1, sd2, normalize(td1[i])) + magd2 = magnitude(d2) if magd2 > 0.0: - scale = vector.magnitude(td1[i])/magd2 + scale = magnitude(td1[i])/magd2 td2[i] = [ d*scale for d in d2 ] else: td2[i] = [ 0.0, 0.0, 0.0 ] @@ -361,8 +360,8 @@ def generateBaseMesh(cls, region, options): ################# centre = [ 0.0, 0.0, 0.0 ] - poleAxis = vector.setMagnitude([ 0.0, 0.0, -1.0 ], lvOuterHeightAxisLength) - sideAxis = vector.setMagnitude([ -1.0, 0.0, 0.0 ], lvOuterRadius) + poleAxis = set_magnitude([ 0.0, 0.0, -1.0 ], lvOuterHeightAxisLength) + sideAxis = set_magnitude([ -1.0, 0.0, 0.0 ], lvOuterRadius) elementsCountAroundLVTrackSurface = 16 elementsCountUpLVTrackSurface = 8 @@ -379,11 +378,11 @@ def generateBaseMesh(cls, region, options): rvApexProportion2 = baseProportionUp*(1.0 - rvProportionUpLV) rvApexPosition = lvTrackSurface.createPositionProportion(rvApexProportion1, rvApexProportion2) rvApexCuspCoordinates, sd1, sd2 = lvTrackSurface.evaluateCoordinates(rvApexPosition, derivatives = True) - rvApexLengthDerivative = vector.setMagnitude(sd1, rvApexLengthFactor) - rvApexCuspDirection = vector.setMagnitude(sd2, -1.0) + rvApexLengthDerivative = set_magnitude(sd1, rvApexLengthFactor) + rvApexCuspDirection = set_magnitude(sd2, -1.0) rvax, rvad1, rvad2, rvad3 = cls.getRVEdgePoints(lvTrackSurface, rvApexPosition, rvApexCuspDirection, rvApexCuspAngleRadians, rvFreeWallThickness) - elementLengthUpLV = vector.magnitude(lvTrackSurface.createHermiteCurvePoints(0.0, 0.0, 0.0, baseProportionUp, elementsCountUpLVFreeWall)[2][-1]) + elementLengthUpLV = magnitude(lvTrackSurface.createHermiteCurvePoints(0.0, 0.0, 0.0, baseProportionUp, elementsCountUpLVFreeWall)[2][-1]) # Get RV inlet, outlet # sample curves around RV cusp @@ -395,7 +394,7 @@ def generateBaseMesh(cls, region, options): inletPosition = lvTrackSurface.createPositionProportion(rvInletPositionAroundLV, baseProportionUp) rvInletCuspCoordinates, sd1, sd2 = lvTrackSurface.evaluateCoordinates(inletPosition, derivatives = True) - td1, td2, rvInletSurfaceNormal = calculate_surface_axes(sd1, sd2, vector.normalise(sd1)) + td1, td2, rvInletSurfaceNormal = calculate_surface_axes(sd1, sd2, normalize(sd1)) cosAngle = math.cos(rvInletAngleRadians) sinAngle = math.sin(rvInletAngleRadians) rvBaseLengthFactor = (1.0 - rvApexLengthFactor) @@ -403,7 +402,7 @@ def generateBaseMesh(cls, region, options): outletPosition = lvTrackSurface.createPositionProportion(rvOutletPositionAroundLV, baseProportionUp) rvOutletCuspCoordinates, sd1, sd2 = lvTrackSurface.evaluateCoordinates(outletPosition, derivatives = True) - td1, td2, rvOutletSurfaceNormal = calculate_surface_axes(sd1, sd2, vector.normalise(sd1)) + td1, td2, rvOutletSurfaceNormal = calculate_surface_axes(sd1, sd2, normalize(sd1)) cosAngle = math.cos(rvOutletAngleRadians) sinAngle = math.sin(rvOutletAngleRadians) rvOutletDerivative = [ rvBaseLengthFactor*(cosAngle*td1[c] + sinAngle*td2[c]) for c in range(3) ] @@ -420,7 +419,7 @@ def generateBaseMesh(cls, region, options): rvix, rvid1, rvid2, rvid3, rviProportions = lvTrackSurface.createHermiteCurvePoints(rvInletPositionAroundLV, baseProportionUp, rvApexProportion1, rvApexProportion2, elementsCountAroundHalf, derivativeStart = rvInletDerivative, derivativeEnd = rvApexLengthDerivative, curveMode = TrackSurface.HermiteCurveMode.UNIFORM_SIZE) rvix, rvid1, rvid2, rvid3, rviProportions = lvTrackSurface.resampleHermiteCurvePointsSmooth(rvix, rvid1, rvid2, rvid3, rviProportions, - derivativeMagnitudeStart=elementLengthUpLV, derivativeMagnitudeEnd=vector.magnitude(rvod1[0])) + derivativeMagnitudeStart=elementLengthUpLV, derivativeMagnitudeEnd=magnitude(rvod1[0])) rvix += rvox [1:] rvid1 += rvod1[1:] rvid2 += rvod2[1:] @@ -436,7 +435,7 @@ def generateBaseMesh(cls, region, options): rvod3 = [] for n in range(elementsCountAroundFull + 1): proportion1, proportion2 = rviProportions[n] - cuspDirection = vector.setMagnitude(rvid2[n], -1.0) + cuspDirection = set_magnitude(rvid2[n], -1.0) nRadians = n*math.pi/elementsCountAroundFull cuspAngleRadians = math.sin(nRadians)*math.sin(nRadians)*rvApexCuspAngleRadians + \ math.cos(nRadians)*math.cos(nRadians)*(rvInletCuspAngleRadians if (n <= elementsCountAroundHalf) else rvOutletCuspAngleRadians) @@ -451,19 +450,19 @@ def generateBaseMesh(cls, region, options): rvod2.append(rd2[1]) rvod3.append(rd3[1]) length = sum(math.sqrt(sum(s*s for s in [ rvox[i + 1][c] - rvox[i][c] for c in range(3) ])) for i in range(elementsCountAroundFull))/elementsCountAroundFull - rvid2 = smoothCubicHermiteDerivativesLine(rvox, [ vector.setMagnitude(rvid2[i], length) for i in range(len(rvid2)) ], fixAllDirections = True) + rvid2 = smoothCubicHermiteDerivativesLine(rvox, [ set_magnitude(rvid2[i], length) for i in range(len(rvid2)) ], fixAllDirections = True) rvod2 = smoothCubicHermiteDerivativesLine(rvox, rvid2, fixAllDirections = True) # get position of centre of RV freewall - rvInletLateralDirection = vector.normalise(vector.crossproduct3(rvInletSurfaceNormal, rvInletDerivative)) - rvOutletLateralDirection = vector.normalise(vector.crossproduct3(rvOutletDerivative, rvOutletSurfaceNormal)) + rvInletLateralDirection = normalize(cross(rvInletSurfaceNormal, rvInletDerivative)) + rvOutletLateralDirection = normalize(cross(rvOutletDerivative, rvOutletSurfaceNormal)) td1 = smoothCubicHermiteDerivativesLine([ rvInletCuspCoordinates, rvOutletCuspCoordinates ], [ rvInletLateralDirection, rvOutletLateralDirection ], fixAllDirections = True) px, pd1, pd2, pd3, pProportions = lvTrackSurface.createHermiteCurvePoints(rvInletPositionAroundLV, baseProportionUp, rvOutletPositionAroundLV, baseProportionUp, elementsCount = 2, derivativeStart = [ 0.5*d for d in td1[0] ], derivativeEnd = [ 0.5*d for d in td1[1] ], curveMode = TrackSurface.HermiteCurveMode.UNIFORM_SIZE) # get regular spacing up centre of RV rbcx = [ (px[1][c] + (rvWidth + rvFreeWallThickness)*pd3[1][c]) for c in range(3) ] rcx = [ rvax[1], rbcx ] - rcd1 = [ vector.setMagnitude(rvad1[1], -1.0), pd2[1] ] + rcd1 = [ set_magnitude(rvad1[1], -1.0), pd2[1] ] #print("rcx ", rcx) #print("rcd1", rcd1) rcd1 = smoothCubicHermiteDerivativesLine(rcx, rcd1, fixAllDirections = True) @@ -475,22 +474,22 @@ def generateBaseMesh(cls, region, options): rscd3 = [] for n in range(len(rscx)): a = [ (rscx[n][c] - centre[c]) for c in range(3) ] - d1 = vector.normalise(vector.crossproduct3(rscd2[n], a)) - d3 = vector.normalise(vector.crossproduct3(d1, rscd2[n])) + d1 = normalize(cross(rscd2[n], a)) + d3 = normalize(cross(d1, rscd2[n])) rscd1.append(d1) rscd3.append(d3) # transition from RV apex to LV apex: lax, lad2, lad1, lad3, laProportions = cls.getTransitionRVtoLV(lvTrackSurface, elementsCountUpLVApex + 1, - [ rvApexProportion1, rvApexProportion2 ], vector.setMagnitude(rvApexCuspDirection, ivSulcusApexTransitionLength), + [ rvApexProportion1, rvApexProportion2 ], set_magnitude(rvApexCuspDirection, ivSulcusApexTransitionLength), [ 0.5, 0.0 ], ivSulcusApexTransitionLength, rvx=rscx[0], rvd1=[ -s for s in rscd2[0] ]) # fix apex - lad3[-1] = vector.normalise(poleAxis) - lad1[-1] = vector.crossproduct3(lad3[-1], lad2[-1]) + lad3[-1] = normalize(poleAxis) + lad1[-1] = cross(lad3[-1], lad2[-1]) # reverse d1 lad1 = [ [-s for s in d ] for d in lad1 ] - lad3 = [ vector.setMagnitude(d, lvFreeWallThickness) for d in lad3 ] + lad3 = [ set_magnitude(d, lvFreeWallThickness) for d in lad3 ] rvShield = ShieldMesh2D(elementsCountAroundRVFreeWall, elementsCountUpRVFreeWall, 0) rvx = rvShield.px @@ -568,7 +567,7 @@ def generateBaseMesh(cls, region, options): for n2 in range(elementsCountUpRVFreeWall + 1): for n1 in range(elementsCountAroundRVFreeWall + 1): if rvd1[1][n2][n1]: - rvd3[0][n2][n1] = rvd3[1][n2][n1] = vector.setMagnitude(vector.crossproduct3(rvd1[1][n2][n1], rvd2[1][n2][n1]), rvFreeWallThickness) + rvd3[0][n2][n1] = rvd3[1][n2][n1] = set_magnitude(cross(rvd1[1][n2][n1], rvd2[1][n2][n1]), rvFreeWallThickness) rvx [0][n2][n1] = [ (rvx [1][n2][n1][c] - rvd3[1][n2][n1][c]) for c in range(3) ] # get inner d1, d2 @@ -621,7 +620,7 @@ def generateBaseMesh(cls, region, options): # sample around top of LV free wall to get centre sulcusDerivativeFactor = 2.0 # GRC fudge factor - td1 = [ vector.setMagnitude(d, sulcusDerivativeFactor*ivSulcusBaseTransitionLength) for d in [ rvOutletLateralDirection, rvInletLateralDirection ] ] + td1 = [ set_magnitude(d, sulcusDerivativeFactor*ivSulcusBaseTransitionLength) for d in [ rvOutletLateralDirection, rvInletLateralDirection ] ] px, pd1, pd2, pd3, pProportions = lvTrackSurface.createHermiteCurvePoints(rvOutletPositionAroundLV, baseProportionUp, rvInletPositionAroundLV + 1.0, baseProportionUp, elementsCount = 2, derivativeStart = td1[0], derivativeEnd = td1[1], curveMode = TrackSurface.HermiteCurveMode.SMOOTH) # sample up centre of LV free wall @@ -632,15 +631,15 @@ def generateBaseMesh(cls, region, options): td2[i] = [ d/elementsCountUpLVFreeWall for d in td2[i] ] lscx, lscd2, lscd1, lscd3, lscProportions = lvTrackSurface.createHermiteCurvePoints(1.0, 0.0, pProportions[1][0], pProportions[1][1], elementsCount = elementsCountUpLVFreeWall, derivativeStart = lad2[-1], derivativeEnd = td2[1], curveMode = TrackSurface.HermiteCurveMode.TRANSITION_START) - lscd3[0] = vector.normalise(vector.crossproduct3(nd2[0], nd2[elementsCountAroundLVTrackSurface*3//4])) - lscd1 = [ vector.crossproduct3(lscd2[i], lscd3[i]) for i in range(len(lscd2)) ] + lscd3[0] = normalize(cross(nd2[0], nd2[elementsCountAroundLVTrackSurface*3//4])) + lscd1 = [ cross(lscd2[i], lscd3[i]) for i in range(len(lscd2)) ] # around regular rows of LV free wall, transition from RV for n2 in range(2 + elementsCountUpLVApex, elementsCountUpLV + 1): ix = elementsCountUpLV - n2 ox = -1 - ix td1 = smoothCubicHermiteDerivativesLine([ lscx[ox], rvix[ix] ], - [ lscd1[ox], vector.setMagnitude(rvid1[ix], -sulcusDerivativeFactor*ivSulcusBaseTransitionLength) ], fixStartDirection=True, fixEndDerivative=True) + [ lscd1[ox], set_magnitude(rvid1[ix], -sulcusDerivativeFactor*ivSulcusBaseTransitionLength) ], fixStartDirection=True, fixEndDerivative=True) td1 = [ [ 0.25*s for s in d ] for d in td1 ] px, pd1, pd2, pd3, pProportions = lvTrackSurface.createHermiteCurvePoints( lscProportions[ox][0], lscProportions[ox][1], @@ -650,7 +649,7 @@ def generateBaseMesh(cls, region, options): qx, qd1, qd2, qd3, qProportions = lvTrackSurface.createHermiteCurvePoints( rviProportions[ox][0], rviProportions[ox][1], lscProportions[ox][0], lscProportions[ox][1], - elementsCount = 4, derivativeStart = vector.setMagnitude(rvid1[ox], 0.25*sulcusDerivativeFactor*ivSulcusBaseTransitionLength), derivativeEnd = td1[0], + elementsCount = 4, derivativeStart = set_magnitude(rvid1[ox], 0.25*sulcusDerivativeFactor*ivSulcusBaseTransitionLength), derivativeEnd = td1[0], curveMode = TrackSurface.HermiteCurveMode.UNIFORM_SIZE) qx += px [1:] qd1 += pd1[1:] @@ -680,8 +679,8 @@ def generateBaseMesh(cls, region, options): rx [-1] = rvx[1][0][hx] rd1[-1] = rvd1[1][0][hx] tx, td1 = sampleCubicHermiteCurves(rx, rd1, elementsCountAroundLVFreeWall + 2, - addLengthStart=0.5*vector.magnitude(rd1[0]), lengthFractionStart=0.5, - addLengthEnd=0.5*vector.magnitude(rd1[-1]), lengthFractionEnd=0.5, + addLengthStart=0.5*magnitude(rd1[0]), lengthFractionStart=0.5, + addLengthEnd=0.5*magnitude(rd1[-1]), lengthFractionEnd=0.5, arcLengthDerivatives=True)[0:2] td2 = [ None ]*(elementsCountAroundLVFreeWall + 3) td3 = [ None ]*(elementsCountAroundLVFreeWall + 3) @@ -690,9 +689,9 @@ def generateBaseMesh(cls, region, options): for i in range(1, elementsCountAroundLVFreeWall + 2): p1 = lvTrackSurface.findNearestPosition(tx[i], p1) x, sd1, sd2 = lvTrackSurface.evaluateCoordinates(p1, derivatives = True) - d1, d2, d3 = calculate_surface_axes(sd1, sd2, vector.normalise(td1[i])) - td2[i] = vector.setMagnitude(d2, vector.magnitude(td1[i])) - td3[i] = vector.setMagnitude(d3, lvFreeWallThickness) + d1, d2, d3 = calculate_surface_axes(sd1, sd2, normalize(td1[i])) + td2[i] = set_magnitude(d2, magnitude(td1[i])) + td3[i] = set_magnitude(d3, lvFreeWallThickness) tProportions[i] = lvTrackSurface.getProportion(p1) lvx [1][ox] = tx [1:-1] lvd1[1][ox] = td1[1:-1] @@ -730,8 +729,8 @@ def generateBaseMesh(cls, region, options): derivativeStart=[ -d for d in lvd2[1][n2reg][n1reg] ], derivativeEnd=lad1[lan]) tx, td1, td2, td3, tProportions = lvTrackSurface.resampleHermiteCurvePointsSmooth(tx, td1, td2, td3, tProportions, - derivativeMagnitudeStart=vector.magnitude(lvd2[1][n2reg][n1reg]), - derivativeMagnitudeEnd=None) # vector.magnitude(lad1[lan])) + derivativeMagnitudeStart=magnitude(lvd2[1][n2reg][n1reg]), + derivativeMagnitudeEnd=None) # magnitude(lad1[lan])) # substitute apex derivatives: td2[-1] = lad2[lan] td3[-1] = lad3[lan] @@ -748,12 +747,12 @@ def generateBaseMesh(cls, region, options): derivativeStart=lad1[lan], derivativeEnd=lvd2[1][n2reg][n1reg]) ux, ud1, ud2, ud3, uProportions = lvTrackSurface.resampleHermiteCurvePointsSmooth(ux, ud1, ud2, ud3, uProportions, - derivativeMagnitudeStart=vector.magnitude(td1[-1]), # vector.magnitude(lad1[lan]), - derivativeMagnitudeEnd=vector.magnitude(lvd2[1][n2reg][n1reg])) + derivativeMagnitudeStart=magnitude(td1[-1]), # magnitude(lad1[lan]), + derivativeMagnitudeEnd=magnitude(lvd2[1][n2reg][n1reg])) lvx [1][n2][n1b:m1a] = tx [1:] + ux [1:-1] lvd1[1][n2][n1b:m1a] = td1[1:] + ud1[1:-1] lvd2[1][n2][n1b:m1a] = td2[1:] + ud2[1:-1] - lvd3[1][n2][n1b:m1a] = [ vector.setMagnitude(d, lvFreeWallThickness) for d in (td3[1:] + ud3[1:-1]) ] + lvd3[1][n2][n1b:m1a] = [ set_magnitude(d, lvFreeWallThickness) for d in (td3[1:] + ud3[1:-1]) ] lvProportions[n2][n1b:m1a] = tProportions[1:] + uProportions[1:-1] # up regular columns of LV @@ -767,7 +766,7 @@ def generateBaseMesh(cls, region, options): lvx [1][n2b][n1] = tx [1] lvd1[1][n2b][n1] = [ -d for d in td1[1] ] lvd2[1][n2b][n1] = td2[1] - lvd3[1][n2b][n1] = vector.setMagnitude(td3[1], lvFreeWallThickness) + lvd3[1][n2b][n1] = set_magnitude(td3[1], lvFreeWallThickness) lvProportions[n2b][n1] = tProportions[1] lvShield.getTriplePoints(n3=1) @@ -795,7 +794,7 @@ def generateBaseMesh(cls, region, options): for n1 in range(elementsCountAroundLVFreeWall + 1): if lvd1[1][n2][n1]: d3 = lvd3[1][n2][n1] - lvd3[0][n2][n1] = lvd3[1][n2][n1] = vector.setMagnitude(d3, lvFreeWallThickness) + lvd3[0][n2][n1] = lvd3[1][n2][n1] = set_magnitude(d3, lvFreeWallThickness) lvx [0][n2][n1] = [ (lvx [1][n2][n1][c] - lvd3[1][n2][n1][c]) for c in range(3) ] # get inner d1 (-/+d2 on regular rows) around full rows of rim up LV apex diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_heartventriclesbase1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_heartventriclesbase1.py index 2488ff28..73c263cc 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_heartventriclesbase1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_heartventriclesbase1.py @@ -21,7 +21,6 @@ from scaffoldmaker.meshtypes.meshtype_3d_heartventricles1 import MeshType_3d_heartventricles1 from scaffoldmaker.meshtypes.scaffold_base import Scaffold_base from scaffoldmaker.utils import interpolation as interp -from scaffoldmaker.utils import vector from scaffoldmaker.utils.eft_utils import remapEftLocalNodes, remapEftNodeValueLabel, scaleEftNodeValueLabels, setEftScaleFactorIds from scaffoldmaker.utils.eftfactory_tricubichermite import eftfactory_tricubichermite from scaffoldmaker.utils.geometry import createCirclePoints @@ -410,9 +409,9 @@ def generateBaseMesh(cls, region, options): axis1 = [ 0.0, -cosLvOutletFrontInclineRadians, sinLvOutletFrontInclineRadians ] axis2 = [ 1.0, 0.0, 0.0 ] lvOutletInnerx, lvOutletInnerd1 = createCirclePoints(lvOutletCentre, - vector.setMagnitude(axis1, lvOutletInnerRadius), vector.setMagnitude(axis2, lvOutletInnerRadius), elementsCountAroundOutlet) + set_magnitude(axis1, lvOutletInnerRadius), set_magnitude(axis2, lvOutletInnerRadius), elementsCountAroundOutlet) lvOutletOuterx, lvOutletOuterd1 = createCirclePoints(lvOutletCentre, - vector.setMagnitude(axis1, lvOutletOuterRadius), vector.setMagnitude(axis2, lvOutletOuterRadius), elementsCountAroundOutlet) + set_magnitude(axis1, lvOutletOuterRadius), set_magnitude(axis2, lvOutletOuterRadius), elementsCountAroundOutlet) lvOutletd2 = [ 0.0, vOutletElementLength*sinLvOutletFrontInclineRadians, vOutletElementLength*cosLvOutletFrontInclineRadians ] zero = [ 0.0, 0.0, 0.0 ] lvOutletOuterd3 = [ None ]*elementsCountAroundOutlet @@ -424,12 +423,12 @@ def generateBaseMesh(cls, region, options): dz = vOutletSpacingz axis1 = [ 0.0, -cosLvOutletFrontInclineRadians, sinLvOutletFrontInclineRadians ] axis2 = [ cosRvOutletLeftInclineRadians, sinRvOutletLeftInclineRadians*sinLvOutletFrontInclineRadians, sinRvOutletLeftInclineRadians*cosLvOutletFrontInclineRadians ] - axis3 = vector.crossproduct3(axis1, axis2) + axis3 = cross(axis1, axis2) rvOutletCentre = [ (lvOutletOuterx[3][c] - dy*axis1[c] + dz*axis3[c]) for c in range(3) ] rvOutletInnerx, rvOutletInnerd1 = createCirclePoints(rvOutletCentre, - vector.setMagnitude(axis1, rvOutletInnerRadius), vector.setMagnitude(axis2, rvOutletInnerRadius), elementsCountAroundOutlet) + set_magnitude(axis1, rvOutletInnerRadius), set_magnitude(axis2, rvOutletInnerRadius), elementsCountAroundOutlet) rvOutletOuterx, rvOutletOuterd1 = createCirclePoints(rvOutletCentre, - vector.setMagnitude(axis1, rvOutletOuterRadius), vector.setMagnitude(axis2, rvOutletOuterRadius), elementsCountAroundOutlet) + set_magnitude(axis1, rvOutletOuterRadius), set_magnitude(axis2, rvOutletOuterRadius), elementsCountAroundOutlet) rvOutletd2 = [ vOutletElementLength*axis3[c] for c in range(3) ] # fix derivative 3 on lv outlet adjacent to rv outlet @@ -543,15 +542,15 @@ def generateBaseMesh(cls, region, options): noa = elementsCountAroundRightAtriumFreeWall - 1 nov = -elementsCountAroundAtrialSeptum - 1 mag = baseHeight + baseThickness - d2 = vector.setMagnitude(vector.crossproduct3(ravd3[0][0][noa], ravd1[0][0][noa]), mag) + d2 = set_magnitude(cross(ravd3[0][0][noa], ravd1[0][0][noa]), mag) pd2 = interp.smoothCubicHermiteDerivativesLine([ rvInnerx[nov], ravx[0][0][noa]], [ rvInnerd2[nov], d2 ], fixStartDerivative=True, fixEndDirection=True) ravd2[0][0][noa] = pd2[1] # set d2 at ra node mid supraventricular crest to be normal to surface; smooth to get final magnitude later ravsvcn1 = elementsCountAroundRightAtriumFreeWall - 2 mag = baseHeight + baseThickness - ravd2[0][0][ravsvcn1] = vector.setMagnitude(vector.crossproduct3(ravd3[0][0][ravsvcn1], ravd1[0][0][ravsvcn1]), mag) - ravd2[1][0][ravsvcn1] = vector.setMagnitude(vector.crossproduct3(ravd3[1][0][ravsvcn1], ravd1[1][0][ravsvcn1]), mag) + ravd2[0][0][ravsvcn1] = set_magnitude(cross(ravd3[0][0][ravsvcn1], ravd1[0][0][ravsvcn1]), mag) + ravd2[1][0][ravsvcn1] = set_magnitude(cross(ravd3[1][0][ravsvcn1], ravd1[1][0][ravsvcn1]), mag) ravsvcn2 = elementsCountAroundRightAtriumFreeWall - 3 # copy derivative 3 from av points to LV outlet at centre, left and right cfb; negate as d1 is reversed: @@ -562,8 +561,8 @@ def generateBaseMesh(cls, region, options): # create point above anterior ventricular septum end xi = 0.5 fd2 = [ (lvOutletOuterx[4][c] - vOuterx[1][c]) for c in range(3) ] - mag = 1.0*vector.magnitude(fd2) - fd2 = vector.setMagnitude([ fd2[0], fd2[1], 0.0 ], mag) + mag = 1.0*magnitude(fd2) + fd2 = set_magnitude([ fd2[0], fd2[1], 0.0 ], mag) x = interp.interpolateCubicHermite(vOuterx[0], vOuterd2[0], lvOutletOuterx[4], fd2, xi) d2 = interp.interpolateCubicHermiteDerivative(vOuterx[0], vOuterd2[0], lvOutletOuterx[4], fd2, xi) pd2 = interp.smoothCubicHermiteDerivativesLine([ vOuterx[0], x, lvOutletOuterx[4] ], [ vOuterd2[0], d2, lvOutletd2 ], fixAllDirections=True, fixStartDerivative=True, fixEndDerivative=True) @@ -581,11 +580,11 @@ def generateBaseMesh(cls, region, options): sx = vOuterx[ns] sd2 = vOuterd2[ns] fx = lvOutletOuterx[nf] - fd2 = vector.setMagnitude([lvOutletInnerx[nf][c] - lvOutletOuterx[nf][c] for c in range(3)], - vector.magnitude(sd2)) + fd2 = set_magnitude([lvOutletInnerx[nf][c] - lvOutletOuterx[nf][c] for c in range(3)], + magnitude(sd2)) scale = interp.computeCubicHermiteDerivativeScaling(sx, sd2, fx, fd2) px, pd2 = interp.sampleCubicHermiteCurvesSmooth([sx, fx], [[d*scale for d in sd2], [d*scale for d in fd2]], 2, - derivativeMagnitudeStart=vector.magnitude(sd2))[0:2] + derivativeMagnitudeStart=magnitude(sd2))[0:2] svcox = px[1] sd1 = [-d for d in ravd2[1][0][ravsvcn2]] fd1 = [rvOutletOuterd1[2][c] + rvOutletd2[c] for c in range(3)] @@ -594,7 +593,7 @@ def generateBaseMesh(cls, region, options): fixStartDerivative=True, fixEndDerivative=True) svcod1 = pd1[1] svcod2 = pd2[1] - svcod3 = vector.setMagnitude(vector.crossproduct3(svcod1, svcod2), baseThickness) + svcod3 = set_magnitude(cross(svcod1, svcod2), baseThickness) lvOutletOuterd3[nf] = [-d for d in pd2[2]] # set a reasonable value for next d2 up on right fibrous ring by supraventricular crest 1 sd12 = [(svcod2[c] - svcod1[c]) for c in range(3)] diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_heartventriclesbase2.py b/src/scaffoldmaker/meshtypes/meshtype_3d_heartventriclesbase2.py index e01b35e7..0c58a117 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_heartventriclesbase2.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_heartventriclesbase2.py @@ -18,7 +18,6 @@ from scaffoldmaker.meshtypes.meshtype_3d_heartventricles2 import MeshType_3d_heartventricles2 from scaffoldmaker.meshtypes.scaffold_base import Scaffold_base from scaffoldmaker.utils import interpolation as interp -from scaffoldmaker.utils import vector from scaffoldmaker.utils.eft_utils import remapEftLocalNodes, remapEftNodeValueLabel, scaleEftNodeValueLabels, setEftScaleFactorIds from scaffoldmaker.utils.eftfactory_bicubichermitelinear import eftfactory_bicubichermitelinear from scaffoldmaker.utils.eftfactory_tricubichermite import eftfactory_tricubichermite @@ -253,7 +252,7 @@ def generateBaseMesh(cls, region, options): bx = [ 0.5*(ax[i] + bx[i]) for i in range(2) ] bx.append(ax[2]) ax = [ (bx[c] - px[c]) for c in range(3) ] - ax = vector.normalise(ax) + ax = normalize(ax) baseRotationRadians = math.atan2(ax[1], ax[0]) # get crux location outletSpacingRadians = 0.25*math.pi # GRC make option? @@ -284,8 +283,8 @@ def generateBaseMesh(cls, region, options): radius = lvOutletInnerRadius if (n3 == 0) else lvOutletOuterRadius loAxis1 = [ radius*ax[c] for c in range(3) ] loAxis2 = [ -loAxis1[1]*cosOutletInclineRadians, loAxis1[0]*cosOutletInclineRadians, -radius*sinOutletInclineRadians ] - loAxis3 = vector.crossproduct3(loAxis1, loAxis2) - scale = vOutletElementLength/vector.magnitude(loAxis3) + loAxis3 = cross(loAxis1, loAxis2) + scale = vOutletElementLength/magnitude(loAxis3) dx_ds2 = [ v*scale for v in loAxis3 ] outletNodeId = [] for n1 in range(elementsCountAroundOutlet): @@ -309,7 +308,7 @@ def generateBaseMesh(cls, region, options): dx_ds3[2] = -dx_ds3[2] else: dx_ds3[2] = -2.0*dx_ds3[2] - scale = radiansPerElementAroundOutlet*rvOutletOuterRadius/vector.magnitude(dx_ds3) + scale = radiansPerElementAroundOutlet*rvOutletOuterRadius/magnitude(dx_ds3) dx_ds3 = [ d*scale for d in dx_ds3 ] coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS3, 1, dx_ds3) if n1 == 0: @@ -330,7 +329,7 @@ def generateBaseMesh(cls, region, options): outletCentreSpacing = lvOutletOuterRadius + outletSpacingHorizontal + rvOutletOuterRadius rvOutletCentre = [ (lvOutletCentre[c] - outletCentreSpacing*ax[c]) for c in range(3) ] # add outletSpacingVertical rotated by vOutletInclineRadians - unitCrossX = vector.normalise([-ax[1], ax[0]]) + unitCrossX = normalize([-ax[1], ax[0]]) rvOutletCentre[0] -= outletSpacingVertical*sinOutletInclineRadians*unitCrossX[0] rvOutletCentre[1] -= outletSpacingVertical*sinOutletInclineRadians*unitCrossX[1] rvOutletCentre[2] += outletSpacingVertical*cosOutletInclineRadians @@ -340,8 +339,8 @@ def generateBaseMesh(cls, region, options): radius = rvOutletInnerRadius if (n3 == 0) else rvOutletOuterRadius roAxis1 = [ radius*ax[c] for c in range(3) ] roAxis2 = [ -roAxis1[1]*cosOutletInclineRadians, roAxis1[0]*cosOutletInclineRadians, radius*sinOutletInclineRadians ] - roAxis3 = vector.crossproduct3(roAxis1, roAxis2) - scale = vOutletElementLength/vector.magnitude(roAxis3) + roAxis3 = cross(roAxis1, roAxis2) + scale = vOutletElementLength/magnitude(roAxis3) dx_ds2 = [ v*scale for v in roAxis3 ] outletNodeId = [] for n1 in range(elementsCountAroundOutlet): @@ -368,7 +367,7 @@ def generateBaseMesh(cls, region, options): dx_ds3[2] = 4.0*dx_ds3[2] dx_ds3 = [ (dx_ds1[c] + dx_ds3[c]) for c in range(3) ] mag3 = radiansPerElementAroundOutlet*rvOutletOuterRadius - scale = mag3/vector.magnitude(dx_ds3) + scale = mag3/magnitude(dx_ds3) dx_ds3 = [ d*scale for d in dx_ds3 ] coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS3, 1, dx_ds3) if n1 == 0: @@ -542,7 +541,7 @@ def generateBaseMesh(cls, region, options): -sinRadiansAround*laOuterMajor[1] + cosRadiansAround*laOuterMinor[1], 0.0 ] scale1 = laOuterDerivatives[n1] - scale1 /= vector.magnitude(dx_ds1) + scale1 /= magnitude(dx_ds1) dx_ds1 = [ d*scale1 for d in dx_ds1 ] dx_ds3 = [ outer[0] - inner[0], outer[1] - inner[1], outer[2] - inner[2] ] if (n1 < lan1CruxLimit) or (n1 > lan1SeptumLimit): @@ -559,7 +558,7 @@ def generateBaseMesh(cls, region, options): else: # GRC check scaling here: mag2 = inner[2] - scale2 = mag2/vector.magnitude(dx_ds2) + scale2 = mag2/magnitude(dx_ds2) dx_ds2 = [ d*scale2 for d in dx_ds2 ] coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 1, dx_ds1) coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, dx_ds2) @@ -623,7 +622,7 @@ def generateBaseMesh(cls, region, options): -sinRadiansAround*raOuterMajor[1] + cosRadiansAround*raOuterMinor[1], 0.0 ] scale1 = raOuterDerivatives[n1] - scale1 /= vector.magnitude(dx_ds1) + scale1 /= magnitude(dx_ds1) dx_ds1 = [ d*scale1 for d in dx_ds1 ] dx_ds3 = [ outer[0] - inner[0], outer[1] - inner[1], outer[2] - inner[2] ] if (n1 <= ran1SeptumLimit) or (n1 >= ran1CruxLimit): @@ -640,7 +639,7 @@ def generateBaseMesh(cls, region, options): else: # GRC check scaling here: mag2 = inner[2] - scale2 = mag2/vector.magnitude(dx_ds2) + scale2 = mag2/magnitude(dx_ds2) dx_ds2 = [ d*scale2 for d in dx_ds2 ] coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 1, dx_ds1) coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, dx_ds2) @@ -705,7 +704,7 @@ def generateBaseMesh(cls, region, options): d1 = [ (x1[c] - xc[c]) for c in range(3) ] d2 = [ (x2[c] - xc[c]) for c in range(3) ] dx_ds3 = [ d1[0] + d2[0], d1[1] + d2[1], d1[2] + d2[2] ] - scale = vector.magnitude(d1)/vector.magnitude(dx_ds3) + scale = magnitude(d1)/magnitude(dx_ds3) dx_ds3 = [ d*scale for d in dx_ds3 ] result = coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS3, 1, dx_ds3 ) @@ -713,8 +712,8 @@ def generateBaseMesh(cls, region, options): cache.setNode(nodes.findNodeByIdentifier(laNodeId[0][2])) result, dx_ds1 = coordinates.getNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 1, 3 ) result, dx_ds3 = coordinates.getNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS3, 1, 3 ) - dx_ds2 = vector.crossproduct3(dx_ds3, dx_ds1) - scale2 = 0.5*vector.magnitude(dx_ds3)/vector.magnitude(dx_ds2) + dx_ds2 = cross(dx_ds3, dx_ds1) + scale2 = 0.5*magnitude(dx_ds3)/magnitude(dx_ds2) dx_ds2 = [ scale2*d for d in dx_ds2 ] result = coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, dx_ds2 ) @@ -752,7 +751,7 @@ def generateBaseMesh(cls, region, options): result, d1b = coordinates.getNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 1, 3 ) result, d2b = coordinates.getNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS3, 1, 3 ) d2b = [ -2.0*d for d in d2b ] - scale = 4.0*(baseHeight + baseThickness)/vector.magnitude(d2a) + scale = 4.0*(baseHeight + baseThickness)/magnitude(d2a) d2a = [ scale*d for d in d2a ] xi = 0.5 xr = 1.0 - xi @@ -760,7 +759,7 @@ def generateBaseMesh(cls, region, options): dx_ds1 = [ (xr*d1a[c] + xi*d1b[c]) for c in range(3) ] dx_ds2 = interp.interpolateCubicHermiteDerivative(xa, d2a, xb, d2b, xi) dx_ds2 = [ xr*d for d in dx_ds2 ] - radialVector = vector.normalise(vector.crossproduct3(dx_ds1, dx_ds2)) + radialVector = normalize(cross(dx_ds1, dx_ds2)) dx_ds3 = [ baseThickness*d for d in radialVector ] x_inner = [ (x[c] - dx_ds3[c]) for c in range(3) ] diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_ostium1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_ostium1.py index 5d058351..4d14e747 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_ostium1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_ostium1.py @@ -12,9 +12,7 @@ from cmlibs.zinc.node import Node from scaffoldmaker.meshtypes.scaffold_base import Scaffold_base from scaffoldmaker.utils import interpolation as interp -from scaffoldmaker.utils import vector from scaffoldmaker.utils.annulusmesh import createAnnulusMesh3d -from scaffoldmaker.utils.eftfactory_tricubichermite import eftfactory_tricubichermite from scaffoldmaker.utils.geometry import createCirclePoints, getCircleProjectionAxes from scaffoldmaker.utils.meshrefinement import MeshRefinement from scaffoldmaker.utils.tracksurface import TrackSurface, TrackSurfacePosition, calculate_surface_axes @@ -388,9 +386,9 @@ def generateOstiumMesh(region, options, trackSurface, centrePosition, axis1, sta pd2, pd1, pd3 = calculate_surface_axes(d1, d2, sideDirection) # get outer coordinates opx = px - opd1 = vector.setMagnitude([-d for d in pd1], elementLengthAroundOstiumEnd) - opd2 = vector.setMagnitude(pd2, elementLengthAroundOstiumEnd) # smoothed later - opd3 = vector.setMagnitude(pd3, ostiumWallThickness) + opd1 = set_magnitude([-d for d in pd1], elementLengthAroundOstiumEnd) + opd2 = set_magnitude(pd2, elementLengthAroundOstiumEnd) # smoothed later + opd3 = set_magnitude(pd3, ostiumWallThickness) ostiumTotalXi3 = 0.0 vesselTotalXi3 = 0.0 @@ -459,13 +457,13 @@ def generateOstiumMesh(region, options, trackSurface, centrePosition, axis1, sta mx, d1, d2 = trackSurface.evaluateCoordinates(position, derivatives=True) md1, md2, md3 = calculate_surface_axes(d1, d2, trackDirection1) nx .insert(1, [(mx[c] + interVesselHeight * md3[c]) for c in range(3)]) - nd1.insert(1, vector.setMagnitude(md1, elementLengthAroundOstiumMid if (0 < iv < (vesselsCount - 2)) else + nd1.insert(1, set_magnitude(md1, elementLengthAroundOstiumMid if (0 < iv < (vesselsCount - 2)) else elementLengthAroundOstiumEnd)) - nd2.insert(1, vector.setMagnitude(md2, ostiumRadius)) + nd2.insert(1, set_magnitude(md2, ostiumRadius)) nd2 = interp.smoothCubicHermiteDerivativesLine(nx, nd2, fixAllDirections=True) px, pd2, pe, pxi = interp.sampleCubicHermiteCurves(nx, nd2, elementsCountAcross)[0:4] pd1 = interp.interpolateSampleLinear(nd1, pe, pxi) - pd3 = [vector.setMagnitude(vector.crossproduct3(pd1[n2], pd2[n2]), commonOstiumWallThickness) + pd3 = [set_magnitude(cross(pd1[n2], pd2[n2]), commonOstiumWallThickness) for n2 in range(elementsCountAcross + 1)] for n3 in range(elementsCountThroughWall + 1): xi3 = 1 - commonOstiumWallThicknessXi3List[n3] @@ -478,7 +476,7 @@ def generateOstiumMesh(region, options, trackSurface, centrePosition, axis1, sta od2[n3][oa] = [-d for d in ld2[0]] od2[n3][ob] = ld2[-1] if useCubicHermiteThroughOstiumWall: - pd3Element = [vector.setMagnitude(vector.crossproduct3(pd1[n2], pd2[n2]), + pd3Element = [set_magnitude(cross(pd1[n2], pd2[n2]), commonOstiumWallThickness * commonOstiumWallThicknessProportions[n3]) for n2 in range(elementsCountAcross + 1)] xd3[iv][n3] = copy.deepcopy(pd3Element[1:elementsCountAcross]) @@ -513,8 +511,8 @@ def generateOstiumMesh(region, options, trackSurface, centrePosition, axis1, sta vod3.append([]) for n3 in range(elementsCountThroughWall + 1): radius = vesselInnerRadius + vesselWallThicknessXi3List[n3] * vesselWallThickness - vAxis1 = vector.setMagnitude(vd1, radius) - vAxis2 = vector.setMagnitude(vd2, radius) + vAxis1 = set_magnitude(vd1, radius) + vAxis2 = set_magnitude(vd2, radius) if vesselsCount == 1: startRadians = 0.5 * math.pi else: @@ -526,7 +524,7 @@ def generateOstiumMesh(region, options, trackSurface, centrePosition, axis1, sta vod1[-1].append(pd1) vod2[-1].append([vd3] * elementsCountAroundVessel) if useCubicHermiteThroughVesselWall: - vod3[-1].append([vector.setMagnitude(vector.crossproduct3(d1, vd3), + vod3[-1].append([set_magnitude(cross(d1, vd3), vesselWallThickness * vesselWallThicknessProportions[n3]) for d1 in pd1]) @@ -668,28 +666,28 @@ def generateOstiumMesh(region, options, trackSurface, centrePosition, axis1, sta for c in range(3): mvPointsd1[v][n3][n1][c] = 0.5 * (mvPointsd1[v][n3][n1][c] + sf1 * ndf[c]) else: - # print('v', v, 'n3', n3, 'n1', n1, ':', vector.magnitude(ndf), 'vs.', vector.magnitude(nd2[1]), + # print('v', v, 'n3', n3, 'n1', n1, ':', magnitude(ndf), 'vs.', magnitude(nd2[1]), # 'd2Map', d2Map) pass # calculate inner d2 derivatives around ostium from outer using track surface curvature factor = 1.0 for n1 in range(elementsCountAroundOstium): - trackDirection = vector.normalise(od2[elementsCountThroughWall][n1]) - trackDistance = factor * vector.magnitude(od2[elementsCountThroughWall][n1]) + trackDirection = normalize(od2[elementsCountThroughWall][n1]) + trackDistance = factor * magnitude(od2[elementsCountThroughWall][n1]) tx = [None, ox[elementsCountThroughWall][n1], None] - td1 = [None, vector.setMagnitude(od2[elementsCountThroughWall][n1], trackDistance), None] - td2 = [None, vector.setMagnitude(od1[elementsCountThroughWall][n1], -trackDistance), None] + td1 = [None, set_magnitude(od2[elementsCountThroughWall][n1], trackDistance), None] + td2 = [None, set_magnitude(od1[elementsCountThroughWall][n1], -trackDistance), None] positionBackward = trackSurface.trackVector(oPositions[n1], trackDirection, -trackDistance) tx[0], d1, d2 = trackSurface.evaluateCoordinates(positionBackward, derivatives=True) sd1, sd2, sd3 = calculate_surface_axes(d1, d2, trackDirection) - td1[0] = vector.setMagnitude(sd1, trackDistance) - td2[0] = vector.setMagnitude(sd2, trackDistance) + td1[0] = set_magnitude(sd1, trackDistance) + td2[0] = set_magnitude(sd2, trackDistance) positionForward = trackSurface.trackVector(oPositions[n1], trackDirection, trackDistance) tx[2], d1, d2 = trackSurface.evaluateCoordinates(positionForward, derivatives=True) sd1, sd2, sd3 = calculate_surface_axes(d1, d2, trackDirection) - td1[2] = vector.setMagnitude(sd1, trackDistance) - td2[2] = vector.setMagnitude(sd2, trackDistance) + td1[2] = set_magnitude(sd1, trackDistance) + td2[2] = set_magnitude(sd2, trackDistance) for n3 in range(elementsCountThroughWall): xi3 = 1 - ostiumWallThicknessXi3List[n3] newd2 = interp.projectHermiteCurvesThroughWall(tx, td1, td2, 1, -ostiumWallThickness * xi3)[1] diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_ostium2.py b/src/scaffoldmaker/meshtypes/meshtype_3d_ostium2.py index 72559d98..57110005 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_ostium2.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_ostium2.py @@ -15,7 +15,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 vector from scaffoldmaker.utils.annulusmesh import createAnnulusMesh3d from scaffoldmaker.utils.eftfactory_tricubichermite import eftfactory_tricubichermite from scaffoldmaker.utils.geometry import sampleEllipsePoints, getApproximateEllipsePerimeter @@ -206,19 +205,19 @@ def generateBaseMesh(cls, region, options): ox = centralPath.cxPath[-1] od2 = centralPath.cd2Path[-1] od3 = centralPath.cd3Path[-1] - scale = 4.0 * vector.magnitude([od2[c] + od3[c] for c in range(3)]) + scale = 4.0 * magnitude([od2[c] + od3[c] for c in range(3)]) - nxOffset = [vector.setMagnitude([-od3[c] - od2[c] for c in range(3)], scale), - vector.setMagnitude([od3[c] - od2[c] for c in range(3)], scale), - vector.setMagnitude([-od3[c] + od2[c] for c in range(3)], scale), - vector.setMagnitude([od3[c] + od2[c] for c in range(3)], scale)] + nxOffset = [set_magnitude([-od3[c] - od2[c] for c in range(3)], scale), + set_magnitude([od3[c] - od2[c] for c in range(3)], scale), + set_magnitude([-od3[c] + od2[c] for c in range(3)], scale), + set_magnitude([od3[c] + od2[c] for c in range(3)], scale)] nx = [] for n in range(len(nxOffset)): nx.append([ox[c] + nxOffset[n][c] for c in range(3)]) - nd1 = [vector.setMagnitude(od3, 2.0 * scale)] * 4 - nd2 = [vector.setMagnitude(od2, 2.0 * scale)] * 4 + nd1 = [set_magnitude(od3, 2.0 * scale)] * 4 + nd2 = [set_magnitude(od2, 2.0 * scale)] * 4 trackSurface = TrackSurface(1, 1, nx, nd1, nd2) generateOstiumMesh(region, options, trackSurface, centralPath) @@ -319,7 +318,7 @@ def generateOstiumMesh(region, options, trackSurface, centralPath, startNodeIden vesselWallThicknessProportionsUI = copy.deepcopy(options['Vessel wall relative thicknesses']) useCubicHermiteThroughVesselWall = not(options['Use linear through vessel wall']) useCrossDerivatives = False # options['Use cross derivatives'] # not implemented - # ostiumRadius = vector.magnitude(centralPath.cd2Path[-1]) + # ostiumRadius = magnitude(centralPath.cd2Path[-1]) arcLength = 0.0 for e in range(len(centralPath.cxPath) - 1): @@ -356,8 +355,8 @@ def generateOstiumMesh(region, options, trackSurface, centralPath, startNodeIden cx, cd1, cd2 = trackSurface.evaluateCoordinates(centrePosition, True) trackDirection1, trackDirection2, centreNormal = calculate_surface_axes(cd1, cd2, cd1) - ellipsePerimeter = getApproximateEllipsePerimeter(vector.magnitude(centralPath.cd2Path[-1]), - vector.magnitude(centralPath.cd3Path[-1])) + ellipsePerimeter = getApproximateEllipsePerimeter(magnitude(centralPath.cd2Path[-1]), + magnitude(centralPath.cd3Path[-1])) # halfCircumference = math.pi * ostiumRadius # circumference = 2.0 * halfCircumference distance = 0.0 @@ -457,10 +456,10 @@ def generateOstiumMesh(region, options, trackSurface, centralPath, startNodeIden px, d1, d2 = trackSurface.evaluateCoordinates(position, True) angleRadians = 2.0 * math.pi / elementsCountAroundOstium * n1 - cosTheta1 = vector.dotproduct(d1, centralPath.cd3Path[-1])/\ - (vector.magnitude(d1) * vector.magnitude(centralPath.cd3Path[-1])) - cosTheta2 = vector.dotproduct(d2, centralPath.cd2Path[-1]) / \ - (vector.magnitude(d2) * vector.magnitude(centralPath.cd2Path[-1])) + cosTheta1 = dot(d1, centralPath.cd3Path[-1])/\ + (magnitude(d1) * magnitude(centralPath.cd3Path[-1])) + cosTheta2 = dot(d2, centralPath.cd2Path[-1]) / \ + (magnitude(d2) * magnitude(centralPath.cd2Path[-1])) signed1 = cosTheta1 / abs(cosTheta1) signed2 = cosTheta2 / abs(cosTheta2) @@ -471,9 +470,9 @@ def generateOstiumMesh(region, options, trackSurface, centralPath, startNodeIden # get inner coordinates opx = px - opd1 = vector.setMagnitude([-d for d in pd1], elementLengthAroundOstiumEnd) - opd2 = vector.setMagnitude(pd2, vector.magnitude(centralPath.cd2Path[0])) # smoothed later - opd3 = vector.setMagnitude(pd3, ostiumWallThickness) + opd1 = set_magnitude([-d for d in pd1], elementLengthAroundOstiumEnd) + opd2 = set_magnitude(pd2, magnitude(centralPath.cd2Path[0])) # smoothed later + opd3 = set_magnitude(pd3, ostiumWallThickness) # node = nodes.createNode(nodeIdentifier, nodetemplate) # cache.setNode(node) @@ -560,13 +559,13 @@ def generateOstiumMesh(region, options, trackSurface, centralPath, startNodeIden # mx, d1, d2 = trackSurface.evaluateCoordinates(position, derivatives=True) # md1, md2, md3 = calculate_surface_axes(d1, d2, trackDirection1) # nx .insert(1, [(mx[c] + interVesselHeight * md3[c]) for c in range(3)]) - # nd1.insert(1, vector.setMagnitude(md1, elementLengthAroundOstiumMid if (0 < iv < (vesselsCount - 2)) else + # nd1.insert(1, set_magnitude(md1, elementLengthAroundOstiumMid if (0 < iv < (vesselsCount - 2)) else # elementLengthAroundOstiumEnd)) - # nd2.insert(1, vector.setMagnitude(md2, ostiumRadius)) + # nd2.insert(1, set_magnitude(md2, ostiumRadius)) # nd2 = interp.smoothCubicHermiteDerivativesLine(nx, nd2, fixAllDirections=True) # px, pd2, pe, pxi = interp.sampleCubicHermiteCurves(nx, nd2, elementsCountAcross)[0:4] # pd1 = interp.interpolateSampleLinear(nd1, pe, pxi) - # pd3 = [vector.setMagnitude(vector.crossproduct3(pd1[n2], pd2[n2]), commonOstiumWallThickness) + # pd3 = [set_magnitude(cross(pd1[n2], pd2[n2]), commonOstiumWallThickness) # for n2 in range(elementsCountAcross + 1)] # for n3 in range(elementsCountThroughWall + 1): # xi3 = 1 - commonOstiumWallThicknessXi3List[n3] @@ -579,7 +578,7 @@ def generateOstiumMesh(region, options, trackSurface, centralPath, startNodeIden # od2[n3][oa] = [-d for d in ld2[0]] # od2[n3][ob] = ld2[-1] # if useCubicHermiteThroughOstiumWall: - # pd3Element = [vector.setMagnitude(vector.crossproduct3(pd1[n2], pd2[n2]), + # pd3Element = [set_magnitude(cross(pd1[n2], pd2[n2]), # commonOstiumWallThickness * commonOstiumWallThicknessProportions[n3]) # for n2 in range(elementsCountAcross + 1)] # xd3[iv][n3] = copy.deepcopy(pd3Element[1:elementsCountAcross]) @@ -630,10 +629,10 @@ def generateOstiumMesh(region, options, trackSurface, centralPath, startNodeIden vod2[-1].append(d2Proximal[n3]) vod3[-1].append(d3Proximal[n3]) else: - radius1 = vector.magnitude(centralPath.cd2Path[0]) - (1.0 - vesselWallThicknessXi3List[n3]) * vesselWallThickness - radius2 = vector.magnitude(centralPath.cd3Path[0]) - (1.0 - vesselWallThicknessXi3List[n3]) * vesselWallThickness - vAxis1 = vector.setMagnitude(centralPath.cd2Path[0], radius1) - vAxis2 = vector.setMagnitude(centralPath.cd3Path[0], radius2) + radius1 = magnitude(centralPath.cd2Path[0]) - (1.0 - vesselWallThicknessXi3List[n3]) * vesselWallThickness + radius2 = magnitude(centralPath.cd3Path[0]) - (1.0 - vesselWallThicknessXi3List[n3]) * vesselWallThickness + vAxis1 = set_magnitude(centralPath.cd2Path[0], radius1) + vAxis2 = set_magnitude(centralPath.cd3Path[0], radius2) # print(centralPath.cd2Path[0], centralPath.cd3Path[0]) # if vesselsCount == 1: startRadians = 0.0 # 0.5 * math.pi @@ -646,10 +645,10 @@ def generateOstiumMesh(region, options, trackSurface, centralPath, startNodeIden vox[-1].append(px) vod1[-1].append(pd1) - vesselEndDerivative = vector.setMagnitude(centralPath.cd1Path[0], arcLength/elementsCountAlong) + vesselEndDerivative = set_magnitude(centralPath.cd1Path[0], arcLength/elementsCountAlong) vod2[-1].append([vesselEndDerivative] * elementsCountAroundVessel) # need to scale with bending? if useCubicHermiteThroughVesselWall: - vod3[-1].append([vector.setMagnitude(vector.crossproduct3(d1, vesselEndDerivative), #cd1Path[-1]), + vod3[-1].append([set_magnitude(cross(d1, vesselEndDerivative), #cd1Path[-1]), vesselWallThickness * vesselWallThicknessProportions[n3]) for d1 in pd1]) @@ -802,28 +801,28 @@ def generateOstiumMesh(region, options, trackSurface, centralPath, startNodeIden for c in range(3): mvPointsd1[v][n3][n1][c] = 0.5 * (mvPointsd1[v][n3][n1][c] + sf1 * ndf[c]) else: - # print('v', v, 'n3', n3, 'n1', n1, ':', vector.magnitude(ndf), 'vs.', vector.magnitude(nd2[1]), + # print('v', v, 'n3', n3, 'n1', n1, ':', magnitude(ndf), 'vs.', magnitude(nd2[1]), # 'd2Map', d2Map) pass # calculate inner d2 derivatives around ostium from outer using track surface curvature factor = 1.0 for n1 in range(elementsCountAroundOstium): - trackDirection = vector.normalise(od2[elementsCountThroughWall][n1]) - trackDistance = factor * vector.magnitude(od2[elementsCountThroughWall][n1]) + trackDirection = normalize(od2[elementsCountThroughWall][n1]) + trackDistance = factor * magnitude(od2[elementsCountThroughWall][n1]) tx = [None, ox[elementsCountThroughWall][n1], None] - td1 = [None, vector.setMagnitude(od2[elementsCountThroughWall][n1], trackDistance), None] - td2 = [None, vector.setMagnitude(od1[elementsCountThroughWall][n1], -trackDistance), None] + td1 = [None, set_magnitude(od2[elementsCountThroughWall][n1], trackDistance), None] + td2 = [None, set_magnitude(od1[elementsCountThroughWall][n1], -trackDistance), None] positionBackward = trackSurface.trackVector(oPositions[n1], trackDirection, -trackDistance) tx[0], d1, d2 = trackSurface.evaluateCoordinates(positionBackward, derivatives=True) sd1, sd2, sd3 = calculate_surface_axes(d1, d2, trackDirection) - td1[0] = vector.setMagnitude(sd1, trackDistance) - td2[0] = vector.setMagnitude(sd2, trackDistance) + td1[0] = set_magnitude(sd1, trackDistance) + td2[0] = set_magnitude(sd2, trackDistance) positionForward = trackSurface.trackVector(oPositions[n1], trackDirection, trackDistance) tx[2], d1, d2 = trackSurface.evaluateCoordinates(positionForward, derivatives=True) sd1, sd2, sd3 = calculate_surface_axes(d1, d2, trackDirection) - td1[2] = vector.setMagnitude(sd1, trackDistance) - td2[2] = vector.setMagnitude(sd2, trackDistance) + td1[2] = set_magnitude(sd1, trackDistance) + td2[2] = set_magnitude(sd2, trackDistance) for n3 in range(elementsCountThroughWall): xi3 = 1 - ostiumWallThicknessXi3List[n3] newd2 = interp.projectHermiteCurvesThroughWall(tx, td1, td2, 1, -ostiumWallThickness * xi3)[1] diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_smallintestine1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_smallintestine1.py index 3ab981b7..83fb7d83 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_smallintestine1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_smallintestine1.py @@ -17,7 +17,6 @@ from scaffoldmaker.scaffoldpackage import ScaffoldPackage from scaffoldmaker.utils import interpolation as interp from scaffoldmaker.utils import tubemesh -from scaffoldmaker.utils import vector from scaffoldmaker.utils.tubemesh import CylindricalSegmentTubeMeshOuterPoints from scaffoldmaker.utils.zinc_utils import exnode_string_from_nodeset_field_parameters, \ get_nodeset_path_field_parameters @@ -1178,7 +1177,7 @@ def createSmallIntestineMesh3d(region, options, networkLayout, nextNodeIdentifie lengthFractionStart=lengthFraction) sd2, sd12 = interp.interpolateSampleCubicHermite(cd2, cd12, se, sxi, ssf) - outerRadiusListCP = [vector.magnitude(c) for c in cd2] + outerRadiusListCP = [magnitude(c) for c in cd2] dOuterRadiusListCP = [] for n in range(len(outerRadiusListCP) - 1): dOuterRadiusListCP.append(outerRadiusListCP[n + 1] - outerRadiusListCP[n]) @@ -1304,8 +1303,8 @@ def createSmallIntestineMesh3d(region, options, networkLayout, nextNodeIdentifie fixEndDerivative = True)[1] d2Extrude.append(d2) d3Unit = \ - vector.normalise(vector.crossproduct3(vector.normalise(d1LastTwoFaces[n1 + elementsCountAround]), - vector.normalise(d2))) + normalize(cross(normalize(d1LastTwoFaces[n1 + elementsCountAround]), + normalize(d2))) d3UnitExtrude.append(d3Unit) d2Extrude = d2Extrude + \ (d2WarpedList[elementsCountAround:-elementsCountAround] if nSegment < segmentCount - 1 else diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_solidsphere1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_solidsphere1.py index 3dbb0b8d..cda172ad 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_solidsphere1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_solidsphere1.py @@ -13,7 +13,6 @@ from cmlibs.zinc.node import Node from scaffoldmaker.meshtypes.scaffold_base import Scaffold_base from scaffoldmaker.utils import interpolation as interp -from scaffoldmaker.utils import vector from scaffoldmaker.utils.eftfactory_tricubichermite import eftfactory_tricubichermite from scaffoldmaker.utils.meshrefinement import MeshRefinement @@ -209,7 +208,7 @@ def generateBaseMesh(cls, region, options): cubicArcLength = cubicArcLengthList[n2] v1 = [0.0, 0.0, -radius+n2*2.0*radius/elementsCountUp] d1 = [math.cos(math.pi/2.0), math.sin(math.pi/2.0), 0.0] - d1 = vector.normalise(d1) + d1 = normalize(d1) d1 = [d*cubicArcLength for d in d1] v2 = [ radius*math.cos(math.pi/2.0)*sinRadiansUp, @@ -217,7 +216,7 @@ def generateBaseMesh(cls, region, options): -radius*cosRadiansUp ] d2 = [math.cos(math.pi/2.0)*sinRadiansUp,math.sin(math.pi/2.0)*sinRadiansUp,-cosRadiansUp] - d2 = vector.normalise(d2) + d2 = normalize(d2) d2 = [d*cubicArcLength for d in d2] x = interp.interpolateCubicHermite(v1, d1, v2, d2, xi) @@ -238,7 +237,7 @@ def generateBaseMesh(cls, region, options): # Calculate node coordinates on arc using cubic hermite interpolation v1 = [0.0, 0.0, -radius+n2*2.0*radius/elementsCountUp] d1 = [cosRadiansAround, sinRadiansAround, 0.0] - d1 = vector.normalise(d1) + d1 = normalize(d1) d1 = [d*cubicArcLength for d in d1] v2 = [ radius*cosRadiansAround*sinRadiansUp, @@ -246,12 +245,12 @@ def generateBaseMesh(cls, region, options): -radius*cosRadiansUp ] d2 = [cosRadiansAround*sinRadiansUp,sinRadiansAround*sinRadiansUp,-cosRadiansUp] - d2 = vector.normalise(d2) + d2 = normalize(d2) d2 = [d*cubicArcLength for d in d2] x = interp.interpolateCubicHermite(v1, d1, v2, d2, xi) # For dx_ds1 - Calculate radius wrt origin where interpolated points lie on - orthoRadius = vector.magnitude(x) + orthoRadius = magnitude(x) orthoRadiansUp = math.pi - math.acos(x[2]/orthoRadius) sinOrthoRadiansUp = math.sin(orthoRadiansUp) cosOrthoRadiansUp = math.cos(orthoRadiansUp) @@ -275,7 +274,7 @@ def generateBaseMesh(cls, region, options): ] dx_ds3 = interp.interpolateCubicHermiteDerivative(v1, d1, v2, d2, xi) - dx_ds3 = vector.normalise(dx_ds3) + dx_ds3 = normalize(dx_ds3) dx_ds3 = [d*cubicArcLength/elementsCountRadial for d in dx_ds3] node = nodes.createNode(nodeIdentifier, nodetemplate) diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_solidsphere2.py b/src/scaffoldmaker/meshtypes/meshtype_3d_solidsphere2.py index fb88a5f5..0b12caf3 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_solidsphere2.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_solidsphere2.py @@ -5,10 +5,10 @@ from __future__ import division +from cmlibs.maths.vectorops import identity_matrix, matrix_vector_mult from cmlibs.utils.zinc.field import findOrCreateFieldCoordinates from scaffoldmaker.annotation.annotationgroup import AnnotationGroup from scaffoldmaker.meshtypes.scaffold_base import Scaffold_base -from scaffoldmaker.utils import vector from scaffoldmaker.utils.meshrefinement import MeshRefinement from scaffoldmaker.utils.spheremesh import SphereMesh, SphereShape @@ -173,12 +173,7 @@ def generateBaseMesh(region, options): annotationGroups = [boxGroup, transitionGroup] centre = [0.0, 0.0, 0.0] - axis1 = [1.0, 0.0, 0.0] - axis2 = [0.0, 1.0, 0.0] - axis3 = [0.0, 0.0, 1.0] - axes = [vector.scaleVector(axis1, radius[0]), - vector.scaleVector(axis2, radius[1]), - vector.scaleVector(axis3, radius[2])] + axes = matrix_vector_mult(identity_matrix(3), radius) elementsCountAcross = [elementsCountAcrossAxis1, elementsCountAcrossAxis2, elementsCountAcrossAxis3] sphere1 = SphereMesh(fm, coordinates, centre, axes, elementsCountAcross, diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_stomach1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_stomach1.py index 073468ad..de11243d 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_stomach1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_stomach1.py @@ -28,7 +28,6 @@ from scaffoldmaker.scaffoldpackage import ScaffoldPackage from scaffoldmaker.utils import interpolation as interp from scaffoldmaker.utils import matrix -from scaffoldmaker.utils import vector from scaffoldmaker.utils.annulusmesh import createAnnulusMesh3d from scaffoldmaker.utils.eftfactory_bicubichermitelinear import eftfactory_bicubichermitelinear from scaffoldmaker.utils.eftfactory_tricubichermite import eftfactory_tricubichermite @@ -1318,7 +1317,7 @@ def createStomachMesh3d(region, fm, coordinates, stomachTermsAlong, allAnnotatio GEJSettings = updateOstiumOptions(options, ostiumOptions) elementsAlongEsophagus = GEJSettings['Number of elements along'] elementsThroughEsophagusWall = GEJSettings['Number of elements through wall'] - ostiumRadius = vector.magnitude(networkLayout.cd2Groups[-1][1]) + ostiumRadius = magnitude(networkLayout.cd2Groups[-1][1]) limitingRidge = options['Limiting ridge'] wallThickness = options['Wall thickness'] GEJSettings['Use linear through ostium wall'] = options['Use linear through wall'] @@ -1454,8 +1453,8 @@ def createStomachMesh3d(region, fm, coordinates, stomachTermsAlong, allAnnotatio elementsOutSection = math.ceil(arcLengthRatioForGroupsFromFundusApex[i]/targetLengthTS) cxSection, cd1Section, pe, pxi, psf = interp.sampleCubicHermiteCurvesSmooth( cxGroup, cd1Group, elementsOutSection, - derivativeMagnitudeStart=None if (i == 1) else 0.0 if (i == 0) else vector.magnitude(cd1Sections[-1][-1]), - derivativeMagnitudeEnd=None if (i != 0) else vector.magnitude(cd1Sections[0][0])) + derivativeMagnitudeStart=None if (i == 1) else 0.0 if (i == 0) else magnitude(cd1Sections[-1][-1]), + derivativeMagnitudeEnd=None if (i != 0) else magnitude(cd1Sections[0][0])) cd2Section = interp.interpolateSampleCubicHermite(cd2Group, cd12Group, pe, pxi, psf)[0] cd3Section = interp.interpolateSampleCubicHermite(cd3Group, cd13Group, pe, pxi, psf)[0] @@ -1499,7 +1498,7 @@ def createStomachMesh3d(region, fm, coordinates, stomachTermsAlong, allAnnotatio d2Apex = cd2PlusSections[0][0] d3Apex = cd3PlusSections[0][0] - rotAxisApex = vector.normalise(vector.crossproduct3(d3Apex, d2Apex)) + rotAxisApex = normalize(cross(d3Apex, d2Apex)) px = sampleEllipsePoints(cxPlusSections[0][0], cd2PlusSections[0][0], cd3PlusSections[0][0], 0.0, math.pi * 2.0, elementsCountAroundDuod)[0] @@ -1553,12 +1552,12 @@ def createStomachMesh3d(region, fm, coordinates, stomachTermsAlong, allAnnotatio # Scale d1 and d2 at apex for n in range(len(xEllipseAroundAll[0])): d1EllipseAroundAll[0][n] = \ - vector.setMagnitude(d1EllipseAroundAll[0][n], + set_magnitude(d1EllipseAroundAll[0][n], interp.computeCubicHermiteArcLength(xEllipseAroundAll[0][n], d2EllipseAroundAll[0][n], xEllipseAroundAll[1][n], d2EllipseAroundAll[1][n], True)) d2EllipseAroundAll[0][n] = \ - vector.setMagnitude(d2EllipseAroundAll[0][n], + set_magnitude(d2EllipseAroundAll[0][n], interp.computeCubicHermiteArcLength(xEllipseAroundAll[0][n], d2EllipseAroundAll[0][n], xEllipseAroundAll[1][n], d2EllipseAroundAll[1][n], True)) @@ -1619,17 +1618,17 @@ def createStomachMesh3d(region, fm, coordinates, stomachTermsAlong, allAnnotatio arcEnd = networkLayout.arcLengthOfGroupsAlong[-1] nearestPosition = trackSurfaceStomach.findNearestPosition(cxEso[0]) xNearestStart = trackSurfaceStomach.evaluateCoordinates(nearestPosition, derivatives=False) - distStart = vector.magnitude([cxEso[0][c] - xNearestStart[c] for c in range(3)]) + distStart = magnitude([cxEso[0][c] - xNearestStart[c] for c in range(3)]) nearestPosition = trackSurfaceStomach.findNearestPosition(cxEso[-1]) xNearestEnd = trackSurfaceStomach.evaluateCoordinates(nearestPosition, derivatives=False) - distEnd = vector.magnitude([cxEso[-1][c] - xNearestEnd[c] for c in range(3)]) + distEnd = magnitude([cxEso[-1][c] - xNearestEnd[c] for c in range(3)]) for iter in range(100): arcDistance = (arcStart + arcEnd) * 0.5 x, d1 = interp.getCubicHermiteCurvesPointAtArcDistance(cxEso, cd1Eso, arcDistance)[0:2] nearestPosition = trackSurfaceStomach.findNearestPosition(x) xNearest = trackSurfaceStomach.evaluateCoordinates(nearestPosition, derivatives=False) - dist = vector.magnitude([x[c] - xNearest[c] for c in range(3)]) + dist = magnitude([x[c] - xNearest[c] for c in range(3)]) if abs(distStart - distEnd) > xTol: if distStart < distEnd: @@ -1641,7 +1640,7 @@ def createStomachMesh3d(region, fm, coordinates, stomachTermsAlong, allAnnotatio else: xCentre, d1Centre, d2Centre = trackSurfaceStomach.evaluateCoordinates(nearestPosition, derivatives=True) - normAxis = vector.normalise([-d for d in d1]) + normAxis = normalize([-d for d in d1]) eIdx = interp.getNearestPointIndex(cxEso, xCentre) - 1 arcLenghtSum = 0.0 for e in range(eIdx): @@ -1686,9 +1685,9 @@ def createStomachMesh3d(region, fm, coordinates, stomachTermsAlong, allAnnotatio d1Path = [cd1Eso[0], [-d for d in normAxis]] ostiumLength = interp.computeCubicHermiteArcLength(xPath[0], d1Path[0], xPath[1], d1Path[1], rescaleDerivatives=True) - d1Path[1] = vector.setMagnitude(d1Path[1], ostiumLength*0.1) + d1Path[1] = set_magnitude(d1Path[1], ostiumLength*0.1) d2Path = [cd2Eso[0], d2Centre] - d3Path = [cd3Eso[0], vector.setMagnitude([-d for d in d1Centre], vector.magnitude(d2Centre))] + d3Path = [cd3Eso[0], set_magnitude([-d for d in d1Centre], magnitude(d2Centre))] d12Path = [cd2Eso[0], [0.0, 0.0, 0.0]] d13Path = [cd3Eso[0], [0.0, 0.0, 0.0]] @@ -1720,9 +1719,9 @@ def createStomachMesh3d(region, fm, coordinates, stomachTermsAlong, allAnnotatio d2AnnulusNorm = [] d2AnnulusOuter = [] for n1 in range(elementsCountAroundEso): - normD2 = vector.normalise(o1_d2[-1][n1]) + normD2 = normalize(o1_d2[-1][n1]) d2AnnulusNorm.append(normD2) - d2AnnulusOuter.append(vector.setMagnitude(o1_d2[-1][n1], sf)) + d2AnnulusOuter.append(set_magnitude(o1_d2[-1][n1], sf)) x = [o1_x[-1][n1][c] + sf * normD2[c] for c in range(3)] nearestPosition = trackSurfaceStomach.findNearestPosition(x, startPosition=o1_Positions[n1]) xAnnulusOuterPosition[n1] = nearestPosition @@ -1734,13 +1733,13 @@ def createStomachMesh3d(region, fm, coordinates, stomachTermsAlong, allAnnotatio v2 = xAnnulusOuter[(n + 1) % elementsCountAroundEso] d = [v2[c] - v1[c] for c in range(3)] arcLengthAround = interp.computeCubicHermiteArcLength(v1, d, v2, d, True) - d1 = [c * arcLengthAround for c in vector.normalise(d)] + d1 = [c * arcLengthAround for c in normalize(d)] d1AnnulusOuter.append(d1) d1AnnulusOuter = interp.smoothCubicHermiteDerivativesLoop(xAnnulusOuter, d2AnnulusOuter) d3Annulus = [] for n in range(elementsCountAroundEso): - d3 = vector.normalise(vector.crossproduct3(vector.normalise(d1AnnulusOuter[n]), d2AnnulusNorm[n])) + d3 = normalize(cross(normalize(d1AnnulusOuter[n]), d2AnnulusNorm[n])) d3Annulus.append(d3) d1AnnulusCurvatureOuter = findCurvatureAroundLoop(xAnnulusOuter, d1AnnulusOuter, d3Annulus) @@ -1817,7 +1816,7 @@ def createStomachMesh3d(region, fm, coordinates, stomachTermsAlong, allAnnotatio nxR, nd1R, nd2R, nd3R = \ trackSurfaceStomach.resampleHermiteCurvePointsSmooth( nx, nd1, nd2, nd3, proportions, derivativeMagnitudeStart= - vector.magnitude(d1EllipseAroundAll[sectionIdx[1]][elementsAroundQuarterDuod]))[0:-1] + magnitude(d1EllipseAroundAll[sectionIdx[1]][elementsAroundQuarterDuod]))[0:-1] # Replace the values in xEllipseAroundAll at quadrants for n in range(len(nxR)): @@ -1846,7 +1845,7 @@ def createStomachMesh3d(region, fm, coordinates, stomachTermsAlong, allAnnotatio nxR, nd1R, nd2R, nd3R = \ trackSurfaceStomach.resampleHermiteCurvePointsSmooth( nx, nd1, nd2, nd3, proportions, derivativeMagnitudeEnd= - vector.magnitude( + magnitude( d1EllipseAroundAll[sectionIdx[1]][elementsAroundQuarterDuod + elementsAroundHalfDuod]))[0:-1] for n in range(len(nxR)): @@ -1908,7 +1907,7 @@ def createStomachMesh3d(region, fm, coordinates, stomachTermsAlong, allAnnotatio d2Uniform = \ trackSurfaceStomach.resampleHermiteCurvePointsSmooth(nx, nd1, nd2, nd3, proportions)[1] startDerivative = d2Uniform[0] - startDerivativeMag = vector.magnitude(startDerivative) + startDerivativeMag = magnitude(startDerivative) # Sample from apex to annulus bPosition = xAnnulusOuterPosition[annulusIdxAtBodyStartIdxMinusOne[count]] @@ -2063,7 +2062,7 @@ def createStomachMesh3d(region, fm, coordinates, stomachTermsAlong, allAnnotatio nd2.append(d2SampledAroundAlong[n2][n1]) d2SmoothedAlongGC = interp.smoothCubicHermiteDerivativesLine(nx, nd2, fixAllDirections=True) d2SmoothedAlongGCB4Change = copy.deepcopy(d2SmoothedAlongGC) - d2SmoothedAlongGC[-1] = vector.setMagnitude(d2AnnulusOuter[0], vector.magnitude(d2SmoothedAlongGC[-1])) + d2SmoothedAlongGC[-1] = set_magnitude(d2AnnulusOuter[0], magnitude(d2SmoothedAlongGC[-1])) d2AnnulusNew[0] = d2SmoothedAlongGC[-1] nx = [] @@ -2073,8 +2072,8 @@ def createStomachMesh3d(region, fm, coordinates, stomachTermsAlong, allAnnotatio nd2.append(d2SampledAroundAlong[n2][n1]) d2SmoothedAlongLC = interp.smoothCubicHermiteDerivativesLine(nx, nd2, fixAllDirections=True) d2SmoothedAlongLCB4Change = copy.deepcopy(d2SmoothedAlongLC) - d2SmoothedAlongLC[0] = vector.setMagnitude(d2AnnulusOuter[elementsAroundHalfEso], - vector.magnitude(d2SmoothedAlongLC[0])) + d2SmoothedAlongLC[0] = set_magnitude(d2AnnulusOuter[elementsAroundHalfEso], + magnitude(d2SmoothedAlongLC[0])) d2AnnulusNew[elementsAroundHalfEso] = d2SmoothedAlongLC[0] d2Smoothed = d2SmoothedAlongGC + \ [[0.0, 1.0, 0.0] for n in range(2 * (elementsAroundQuarterEso - 2) + 1)] + \ @@ -2092,36 +2091,36 @@ def createStomachMesh3d(region, fm, coordinates, stomachTermsAlong, allAnnotatio if n1 == elementsAroundHalfDuod - 1: d2Smoothed[annulusFundusOpenRingIdx - 1] = \ - vector.setMagnitude(d2AnnulusOuter[1], vector.magnitude(nd2[annulusFundusOpenRingIdx - 1])) + set_magnitude(d2AnnulusOuter[1], magnitude(nd2[annulusFundusOpenRingIdx - 1])) d2AnnulusNew[1] = d2Smoothed[annulusFundusOpenRingIdx - 1] for m in range(2 * (elementsAroundQuarterEso - 2) + 1): annulusIdx = m + 2 d2Smoothed[annulusFundusOpenRingIdx + m] = \ - vector.setMagnitude(d2AnnulusOuter[annulusIdx], - vector.magnitude(d1SampledAroundAlong[annulusFundusOpenRingIdx + m][n1])) + set_magnitude(d2AnnulusOuter[annulusIdx], + magnitude(d1SampledAroundAlong[annulusFundusOpenRingIdx + m][n1])) d2AnnulusNew[annulusIdx] = d2Smoothed[annulusFundusOpenRingIdx + m] d2Smoothed[annulusBodyOpenRingIdx + 1] = \ - vector.setMagnitude(d2AnnulusOuter[elementsAroundHalfEso - 1], - vector.magnitude(nd2[annulusBodyOpenRingIdx + 1])) + set_magnitude(d2AnnulusOuter[elementsAroundHalfEso - 1], + magnitude(nd2[annulusBodyOpenRingIdx + 1])) d2AnnulusNew[elementsAroundHalfEso - 1] = d2Smoothed[annulusBodyOpenRingIdx + 1] if n1 == elementsAroundHalfDuod + 1: d2Smoothed[annulusFundusOpenRingIdx - 1] = \ - vector.setMagnitude(d2AnnulusOuter[-1], vector.magnitude(nd2[annulusFundusOpenRingIdx - 1])) + set_magnitude(d2AnnulusOuter[-1], magnitude(nd2[annulusFundusOpenRingIdx - 1])) d2AnnulusNew[-1] = d2Smoothed[annulusFundusOpenRingIdx - 1] for m in range(2 * (elementsAroundQuarterEso - 2) + 1): annulusIdx = -(m + 2) d2Smoothed[annulusFundusOpenRingIdx + m] = \ - vector.setMagnitude(d2AnnulusOuter[annulusIdx], - vector.magnitude(d1SampledAroundAlong[annulusFundusOpenRingIdx + m][n1])) + set_magnitude(d2AnnulusOuter[annulusIdx], + magnitude(d1SampledAroundAlong[annulusFundusOpenRingIdx + m][n1])) d2AnnulusNew[annulusIdx] = d2Smoothed[annulusFundusOpenRingIdx + m] d2Smoothed[annulusBodyOpenRingIdx + 1] = \ - vector.setMagnitude(d2AnnulusOuter[elementsAroundHalfEso + 1], - vector.magnitude(nd2[annulusBodyOpenRingIdx + 1])) + set_magnitude(d2AnnulusOuter[elementsAroundHalfEso + 1], + magnitude(nd2[annulusBodyOpenRingIdx + 1])) d2AnnulusNew[elementsAroundHalfEso + 1] = d2Smoothed[annulusBodyOpenRingIdx + 1] for n2 in range(len(d2Smoothed)): @@ -2153,8 +2152,8 @@ def createStomachMesh3d(region, fm, coordinates, stomachTermsAlong, allAnnotatio # calculate d3 for n2 in range(len(xSampledAroundAlong)): for n1 in range(len(xSampledAroundAlong[n2])): - d3SampledAroundAlong[n2][n1] = vector.normalise(vector.crossproduct3( - vector.normalise(d1SampledAroundAlong[n2][n1]), vector.normalise(d2SampledAroundAlong[n2][n1]))) + d3SampledAroundAlong[n2][n1] = normalize(cross( + normalize(d1SampledAroundAlong[n2][n1]), normalize(d2SampledAroundAlong[n2][n1]))) # for n2 in range(len(xSampledAroundAlong)): # for n1 in range(len(xSampledAroundAlong[n2])): @@ -2202,7 +2201,7 @@ def createStomachMesh3d(region, fm, coordinates, stomachTermsAlong, allAnnotatio # Calculate curvature along d2AnnulusCurvature = [] for n in range(elementsCountAroundEso): - d2AnnulusCurvature.append(interp.getCubicHermiteCurvature(o1_x[-1][n], vector.setMagnitude(o1_d2[-1][n], sf), + d2AnnulusCurvature.append(interp.getCubicHermiteCurvature(o1_x[-1][n], set_magnitude(o1_d2[-1][n], sf), xAnnulusOuter[n], d2AnnulusNew[n], d3Annulus[n], 1.0)) d2CurvatureAroundAlong = [[[] for n1 in range(len(xSampledAroundAlong[n2]))] @@ -2288,8 +2287,8 @@ def createStomachMesh3d(region, fm, coordinates, stomachTermsAlong, allAnnotatio xSampledAroundAlong[1][elementsAroundQuarterDuod], d2SampledAroundAlong[1][elementsAroundQuarterDuod], rescaleDerivatives=True) - d1SampledAroundAlong[0][0] = vector.setMagnitude(d1SampledAroundAlong[0][0], arcLength) - d2SampledAroundAlong[0][0] = vector.setMagnitude(d2SampledAroundAlong[0][0], arcLength) + d1SampledAroundAlong[0][0] = set_magnitude(d1SampledAroundAlong[0][0], arcLength) + d2SampledAroundAlong[0][0] = set_magnitude(d2SampledAroundAlong[0][0], arcLength) # Replace d1Curvature with d2Curvature d1CurvatureAroundAlong[0][0] = d2CurvatureAroundAlong[0][0] @@ -2345,7 +2344,7 @@ def createStomachMesh3d(region, fm, coordinates, stomachTermsAlong, allAnnotatio for n1 in range(len(xSampledAroundAlong[n2])): # Coordinates - norm = vector.normalise(d3SampledAroundAlong[n2][n1]) + norm = normalize(d3SampledAroundAlong[n2][n1]) xOut = xSampledAroundAlong[n2][n1] xIn = [xOut[i] - norm[i] * wallThickness for i in range(3)] dWall = [wallThickness * c for c in norm] diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_uterus1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_uterus1.py index 03d4b730..5cb86196 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_uterus1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_uterus1.py @@ -21,7 +21,6 @@ from scaffoldmaker.scaffoldpackage import ScaffoldPackage from scaffoldmaker.utils import interpolation as interp from scaffoldmaker.utils import tubemesh -from scaffoldmaker.utils import vector from scaffoldmaker.utils.eftfactory_bicubichermitelinear import eftfactory_bicubichermitelinear from scaffoldmaker.utils.geometry import createEllipsePoints from scaffoldmaker.utils.interpolation import ( @@ -549,7 +548,7 @@ def findDerivativeBetweenPoints(v1, v2): """ d = [v2[c] - v1[c] for c in range(3)] arcLengthAround = interp.computeCubicHermiteArcLength(v1, d, v2, d, True) - d = [c * arcLengthAround for c in vector.normalise(d)] + d = [c * arcLengthAround for c in normalize(d)] return d @@ -1143,8 +1142,8 @@ def getTubeNodes(cx_group, elementsCountAround, elementsCountAlongTube, elements for n2 in range(elementsCountAlongTube + 1): d3Around = [] for n1 in range(elementsCountAround): - d3Around.append(vector.normalise( - vector.crossproduct3(vector.normalise(d1SampledTube[n2][n1]), vector.normalise(d2SampledTube[n2][n1])))) + d3Around.append(normalize( + cross(normalize(d1SampledTube[n2][n1]), normalize(d2SampledTube[n2][n1])))) d3Tube.append(d3Around) xInner = [] @@ -1802,25 +1801,25 @@ def getDoubleTubeNodes(cx_tube_group, elementsCountAlong, elementsCountAround, e d3lList = [] for n in range(len(cx_tube_group[0])): x = cx_tube_group[0][n] - v2r = vector.normalise(cx_tube_group[2][n]) - v3r = vector.normalise(cx_tube_group[4][n]) - d2Mag = vector.magnitude(cx_tube_group[2][n]) + v2r = normalize(cx_tube_group[2][n]) + v3r = normalize(cx_tube_group[4][n]) + d2Mag = magnitude(cx_tube_group[2][n]) r2 = (d2Mag - distance - wallThickness) / 2 - # d3Mag = vector.magnitude(cx_tube_group[4][n]) + # d3Mag = magnitude(cx_tube_group[4][n]) # r3 = (d3Mag - wallThickness) - v_trans = vector.setMagnitude(v2r, distance + r2) + v_trans = set_magnitude(v2r, distance + r2) x_right = [x[c] + v_trans[c] for c in range(3)] xrList.append(x_right) - d2rList.append(vector.setMagnitude(v2r, r2)) - d3rList.append(vector.setMagnitude(v3r, r2)) + d2rList.append(set_magnitude(v2r, r2)) + d3rList.append(set_magnitude(v3r, r2)) x_left = [x[c] - v_trans[c] for c in range(3)] v2l = [-v2r[c] for c in range(3)] v3l = [-v3r[c] for c in range(3)] xlList.append(x_left) - d2lList.append(vector.setMagnitude(v2l, r2)) - d3lList.append(vector.setMagnitude(v3l, r2)) + d2lList.append(set_magnitude(v2l, r2)) + d3lList.append(set_magnitude(v3l, r2)) cx_tube_group_right = [xrList, cx_tube_group[1], d2rList, [], d3rList] cx_tube_group_left = [xlList, cx_tube_group[1], d2lList, [], d3lList] diff --git a/src/scaffoldmaker/scaffoldpackage.py b/src/scaffoldmaker/scaffoldpackage.py index 9aedfef4..4ec401ce 100644 --- a/src/scaffoldmaker/scaffoldpackage.py +++ b/src/scaffoldmaker/scaffoldpackage.py @@ -6,7 +6,7 @@ import copy import math -from cmlibs.maths.vectorops import euler_to_rotation_matrix +from cmlibs.maths.vectorops import euler_to_rotation_matrix, magnitude from cmlibs.utils.zinc.field import createFieldEulerAnglesRotationMatrix from cmlibs.utils.zinc.finiteelement import get_highest_dimension_mesh, get_maximum_node_identifier from cmlibs.utils.zinc.general import ChangeManager @@ -14,7 +14,6 @@ from scaffoldmaker.annotation.annotationgroup import AnnotationGroup, findAnnotationGroupByName, \ getAnnotationMarkerLocationField # , getAnnotationMarkerNameField from scaffoldmaker.meshtypes.scaffold_base import Scaffold_base -from scaffoldmaker.utils import vector class ScaffoldPackage: @@ -348,7 +347,7 @@ def deleteElementsInRanges(self, region, deleteElementRanges): diff = [abs(evaluatedMaterialCoordinates[c] - materialCoordinates[c]) for c in range(3)] # threshold designed for material coordinates of nominally unit scale - if vector.magnitude(diff) < 1e-03: + if magnitude(diff) < 1e-03: destroyNodes.removeNode(annotationGroup.getMarkerNode()) removeMarkerGroup = False diff --git a/src/scaffoldmaker/utils/annulusmesh.py b/src/scaffoldmaker/utils/annulusmesh.py index 3c4e5820..dc92f5ea 100644 --- a/src/scaffoldmaker/utils/annulusmesh.py +++ b/src/scaffoldmaker/utils/annulusmesh.py @@ -6,11 +6,11 @@ import copy from collections.abc import Sequence +from cmlibs.maths.vectorops import magnitude, set_magnitude, normalize, cross, dot from cmlibs.utils.zinc.field import findOrCreateFieldCoordinates from cmlibs.zinc.element import Element from cmlibs.zinc.node import Node from scaffoldmaker.utils import interpolation as interp -from scaffoldmaker.utils import vector from scaffoldmaker.utils.eft_utils import remapEftNodeValueLabel, setEftScaleFactorIds from scaffoldmaker.utils.eftfactory_bicubichermitelinear import eftfactory_bicubichermitelinear from scaffoldmaker.utils.eftfactory_tricubichermite import eftfactory_tricubichermite @@ -206,11 +206,11 @@ def createAnnulusMesh3d(nodes, mesh, nextNodeIdentifier, nextElementIdentifier, # Find total wall thickness thicknessProportions = [] thicknesses = [] - thicknesses.append([vector.magnitude([(startPointsx[nodesCountWall - 1][n1][c] - startPointsx[0][n1][c]) + thicknesses.append([magnitude([(startPointsx[nodesCountWall - 1][n1][c] - startPointsx[0][n1][c]) for c in range(3)]) for n1 in range(nodesCountAround)]) for n2 in range(1, elementsCountRadial): thicknesses.append([None] * nodesCountAround) - thicknesses.append([vector.magnitude([(endPointsx[nodesCountWall - 1][n1][c] - endPointsx[0][n1][c]) + thicknesses.append([magnitude([(endPointsx[nodesCountWall - 1][n1][c] - endPointsx[0][n1][c]) for c in range(3)]) for n1 in range(nodesCountAround)]) for n3 in range(nodesCountWall): @@ -221,10 +221,10 @@ def createAnnulusMesh3d(nodes, mesh, nextNodeIdentifier, nextElementIdentifier, endPointsd3[n3] if (endPointsd3 is not None) else None] startThicknessList = \ - [vector.magnitude([(startPointsx[n3][n1][c] - startPointsx[n3 - (1 if n3 > 0 else 0)][n1][c]) + [magnitude([(startPointsx[n3][n1][c] - startPointsx[n3 - (1 if n3 > 0 else 0)][n1][c]) for c in range(3)]) for n1 in range(len(startPointsx[n3]))] endThicknessList = \ - [vector.magnitude([(endPointsx[n3][n1][c] - endPointsx[n3 - (1 if n3 > 0 else 0)][n1][c]) + [magnitude([(endPointsx[n3][n1][c] - endPointsx[n3 - (1 if n3 > 0 else 0)][n1][c]) for c in range(3)]) for n1 in range(len(endPointsx[n3]))] thicknessList = [startThicknessList, endThicknessList] # thickness of each layer @@ -251,14 +251,14 @@ def createAnnulusMesh3d(nodes, mesh, nextNodeIdentifier, nextElementIdentifier, bd1, bd2 = getMappedD1D2([endPointsd1[n3][n1], endPointsd2[n3][n1]] + ([endPointsd3[n3][n1]] if endPointsd3 else []), endDerivativesMap[n3][n1] if endDerivativesMap else None) - bd2 = vector.setMagnitude(bd2, vector.magnitude(ad2)) + bd2 = set_magnitude(bd2, magnitude(ad2)) scaling = interp.computeCubicHermiteDerivativeScaling(ax, ad2, bx, bd2) if fixMinimumStart: - mag = vector.magnitude(ad2) * scaling + mag = magnitude(ad2) * scaling if (fixStartMagnitude is None) or (mag < fixStartMagnitude): fixStartMagnitude = mag else: # fixMinimumEnd: - mag = vector.magnitude(bd2) * scaling + mag = magnitude(bd2) * scaling if (fixEndMagnitude is None) or (mag < fixEndMagnitude): fixEndMagnitude = mag @@ -290,17 +290,17 @@ def createAnnulusMesh3d(nodes, mesh, nextNodeIdentifier, nextElementIdentifier, # sample between start and end points and derivatives # scaling end derivatives to arc length gives even curvature along the curve - aMag = vector.magnitude(ad2) - bMag = vector.magnitude(bd2) + aMag = magnitude(ad2) + bMag = magnitude(bd2) ad2mag = 0.5 * ((1.0 + sampleBlend) * aMag + (1.0 - sampleBlend) * bMag) - ad2Scaled = vector.setMagnitude(ad2, ad2mag) if (aMag > 0.0) else [0.0, 0.0, 0.0] + ad2Scaled = set_magnitude(ad2, ad2mag) if (aMag > 0.0) else [0.0, 0.0, 0.0] bd2mag = 0.5 * ((1.0 + sampleBlend) * bMag + (1.0 - sampleBlend) * aMag) - bd2Scaled = vector.setMagnitude(bd2, bd2mag) if (bMag > 0.0) else [0.0, 0.0, 0.0] + bd2Scaled = set_magnitude(bd2, bd2mag) if (bMag > 0.0) else [0.0, 0.0, 0.0] scaling = interp.computeCubicHermiteDerivativeScaling(ax, ad2Scaled, bx, bd2Scaled) ad2Scaled = [d * scaling for d in ad2Scaled] bd2Scaled = [d * scaling for d in bd2Scaled] - derivativeMagnitudeStart = None if rescaleStartDerivatives else vector.magnitude(ad2) - derivativeMagnitudeEnd = None if rescaleEndDerivatives else vector.magnitude(bd2) + derivativeMagnitudeStart = None if rescaleStartDerivatives else magnitude(ad2) + derivativeMagnitudeEnd = None if rescaleEndDerivatives else magnitude(bd2) if fixStartMagnitude: derivativeMagnitudeStart = fixStartMagnitude / elementsCountRadial elif fixEndMagnitude: @@ -345,10 +345,10 @@ def createAnnulusMesh3d(nodes, mesh, nextNodeIdentifier, nextElementIdentifier, # set scalefactors if rescaling, make same on inside for now if rescaleStartDerivatives: - scaleFactor = vector.magnitude(md2[0]) / vector.magnitude(ad2) + scaleFactor = magnitude(md2[0]) / magnitude(ad2) scaleFactorMapStart[n3].append(scaleFactor) if rescaleEndDerivatives: - scaleFactor = vector.magnitude(md2[-1]) / vector.magnitude(bd2) + scaleFactor = magnitude(md2[-1]) / magnitude(bd2) scaleFactorMapEnd[n3].append(scaleFactor) for n2 in range(1, elementsCountRadial): @@ -378,7 +378,7 @@ def createAnnulusMesh3d(nodes, mesh, nextNodeIdentifier, nextElementIdentifier, for n3 in range(0, nodesCountWall - 1): for n1 in range(nodesCountAround): xi3 = 1 - xi3List[n3][n2][n1] - normal = vector.normalise(vector.crossproduct3(pd1[-1][n2][n1], pd2[-1][n2][n1])) + normal = normalize(cross(pd1[-1][n2][n1], pd2[-1][n2][n1])) thickness = thicknesses[n2][n1] * xi3 d3 = [d * thickness for d in normal] px[n3][n2][n1] = [(px[-1][n2][n1][c] - d3[c]) for c in range(3)] @@ -403,7 +403,7 @@ def createAnnulusMesh3d(nodes, mesh, nextNodeIdentifier, nextElementIdentifier, factor = 1.0 + curvature * thickness pd2[n3][n2][n1] = [factor * d for d in pd2[-1][n2][n1]] d2Scaled = [factor * d for d in pd2[-1][n2][n1]] - if vector.dotproduct(vector.normalise(pd2[-1][n2][n1]), vector.normalise(d2Scaled)) == -1: + if dot(normalize(pd2[-1][n2][n1]), normalize(d2Scaled)) == -1: pd2[n3][n2][n1] = [-factor * d for d in pd2[-1][n2][n1]] if not midLinearXi3: pd3[n3][n2][n1] = pd3[-1][n2][n1] = \ @@ -435,10 +435,10 @@ def createAnnulusMesh3d(nodes, mesh, nextNodeIdentifier, nextElementIdentifier, magnitudeScalingMode=interp.DerivativeScalingMode. HARMONIC_MEAN) if rescaleStartDerivatives: - scaleFactor = vector.magnitude(sd2[0]) / vector.magnitude(md2[0]) + scaleFactor = magnitude(sd2[0]) / magnitude(md2[0]) scaleFactorMapStart[n3].append(scaleFactor) if rescaleEndDerivatives: - scaleFactor = vector.magnitude(sd2[-1]) / vector.magnitude(md2[-1]) + scaleFactor = magnitude(sd2[-1]) / magnitude(md2[-1]) scaleFactorMapEnd[n3].append(scaleFactor) for n2 in range(1, elementsCountRadial): diff --git a/src/scaffoldmaker/utils/cylindermesh.py b/src/scaffoldmaker/utils/cylindermesh.py index 3a4782a6..bc86b514 100644 --- a/src/scaffoldmaker/utils/cylindermesh.py +++ b/src/scaffoldmaker/utils/cylindermesh.py @@ -7,10 +7,11 @@ import math from enum import Enum +from cmlibs.maths.vectorops import set_magnitude, cross, magnitude, normalize, dot from cmlibs.utils.zinc.finiteelement import getMaximumNodeIdentifier, getMaximumElementIdentifier from cmlibs.zinc.field import Field from cmlibs.zinc.node import Node -from scaffoldmaker.utils import vector, geometry +from scaffoldmaker.utils import geometry from scaffoldmaker.utils.interpolation import sampleCubicHermiteCurves, interpolateSampleCubicHermite, \ smoothCubicHermiteDerivativesLine from scaffoldmaker.utils.mirror import Mirror @@ -63,13 +64,13 @@ def __init__(self, elementsCountAcrossMajor, elementsCountAcrossMinor, elementsC self._majorAxis = majorAxis self._minorRadius = minorRadius if alongAxis: - self._minorAxis = vector.setMagnitude(vector.crossproduct3(alongAxis, majorAxis), minorRadius) + self._minorAxis = set_magnitude(cross(alongAxis, majorAxis), minorRadius) self._elementsCountAcrossMinor = elementsCountAcrossMinor self._elementsCountAcrossMajor = elementsCountAcrossMajor self._elementsCountAcrossShell = elementsCountAcrossShell self._elementsCountAcrossTransition = elementsCountAcrossTransition self._shellProportion = shellProportion - self._majorRadius = vector.magnitude(majorAxis) + self._majorRadius = magnitude(majorAxis) self.px = None self.pd1 = None self.pd2 = None @@ -123,9 +124,9 @@ def __init__(self, region, centralPath, elementsCount): sd3, sd13 = interpolateSampleCubicHermite(cd3, cd13, se, sxi, ssf) self.centres = sx - self.majorRadii = [vector.magnitude(a) for a in sd2] + self.majorRadii = [magnitude(a) for a in sd2] self.majorAxis = sd2 - self.minorRadii = [vector.magnitude(a) for a in sd3] + self.minorRadii = [magnitude(a) for a in sd3] self.minorAxis = sd3 self.alongAxis = sd1 @@ -199,7 +200,7 @@ def __init__(self, fieldModule, coordinates, elementsCountAlong, base=None, end= base._shellProportion, self._centres[0], None, self._majorAxis[0], self._minorRadii[0]) else: - self._length = vector.magnitude(base._alongAxis) + self._length = magnitude(base._alongAxis) arcLengthAlong = self._length / elementsCountAlong self.calculateEllipseParams(arcLengthAlong, cylinderCentralPath=self._cylinderCentralPath) @@ -256,13 +257,13 @@ def createCylinderMesh3d(self, fieldModule, coordinates): # The other ellipses for a straight cylinder. if self._cylinderType == CylinderType.CYLINDER_STRAIGHT: - arcLengthAlong = vector.magnitude(self._base._alongAxis) / self._elementsCountAlong + arcLengthAlong = magnitude(self._base._alongAxis) / self._elementsCountAlong for n2 in range(self._elementsCountUp + 1): for n3 in range(self._elementsCountAlong + 1): for n1 in range(self._elementsCountAcrossMinor + 1): if self._shield.px[0][n2][n1]: temx = [self._shield.px[0][n2][n1][c] + n3 * arcLengthAlong * - vector.normalise(self._base._alongAxis)[c] for c in range(3)] + normalize(self._base._alongAxis)[c] for c in range(3)] self._shield.px[n3][n2][n1] = temx self._shield.pd1[n3][n2][n1] = self._shield.pd1[0][n2][n1] self._shield.pd2[n3][n2][n1] = self._shield.pd2[0][n2][n1] @@ -276,7 +277,7 @@ def createCylinderMesh3d(self, fieldModule, coordinates): self._elementsCountAcrossShell, self._elementsCountAcrossTransition, self._shellProportion, self._centres[-1], self._shield.pd2[-1][0][1], - vector.setMagnitude(self._base._majorAxis, self._majorRadii[-1]), + set_magnitude(self._base._majorAxis, self._majorRadii[-1]), self._minorRadii[-1]) self.setEndsNodes() @@ -437,8 +438,8 @@ def __init__(self, centre, majorAxis, minorAxis, self.centre = centre self.majorAxis = majorAxis self.minorAxis = minorAxis - self.majorRadius = vector.magnitude(majorAxis) - self.minorRadius = vector.magnitude(minorAxis) + self.majorRadius = magnitude(majorAxis) + self.minorRadius = magnitude(minorAxis) self.coreMajorRadius = coreMajorRadius self.coreMinorRadius = coreMinorRadius elementsCountRim = elementsCountAcrossShell + elementsCountAcrossTransition - 1 @@ -495,7 +496,7 @@ def generateBase1DMesh(self, rx): ratio = rx/self.elementsCountAcrossShell if self.elementsCountAcrossShell > 0 else 0 majorAxis = [d*(1 - ratio*(1-self.coreMajorRadius/self.majorRadius)) for d in self.majorAxis] minorAxis = [d*(1 - ratio*(1-self.coreMinorRadius/self.minorRadius)) for d in self.minorAxis] - majorRadius = vector.magnitude(majorAxis) + majorRadius = magnitude(majorAxis) nx, nd1 = createEllipsePerimeter( self.centre, majorAxis, minorAxis, self.elementsCountAround, majorRadius) @@ -506,7 +507,7 @@ def generateBase1DMesh(self, rx): tbx.append(nx[n]) tbd1.append(nd1[n]) tbd2.append(nte) - tbd3.append(vector.normalise(vector.crossproduct3(tbd1[n], nte))) + tbd3.append(normalize(cross(tbd1[n], nte))) for n in range(self.elementsCountAround + 1): n1, n2 = self.__shield.convertRimIndex(n, rx) @@ -562,17 +563,17 @@ def createMirrorCurve(self): [c]) for c in range(3)] rcx.append(tmdx) rcx.append(tmux) - rcd3 = [vector.setMagnitude(tmdd3, -1), vector.setMagnitude(tmdd3, -1)] + rcd3 = [set_magnitude(tmdd3, -1), set_magnitude(tmdd3, -1)] rscx, rscd1 = sampleCubicHermiteCurves(rcx, rcd3, self.elementsCountUp - n2a, arcLengthDerivatives=True)[0:2] # get d2, d3 rscd2 = [] rscd3 = [] for n in range(len(rscx)): - d3 = vector.normalise( + d3 = normalize( [btx[self.elementsCountUp][self.elementsCountAcrossMinor][c] - btx[self.elementsCountUp][0] [c] for c in range(3)]) - d2 = vector.normalise(vector.crossproduct3(d3, rscd1[n])) + d2 = normalize(cross(d3, rscd1[n])) rscd2.append(d2) rscd3.append(d3) @@ -601,7 +602,7 @@ def createRegularRowCurves(self, rscx, rscd1, rscd3): for n2 in range(n2d, n2m + 1): txm, td3m, pe, pxi, psf = sampleCubicHermiteCurves( [btx[n2][n1a], rscx[n2 - n2a], btx[n2][n1z]], - [vector.setMagnitude(btd3[n2][n1a], -1.0), rscd3[n2 - n2a], btd3[n2][n1z]], + [set_magnitude(btd3[n2][n1a], -1.0), rscd3[n2 - n2a], btd3[n2][n1z]], self.elementsCountAcrossMinor-2*self.elementsCountAcrossShell, arcLengthDerivatives=True) td1m = interpolateSampleCubicHermite([[-btd1[n2][n1a][c] for c in range(3)], rscd1[n2 - n2a], btd1[n2][n1z]], [[0.0, 0.0, 0.0]] * 3, pe, pxi, psf)[0] @@ -628,9 +629,9 @@ def createRegularRowCurves(self, rscx, rscd1, rscd3): btx[n2][n1] = tx[n1] if n2 == n2m: if n1 <= n1b: - btd1[n2][n1] = vector.setMagnitude(btd1[n2m][-1], -1.0) + btd1[n2][n1] = set_magnitude(btd1[n2m][-1], -1.0) else: - btd1[n2][n1] = vector.setMagnitude(btd1[n2m][-1], 1.0) + btd1[n2][n1] = set_magnitude(btd1[n2m][-1], 1.0) else: if n1 <= n1b: btd1[n2][n1] = [-d for d in td1m[n1 - n1a]] @@ -756,7 +757,7 @@ def generateNodesForUpperHalf(self): It keeps the d1 direction. It uses mirrorPlane: plane ax+by+cz=d in form of [a,b,c,d] """ - mirrorPlane = [-d for d in self.majorAxis] + [-vector.dotproduct(self.majorAxis, self.centre)] + mirrorPlane = [-d for d in self.majorAxis] + [-dot(self.majorAxis, self.centre)] mirror = Mirror(mirrorPlane) for n2 in range(self.elementsCountUp): for n1 in range(self.elementsCountAcrossMinor + 1): @@ -798,10 +799,10 @@ def createEllipsePerimeter(centre, majorAxis, minorAxis, elementsCountAround, he """ nx = [] nd1 = [] - magMajorAxis = vector.magnitude(majorAxis) - magMinorAxis = vector.magnitude(minorAxis) - unitMajorAxis = vector.normalise(majorAxis) - unitMinorAxis = vector.normalise(minorAxis) + magMajorAxis = magnitude(majorAxis) + magMinorAxis = magnitude(minorAxis) + unitMajorAxis = normalize(majorAxis) + unitMinorAxis = normalize(minorAxis) useHeight = min(max(0.0, height), 2.0 * magMajorAxis) totalRadians = math.acos((magMajorAxis - useHeight)/magMajorAxis) @@ -818,7 +819,7 @@ def createEllipsePerimeter(centre, majorAxis, minorAxis, elementsCountAround, he nx.append( [(centre[c] + cosRadians * majorAxis[c] + sinRadians * minorAxis[c]) for c in range(3)]) - ndab = vector.setMagnitude([-sinRadians * magMajorAxis, cosRadians * magMinorAxis], elementArcLength) + ndab = set_magnitude([-sinRadians * magMajorAxis, cosRadians * magMinorAxis], elementArcLength) nd1.append( [(ndab[0] * unitMajorAxis[c] + ndab[1] * unitMinorAxis[c]) for c in range(3)]) radians = geometry.updateEllipseAngleByArcLength(magMajorAxis, magMinorAxis, radians, elementArcLength, @@ -833,7 +834,7 @@ def normalToEllipse(v1, v2): :param v2: vector 2. :return: """ - nte = vector.normalise(vector.crossproduct3(v1, v2)) + nte = normalize(cross(v1, v2)) return nte @@ -850,7 +851,7 @@ def computeNextRadius(radius, axis, ratio, progression): radius = radius * ratio elif progression == ConeBaseProgression.ARITHMETIC_PROGRESSION: radius += ratio - axis = vector.setMagnitude(axis, radius) + axis = set_magnitude(axis, radius) return radius, axis @@ -862,5 +863,5 @@ def computeNextCentre(centre, arcLength, axis): :param centre: the start centre. :return: next centre coordinates. """ - centre = [centre[c]+arcLength * vector.normalise(axis)[c] for c in range(3)] + centre = [centre[c]+arcLength * normalize(axis)[c] for c in range(3)] return centre diff --git a/src/scaffoldmaker/utils/eftfactory_tricubichermite.py b/src/scaffoldmaker/utils/eftfactory_tricubichermite.py index 0adf3ef5..c3a01b57 100644 --- a/src/scaffoldmaker/utils/eftfactory_tricubichermite.py +++ b/src/scaffoldmaker/utils/eftfactory_tricubichermite.py @@ -3,12 +3,12 @@ ''' import math +from cmlibs.maths.vectorops import normalize, cross from cmlibs.utils.zinc.field import findOrCreateFieldCoordinates from cmlibs.utils.zinc.finiteelement import getElementNodeIdentifiers from cmlibs.zinc.element import Element, Elementbasis, Elementfieldtemplate from cmlibs.zinc.field import Field from cmlibs.zinc.node import Node -from scaffoldmaker.utils import vector from scaffoldmaker.utils.eft_utils import mapEftFunction1Node1Term, remapEftLocalNodes, remapEftNodeValueLabel, scaleEftNodeValueLabels, setEftScaleFactorIds @@ -1497,19 +1497,19 @@ def replaceElementWithInlet4(self, origElement, startElementId, nodetemplate, st result, fc = coordinates.evaluateReal(cache, 3) resulta, a = coordinates.evaluateDerivative(diff1, cache, 3) resultb, b = coordinates.evaluateDerivative(diff2, cache, 3) - n = vector.normalise(vector.crossproduct3(a, b)) - t = vector.normalise(a if (inclineAxis == 1) else b) + n = normalize(cross(a, b)) + t = normalize(a if (inclineAxis == 1) else b) n = [ (math.cos(inclineRadians)*n[c] + math.sin(inclineRadians)*t[c]) for c in range(3) ] #print(resulta, 'a =', a, ',', resultb, ' b=', b, ' fc=', fc, ' n=',n) ic = [ (fc[i] + tubeLength*n[i]) for i in range(3) ] if inclineAxis == 2: - na = vector.normalise(a) - nb = vector.normalise(vector.crossproduct3(n, a)) + na = normalize(a) + nb = normalize(cross(n, a)) else: - nb = vector.normalise(b) - na = vector.normalise(vector.crossproduct3(b, n)) - a = vector.normalise([ -(na[i] + nb[i]) for i in range(3) ]) - b = vector.normalise(vector.crossproduct3(a, n)) + nb = normalize(b) + na = normalize(cross(b, n)) + a = normalize([ -(na[i] + nb[i]) for i in range(3) ]) + b = normalize(cross(a, n)) zero = [ 0.0, 0.0, 0.0 ] nodeIdentifier = startNodeId @@ -1615,8 +1615,8 @@ def replaceTwoElementWithInlet6(self, origElement1, origElement2, startElementId fm.beginChange() cache = fm.createFieldcache() coordinates = findOrCreateFieldCoordinates(fm) - a = vector.normalise(inletSide) - b = vector.normalise(vector.crossproduct3(inletAxis, inletSide)) + a = normalize(inletSide) + b = normalize(cross(inletAxis, inletSide)) zero = [ 0.0, 0.0, 0.0 ] nodeIdentifier = startNodeId diff --git a/src/scaffoldmaker/utils/geometry.py b/src/scaffoldmaker/utils/geometry.py index bdad3abf..5e80ba75 100644 --- a/src/scaffoldmaker/utils/geometry.py +++ b/src/scaffoldmaker/utils/geometry.py @@ -7,8 +7,7 @@ import copy import math -from cmlibs.maths.vectorops import magnitude, mult -from scaffoldmaker.utils import vector +from cmlibs.maths.vectorops import magnitude, mult, normalize, cross, set_magnitude from scaffoldmaker.utils.tracksurface import calculate_surface_delta_xi @@ -216,11 +215,11 @@ def createEllipsoidPoints(centre, poleAxis, sideAxis, elementsCountAround, eleme nx = [] nd1 = [] nd2 = [] - magPoleAxis = vector.magnitude(poleAxis) - magSideAxis = vector.magnitude(sideAxis) - unitPoleAxis = vector.normalise(poleAxis) - unitSideAxis1 = vector.normalise(sideAxis) - unitSideAxis2 = vector.normalise(vector.crossproduct3(sideAxis, poleAxis)) + magPoleAxis = magnitude(poleAxis) + magSideAxis = magnitude(sideAxis) + unitPoleAxis = normalize(poleAxis) + unitSideAxis1 = normalize(sideAxis) + unitSideAxis2 = normalize(cross(sideAxis, poleAxis)) useHeight = min(max(0.0, height), 2.0*magPoleAxis) totalRadiansUp = getEllipseRadiansToX(magPoleAxis, 0.0, magPoleAxis - useHeight, initialTheta = 0.5*math.pi*useHeight/magPoleAxis) radiansUp = 0.0 @@ -231,7 +230,7 @@ def createEllipsoidPoints(centre, poleAxis, sideAxis, elementsCountAround, eleme cosRadiansUp = math.cos(radiansUp) sinRadiansUp = math.sin(radiansUp) radius = sinRadiansUp*magSideAxis - d2r, d2z = vector.setMagnitude([ cosRadiansUp*magSideAxis, sinRadiansUp*magPoleAxis ], elementLengthUp) + d2r, d2z = set_magnitude([ cosRadiansUp*magSideAxis, sinRadiansUp*magPoleAxis ], elementLengthUp) cx = [ (centre[c] + cosRadiansUp*poleAxis[c]) for c in range(3) ] elementLengthAround = radius*radiansPerElementAround radiansAround = 0.0 @@ -442,9 +441,9 @@ def getCircleProjectionAxes(ax, ad1, ad2, ad3, length, angle1radians, angle2radi w2 = br*math.sin(angleAroundRadians) # f2/fh w3 = arcRadius*math.sin(arcAngleRadians) bx = [ (ax[c] + w1*ad1[c] + w2*ad2[c] + w3*ad3[c]) for c in range(3) ] - bd3 = vector.normalise([ (f1*ad1[c] + f2*ad2[c] + f3*ad3[c]) for c in range(3) ]) - bd1 = vector.normalise(vector.crossproduct3(ad2, bd3)) - bd2 = vector.crossproduct3(bd3, bd1) + bd3 = normalize([ (f1*ad1[c] + f2*ad2[c] + f3*ad3[c]) for c in range(3) ]) + bd1 = normalize(cross(ad2, bd3)) + bd2 = cross(bd3, bd1) if angle3radians: cosAngle3 = math.cos(angle3radians) sinAngle3 = math.sin(angle3radians) @@ -470,8 +469,8 @@ def getSurfaceProjectionAxes(ax, ad1, ad2, ad3, angle1radians, angle2radians, le f3 = cosAngle1*cosAngle2 bd3 = [ (f1*ad1[c] + f2*ad2[c] + f3*ad3[c]) for c in range (3) ] bx = [ (ax[c] + length*bd3[c]) for c in range(3) ] - bd1 = vector.crossproduct3(ad2, bd3) - bd2 = vector.crossproduct3(bd3, bd1) + bd1 = cross(ad2, bd3) + bd2 = cross(bd3, bd1) return bx, bd1, bd2, bd3 def createEllipsePoints(cx, radian, axis1, axis2, elementsCountAround, startRadians = 0.0): @@ -488,16 +487,16 @@ def createEllipsePoints(cx, radian, axis1, axis2, elementsCountAround, startRadi ''' px = [] pd1 = [] - magAxis1 = vector.magnitude(axis1) - magAxis2 = vector.magnitude(axis2) + magAxis1 = magnitude(axis1) + magAxis2 = magnitude(axis2) totalEllipsePerimeter = getApproximateEllipsePerimeter(magAxis1, magAxis2) partOfEllipsePerimeter = radian * totalEllipsePerimeter / (2 * math.pi) elementLength = partOfEllipsePerimeter / elementsCountAround if radian != 2 * math.pi: elementsCountAround = elementsCountAround + 1 - unitSideAxis1 = vector.normalise(axis1) - unitSideAxis2 = vector.normalise(axis2) + unitSideAxis1 = normalize(axis1) + unitSideAxis2 = normalize(axis2) for n in range(elementsCountAround): angle = startRadians arcLength = n * elementLength @@ -511,7 +510,7 @@ def createEllipsePoints(cx, radian, axis1, axis2, elementsCountAround, startRadi ] px.append(x) rd1 = [magAxis1 * (-sinRadiansAround * unitSideAxis1[c]) + magAxis2 * (cosRadiansAround * unitSideAxis2[c]) for c in range(3)] - rd1Norm = vector.normalise(rd1) + rd1Norm = normalize(rd1) pd1.append([elementLength * (rd1Norm[c])for c in range(3)]) return px, pd1 diff --git a/src/scaffoldmaker/utils/interpolation.py b/src/scaffoldmaker/utils/interpolation.py index 54d05de5..b5548976 100644 --- a/src/scaffoldmaker/utils/interpolation.py +++ b/src/scaffoldmaker/utils/interpolation.py @@ -2,13 +2,12 @@ Interpolation functions shared by mesh generators. """ -from cmlibs.maths.vectorops import add, cross, dot, magnitude, mult, normalize, sub +from cmlibs.maths.vectorops import add, cross, dot, magnitude, mult, normalize, sub, set_magnitude import copy import math from collections.abc import Sequence from enum import Enum -from scaffoldmaker.utils import vector gaussXi3 = ( (-math.sqrt(0.6)+1.0)/2.0, 0.5, (+math.sqrt(0.6)+1.0)/2.0 ) gaussWt3 = ( 5.0/18.0, 4.0/9.0, 5.0/18.0 ) @@ -105,8 +104,8 @@ def computeCubicHermiteArcLength(v1, d1, v2, d2, rescaleDerivatives): lastArcLength = math.sqrt(sum((v2[i] - v1[i])*(v2[i] - v1[i]) for i in range(len(v1)))) else: lastArcLength = getCubicHermiteArcLength(v1, d1, v2, d2) - d1 = vector.normalise(d1) - d2 = vector.normalise(d2) + d1 = normalize(d1) + d2 = normalize(d2) tol = 1.0E-6 for iters in range(100): #print('iter',iters,'=',lastArcLength) @@ -237,9 +236,9 @@ def getCubicHermiteCurvature(v1, d1, v2, d2, radialVector, xi): """ tangent = interpolateCubicHermiteDerivative(v1, d1, v2, d2, xi) dTangent = interpolateCubicHermiteSecondDerivative(v1, d1, v2, d2, xi) - #tangentVector = vector.normalise(tangent) - #tangentCurvature = vector.dotproduct(dTangent, tangentVector) - radialCurvature = vector.dotproduct(dTangent, radialVector) + #tangentVector = normalize(tangent) + #tangentCurvature = dot(dTangent, tangentVector) + radialCurvature = dot(dTangent, radialVector) magTangent = magnitude(tangent) curvature = radialCurvature/(magTangent*magTangent) @@ -256,7 +255,7 @@ def getCubicHermiteCurvatureSimple(v1, d1, v2, d2, xi): mag_tangent = magnitude(tangent) if mag_tangent > 0.0: dTangent = interpolateCubicHermiteSecondDerivative(v1, d1, v2, d2, xi) - cp = vector.crossproduct3(tangent, dTangent) + cp = cross(tangent, dTangent) curvature = magnitude(cp) / (mag_tangent * mag_tangent * mag_tangent) else: curvature = 0.0 @@ -331,7 +330,7 @@ def projectHermiteCurvesThroughWall(nx, nd1, nd2, n, wallThickness, loop=False): """ maxPointIndex = len(nx) - 1 assert (0 <= n <= maxPointIndex), 'projectHermiteCurvesThroughWall. Invalid index' - unitNormal = vector.normalise(vector.crossproduct3(nd1[n], nd2[n])) + unitNormal = normalize(cross(nd1[n], nd2[n])) x = [ (nx[n][c] + wallThickness*unitNormal[c]) for c in range(3) ] # calculate inner d1 from curvature around curvature = 0.0 @@ -346,7 +345,7 @@ def projectHermiteCurvesThroughWall(nx, nd1, nd2, n, wallThickness, loop=False): factor = 1.0 - curvature*wallThickness d1 = [ factor*c for c in nd1[n] ] d2 = copy.deepcopy(nd2[n]) # magnitude can't be determined here - d3 = vector.setMagnitude(unitNormal, math.fabs(wallThickness)) + d3 = set_magnitude(unitNormal, math.fabs(wallThickness)) return x, d1, d2, d3 @@ -383,8 +382,8 @@ def sampleCubicHermiteCurves(nx, nd1, elementsCountOut, for e in range(elementsCountIn): if arcLengthDerivatives: arcLength = computeCubicHermiteArcLength(nx[e], nd1[e], nx[e + 1], nd1[e + 1], rescaleDerivatives = True) - nd1a.append(vector.setMagnitude(nd1[e], arcLength)) - nd1b.append(vector.setMagnitude(nd1[e + 1], arcLength)) + nd1a.append(set_magnitude(nd1[e], arcLength)) + nd1b.append(set_magnitude(nd1[e + 1], arcLength)) else: arcLength = getCubicHermiteArcLength(nx[e], nd1[e], nx[e + 1], nd1[e + 1]) length += arcLength @@ -733,7 +732,7 @@ def smoothCubicHermiteDerivativesLine(nx, nd1, if fixAllDirections or (fixStartDirection and fixEndDirection): # fixed directions, equal magnitude arcLength = computeCubicHermiteArcLength(nx[0], nd1[0], nx[1], nd1[1], rescaleDerivatives=True) - return [ vector.setMagnitude(nd1[0], arcLength), vector.setMagnitude(nd1[1], arcLength) ] + return [ set_magnitude(nd1[0], arcLength), set_magnitude(nd1[1], arcLength) ] # fixStartMag = magnitude(md1[0]) if fixStartDerivative else 0.0 # fixEndMag = magnitude(md1[-1]) if fixEndDerivative else 0.0 tol = 1.0E-6 @@ -746,7 +745,7 @@ def smoothCubicHermiteDerivativesLine(nx, nd1, if not fixStartDerivative: if fixAllDirections or fixStartDirection: mag = 2.0*arcLengths[0] - magnitude(lastmd1[1]) - md1[0] = vector.setMagnitude(nd1[0], mag) if (mag > 0.0) else [ 0.0, 0.0, 0.0 ] + md1[0] = set_magnitude(nd1[0], mag) if (mag > 0.0) else [ 0.0, 0.0, 0.0 ] else: md1[0] = interpolateLagrangeHermiteDerivative(nx[0], nx[1], lastmd1[1], 0.0) # middle @@ -773,12 +772,12 @@ def smoothCubicHermiteDerivativesLine(nx, nd1, mag = 0.5 * (dnm + dn) else: # harmonicMeanMagnitude mag = 2.0 / (1.0 / dnm + 1.0 / dn) - md1[n] = vector.setMagnitude(md1[n], mag) + md1[n] = set_magnitude(md1[n], mag) # end if not fixEndDerivative: if fixAllDirections or fixEndDirection: mag = 2.0*arcLengths[-1] - magnitude(lastmd1[-2]) - md1[-1] = vector.setMagnitude(nd1[-1], mag) if (mag > 0.0) else [ 0.0, 0.0, 0.0 ] + md1[-1] = set_magnitude(nd1[-1], mag) if (mag > 0.0) else [ 0.0, 0.0, 0.0 ] else: md1[-1] = interpolateHermiteLagrangeDerivative(nx[-2], lastmd1[-2], nx[-1], 1.0) if instrument: @@ -1002,7 +1001,7 @@ def smoothCubicHermiteDerivativesLoop(nx, nd1, mag = 0.5*(arcLengths[nm] + arcLengths[n]) else: # harmonicMeanMagnitude mag = 2.0/(1.0/arcLengths[nm] + 1.0/arcLengths[n]) - md1[n] = vector.setMagnitude(md1[n], mag) + md1[n] = set_magnitude(md1[n], mag) if instrument: print('iter', iter + 1, md1) dtol = tol*sum(arcLengths)/len(arcLengths) @@ -1037,7 +1036,7 @@ def getDoubleCubicHermiteCurvesMidDerivative(ax, ad1, mx, bx, bd1): maga = magnitude(ad1) magb = magnitude(bd1) magm = arcLengtha + arcLengthb - 0.5*(maga + magb) - return vector.setMagnitude(md1, magm) + return set_magnitude(md1, magm) def sampleParameterAlongLine(lengthList, paramList, elementsCountOut): """ diff --git a/src/scaffoldmaker/utils/mirror.py b/src/scaffoldmaker/utils/mirror.py index 94801e94..6c6d166c 100644 --- a/src/scaffoldmaker/utils/mirror.py +++ b/src/scaffoldmaker/utils/mirror.py @@ -4,20 +4,20 @@ from __future__ import division -from scaffoldmaker.utils import vector +from cmlibs.maths.vectorops import magnitude, normalize, dot class Mirror: - ''' + """ Utilities for getting a mirror image of a scaffold. - ''' + """ def __init__(self, mirrorPlane): - ''' + """ :param mirrorPlane: Plane ax+by+cd=d defined as a list p=[a,b,c,d]. - ''' + """ self._plane = mirrorPlane self._planeUnitNormalVector = self.getPlaneNormalVector() - self._magNormal = vector.magnitude(mirrorPlane[:3]) + self._magNormal = magnitude(mirrorPlane[:3]) self._planeDParam = mirrorPlane[3] def getPointDistanceFromPlane(self, x): @@ -29,7 +29,7 @@ def getPointDistanceFromPlane(self, x): n = self._planeUnitNormalVector ma = self._magNormal d = self._planeDParam - return vector.dotproduct(x, n) - d / ma + return dot(x, n) - d / ma def getPlaneNormalVector(self): """ @@ -37,7 +37,7 @@ def getPlaneNormalVector(self): :return: Unit normal vector. """ p = self._plane - return vector.normalise([p[0], p[1], p[2]]) + return normalize([p[0], p[1], p[2]]) def pointProjectionToPlane(self, x): """ @@ -66,7 +66,7 @@ def mirrorVector(self,v): :return: image vector. """ n = self._planeUnitNormalVector - return [(v[c] - 2 * vector.dotproduct(v, n) * n[c]) for c in range(3)] + return [(v[c] - 2 * dot(v, n) * n[c]) for c in range(3)] def reverseMirrorVector(self,v): """ @@ -75,4 +75,4 @@ def reverseMirrorVector(self,v): :return: image vector. """ n = self._planeUnitNormalVector - return [-v[c] + 2 * vector.dotproduct(v, n) * n[c] for c in range(3)] + return [-v[c] + 2 * dot(v, n) * n[c] for c in range(3)] diff --git a/src/scaffoldmaker/utils/shieldmesh.py b/src/scaffoldmaker/utils/shieldmesh.py index 102e50db..c469f2b8 100644 --- a/src/scaffoldmaker/utils/shieldmesh.py +++ b/src/scaffoldmaker/utils/shieldmesh.py @@ -8,10 +8,10 @@ import copy from enum import Enum +from cmlibs.maths.vectorops import cross, normalize from cmlibs.zinc.element import Element from cmlibs.zinc.field import Field from cmlibs.zinc.node import Node -from scaffoldmaker.utils import vector from scaffoldmaker.utils.eft_utils import remapEftNodeValueLabel, setEftScaleFactorIds from scaffoldmaker.utils.eftfactory_tricubichermite import eftfactory_tricubichermite from scaffoldmaker.utils.interpolation import DerivativeScalingMode, sampleCubicHermiteCurves, \ @@ -189,7 +189,7 @@ def getTriplePoints(self, n3): self.pProportions[n2b][n1c]))) self.pProportions[n2b][n1b] = self.trackSurface.getProportion(p) x, sd1, sd2 = self.trackSurface.evaluateCoordinates(p, derivatives=True) - d1, d2, d3 = calculate_surface_axes(sd1, sd2, vector.normalise(sd1)) + d1, d2, d3 = calculate_surface_axes(sd1, sd2, normalize(sd1)) self.pd3[n3][n2b][n1b] = d3 self.px[n3][n2b][n1b] = x if self._type == ShieldRimDerivativeMode.SHIELD_RIM_DERIVATIVE_MODE_AROUND: @@ -200,10 +200,10 @@ def getTriplePoints(self, n3): self.pd2[n3][n2b][n1b] = [(self.px[n3][n2c][n1b][c] - self.px[n3][n2b][n1b][c]) for c in range(3)] if not self.trackSurface: if self._type == ShieldRimDerivativeMode.SHIELD_RIM_DERIVATIVE_MODE_AROUND: - self.pd2[n3][n2b][n1b] = vector.normalise(vector.crossproduct3(self.pd3[n3][n2b][n1b], + self.pd2[n3][n2b][n1b] = normalize(cross(self.pd3[n3][n2b][n1b], self.pd1[n3][n2b][n1b])) elif self._type == ShieldRimDerivativeMode.SHIELD_RIM_DERIVATIVE_MODE_REGULAR: - self.pd3[n3][n2b][n1b] = vector.normalise(vector.crossproduct3(self.pd1[n3][n2b][n1b], + self.pd3[n3][n2b][n1b] = normalize(cross(self.pd1[n3][n2b][n1b], self.pd2[n3][n2b][n1b])) # right rtx = [] @@ -246,7 +246,7 @@ def getTriplePoints(self, n3): self.pProportions[n2b][m1c]))) self.pProportions[n2b][m1b] = self.trackSurface.getProportion(p) x, sd1, sd2 = self.trackSurface.evaluateCoordinates(p, derivatives=True) - d1, d2, d3 = calculate_surface_axes(sd1, sd2, vector.normalise(sd1)) + d1, d2, d3 = calculate_surface_axes(sd1, sd2, normalize(sd1)) self.pd3[n3][n2b][m1b] = d3 self.px[n3][n2b][m1b] = x if self._type == ShieldRimDerivativeMode.SHIELD_RIM_DERIVATIVE_MODE_AROUND: @@ -257,10 +257,10 @@ def getTriplePoints(self, n3): self.pd2[n3][n2b][m1b] = [(self.px[n3][n2c][m1b][c] - self.px[n3][n2b][m1b][c]) for c in range(3)] if not self.trackSurface: if self._type == ShieldRimDerivativeMode.SHIELD_RIM_DERIVATIVE_MODE_AROUND: - self.pd2[n3][n2b][m1b] = vector.normalise(vector.crossproduct3(self.pd3[n3][n2b][m1b], + self.pd2[n3][n2b][m1b] = normalize(cross(self.pd3[n3][n2b][m1b], self.pd1[n3][n2b][m1b])) elif self._type == ShieldRimDerivativeMode.SHIELD_RIM_DERIVATIVE_MODE_REGULAR: - self.pd3[n3][n2b][m1b] = vector.normalise(vector.crossproduct3(self.pd1[n3][n2b][m1b], + self.pd3[n3][n2b][m1b] = normalize(cross(self.pd1[n3][n2b][m1b], self.pd2[n3][n2b][m1b])) def smoothDerivativesToTriplePoints(self, n3, fixAllDirections=False): diff --git a/src/scaffoldmaker/utils/spheremesh.py b/src/scaffoldmaker/utils/spheremesh.py index 896478a8..18f3cb0a 100644 --- a/src/scaffoldmaker/utils/spheremesh.py +++ b/src/scaffoldmaker/utils/spheremesh.py @@ -6,9 +6,9 @@ import math +from cmlibs.maths.vectorops import magnitude, mult, normalize, add_vectors, set_magnitude, rejection, cross, angle, rotate_vector_around_vector from cmlibs.utils.zinc.finiteelement import getMaximumNodeIdentifier, getMaximumElementIdentifier from cmlibs.zinc.field import Field -from scaffoldmaker.utils import vector from scaffoldmaker.utils.interpolation import sampleCubicHermiteCurves, smoothCubicHermiteDerivativesLine from scaffoldmaker.utils.cylindermesh import Ellipse2D, EllipseShape from scaffoldmaker.utils.shieldmesh import ShieldMesh3D @@ -57,7 +57,7 @@ def __init__(self, fieldModule, coordinates, centre, axes, elementsCountAcross, """ self._axes = axes - self._radius = [vector.magnitude(axis) for axis in axes] + self._radius = [magnitude(axis) for axis in axes] self._coreRadius = [] self._shield = None self._elementsCount = [ec - 2 * elementsCountAcrossShell for ec in elementsCountAcross] @@ -198,11 +198,11 @@ def get_octant_axes_and_elements_count(self, octant_shape): if octant_number in [2, 4, 5, 7]: boxDerivatives = [boxDerivativesV2[c] * signs[c] for c in range(3)] elementsCountAcross = elementsCountAcrossV2 - axes = [vector.scaleVector(axesV2[c], axesSignsV2[c]) for c in range(3)] + axes = [mult(axesV2[c], axesSignsV2[c]) for c in range(3)] else: boxDerivatives = [boxDerivativesV1[c] * signs[c] for c in range(3)] elementsCountAcross = elementsCountAcrossV1 - axes = [vector.scaleVector(axesV1[c], axesSignsV1[c]) for c in range(3)] + axes = [mult(axesV1[c], axesSignsV1[c]) for c in range(3)] if self._sphereShape == SphereShape.SPHERE_SHAPE_HALF_AAP: elementsCountAcross[0] = elementsCountAcross[0]//2 @@ -212,7 +212,7 @@ def get_octant_axes_and_elements_count(self, octant_shape): elementsCountAcross[1] = elementsCountAcross[1] // 2 elementsCountAcross[2] = elementsCountAcross[2] // 2 - axes = [vector.normalise(v) for v in axes] + axes = [normalize(v) for v in axes] return axes, elementsCountAcross, boxDerivatives def copy_octant_nodes_to_sphere_shield(self, octant, octant_shape): @@ -474,7 +474,7 @@ def __init__(self, centre, axes, elementsCountAcross, example, is the octant in positive axis1, positive axis2 and positive axis3. """ self._axes = axes - self._radius = [vector.magnitude(axis) for axis in axes] + self._radius = [magnitude(axis) for axis in axes] self._coreRadius = [] self._shield = None self._elementsCount = elementsCountAcross @@ -815,13 +815,13 @@ def sample_curves_between_two_nodes_on_sphere(self, id1, id2, elementsOut, dStar if ni < len(dStart): if dStart[ni]: - btd[dStart[ni]][idi[0]][idi[1]][idi[2]] = nd1[nit] if dStart[ni] > 0 else vector.scaleVector( + btd[dStart[ni]][idi[0]][idi[1]][idi[2]] = nd1[nit] if dStart[ni] > 0 else mult( nd1[nit], -1) nit += 1 elif ni > elementsCount - len(dEnd): nie = ni - elementsCount + len(dEnd) - 1 if dEnd[nie]: - btd[abs(dEnd[nie])][idi[0]][idi[1]][idi[2]] = nd1[nit] if dEnd[nie] > 0 else vector.scaleVector( + btd[abs(dEnd[nie])][idi[0]][idi[1]][idi[2]] = nd1[nit] if dEnd[nie] > 0 else mult( nd1[nit], -1) nit += 1 else: @@ -832,7 +832,7 @@ def sample_curves_between_two_nodes_on_sphere(self, id1, id2, elementsOut, dStar btd2[idi[0]][idi[1]][idi[2]] = a2 # initialise btd3[idi[0]][idi[1]][idi[2]] = a3 # initialise - btd[abs(dbetween[0])][idi[0]][idi[1]][idi[2]] = nd1[nit] if dbetween[0] > 0 else vector.scaleVector( + btd[abs(dbetween[0])][idi[0]][idi[1]][idi[2]] = nd1[nit] if dbetween[0] > 0 else mult( nd1[nit], -1) nit += 1 @@ -880,7 +880,7 @@ def smooth_derivatives_regular_surface_curve(self, constant_index, nc, dStart, d if dStart[0][ni] > 0: td.append(btd[abs(dStart[se][ni])][ids[0]][ids[1]][ids[2]]) else: - td.append(vector.scaleVector(btd[abs(dStart[se][ni])][ids[0]][ids[1]][ids[2]], -1)) + td.append(mult(btd[abs(dStart[se][ni])][ids[0]][ids[1]][ids[2]], -1)) elif ni > elementsCount[se] - len(dEnd[se]): nie = ni - elementsCount[se] + len(dEnd[se]) - 1 if dEnd[se][nie]: @@ -888,13 +888,13 @@ def smooth_derivatives_regular_surface_curve(self, constant_index, nc, dStart, d if dEnd[se][nie] > 0: td.append(btd[abs(dEnd[se][nie])][ids[0]][ids[1]][ids[2]]) else: - td.append(vector.scaleVector(btd[abs(dEnd[se][nie])][ids[0]][ids[1]][ids[2]], -1)) + td.append(mult(btd[abs(dEnd[se][nie])][ids[0]][ids[1]][ids[2]], -1)) else: tx.append(btx[ids[0]][ids[1]][ids[2]]) if dBetween[se][0] > 0: td.append(btd[abs(dBetween[se][0])][ids[0]][ids[1]][ids[2]]) else: - td.append(vector.scaleVector(btd[abs(dBetween[se][0])][ids[0]][ids[1]][ids[2]], -1)) + td.append(mult(btd[abs(dBetween[se][0])][ids[0]][ids[1]][ids[2]], -1)) td = smoothCubicHermiteDerivativesLine(tx, td, fixStartDirection=True, fixEndDirection=True) @@ -913,7 +913,7 @@ def smooth_derivatives_regular_surface_curve(self, constant_index, nc, dStart, d if dStart[0][ni] > 0: btd[abs(dStart[se][ni])][ids[0]][ids[1]][ids[2]] = td[nit] else: - btd[abs(dStart[se][ni])][ids[0]][ids[1]][ids[2]] = vector.scaleVector(td[nit], -1) + btd[abs(dStart[se][ni])][ids[0]][ids[1]][ids[2]] = mult(td[nit], -1) nit += 1 elif ni > elementsCount[se] - len(dEnd[se]): nie = ni - elementsCount[se] + len(dEnd[se]) - 1 @@ -921,13 +921,13 @@ def smooth_derivatives_regular_surface_curve(self, constant_index, nc, dStart, d if dEnd[se][nie] > 0: btd[abs(dEnd[se][nie])][ids[0]][ids[1]][ids[2]] = td[nit] else: - btd[abs(dEnd[se][nie])][ids[0]][ids[1]][ids[2]] = vector.scaleVector(td[nit], -1) + btd[abs(dEnd[se][nie])][ids[0]][ids[1]][ids[2]] = mult(td[nit], -1) nit += 1 else: if dBetween[se][0] > 0: btd[abs(dBetween[se][0])][ids[0]][ids[1]][ids[2]] = td[nit] else: - btd[abs(dBetween[se][0])][ids[0]][ids[1]][ids[2]] = vector.scaleVector(td[nit], -1) + btd[abs(dBetween[se][0])][ids[0]][ids[1]][ids[2]] = mult(td[nit], -1) nit += 1 def calculate_interior_quadruple_point(self): @@ -955,11 +955,11 @@ def calculate_interior_quadruple_point(self): n3r0, n2r0, n1r0 = self.get_triple_curves_end_node_parameters(0, cx=cx, index_output=True) n3r, n2r, n1r = self.get_triple_curves_end_node_parameters(1, cx=cx, index_output=True) - ts = vector.magnitude(vector.addVectors([btx[n3r0][n2r0][n1r0], btx[n3r][n2r][n1r]], [1, -1])) - ra = vector.addVectors([btx[n3z][0][n1z], self._centre], [1, -1]) - radius = vector.magnitude(ra) - local_x = vector.scaleVector(ra, (1 - ts/radius)) - x = vector.addVectors([local_x, self._centre], [1, 1]) + ts = magnitude(add_vectors([btx[n3r0][n2r0][n1r0], btx[n3r][n2r][n1r]], [1, -1])) + ra = add_vectors([btx[n3z][0][n1z], self._centre], [1, -1]) + radius = magnitude(ra) + local_x = mult(ra, (1 - ts/radius)) + x = add_vectors([local_x, self._centre], [1, 1]) n3r0, n2r0, n1r0 = self.get_triple_curves_end_node_parameters(0, index_output=True) n3r1, n2r1, n1r1 = self.get_triple_curves_end_node_parameters(1, index_output=True) btx[n3r0][n2r0][n1r0] = x @@ -1007,8 +1007,8 @@ def sample_interior_curves(self): btd1[n3y][n2][n1] = [btx[n3y][1][n1][0] - btx[n3z][0][n1][0], 0.0, 0.0] btd2[n3y][n2][n1] = [0.0, 0.0, -btx[n3y][1][n1][2] + btx[n3z][0][n1][2]] else: - btd1[n3y][n2][n1] = vector.addVectors([btx[n3y][n2][n1], btx[n3y][n2+1][n1]], [-1, 1]) - btd2[n3y][n2][n1] = vector.addVectors([btx[n3y][n2][n1], btx[0][n2][n1]], [1, -1]) + btd1[n3y][n2][n1] = add_vectors([btx[n3y][n2][n1], btx[n3y][n2+1][n1]], [-1, 1]) + btd2[n3y][n2][n1] = add_vectors([btx[n3y][n2][n1], btx[0][n2][n1]], [1, -1]) # sample along curve0_3 for n2 in range(1, self._elementsCount[0]): @@ -1027,8 +1027,8 @@ def sample_interior_curves(self): btd1[n3][n2][n1] = [btx[n3][n2][n1][0] - btx[n3][n2-1][n1+1][0], 0.0, 0.0] btd3[n3][n2][n1] = [0.0, btx[n3][n2-1][n1+1][1] - btx[n3][n2][n1][1], 0.0] else: - btd1[n3][n2][n1] = vector.addVectors([btx[n3][n2+1][n1], btx[n3][n2][n1]], [1, -1]) - btd3[n3][n2][n1] = vector.addVectors([btx[n3][n2][n1+1], btx[n3][n2][n1]], [1, -1]) + btd1[n3][n2][n1] = add_vectors([btx[n3][n2+1][n1], btx[n3][n2][n1]], [1, -1]) + btd3[n3][n2][n1] = add_vectors([btx[n3][n2][n1+1], btx[n3][n2][n1]], [1, -1]) def smooth_regular_interior_curves(self): """ @@ -1061,8 +1061,8 @@ def smooth_regular_interior_curves(self): else: for n3 in range(1, self._elementsCount[2]): for n1 in range(1, self._elementsCount[1]): - btd1[n3][1][n1] = vector.addVectors([btx[n3][2][n1], btx[n3][1][n1]], [1, -1]) - btd1[n3][2][n1] = vector.setMagnitude(btd1[n3][2][n1], vector.magnitude(btd1[n3][1][n1])) + btd1[n3][1][n1] = add_vectors([btx[n3][2][n1], btx[n3][1][n1]], [1, -1]) + btd1[n3][2][n1] = set_magnitude(btd1[n3][2][n1], magnitude(btd1[n3][1][n1])) # smooth d3 in regular if self._elementsCount[1] >= 3: @@ -1081,8 +1081,8 @@ def smooth_regular_interior_curves(self): else: for n3 in range(1, self._elementsCount[2]): for n2 in range(1, self._elementsCount[0]): - btd3[n3][n2][1] = vector.addVectors([btx[n3][n2][1], btx[n3][n2][0]], [1, -1]) - btd3[n3][n2][0] = vector.setMagnitude(btd3[n3][n2][0], vector.magnitude(btd3[n3][n2][1])) + btd3[n3][n2][1] = add_vectors([btx[n3][n2][1], btx[n3][n2][0]], [1, -1]) + btd3[n3][n2][0] = set_magnitude(btd3[n3][n2][0], magnitude(btd3[n3][n2][1])) # regular curves d2 for n2 in range(1, self._elementsCount[0]): @@ -1097,8 +1097,8 @@ def smooth_regular_interior_curves(self): for n3 in range(self._elementsCount[2]): btd2[n3][n2][n1] = td2[n3] else: - btd2[1][n2][n1] = vector.addVectors([btx[1][n2][n1], btx[0][n2][n1]], [1, -1]) - btd2[0][n2][n1] = vector.setMagnitude(btd2[0][n2][n1], vector.magnitude(btd2[1][n2][n1])) + btd2[1][n2][n1] = add_vectors([btx[1][n2][n1], btx[0][n2][n1]], [1, -1]) + btd2[0][n2][n1] = set_magnitude(btd2[0][n2][n1], magnitude(btd2[1][n2][n1])) def get_triple_curves_end_node_parameters(self, rx, cx=None, index_output=False): """ @@ -1182,11 +1182,11 @@ def local_orthogonal_unit_vectors(x, axis3, origin): :param axis3: The third axis in Cartesian coordinate system (axis1, axis2, axis3) :return: e1, e2, e3. Unit vectors. e3 is normal to the boundary, e2 is in (e3, axis3) plane and e1 normal to them. """ - r = vector.addVectors([x, origin], [1, -1]) - e3 = vector.normalise(r) - e2 = vector.vectorRejection(axis3, e3) - e2 = vector.normalise(e2) - e1 = vector.crossproduct3(e2, e3) + r = add_vectors([x, origin], [1, -1]) + e3 = normalize(r) + e2 = rejection(axis3, e3) + e2 = normalize(e2) + e1 = cross(e2, e3) return e1, e2, e3 @@ -1197,11 +1197,11 @@ def calculate_arc_length(x1, x2, origin): :param x1, x2: points coordinates. :return: arc length """ - r1 = vector.addVectors([x1, origin], [1, -1]) - r2 = vector.addVectors([x2, origin], [1, -1]) - radius = vector.magnitude(r1) - angle = vector.angleBetweenVectors(r1, r2) - return radius * angle + r1 = add_vectors([x1, origin], [1, -1]) + r2 = add_vectors([x2, origin], [1, -1]) + radius = magnitude(r1) + theta = angle(r1, r2) + return radius * theta def sample_curves_on_sphere(x1, x2, origin, elementsOut): @@ -1211,11 +1211,11 @@ def sample_curves_on_sphere(x1, x2, origin, elementsOut): :param elementsOut: :return: """ - r1 = vector.addVectors([x1, origin], [1, -1]) - r2 = vector.addVectors([x2, origin], [1, -1]) - deltax = vector.addVectors([r1, r2], [-1, 1]) - normal = vector.crossproduct3(r1, deltax) - angle = vector.angleBetweenVectors(r1, r2) + r1 = add_vectors([x1, origin], [1, -1]) + r2 = add_vectors([x2, origin], [1, -1]) + deltax = add_vectors([r1, r2], [-1, 1]) + normal = cross(r1, deltax) + angle = angle(r1, r2) anglePerElement = angle/elementsOut arcLengthPerElement = calculate_arc_length(x1, x2, origin)/elementsOut @@ -1223,9 +1223,9 @@ def sample_curves_on_sphere(x1, x2, origin, elementsOut): nd1 = [] for n1 in range(elementsOut + 1): radiansAcross = n1 * anglePerElement - r = vector.rotateVectorAroundVector(r1, normal, radiansAcross) - x = vector.addVectors([r, origin], [1, 1]) - d1 = vector.setMagnitude(vector.crossproduct3(normal, r), arcLengthPerElement) + r = rotate_vector_around_vector(r1, normal, radiansAcross) + x = add_vectors([r, origin], [1, 1]) + d1 = set_magnitude(cross(normal, r), arcLengthPerElement) nx.append(x) nd1.append(d1) @@ -1246,7 +1246,7 @@ def cartesian_to_spherical(x): """ :return: [r, theta, phi]. """ - r = vector.magnitude(x) + r = magnitude(x) theta = math.atan2(x[1], x[0]) phi = math.acos(x[2]/r) return r, theta, phi @@ -1264,8 +1264,8 @@ def local_to_global_coordinates(local_x, local_axes, local_origin=None): """ if local_origin is None: local_origin = [0.0, 0.0, 0.0] - normalised_axes = [vector.normalise(v) for v in local_axes] - return vector.addVectors([vector.addVectors(normalised_axes, local_x), local_origin]) + normalised_axes = [normalize(v) for v in local_axes] + return add_vectors([add_vectors(normalised_axes, local_x), local_origin]) def intersection_of_two_great_circles_on_sphere(p1, q1, p2, q2): @@ -1274,19 +1274,19 @@ def intersection_of_two_great_circles_on_sphere(p1, q1, p2, q2): :param p1, q1, p2, q2: arcs extremities coordinates. :return: Point Sx, intersection between the arcs. """ - normal_to_plane_OP1Q1 = vector.crossproduct3(p1, q1) - normal_to_plane_OP2Q2 = vector.crossproduct3(p2, q2) + normal_to_plane_OP1Q1 = cross(p1, q1) + normal_to_plane_OP2Q2 = cross(p2, q2) - planes_intersection_vector = vector.crossproduct3(normal_to_plane_OP1Q1, normal_to_plane_OP2Q2) - if vector.magnitude(planes_intersection_vector) == 0: + planes_intersection_vector = cross(normal_to_plane_OP1Q1, normal_to_plane_OP2Q2) + if magnitude(planes_intersection_vector) == 0: sx = None else: - sx = vector.setMagnitude(planes_intersection_vector, vector.magnitude(p1)) - p1q1_angle = vector.angleBetweenVectors(p1, q1) - p1s_angle = vector.angleBetweenVectors(p1, sx) - p2s_angle = vector.angleBetweenVectors(p2, sx) + sx = set_magnitude(planes_intersection_vector, magnitude(p1)) + p1q1_angle = angle(p1, q1) + p1s_angle = angle(p1, sx) + p2s_angle = angle(p2, sx) if p1s_angle > p1q1_angle or p2s_angle > p1q1_angle: - sx = vector.scaleVector(sx, -1) + sx = mult(sx, -1) return sx @@ -1297,4 +1297,4 @@ def point_projection_on_sphere(p1, radius): :param p1: point. :return: """ - return vector.setMagnitude(p1, radius) + return set_magnitude(p1, radius) diff --git a/src/scaffoldmaker/utils/tracksurface.py b/src/scaffoldmaker/utils/tracksurface.py index 936b7c2b..e4f93bb2 100644 --- a/src/scaffoldmaker/utils/tracksurface.py +++ b/src/scaffoldmaker/utils/tracksurface.py @@ -5,7 +5,7 @@ import copy from enum import Enum import math -from cmlibs.maths.vectorops import add, cross, dot, magnitude, mult, normalize, sub +from cmlibs.maths.vectorops import add, cross, dot, magnitude, mult, normalize, sub, set_magnitude from cmlibs.utils.zinc.general import ChangeManager from cmlibs.utils.zinc.field import find_or_create_field_coordinates, find_or_create_field_group from cmlibs.utils.zinc.finiteelement import get_maximum_element_identifier, get_maximum_node_identifier @@ -18,7 +18,6 @@ interpolateHermiteLagrangeDerivative, interpolateLagrangeHermiteDerivative, sampleCubicHermiteCurves, \ sampleCubicHermiteCurvesSmooth, smoothCubicHermiteDerivativesLine, smoothCubicHermiteDerivativesLoop, \ updateCurveLocationToFaceNumber -from scaffoldmaker.utils import vector class TrackSurfacePosition: @@ -316,12 +315,12 @@ def createHermiteCurvePoints(self, aProportion1, aProportion2, bProportion1, bPr f1 = dproportions[n][0] * self._elementsCount1 f2 = dproportions[n][1] * self._elementsCount2 d1 = [(f1*sd1[c] + f2*sd2[c]) for c in range(3)] - d3 = vector.crossproduct3(sd1, sd2) + d3 = cross(sd1, sd2) # handle zero magnitude of d3 mag = math.sqrt(sum(d3[c]*d3[c] for c in range(3))) if mag > 0.0: d3 = [(d3[c]/mag) for c in range(3)] - d2 = vector.crossproduct3(d3, d1) + d2 = cross(d3, d1) nx .append(x) nd2.append(d2) nd1.append(d1) @@ -339,19 +338,19 @@ def resampleHermiteCurvePointsSmooth(self, nx, nd1, nd2, nd3, nProportions, # print(nx, nd1, elementsCount, derivativeMagnitudeStart, derivativeMagnitudeEnd) nx, nd1 = sampleCubicHermiteCurvesSmooth( nx, nd1, elementsCount, derivativeMagnitudeStart, derivativeMagnitudeEnd)[0:2] - mag2 = vector.magnitude(nd2[0]) + mag2 = magnitude(nd2[0]) if mag2 > 0.0: - nd2[0] = vector.setMagnitude(nd2[0], vector.magnitude(nd1[0])) + nd2[0] = set_magnitude(nd2[0], magnitude(nd1[0])) for n in range(1, elementsCount): p = self.findNearestPosition(nx[n], self.createPositionProportion(*nProportions[n])) nProportions[n] = self.getProportion(p) _, sd1, sd2 = self.evaluateCoordinates(p, derivatives=True) - _, d2, d3 = calculate_surface_axes(sd1, sd2, vector.normalise(nd1[n])) - nd2[n] = vector.setMagnitude(d2, vector.magnitude(nd1[n])) + _, d2, d3 = calculate_surface_axes(sd1, sd2, normalize(nd1[n])) + nd2[n] = set_magnitude(d2, magnitude(nd1[n])) nd3[n] = d3 - mag2 = vector.magnitude(nd2[-1]) + mag2 = magnitude(nd2[-1]) if mag2 > 0.0: - nd2[-1] = vector.setMagnitude(nd2[-1], vector.magnitude(nd1[-1])) + nd2[-1] = set_magnitude(nd2[-1], magnitude(nd1[-1])) return nx, nd1, nd2, nd3, nProportions def positionOnBoundary(self, position: TrackSurfacePosition): @@ -1293,8 +1292,8 @@ def trackVector(self, startPosition, direction, trackDistance): ad = [proportion*(adxi1*ad1[c] + adxi2*ad2[c]) for c in range(3)] bd = [proportion*(bdxi1*bd1[c] + bdxi2*bd2[c]) for c in range(3)] arcLength = getCubicHermiteArcLength(ax, ad, bx, bd) - # print('scales', vector.magnitude([ (bx[c] - ax[c]) for c in range(3) ]), - # vector.magnitude(ad), vector.magnitude(bd), 'arc length', arcLength) + # print('scales', magnitude([ (bx[c] - ax[c]) for c in range(3) ]), + # magnitude(ad), magnitude(bd), 'arc length', arcLength) if (distance + arcLength) >= distanceLimit: # limit to useTrackDistance, approximately, and finish r = proportion*(useTrackDistance - distance)/arcLength @@ -1475,17 +1474,17 @@ def calculate_surface_delta_xi(d1, d2, direction): delta_xi2 = inva[1][0]*b[0] + inva[1][1]*b[1] else: # at pole: assume direction is inline with d1 or d2 and other is zero - delta_xi2 = vector.dotproduct(d2, direction) + delta_xi2 = dot(d2, direction) if math.fabs(delta_xi2) > 0.0: delta_xi1 = 0.0 - delta_xi2 = (1.0 if (delta_xi2 > 0.0) else -1.0)*vector.magnitude(direction)/vector.magnitude(d2) + delta_xi2 = (1.0 if (delta_xi2 > 0.0) else -1.0)*magnitude(direction)/magnitude(d2) else: - delta_xi1 = vector.dotproduct(d1, direction) + delta_xi1 = dot(d1, direction) if math.fabs(delta_xi1) > 0.0: - delta_xi1 = (1.0 if (delta_xi1 > 0.0) else -1.0)*vector.magnitude(direction)/vector.magnitude(d1) + delta_xi1 = (1.0 if (delta_xi1 > 0.0) else -1.0)*magnitude(direction)/magnitude(d1) delta_xi2 = 0.0 # delx = [ (delta_xi1*d1[c] + delta_xi2*d2[c]) for c in range(3) ] - # print('delx', delx, 'dir', direction, 'diff', vector.magnitude([ (delx[c] - direction[c]) for c in range(3) ])) + # print('delx', delx, 'dir', direction, 'diff', magnitude([ (delx[c] - direction[c]) for c in range(3) ])) return delta_xi1, delta_xi2 @@ -1496,12 +1495,12 @@ def calculate_surface_axes(d1, d2, direction): Vectors all have unit magnitude. """ delta_xi1, delta_xi2 = calculate_surface_delta_xi(d1, d2, direction) - ax1 = vector.normalise([delta_xi1*d1[c] + delta_xi2*d2[c] for c in range(3)]) - ax3 = vector.crossproduct3(d1, d2) - mag3 = vector.magnitude(ax3) + ax1 = normalize([delta_xi1*d1[c] + delta_xi2*d2[c] for c in range(3)]) + ax3 = cross(d1, d2) + mag3 = magnitude(ax3) if mag3 > 0.0: ax3 = [s/mag3 for s in ax3] - ax2 = vector.normalise(vector.crossproduct3(ax3, ax1)) + ax2 = normalize(cross(ax3, ax1)) else: ax3 = [0.0, 0.0, 0.0] ax2 = [0.0, 0.0, 0.0] diff --git a/src/scaffoldmaker/utils/tubemesh.py b/src/scaffoldmaker/utils/tubemesh.py index e8e356dd..e8c04ccc 100644 --- a/src/scaffoldmaker/utils/tubemesh.py +++ b/src/scaffoldmaker/utils/tubemesh.py @@ -6,6 +6,7 @@ import math +from cmlibs.maths.vectorops import normalize, dot, cross, magnitude, set_magnitude from cmlibs.utils.zinc.field import findOrCreateFieldCoordinates from cmlibs.zinc.element import Element from cmlibs.zinc.field import Field @@ -13,7 +14,6 @@ from scaffoldmaker.annotation.annotationgroup import mergeAnnotationGroups from scaffoldmaker.utils import interpolation as interp from scaffoldmaker.utils import matrix -from scaffoldmaker.utils import vector from scaffoldmaker.utils.eftfactory_bicubichermitelinear import eftfactory_bicubichermitelinear from scaffoldmaker.utils.eftfactory_tricubichermite import eftfactory_tricubichermite from scaffoldmaker.utils.eft_utils import remapEftNodeValueLabelsVersion @@ -77,10 +77,10 @@ def getPlaneProjectionOnCentralPath(x, elementsCountAround, elementsCountAlong, sd2ProjectedListRef = [] for n in range(len(sd2RefList)): - sd1Normalised = vector.normalise(sd1RefList[n]) - dp = vector.dotproduct(sd2RefList[n], sd1Normalised) + sd1Normalised = normalize(sd1RefList[n]) + dp = dot(sd2RefList[n], sd1Normalised) dpScaled = [dp * c for c in sd1Normalised] - sd2Projected = vector.normalise([sd2RefList[n][c] - dpScaled[c] for c in range(3)]) + sd2Projected = normalize([sd2RefList[n][c] - dpScaled[c] for c in range(3)]) sd2ProjectedListRef.append(sd2Projected) return sxRefList, sd1RefList, sd2ProjectedListRef, zRefList @@ -118,12 +118,12 @@ def warpSegmentPoints(xList, d1List, d2List, segmentAxis, sx, sd1, sd2, elements centroid = [0.0, 0.0, refPointZ[nAlongSegment]] # Rotate to align segment axis with tangent of central line - unitTangent = vector.normalise(sd1[nAlongSegment]) - cp = vector.crossproduct3(segmentAxis, unitTangent) - dp = vector.dotproduct(segmentAxis, unitTangent) - if vector.magnitude(cp)> 0.0: # path tangent not parallel to segment axis - axisRot = vector.normalise(cp) - thetaRot = math.acos(vector.dotproduct(segmentAxis, unitTangent)) + unitTangent = normalize(sd1[nAlongSegment]) + cp = cross(segmentAxis, unitTangent) + dp = dot(segmentAxis, unitTangent) + 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) centroidRot = [rotFrame[j][0]*centroid[0] + rotFrame[j][1]*centroid[1] + rotFrame[j][2]*centroid[2] for j in range(3)] @@ -145,7 +145,7 @@ def warpSegmentPoints(xList, d1List, d2List, segmentAxis, sx, sd1, sd2, elements d1 = d1ElementAlongSegment[n1] d2 = d2ElementAlongSegment[n1] - if vector.magnitude(cp)> 0.0: # path tangent not parallel to segment axis + if magnitude(cp)> 0.0: # path tangent not parallel to segment axis xRot1 = [rotFrame[j][0]*x[0] + rotFrame[j][1]*x[1] + rotFrame[j][2]*x[2] for j in range(3)] d1Rot1 = [rotFrame[j][0]*d1[0] + rotFrame[j][1]*d1[1] + rotFrame[j][2]*d1[2] for j in range(3)] d2Rot1 = [rotFrame[j][0]*d2[0] + rotFrame[j][1]*d2[1] + rotFrame[j][2]*d2[2] for j in range(3)] @@ -159,13 +159,13 @@ def warpSegmentPoints(xList, d1List, d2List, segmentAxis, sx, sd1, sd2, elements if n1 == 0: # Find angle between xCentroidRot and first node in the face vectorToFirstNode = [xRot1[c] - centroidRot[c] for c in range(3)] - if vector.magnitude(vectorToFirstNode) > 0.0: - cp = vector.crossproduct3(vector.normalise(vectorToFirstNode), vector.normalise(sd2[nAlongSegment])) - if vector.magnitude(cp) > 1e-7: - cp = vector.normalise(cp) - signThetaRot2 = vector.dotproduct(unitTangent, cp) + if magnitude(vectorToFirstNode) > 0.0: + cp = cross(normalize(vectorToFirstNode), normalize(sd2[nAlongSegment])) + if magnitude(cp) > 1e-7: + cp = normalize(cp) + signThetaRot2 = dot(unitTangent, cp) thetaRot2 = math.acos( - vector.dotproduct(vector.normalise(vectorToFirstNode), vector.normalise(sd2[nAlongSegment]))) + dot(normalize(vectorToFirstNode), normalize(sd2[nAlongSegment]))) axisRot2 = unitTangent rotFrame2 = matrix.getRotationMatrixFromAxisAngle(axisRot2, signThetaRot2*thetaRot2) else: @@ -188,13 +188,13 @@ def warpSegmentPoints(xList, d1List, d2List, segmentAxis, sx, sd1, sd2, elements for n1 in range(elementsCountAround): n = nAlongSegment * elementsCountAround + n1 # Calculate norm - sd1Normalised = vector.normalise(sd1[nAlongSegment]) + sd1Normalised = normalize(sd1[nAlongSegment]) v = [xWarpedList[n][c] - sx[nAlongSegment][c] for c in range(3)] - dp = vector.dotproduct(v, sd1Normalised) + dp = dot(v, sd1Normalised) dpScaled = [dp * c for c in sd1Normalised] vProjected = [v[c] - dpScaled[c] for c in range(3)] - if vector.magnitude(vProjected) > 0.0: - vProjectedNormalised = vector.normalise(vProjected) + if magnitude(vProjected) > 0.0: + vProjectedNormalised = normalize(vProjected) else: vProjectedNormalised = [0.0, 0.0, 0.0] @@ -212,7 +212,7 @@ def warpSegmentPoints(xList, d1List, d2List, segmentAxis, sx, sd1, sd2, elements vProjectedNormalised, 0.0)) # Scale if nAlongSegment < elementsCountAlongSegment: - factor = 1.0 - curvature * vector.magnitude(v) + factor = 1.0 - curvature * magnitude(v) d2 = [factor * c for c in d2WarpedList[n]] d2WarpedListScaled.append(d2) else: @@ -237,8 +237,8 @@ def warpSegmentPoints(xList, d1List, d2List, segmentAxis, sx, sd1, sd2, elements # Calculate unit d3 for n in range(len(xWarpedList)): - d3Unit = vector.normalise(vector.crossproduct3(vector.normalise(d1WarpedList[n]), - vector.normalise(d2WarpedListFinal[n]))) + d3Unit = normalize(cross(normalize(d1WarpedList[n]), + normalize(d2WarpedListFinal[n]))) d3WarpedUnitList.append(d3Unit) return xWarpedList, d1WarpedList, d2WarpedListFinal, d3WarpedUnitList @@ -300,9 +300,9 @@ def extrudeSurfaceCoordinates(xSurf, d1Surf, d2Surf, d3Surf, wallThicknessList, prevIdx = n - 1 if (n1 != 0) else (n2 + 1)*elementsCountAround - 1 nextIdx = n + 1 if (n1 < elementsCountAround - 1) else n2*elementsCountAround kappam = interp.getCubicHermiteCurvature(xSurf[prevIdx], d1Surf[prevIdx], xSurf[n], d1Surf[n], - vector.normalise(d3Surf[n]), 1.0) + normalize(d3Surf[n]), 1.0) kappap = interp.getCubicHermiteCurvature(xSurf[n], d1Surf[n], xSurf[nextIdx], d1Surf[nextIdx], - vector.normalise(d3Surf[n]), 0.0) + normalize(d3Surf[n]), 0.0) if not transitElementList[n1] and not transitElementList[(n1-1)%elementsCountAround]: curvatureAround = 0.5*(kappam + kappap) elif transitElementList[n1]: @@ -315,18 +315,18 @@ def extrudeSurfaceCoordinates(xSurf, d1Surf, d2Surf, d3Surf, wallThicknessList, if n2 == 0: curvature = interp.getCubicHermiteCurvature(xSurf[n], d2Surf[n], xSurf[n + elementsCountAround], d2Surf[n + elementsCountAround], - vector.normalise(d3Surf[n]), 0.0) + normalize(d3Surf[n]), 0.0) elif n2 == elementsCountAlong: curvature = interp.getCubicHermiteCurvature(xSurf[n - elementsCountAround], d2Surf[n - elementsCountAround], - xSurf[n], d2Surf[n], vector.normalise(d3Surf[n]), 1.0) + xSurf[n], d2Surf[n], normalize(d3Surf[n]), 1.0) else: curvature = 0.5*( interp.getCubicHermiteCurvature(xSurf[n - elementsCountAround], d2Surf[n - elementsCountAround], - xSurf[n], d2Surf[n], vector.normalise(d3Surf[n]), 1.0) + + xSurf[n], d2Surf[n], normalize(d3Surf[n]), 1.0) + interp.getCubicHermiteCurvature(xSurf[n], d2Surf[n], xSurf[n + elementsCountAround], d2Surf[n + elementsCountAround], - vector.normalise(d3Surf[n]), 0.0)) + normalize(d3Surf[n]), 0.0)) curvatureAlong.append(curvature) for n3 in range(elementsCountThroughWall + 1): @@ -451,12 +451,13 @@ def createFlatCoordinates(xiList, lengthAroundList, totalLengthAlong, wallThickn v2 = xFlatList[nodeIdx] d1 = d2 = [v1[i] - v2[i] for i in range(3)] arclength = interp.computeCubicHermiteArcLength(v1, d1, v2, d2, True) - d2Flat = vector.setMagnitude(d1, arclength) + d2Flat = set_magnitude(d1, arclength) d2FlatList.append(d2Flat) d2FlatList = d2FlatList + d2FlatList[-(elementsCountAround+1)*(elementsCountThroughWall+1):] return xFlatList, d1FlatList, d2FlatList + def createOrganCoordinates(xiList, relativeThicknessList, lengthToDiameterRatio, wallThicknessToDiameterRatio, elementsCountAround, elementsCountAlong, elementsCountThroughWall, transitElementList): """ @@ -507,8 +508,8 @@ def createOrganCoordinates(xiList, relativeThicknessList, lengthToDiameterRatio, # To modify derivative along transition elements for i in range(len(transitElementList)): if transitElementList[i]: - d1List[i] = vector.setMagnitude(d1List[i], vector.magnitude(d1List[i - 1])) - d1List[i + 1] = vector.setMagnitude(d1List[i+ 1], vector.magnitude(d1List[(i + 2) % elementsCountAround])) + d1List[i] = set_magnitude(d1List[i], magnitude(d1List[i - 1])) + d1List[i + 1] = set_magnitude(d1List[i+ 1], magnitude(d1List[(i + 2) % elementsCountAround])) d1OrganList += d1List diff --git a/src/scaffoldmaker/utils/vector.py b/src/scaffoldmaker/utils/vector.py deleted file mode 100644 index 4b794a66..00000000 --- a/src/scaffoldmaker/utils/vector.py +++ /dev/null @@ -1,119 +0,0 @@ -''' -Utility functions for vectors. -''' - -import math - -def crossproduct3(a, b): - ''' - :return: vector 3-D cross product of a and b - ''' - return [ a[1]*b[2] - a[2]*b[1], a[2]*b[0] - a[0]*b[2], a[0]*b[1] - a[1]*b[0] ] - -def dotproduct(a, b): - ''' - :return: vector dot (inner) product of a and b - ''' - return sum(a[i]*b[i] for i in range(len(a))) - -def magnitude(v): - ''' - return: scalar magnitude of vector v - ''' - return math.sqrt(sum(c*c for c in v)) - -def normalise(v): - ''' - :return: vector v normalised to unit length - ''' - mag = math.sqrt(sum(c*c for c in v)) - return [ c/mag for c in v ] - - -def setMagnitude(v, mag): - ''' - return: Vector v with magnitude set to mag. - ''' - oldMag = magnitude(v) - scale = mag / oldMag if (oldMag != 0.0) else 0.0 - return [c *scale for c in v] - - -def addVectors(vectors, scalars = None): - ''' - returns s1*v1+s2*v2+... where scalars = [s1, s2, ...] and vectors=[v1, v2, ...]. - :return: Resultant vector - ''' - if not scalars: - scalars = [1]*len(vectors) - else: - assert len(vectors) == len(scalars) - - resultant = [0, 0, 0] - for i in range(len(vectors)): - resultant = [resultant[c] + scalars[i] * vectors[i][c] for c in range(3)] - return resultant - - -def scalarProjection(v1, v2): - """ - :return: Scalar projection of v1 onto v2. - """ - return dotproduct(v1, normalise(v2)) - - -def vectorProjection(v1, v2): - """ - Calculate vector projection of v1 on v2 - :return: A projection vector. - """ - s1 = scalarProjection(v1, v2) - return scaleVector(normalise(v2), s1) - - -def vectorRejection(v1, v2): - """ - Calculate vector rejection of v1 on v2 - :return: A rejection vector. - """ - v1p = vectorProjection(v1, v2) - return addVectors([v1, v1p], [1.0, -1.0]) - - -def scaleVector(v, s): - """ - Calculate s * v - :param v: Vector. - :param s: Scalar. - :return: - """ - return [s * c for c in v] - - -def parallelVectors(v1, v2): - """ - :return: True if the vectors are parallel. - """ - assert (len(v2) == len(v1)), 'Vectors lengths are not the same.' - TOL = 1.0e-6/2.0 - if magnitude(crossproduct3(v1, v2)) < TOL * (magnitude(v1)+magnitude(v2)): - return True - return False - - -def angleBetweenVectors(v1, v2): - """ - :return: Angle between vectors v1 and v2 in radians - """ - return math.acos(dotproduct(normalise(v1), normalise(v2))) - - -def rotateVectorAroundVector(v, k, a): - """ - Rotate vector v, by an angle a (right-hand rule) in radians around vector k. - :return: rotated vector. - """ - k = normalise(k) - vperp = addVectors([v, crossproduct3(k, v)], [math.cos(a), math.sin(a)]) - vparal = scaleVector(k, dotproduct(k, v)*(1 - math.cos(a))) - return addVectors([vperp, vparal]) diff --git a/src/scaffoldmaker/utils/zinc_utils.py b/src/scaffoldmaker/utils/zinc_utils.py index 66e33f40..9ad92f94 100644 --- a/src/scaffoldmaker/utils/zinc_utils.py +++ b/src/scaffoldmaker/utils/zinc_utils.py @@ -1,6 +1,7 @@ """ Utility functions for easing use of Zinc API. """ +from cmlibs.maths.vectorops import rejection, set_magnitude, magnitude, cross, normalize from cmlibs.utils.zinc.field import find_or_create_field_coordinates, find_or_create_field_group from cmlibs.utils.zinc.finiteelement import get_maximum_element_identifier, get_maximum_node_identifier from cmlibs.utils.zinc.general import ChangeManager, HierarchicalChangeManager @@ -11,7 +12,6 @@ from cmlibs.zinc.node import Node, NodesetGroup from cmlibs.zinc.result import RESULT_OK from scaffoldmaker.utils import interpolation as interp -from scaffoldmaker.utils import vector def interpolateNodesCubicHermite(cache, coordinates, xi, normal_scale, @@ -45,9 +45,9 @@ def interpolateNodesCubicHermite(cache, coordinates, xi, normal_scale, d2c = [ cross_scale2*d for d in d2c ] arcLength = interp.computeCubicHermiteArcLength(v1, d1, v2, d2, True) - mag = arcLength/vector.magnitude(d1) + mag = arcLength/magnitude(d1) d1 = [ mag*d for d in d1 ] - mag = arcLength/vector.magnitude(d2) + mag = arcLength/magnitude(d2) d2 = [ mag*d for d in d2 ] xr = 1.0 - xi @@ -57,7 +57,7 @@ def interpolateNodesCubicHermite(cache, coordinates, xi, normal_scale, dx_ds = [ scale*d for d in dx_ds ] dx_ds_cross = [ (xr*d1c[c] + xi*d2c[c]) for c in range(3) ] - radialVector = vector.normalise(vector.crossproduct3(dx_ds_cross, dx_ds)) + radialVector = normalize(cross(dx_ds_cross, dx_ds)) dx_ds_normal = [ normal_scale*d for d in radialVector ] return x, dx_ds, dx_ds_cross, dx_ds_normal @@ -261,16 +261,16 @@ def make_nodeset_derivatives_orthogonal(nodeset, field, make_d2_normal: bool=Tru d1 = node_parameters[d1_index][v] d2 = node_parameters[d2_index][v] if (d2_index is not None) else None if make_d2_normal: - td2 = vector.vectorRejection(d2, d1) - td2 = vector.setMagnitude(td2, vector.magnitude(d2)) + td2 = rejection(d2, d1) + td2 = set_magnitude(td2, magnitude(d2)) d2 = node_parameters[d2_index][v] = td2 if make_d3_normal: d3 = node_parameters[d3_index][v] if d2: - td3 = vector.crossproduct3(d1, d2) + td3 = cross(d1, d2) else: - td3 = vector.vectorRejection(d3, d1) - td3 = vector.setMagnitude(td3, vector.magnitude(d3)) + td3 = rejection(d3, d1) + td3 = set_magnitude(td3, magnitude(d3)) d3 = node_parameters[d3_index][v] = td3 remove_indexes = [d1_index] if (d2_index is not None) and (not make_d2_normal):