Merge pull request #79 from mlin865/pigColon

Pig colon

src/scaffoldmaker/annotation/colon_terms.py

@@ -4,13 +4,23 @@
 
 # convention: preferred name, preferred id, followed by any other ids and alternative names
 colon_terms = [
-    ( "colon", "UBERON:0001155", "FMA:14543" ),
-    ( "mesenteric zone", None ),
-    ( "non-mesenteric zone", None ),
-    ( "taenia coli", "UBERON:0012419", "FMA:15041" ),
-    ( "tenia libera", None ),
-    ( "tenia mesocolica", None ),
-    ( "tenia omentalis", None )
+    ( "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"),
+    ( "taenia omentalis", "ILX:0739286"),
+    ( "transverse colon", "UBERON:0001157", "FMA:14546", "ILX:0728767")
     ]
 
 def get_colon_term(name : str):

src/scaffoldmaker/annotation/smallintestine_terms.py

@@ -0,0 +1,22 @@
+"""
+Common resource for small intestine annotation terms.
+"""
+
+# convention: preferred name, preferred id, followed by any other ids and alternative names
+smallintestine_terms = [
+    ( "duodenum", "UBERON:0002114", "FMA:7206", "ILX:0726125"),
+    ( "ileum", "UBERON:0002116", "FMA:7208", "ILX:0728151"),
+    ( "jejunum", "UBERON:0002115", "FMA:7207", "ILX:0724224"),
+    ( "small intestine", "UBERON:0002108", "FMA:7200", "ILX:0726770")
+    ]
+
+def get_smallintestine_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 smallintestine_terms:
+        if name in term:
+            return ( term[0], term[1] )
+    raise NameError("Small intestine annotation term '" + name + "' not found.")

src/scaffoldmaker/meshtypes/meshtype_3d_smallintestine1.py

