Add meshtype, terms and unit test for esophagus

src/scaffoldmaker/annotation/esophagus_terms.py

@@ -0,0 +1,31 @@
+"""
+Common resource for esophagus annotation terms.
+"""
+
+# convention: preferred name, preferred id, followed by any other ids and alternative names
+esophagus_terms = [
+    ("abdominal part of esophagus", "UBERON:0035177", "FMA:9397", "ILX:0735274"),
+    ("cervical part of esophagus", "UBERON:0035450", "FMA:9395", "ILX:0734725"),
+    ("distal point of lower esophageal sphincter serosa on the greater curvature of stomach", "None"),
+    ("distal point of lower esophageal sphincter serosa on the lesser curvature of stomach", "None"),
+    ("esophagus", "UBERON:0001043", "FMA:7131", "ILX:0735017"),
+    ("esophagus mucosa", "UBERON:0002469", "FMA:62996", "ILX:0725079"),
+    ("esophagus smooth muscle circular layer", "UBERON:0009960", "FMA:67605", "ILX:0727608"),
+    ("esophagus smooth muscle longitudinal layer", "UBERON:0009961", "FMA:63573", "ILX:0735992"),
+    ("proximodorsal midpoint on serosa of upper esophageal sphincter", "None"),
+    ("proximoventral midpoint on serosa of upper esophageal sphincter", "None"),
+    ("serosa of esophagus", "UBERON:0001975", "FMA:63057", "ILX:0725745"),
+    ("submucosa of esophagus", "UBERON:0001972", "FMA:62997", "ILX:0728662"),
+    ("thoracic part of esophagus", "UBERON:0035216", "FMA:9396", "ILX:0732442"),
+    ]
+
+def get_esophagus_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 esophagus_terms:
+        if name in term:
+            return ( term[0], term[1] )
+    raise NameError("Esophagus annotation term '" + name + "' not found.")

src/scaffoldmaker/annotation/stomach_terms.py

@@ -22,6 +22,8 @@
     ("circular-longitudinal muscle interface of pyloric canal along the lesser curvature", "ILX:0793138"),
     ("circular-longitudinal muscle interface of stomach", "ILX:0793096"),
     ("circular-longitudinal muscle interface of ventral stomach", "ILX:0793097"),
+    ("distal point of lower esophageal sphincter serosa on the greater curvature of stomach", "None"),
+    ("distal point of lower esophageal sphincter serosa on the lesser curvature of stomach", "None"),
     ("dorsal stomach", "ILX:0793086"),
     ("duodenum", "UBERON:0002114", " FMA:7206", "ILX:0726125"),
     ("esophagogastric junction", "UBERON:0007650", "FMA: 9434", "ILX:0733910"),

src/scaffoldmaker/meshtypes/meshtype_3d_stomach1.py

@@ -2765,7 +2765,9 @@ def generateBaseMesh(cls, region, options):
                             "gastroduodenal junction along the lesser curvature on serosa",
                             "body-antrum junction along the greater curvature on serosa",
                             "limiting ridge at the greater curvature on serosa" if limitingRidge else
