Allow for annotation of elements through wall

src/scaffoldmaker/annotation/colon_terms.py

@@ -6,12 +6,16 @@
 colon_terms = [
     ( "ascending colon", "UBERON:0001156", "FMA:14545", "ILX:0734393"),
     ( "descending colon", "UBERON:0001158", "FMA:14547", "ILX:0724444"),
+    ( "caecum", "UBERON:0001153", "FMA:14541", "ILX:0732270"),
     ( "colon", "UBERON:0001155", "FMA:14543", "ILX:0736005"),
+    ( "colonic mucosa", "UBERON:0000317", "FMA:14984", "ILX:0731046"),
     ( "distal colon", "UBERON:0008971", "ILX:0727523"),
     ( "mesenteric zone", None),
     ( "non-mesenteric zone", None),
     ( "proximal colon", "UBERON:0008972", "ILX:0733240"),
+    ( "serosa of colon", "UBERON:0003335", "FMA:14990", "ILX:0736932"),
     ( "spiral colon", "UBERON:0010239", "ILX:0735018"),
+    ( "submucosa of colon", "UBERON:0003331", "FMA:14985", "ILX:0728041"),
     ( "taenia coli", "UBERON:0012419", "FMA:15041", "ILX:0731555"),
     ( "taenia libera", "ILX:0739285"),
     ( "taenia mesocolica", "ILX:0739284"),

src/scaffoldmaker/meshtypes/meshtype_3d_cecum1.py

@@ -6,6 +6,8 @@
 
 import copy
 import math
+from scaffoldmaker.annotation.annotationgroup import AnnotationGroup
+from scaffoldmaker.annotation.colon_terms import get_colon_term
 from scaffoldmaker.meshtypes.meshtype_1d_path1 import MeshType_1d_path1, extractPathParametersFromRegion
 from scaffoldmaker.meshtypes.meshtype_3d_colonsegment1 import ColonSegmentTubeMeshInnerPoints, \
     getFullProfileFromHalfHaustrum, getTeniaColi, createNodesAndElementsTeniaColi
@@ -364,11 +366,16 @@ def generateBaseMesh(cls, region, options):
             segmentLength, wallThickness, cornerInnerRadiusFactor, haustrumInnerRadiusFactorAlongCecum,
             innerRadiusAlongCecum, dInnerRadiusAlongCecum, tcWidthAlongCecum, startPhase)
 
-        annotationArrayAlong = [''] * elementsCountAlongSegment * segmentCount
+        # Create annotation
+        cecumGroup = AnnotationGroup(region, get_colon_term("caecum"))
+        annotationGroups = [cecumGroup]
+        annotationArrayAlong = (['caecum'] * elementsCountAlong)
+
+        annotationArrayThroughWall = [''] * elementsCountThroughWall
 
         for nSegment in range(segmentCount):
             # Make regular segments
-            xInner, d1Inner, d2Inner, transitElementList, segmentAxis, annotationGroups, annotationArrayAround \
+            xInner, d1Inner, d2Inner, transitElementList, segmentAxis, annotationGroupsAround, annotationArrayAround \
                 = colonSegmentTubeMeshInnerPoints.getColonSegmentTubeMeshInnerPoints(nSegment)
 
             # Replace first half of first segment with apex and sample along apex and second half of segment
@@ -416,6 +423,7 @@ def generateBaseMesh(cls, region, options):
                                        zRefList, innerRadiusAlongCecum, closedProximalEnd)
 
         # Create coordinates and derivatives
