Merge pull request #30 from rchristie/heart1

Add heart fiducial marker annotation points

scaffoldmaker/meshtypes/meshtype_3d_heart1.py

@@ -12,8 +12,6 @@
 from scaffoldmaker.utils.zinc_utils import *
 from scaffoldmaker.utils.meshrefinement import MeshRefinement
 from opencmiss.zinc.element import Element
-#from opencmiss.zinc.field import Field
-#from opencmiss.zinc.node import Node
 
 class MeshType_3d_heart1(object):
     '''
@@ -100,6 +98,16 @@ def generateBaseMesh(cls, region, options):
         rFibrousRingGroup = AnnotationGroup(region, 'right fibrous ring', FMANumber = 77125, lyphID = 'Lyph ID unknown')
         annotationGroups += [ lFibrousRingGroup, rFibrousRingGroup ]
 
+        # annotation points
+        #dataCoordinates = getOrCreateCoordinateField(fm, 'data_coordinates')
+        dataLabel = getOrCreateLabelField(fm, 'data_label')
+        dataElementXi = getOrCreateElementXiField(fm, 'data_element_xi')
+
+        datapoints = fm.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_DATAPOINTS)
+        datapointTemplateInternal = datapoints.createNodetemplate()
+        datapointTemplateInternal.defineField(dataLabel)
+        datapointTemplateInternal.defineField(dataElementXi)
+
         ##############
         # Create nodes
         ##############
@@ -176,6 +184,7 @@ def generateBaseMesh(cls, region, options):
         eftFibrousRing = bicubichermitelinear.createEftBasic()
 
         # left fibrous ring, starting at crux / collapsed posterior interatrial sulcus
+        cruxElementId = None
         for e in range(-1, elementsCountAroundAtrialFreeWall):
             eft1 = eftFibrousRing
             n1 = e
@@ -270,6 +279,7 @@ def generateBaseMesh(cls, region, options):
                 ln_map = [ 1, 2, 3, 4, 5, 5, 6, 6 ]
                 remapEftLocalNodes(eft1, 6, ln_map)
                 meshGroups += [ lFibrousRingMeshGroup ]
+                cruxElementId = elementIdentifier
             elif e == 0:
                 # general linear map d3 adjacent to collapsed crux/posterior sulcus
                 eft1 = bicubichermitelinear.createEftNoCrossDerivatives()
@@ -349,6 +359,14 @@ def generateBaseMesh(cls, region, options):
             for meshGroup in meshGroups:
                 meshGroup.addElement(element)
 
+        # annotation points
+        cruxElement = mesh.findElementByIdentifier(cruxElementId)
+        datapoint = datapoints.createNode(-1, datapointTemplateInternal)
+        cache.setNode(datapoint)
+        #dataCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, vApexInnerx)
+        dataLabel.assignString(cache, 'crux')
+        dataElementXi.assignMeshLocation(cache, cruxElement, [ 0.0, 0.5, 1.0 ])
+
         fm.endChange()
         return annotationGroups
 

scaffoldmaker/meshtypes/meshtype_3d_heartventricles1.py

@@ -182,6 +182,16 @@ def generateBaseMesh(cls, region, options):
         vSeptumGroup = AnnotationGroup(region, 'interventricular septum', FMANumber = 7133, lyphID = 'Lyph ID unknown')
         annotationGroups = [ lvGroup, rvGroup, vSeptumGroup ]
 
+        # annotation points
+        #dataCoordinates = getOrCreateCoordinateField(fm, 'data_coordinates')
+        dataLabel = getOrCreateLabelField(fm, 'data_label')
+        dataElementXi = getOrCreateElementXiField(fm, 'data_element_xi')
+
+        datapoints = fm.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_DATAPOINTS)
+        datapointTemplateInternal = datapoints.createNodetemplate()
+        datapointTemplateInternal.defineField(dataLabel)
+        datapointTemplateInternal.defineField(dataElementXi)
+
         #################
         # Create nodes
         #################
@@ -954,6 +964,17 @@ def generateBaseMesh(cls, region, options):
                     for meshGroup in meshGroups:
                         meshGroup.addElement(element)
 
+        # apex annotation points
+        element1 = mesh.findElementByIdentifier(1)
+        datapoint = datapoints.createNode(-1, datapointTemplateInternal)
+        cache.setNode(datapoint)
+        dataLabel.assignString(cache, 'apex endo')
+        dataElementXi.assignMeshLocation(cache, element1, [ 0.0, 0.0, 0.0 ])
+        datapoint = datapoints.createNode(-1, datapointTemplateInternal)
+        cache.setNode(datapoint)
+        dataLabel.assignString(cache, 'apex epi')
+        dataElementXi.assignMeshLocation(cache, element1, [ 0.0, 0.0, 1.0 ])
+
         fm.endChange()
         return annotationGroups
 

scaffoldmaker/meshtypes/meshtype_3d_heartventriclesbase1.py

@@ -234,6 +234,16 @@ def generateBaseMesh(cls, region, options):
         lFibrousRingGroup = AnnotationGroup(region, 'left fibrous ring', FMANumber = 77124, lyphID = 'Lyph ID unknown')
         rFibrousRingGroup = AnnotationGroup(region, 'right fibrous ring', FMANumber = 77125, lyphID = 'Lyph ID unknown')
 
+        # annotation points
+        dataCoordinates = getOrCreateCoordinateField(fm, 'data_coordinates')
+        dataLabel = getOrCreateLabelField(fm, 'data_label')
+        #dataElementXi = getOrCreateElementXiField(fm, 'data_element_xi')
+
+        datapoints = fm.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_DATAPOINTS)
+        datapointTemplateExternal = datapoints.createNodetemplate()
+        datapointTemplateExternal.defineField(dataCoordinates)
+        datapointTemplateExternal.defineField(dataLabel)
+
         #################
         # Create nodes
         #################
@@ -315,6 +325,11 @@ def generateBaseMesh(cls, region, options):
         zero = [ 0.0, 0.0, 0.0 ]
         lvOutletOuterd3 = [ None ]*elementsCountAroundOutlet
 
+        datapoint = datapoints.createNode(-1, datapointTemplateExternal)
+        cache.setNode(datapoint)
+        dataCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, lvOutletCentre)
+        dataLabel.assignString(cache, 'aortic valve ctr')
+
         # RV outlet points
         cosRvOutletLeftInclineRadians = math.cos(rvOutletLeftInclineRadians)
         sinRvOutletLeftInclineRadians = math.sin(rvOutletLeftInclineRadians)
@@ -330,6 +345,11 @@ def generateBaseMesh(cls, region, options):
             vector.setMagnitude(axis1, rvOutletOuterRadius), vector.setMagnitude(axis2, rvOutletOuterRadius), elementsCountAroundOutlet)
         rvOutletd2 = [ vOutletElementLength*axis3[c] for c in range(3) ]
 
+        datapoint = datapoints.createNode(-1, datapointTemplateExternal)
+        cache.setNode(datapoint)
+        dataCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, rvOutletCentre)
+        dataLabel.assignString(cache, 'pulmonary valve ctr')
+
         # fix derivative 3 on lv outlet adjacent to rv outlet
         n1 = elementsCountAroundOutlet//2
         lvOutletOuterd3[n1] = interpolateLagrangeHermiteDerivative(lvOutletOuterx[n1], rvOutletOuterx[0], rvOutletd2, 0.0)
@@ -344,6 +364,8 @@ def generateBaseMesh(cls, region, options):
                 aBaseInnerMajorMag, aBaseInnerMinorMag, aMajorAxisRadians,
                 aBaseWallThickness, aBaseSlopeHeight, aBaseSlopeLength, aSeptumThickness,
                 aortaOuterRadius, aBaseFrontInclineRadians, aBaseSideInclineRadians = 0.0, aBaseBackInclineRadians = 0.0)
+        laCentre[2] -= (aBaseSlopeHeight + fibrousRingThickness)
+
         # get d3 from inner-outer difference
         laBaseInnerd3 = []
         laBaseOuterd3 = []
@@ -468,6 +490,16 @@ def generateBaseMesh(cls, region, options):
         pd2 = smoothCubicHermiteDerivativesLine([ rvInnerx[nov], ravInnerx[0][noa]], [ rvInnerd2[nov], d2 ], fixStartDerivative=True, fixEndDirection=True)
         ravInnerd2[0][noa] = pd2[1]
 
+        datapoint = datapoints.createNode(-1, datapointTemplateExternal)
+        cache.setNode(datapoint)
+        dataCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, laCentre)
+        dataLabel.assignString(cache, 'mitral valve ctr')
+
+        datapoint = datapoints.createNode(-1, datapointTemplateExternal)
+        cache.setNode(datapoint)
+        dataCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, [ -laCentre[0], laCentre[1], laCentre[2] ])
+        dataLabel.assignString(cache, 'tricuspid valve ctr')
+
         # set d2 at ra node mid supraventricular crest to be normal to surface; smooth to get final magnitude later
         ravsvcn1 = elementsCountAroundAtria - 3
         mag = baseHeight + baseThickness
@@ -520,8 +552,6 @@ def generateBaseMesh(cls, region, options):
             fixStartDerivative=True, fixEndDirection = True)
         ravInnerd2[0][ravsvcn1] = pd2[1]
 
-        #ravInnerd2[0][ravsvcn1] = [ -d for d in pd1[0] ]
-
         mx = [ (mx[c] + svcid3[c]) for c in range(3) ]
         md2 = svcid2
         nf = 2

scaffoldmaker/utils/zinc_utils.py

@@ -38,6 +38,50 @@ def getOrCreateCoordinateField(fieldmodule, name='coordinates', componentsCount=
     fieldmodule.endChange()
     return coordinates
 
+def getOrCreateElementXiField(fieldmodule, name='label', mesh=None):
+    '''
+    Finds or creates a stored mesh location field for storing locations in the
+    supplied mesh e.g. for defining on annotation points with mesh locations.
+    Raises exception if existing field of name is not stored mesh location type.
+    Note can't currently verify existing field stores locations in the supplied mesh.
+    :param fieldmodule:  Zinc fieldmodule to find or create field in.
+    :param name:  Name of field to find or create.
+    :param mesh:  Mesh to store locations in.
+    '''
+    if mesh is None:
+        mesh = fieldmodule.findMeshByDimension(3)
+    assert mesh.isValid(), 'getOrCreateElementXiField.  Invalid mesh'
+    elementXiField = fieldmodule.findFieldByName(name)
+    if elementXiField.isValid():
+        elementXiField = elementXiField.castStoredMeshLocation()
+        assert elementXiField.isValid(), 'getOrCreateElementXiField.  Existing field \'' + name + '\' is not stored mesh location type'
+        return elementXiField
+    fieldmodule.beginChange()
+    elementXiField = fieldmodule.createFieldStoredMeshLocation(mesh)
+    elementXiField.setName(name)
+    elementXiField.setManaged(True)
+    fieldmodule.endChange()
+    return elementXiField
+
+def getOrCreateLabelField(fieldmodule, name='label'):
+    '''
+    Finds or creates a stored string field for defining labels on nodes, e.g. annotation points.
+    Raises exception if existing field of name is not string-valued.
+    Note can't currently distinguish stored string from constant string fields.
+    :param fieldmodule:  Zinc fieldmodule to find or create field in.
+    :param name:  Name of field to find or create.
+    '''
+    labelField = fieldmodule.findFieldByName(name)
+    if labelField.isValid():
+        assert labelField.getValueType() == Field.VALUE_TYPE_STRING, 'getOrCreateLabelField.  Existing field \'' + name + '\' is not string valued'
+        return labelField
+    fieldmodule.beginChange()
+    labelField = fieldmodule.createFieldStoredString()
+    labelField.setName(name)
+    labelField.setManaged(True)
+    fieldmodule.endChange()
+    return labelField
+
 def getElementNodeIdentifiers(element, eft):
     '''
     Get identifiers of all nodes used by eft in element.