-                            "fundus-body junction along the greater curvature on serosa"])
+                            "fundus-body junction along the greater curvature on serosa",
+                            "distal point of lower esophageal sphincter serosa on the greater curvature of stomach",
+                            "distal point of lower esophageal sphincter serosa on the lesser curvature of stomach"])
 
         markerInnerElementIdentifiers = [stomachStartElement - elementsCountThroughWall * elementsCountAroundEso,
                                          stomachStartElement - (elementsCountThroughWall - 1) * elementsCountAroundEso -
@@ -2777,7 +2779,9 @@ def generateBaseMesh(cls, region, options):
                                          elementsAroundHalfDuod,
                                          lastDuodenumElementIdentifier - elementsCountThroughWall *
                                          elementsCountAroundDuod * (sum(elementsCountAlongGroups[-3:]) + 1),
-                                         fundusBodyJunctionInnerElementIdentifier]
+                                         fundusBodyJunctionInnerElementIdentifier,
+                                         elementsCountAroundEso * (elementsCountThroughWall - 1) + 1,
+                                         elementsCountAroundEso * elementsCountThroughWall - elementsAroundHalfEso + 1]
 
         elementsCountAroundLayer = [elementsCountAroundEso, elementsCountAroundEso,
                                     elementsCountAroundDuod, elementsCountAroundDuod,
@@ -2787,13 +2791,18 @@ def generateBaseMesh(cls, region, options):
             for n in range(len(markerNames[n3])):
                 markerGroup = findOrCreateAnnotationGroupForTerm(allAnnotationGroups, region,
                                                                  get_stomach_term(markerNames[n3][n]))
-                markerElementIdentifier = \
-                    markerInnerElementIdentifiers[n] + \
-                    (0 if n3 == 0 or elementsCountThroughWall == 1 else elementsCountAroundLayer[n] *
-                                                                        (elementsCountThroughWall - 1))
+                if n < 6:
+                    markerElementIdentifier = \
+                        markerInnerElementIdentifiers[n] + \
+                        (0 if n3 == 0 or elementsCountThroughWall == 1 else elementsCountAroundLayer[n] *
+                                                                            (elementsCountThroughWall - 1))
+                    markerXi = [0.0, 1.0, 0.0 if n3 != len(markerNames) - 1 else 1.0] if n < len(markerNames[n3]) - 1 \
+                        else [0.0, 0.0 if limitingRidge else 1.0, 0.0 if n3 != len(markerNames) - 1 else 1.0]
+                else:
+                    markerElementIdentifier = markerInnerElementIdentifiers[n]
+                    markerXi = [0.0, 0.0, 1.0]
                 markerElement = mesh.findElementByIdentifier(markerElementIdentifier)
-                markerXi = [0.0, 1.0, 0.0 if n3 != len(markerNames) - 1 else 1.0] if n < len(markerNames[n3]) - 1 else \
-                    [0.0, 0.0 if limitingRidge else 1.0, 0.0 if n3 != len(markerNames) - 1 else 1.0]
+
                 cache.setMeshLocation(markerElement, markerXi)
                 markerPoint = markerPoints.createNode(nodeIdentifier, markerTemplateInternal)
                 nodeIdentifier += 1

src/scaffoldmaker/scaffolds.py

@@ -19,6 +19,7 @@
 from scaffoldmaker.meshtypes.meshtype_3d_cecum1 import MeshType_3d_cecum1
 from scaffoldmaker.meshtypes.meshtype_3d_colon1 import MeshType_3d_colon1
 from scaffoldmaker.meshtypes.meshtype_3d_colonsegment1 import MeshType_3d_colonsegment1
+from scaffoldmaker.meshtypes.meshtype_3d_esophagus1 import MeshType_3d_esophagus1
 from scaffoldmaker.meshtypes.meshtype_3d_heart1 import MeshType_3d_heart1
 from scaffoldmaker.meshtypes.meshtype_3d_heart2 import MeshType_3d_heart2
 from scaffoldmaker.meshtypes.meshtype_3d_heartarterialroot1 import MeshType_3d_heartarterialroot1
@@ -68,6 +69,7 @@ def __init__(self):
             MeshType_3d_cecum1,
             MeshType_3d_colon1,
             MeshType_3d_colonsegment1,
+            MeshType_3d_esophagus1,
             MeshType_3d_heart1,
             MeshType_3d_heart2,
             MeshType_3d_heartarterialroot1,

tests/test_esophagus.py

@@ -0,0 +1,162 @@
+import copy
+import unittest
+
+from opencmiss.utils.zinc.finiteelement import evaluateFieldNodesetRange, findNodeWithName
+from opencmiss.utils.zinc.general import ChangeManager
+from opencmiss.zinc.context import Context
+from opencmiss.zinc.element import Element
+from opencmiss.zinc.field import Field
+from opencmiss.zinc.result import RESULT_OK
+from scaffoldmaker.annotation.annotationgroup import getAnnotationGroupForTerm
+from scaffoldmaker.annotation.esophagus_terms import get_esophagus_term
+from scaffoldmaker.meshtypes.meshtype_3d_esophagus1 import MeshType_3d_esophagus1
+from scaffoldmaker.utils.meshrefinement import MeshRefinement
+from scaffoldmaker.utils.zinc_utils import createFaceMeshGroupExteriorOnFace
+
+from testutils import assertAlmostEqualList
+
+
+class EsophagusScaffoldTestCase(unittest.TestCase):
+
+    def test_esophagus1(self):
+        """
+        Test creation of esophagus scaffold.
+        """
+        scaffold = MeshType_3d_esophagus1
+        parameterSetNames = scaffold.getParameterSetNames()
+        self.assertEqual(parameterSetNames, ["Default", "Human 1"])
+        options = scaffold.getDefaultOptions("Human 1")
+        self.assertEqual(15, len(options))
+        self.assertEqual(8, options.get("Number of elements around"))
+        self.assertEqual(12, options.get("Number of elements along"))
+        self.assertEqual(4, options.get("Number of elements through wall"))
+        self.assertEqual(5.0, options.get("Wall thickness"))
+        self.assertEqual(0.35, options.get("Mucosa relative thickness"))
+        self.assertEqual(0.15, options.get("Submucosa relative thickness"))
+        self.assertEqual(0.25, options.get("Circular muscle layer relative thickness"))
+        self.assertEqual(0.25, options.get("Longitudinal muscle layer relative thickness"))
+
+        context = Context("Test")
+        region = context.getDefaultRegion()
+        self.assertTrue(region.isValid())
+        annotationGroups = scaffold.generateBaseMesh(region, options)
+        self.assertEqual(12, len(annotationGroups))
+
+        fieldmodule = region.getFieldmodule()
+        self.assertEqual(RESULT_OK, fieldmodule.defineAllFaces())
+        mesh3d = fieldmodule.findMeshByDimension(3)
+        self.assertEqual(384, mesh3d.getSize())
+        mesh2d = fieldmodule.findMeshByDimension(2)
+        self.assertEqual(1280, mesh2d.getSize())
+        mesh1d = fieldmodule.findMeshByDimension(1)
+        self.assertEqual(1416, mesh1d.getSize())
+        nodes = fieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES)
+        self.assertEqual(524, nodes.getSize())
+        datapoints = fieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_DATAPOINTS)
+        self.assertEqual(0, datapoints.getSize())
+
+        coordinates = fieldmodule.findFieldByName("coordinates").castFiniteElement()
+        self.assertTrue(coordinates.isValid())
+        minimums, maximums = evaluateFieldNodesetRange(coordinates, nodes)
+        assertAlmostEqualList(self, minimums, [-17.57777917463216, -9.777573632674187, -14.999995678650333], 1.0E-6)
+        assertAlmostEqualList(self, maximums, [40.16077850310996, 291.0400575669104, 14.999999571130854], 1.0E-6)
+
+        with ChangeManager(fieldmodule):
+            one = fieldmodule.createFieldConstant(1.0)
+            faceMeshGroup = createFaceMeshGroupExteriorOnFace(fieldmodule, Element.FACE_TYPE_XI3_1)
+            surfaceAreaField = fieldmodule.createFieldMeshIntegral(one, coordinates, faceMeshGroup)
+            surfaceAreaField.setNumbersOfPoints(4)
+            volumeField = fieldmodule.createFieldMeshIntegral(one, coordinates, mesh3d)
+            volumeField.setNumbersOfPoints(3)
+        fieldcache = fieldmodule.createFieldcache()
+        result, surfaceArea = surfaceAreaField.evaluateReal(fieldcache, 1)
+        self.assertEqual(result, RESULT_OK)
+        self.assertAlmostEqual(surfaceArea, 27931.063402890817, delta=1.0E-6)
+        result, volume = volumeField.evaluateReal(fieldcache, 1)
+        self.assertEqual(result, RESULT_OK)
+        self.assertAlmostEqual(volume, 116538.12918074048, delta=1.0E-6)
+
+        # check some annotationGroups:
+        expectedSizes3d = {
+            "abdominal part of esophagus": 64,
+            "cervical part of esophagus": 64,
+            "esophagus": 384,
+            "esophagus mucosa": 96,
+            "esophagus smooth muscle circular layer": 96,
+            "esophagus smooth muscle longitudinal layer": 96,
+            "submucosa of esophagus": 96,
+            "thoracic part of esophagus": 256
+            }
+        for name in expectedSizes3d:
+            group = getAnnotationGroupForTerm(annotationGroups, get_esophagus_term(name))
+            size = group.getMeshGroup(mesh3d).getSize()
+            self.assertEqual(expectedSizes3d[name], size, name)
+
+        # refine 4x4x4 and check result
+        # first remove any surface annotation groups as they are re-added by defineFaceAnnotations
+        removeAnnotationGroups = []
+        for annotationGroup in annotationGroups:
+            if (not annotationGroup.hasMeshGroup(mesh3d)) and \
+                    (annotationGroup.hasMeshGroup(mesh2d) or annotationGroup.hasMeshGroup(mesh1d)):
+                removeAnnotationGroups.append(annotationGroup)
+
+        for annotationGroup in removeAnnotationGroups:
+            annotationGroups.remove(annotationGroup)
+        self.assertEqual(12, len(annotationGroups))
+
+        refineRegion = region.createRegion()
+        refineFieldmodule = refineRegion.getFieldmodule()
+        options['Refine number of elements around'] = 4
+        options['Refine number of elements along'] = 4
+        options['Refine number of elements through wall'] = 4
+        meshrefinement = MeshRefinement(region, refineRegion, annotationGroups)
+        scaffold.refineMesh(meshrefinement, options)
+        annotationGroups = meshrefinement.getAnnotationGroups()
+
+        refineFieldmodule.defineAllFaces()
+        oldAnnotationGroups = copy.copy(annotationGroups)
+        for annotationGroup in annotationGroups:
+            annotationGroup.addSubelements()
+        scaffold.defineFaceAnnotations(refineRegion, options, annotationGroups)
+        for annotation in annotationGroups:
+            if annotation not in oldAnnotationGroups:
+                annotationGroup.addSubelements()
+        self.assertEqual(13, len(annotationGroups))
+#
+        mesh3d = refineFieldmodule.findMeshByDimension(3)
+        self.assertEqual(24576, mesh3d.getSize())
+        mesh2d = refineFieldmodule.findMeshByDimension(2)
+        self.assertEqual(75776, mesh2d.getSize())
+        mesh1d = refineFieldmodule.findMeshByDimension(1)
+        self.assertEqual(77856, mesh1d.getSize())
+        nodes = refineFieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES)
+        self.assertEqual(26660, nodes.getSize())
+        datapoints = refineFieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_DATAPOINTS)
+        self.assertEqual(0, datapoints.getSize())
+
+        # check some refined annotationGroups:
+        for name in expectedSizes3d:
+            group = getAnnotationGroupForTerm(annotationGroups, get_esophagus_term(name))
+            size = group.getMeshGroup(mesh3d).getSize()
+            self.assertEqual(expectedSizes3d[name]*64, size, name)
+
+        # test finding a marker in refined scaffold
+        markerGroup = refineFieldmodule.findFieldByName("marker").castGroup()
+        refinedNodes = refineFieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES)
+        markerNodes = markerGroup.getFieldNodeGroup(refinedNodes).getNodesetGroup()
+        self.assertEqual(4, markerNodes.getSize())
+        markerName = refineFieldmodule.findFieldByName("marker_name")
+        self.assertTrue(markerName.isValid())
+        markerLocation = refineFieldmodule.findFieldByName("marker_location")
+        self.assertTrue(markerLocation.isValid())
+        cache = refineFieldmodule.createFieldcache()
+        node = findNodeWithName(markerNodes, markerName,
+                                "distal point of lower esophageal sphincter serosa on the greater curvature of stomach")
+        self.assertTrue(node.isValid())
+        cache.setNode(node)
+        element, xi = markerLocation.evaluateMeshLocation(cache, 3)
+        self.assertEqual(24125, element.getIdentifier())
+        assertAlmostEqualList(self, xi, [0.0, 1.0, 1.0], 1.0E-10)
+
+if __name__ == "__main__":
+    unittest.main()