+        annotationGroups += annotationGroupsAround
         wallThicknessList = [wallThickness] * (elementsCountAlong + 1)
 
         xList, d1List, d2List, d3List, curvatureList = tubemesh.getCoordinatesFromInner(xWarpedList, d1WarpedList,
@@ -466,15 +474,17 @@ def generateBaseMesh(cls, region, options):
             nextNodeIdentifier, nextElementIdentifier, annotationGroups = createNodesAndElementsTeniaColi(
                     region, xCecum, d1Cecum, d2Cecum, d3Cecum, xFlat, d1Flat, d2Flat, xTexture, d1Texture, d2Texture,
                     elementsCountAroundTC, elementsCountAroundHaustrum, elementsCountAlong, elementsCountThroughWall,
-                    tcCount, annotationGroups, annotationArrayAround, annotationArrayAlong, firstNodeIdentifier,
-                    firstElementIdentifier, useCubicHermiteThroughWall, useCrossDerivatives, closedProximalEnd)
+                    tcCount, annotationGroups, annotationArrayAround, annotationArrayAlong, annotationArrayThroughWall,
+                    firstNodeIdentifier, firstElementIdentifier, useCubicHermiteThroughWall, useCrossDerivatives,
+                    closedProximalEnd)
 
         else:
             nextNodeIdentifier, nextElementIdentifier, annotationGroups = tubemesh.createNodesAndElements(
                 region, xCecum, d1Cecum, d2Cecum, d3Cecum, xFlat, d1Flat, d2Flat, xTexture, d1Texture, d2Texture,
                 elementsCountAround, elementsCountAlong, elementsCountThroughWall,
-                annotationGroups, annotationArrayAround, annotationArrayAlong, firstNodeIdentifier,
-                firstElementIdentifier, useCubicHermiteThroughWall, useCrossDerivatives, closedProximalEnd)
+                annotationGroups, annotationArrayAround, annotationArrayAlong, annotationArrayThroughWall,
+                firstNodeIdentifier, firstElementIdentifier, useCubicHermiteThroughWall, useCrossDerivatives,
+                closedProximalEnd)
 
         # Add ostium on track surface between two tenia on the last segment
         elementsAroundTrackSurface = elementsCountAroundHaustrum
@@ -665,14 +675,16 @@ def generateBaseMesh(cls, region, options):
 
         ostiumSettings['Number of elements around ostium'] = len(innerEndPoints_Id)
 
-        nextNodeIdentifier, nextElementIdentifier, (o1_x, o1_d1, o1_d2, o1_d3, o1_NodeId, o1_Positions) = \
-            generateOstiumMesh(region, ostiumSettings, trackSurfaceOstium, centrePosition, axis1,
-                               nextNodeIdentifier, nextElementIdentifier)
-
         fm = region.getFieldmodule()
         mesh = fm.findMeshByDimension(3)
         nodes = fm.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES)
 
+        cecumMeshGroup = cecumGroup.getMeshGroup(mesh)
+
+        nextNodeIdentifier, nextElementIdentifier, (o1_x, o1_d1, o1_d2, o1_d3, o1_NodeId, o1_Positions) = \
+            generateOstiumMesh(region, ostiumSettings, trackSurfaceOstium, centrePosition, axis1,
+                               nextNodeIdentifier, nextElementIdentifier, ostiumMeshGroups= [cecumMeshGroup] )
+
         startProportions = []
         for n in range(len(innerEndPoints_Id)):
             startProportions.append(trackSurfaceOstium.getProportion(o1_Positions[n]))
@@ -681,8 +693,8 @@ def generateBaseMesh(cls, region, options):
             nodes, mesh, nextNodeIdentifier, nextElementIdentifier,
             o1_x, o1_d1, o1_d2, None, o1_NodeId, None,
             endPoints_x, endPoints_d1, endPoints_d2, None, endPoints_Id, endDerivativesMap,
-            elementsCountRadial = 2, tracksurface = trackSurfaceOstium, startProportions = startProportions,
-            endProportions = endProportions)
+            elementsCountRadial = 2, meshGroups = [cecumMeshGroup], tracksurface = trackSurfaceOstium,
+            startProportions = startProportions, endProportions = endProportions)
 
         # Delete elements under annulus mesh
         deleteElementsAndNodesUnderAnnulusMesh(fm, nodes, mesh, deleteElementIdentifier, deleteNodeIdentifier)

src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py

@@ -465,6 +465,8 @@ def generateBaseMesh(cls, region, options):
                                     ['transverse colon'] * elementsAlongInTransverse +
                                     ['descending colon'] * elementsAlongInDistal)
 
+        annotationArrayThroughWall = ([''] * elementsCountThroughWall)
+
         xExtrude = []
         d1Extrude = []
         d2Extrude = []