@@ -5,10 +5,11 @@
 """
 
 import copy
+from scaffoldmaker.annotation.annotationgroup import AnnotationGroup
+from scaffoldmaker.annotation.smallintestine_terms import get_smallintestine_term
 from scaffoldmaker.meshtypes.meshtype_1d_path1 import MeshType_1d_path1, extractPathParametersFromRegion
 from scaffoldmaker.meshtypes.scaffold_base import Scaffold_base
 from scaffoldmaker.scaffoldpackage import ScaffoldPackage
-from scaffoldmaker.utils.meshrefinement import MeshRefinement
 from scaffoldmaker.utils import interpolation as interp
 from scaffoldmaker.utils import tubemesh
 from scaffoldmaker.utils.tubemesh import CylindricalSegmentTubeMeshInnerPoints
@@ -245,6 +246,7 @@ def generateBaseMesh(cls, region, options):
             # print(e+1, arcLength)
             length += arcLength
         segmentLength = length / segmentCount
+        elementAlongLength = length / elementsCountAlong
         # print('Length = ', length)
 
         # Sample central path
@@ -257,6 +259,21 @@ def generateBaseMesh(cls, region, options):
         innerRadiusSegmentList, dInnerRadiusSegmentList = interp.sampleParameterAlongLine(lengthList, innerRadiusList,
                                                                                           segmentCount)
 
+        # Create annotation groups for small intestine sections
+        elementsAlongDuodenum = round(duodenumLength / elementAlongLength)
+        elementsAlongJejunum = round(jejunumLength / elementAlongLength)
+        elementsAlongIleum = elementsCountAlong - elementsAlongDuodenum - elementsAlongJejunum
+
+        duodenumGroup = AnnotationGroup(region, get_smallintestine_term("duodenum"))
+        jejunumGroup = AnnotationGroup(region, get_smallintestine_term("jejunum"))
+        ileumGroup = AnnotationGroup(region, get_smallintestine_term("ileum"))
+        annotationGroups = [duodenumGroup, jejunumGroup, ileumGroup]
+        annotationArrayAlong = (['duodenum'] * elementsAlongDuodenum +
+                                ['jejunum'] * elementsAlongJejunum +
+                                ['ileum'] * elementsAlongIleum)
+        annotationArrayAround = [''] * (elementsCountAround)
+        annotationArrayThroughWall = [''] * elementsCountThroughWall
+
         xExtrude = []
         d1Extrude = []
         d2Extrude = []
@@ -328,16 +345,13 @@ def generateBaseMesh(cls, region, options):
             xiList, flatWidthList, length, wallThickness, elementsCountAround,
             elementsCountAlong, elementsCountThroughWall, transitElementList)
 
-        # Create annotation groups
-        annotationGroups = []
-        annotationArray = [''] * (elementsCountAround)
-
         # Create nodes and elements
         nextNodeIdentifier, nextElementIdentifier, annotationGroups = tubemesh.createNodesAndElements(
             region, xList, d1List, d2List, d3List, xFlat, d1Flat, d2Flat, xTexture, d1Texture, d2Texture,
             elementsCountAround, elementsCountAlong, elementsCountThroughWall,
-            annotationGroups, annotationArray, firstNodeIdentifier, firstElementIdentifier,
-            useCubicHermiteThroughWall, useCrossDerivatives, closedProximalEnd=False)
+            annotationGroups, annotationArrayAround, annotationArrayAlong, annotationArrayThroughWall,
+            firstNodeIdentifier, firstElementIdentifier, useCubicHermiteThroughWall, useCrossDerivatives,
+            closedProximalEnd=False)
 
         return annotationGroups
 

src/scaffoldmaker/utils/eftfactory_bicubichermitelinear.py

@@ -268,3 +268,166 @@ def createEftWedgeXi1ZeroOpenTube(self):
         remapEftLocalNodes(eft, 6, ln_map)
         assert eft.validate(), 'eftfactory_tricubichermite.createEftWedgeXi1ZeroOpenTube:  Failed to validate eft'
         return eft
+
+    def createEftTetrahedronXi1One(self, nodeScaleFactorOffset0, nodeScaleFactorOffset1):
+        '''
+        Create a bicubic hermite linear element field for a solid tetrahedron for the apex of cecum,
+        with xi1 and xi3 collapsed on xi2 = 0, and xi3 collapsed on xi1 = 1 and xi2 = 1.
+        Each collapsed node on xi2 = 0 has 3 scale factors giving the cos, sin coefficients of
+        the radial line from global derivatives, plus the arc subtended by the element in radians,
+        so the circumferential direction is rounded.
+        Need to create a new template for each sector around axis giving common
+        nodeScaleFactorOffset values on common faces. Suggestion is to start at 0 and
+        add 10000 for each radial line around axis.
+        :param nodeScaleFactorOffset0: offset of node scale factors at axis on xi1=0
+        :param nodeScaleFactorOffset1: offset of node scale factors at axis on xi1=1
+        :return: Element field template
+        '''
+        # start with full bicubic hermite linear
+        eft = self._mesh.createElementfieldtemplate(self._basis)
+
+        for n in [ 2, 3, 6, 7 ]:
+            eft.setFunctionNumberOfTerms(n * 4 + 4, 0)
+
+        # GRC: allow scale factor identifier for global -1.0 to be prescribed
+        setEftScaleFactorIds(eft, [1], [
+            nodeScaleFactorOffset0 + 1, nodeScaleFactorOffset0 + 2, nodeScaleFactorOffset0 + 3,
+            nodeScaleFactorOffset1 + 1, nodeScaleFactorOffset1 + 2, nodeScaleFactorOffset1 + 3 ] )
+
+        # remap parameters on xi2 = 0 before collapsing nodes
+        remapEftNodeValueLabel(eft, [ 1, 2, 5, 6 ], Node.VALUE_LABEL_D_DS1, [])
+
+        for layer in range(2):
+            soAround = 1
+            ln = layer * 4 + 1
+            # 2 terms for d/dxi2 via general linear map:
+            remapEftNodeValueLabel(eft, [ln], Node.VALUE_LABEL_D_DS2,
+                                   [(Node.VALUE_LABEL_D_DS1, [soAround + 1]),
+                                   (Node.VALUE_LABEL_D_DS2, [soAround + 2])])
+            # 2 terms for cross derivative 1 2 to correct circular apex: cos(theta).phi, -sin(theta).phi
+            remapEftNodeValueLabel(eft, [ln], Node.VALUE_LABEL_D2_DS1DS2,
+                                   [(Node.VALUE_LABEL_D_DS1, [soAround + 2, soAround + 3]),
+                                   (Node.VALUE_LABEL_D_DS2, [1, soAround + 1, soAround + 3])])
+
+            ln = layer * 4 + 2
+            # 2 terms for d/dxi2 via general linear map:
+            remapEftNodeValueLabel(eft, [ln], Node.VALUE_LABEL_D_DS2,
+                                   [(Node.VALUE_LABEL_D_DS1, [soAround + 4]),
+                                   (Node.VALUE_LABEL_D_DS2, [soAround + 5])])
+            # 2 terms for cross derivative 1 2 to correct circular apex: cos(theta).phi, -sin(theta).phi
+            remapEftNodeValueLabel(eft, [ln], Node.VALUE_LABEL_D2_DS1DS2,
+                                   [(Node.VALUE_LABEL_D_DS1, [soAround + 5, soAround + 6]),
+                                    (Node.VALUE_LABEL_D_DS2, [1, soAround + 4, soAround + 6])])
+
+        ln_map = [ 1, 1, 2, 3, 1, 1, 4, 3]
+        remapEftLocalNodes(eft, 4, ln_map)
+
+        assert eft.validate(), 'eftfactory_bicubichermitelinear.createEftTetrahedronXi1One:  Failed to validate eft'
+        return eft
+
+    def createEftTetrahedronXi1Zero(self, nodeScaleFactorOffset0, nodeScaleFactorOffset1):
+        '''
+        Create a bicubic hermite linear element field for a solid tetrahedron for the apex of cecum,
+        with xi1 and xi3 collapsed on xi2 = 0, and xi3 collapsed on xi1 = 0, xi2 = 1.
+        Each collapsed node on xi2 = 0 has 3 scale factors giving the cos, sin coefficients of
+        the radial line from global derivatives, plus the arc subtended by the element in radians,
+        so the circumferential direction is rounded.
+        Need to create a new template for each sector around axis giving common
+        nodeScaleFactorOffset values on common faces. Suggestion is to start at 0 and
+        add 10000 for each radial line around axis.
+        :param nodeScaleFactorOffset0: offset of node scale factors at axis on xi1=0
+        :param nodeScaleFactorOffset1: offset of node scale factors at axis on xi1=1
+        :return: Element field template
+        '''
+        # start with full bicubic hermite linear
+        eft = self._mesh.createElementfieldtemplate(self._basis)
+        for n in [ 2, 3, 6, 7 ]:
+            eft.setFunctionNumberOfTerms(n * 4 + 4, 0)
+
+        # GRC: allow scale factor identifier for global -1.0 to be prescribed
+        setEftScaleFactorIds(eft, [1], [
+            nodeScaleFactorOffset0 + 1, nodeScaleFactorOffset0 + 2, nodeScaleFactorOffset0 + 3,
+            nodeScaleFactorOffset1 + 1, nodeScaleFactorOffset1 + 2, nodeScaleFactorOffset1 + 3 ])
+
+        # remap parameters on xi2 = 0 before collapsing nodes
+        remapEftNodeValueLabel(eft, [ 1, 2, 5, 6 ], Node.VALUE_LABEL_D_DS1, [])
+
+        for layer in range(2):
+            soAround = 1
+            ln = layer * 4 + 1
+            # 2 terms for d/dxi2 via general linear map:
+            remapEftNodeValueLabel(eft, [ln], Node.VALUE_LABEL_D_DS2,
+                                   [(Node.VALUE_LABEL_D_DS1, [soAround + 1]),
+                                   (Node.VALUE_LABEL_D_DS2, [soAround + 2])])
+            # 2 terms for cross derivative 1 2 to correct circular apex: cos(theta).phi, -sin(theta).phi
+            remapEftNodeValueLabel(eft, [ln], Node.VALUE_LABEL_D2_DS1DS2,
+                                   [(Node.VALUE_LABEL_D_DS1, [soAround + 2, soAround + 3]),
+                                   (Node.VALUE_LABEL_D_DS2, [1, soAround + 1, soAround + 3])])
+
+            ln = layer * 4 + 2
+            # 2 terms for d/dxi2 via general linear map:
+            remapEftNodeValueLabel(eft, [ln], Node.VALUE_LABEL_D_DS2,
+                                   [(Node.VALUE_LABEL_D_DS1, [soAround + 4]),
+                                   (Node.VALUE_LABEL_D_DS2, [soAround + 5])])
+            # 2 terms for cross derivative 1 2 to correct circular apex: cos(theta).phi, -sin(theta).phi
+            remapEftNodeValueLabel(eft, [ln], Node.VALUE_LABEL_D2_DS1DS2,
+                                   [(Node.VALUE_LABEL_D_DS1, [soAround + 5, soAround + 6]),
+                                    (Node.VALUE_LABEL_D_DS2, [1, soAround + 4, soAround + 6])])
+
+        ln_map = [ 1, 1, 2, 3, 1, 1, 2, 4]
+        remapEftLocalNodes(eft, 4, ln_map)
+
+        assert eft.validate(), 'eftfactory_bicubichermitelinear.createEftTetrahedronXi1Zero:  Failed to validate eft'
+        return eft
+
+    def createEftPyramidBottomSimple(self, nodeScaleFactorOffset0, nodeScaleFactorOffset1):
+        '''
+        Create a bicubic hermite linear element field for a solid pyramid for elements within
+        a tenia coli joining to the cecal apex, with xi1 and xi3 collapsed on xi2 = 0.
+        Each collapsed node has 3 scale factors giving the cos, sin coefficients of the
+        radial line from global derivatives, plus the arc subtended by the element in radians,
+        so the circumferential direction is rounded. Need to create a new template for each
+        sector around axis giving common nodeScaleFactorOffset values on common faces.
+        Suggestion is to start at 0 and add 10000 for each radial line around axis.
+        :param nodeScaleFactorOffset0: offset of node scale factors at axis on xi1=0
+        :param nodeScaleFactorOffset1: offset of node scale factors at axis on xi1=1
+        :return: Element field template
+        '''
+        # start with full bicubic hermite linear
+        eft = self._mesh.createElementfieldtemplate(self._basis)
+        for n in [ 2, 3, 6, 7 ]:
+            eft.setFunctionNumberOfTerms(n * 4 + 4, 0)
+
+        # GRC: allow scale factor identifier for global -1.0 to be prescribed
+        setEftScaleFactorIds(eft, [1], [
+            nodeScaleFactorOffset0 + 1, nodeScaleFactorOffset0 + 2, nodeScaleFactorOffset0 + 3,
+            nodeScaleFactorOffset1 + 1, nodeScaleFactorOffset1 + 2, nodeScaleFactorOffset1 + 3])
+
+        # remap parameters on xi2 = 0 before collapsing nodes
+        remapEftNodeValueLabel(eft, [ 1, 2, 5, 6 ], Node.VALUE_LABEL_D_DS1, [])
+
+        for layer in range(2):
+            soAround = 1
+            ln = layer * 4 + 1
+            # 2 terms for d/dxi2 via general linear map:
+            remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D_DS2,
+                                   [(Node.VALUE_LABEL_D_DS1, [soAround + 1]), (Node.VALUE_LABEL_D_DS2, [soAround + 2])])
+            # 2 terms for cross derivative 1 2 to correct circular apex: cos(theta).phi, -sin(theta).phi
+            remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS2,
+                                   [(Node.VALUE_LABEL_D_DS1, [soAround + 2, soAround + 3]),
+                                    (Node.VALUE_LABEL_D_DS2, [1, soAround + 1, soAround + 3])])
+
+            ln = layer * 4 + 2
+            # 2 terms for d/dxi2 via general linear map:
+            remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D_DS2,
+                                   [(Node.VALUE_LABEL_D_DS1, [soAround + 4]), (Node.VALUE_LABEL_D_DS2, [soAround + 5])])
+            # 2 terms for cross derivative 1 2 to correct circular apex: cos(theta).phi, -sin(theta).phi
+            remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS2,
+                                   [(Node.VALUE_LABEL_D_DS1, [soAround + 5, soAround + 6]),
+                                    (Node.VALUE_LABEL_D_DS2, [1, soAround + 4, soAround + 6])])
+
+        ln_map = [ 1, 1, 2, 3, 1, 1, 4, 5 ]
+        remapEftLocalNodes(eft, 5, ln_map)
+
+        assert eft.validate(), 'eftfactory_bicubichermitelinear.createEftPyramidBottomSimple:  Failed to validate eft'
+        return eft