tests/test_stomach.py

@@ -53,7 +53,7 @@ def test_stomach1(self):
         region = context.getDefaultRegion()
         self.assertTrue(region.isValid())
         annotationGroups = scaffold.generateBaseMesh(region, options)
-        self.assertEqual(38, len(annotationGroups))
+        self.assertEqual(40, len(annotationGroups))
 
         fieldmodule = region.getFieldmodule()
         self.assertEqual(RESULT_OK, fieldmodule.defineAllFaces())
@@ -64,7 +64,7 @@ def test_stomach1(self):
         mesh1d = fieldmodule.findMeshByDimension(1)
         self.assertEqual(2195, mesh1d.getSize())
         nodes = fieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES)
-        self.assertEqual(830, nodes.getSize())
+        self.assertEqual(832, nodes.getSize())
         datapoints = fieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_DATAPOINTS)
         self.assertEqual(0, datapoints.getSize())
 
@@ -115,7 +115,7 @@ def test_stomach1(self):
 
         for annotationGroup in removeAnnotationGroups:
             annotationGroups.remove(annotationGroup)
-        self.assertEqual(38, len(annotationGroups))
+        self.assertEqual(40, len(annotationGroups))
 
         refineRegion = region.createRegion()
         refineFieldmodule = refineRegion.getFieldmodule()