@@ -536,8 +538,9 @@ def generateBaseMesh(cls, region, options):
             nextNodeIdentifier, nextElementIdentifier, annotationGroups = createNodesAndElementsTeniaColi(
                 region, xList, d1List, d2List, d3List, xFlat, d1Flat, d2Flat, xTexture, d1Texture, d2Texture,
                 elementsCountAroundTC, elementsCountAroundHaustrum, elementsCountAlong, elementsCountThroughWall,
-                tcCount, annotationGroups, annotationArrayAround, annotationArrayAlong, firstNodeIdentifier,
-                firstElementIdentifier, useCubicHermiteThroughWall, useCrossDerivatives, closedProximalEnd)
+                tcCount, annotationGroups, annotationArrayAround, annotationArrayAlong, annotationArrayThroughWall,
+                firstNodeIdentifier, firstElementIdentifier, useCubicHermiteThroughWall, useCrossDerivatives,
+                closedProximalEnd)
 
         else:
             # Create flat and texture coordinates
@@ -549,8 +552,9 @@ def generateBaseMesh(cls, region, options):
             nextNodeIdentifier, nextElementIdentifier, annotationGroups = tubemesh.createNodesAndElements(
                 region, xList, d1List, d2List, d3List, xFlat, d1Flat, d2Flat, xTexture, d1Texture, d2Texture,
                 elementsCountAround, elementsCountAlong, elementsCountThroughWall,
-                annotationGroups, annotationArrayAround, annotationArrayAlong, firstNodeIdentifier,
-                firstElementIdentifier, useCubicHermiteThroughWall, useCrossDerivatives, closedProximalEnd=False)
+                annotationGroups, annotationArrayAround, annotationArrayAlong, annotationArrayThroughWall,
+                firstNodeIdentifier, firstElementIdentifier, useCubicHermiteThroughWall, useCrossDerivatives,
+                closedProximalEnd)
 
         return annotationGroups
 

src/scaffoldmaker/meshtypes/meshtype_3d_colonsegment1.py

@@ -260,14 +260,26 @@ def generateBaseMesh(cls, region, options):
             segmentLength, wallThickness, cornerInnerRadiusFactor, haustrumInnerRadiusFactorAlongSegment,
             radiusAlongSegment, dRadiusAlongSegment, tcWidthAlongSegment, startPhase)
 
+        # Create annotation
+        annotationGroups = []
         annotationArrayAlong = [''] * elementsCountAlongSegment
+        # serosaGroup = AnnotationGroup(region, get_colon_term("serosa of colon"))
+        # submucosaGroup = AnnotationGroup(region, get_colon_term("submucosa of colon"))
+        # mucosaGroup = AnnotationGroup(region, get_colon_term("colonic mucosa"))
+        # annotationGroupsThroughWall = [serosaGroup, submucosaGroup, mucosaGroup]
+        # annotationGroups += annotationGroupsThroughWall
+        # annotationArrayThroughWall = ['colonic mucosa', 'submucosa of colon', 'serosa of colon']
+        annotationArrayThroughWall = [''] * elementsCountThroughWall
 
         # Create inner points
         nSegment = 0
+        closedProximalEnd = False
 
-        xInner, d1Inner, d2Inner, transitElementList, segmentAxis, annotationGroups, annotationArrayAround = \
+        xInner, d1Inner, d2Inner, transitElementList, segmentAxis, annotationGroupsAround, annotationArrayAround = \
             colonSegmentTubeMeshInnerPoints.getColonSegmentTubeMeshInnerPoints(nSegment)
 
+        annotationGroups += annotationGroupsAround
+
         # Project reference point for warping onto central path
         sxRefList, sd1RefList, sd2ProjectedListRef, zRefList = \
             tubemesh.getPlaneProjectionOnCentralPath(xInner, elementsCountAround, elementsCountAlongSegment,
@@ -277,7 +289,7 @@ def generateBaseMesh(cls, region, options):
         xWarpedList, d1WarpedList, d2WarpedList, d3WarpedUnitList = tubemesh.warpSegmentPoints(
             xInner, d1Inner, d2Inner, segmentAxis, sxRefList, sd1RefList, sd2ProjectedListRef,
             elementsCountAround, elementsCountAlongSegment, zRefList, radiusAlongSegment,
-            closedProximalEnd=False)
+            closedProximalEnd)
 
         contractedWallThicknessList = colonSegmentTubeMeshInnerPoints.getContractedWallThicknessList()
 
