From ad4e5065eef4a9906552b6655cca0389af88342f Mon Sep 17 00:00:00 2001 From: Chang-Joon Lee Date: Mon, 16 Dec 2024 12:23:13 +1300 Subject: [PATCH 1/2] Add renal pelvis scaffold --- src/scaffoldmaker/annotation/kidney_terms.py | 24 + src/scaffoldmaker/annotation/ureter_terms.py | 19 + .../meshtypes/meshtype_3d_renal_pelvis1.py | 683 ++++++++++++++++++ src/scaffoldmaker/scaffolds.py | 6 +- src/scaffoldmaker/utils/tracksurface.py | 1 + src/scaffoldmaker/utils/tubenetworkmesh.py | 40 + tests/test_renalpelvis.py | 131 ++++ 7 files changed, 903 insertions(+), 1 deletion(-) create mode 100644 src/scaffoldmaker/annotation/kidney_terms.py create mode 100644 src/scaffoldmaker/annotation/ureter_terms.py create mode 100644 src/scaffoldmaker/meshtypes/meshtype_3d_renal_pelvis1.py create mode 100644 tests/test_renalpelvis.py diff --git a/src/scaffoldmaker/annotation/kidney_terms.py b/src/scaffoldmaker/annotation/kidney_terms.py new file mode 100644 index 00000000..a0b65cde --- /dev/null +++ b/src/scaffoldmaker/annotation/kidney_terms.py @@ -0,0 +1,24 @@ +""" +Common resource for kidney annotation terms. +""" + +# convention: preferred name, preferred id, followed by any other ids and alternative names +kidney_terms = [ + ("core", ""), + ("renal pelvis", "UBERON:0001224", "ILX:0723968"), + ("major calyx", "UBERON:0001226", "ILX:0730785"), + ("minor calyx", "UBERON:0001227", "ILX:0730473"), + ("renal pyramid", "UBERON:0004200", "ILX:0727514"), + ("shell", "") + ] + +def get_kidney_term(name : str): + """ + Find term by matching name to any identifier held for a term. + Raise exception if name not found. + :return ( preferred name, preferred id ) + """ + for term in kidney_terms: + if name in term: + return (term[0], term[1]) + raise NameError("Kidney annotation term '" + name + "' not found.") diff --git a/src/scaffoldmaker/annotation/ureter_terms.py b/src/scaffoldmaker/annotation/ureter_terms.py new file mode 100644 index 00000000..5c052657 --- /dev/null +++ b/src/scaffoldmaker/annotation/ureter_terms.py @@ -0,0 +1,19 @@ +""" +Common resource for ureter annotation terms. +""" + +# convention: preferred name, preferred id, followed by any other ids and alternative names +ureter_terms = [ + ("ureter", "UBERON:0000056", "ILX:0728080") + ] + +def get_ureter_term(name : str): + """ + Find term by matching name to any identifier held for a term. + Raise exception if name not found. + :return ( preferred name, preferred id ) + """ + for term in ureter_terms: + if name in term: + return (term[0], term[1]) + raise NameError("Body annotation term '" + name + "' not found.") \ No newline at end of file diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_renal_pelvis1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_renal_pelvis1.py new file mode 100644 index 00000000..a379d8d5 --- /dev/null +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_renal_pelvis1.py @@ -0,0 +1,683 @@ +""" +Generates a 3D renal pelvis using tube network mesh. +""" +import math + +from cmlibs.maths.vectorops import mult, cross, add, sub, set_magnitude +from cmlibs.utils.zinc.field import find_or_create_field_coordinates +from cmlibs.zinc.field import Field + +from scaffoldmaker.annotation.annotationgroup import AnnotationGroup +from scaffoldmaker.annotation.kidney_terms import get_kidney_term +from scaffoldmaker.annotation.ureter_terms import get_ureter_term +from scaffoldmaker.meshtypes.meshtype_1d_network_layout1 import MeshType_1d_network_layout1 +from scaffoldmaker.meshtypes.scaffold_base import Scaffold_base +from scaffoldmaker.scaffoldpackage import ScaffoldPackage +from scaffoldmaker.utils.interpolation import sampleCubicHermiteCurves, smoothCubicHermiteDerivativesLine +from scaffoldmaker.utils.networkmesh import NetworkMesh +from scaffoldmaker.utils.tubenetworkmesh import TubeNetworkMeshGenerateData, RenalPelvisTubeNetworkMeshBuilder +from cmlibs.zinc.node import Node + + +class MeshType_1d_renal_pelvis_network_layout1(MeshType_1d_network_layout1): + """ + Defines renal pelvis network layout. + """ + + + @classmethod + def getName(cls): + return "1D Renal Pelvis Network Layout 1" + + @classmethod + def getParameterSetNames(cls): + return ["Default"] + + @classmethod + def getDefaultOptions(cls, parameterSetName="Default"): + options = {} + options["Base parameter set"] = "Human 1" if (parameterSetName == "Default") else parameterSetName + options["Structure"] = ( + "1-2, 2-3.1,3.2-4,3.3-5,3.4-6," + "4.2-7, 4.3-8, 6.2-9, 6.3-10," + "7.2-11, 7.3-12, 8.2-13, 8.3-14, 5.2-15, 5.3-16, 9.2-17, 9.3-18, 10.2-19, 10.3-20," + "11-21-22, 12-23-24, 13-25-26, 14-27-28, 15-29-30, 16-31-32, 17-33-34, 18-35-36, 19-37-38, 20-39-40") + options["Define inner coordinates"] = True + options["Ureter length"] = 3.0 + options["Ureter radius"] = 0.1 + options["Ureter bend angle degrees"] = 45 + options["Major calyx length"] = 0.6 + options["Major calyx radius"] = 0.1 + options["Major calyx angle degrees"] = 170 + options["Middle major calyx length"] = 0.4 + options["Major to bottom/top minor calyx length"] = 0.3 + options["Major to lower/upper minor calyx length"] = 0.3 + options["Bottom/top minor calyx length"] = 0.2 + options["Lower/upper minor calyx length"] = 0.2 + options["Minor calyx radius"] = 0.1 + options["Bottom/top minor calyx bifurcation angle degrees"] = 90 + options["Lower/upper minor calyx bifurcation angle degrees"] = 90 + options["Lower/upper minor calyx bend angle degrees"] = 10 + options["Renal pyramid length"] = 0.5 + options["Renal pyramid width"] = 0.5 + options["Inner proportion default"] = 0.8 + options["Inner proportion ureter"] = 0.7 + return options + + @classmethod + def getOrderedOptionNames(cls): + return [ + "Ureter length", + "Ureter radius", + "Ureter bend angle degrees", + "Major calyx length", + "Major calyx radius", + "Major calyx angle degrees", + "Middle major calyx length", + "Major to bottom/top minor calyx length", + "Major to lower/upper minor calyx length", + "Bottom/top minor calyx length", + "Lower/upper minor calyx length", + "Minor calyx radius", + "Bottom/top minor calyx bifurcation angle degrees", + "Lower/upper minor calyx bifurcation angle degrees", + "Lower/upper minor calyx bend angle degrees", + "Renal pyramid length", + "Renal pyramid width", + "Inner proportion default", + "Inner proportion ureter" + ] + + @classmethod + def checkOptions(cls, options): + dependentChanges = False + for key in [ + "Ureter length", + "Ureter radius", + "Ureter bend angle degrees", + "Major calyx length", + "Major calyx radius", + "Major calyx angle degrees", + "Middle major calyx length", + "Major to bottom/top minor calyx length", + "Major to lower/upper minor calyx length", + "Bottom/top minor calyx length", + "Lower/upper minor calyx length", + "Minor calyx radius", + "Bottom/top minor calyx bifurcation angle degrees", + "Lower/upper minor calyx bifurcation angle degrees", + "Lower/upper minor calyx bend angle degrees", + "Renal pyramid length", + "Renal pyramid width", + "Inner proportion default", + "Inner proportion ureter" + ]: + if options[key] < 0.1: + options[key] = 0.1 # check again + for key in [ + "Inner proportion default", + "Inner proportion ureter" + ]: + if options[key] < 0.1: + options[key] = 0.1 + elif options[key] > 0.9: + options[key] = 0.9 + for key, angleRange in { + "Ureter bend angle degrees": (0.0, 45.0), + "Major calyx angle degrees": (130.0, 170.0), + "Bottom/top minor calyx bifurcation angle degrees": (60.0, 120.0), + "Lower/upper minor calyx bifurcation angle degrees": (60.0, 120.0), + "Lower/upper minor calyx bend angle degrees": (0.0, 10.0) + }.items(): + if options[key] < angleRange[0]: + options[key] = angleRange[0] + elif options[key] > angleRange[1]: + options[key] = angleRange[1] + + return dependentChanges + + @classmethod + def generateBaseMesh(cls, region, options): + """ + Generate the unrefined mesh. + :param region: Zinc region to define model in. Must be empty. + :param options: Dict containing options. See getDefaultOptions(). + :return [] empty list of AnnotationGroup, NetworkMesh + """ + # parameters + structure = options["Structure"] + ureterLength = options["Ureter length"] + ureterRadius = options["Ureter radius"] + ureterBendAngle = options["Ureter bend angle degrees"] + majorCalyxLength = options["Major calyx length"] + majorCalyxRadius = options["Major calyx radius"] + majorCalyxAngle = options["Major calyx angle degrees"] + majorToBottomMinorCalyxLength = options["Major to bottom/top minor calyx length"] + majorToLowerMinorCalyxLength = options["Major to lower/upper minor calyx length"] + middleMajorLength = options["Middle major calyx length"] + bottomMinorCalyxLength = options["Bottom/top minor calyx length"] + minorCalyxRadius = options["Minor calyx radius"] + lowerMinorCalyxLength = options["Lower/upper minor calyx length"] + bottomMinorCalyxAngle = options["Bottom/top minor calyx bifurcation angle degrees"] + lowerMinorCalyxAngle = options["Lower/upper minor calyx bifurcation angle degrees"] + minorCalyxBendAngle = options["Lower/upper minor calyx bend angle degrees"] + pyramidLength = options["Renal pyramid length"] + pyramidWidth = options["Renal pyramid width"] + innerProportionDefault = options["Inner proportion default"] + innerProportionUreter = options["Inner proportion ureter"] + + networkMesh = NetworkMesh(structure) + networkMesh.create1DLayoutMesh(region) + + fieldmodule = region.getFieldmodule() + mesh = fieldmodule.findMeshByDimension(1) + + # set up element annotations + renalPelvisGroup = AnnotationGroup(region, get_kidney_term("renal pelvis")) + ureterGroup = AnnotationGroup(region, get_ureter_term("ureter")) + majorCalyxGroup = AnnotationGroup(region, get_kidney_term("major calyx")) + minorCalyxGroup = AnnotationGroup(region, get_kidney_term("minor calyx")) + renalPyramidGroup = AnnotationGroup(region, get_kidney_term("renal pyramid")) + annotationGroups = [renalPelvisGroup, ureterGroup, majorCalyxGroup, minorCalyxGroup, renalPyramidGroup] + + renalPelvisMeshGroup = renalPelvisGroup.getMeshGroup(mesh) + elementIdentifier = 1 + ureterElementsCount = 2 + meshGroups = [renalPelvisMeshGroup, ureterGroup.getMeshGroup(mesh)] + for e in range(ureterElementsCount): + element = mesh.findElementByIdentifier(elementIdentifier) + for meshGroup in meshGroups: + meshGroup.addElement(element) + elementIdentifier += 1 + + majorCalyxElementsCount = 3 + meshGroups = [renalPelvisMeshGroup, majorCalyxGroup.getMeshGroup(mesh)] + bottomMajor, middleMajor, topMajor = 0, 1, 2 + for e in range(majorCalyxElementsCount): + element = mesh.findElementByIdentifier(elementIdentifier) + if e == middleMajor: + middleMajorCalyxIdentifier = elementIdentifier + for meshGroup in meshGroups: + meshGroup.addElement(element) + elementIdentifier += 1 + + meshGroups = [renalPelvisMeshGroup, minorCalyxGroup.getMeshGroup(mesh)] + bottomMinor, lowerMinor, middleMajor, topMinor, upperMinor = 0, 1, 2, 3, 4 + for calyx in (bottomMinor, lowerMinor, middleMajor, topMinor, upperMinor): + cElementIdentifier = middleMajorCalyxIdentifier if calyx == middleMajor else elementIdentifier + element = mesh.findElementByIdentifier(cElementIdentifier) + for meshGroup in meshGroups: + meshGroup.addElement(element) + elementIdentifier = elementIdentifier if calyx == middleMajor else elementIdentifier + 1 + + renalPyramidMeshGroup = renalPyramidGroup.getMeshGroup(mesh) + minorCalyxElementsCount = 1 + renalPyramidElementsCount = 2 + anterior, posterior = 0, 1 + bottomMinor, lowerMinor, middleMajor, topMinor, upperMinor = 0, 1, 2, 3, 4 + for calyx in (bottomMinor, lowerMinor, middleMajor, topMinor, upperMinor): + for side in (anterior, posterior): + for count, groupCount in [(minorCalyxElementsCount, renalPyramidElementsCount)]: + for e in range(count): + element = mesh.findElementByIdentifier(elementIdentifier) + renalPyramidMeshGroup.addElement(element) + elementIdentifier += 1 + + # set coordinates (outer) + fieldcache = fieldmodule.createFieldcache() + coordinates = find_or_create_field_coordinates(fieldmodule) + # need to ensure inner coordinates are at least defined: + cls.defineInnerCoordinates(region, coordinates, options, networkMesh, innerProportion=innerProportionDefault) + innerCoordinates = find_or_create_field_coordinates(fieldmodule, "inner coordinates") + nodes = fieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES) + + # ureter + nodeIdentifier = 1 + ureterScale = ureterLength / ureterElementsCount + ureterBendAngleRadians = math.radians(ureterBendAngle) + sinUreterAngle = math.sin(ureterBendAngleRadians) + cosUreterAngle = math.cos(ureterBendAngleRadians) + endX = [ureterLength, 0.0, 0.0] + tx = endX[0] - ureterLength * cosUreterAngle + ty = endX[1] - ureterLength * sinUreterAngle + startX = [tx, ty, 0.0] + d1 = [ureterScale, 0.0, 0.0] + d3 = [0.0, 0.0, ureterRadius] + id3 = mult(d3, innerProportionUreter) + td1 = [0.0, ureterScale, 0.0] + sx, sd1 = sampleCubicHermiteCurves([startX, endX], [td1, d1], ureterElementsCount, arcLengthDerivatives=True)[0:2] + sd1 = smoothCubicHermiteDerivativesLine(sx, sd1, fixEndDirection=True) + for e in range(ureterElementsCount + 1): + node = nodes.findNodeByIdentifier(nodeIdentifier) + fieldcache.setNode(node) + sd2 = set_magnitude(cross(d3, sd1[e]), ureterRadius) + sid2 = mult(sd2, innerProportionUreter) + for field, derivatives in ((coordinates, (sd1[e], sd2, d3)), (innerCoordinates, (sd1[e], sid2, id3))): + setNodeFieldParameters(field, fieldcache, sx[e], *derivatives) + nodeIdentifier += 1 + majorCalyxJunctionNodeIdentifier = nodeIdentifier - 1 + majorCalyxStartX = endX + + # major calyx + majorCalyxAngleRadians = math.radians(majorCalyxAngle / 2) + sx = majorCalyxStartX + bottomMajor, middleMajor, topMajor = 0, 1, 2 + majorCalyxNodeIdentifiers, majorCalyxXList = [], [] + for side in (bottomMajor, middleMajor, topMajor): + majorCalyxNodeIdentifiers.append(nodeIdentifier) + calyxLength = middleMajorLength if side == middleMajor else majorCalyxLength + majorCalyxAngle = 0 if side == middleMajor else majorCalyxAngleRadians + cosMajorCalyxAngle = math.cos(majorCalyxAngle) + sinMajorCalyxAngle = math.sin(-majorCalyxAngle if side == bottomMajor else majorCalyxAngle) + majorCalyxEndX = sx[0] + calyxLength * cosMajorCalyxAngle + majorCalyxEndY = sx[1] + calyxLength * sinMajorCalyxAngle + x = [majorCalyxEndX, majorCalyxEndY, 0.0] + majorCalyxXList.append(x) + sd1 = sub(x, sx) + sd2_list, sd3_list, sNodeIdentifiers = [], [], [] + for i in range(2): + isJunctionNode = i == 0 + nodeId = majorCalyxJunctionNodeIdentifier if isJunctionNode else nodeIdentifier + sNodeIdentifiers.append(nodeId) + node = nodes.findNodeByIdentifier(nodeId) + fieldcache.setNode(node) + version = {middleMajor: 3, bottomMajor: 2, topMajor: 4}[side] + sd3 = [0.0, 0.0, majorCalyxRadius] + sid3 = mult(sd3, innerProportionDefault) + sd2 = set_magnitude(cross(sd3, sd1), majorCalyxRadius) + sid2 = mult(sd2, innerProportionDefault) + sd2_list.append(sd2) + sd3_list.append(sd3) + if not isJunctionNode: + for field, derivatives in ( + (coordinates, [x, sd1, sd2, sd3]), + (innerCoordinates, [x, sd1, sid2, sid3]) + ): + for valueLabel, value in zip( + (Node.VALUE_LABEL_VALUE, Node.VALUE_LABEL_D_DS1, Node.VALUE_LABEL_D_DS2, + Node.VALUE_LABEL_D_DS3), + derivatives + ): + field.setNodeParameters(fieldcache, -1, valueLabel, 1, value) + nodeIdentifier += 1 + setNodeFieldVersionDerivatives(coordinates, fieldcache, version, sd1, sd2, sd3) + setNodeFieldVersionDerivatives(innerCoordinates, fieldcache, version, sd1, sid2, sid3) + + # top and bottom major calyx to minor calyx + lowerMinor, bottomMinor = 1, 0 + topMinor, upperMinor = 0, 1 + minorCalyxNodeIdentifiers, minorCalyxXList = [], [] + for calyx in [bottomMajor, topMajor]: + cNodeIdentifier = majorCalyxNodeIdentifiers[calyx] + sides = [bottomMinor, lowerMinor] if calyx == bottomMajor else [topMinor, upperMinor] + signValue = -1.0 if calyx == bottomMajor else 1.0 + sx = majorCalyxXList[calyx] + for side in sides: + calyxLength = majorToBottomMinorCalyxLength if side in (bottomMinor, topMinor) else majorToLowerMinorCalyxLength + x = [sx[0], sx[1] + calyxLength * signValue, sx[2]] if side in (bottomMinor, topMinor) else \ + [sx[0] + calyxLength, sx[1], sx[2]] + if side in (lowerMinor, upperMinor): + theta = math.radians(-minorCalyxBendAngle) if calyx == bottomMajor else math.radians(minorCalyxBendAngle) + rx = sx[0] + (x[0] - sx[0]) * math.cos(theta) - (x[1] - sx[1]) * math.sin(theta) + ry = sx[1] + (x[0] - sx[0]) * math.sin(theta) + (x[1] - sx[1]) * math.cos(theta) + x = [rx, ry, 0.0] + sd1 = sub(x, sx) + minorCalyxNodeIdentifiers.append(nodeIdentifier) + minorCalyxXList.append(x) + minorCalyx_sd2_list, minorCalyx_sd3_list, sNodeIdentifiers = [], [], [] + for i in range(2): + isJunctionNode = i == 0 + nodeId = cNodeIdentifier if isJunctionNode else nodeIdentifier + sNodeIdentifiers.append(nodeId) + node = nodes.findNodeByIdentifier(nodeId) + fieldcache.setNode(node) + + version = 2 if side == 0 else 3 + sd3 = [0.0, 0.0, minorCalyxRadius] + sd2 = [minorCalyxRadius, 0.0, 0.0] if (calyx == bottomMajor and version == 2) \ + else set_magnitude(cross(sd3, sd1), minorCalyxRadius) + sid2 = mult(sd2, innerProportionDefault) + sid3 = mult(sd3, innerProportionDefault) + minorCalyx_sd2_list.append(sd2) + minorCalyx_sd3_list.append(sd3) + if not isJunctionNode: + for field, derivatives in ( + (coordinates, [x, sd1, sd2, sd3]), + (innerCoordinates, [x, sd1, sid2, sid3]) + ): + for valueLabel, value in zip( + (Node.VALUE_LABEL_VALUE, Node.VALUE_LABEL_D_DS1, Node.VALUE_LABEL_D_DS2, + Node.VALUE_LABEL_D_DS3), + derivatives + ): + field.setNodeParameters(fieldcache, -1, valueLabel, 1, value) + nodeIdentifier += 1 + setNodeFieldVersionDerivatives(coordinates, fieldcache, version, sd1, sd2, sd3) + setNodeFieldVersionDerivatives(innerCoordinates, fieldcache, version, sd1, sid2, sid3) + + # minor calyx to renal pyramid connection + anterior, posterior = 0, 1 + bottomMinor, lowerMinor, topMinor, upperMinor, middleMajor = 0, 1, 2, 3, 4 + renalPyramidStartX, renalPyramidNodeIdentifiers = [], [] + connection_sd2_list, connection_sd3_list = [], [] + for calyx in (bottomMinor, lowerMinor, middleMajor, topMinor, upperMinor): + cNodeIdentifier = majorCalyxNodeIdentifiers[1] if calyx == middleMajor else minorCalyxNodeIdentifiers[calyx] + + sx = majorCalyxXList[1] if calyx == middleMajor else minorCalyxXList[calyx] + minorCalyxLength = bottomMinorCalyxLength if calyx in (bottomMinor, topMinor) else lowerMinorCalyxLength + minorCalyxAngle = bottomMinorCalyxAngle if calyx in (bottomMinor, topMinor) else lowerMinorCalyxAngle + minorCalyxHalfAngle = 0.5 * minorCalyxAngle + minorCalyxHalfAngleRadians = math.radians(minorCalyxHalfAngle) + sinMinorCalyxAngle = math.sin(minorCalyxHalfAngleRadians) + cosMinorCalyxAngle = math.cos(minorCalyxHalfAngleRadians) + + renalPyramidNodeIdentifiers.append([]) + renalPyramidStartX.append([]) + connection_sd2_list.append([]) + connection_sd3_list.append([]) + for side in [anterior, posterior]: + if calyx in [bottomMinor, topMinor]: + nx = sx[0] + minorCalyxLength * (-sinMinorCalyxAngle if side == anterior else sinMinorCalyxAngle) + ny = sx[1] + minorCalyxLength * (-cosMinorCalyxAngle if calyx == bottomMinor else cosMinorCalyxAngle) + x = [nx, ny, 0.0] + elif calyx in [lowerMinor, upperMinor]: + signValue = -1 if calyx == lowerMinor else 1 + nx = sx[0] + minorCalyxLength * cosMinorCalyxAngle + ny = sx[1] + minorCalyxLength * math.sin(math.radians(minorCalyxBendAngle)) * signValue + nz = sx[2] + minorCalyxLength * (sinMinorCalyxAngle if side == anterior else -sinMinorCalyxAngle) + x = [nx, ny, nz] + else: + nx = sx[0] + minorCalyxLength * cosMinorCalyxAngle + nz = sx[2] + minorCalyxLength * (sinMinorCalyxAngle if side == anterior else -sinMinorCalyxAngle) + x = [nx, sx[1], nz] + renalPyramidNodeIdentifiers[-1].append(nodeIdentifier) + renalPyramidStartX[-1].append(x) + sd1 = sub(x, sx) + if calyx in [bottomMinor, topMinor]: + sd3 = [0.0, 0.0, minorCalyxRadius] + sd2 = set_magnitude(cross(sd3, sd1), minorCalyxRadius) + elif calyx in [lowerMinor, upperMinor]: + sd2 = set_magnitude(cross([0.0, 0.0, 1.0], sd1), minorCalyxRadius) + sd3 = set_magnitude(cross(sd1, sd2), minorCalyxRadius) + else: + sd2 = [0.0, minorCalyxRadius, 0.0] + sd3 = set_magnitude(cross(sd1, sd2), minorCalyxRadius) + connection_sd2_list[-1].append(sd2) + connection_sd3_list[-1].append(sd3) + + sid2 = mult(sd2, innerProportionDefault) + sid3 = mult(sd3, innerProportionDefault) + sNodeIdentifiers = [] + version = 2 if side == anterior else 3 + for i in range(2): + isJunctionNode = i == 0 + nodeId = cNodeIdentifier if isJunctionNode else nodeIdentifier + sNodeIdentifiers.append(nodeId) + node = nodes.findNodeByIdentifier(nodeId) + fieldcache.setNode(node) + if not isJunctionNode: + for field, derivatives in ( + (coordinates, [x, sd1, sd2, sd3]), + (innerCoordinates, [x, sd1, sid2, sid3]) + ): + for valueLabel, value in zip( + (Node.VALUE_LABEL_VALUE, Node.VALUE_LABEL_D_DS1, Node.VALUE_LABEL_D_DS2, + Node.VALUE_LABEL_D_DS3), + derivatives + ): + field.setNodeParameters(fieldcache, -1, valueLabel, 1, value) + nodeIdentifier += 1 + setNodeFieldVersionDerivatives(coordinates, fieldcache, version, sd1, sd2, sd3) + setNodeFieldVersionDerivatives(innerCoordinates, fieldcache, version, sd1, sid2, sid3) + + # minor calyx to renal pyramids + pyramidElementsCount = 2 + pyramidHalfWidth = 0.5 * pyramidWidth + bottomMinor, lowerMinor, middleMajor, topMinor, upperMinor = 0, 1, 2, 3, 4 + for calyx in [bottomMinor, lowerMinor, middleMajor, topMinor, upperMinor]: + minorCalyxAngle = bottomMinorCalyxAngle if calyx in (bottomMinor, topMinor) else lowerMinorCalyxAngle + minorCalyxHalfAngle = 0.5 * minorCalyxAngle + minorCalyxHalfAngleRadians = math.radians(minorCalyxHalfAngle) + sinMinorCalyxAngle = math.sin(minorCalyxHalfAngleRadians) + cosMinorCalyxAngle = math.cos(minorCalyxHalfAngleRadians) + for side in [anterior, posterior]: + cNodeIdentifier = renalPyramidNodeIdentifiers[calyx][side] + sx = renalPyramidStartX[calyx][side] + if calyx in [bottomMinor, topMinor]: + nx = -sinMinorCalyxAngle if side == anterior else sinMinorCalyxAngle + ny = -cosMinorCalyxAngle if calyx == bottomMinor else cosMinorCalyxAngle + pyramidDirn = [nx, ny, 0.0] + elif calyx in [lowerMinor, upperMinor]: + tx = [1.0, 0.0, 0.0] + theta = math.radians(-minorCalyxBendAngle) if calyx == lowerMinor else math.radians(minorCalyxBendAngle) + rx = tx[0] * math.cos(theta) - tx[1] * math.sin(theta) + ry = tx[0] * math.sin(theta) + tx[1] * math.cos(theta) + pyramidDirn = [rx, ry, 0.0] + else: + pyramidDirn = [1.0, 0.0, 0.0] + xList = [] + for e in range(pyramidElementsCount): + pyramidLengthScale = 0.7 * pyramidLength if e == 0 else pyramidLength + x = add(sx, mult(pyramidDirn, (pyramidLengthScale))) + xList.append(x) + pyramid_sd2_list = [] + pyramid_sd3_list = [] + sNodeIdentifiers = [] + for e in range(pyramidElementsCount): + sNodeIdentifiers.append(nodeIdentifier) + node = nodes.findNodeByIdentifier(nodeIdentifier) + fieldcache.setNode(node) + + sd1 = sub(xList[1], xList[0]) + pyramidWidthScale = pyramidHalfWidth if e == 0 else 0.9 * pyramidHalfWidth + pyramidThickness = 1.1 * minorCalyxRadius if e == 0 else minorCalyxRadius + if calyx in [bottomMinor, topMinor]: + sd3 = [0.0, 0.0, pyramidThickness] + sd2 = set_magnitude(cross(sd3, sd1), pyramidWidthScale) + elif calyx in [lowerMinor, upperMinor]: + sd3 = [0.0, 0.0, pyramidThickness] + sd2 = set_magnitude(cross(sd3, sd1), pyramidWidthScale) + else: + sd2 = [0.0, pyramidWidthScale, 0.0] + sd3 = set_magnitude(cross(sd1, sd2), pyramidThickness) + pyramid_sd2_list.append(sd2) + pyramid_sd3_list.append(sd3) + sid2 = mult(sd2, innerProportionDefault) + sid3 = mult(sd3, innerProportionDefault) + for field, derivatives in ( + (coordinates, [xList[e], sd1, sd2, sd3]), + (innerCoordinates, [xList[e], sd1, sid2, sid3]) + ): + for valueLabel, value in zip( + (Node.VALUE_LABEL_VALUE, Node.VALUE_LABEL_D_DS1, Node.VALUE_LABEL_D_DS2, + Node.VALUE_LABEL_D_DS3), + derivatives + ): + field.setNodeParameters(fieldcache, -1, valueLabel, 1, value) + nodeIdentifier += 1 + pyramid_sd2_list.append(connection_sd2_list[calyx][side]) + pyramid_sd3_list.append(connection_sd3_list[calyx][side]) + for e in range(pyramidElementsCount): + node = nodes.findNodeByIdentifier(sNodeIdentifiers[e]) + fieldcache.setNode(node) + sd12 = sub(pyramid_sd2_list[e + 1], pyramid_sd2_list[e]) + sd13 = sub(pyramid_sd3_list[e + 1], pyramid_sd3_list[e]) + coordinates.setNodeParameters(fieldcache, -1, Node.VALUE_LABEL_D2_DS1DS2, 1, sd12) + coordinates.setNodeParameters(fieldcache, -1, Node.VALUE_LABEL_D2_DS1DS3, 1, sd13) + sid12 = mult(sd12, innerProportionDefault) + sid13 = mult(sd13, innerProportionDefault) + innerCoordinates.setNodeParameters(fieldcache, -1, Node.VALUE_LABEL_D2_DS1DS2, 1, sid12) + innerCoordinates.setNodeParameters(fieldcache, -1, Node.VALUE_LABEL_D2_DS1DS3, 1, sid13) + + return annotationGroups, networkMesh + + @classmethod + def getInteractiveFunctions(cls): + """ + Edit base class list to include only valid functions. + """ + interactiveFunctions = super(MeshType_1d_renal_pelvis_network_layout1, cls).getInteractiveFunctions() + for interactiveFunction in interactiveFunctions: + if interactiveFunction[0] == "Edit structure...": + interactiveFunctions.remove(interactiveFunction) + break + return interactiveFunctions + + +class MeshType_3d_renal_pelvis1(Scaffold_base): + """ + Generates a 3-D renal pelvis. + """ + + @classmethod + def getName(cls): + return "3D Renal Pelvis 1" + + @classmethod + def getParameterSetNames(cls): + return [ + "Default", + "Human 1" + ] + + @classmethod + def getDefaultOptions(cls, parameterSetName='Default'): + options = {} + useParameterSetName = "Human 1" if (parameterSetName == "Default") else parameterSetName + options["Base parameter set"] = useParameterSetName + options["Network layout"] = ScaffoldPackage(MeshType_1d_renal_pelvis_network_layout1) + options["Elements count around"] = 8 + options["Elements count through shell"] = 1 + options["Annotation elements counts around"] = [0] + options["Target element density along longest segment"] = 4.0 + options["Use linear through shell"] = False + options["Use outer trim surfaces"] = True + options["Show trim surfaces"] = False + return options + + @classmethod + def getOrderedOptionNames(cls): + return [ + "Network layout", + "Elements count around", + "Elements count through shell", + "Annotation elements counts around", + "Target element density along longest segment", + "Use linear through shell", + "Use outer trim surfaces", + "Show trim surfaces" + ] + + @classmethod + def getOptionValidScaffoldTypes(cls, optionName): + if optionName == "Network layout": + return [MeshType_1d_renal_pelvis_network_layout1] + return [] + + @classmethod + def getOptionScaffoldTypeParameterSetNames(cls, optionName, scaffoldType): + assert scaffoldType in cls.getOptionValidScaffoldTypes(optionName), \ + cls.__name__ + ".getOptionScaffoldTypeParameterSetNames. " + \ + "Invalid option \"" + optionName + "\" scaffold type " + scaffoldType.getName() + return scaffoldType.getParameterSetNames() # use the defaults from the network layout scaffold + + @classmethod + def getOptionScaffoldPackage(cls, optionName, scaffoldType, parameterSetName=None): + """ + :param parameterSetName: Name of valid parameter set for option Scaffold, or None for default. + :return: ScaffoldPackage. + """ + if parameterSetName: + assert parameterSetName in cls.getOptionScaffoldTypeParameterSetNames(optionName, scaffoldType), \ + "Invalid parameter set " + str(parameterSetName) + " for scaffold " + str(scaffoldType.getName()) + \ + " in option " + str(optionName) + " of scaffold " + cls.getName() + if optionName == "Network layout": + if not parameterSetName: + parameterSetName = "Default" + return ScaffoldPackage(MeshType_1d_renal_pelvis_network_layout1, defaultParameterSetName=parameterSetName) + assert False, cls.__name__ + ".getOptionScaffoldPackage: Option " + optionName + " is not a scaffold" + + @classmethod + def checkOptions(cls, options): + dependentChanges = False + if (options["Network layout"].getScaffoldType() not in + cls.getOptionValidScaffoldTypes("Network layout")): + options["Body network layout"] = ScaffoldPackage(MeshType_1d_renal_pelvis_network_layout1) + + return dependentChanges + + @classmethod + def generateBaseMesh(cls, region, options): + """ + Generate the base hermite-bilinear mesh. See also generateMesh(). + :param region: Zinc region to define model in. Must be empty. + :param options: Dict containing options. See getDefaultOptions(). + :return: list of AnnotationGroup, None + """ + layoutRegion = region.createRegion() + networkLayout = options["Network layout"] + networkLayout.generate(layoutRegion) # ask scaffold to generate to get user-edited parameters + layoutAnnotationGroups = networkLayout.getAnnotationGroups() + networkMesh = networkLayout.getConstructionObject() + + tubeNetworkMeshBuilder = RenalPelvisTubeNetworkMeshBuilder( + networkMesh, + targetElementDensityAlongLongestSegment=options["Target element density along longest segment"], + defaultElementsCountAround=options["Elements count around"], + elementsCountThroughShell=options["Elements count through shell"], + layoutAnnotationGroups=layoutAnnotationGroups, + annotationElementsCountsAround=options["Annotation elements counts around"]) + + tubeNetworkMeshBuilder.build() + generateData = TubeNetworkMeshGenerateData( + region, 3, + isLinearThroughShell=False, + isShowTrimSurfaces=options["Show trim surfaces"]) + tubeNetworkMeshBuilder.generateMesh(generateData) + annotationGroups = generateData.getAnnotationGroups() + + return annotationGroups, None + +def setNodeFieldParameters(field, fieldcache, x, d1, d2, d3, d12=None, d13=None): + """ + Assign node field parameters x, d1, d2, d3 of field. + :param field: Field parameters to assign. + :param fieldcache: Fieldcache with node set. + :param x: Parameters to set for Node.VALUE_LABEL_VALUE. + :param d1: Parameters to set for Node.VALUE_LABEL_D_DS1. + :param d2: Parameters to set for Node.VALUE_LABEL_D_DS2. + :param d3: Parameters to set for Node.VALUE_LABEL_D_DS3. + :param d12: Optional parameters to set for Node.VALUE_LABEL_D2_DS1DS2. + :param d13: Optional parameters to set for Node.VALUE_LABEL_D2_DS1DS3. + :return: + """ + field.setNodeParameters(fieldcache, -1, Node.VALUE_LABEL_VALUE, 1, x) + field.setNodeParameters(fieldcache, -1, Node.VALUE_LABEL_D_DS1, 1, d1) + field.setNodeParameters(fieldcache, -1, Node.VALUE_LABEL_D_DS2, 1, d2) + field.setNodeParameters(fieldcache, -1, Node.VALUE_LABEL_D_DS3, 1, d3) + if d12: + field.setNodeParameters(fieldcache, -1, Node.VALUE_LABEL_D2_DS1DS2, 1, d12) + if d13: + field.setNodeParameters(fieldcache, -1, Node.VALUE_LABEL_D2_DS1DS3, 1, d13) + + +def setNodeFieldVersionDerivatives(field, fieldcache, version, d1, d2, d3, d12=None, d13=None): + """ + Assign node field parameters d1, d2, d3 of field. + :param field: Field to assign parameters of. + :param fieldcache: Fieldcache with node set. + :param version: Version of d1, d2, d3 >= 1. + :param d1: Parameters to set for Node.VALUE_LABEL_D_DS1. + :param d2: Parameters to set for Node.VALUE_LABEL_D_DS2. + :param d3: Parameters to set for Node.VALUE_LABEL_D_DS3. + :param d12: Optional parameters to set for Node.VALUE_LABEL_D2_DS1DS2. + :param d13: Optional parameters to set for Node.VALUE_LABEL_D2_DS1DS3. + :return: + """ + field.setNodeParameters(fieldcache, -1, Node.VALUE_LABEL_D_DS1, version, d1) + field.setNodeParameters(fieldcache, -1, Node.VALUE_LABEL_D_DS2, version, d2) + field.setNodeParameters(fieldcache, -1, Node.VALUE_LABEL_D_DS3, version, d3) + if d12: + field.setNodeParameters(fieldcache, -1, Node.VALUE_LABEL_D2_DS1DS2, version, d12) + if d13: + field.setNodeParameters(fieldcache, -1, Node.VALUE_LABEL_D2_DS1DS3, version, d13) diff --git a/src/scaffoldmaker/scaffolds.py b/src/scaffoldmaker/scaffolds.py index 6c35deb2..7ac6e0a2 100644 --- a/src/scaffoldmaker/scaffolds.py +++ b/src/scaffoldmaker/scaffolds.py @@ -41,6 +41,8 @@ from scaffoldmaker.meshtypes.meshtype_3d_musclefusiform1 import MeshType_3d_musclefusiform1 from scaffoldmaker.meshtypes.meshtype_3d_ostium1 import MeshType_3d_ostium1 from scaffoldmaker.meshtypes.meshtype_3d_ostium2 import MeshType_3d_ostium2 +from scaffoldmaker.meshtypes.meshtype_3d_renal_pelvis1 import MeshType_3d_renal_pelvis1, \ + MeshType_1d_renal_pelvis_network_layout1 from scaffoldmaker.meshtypes.meshtype_3d_smallintestine1 import MeshType_3d_smallintestine1 from scaffoldmaker.meshtypes.meshtype_3d_solidcylinder1 import MeshType_3d_solidcylinder1 from scaffoldmaker.meshtypes.meshtype_3d_solidsphere1 import MeshType_3d_solidsphere1 @@ -102,6 +104,7 @@ def __init__(self): MeshType_3d_musclefusiform1, MeshType_3d_ostium1, MeshType_3d_ostium2, + MeshType_3d_renal_pelvis1, MeshType_3d_smallintestine1, MeshType_3d_solidcylinder1, MeshType_3d_solidsphere1, @@ -120,7 +123,8 @@ def __init__(self): MeshType_3d_wholebody2 ] self._allPrivateScaffoldTypes = [ - MeshType_1d_human_body_network_layout1 + MeshType_1d_human_body_network_layout1, + MeshType_1d_renal_pelvis_network_layout1 ] def findScaffoldTypeByName(self, name): diff --git a/src/scaffoldmaker/utils/tracksurface.py b/src/scaffoldmaker/utils/tracksurface.py index 09bcc78f..17029b0b 100644 --- a/src/scaffoldmaker/utils/tracksurface.py +++ b/src/scaffoldmaker/utils/tracksurface.py @@ -1163,6 +1163,7 @@ def findNearestPositionOnCurve(self, cx, cd1, loop=False, startCurveLocation=Non else: # add out-of-plane slope component if it < 10: + mag_ri = MAX_SLOPE_FACTOR if mag_ri == 0.0 else mag_ri slope_factor = mag_r * mag_r / (mag_ri * mag_ri) else: slope_factor = 1.0 + r_dot_n / mag_r # wrong, but more reliable diff --git a/src/scaffoldmaker/utils/tubenetworkmesh.py b/src/scaffoldmaker/utils/tubenetworkmesh.py index f41e6770..cbe5d3fd 100644 --- a/src/scaffoldmaker/utils/tubenetworkmesh.py +++ b/src/scaffoldmaker/utils/tubenetworkmesh.py @@ -3152,6 +3152,18 @@ def __init__(self, networkMesh: NetworkMesh, targetElementDensityAlongLongestSeg self._layoutInnerCoordinates = None self._useOuterTrimSurfaces = useOuterTrimSurfaces if self._layoutInnerCoordinates else False + def checkSegmentCore(self, networkSegment): + """ + Checks whether a segment should have a core locally, depending on its annotation. Currently, the annotation + requires to include "core" in its name to indicate the segment has a core. + If a segment requires a core locally, it overrides the global value of self._isCore variable. + """ + for layoutAnnotationGroup in self._layoutAnnotationGroups: + if networkSegment.hasLayoutElementsInMeshGroup(layoutAnnotationGroup.getMeshGroup(self._layoutMesh)): + if "core" in layoutAnnotationGroup.getTerm(): + self._isCore = True + return self._isCore + def createSegment(self, networkSegment): pathParametersList = [get_nodeset_path_ordered_field_parameters( self._layoutNodes, self._layoutCoordinates, pathValueLabels, @@ -3173,6 +3185,8 @@ def createSegment(self, networkSegment): elementsCountAround = self._annotationElementsCountsAround[i] break i += 1 + + self.checkSegmentCore(networkSegment) if self._isCore: annotationElementsCountAcrossMinor = [] i = 0 @@ -3253,6 +3267,32 @@ def generateMesh(self, generateData): segment.addSideD3ElementsToMeshGroup(True, dorsalMeshGroup) +class RenalPelvisTubeNetworkMeshBuilder(TubeNetworkMeshBuilder): + """ + Specialization of TubeNetworkMeshBuilder adding annotations for the renal pelvis. + Requires network layout to follow these conventions: + - +y-axis is top, and -y-axis is bottom. + - +d3 direction is anterior, -d3 is posterior. + - naming of major calyxes: top, middle, bottom + - naming of minor calyxes: top, upper, middle, lower, bottom + """ + + def checkSegmentCore(self, networkSegment): + super(RenalPelvisTubeNetworkMeshBuilder, self).checkSegmentCore(networkSegment) + """ + Checks whether a segment should have a core locally, depending on its annotation. + For the renal pelvis scaffold, the annotation term "renal pyramid" indicates the segment has a core. + """ + for layoutAnnotationGroup in self._layoutAnnotationGroups: + if networkSegment.hasLayoutElementsInMeshGroup(layoutAnnotationGroup.getMeshGroup(self._layoutMesh)): + if "renal pyramid" in layoutAnnotationGroup.getTerm(): + self._isCore = True + return self._isCore + + def generateMesh(self, generateData): + super(RenalPelvisTubeNetworkMeshBuilder, self).generateMesh(generateData) + + class TubeEllipseGenerator: """ Generates tube ellipse curves with even-sized elements with specified radius, phase angle, diff --git a/tests/test_renalpelvis.py b/tests/test_renalpelvis.py new file mode 100644 index 00000000..df7e3b3d --- /dev/null +++ b/tests/test_renalpelvis.py @@ -0,0 +1,131 @@ +import math +import unittest + +from cmlibs.utils.zinc.finiteelement import evaluateFieldNodesetRange +from cmlibs.utils.zinc.general import ChangeManager + +from cmlibs.zinc.context import Context +from cmlibs.zinc.element import Element +from cmlibs.zinc.field import Field +from cmlibs.zinc.result import RESULT_OK + +from scaffoldmaker.annotation.annotationgroup import getAnnotationGroupForTerm +from scaffoldmaker.annotation.kidney_terms import get_kidney_term +from scaffoldmaker.annotation.ureter_terms import get_ureter_term +from scaffoldmaker.meshtypes.meshtype_3d_renal_pelvis1 import MeshType_3d_renal_pelvis1 + + +from testutils import assertAlmostEqualList + + +class RenalPelviscaffoldTestCase(unittest.TestCase): + + def test_renalpelvis(self): + """ + Test creation of renal pelvis scaffold. + """ + scaffold = MeshType_3d_renal_pelvis1 + parameterSetNames = scaffold.getParameterSetNames() + self.assertEqual(parameterSetNames, ["Default", "Human 1"]) + options = scaffold.getDefaultOptions("Human 1") + + self.assertEqual(9, len(options)) + self.assertEqual(8, options["Elements count around"]) + self.assertEqual(1, options["Elements count through shell"]) + self.assertEqual([0], options["Annotation elements counts around"]) + self.assertEqual(4.0, options["Target element density along longest segment"]) + self.assertEqual(False, options["Use linear through shell"]) + self.assertEqual(True, options["Use outer trim surfaces"]) + self.assertEqual(False, options["Show trim surfaces"]) + + context = Context("Test") + region = context.getDefaultRegion() + self.assertTrue(region.isValid()) + annotationGroups = scaffold.generateMesh(region, options)[0] + self.assertEqual(7, len(annotationGroups)) + + fieldmodule = region.getFieldmodule() + self.assertEqual(RESULT_OK, fieldmodule.defineAllFaces()) + mesh3d = fieldmodule.findMeshByDimension(3) + self.assertEqual(656, mesh3d.getSize()) + mesh2d = fieldmodule.findMeshByDimension(2) + self.assertEqual(2476, mesh2d.getSize()) + mesh1d = fieldmodule.findMeshByDimension(1) + self.assertEqual(2991, mesh1d.getSize()) + nodes = fieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES) + self.assertEqual(1172, nodes.getSize()) + datapoints = fieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_DATAPOINTS) + self.assertEqual(0, datapoints.getSize()) + + # Check coordinates range, volume + coordinates = fieldmodule.findFieldByName("coordinates").castFiniteElement() + self.assertTrue(coordinates.isValid()) + minimums, maximums = evaluateFieldNodesetRange(coordinates, nodes) + tol = 1.0E-4 + assertAlmostEqualList(self, minimums, [0.7789485582121838, -2.128648920925165, -0.25340433034411736], tol) + assertAlmostEqualList(self, maximums, [4.04142135623731, 1.5517905914526038, 0.2534043303441174], tol) + + with ChangeManager(fieldmodule): + one = fieldmodule.createFieldConstant(1.0) + isExterior = fieldmodule.createFieldIsExterior() + mesh2d = fieldmodule.findMeshByDimension(2) + fieldcache = fieldmodule.createFieldcache() + + volumeField = fieldmodule.createFieldMeshIntegral(one, coordinates, mesh3d) + volumeField.setNumbersOfPoints(3) + result, volume = volumeField.evaluateReal(fieldcache, 1) + self.assertEqual(result, RESULT_OK) + + surfaceAreaField = fieldmodule.createFieldMeshIntegral(isExterior, coordinates, mesh2d) + surfaceAreaField.setNumbersOfPoints(4) + result, surfaceArea = surfaceAreaField.evaluateReal(fieldcache, 1) + self.assertEqual(result, RESULT_OK) + + self.assertAlmostEqual(volume, 0.4091337158319937, delta=tol) + self.assertAlmostEqual(surfaceArea, 13.670578210256929, delta=tol) + + # check some annotation groups: + + expectedSizes3d = { + "core": (320, 0.2200233669825981), + "major calyx": (48, 0.015284024590047216), + "minor calyx": (80, 0.014873238147161206), + "renal pelvis": (176, 0.07598967760825394), + "renal pyramid": (80, 0.017177432136007482), + "shell": (336, 0.18913226215956122), + "ureter": (64, 0.049524462117630445) + } + for name in expectedSizes3d: + term = get_ureter_term(name) if name == "ureter" else get_kidney_term(name) + annotationGroup = getAnnotationGroupForTerm(annotationGroups, term) + size = annotationGroup.getMeshGroup(mesh3d).getSize() + self.assertEqual(expectedSizes3d[name][0], size, name) + volumeMeshGroup = annotationGroup.getMeshGroup(mesh3d) + volumeField = fieldmodule.createFieldMeshIntegral(one, coordinates, volumeMeshGroup) + volumeField.setNumbersOfPoints(4) + fieldcache = fieldmodule.createFieldcache() + result, volume = volumeField.evaluateReal(fieldcache, 1) + self.assertEqual(result, RESULT_OK) + self.assertAlmostEqual(volume, expectedSizes3d[name][1], delta=tol) + + expectedSizes2d = { + "core": (1302, 12.615540404295311), + "major calyx": (208, 1.837248623082175), + "minor calyx": (358, 1.8841596431915135), + "renal pelvis": (734, 7.5061456904809285), + "renal pyramid": (382, 2.1844517878410996), + "shell": (1454, 18.19860793812061), + "ureter": (264, 4.292647846731792) + } + for name in expectedSizes2d: + term = get_ureter_term(name) if name == "ureter" else get_kidney_term(name) + annotationGroup = getAnnotationGroupForTerm(annotationGroups, term) + size = annotationGroup.getMeshGroup(mesh2d).getSize() + self.assertEqual(expectedSizes2d[name][0], size, name) + surfaceMeshGroup = annotationGroup.getMeshGroup(mesh2d) + surfaceAreaField = fieldmodule.createFieldMeshIntegral(one, coordinates, surfaceMeshGroup) + surfaceAreaField.setNumbersOfPoints(4) + fieldcache = fieldmodule.createFieldcache() + result, surfaceArea = surfaceAreaField.evaluateReal(fieldcache, 1) + self.assertEqual(result, RESULT_OK) + self.assertAlmostEqual(surfaceArea, expectedSizes2d[name][1], delta=tol) From cc32301e11e41c45230886d206a672b309d9f7fa Mon Sep 17 00:00:00 2001 From: Chang-Joon Lee Date: Thu, 19 Dec 2024 11:30:12 +1300 Subject: [PATCH 2/2] Fix bending of ureter --- src/scaffoldmaker/meshtypes/meshtype_3d_renal_pelvis1.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_renal_pelvis1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_renal_pelvis1.py index a379d8d5..0e0263b1 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_renal_pelvis1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_renal_pelvis1.py @@ -3,7 +3,7 @@ """ import math -from cmlibs.maths.vectorops import mult, cross, add, sub, set_magnitude +from cmlibs.maths.vectorops import mult, cross, add, sub, set_magnitude, rotate_about_z_axis from cmlibs.utils.zinc.field import find_or_create_field_coordinates from cmlibs.zinc.field import Field @@ -244,7 +244,8 @@ def generateBaseMesh(cls, region, options): d1 = [ureterScale, 0.0, 0.0] d3 = [0.0, 0.0, ureterRadius] id3 = mult(d3, innerProportionUreter) - td1 = [0.0, ureterScale, 0.0] + # td1 = [0.0, ureterScale, 0.0] + td1 = rotate_about_z_axis(d1, 2 * ureterBendAngleRadians) sx, sd1 = sampleCubicHermiteCurves([startX, endX], [td1, d1], ureterElementsCount, arcLengthDerivatives=True)[0:2] sd1 = smoothCubicHermiteDerivativesLine(sx, sd1, fixEndDirection=True) for e in range(ureterElementsCount + 1):