Skip to content

Commit

Permalink
Switch over to using cmlibs.maths 0.6.0 for math functions.
Browse files Browse the repository at this point in the history
  • Loading branch information
hsorby committed Jul 23, 2024
1 parent e15eeec commit 3f60540
Show file tree
Hide file tree
Showing 37 changed files with 656 additions and 807 deletions.
2 changes: 1 addition & 1 deletion setup.py
Original file line number Diff line number Diff line change
Expand Up @@ -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",
Expand Down
38 changes: 19 additions & 19 deletions src/scaffoldmaker/meshtypes/meshtype_3d_bladder1.py
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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])
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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

Expand All @@ -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)
Expand Down Expand Up @@ -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]
Expand Down Expand Up @@ -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]
Expand Down Expand Up @@ -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 = []
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down
23 changes: 11 additions & 12 deletions src/scaffoldmaker/meshtypes/meshtype_3d_bladderurethra1.py
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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 = []
Expand Down Expand Up @@ -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,
Expand All @@ -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,
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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 = []
Expand Down
14 changes: 5 additions & 9 deletions src/scaffoldmaker/meshtypes/meshtype_3d_bone1.py
Original file line number Diff line number Diff line change
Expand Up @@ -7,21 +7,19 @@

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
from cmlibs.zinc.node import Node
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):
Expand Down Expand Up @@ -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]
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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]
Expand Down Expand Up @@ -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
Expand Down
Loading

0 comments on commit 3f60540

Please sign in to comment.