@@ -288,8 +300,6 @@ def generateBaseMesh(cls, region, options):
 
         relaxedLengthList, xiList = colonSegmentTubeMeshInnerPoints.getRelaxedLengthAndXiList()
 
-        closedProximalEnd = False
-
         if tcThickness > 0:
             tubeTCWidthList = colonSegmentTubeMeshInnerPoints.getTubeTCWidthList()
             xList, d1List, d2List, d3List, annotationGroups, annotationArrayAround = getTeniaColi(
@@ -307,8 +317,9 @@ def generateBaseMesh(cls, region, options):
             nextNodeIdentifier, nextElementIdentifier, annotationGroups = createNodesAndElementsTeniaColi(
                 region, xList, d1List, d2List, d3List, xFlat, d1Flat, d2Flat, xTexture, d1Texture, d2Texture,
                 elementsCountAroundTC, elementsCountAroundHaustrum, elementsCountAlongSegment, elementsCountThroughWall,
-                tcCount, annotationGroups, annotationArrayAround, annotationArrayAlong, firstNodeIdentifier,
-                firstElementIdentifier, useCubicHermiteThroughWall, useCrossDerivatives, closedProximalEnd)
+                tcCount, annotationGroups, annotationArrayAround, annotationArrayAlong, annotationArrayThroughWall,
+                firstNodeIdentifier, firstElementIdentifier, useCubicHermiteThroughWall, useCrossDerivatives,
+                closedProximalEnd)
         else:
             # Create flat and texture coordinates
             xFlat, d1Flat, d2Flat, xTexture, d1Texture, d2Texture = tubemesh.createFlatAndTextureCoordinates(
@@ -319,8 +330,9 @@ def generateBaseMesh(cls, region, options):
             nextNodeIdentifier, nextElementIdentifier, annotationGroups = tubemesh.createNodesAndElements(
                 region, xList, d1List, d2List, d3List, xFlat, d1Flat, d2Flat, xTexture, d1Texture, d2Texture,
                 elementsCountAround, elementsCountAlongSegment, elementsCountThroughWall,
-                annotationGroups, annotationArrayAround, annotationArrayAlong, firstNodeIdentifier,
-                firstElementIdentifier, useCubicHermiteThroughWall, useCrossDerivatives, closedProximalEnd)
+                annotationGroups, annotationArrayAround, annotationArrayAlong, annotationArrayThroughWall,
+                firstNodeIdentifier, firstElementIdentifier, useCubicHermiteThroughWall, useCrossDerivatives,
+                closedProximalEnd)
 
         return annotationGroups
 
@@ -384,9 +396,9 @@ def getColonSegmentTubeMeshInnerPoints(self, nSegment):
                                                nSegment * self._elementsCountAlongSegment:
                                                (nSegment + 1) * self._elementsCountAlongSegment + 1]
 