@@ -133,7 +133,7 @@ def test_stomach1(self):
         for annotation in annotationGroups:
             if annotation not in oldAnnotationGroups:
                 annotationGroup.addSubelements()
-        self.assertEqual(72, len(annotationGroups))
+        self.assertEqual(74, len(annotationGroups))
 #
         mesh3d = refineFieldmodule.findMeshByDimension(3)
         self.assertEqual(37248, mesh3d.getSize())
@@ -142,7 +142,7 @@ def test_stomach1(self):
         mesh1d = refineFieldmodule.findMeshByDimension(1)
         self.assertEqual(118796, mesh1d.getSize())
         nodes = refineFieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES)
-        self.assertEqual(40814, nodes.getSize())
+        self.assertEqual(40816, nodes.getSize())
         datapoints = refineFieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_DATAPOINTS)
         self.assertEqual(0, datapoints.getSize())
 
@@ -156,7 +156,7 @@ def test_stomach1(self):
         markerGroup = refineFieldmodule.findFieldByName("marker").castGroup()
         refinedNodes = refineFieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES)
         markerNodes = markerGroup.getFieldNodeGroup(refinedNodes).getNodesetGroup()
-        self.assertEqual(18, markerNodes.getSize())
+        self.assertEqual(20, markerNodes.getSize())
         markerName = refineFieldmodule.findFieldByName("marker_name")
         self.assertTrue(markerName.isValid())
         markerLocation = refineFieldmodule.findFieldByName("marker_location")