Renal pelvis scaffold

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.")

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.")

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):

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

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,

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)