-        xInner, d1Inner, d2Inner, transitElementList, xiSegment, relaxedLengthSegment, \
-            contractedWallThicknessSegment, segmentAxis, annotationGroups, annotationArray\
-            = getColonSegmentInnerPoints(self._region,
+        xInner, d1Inner, d2Inner, transitElementList, xiSegment, relaxedLengthSegment, contractedWallThicknessSegment, \
+        segmentAxis, annotationGroupsAround, annotationArrayAround = \
+            getColonSegmentInnerPoints(self._region,
                 self._elementsCountAroundTC, self._elementsCountAroundHaustrum, self._elementsCountAlongSegment,
                 self._tcCount, self._segmentLengthEndDerivativeFactor, self._segmentLengthMidDerivativeFactor,
                 self._segmentLength, self._wallThickness,
@@ -407,7 +419,8 @@ def getColonSegmentTubeMeshInnerPoints(self, nSegment):
         contractedWallThickness = contractedWallThicknessSegment[startIdx:self._elementsCountAlongSegment + 1]
         self._contractedWallThicknessList += contractedWallThickness
 
-        return xInner, d1Inner, d2Inner, transitElementList, segmentAxis, annotationGroups, annotationArray
+        return xInner, d1Inner, d2Inner, transitElementList, segmentAxis, \
+               annotationGroupsAround, annotationArrayAround
 
     def getTubeTCWidthList(self):
         return self._tubeTCWidthList
@@ -471,7 +484,7 @@ def getColonSegmentInnerPoints(region, elementsCountAroundTC,
     along colon segment. Assume incompressiblity and a shortened length around will
     result in a thicker wall and vice-versa.
     :return segmentAxis: Axis of segment.
-    :return annotationGroups, annotationArray: annotationArray stores annotation
+    :return annotationGroupsAround, annotationArrayAround: annotationArray stores annotation
     names of elements around.
     """
 
@@ -545,8 +558,8 @@ def getColonSegmentInnerPoints(region, elementsCountAroundTC,
         # Create annotation groups for mouse colon
         mzGroup = AnnotationGroup(region, get_colon_term("mesenteric zone"))
         nonmzGroup = AnnotationGroup(region, get_colon_term("non-mesenteric zone"))
-        annotationGroups = [mzGroup, nonmzGroup]
-        annotationArray = (['mesenteric zone']*int(elementsCountAroundTC*0.5) +
+        annotationGroupsAround = [mzGroup, nonmzGroup]
+        annotationArrayAround = (['mesenteric zone']*int(elementsCountAroundTC*0.5) +
                            ['non-mesenteric zone']*elementsCountAroundHaustrum +
                            ['mesenteric zone']*int(elementsCountAroundTC*0.5))
 
@@ -746,11 +759,11 @@ def getColonSegmentInnerPoints(region, elementsCountAroundTC,
             d2Final = d2Final + d2AlongList
 
         # Create annotation groups
-        annotationGroups = []
-        annotationArray = ['']*(elementsCountAround)
+        annotationGroupsAround = []
+        annotationArrayAround = ['']*(elementsCountAround)
 
     return xFinal, d1Final, d2Final, transitElementList, xiList, relaxedLengthList, contractedWallThicknessList, \
-           segmentAxis, annotationGroups, annotationArray
+           segmentAxis, annotationGroupsAround, annotationArrayAround
 
 def createHalfSetInterHaustralSegment(elementsCountAroundTC, elementsCountAroundHaustrum,
     tcCount, tcWidth, radius, cornerInnerRadiusFactor, sampleElementOut):
@@ -1540,7 +1553,7 @@ def createNodesAndElementsTeniaColi(region,
     xTexture, d1Texture, d2Texture,
     elementsCountAroundTC, elementsCountAroundHaustrum,
     elementsCountAlong, elementsCountThroughWall, tcCount,
-    annotationGroups, annotationArrayAround, annotationArrayAlong,
+    annotationGroups, annotationArrayAround, annotationArrayAlong, annotationArrayThroughWall,
     firstNodeIdentifier, firstElementIdentifier,
     useCubicHermiteThroughWall, useCrossDerivatives, closedProximalEnd):
     """
@@ -1560,6 +1573,7 @@ def createNodesAndElementsTeniaColi(region,
     :param annotationGroups: stores information about annotation groups.
     :param annotationArrayAround: stores annotation names of elements around.
     :param annotationArrayAlong: stores annotation names of elements along.
+    :param annotationArrayThroughWall: stores annotation names of elements through wall.
     :param firstNodeIdentifier, firstElementIdentifier: first node and element
     identifier to use.
     :param useCubicHermiteThroughWall: use linear when false.
@@ -1822,7 +1836,8 @@ def createNodesAndElementsTeniaColi(region,
                 elementIdentifier = elementIdentifier + 1
                 for annotationGroup in annotationGroups:
                     if annotationArrayAround[e1] == annotationGroup._name or \
-                            annotationArrayAlong[0] == annotationGroup._name:
+                            annotationArrayAlong[0] == annotationGroup._name or \
+                            annotationArrayThroughWall[e3] == annotationGroup._name:
                         meshGroup = annotationGroup.getMeshGroup(mesh)
                         meshGroup.addElement(element)
 
@@ -1994,7 +2009,8 @@ def createNodesAndElementsTeniaColi(region,
                 if annotationGroups:
                     for annotationGroup in annotationGroups:
                         if annotationArrayAround[e1] == annotationGroup._name or \
-                                annotationArrayAlong[e2] == annotationGroup._name:
+                                annotationArrayAlong[e2] == annotationGroup._name or \
+                                annotationArrayThroughWall[e3] == annotationGroup._name:
                             meshGroup = annotationGroup.getMeshGroup(mesh)
                             meshGroup.addElement(element)