From 3edf2a0aa0b6034991f6ade8a77a65ec7d659f32 Mon Sep 17 00:00:00 2001 From: Richard Christie Date: Fri, 10 Apr 2020 16:52:17 +1200 Subject: [PATCH 01/15] Fix atria annotation group names --- src/scaffoldmaker/meshtypes/meshtype_3d_heartatria1.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_heartatria1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_heartatria1.py index 1e06cab9..157c7262 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_heartatria1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_heartatria1.py @@ -3218,9 +3218,9 @@ def refineMesh(cls, meshrefinement, options): sourceFm = meshrefinement._sourceFm annotationGroups = meshrefinement._sourceAnnotationGroups - laGroup = findAnnotationGroupByName(annotationGroups, 'left atrium') + laGroup = findAnnotationGroupByName(annotationGroups, 'left atrium myocardium') laElementGroupField = laGroup.getFieldElementGroup(meshrefinement._sourceMesh) - raGroup = findAnnotationGroupByName(annotationGroups, 'right atrium') + raGroup = findAnnotationGroupByName(annotationGroups, 'right atrium myocardium') raElementGroupField = raGroup.getFieldElementGroup(meshrefinement._sourceMesh) aSeptumGroup = findAnnotationGroupByName(annotationGroups, 'interatrial septum') aSeptumElementGroupField = aSeptumGroup.getFieldElementGroup(meshrefinement._sourceMesh) From 80d485ef53a90a75efdca47d6332df333301902c Mon Sep 17 00:00:00 2001 From: Richard Christie Date: Thu, 16 Apr 2020 20:32:20 +1200 Subject: [PATCH 02/15] Transfer embedded marker points to refined mesh --- .../meshtypes/meshtype_3d_heart1.py | 8 +-- .../meshtypes/meshtype_3d_heartventricles1.py | 8 +-- src/scaffoldmaker/utils/exportvtk.py | 3 +- src/scaffoldmaker/utils/meshrefinement.py | 66 ++++++++++++++++++- 4 files changed, 74 insertions(+), 11 deletions(-) diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_heart1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_heart1.py index f01dbf7e..7acaae84 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_heart1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_heart1.py @@ -136,14 +136,14 @@ def generateBaseMesh(cls, region, options): # annotation fiducial points markerGroup = findOrCreateFieldGroup(fm, "marker") - markerCoordinates = findOrCreateFieldCoordinates(fm, "marker_coordinates") + #markerCoordinates = findOrCreateFieldCoordinates(fm, "marker_coordinates") markerName = findOrCreateFieldStoredString(fm, name="marker_name") markerLocation = findOrCreateFieldStoredMeshLocation(fm, mesh, name="marker_location") nodes = fm.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES) markerPoints = findOrCreateFieldNodeGroup(markerGroup, nodes).getNodesetGroup() markerTemplateInternal = nodes.createNodetemplate() - markerTemplateInternal.defineField(markerCoordinates) + #markerTemplateInternal.defineField(markerCoordinates) markerTemplateInternal.defineField(markerName) markerTemplateInternal.defineField(markerLocation) @@ -401,13 +401,13 @@ def generateBaseMesh(cls, region, options): # annotation fiducial points cruxElement = mesh.findElementByIdentifier(cruxElementId) - cruxXi = [ 0.0, 0.5, 1.0 ] + cruxXi = [ 0.5, 0.5, 1.0 ] cache.setMeshLocation(cruxElement, cruxXi) result, cruxCoordinates = coordinates.evaluateReal(cache, 3) markerPoint = markerPoints.createNode(nodeIdentifier, markerTemplateInternal) nodeIdentifier += 1 cache.setNode(markerPoint) - markerCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, cruxCoordinates) + #markerCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, cruxCoordinates) markerName.assignString(cache, 'crux') markerLocation.assignMeshLocation(cache, cruxElement, cruxXi) diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_heartventricles1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_heartventricles1.py index d3864d87..114a00bf 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_heartventricles1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_heartventricles1.py @@ -256,14 +256,14 @@ def generateBaseMesh(cls, region, options): # annotation fiducial points markerGroup = findOrCreateFieldGroup(fm, "marker") - markerCoordinates = findOrCreateFieldCoordinates(fm, "marker_coordinates") + #markerCoordinates = findOrCreateFieldCoordinates(fm, "marker_coordinates") markerName = findOrCreateFieldStoredString(fm, name="marker_name") markerLocation = findOrCreateFieldStoredMeshLocation(fm, mesh, name="marker_location") nodes = fm.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES) markerPoints = findOrCreateFieldNodeGroup(markerGroup, nodes).getNodesetGroup() markerTemplateInternal = nodes.createNodetemplate() - markerTemplateInternal.defineField(markerCoordinates) + #markerTemplateInternal.defineField(markerCoordinates) markerTemplateInternal.defineField(markerName) markerTemplateInternal.defineField(markerLocation) @@ -1073,13 +1073,13 @@ def generateBaseMesh(cls, region, options): markerPoint = markerPoints.createNode(nodeIdentifier, markerTemplateInternal) nodeIdentifier += 1 cache.setNode(markerPoint) - markerCoordinates.assignReal(cache, lvApexInnerx) + #markerCoordinates.assignReal(cache, lvApexInnerx) markerName.assignString(cache, 'apex endo') markerLocation.assignMeshLocation(cache, element1, [ 0.0, 0.0, 0.0 ]) markerPoint = markerPoints.createNode(nodeIdentifier, markerTemplateInternal) nodeIdentifier += 1 cache.setNode(markerPoint) - markerCoordinates.assignReal(cache, lvApexOuterx) + #markerCoordinates.assignReal(cache, lvApexOuterx) markerName.assignString(cache, 'apex epi') markerLocation.assignMeshLocation(cache, element1, [ 0.0, 0.0, 1.0 ]) diff --git a/src/scaffoldmaker/utils/exportvtk.py b/src/scaffoldmaker/utils/exportvtk.py index 03b5d058..a72d80f1 100644 --- a/src/scaffoldmaker/utils/exportvtk.py +++ b/src/scaffoldmaker/utils/exportvtk.py @@ -5,6 +5,7 @@ import io from sys import version_info from opencmiss.zinc.field import Field +from opencmiss.utils.zinc.field import findOrCreateFieldCoordinates from opencmiss.utils.zinc.finiteelement import getElementNodeIdentifiersBasisOrder @@ -26,7 +27,7 @@ def __init__(self, region, description, annotationGroups = None): self._mesh = self._fieldmodule.findMeshByDimension(dimension) if self._mesh.getSize() > 0: break - self._coordinates = self._fieldmodule.findFieldByName('coordinates') + self._coordinates = findOrCreateFieldCoordinates(self._fieldmodule) self._description = description self._annotationGroups = annotationGroups if annotationGroups else [] diff --git a/src/scaffoldmaker/utils/meshrefinement.py b/src/scaffoldmaker/utils/meshrefinement.py index a35728cf..38383994 100644 --- a/src/scaffoldmaker/utils/meshrefinement.py +++ b/src/scaffoldmaker/utils/meshrefinement.py @@ -3,7 +3,8 @@ ''' from __future__ import division import math -from opencmiss.utils.zinc.field import findOrCreateFieldCoordinates +from opencmiss.utils.zinc.field import findOrCreateFieldCoordinates, findOrCreateFieldGroup, findOrCreateFieldNodeGroup, \ + findOrCreateFieldStoredMeshLocation, findOrCreateFieldStoredString from opencmiss.zinc.element import Element, Elementbasis from opencmiss.zinc.field import Field from opencmiss.zinc.node import Node @@ -67,7 +68,7 @@ def __init__(self, sourceRegion, targetRegion, sourceAnnotationGroups = []): self._nodeIdentifier = 1 self._elementIdentifier = 1 - + # prepare annotation group map self._sourceAnnotationGroups = sourceAnnotationGroups self._annotationGroups = [] self._sourceAndTargetMeshGroups = [] @@ -78,14 +79,45 @@ def __init__(self, sourceRegion, targetRegion, sourceAnnotationGroups = []): targetMeshGroup = targetAnnotationGroup.getMeshGroup(self._targetMesh) self._annotationGroups.append(targetAnnotationGroup) self._sourceAndTargetMeshGroups.append( ( sourceMeshGroup, targetMeshGroup) ) + # prepare element -> marker point list map + self.elementMarkerMap = {} + sourceMarkerGroup = findOrCreateFieldGroup(self._sourceFm, "marker") + sourceMarkerName = findOrCreateFieldStoredString(self._sourceFm, name="marker_name") + sourceMarkerLocation = findOrCreateFieldStoredMeshLocation(self._sourceFm, self._sourceMesh, name="marker_location") + sourceMarkerNodes = findOrCreateFieldNodeGroup(sourceMarkerGroup, sourceNodes).getNodesetGroup() + nodeIter = sourceMarkerNodes.createNodeiterator() + node = nodeIter.next() + while node.isValid(): + self._sourceCache.setNode(node) + element, xi = sourceMarkerLocation.evaluateMeshLocation(self._sourceCache, self._sourceMesh.getDimension()) + if element.isValid(): + elementIdentifier = element.getIdentifier() + markerName = sourceMarkerName.evaluateString(self._sourceCache) + markerList = self.elementMarkerMap.get(elementIdentifier) + if not markerList: + markerList = [] + self.elementMarkerMap[elementIdentifier] = markerList + markerList.append( ( markerName, xi ) ) + node = nodeIter.next() + if self.elementMarkerMap: + self._targetMarkerGroup = findOrCreateFieldGroup(self._targetFm, "marker") + self._targetMarkerName = findOrCreateFieldStoredString(self._targetFm, name="marker_name") + self._targetMarkerLocation = findOrCreateFieldStoredMeshLocation(self._targetFm, self._targetMesh, name="marker_location") + self._targetMarkerNodes = findOrCreateFieldNodeGroup(self._targetMarkerGroup, self._targetNodes).getNodesetGroup() + self._targetMarkerTemplate = self._targetMarkerNodes.createNodetemplate() + self._targetMarkerTemplate.defineField(self._targetMarkerName) + self._targetMarkerTemplate.defineField(self._targetMarkerLocation) + def __del__(self): self._sourceFm.endChange() self._targetFm.endChange() + def getAnnotationGroups(self): return self._annotationGroups + def refineElementCubeStandard3d(self, sourceElement, numberInXi1, numberInXi2, numberInXi3, addNewNodesToOctree=True, shareNodeIds=None, shareNodeCoordinates=None): ''' @@ -146,6 +178,7 @@ def refineElementCubeStandard3d(self, sourceElement, numberInXi1, numberInXi2, n nids.append(nodeId) nx.append(x) # create elements + startElementIdentifier = self._elementIdentifier for k in range(numberInXi3): ok = (numberInXi2 + 1)*(numberInXi1 + 1) for j in range(numberInXi2): @@ -162,8 +195,37 @@ def refineElementCubeStandard3d(self, sourceElement, numberInXi1, numberInXi2, n for meshGroup in meshGroups: meshGroup.addElement(element) + + # re-map any markers embedded in the source element + if self.elementMarkerMap: + markerList = self.elementMarkerMap.get(sourceElement.getIdentifier()) + if markerList: + numberInXi = [ numberInXi1, numberInXi2, numberInXi3 ] + elementOffset = [ 1, numberInXi1, numberInXi1*numberInXi2 ] + targetXi = [ 0.0 ]*3 + for marker in markerList: + markerName, sourceXi = marker + node = self._targetMarkerNodes.createNode(self._nodeIdentifier, self._targetMarkerTemplate) + self._targetCache.setNode(node) + self._targetMarkerName.assignString(self._targetCache, markerName) + # determine which sub-element, targetXi that sourceXi maps to + targetElementIdentifier = startElementIdentifier + for i in range(3): + targetXi[i] = sourceXi[i]*numberInXi[i] + el = int(targetXi[i]) + if el < numberInXi[i]: + targetXi[i] -= el + else: + el = numberInXi[i] - 1 + targetXi[i] = 1.0 + targetElementIdentifier += el*elementOffset[i] + targetElement = self._targetMesh.findElementByIdentifier(targetElementIdentifier) + result = self._targetMarkerLocation.assignMeshLocation(self._targetCache, targetElement, targetXi) + self._nodeIdentifier += 1 + return nids, nx + def refineAllElementsCubeStandard3d(self, numberInXi1, numberInXi2, numberInXi3): element = self._sourceElementiterator.next() while element.isValid(): From e424b87455e956cacadc9e14a448bf4a96d7bc84 Mon Sep 17 00:00:00 2001 From: Richard Christie Date: Mon, 20 Apr 2020 17:42:08 +1200 Subject: [PATCH 03/15] Make shared annotation terms for heart, colon, bladder Make heart fiducial markers match what is being digitised. --- .../annotation/annotationgroup.py | 71 +++++++++++---- src/scaffoldmaker/annotation/bladder_terms.py | 21 +++++ src/scaffoldmaker/annotation/colon_terms.py | 25 ++++++ src/scaffoldmaker/annotation/heart_terms.py | 63 +++++++++++++ .../meshtypes/meshtype_3d_bladder1.py | 5 +- .../meshtypes/meshtype_3d_colonsegment1.py | 15 ++-- .../meshtypes/meshtype_3d_heart1.py | 14 ++- .../meshtypes/meshtype_3d_heart2.py | 2 +- .../meshtype_3d_heartarterialroot1.py | 19 ++-- .../meshtypes/meshtype_3d_heartatria1.py | 88 +++++++------------ .../meshtypes/meshtype_3d_heartatria2.py | 21 ++--- .../meshtypes/meshtype_3d_heartventricles1.py | 31 +++---- .../meshtypes/meshtype_3d_heartventricles2.py | 7 +- .../meshtype_3d_heartventriclesbase1.py | 46 +++++++--- .../meshtype_3d_heartventriclesbase2.py | 18 ++-- src/scaffoldmaker/utils/meshrefinement.py | 3 +- 16 files changed, 304 insertions(+), 145 deletions(-) create mode 100644 src/scaffoldmaker/annotation/bladder_terms.py create mode 100644 src/scaffoldmaker/annotation/colon_terms.py create mode 100644 src/scaffoldmaker/annotation/heart_terms.py diff --git a/src/scaffoldmaker/annotation/annotationgroup.py b/src/scaffoldmaker/annotation/annotationgroup.py index 96b08fd2..96b83797 100644 --- a/src/scaffoldmaker/annotation/annotationgroup.py +++ b/src/scaffoldmaker/annotation/annotationgroup.py @@ -9,35 +9,41 @@ class AnnotationGroup(object): Describes subdomains of a scaffold with attached names and terms. ''' - def __init__(self, region, name, FMANumber, lyphID): + def __init__(self, region, term): ''' :param region: The Zinc region the AnnotationGroup is to be made for. - :param name: The name of an annotation group e.g. common medical term. - :param FMANumber: The FMA Number of the group. - :param lyphID: The Apinatomy Lyph ID for the group. + :param term: Identifier for anatomical term, currently a tuple of name, id. + e.g. ("heart", "FMA:7088") ''' - self._name = name - self._FMANumber = FMANumber - self._lyphID = lyphID + self._name = term[0] + self._id = term[1] fm = region.getFieldmodule() - field = fm.findFieldByName(name) + field = fm.findFieldByName(self._name) if field.isValid(): self._group = field.castGroup() - assert self._group.isValid(), 'AnnotationGroup found existing non-group field called ' + name + assert self._group.isValid(), 'AnnotationGroup found existing non-group field called ' + self._name else: # assume client is calling between fm.begin/endChange() self._group = fm.createFieldGroup() - self._group.setName(name) + self._group.setName(self._name) self._group.setManaged(True) def getName(self): return self._name + def getId(self): + return self._id + def getFMANumber(self): - return self._FMANumber + """ + :return: Integer FMA number or None. + """ + if self._id and (self.id[:4] == "FMA:"): + return int(self._id[4:]) + return None - def getLyphID(self): - return self._lyphID + def getTerm(self): + return ( self._name, self._id ) def getGroup(self): return self._group @@ -92,13 +98,48 @@ def addSubelements(self): meshGroup.addElementsConditional(elementGroup) # use FieldElementGroup as conditional field -def findAnnotationGroupByName(annotationGroups, name): +def findAnnotationGroupByName(annotationGroups: list, name: str): ''' + Find existing annotation group for name. :param annotationGroups: list(AnnotationGroup) :param name: Name of group. - :return: AnnotationGroup or None. + :return: AnnotationGroup or None if not found. ''' for annotationGroup in annotationGroups: if annotationGroup._name == name: return annotationGroup return None + + +def findOrCreateAnnotationGroupForTerm(annotationGroups: list, region, term) -> AnnotationGroup: + ''' + Find existing annotation group for term, or create it for region if not gound. + If annotation group created here, append it to annotationGroups. + :param annotationGroups: list(AnnotationGroup) + :param region: Zinc region to create group for. + :param term: Identifier for anatomical term, currently a tuple of name, id. + :return: AnnotationGroup. + ''' + name = term[0] + annotationGroup = findAnnotationGroupByName(annotationGroups, name) + if annotationGroup: + assert annotationGroup._id == term[1], "Annotation group '" + name + "' id '" + term[1] + "' does not match existing id '" + annotationGroup._id + "'" + else: + annotationGroup = AnnotationGroup(region, term) + annotationGroups.append(annotationGroup) + return annotationGroup + + +def getAnnotationGroupForTerm(annotationGroups: list, term) -> AnnotationGroup: + ''' + Get existing annotation group for term. Raise exception if not found. + :param annotationGroups: list(AnnotationGroup) + :param term: Identifier for anatomical term, currently a tuple of name, id. + :return: AnnotationGroup. + ''' + name = term[0] + annotationGroup = findAnnotationGroupByName(annotationGroups, name) + if annotationGroup: + assert annotationGroup._id == term[1], "Annotation group '" + name + "' id '" + term[1] + "' does not match existing id '" + annotationGroup._id + "'" + return annotationGroup + raise NameError("Annotation group '" + name + "' not found.") diff --git a/src/scaffoldmaker/annotation/bladder_terms.py b/src/scaffoldmaker/annotation/bladder_terms.py new file mode 100644 index 00000000..0e7d6bb9 --- /dev/null +++ b/src/scaffoldmaker/annotation/bladder_terms.py @@ -0,0 +1,21 @@ +""" +Common resource for bladder annotation terms. +""" + +# convention: preferred name, preferred id, followed by any other ids and alternative names +bladder_terms = [ + ( "urinary bladder", "FMA:15900", "UBERON:0001255" ), + ( "neck of urinary bladder", "FMA:15912", "UBERON:0001258"), + ( "body of urinary bladder", None) # needs to be replaced with an actual term e.g. urinary bladder (which includes neck) + ] + +def get_bladder_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 bladder_terms: + if name in term: + return ( term[0], term[1] ) + raise NameError("Bladder annotation term '" + name + "' not found.") diff --git a/src/scaffoldmaker/annotation/colon_terms.py b/src/scaffoldmaker/annotation/colon_terms.py new file mode 100644 index 00000000..f59e6e3e --- /dev/null +++ b/src/scaffoldmaker/annotation/colon_terms.py @@ -0,0 +1,25 @@ +""" +Common resource for colon annotation terms. +""" + +# 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 ) + ] + +def get_colon_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 colon_terms: + if name in term: + return ( term[0], term[1] ) + raise NameError("Colon annotation term '" + name + "' not found.") diff --git a/src/scaffoldmaker/annotation/heart_terms.py b/src/scaffoldmaker/annotation/heart_terms.py new file mode 100644 index 00000000..fd49ed50 --- /dev/null +++ b/src/scaffoldmaker/annotation/heart_terms.py @@ -0,0 +1,63 @@ +""" +Common resource for heart annotation terms. +""" + +# convention: preferred name, preferred id, followed by any other ids and alternative names +heart_terms = [ + ( "heart", "FMA:7088", "UBERON:0000948" ), + # ventricles + ( "left ventricle myocardium", "FMA:9558" ), + ( "right ventricle myocardium", "FMA:9535" ), + ( "interventricular septum", "FMA:7133" ), + ( "endocardium of left ventricle", "FMA:9559" ), + ( "endocardium of right ventricle", "FMA:9536" ), + ( "epicardium", "FMA:9461", "UBERON:0002348"), + #( "epicardium of ventricle", "FMA:12150", "UBERON:0001082" ), + # ventricles with base + ( "conus arteriosus", "UBERON:0003983" ), + ( "left fibrous ring", "FMA:77124" ), + ( "right fibrous ring", "FMA:77125" ), + # atria + ( "left atrium myocardium", "FMA:7285" ), + ( "right atrium myocardium", "FMA:7282" ), + ( "endocardium of left atrium", "FMA:7286", "UBERON:0034903" ), + ( "endocardium of right atrium", "FMA:7281", "UBERON:0009129" ), + ( "interatrial septum", "FMA:7108" ), + ( "fossa ovalis", "FMA:9246" ), + ( "left atrial appendage", "FMA:7219" ), # GRC rename auricle + ( "right atrial appendage", "FMA:7218" ), # GRC rename auricle + ( "pulmonary vein", "FMA:66643", "UBERON:0002016" ), + ( "left pulmonary vein", "UBERON:0009030" ), + ( "left inferior pulmonary vein", "FMA:49913" ), + ( "left superior pulmonary vein", "FMA:49916" ), + ( "middle pulmonary vein", "ILX:0739222" ), # in mouse, rat, rabbit + ( "right pulmonary vein", "UBERON:0009032" ), + ( "right inferior pulmonary vein", "FMA:49911" ), + ( "right superior pulmonary vein", "FMA:49914" ), + ( "inferior vena cava", "FMA:10951", "UBERON:0001072", "posterior vena cava" ), + ( "inferior vena cava inlet", "ILX:0738358" ), + ( "superior vena cava", "FMA:4720", "UBERON:0001585", "anterior vena cava" ), + ( "superior vena cava inlet", "ILX:0738367" ), + # arterial root + ( "root of aorta", "FMA:3740" ), + ( "posterior cusp of aortic valve", "FMA:7253" ), + ( "right cusp of aortic valve", "FMA:7252" ), + ( "left cusp of aortic valve", "FMA:7251" ), + ( "root of pulmonary trunk", "FMA:8612" ), + ( "right cusp of pulmonary valve", "FMA:7250" ), + ( "anterior cusp of pulmonary valve", "FMA:7249" ), + ( "left cusp of pulmonary valve", "FMA:7247" ), + # future: fiducial markers + ( "apex of heart", "FMA:7164", "UBERON:0002098") + ] + +def get_heart_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 heart_terms: + if name in term: + return ( term[0], term[1] ) + raise NameError("Heart annotation term '" + name + "' not found.") diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_bladder1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_bladder1.py index d80c76ba..8e6b607e 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_bladder1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_bladder1.py @@ -10,6 +10,7 @@ from opencmiss.zinc.field import Field from opencmiss.zinc.node import Node from scaffoldmaker.annotation.annotationgroup import AnnotationGroup +from scaffoldmaker.annotation.bladder_terms import get_bladder_term from scaffoldmaker.meshtypes.meshtype_3d_ostium1 import MeshType_3d_ostium1, generateOstiumMesh from scaffoldmaker.meshtypes.scaffold_base import Scaffold_base from scaffoldmaker.scaffoldpackage import ScaffoldPackage @@ -264,8 +265,8 @@ def generateBaseMesh(region, options): cache = fm.createFieldcache() - neckGroup = AnnotationGroup(region, 'neck of bladder', FMANumber='unknown', lyphID='unknown') - bodyGroup = AnnotationGroup(region, 'body of bladder', FMANumber='unknown', lyphID='unknown') + neckGroup = AnnotationGroup(region, get_bladder_term("neck of urinary bladder")) + bodyGroup = AnnotationGroup(region, get_bladder_term("body of urinary bladder")) annotationGroups = [neckGroup, bodyGroup] neckMeshGroup = neckGroup.getMeshGroup(mesh) diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_colonsegment1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_colonsegment1.py index 4a7a58b1..353889be 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_colonsegment1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_colonsegment1.py @@ -12,6 +12,7 @@ from opencmiss.zinc.field import Field from opencmiss.zinc.node import Node from scaffoldmaker.annotation.annotationgroup import AnnotationGroup +from scaffoldmaker.annotation.colon_terms import get_colon_term from scaffoldmaker.meshtypes.scaffold_base import Scaffold_base from scaffoldmaker.utils.eftfactory_bicubichermitelinear import eftfactory_bicubichermitelinear from scaffoldmaker.utils.eftfactory_tricubichermite import eftfactory_tricubichermite @@ -537,8 +538,8 @@ def getColonSegmentInnerPoints(region, elementsCountAroundTC, d2Final.append(d2Raw[n1][n2]) # Create annotation groups for mouse colon - mzGroup = AnnotationGroup(region, 'mesenteric zone', FMANumber = 'FMANumber unknown', lyphID = 'Lyph ID unknown') - nonmzGroup = AnnotationGroup(region, 'non-mesenteric zone', FMANumber = 'FMANumber unknown', lyphID = 'Lyph ID unknown') + 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) + ['non-mesenteric zone']*elementsCountAroundHaustrum + @@ -1288,14 +1289,14 @@ def getTeniaColi(region, xList, d1List, d2List, d3List, curvatureList, # Update annotation groups if tcCount == 3: - tlGroup = AnnotationGroup(region, 'tenia libera', FMANumber = 'FMANumber unknown', lyphID = 'Lyph ID unknown') - tmGroup = AnnotationGroup(region, 'tenia mesocolica', FMANumber = 'FMANumber unknown', lyphID = 'Lyph ID unknown') - toGroup = AnnotationGroup(region, 'tenia omentalis', FMANumber = 'FMANumber unknown', lyphID = 'Lyph ID unknown') + tlGroup = AnnotationGroup(region, get_colon_term("tenia libera")) + tmGroup = AnnotationGroup(region, get_colon_term("tenia mesocolica")) + toGroup = AnnotationGroup(region, get_colon_term("tenia omentalis")) annotationGroupsTC = [tlGroup, tmGroup, toGroup] - annotationArrayTC = (['tenia omentalis']*int(elementsCountAroundTC*0.5) + + annotationArrayTC = (['tenia omentalis']*(elementsCountAroundTC//2) + ['tenia libera']*elementsCountAroundTC + ['tenia mesocolica']*elementsCountAroundTC + - ['tenia omentalis']*int(elementsCountAroundTC*0.5)) + ['tenia omentalis']*(elementsCountAroundTC//2)) annotationGroups += annotationGroupsTC annotationArray += annotationArrayTC diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_heart1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_heart1.py index 7acaae84..33b890e2 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_heart1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_heart1.py @@ -10,7 +10,8 @@ from opencmiss.zinc.element import Element from opencmiss.zinc.field import Field from opencmiss.zinc.node import Node -from scaffoldmaker.annotation.annotationgroup import AnnotationGroup, findAnnotationGroupByName +from scaffoldmaker.annotation.annotationgroup import AnnotationGroup, findOrCreateAnnotationGroupForTerm, getAnnotationGroupForTerm +from scaffoldmaker.annotation.heart_terms import get_heart_term from scaffoldmaker.meshtypes.scaffold_base import Scaffold_base from scaffoldmaker.meshtypes.meshtype_3d_heartatria1 import MeshType_3d_heartatria1 from scaffoldmaker.meshtypes.meshtype_3d_heartventriclesbase1 import MeshType_3d_heartventriclesbase1 @@ -130,20 +131,18 @@ def generateBaseMesh(cls, region, options): # generate heartventriclesbase1 model and put atria1 on it annotationGroups = MeshType_3d_heartventriclesbase1.generateBaseMesh(region, options) annotationGroups += MeshType_3d_heartatria1.generateBaseMesh(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') - annotationGroups += [ lFibrousRingGroup, rFibrousRingGroup ] + #epiGroup = getAnnotationGroupForTerm(annotationGroups, get_heart_term("epicardium")) + lFibrousRingGroup = findOrCreateAnnotationGroupForTerm(annotationGroups, region, get_heart_term("left fibrous ring")) + rFibrousRingGroup = findOrCreateAnnotationGroupForTerm(annotationGroups, region, get_heart_term("right fibrous ring")) # annotation fiducial points markerGroup = findOrCreateFieldGroup(fm, "marker") - #markerCoordinates = findOrCreateFieldCoordinates(fm, "marker_coordinates") markerName = findOrCreateFieldStoredString(fm, name="marker_name") markerLocation = findOrCreateFieldStoredMeshLocation(fm, mesh, name="marker_location") nodes = fm.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES) markerPoints = findOrCreateFieldNodeGroup(markerGroup, nodes).getNodesetGroup() markerTemplateInternal = nodes.createNodetemplate() - #markerTemplateInternal.defineField(markerCoordinates) markerTemplateInternal.defineField(markerName) markerTemplateInternal.defineField(markerLocation) @@ -407,8 +406,7 @@ def generateBaseMesh(cls, region, options): markerPoint = markerPoints.createNode(nodeIdentifier, markerTemplateInternal) nodeIdentifier += 1 cache.setNode(markerPoint) - #markerCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, cruxCoordinates) - markerName.assignString(cache, 'crux') + markerName.assignString(cache, "crux of heart") markerLocation.assignMeshLocation(cache, cruxElement, cruxXi) fm.endChange() diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_heart2.py b/src/scaffoldmaker/meshtypes/meshtype_3d_heart2.py index 3220183d..7dfd9ae5 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_heart2.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_heart2.py @@ -5,7 +5,7 @@ from __future__ import division import math from opencmiss.utils.zinc.field import findOrCreateFieldCoordinates -from scaffoldmaker.annotation.annotationgroup import AnnotationGroup, findAnnotationGroupByName +from scaffoldmaker.annotation.annotationgroup import AnnotationGroup from scaffoldmaker.meshtypes.meshtype_3d_heartatria2 import MeshType_3d_heartatria2 from scaffoldmaker.meshtypes.meshtype_3d_heartventriclesbase2 import MeshType_3d_heartventriclesbase2 from scaffoldmaker.meshtypes.scaffold_base import Scaffold_base diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_heartarterialroot1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_heartarterialroot1.py index d2eb0488..61674e8f 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_heartarterialroot1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_heartarterialroot1.py @@ -11,7 +11,8 @@ from opencmiss.zinc.element import Element, Elementbasis from opencmiss.zinc.field import Field from opencmiss.zinc.node import Node -from scaffoldmaker.annotation.annotationgroup import AnnotationGroup, findAnnotationGroupByName +from scaffoldmaker.annotation.annotationgroup import AnnotationGroup +from scaffoldmaker.annotation.heart_terms import get_heart_term from scaffoldmaker.meshtypes.scaffold_base import Scaffold_base from scaffoldmaker.utils.eft_utils import remapEftLocalNodes, remapEftNodeValueLabel, setEftScaleFactorIds from scaffoldmaker.utils.geometry import getApproximateEllipsePerimeter, createCirclePoints @@ -121,17 +122,17 @@ def generateBaseMesh(cls, region, options, baseCentre=[ 0.0, 0.0, 0.0 ], axisSid mesh = fm.findMeshByDimension(3) if aorticNotPulmonary: - arterialRootGroup = AnnotationGroup(region, 'root of aorta', FMANumber = 3740, lyphID = 'Lyph ID unknown') + arterialRootGroup = AnnotationGroup(region, get_heart_term("root of aorta")) cuspGroups = [ - AnnotationGroup(region, 'posterior cusp of aortic valve', FMANumber = 7253, lyphID = 'Lyph ID unknown'), - AnnotationGroup(region, 'right cusp of aortic valve', FMANumber = 7252, lyphID = 'Lyph ID unknown'), - AnnotationGroup(region, 'left cusp of aortic valve', FMANumber = 7251, lyphID = 'Lyph ID unknown') ] + AnnotationGroup(region, get_heart_term("posterior cusp of aortic valve")), + AnnotationGroup(region, get_heart_term("right cusp of aortic valve")), + AnnotationGroup(region, get_heart_term("left cusp of aortic valve")) ] else: - arterialRootGroup = AnnotationGroup(region, 'root of pulmonary trunk', FMANumber = 8612, lyphID = 'Lyph ID unknown') + arterialRootGroup = AnnotationGroup(region, get_heart_term("root of pulmonary trunk")) cuspGroups = [ - AnnotationGroup(region, 'right cusp of pulmonary valve', FMANumber = 7250, lyphID = 'Lyph ID unknown'), - AnnotationGroup(region, 'anterior cusp of pulmonary valve', FMANumber = 7249, lyphID = 'Lyph ID unknown'), - AnnotationGroup(region, 'left cusp of pulmonary valve', FMANumber = 7247, lyphID = 'Lyph ID unknown') ] + AnnotationGroup(region, get_heart_term("right cusp of pulmonary valve")), + AnnotationGroup(region, get_heart_term("anterior cusp of pulmonary valve")), + AnnotationGroup(region, get_heart_term("left cusp of pulmonary valve")) ] allGroups = [ arterialRootGroup ] # groups that all elements in scaffold will go in annotationGroups = allGroups + cuspGroups diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_heartatria1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_heartatria1.py index 157c7262..b7a8091c 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_heartatria1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_heartatria1.py @@ -11,7 +11,8 @@ from opencmiss.zinc.element import Element, Elementbasis from opencmiss.zinc.field import Field from opencmiss.zinc.node import Node -from scaffoldmaker.annotation.annotationgroup import AnnotationGroup, findAnnotationGroupByName +from scaffoldmaker.annotation.annotationgroup import AnnotationGroup, findOrCreateAnnotationGroupForTerm, getAnnotationGroupForTerm +from scaffoldmaker.annotation.heart_terms import get_heart_term from scaffoldmaker.meshtypes.scaffold_base import Scaffold_base from scaffoldmaker.meshtypes.meshtype_3d_ostium1 import MeshType_3d_ostium1, generateOstiumMesh from scaffoldmaker.scaffoldpackage import ScaffoldPackage @@ -736,12 +737,12 @@ def generateBaseMesh(cls, region, options): coordinates = findOrCreateFieldCoordinates(fm) cache = fm.createFieldcache() - laGroup = AnnotationGroup(region, "left atrium myocardium", FMANumber = 7285, lyphID = "Lyph ID unknown") - raGroup = AnnotationGroup(region, "right atrium myocardium", FMANumber = 7282, lyphID = "Lyph ID unknown") - aSeptumGroup = AnnotationGroup(region, "interatrial septum", FMANumber = 7108, lyphID = "Lyph ID unknown") - fossaGroup = AnnotationGroup(region, "fossa ovalis", FMANumber = 9246, lyphID = "Lyph ID unknown") - laaGroup = AnnotationGroup(region, 'left atrial appendage', FMANumber = 7219, lyphID = 'Lyph ID unknown') - raaGroup = AnnotationGroup(region, 'right atrial appendage', FMANumber = 7218, lyphID = 'Lyph ID unknown') + laGroup = AnnotationGroup(region, get_heart_term("left atrium myocardium")) + raGroup = AnnotationGroup(region, get_heart_term("right atrium myocardium")) + aSeptumGroup = AnnotationGroup(region, get_heart_term("interatrial septum")) + fossaGroup = AnnotationGroup(region, get_heart_term("fossa ovalis")) + laaGroup = AnnotationGroup(region, get_heart_term("left atrial appendage")) + raaGroup = AnnotationGroup(region, get_heart_term("right atrial appendage")) annotationGroups = [ laGroup, raGroup, aSeptumGroup, fossaGroup, laaGroup, raaGroup ] lpvOstiumSettings = lpvOstium.getScaffoldSettings() @@ -749,52 +750,33 @@ def generateBaseMesh(cls, region, options): if commonLeftRightPvOstium: # use only lpv: if lpvCount == 1: - pvGroup = AnnotationGroup(region, 'pulmonary vein', FMANumber = 'FMA number unknown', lyphID = 'Lyph ID unknown') + pvGroup = AnnotationGroup(region, get_heart_term("pulmonary vein")) lpvGroups = [ pvGroup ] else: - lpvGroup = AnnotationGroup(region, 'left pulmonary vein', FMANumber = 'FMA number unknown', lyphID = 'Lyph ID unknown') + lpvGroup = AnnotationGroup(region, get_heart_term("left pulmonary vein")) + rpvGroup = AnnotationGroup(region, get_heart_term("right pulmonary vein")) if lpvCount == 2: - rpvGroup = AnnotationGroup(region, 'right pulmonary vein', FMANumber = 'FMA number unknown', lyphID = 'Lyph ID unknown') lpvGroups = [ lpvGroup, rpvGroup ] else: - rspvGroup = AnnotationGroup(region, 'right superior pulmonary vein', FMANumber = 49914, lyphID = 'Lyph ID unknown') - ripvGroup = AnnotationGroup(region, 'right inferior pulmonary vein', FMANumber = 49911, lyphID = 'Lyph ID unknown') - lpvGroups = [ lpvGroup, rspvGroup, ripvGroup ] + mpvGroup = AnnotationGroup(region, get_heart_term("middle pulmonary vein")) + lpvGroups = [ lpvGroup, mpvGroup, rpvGroup ] annotationGroups += lpvGroups else: # separate left and right pulmonary vein ostia - if lpvCount == 1: - lpvGroup = AnnotationGroup(region, 'left pulmonary vein', FMANumber = 'FMA number unknown', lyphID = 'Lyph ID unknown') - lpvGroups = [ lpvGroup ] - else: - lipvGroup = AnnotationGroup(region, 'left inferior pulmonary vein', FMANumber = 49913, lyphID = 'Lyph ID unknown') - lspvGroup = AnnotationGroup(region, 'left superior pulmonary vein', FMANumber = 49916, lyphID = 'Lyph ID unknown') - if lpvCount == 2: - lpvGroups = [ lipvGroup, lspvGroup ] - else: # lpvCount == 3: - lmpvGroup = AnnotationGroup(region, 'left middle pulmonary vein', FMANumber = 'FMA number unknown', lyphID = 'Lyph ID unknown') - lpvGroups = [ lipvGroup, lmpvGroup, lspvGroup ] - annotationGroups += lpvGroups + lpvGroup = AnnotationGroup(region, get_heart_term("left pulmonary vein")) + lpvGroups = [ lpvGroup ]*lpvCount + annotationGroups.append(lpvGroup) rpvOstiumSettings = rpvOstium.getScaffoldSettings() rpvCount = rpvOstiumSettings['Number of vessels'] - if rpvCount == 1: - rpvGroup = AnnotationGroup(region, 'right pulmonary vein', FMANumber = 'FMA number unknown', lyphID = 'Lyph ID unknown') - rpvGroups = [ rpvGroup ] - else: - ripvGroup = AnnotationGroup(region, 'right inferior pulmonary vein', FMANumber = 49911, lyphID = 'Lyph ID unknown') - rspvGroup = AnnotationGroup(region, 'right superior pulmonary vein', FMANumber = 49914, lyphID = 'Lyph ID unknown') - if rpvCount == 2: - rpvGroups = [ ripvGroup, rspvGroup ] - else: # rpvCount == 3: - rmpvGroup = AnnotationGroup(region, 'right middle pulmonary vein', FMANumber = 'FMA number unknown', lyphID = 'Lyph ID unknown') - rpvGroups = [ ripvGroup, rmpvGroup, rspvGroup ] - annotationGroups += rpvGroups - - ivcInletGroup = AnnotationGroup(region, 'inferior vena cava inlet', FMANumber = 10951, lyphID = 'Lyph ID unknown') - svcInletGroup = AnnotationGroup(region, 'superior vena cava inlet', FMANumber = 4720, lyphID = 'Lyph ID unknown') + rpvGroup = AnnotationGroup(region, get_heart_term("right pulmonary vein" )) + rpvGroups = [ rpvGroup ]*rpvCount + annotationGroups.append(rpvGroup) + + ivcInletGroup = AnnotationGroup(region, get_heart_term("inferior vena cava inlet")) + svcInletGroup = AnnotationGroup(region, get_heart_term("superior vena cava inlet")) annotationGroups += [ ivcInletGroup, svcInletGroup ] # av boundary nodes are put in left and right fibrous ring groups only so they can be found by heart1 - lFibrousRingGroup = AnnotationGroup(region, 'left fibrous ring', FMANumber = 77124, lyphID = 'Lyph ID unknown') - rFibrousRingGroup = AnnotationGroup(region, 'right fibrous ring', FMANumber = 77125, lyphID = 'Lyph ID unknown') + lFibrousRingGroup = AnnotationGroup(region, get_heart_term("left fibrous ring")) + rFibrousRingGroup = AnnotationGroup(region, get_heart_term("right fibrous ring")) ############## # Create nodes @@ -3180,15 +3162,14 @@ def generateBaseMesh(cls, region, options): fm.createFieldNot(is_la_endo)), fm.createFieldAnd(aSeptumGroup.getFieldElementGroup(mesh2d), is_exterior_face_xi3_1)) is_a_epi = fm.createFieldAnd(fm.createFieldOr(is_la, is_ra), - fm.createFieldAnd(is_exterior_face_xi3_1, - fm.createFieldNot(aSeptumGroup.getFieldElementGroup(mesh2d)))) - laEndoGroup = AnnotationGroup(region, "Endocardium of left atrium", FMANumber = 7286, lyphID = 'Lyph ID unknown') + fm.createFieldAnd(is_exterior_face_xi3_1, + fm.createFieldNot(aSeptumGroup.getFieldElementGroup(mesh2d)))) + epiGroup = findOrCreateAnnotationGroupForTerm(annotationGroups, region, get_heart_term("epicardium")) + epiGroup.getMeshGroup(mesh2d).addElementsConditional(is_a_epi) + laEndoGroup = findOrCreateAnnotationGroupForTerm(annotationGroups, region, get_heart_term("endocardium of left atrium")) laEndoGroup.getMeshGroup(mesh2d).addElementsConditional(is_la_endo) - raEndoGroup = AnnotationGroup(region, "Endocardium of right atrium", FMANumber = 7281, lyphID = 'Lyph ID unknown') + raEndoGroup = findOrCreateAnnotationGroupForTerm(annotationGroups, region, get_heart_term("endocardium of right atrium")) raEndoGroup.getMeshGroup(mesh2d).addElementsConditional(is_ra_endo) - # Note I could not find any epicardium of atria - aEpiGroup = AnnotationGroup(region, "Epicardium of atrium", FMANumber = 0, lyphID = 'Lyph ID unknown') - aEpiGroup.getMeshGroup(mesh2d).addElementsConditional(is_a_epi) del is_exterior del is_exterior_face_xi3_0 del is_exterior_face_xi3_1 @@ -3197,7 +3178,6 @@ def generateBaseMesh(cls, region, options): del is_la_endo del is_ra_endo del is_a_epi - annotationGroups += [ laEndoGroup, raEndoGroup, aEpiGroup ] fm.endChange() return annotationGroups @@ -3218,16 +3198,16 @@ def refineMesh(cls, meshrefinement, options): sourceFm = meshrefinement._sourceFm annotationGroups = meshrefinement._sourceAnnotationGroups - laGroup = findAnnotationGroupByName(annotationGroups, 'left atrium myocardium') + laGroup = getAnnotationGroupForTerm(annotationGroups, get_heart_term("left atrium myocardium")) laElementGroupField = laGroup.getFieldElementGroup(meshrefinement._sourceMesh) - raGroup = findAnnotationGroupByName(annotationGroups, 'right atrium myocardium') + raGroup = getAnnotationGroupForTerm(annotationGroups, get_heart_term("right atrium myocardium")) raElementGroupField = raGroup.getFieldElementGroup(meshrefinement._sourceMesh) - aSeptumGroup = findAnnotationGroupByName(annotationGroups, 'interatrial septum') + aSeptumGroup = getAnnotationGroupForTerm(annotationGroups, get_heart_term("interatrial septum")) aSeptumElementGroupField = aSeptumGroup.getFieldElementGroup(meshrefinement._sourceMesh) isSeptumEdgeWedge = sourceFm.createFieldXor(sourceFm.createFieldAnd(laElementGroupField, raElementGroupField), aSeptumElementGroupField) # last atria element is last element in following group: - lastGroup = findAnnotationGroupByName(annotationGroups, 'right atrial appendage') + lastGroup = getAnnotationGroupForTerm(annotationGroups, get_heart_term("right atrial appendage")) lastMeshGroup = lastGroup.getMeshGroup(meshrefinement._sourceMesh) lastElementIdentifier = -1 elementIter = lastMeshGroup.createElementiterator() diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_heartatria2.py b/src/scaffoldmaker/meshtypes/meshtype_3d_heartatria2.py index 9414e290..61119f40 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_heartatria2.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_heartatria2.py @@ -11,6 +11,7 @@ from opencmiss.zinc.field import Field from opencmiss.zinc.node import Node from scaffoldmaker.annotation.annotationgroup import AnnotationGroup +from scaffoldmaker.annotation.heart_terms import get_heart_term from scaffoldmaker.meshtypes.scaffold_base import Scaffold_base from scaffoldmaker.utils.eft_utils import remapEftLocalNodes, remapEftNodeValueLabel, scaleEftNodeValueLabels, setEftScaleFactorIds from scaffoldmaker.utils.geometry import getApproximateEllipsePerimeter, getEllipseArcLength, getEllipseRadiansToX, updateEllipseAngleByArcLength @@ -190,16 +191,16 @@ def generateBaseMesh(cls, region, options): coordinates = findOrCreateFieldCoordinates(fm) cache = fm.createFieldcache() - laGroup = AnnotationGroup(region, "left atrium myocardium", FMANumber = 7285, lyphID = "Lyph ID unknown") - raGroup = AnnotationGroup(region, "right atrium myocardium", FMANumber = 7282, lyphID = "Lyph ID unknown") - aSeptumGroup = AnnotationGroup(region, "interatrial septum", FMANumber = 7108, lyphID = "Lyph ID unknown") - fossaGroup = AnnotationGroup(region, "fossa ovalis", FMANumber = 9246, lyphID = "Lyph ID unknown") - lipvGroup = AnnotationGroup(region, 'left inferior pulmonary vein', FMANumber = 49913, lyphID = 'Lyph ID unknown') - lspvGroup = AnnotationGroup(region, 'left superior pulmonary vein', FMANumber = 49916, lyphID = 'Lyph ID unknown') - ripvGroup = AnnotationGroup(region, 'right inferior pulmonary vein', FMANumber = 49911, lyphID = 'Lyph ID unknown') - rspvGroup = AnnotationGroup(region, 'right superior pulmonary vein', FMANumber = 49914, lyphID = 'Lyph ID unknown') - ivcInletGroup = AnnotationGroup(region, 'inferior vena cava inlet', FMANumber = 10951, lyphID = 'Lyph ID unknown') - svcInletGroup = AnnotationGroup(region, 'superior vena cava inlet', FMANumber = 4720, lyphID = 'Lyph ID unknown') + laGroup = AnnotationGroup(region, get_heart_term("left atrium myocardium")) + raGroup = AnnotationGroup(region, get_heart_term("right atrium myocardium")) + aSeptumGroup = AnnotationGroup(region, get_heart_term("interatrial septum")) + fossaGroup = AnnotationGroup(region, get_heart_term("fossa ovalis")) + lipvGroup = AnnotationGroup(region, get_heart_term("left inferior pulmonary vein")) + lspvGroup = AnnotationGroup(region, get_heart_term("left superior pulmonary vein")) + ripvGroup = AnnotationGroup(region, get_heart_term("right inferior pulmonary vein")) + rspvGroup = AnnotationGroup(region, get_heart_term("right superior pulmonary vein")) + ivcInletGroup = AnnotationGroup(region, get_heart_term("inferior vena cava inlet")) + svcInletGroup = AnnotationGroup(region, get_heart_term("superior vena cava inlet")) annotationGroups = [ laGroup, raGroup, aSeptumGroup, fossaGroup, lipvGroup, lspvGroup, ripvGroup, rspvGroup, ivcInletGroup, svcInletGroup ] ############## diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_heartventricles1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_heartventricles1.py index 114a00bf..f41d8769 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_heartventricles1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_heartventricles1.py @@ -11,6 +11,7 @@ from opencmiss.zinc.field import Field from opencmiss.zinc.node import Node from scaffoldmaker.annotation.annotationgroup import AnnotationGroup +from scaffoldmaker.annotation.heart_terms import get_heart_term from scaffoldmaker.meshtypes.scaffold_base import Scaffold_base from scaffoldmaker.utils import vector from scaffoldmaker.utils.eft_utils import remapEftLocalNodes, remapEftNodeValueLabel, scaleEftNodeValueLabels, setEftScaleFactorIds @@ -249,9 +250,9 @@ def generateBaseMesh(cls, region, options): mesh = fm.findMeshByDimension(3) - lvGroup = AnnotationGroup(region, "left ventricle myocardium", FMANumber = 9558, lyphID = 'Lyph ID unknown') - rvGroup = AnnotationGroup(region, "right ventricle myocardium", FMANumber = 9535, lyphID = 'Lyph ID unknown') - vSeptumGroup = AnnotationGroup(region, "interventricular septum", FMANumber = 7133, lyphID = 'Lyph ID unknown') + lvGroup = AnnotationGroup(region, get_heart_term("left ventricle myocardium")) + rvGroup = AnnotationGroup(region, get_heart_term("right ventricle myocardium")) + vSeptumGroup = AnnotationGroup(region, get_heart_term("interventricular septum")) annotationGroups = [ lvGroup, rvGroup, vSeptumGroup ] # annotation fiducial points @@ -1052,12 +1053,12 @@ def generateBaseMesh(cls, region, options): fm.createFieldAnd(vSeptumGroup.getFieldElementGroup(mesh2d), is_exterior_face_xi3_1)) is_epi = fm.createFieldAnd(is_exterior_face_xi3_1, fm.createFieldNot(vSeptumGroup.getFieldElementGroup(mesh2d))) - lvEndoGroup = AnnotationGroup(region, "Endocardium of left ventricle", FMANumber = 9559, lyphID = 'Lyph ID unknown') + epiGroup = AnnotationGroup(region, get_heart_term("epicardium")) + epiGroup.getMeshGroup(mesh2d).addElementsConditional(is_epi) + lvEndoGroup = AnnotationGroup(region, get_heart_term("endocardium of left ventricle")) lvEndoGroup.getMeshGroup(mesh2d).addElementsConditional(is_lv_endo) - rvEndoGroup = AnnotationGroup(region, "Endocardium of right ventricle", FMANumber = 9536, lyphID = 'Lyph ID unknown') + rvEndoGroup = AnnotationGroup(region, get_heart_term("endocardium of right ventricle")) rvEndoGroup.getMeshGroup(mesh2d).addElementsConditional(is_rv_endo) - vEpiGroup = AnnotationGroup(region, "Epicardium of ventricle", FMANumber = 12150, lyphID = 'Lyph ID unknown') - vEpiGroup.getMeshGroup(mesh2d).addElementsConditional(is_epi) del is_exterior del is_exterior_face_xi3_0 del is_exterior_face_xi3_1 @@ -1066,21 +1067,21 @@ def generateBaseMesh(cls, region, options): del is_lv_endo del is_rv_endo del is_epi - annotationGroups += [ lvEndoGroup, rvEndoGroup, vEpiGroup ] + annotationGroups += [ lvEndoGroup, rvEndoGroup, epiGroup ] # apex annotation points element1 = mesh.findElementByIdentifier(1) - markerPoint = markerPoints.createNode(nodeIdentifier, markerTemplateInternal) - nodeIdentifier += 1 - cache.setNode(markerPoint) - #markerCoordinates.assignReal(cache, lvApexInnerx) - markerName.assignString(cache, 'apex endo') - markerLocation.assignMeshLocation(cache, element1, [ 0.0, 0.0, 0.0 ]) + #markerPoint = markerPoints.createNode(nodeIdentifier, markerTemplateInternal) + #nodeIdentifier += 1 + #cache.setNode(markerPoint) + ##markerCoordinates.assignReal(cache, lvApexInnerx) + #markerName.assignString(cache, 'apex endo') + #markerLocation.assignMeshLocation(cache, element1, [ 0.0, 0.0, 0.0 ]) markerPoint = markerPoints.createNode(nodeIdentifier, markerTemplateInternal) nodeIdentifier += 1 cache.setNode(markerPoint) #markerCoordinates.assignReal(cache, lvApexOuterx) - markerName.assignString(cache, 'apex epi') + markerName.assignString(cache, 'APEX') # interlex.org has 'apex of heart' markerLocation.assignMeshLocation(cache, element1, [ 0.0, 0.0, 1.0 ]) fm.endChange() diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_heartventricles2.py b/src/scaffoldmaker/meshtypes/meshtype_3d_heartventricles2.py index 86365ec6..f206c2d1 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_heartventricles2.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_heartventricles2.py @@ -9,6 +9,7 @@ from opencmiss.zinc.field import Field from opencmiss.zinc.node import Node from scaffoldmaker.annotation.annotationgroup import AnnotationGroup +from scaffoldmaker.annotation.heart_terms import get_heart_term from scaffoldmaker.meshtypes.scaffold_base import Scaffold_base from scaffoldmaker.utils import vector from scaffoldmaker.utils.eft_utils import remapEftNodeValueLabel, setEftScaleFactorIds @@ -153,9 +154,9 @@ def generateBaseMesh(cls, region, options): coordinates = findOrCreateFieldCoordinates(fm) cache = fm.createFieldcache() - lvGroup = AnnotationGroup(region, "left ventricle myocardium", FMANumber = 9558, lyphID = 'Lyph ID unknown') - rvGroup = AnnotationGroup(region, "right ventricle myocardium", FMANumber = 9535, lyphID = 'Lyph ID unknown') - vSeptumGroup = AnnotationGroup(region, "interventricular septum", FMANumber = 7133, lyphID = 'Lyph ID unknown') + lvGroup = AnnotationGroup(region, get_heart_term("left ventricle myocardium")) + rvGroup = AnnotationGroup(region, get_heart_term("right ventricle myocardium")) + vSeptumGroup = AnnotationGroup(region, get_heart_term("interventricular septum")) annotationGroups = [ lvGroup, rvGroup, vSeptumGroup ] ################# diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_heartventriclesbase1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_heartventriclesbase1.py index b17d5bec..fb502221 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_heartventriclesbase1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_heartventriclesbase1.py @@ -14,7 +14,8 @@ from opencmiss.zinc.field import Field from opencmiss.zinc.node import Node from opencmiss.zinc.result import RESULT_OK as ZINC_OK -from scaffoldmaker.annotation.annotationgroup import AnnotationGroup, findAnnotationGroupByName +from scaffoldmaker.annotation.annotationgroup import AnnotationGroup, getAnnotationGroupForTerm +from scaffoldmaker.annotation.heart_terms import get_heart_term from scaffoldmaker.meshtypes.meshtype_3d_heartatria1 import MeshType_3d_heartatria1, getAtriumBasePoints from scaffoldmaker.meshtypes.meshtype_3d_heartventricles1 import MeshType_3d_heartventricles1 from scaffoldmaker.meshtypes.scaffold_base import Scaffold_base @@ -300,14 +301,14 @@ def generateBaseMesh(cls, region, options): annotationGroups = MeshType_3d_heartventricles1.generateBaseMesh(region, options) # find/add annotation groups - lvGroup = findAnnotationGroupByName(annotationGroups, "left ventricle myocardium") - rvGroup = findAnnotationGroupByName(annotationGroups, "right ventricle myocardium") - vSeptumGroup = findAnnotationGroupByName(annotationGroups, "interventricular septum") - conusArteriosusGroup = AnnotationGroup(region, 'conus arteriosus', FMANumber = 0, lyphID = 'Lyph ID unknown') + lvGroup = getAnnotationGroupForTerm(annotationGroups, get_heart_term("left ventricle myocardium")) + rvGroup = getAnnotationGroupForTerm(annotationGroups, get_heart_term("right ventricle myocardium")) + vSeptumGroup = getAnnotationGroupForTerm(annotationGroups, get_heart_term("interventricular septum")) + conusArteriosusGroup = AnnotationGroup(region, get_heart_term("conus arteriosus")) annotationGroups += [ conusArteriosusGroup ] # av boundary nodes are put in left and right fibrous ring groups only so they can be found by heart1 - lFibrousRingGroup = AnnotationGroup(region, 'left fibrous ring', FMANumber = 77124, lyphID = 'Lyph ID unknown') - rFibrousRingGroup = AnnotationGroup(region, 'right fibrous ring', FMANumber = 77125, lyphID = 'Lyph ID unknown') + lFibrousRingGroup = AnnotationGroup(region, get_heart_term("left fibrous ring")) + rFibrousRingGroup = AnnotationGroup(region, get_heart_term("right fibrous ring")) # annotation fiducial points markerGroup = findOrCreateFieldGroup(fm, "marker") @@ -321,6 +322,10 @@ def generateBaseMesh(cls, region, options): markerTemplateExternal.defineField(markerCoordinates) markerTemplateExternal.defineField(markerName) + markerTemplateInternal = nodes.createNodetemplate() + markerTemplateInternal.defineField(markerName) + markerTemplateInternal.defineField(markerLocation) + ################# # Create nodes ################# @@ -857,6 +862,7 @@ def generateBaseMesh(cls, region, options): nids = None scalefactors = None meshGroups = [ rvMeshGroup ] + addMarker = None noa = e niv = e @@ -976,6 +982,7 @@ def generateBaseMesh(cls, region, options): # outer infundibulum 5, above septum nids = [ rvInnerNodeId[niv - 1], rvInnerNodeId[niv], rvOutletNodeId[0][5], rvOutletNodeId[0][0], avsNodeId, lvOutletNodeId[1][3], rvOutletNodeId[1][5], rvOutletNodeId[1][0] ] + addMarker = { "name" : "Pulmonary valve-RV", "xi" : [ 1.0, 1.0, 0.0 ] } eft1 = tricubichermite.createEftNoCrossDerivatives() setEftScaleFactorIds(eft1, [1], []) scalefactors = [ -1.0 ] @@ -1002,6 +1009,13 @@ def generateBaseMesh(cls, region, options): #print('create element rv base r1', elementIdentifier, result, result2, result3, nids) elementIdentifier += 1 + if addMarker: + markerPoint = markerPoints.createNode(nodeIdentifier, markerTemplateInternal) + nodeIdentifier += 1 + cache.setNode(markerPoint) + markerName.assignString(cache, addMarker["name"]) + markerLocation.assignMeshLocation(cache, element, addMarker["xi"]) + for meshGroup in meshGroups: meshGroup.addElement(element) @@ -1117,6 +1131,7 @@ def generateBaseMesh(cls, region, options): nids = None scalefactors = None meshGroups = [ lvMeshGroup ] + addMarker = None eft1 = tricubichermite.createEftNoCrossDerivatives() setEftScaleFactorIds(eft1, [1], []) @@ -1200,6 +1215,8 @@ def generateBaseMesh(cls, region, options): no = e - 5 ni = elementsCountAroundLVFreeWall + elementsCountAroundAtrialSeptum + no nids = [ lvInnerNodeId[ni], lvInnerNodeId[ni + 1], lvOutletNodeId[0][no], lvOutletNodeId[0][no + 1], lvOutletNodeId[1][no], lvOutletNodeId[1][no + 1] ] + if nids[2] == lvOutletNodeId[0][2]: + addMarker = { "name" : "Aortic Valve-Coronary vessel", "xi" : [ 0.0, 1.0, 0.0 ] } if no == 0: remapEftNodeValueLabel(eft1, [ 1 ], Node.VALUE_LABEL_D_DS2, [ ( Node.VALUE_LABEL_D_DS1, [] ), ( Node.VALUE_LABEL_D_DS2, [] ) ]) remapEftNodeValueLabel(eft1, [ 1 ], Node.VALUE_LABEL_D_DS3, [ ( Node.VALUE_LABEL_D_DS2, [] ), ( Node.VALUE_LABEL_D_DS3, [] ) ]) @@ -1223,6 +1240,13 @@ def generateBaseMesh(cls, region, options): #print('create element lv base r2', elementIdentifier, result, result2, result3, nids) elementIdentifier += 1 + if addMarker: + markerPoint = markerPoints.createNode(nodeIdentifier, markerTemplateInternal) + nodeIdentifier += 1 + cache.setNode(markerPoint) + markerName.assignString(cache, addMarker["name"]) + markerLocation.assignMeshLocation(cache, element, addMarker["xi"]) + for meshGroup in meshGroups: meshGroup.addElement(element) @@ -1378,12 +1402,12 @@ def generateBaseMesh(cls, region, options): fm.createFieldAnd(vSeptumGroup.getFieldElementGroup(mesh2d), is_exterior_face_xi3_1)) is_epi = fm.createFieldAnd(is_exterior_face_xi3_1, fm.createFieldNot(vSeptumGroup.getFieldElementGroup(mesh2d))) - lvEndoGroup = AnnotationGroup(region, "Endocardium of left ventricle", FMANumber = 9559, lyphID = 'Lyph ID unknown') + epiGroup = getAnnotationGroupForTerm(annotationGroups, get_heart_term("epicardium")) + epiGroup.getMeshGroup(mesh2d).addElementsConditional(is_epi) + lvEndoGroup = getAnnotationGroupForTerm(annotationGroups, get_heart_term("endocardium of left ventricle")) lvEndoGroup.getMeshGroup(mesh2d).addElementsConditional(is_lv_endo) - rvEndoGroup = AnnotationGroup(region, "Endocardium of right ventricle", FMANumber = 9536, lyphID = 'Lyph ID unknown') + rvEndoGroup = getAnnotationGroupForTerm(annotationGroups, get_heart_term("endocardium of right ventricle")) rvEndoGroup.getMeshGroup(mesh2d).addElementsConditional(is_rv_endo) - vEpiGroup = AnnotationGroup(region, "Epicardium of ventricle", FMANumber = 12150, lyphID = 'Lyph ID unknown') - vEpiGroup.getMeshGroup(mesh2d).addElementsConditional(is_epi) del is_exterior del is_exterior_face_xi3_0 del is_exterior_face_xi3_1 diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_heartventriclesbase2.py b/src/scaffoldmaker/meshtypes/meshtype_3d_heartventriclesbase2.py index cd46dfd4..b36e91c9 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_heartventriclesbase2.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_heartventriclesbase2.py @@ -11,7 +11,8 @@ from opencmiss.zinc.element import Element, Elementbasis from opencmiss.zinc.field import Field from opencmiss.zinc.node import Node -from scaffoldmaker.annotation.annotationgroup import AnnotationGroup, findAnnotationGroupByName +from scaffoldmaker.annotation.annotationgroup import AnnotationGroup, getAnnotationGroupForTerm +from scaffoldmaker.annotation.heart_terms import get_heart_term from scaffoldmaker.meshtypes.meshtype_3d_heartventricles2 import MeshType_3d_heartventricles2 from scaffoldmaker.meshtypes.scaffold_base import Scaffold_base from scaffoldmaker.utils.eft_utils import remapEftLocalNodes, remapEftNodeValueLabel, scaleEftNodeValueLabels, setEftScaleFactorIds @@ -181,13 +182,14 @@ def generateBaseMesh(cls, region, options): # generate heartventricles2 model to add base plane to annotationGroups = MeshType_3d_heartventricles2.generateBaseMesh(region, options) - lvGroup = findAnnotationGroupByName(annotationGroups, "left ventricle myocardium") - rvGroup = findAnnotationGroupByName(annotationGroups, "right ventricle myocardium") - vSeptumGroup = findAnnotationGroupByName(annotationGroups, "interventricular septum") - conusArteriosusGroup = AnnotationGroup(region, 'conus arteriosus', FMANumber = 0, lyphID = 'Lyph ID unknown') - lFibrousRingGroup = AnnotationGroup(region, 'left fibrous ring', FMANumber = 77124, lyphID = 'Lyph ID unknown') - rFibrousRingGroup = AnnotationGroup(region, 'right fibrous ring', FMANumber = 77125, lyphID = 'Lyph ID unknown') - annotationGroups += [ conusArteriosusGroup, lFibrousRingGroup, rFibrousRingGroup ] + lvGroup = getAnnotationGroupForTerm(annotationGroups, get_heart_term("left ventricle myocardium")) + rvGroup = getAnnotationGroupForTerm(annotationGroups, get_heart_term("right ventricle myocardium")) + vSeptumGroup = getAnnotationGroupForTerm(annotationGroups, get_heart_term("interventricular septum")) + conusArteriosusGroup = AnnotationGroup(region, get_heart_term("conus arteriosus")) + annotationGroups += [ conusArteriosusGroup ] + # av boundary nodes are put in left and right fibrous ring groups only so they can be found by heart1 + lFibrousRingGroup = AnnotationGroup(region, get_heart_term("left fibrous ring")) + rFibrousRingGroup = AnnotationGroup(region, get_heart_term("right fibrous ring")) fm = region.getFieldmodule() fm.beginChange() diff --git a/src/scaffoldmaker/utils/meshrefinement.py b/src/scaffoldmaker/utils/meshrefinement.py index 38383994..08c36e22 100644 --- a/src/scaffoldmaker/utils/meshrefinement.py +++ b/src/scaffoldmaker/utils/meshrefinement.py @@ -74,8 +74,7 @@ def __init__(self, sourceRegion, targetRegion, sourceAnnotationGroups = []): self._sourceAndTargetMeshGroups = [] for sourceAnnotationGroup in sourceAnnotationGroups: sourceMeshGroup = sourceAnnotationGroup.getMeshGroup(self._sourceMesh) - targetAnnotationGroup = AnnotationGroup(self._targetRegion, \ - sourceAnnotationGroup.getName(), sourceAnnotationGroup.getFMANumber(), sourceAnnotationGroup.getLyphID()) + targetAnnotationGroup = AnnotationGroup(self._targetRegion, sourceAnnotationGroup.getTerm()) targetMeshGroup = targetAnnotationGroup.getMeshGroup(self._targetMesh) self._annotationGroups.append(targetAnnotationGroup) self._sourceAndTargetMeshGroups.append( ( sourceMeshGroup, targetMeshGroup) ) From 43074fd182a5b1c1ea4fb8e407d02ac2996e19d2 Mon Sep 17 00:00:00 2001 From: Richard Christie Date: Tue, 21 Apr 2020 10:30:07 +1200 Subject: [PATCH 04/15] Add fibrous ring epicardium --- .../meshtypes/meshtype_3d_heart1.py | 17 ++++++++++++++++- 1 file changed, 16 insertions(+), 1 deletion(-) diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_heart1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_heart1.py index 33b890e2..c36945ed 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_heart1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_heart1.py @@ -131,7 +131,6 @@ def generateBaseMesh(cls, region, options): # generate heartventriclesbase1 model and put atria1 on it annotationGroups = MeshType_3d_heartventriclesbase1.generateBaseMesh(region, options) annotationGroups += MeshType_3d_heartatria1.generateBaseMesh(region, options) - #epiGroup = getAnnotationGroupForTerm(annotationGroups, get_heart_term("epicardium")) lFibrousRingGroup = findOrCreateAnnotationGroupForTerm(annotationGroups, region, get_heart_term("left fibrous ring")) rFibrousRingGroup = findOrCreateAnnotationGroupForTerm(annotationGroups, region, get_heart_term("right fibrous ring")) @@ -409,6 +408,22 @@ def generateBaseMesh(cls, region, options): markerName.assignString(cache, "crux of heart") markerLocation.assignMeshLocation(cache, cruxElement, cruxXi) + # add to epicardium surface group + fm.defineAllFaces() + lFibrousRingGroup.addSubelements() + rFibrousRingGroup.addSubelements() + mesh2d = fm.findMeshByDimension(2) + is_exterior = fm.createFieldIsExterior() + is_exterior_face_xi3_1 = fm.createFieldAnd(is_exterior, fm.createFieldIsOnFace(Element.FACE_TYPE_XI3_1)) + is_fibrous_ring = fm.createFieldOr(lFibrousRingGroup.getFieldElementGroup(mesh2d), rFibrousRingGroup.getFieldElementGroup(mesh2d)) + is_fibrous_ring_epi = fm.createFieldAnd(is_fibrous_ring, is_exterior_face_xi3_1) + epiGroup = getAnnotationGroupForTerm(annotationGroups, get_heart_term("epicardium")) + epiGroup.getMeshGroup(mesh2d).addElementsConditional(is_fibrous_ring_epi) + del is_exterior + del is_exterior_face_xi3_1 + del is_fibrous_ring + del is_fibrous_ring_epi + fm.endChange() return annotationGroups From a436b53a1bc68be56b21ce33286eb3ea01462983 Mon Sep 17 00:00:00 2001 From: Richard Christie Date: Tue, 21 Apr 2020 11:56:35 +1200 Subject: [PATCH 05/15] Merge heart annotation group list without duplicates --- src/scaffoldmaker/annotation/annotationgroup.py | 16 ++++++++++++++++ .../meshtypes/meshtype_3d_heart1.py | 7 ++++--- 2 files changed, 20 insertions(+), 3 deletions(-) diff --git a/src/scaffoldmaker/annotation/annotationgroup.py b/src/scaffoldmaker/annotation/annotationgroup.py index 96b83797..cfa48004 100644 --- a/src/scaffoldmaker/annotation/annotationgroup.py +++ b/src/scaffoldmaker/annotation/annotationgroup.py @@ -143,3 +143,19 @@ def getAnnotationGroupForTerm(annotationGroups: list, term) -> AnnotationGroup: assert annotationGroup._id == term[1], "Annotation group '" + name + "' id '" + term[1] + "' does not match existing id '" + annotationGroup._id + "'" return annotationGroup raise NameError("Annotation group '" + name + "' not found.") + + +def mergeAnnotationGroups(*annotationGroupsIn): + ''' + Merge the supplied sequence of list(annotationGroups) to a single list, + without duplicates. + :param annotationGroupsIn: Variable number of list(AnnotationGroup) to merge. + Groups must be for the same region. + :return: Merged list(AnnotationGroup) + ''' + annotationGroups = [] + for agroups in annotationGroupsIn: + for agroup in agroups: + if not findAnnotationGroupByName(annotationGroups, agroup._name): + annotationGroups.append(agroup) + return annotationGroups diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_heart1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_heart1.py index c36945ed..bb16d599 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_heart1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_heart1.py @@ -10,7 +10,7 @@ from opencmiss.zinc.element import Element from opencmiss.zinc.field import Field from opencmiss.zinc.node import Node -from scaffoldmaker.annotation.annotationgroup import AnnotationGroup, findOrCreateAnnotationGroupForTerm, getAnnotationGroupForTerm +from scaffoldmaker.annotation.annotationgroup import AnnotationGroup, findOrCreateAnnotationGroupForTerm, getAnnotationGroupForTerm, mergeAnnotationGroups from scaffoldmaker.annotation.heart_terms import get_heart_term from scaffoldmaker.meshtypes.scaffold_base import Scaffold_base from scaffoldmaker.meshtypes.meshtype_3d_heartatria1 import MeshType_3d_heartatria1 @@ -129,8 +129,9 @@ def generateBaseMesh(cls, region, options): mesh = fm.findMeshByDimension(3) # generate heartventriclesbase1 model and put atria1 on it - annotationGroups = MeshType_3d_heartventriclesbase1.generateBaseMesh(region, options) - annotationGroups += MeshType_3d_heartatria1.generateBaseMesh(region, options) + ventriclesAnnotationGroups = MeshType_3d_heartventriclesbase1.generateBaseMesh(region, options) + atriaAnnotationGroups = MeshType_3d_heartatria1.generateBaseMesh(region, options) + annotationGroups = mergeAnnotationGroups(ventriclesAnnotationGroups, atriaAnnotationGroups) lFibrousRingGroup = findOrCreateAnnotationGroupForTerm(annotationGroups, region, get_heart_term("left fibrous ring")) rFibrousRingGroup = findOrCreateAnnotationGroupForTerm(annotationGroups, region, get_heart_term("right fibrous ring")) From 7340b7f9ae6353225b9deed1740f93147e6e5c6a Mon Sep 17 00:00:00 2001 From: Richard Christie Date: Tue, 21 Apr 2020 12:27:26 +1200 Subject: [PATCH 06/15] Remove non-embedded marker points --- .../meshtype_3d_heartarterialroot1.py | 18 ------- .../meshtypes/meshtype_3d_heartventricles1.py | 9 ---- .../meshtype_3d_heartventriclesbase1.py | 47 ------------------- 3 files changed, 74 deletions(-) diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_heartarterialroot1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_heartarterialroot1.py index 61674e8f..6e05273b 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_heartarterialroot1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_heartarterialroot1.py @@ -137,17 +137,7 @@ def generateBaseMesh(cls, region, options, baseCentre=[ 0.0, 0.0, 0.0 ], axisSid allGroups = [ arterialRootGroup ] # groups that all elements in scaffold will go in annotationGroups = allGroups + cuspGroups - # annotation fiducial points - markerGroup = findOrCreateFieldGroup(fm, "marker") - markerCoordinates = findOrCreateFieldCoordinates(fm, "marker_coordinates") - markerName = findOrCreateFieldStoredString(fm, name="marker_name") - #markerLocation = findOrCreateFieldStoredMeshLocation(fm, mesh, name="marker_location") - nodes = fm.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES) - markerPoints = findOrCreateFieldNodeGroup(markerGroup, nodes).getNodesetGroup() - markerTemplateExternal = nodes.createNodetemplate() - markerTemplateExternal.defineField(markerCoordinates) - markerTemplateExternal.defineField(markerName) ################# # Create nodes @@ -506,14 +496,6 @@ def generateBaseMesh(cls, region, options, baseCentre=[ 0.0, 0.0, 0.0 ], axisSid for meshGroup in meshGroups: meshGroup.addElement(element) - # create annotation points - - markerPoint = markerPoints.createNode(nodeIdentifier, markerTemplateExternal) - nodeIdentifier += 1 - cache.setNode(markerPoint) - markerCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, noduleCentre) - markerName.assignString(cache, 'aortic valve ctr' if aorticNotPulmonary else 'pulmonary valve ctr') - fm.endChange() return annotationGroups diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_heartventricles1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_heartventricles1.py index f41d8769..9c76fd9c 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_heartventricles1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_heartventricles1.py @@ -257,14 +257,12 @@ def generateBaseMesh(cls, region, options): # annotation fiducial points markerGroup = findOrCreateFieldGroup(fm, "marker") - #markerCoordinates = findOrCreateFieldCoordinates(fm, "marker_coordinates") markerName = findOrCreateFieldStoredString(fm, name="marker_name") markerLocation = findOrCreateFieldStoredMeshLocation(fm, mesh, name="marker_location") nodes = fm.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES) markerPoints = findOrCreateFieldNodeGroup(markerGroup, nodes).getNodesetGroup() markerTemplateInternal = nodes.createNodetemplate() - #markerTemplateInternal.defineField(markerCoordinates) markerTemplateInternal.defineField(markerName) markerTemplateInternal.defineField(markerLocation) @@ -1071,16 +1069,9 @@ def generateBaseMesh(cls, region, options): # apex annotation points element1 = mesh.findElementByIdentifier(1) - #markerPoint = markerPoints.createNode(nodeIdentifier, markerTemplateInternal) - #nodeIdentifier += 1 - #cache.setNode(markerPoint) - ##markerCoordinates.assignReal(cache, lvApexInnerx) - #markerName.assignString(cache, 'apex endo') - #markerLocation.assignMeshLocation(cache, element1, [ 0.0, 0.0, 0.0 ]) markerPoint = markerPoints.createNode(nodeIdentifier, markerTemplateInternal) nodeIdentifier += 1 cache.setNode(markerPoint) - #markerCoordinates.assignReal(cache, lvApexOuterx) markerName.assignString(cache, 'APEX') # interlex.org has 'apex of heart' markerLocation.assignMeshLocation(cache, element1, [ 0.0, 0.0, 1.0 ]) diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_heartventriclesbase1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_heartventriclesbase1.py index fb502221..beb39ca9 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_heartventriclesbase1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_heartventriclesbase1.py @@ -312,16 +312,11 @@ def generateBaseMesh(cls, region, options): # annotation fiducial points markerGroup = findOrCreateFieldGroup(fm, "marker") - markerCoordinates = findOrCreateFieldCoordinates(fm, "marker_coordinates") markerName = findOrCreateFieldStoredString(fm, name="marker_name") markerLocation = findOrCreateFieldStoredMeshLocation(fm, mesh, name="marker_location") nodes = fm.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES) markerPoints = findOrCreateFieldNodeGroup(markerGroup, nodes).getNodesetGroup() - markerTemplateExternal = nodes.createNodetemplate() - markerTemplateExternal.defineField(markerCoordinates) - markerTemplateExternal.defineField(markerName) - markerTemplateInternal = nodes.createNodetemplate() markerTemplateInternal.defineField(markerName) markerTemplateInternal.defineField(markerLocation) @@ -355,24 +350,6 @@ def generateBaseMesh(cls, region, options): fieldassignment = coordinates.createFieldassignment(newCoordinates) fieldassignment.setNodeset(nodes) fieldassignment.assign() - # also transform data point coordinates and/or re-evaluate from embedded locations - newCoordinates = fm.createFieldAdd(fm.createFieldMatrixMultiply(3, rotationMatrix, markerCoordinates), ventriclesOffset) - fieldassignment = markerCoordinates.createFieldassignment(newCoordinates) - fieldassignment.setNodeset(markerPoints) - fieldassignment.assign() - fieldassignment = None - newCoordinates = None - ventriclesOffset = None - fiducialHostCoordinates = fm.createFieldEmbedded(coordinates, markerLocation) - iter = markerPoints.createNodeiterator() - fiducialPoint = iter.next() - while fiducialPoint.isValid(): - cache.setNode(fiducialPoint) - result, fiducialx = fiducialHostCoordinates.evaluateReal(cache, 3) - if result == ZINC_OK: - markerCoordinates.assignReal(cache, fiducialx) - fiducialPoint = iter.next() - fiducialHostCoordinates = None # discover ventricles top LV inner, RV inner, V Outer nodes, coordinates and derivatives startLVInnerNodeId = 2 + (elementsCountUpLV - 1)*elementsCountAroundLV @@ -421,12 +398,6 @@ def generateBaseMesh(cls, region, options): zero = [ 0.0, 0.0, 0.0 ] lvOutletOuterd3 = [ None ]*elementsCountAroundOutlet - markerPoint = markerPoints.createNode(nodeIdentifier, markerTemplateExternal) - nodeIdentifier += 1 - cache.setNode(markerPoint) - markerCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, lvOutletCentre) - markerName.assignString(cache, 'aortic valve ctr') - # RV outlet points cosRvOutletLeftInclineRadians = math.cos(rvOutletLeftInclineRadians) sinRvOutletLeftInclineRadians = math.sin(rvOutletLeftInclineRadians) @@ -442,12 +413,6 @@ def generateBaseMesh(cls, region, options): vector.setMagnitude(axis1, rvOutletOuterRadius), vector.setMagnitude(axis2, rvOutletOuterRadius), elementsCountAroundOutlet) rvOutletd2 = [ vOutletElementLength*axis3[c] for c in range(3) ] - markerPoint = markerPoints.createNode(nodeIdentifier, markerTemplateExternal) - nodeIdentifier += 1 - cache.setNode(markerPoint) - markerCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, rvOutletCentre) - markerName.assignString(cache, 'pulmonary valve ctr') - # fix derivative 3 on lv outlet adjacent to rv outlet n1 = elementsCountAroundOutlet//2 lvOutletOuterd3[n1] = interp.interpolateLagrangeHermiteDerivative(lvOutletOuterx[n1], rvOutletOuterx[0], rvOutletd2, 0.0) @@ -559,18 +524,6 @@ def generateBaseMesh(cls, region, options): pd2 = interp.smoothCubicHermiteDerivativesLine([ rvInnerx[nov], ravx[0][0][noa]], [ rvInnerd2[nov], d2 ], fixStartDerivative=True, fixEndDirection=True) ravd2[0][0][noa] = pd2[1] - markerPoint = markerPoints.createNode(nodeIdentifier, markerTemplateExternal) - nodeIdentifier += 1 - cache.setNode(markerPoint) - markerCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, laCentre) - markerName.assignString(cache, 'mitral valve ctr') - - markerPoint = markerPoints.createNode(nodeIdentifier, markerTemplateExternal) - nodeIdentifier += 1 - cache.setNode(markerPoint) - markerCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, [ -laCentre[0], laCentre[1], laCentre[2] ]) - markerName.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 = elementsCountAroundRightAtriumFreeWall - 2 mag = baseHeight + baseThickness From 2121848f49d80d69a7eae274d9518160d9a205f6 Mon Sep 17 00:00:00 2001 From: Richard Christie Date: Tue, 21 Apr 2020 17:43:01 +1200 Subject: [PATCH 07/15] Make heart surface annotation groups on refinement Put more boilerplate code in Scaffold_base. --- .../meshtypes/meshtype_3d_heart1.py | 47 ++++------ .../meshtypes/meshtype_3d_heartatria1.py | 78 +++++++---------- .../meshtypes/meshtype_3d_heartventricles1.py | 81 +++++++---------- .../meshtype_3d_heartventriclesbase1.py | 53 ++---------- src/scaffoldmaker/meshtypes/scaffold_base.py | 86 +++++++++++++++---- 5 files changed, 158 insertions(+), 187 deletions(-) diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_heart1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_heart1.py index bb16d599..ee05eb30 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_heart1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_heart1.py @@ -122,7 +122,6 @@ def generateBaseMesh(cls, region, options): useCrossDerivatives = False fm = region.getFieldmodule() - fm.beginChange() coordinates = findOrCreateFieldCoordinates(fm) cache = fm.createFieldcache() @@ -409,23 +408,6 @@ def generateBaseMesh(cls, region, options): markerName.assignString(cache, "crux of heart") markerLocation.assignMeshLocation(cache, cruxElement, cruxXi) - # add to epicardium surface group - fm.defineAllFaces() - lFibrousRingGroup.addSubelements() - rFibrousRingGroup.addSubelements() - mesh2d = fm.findMeshByDimension(2) - is_exterior = fm.createFieldIsExterior() - is_exterior_face_xi3_1 = fm.createFieldAnd(is_exterior, fm.createFieldIsOnFace(Element.FACE_TYPE_XI3_1)) - is_fibrous_ring = fm.createFieldOr(lFibrousRingGroup.getFieldElementGroup(mesh2d), rFibrousRingGroup.getFieldElementGroup(mesh2d)) - is_fibrous_ring_epi = fm.createFieldAnd(is_fibrous_ring, is_exterior_face_xi3_1) - epiGroup = getAnnotationGroupForTerm(annotationGroups, get_heart_term("epicardium")) - epiGroup.getMeshGroup(mesh2d).addElementsConditional(is_fibrous_ring_epi) - del is_exterior - del is_exterior_face_xi3_1 - del is_fibrous_ring - del is_fibrous_ring_epi - - fm.endChange() return annotationGroups @classmethod @@ -461,17 +443,24 @@ def refineMesh(cls, meshrefinement, options): i += 1 @classmethod - def generateMesh(cls, region, options): + def defineFaceAnnotations(cls, region, options, annotationGroups): """ - Generate base or refined mesh. - :param region: Zinc region to create mesh in. Must be empty. + Add face annotation groups from the highest dimension mesh. + Must have defined faces and added subelements for highest dimension groups. + :param region: Zinc region containing model. :param options: Dict containing options. See getDefaultOptions(). - :return: list of AnnotationGroup for mesh. + :param annotationGroups: List of annotation groups for top-level elements. + New face annotation groups are appended to this list. """ - if not options['Refine']: - return cls.generateBaseMesh(region, options) - baseRegion = region.createRegion() - baseAnnotationGroups = cls.generateBaseMesh(baseRegion, options) - meshrefinement = MeshRefinement(baseRegion, region, baseAnnotationGroups) - cls.refineMesh(meshrefinement, options) - return meshrefinement.getAnnotationGroups() + # add epicardium of fibrous ring + fm = region.getFieldmodule() + MeshType_3d_heartventriclesbase1.defineFaceAnnotations(region, options, annotationGroups) + MeshType_3d_heartatria1.defineFaceAnnotations(region, options, annotationGroups) + lFibrousRingGroup = getAnnotationGroupForTerm(annotationGroups, get_heart_term("left fibrous ring")) + rFibrousRingGroup = getAnnotationGroupForTerm(annotationGroups, get_heart_term("right fibrous ring")) + epiGroup = getAnnotationGroupForTerm(annotationGroups, get_heart_term("epicardium")) + mesh2d = fm.findMeshByDimension(2) + is_exterior_face_xi3_1 = fm.createFieldAnd(fm.createFieldIsExterior(), fm.createFieldIsOnFace(Element.FACE_TYPE_XI3_1)) + is_non_septal_fibrous_ring = fm.createFieldXor(lFibrousRingGroup.getFieldElementGroup(mesh2d), rFibrousRingGroup.getFieldElementGroup(mesh2d)) + is_non_septal_fibrous_ring_epi = fm.createFieldAnd(is_non_septal_fibrous_ring, is_exterior_face_xi3_1) + epiGroup.getMeshGroup(mesh2d).addElementsConditional(is_non_septal_fibrous_ring_epi) diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_heartatria1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_heartatria1.py index b7a8091c..714a3672 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_heartatria1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_heartatria1.py @@ -733,7 +733,6 @@ def generateBaseMesh(cls, region, options): useCrossDerivatives = options['Use cross derivatives'] fm = region.getFieldmodule() - fm.beginChange() coordinates = findOrCreateFieldCoordinates(fm) cache = fm.createFieldcache() @@ -3146,42 +3145,9 @@ def generateBaseMesh(cls, region, options): nid1 = raTrackSurfaceFirstNodeIdentifier + e2*nodesCount1 + e1 element.setNodesByIdentifier(eft2d, [ nid1, nid1 + 1, nid1 + nodesCount1, nid1 + nodesCount1 + 1 ]) - # create endocardium and epicardium groups - fm.defineAllFaces() - laGroup.addSubelements() - raGroup.addSubelements() - aSeptumGroup.addSubelements() - mesh2d = fm.findMeshByDimension(2) - is_exterior = fm.createFieldIsExterior() - is_exterior_face_xi3_0 = fm.createFieldAnd(is_exterior, fm.createFieldIsOnFace(Element.FACE_TYPE_XI3_0)) - is_exterior_face_xi3_1 = fm.createFieldAnd(is_exterior, fm.createFieldIsOnFace(Element.FACE_TYPE_XI3_1)) - is_la = laGroup.getFieldElementGroup(mesh2d) - is_ra = raGroup.getFieldElementGroup(mesh2d) - is_la_endo = fm.createFieldAnd(is_la, is_exterior_face_xi3_0) - is_ra_endo = fm.createFieldOr(fm.createFieldAnd(fm.createFieldAnd(is_ra, is_exterior_face_xi3_0), - fm.createFieldNot(is_la_endo)), - fm.createFieldAnd(aSeptumGroup.getFieldElementGroup(mesh2d), is_exterior_face_xi3_1)) - is_a_epi = fm.createFieldAnd(fm.createFieldOr(is_la, is_ra), - fm.createFieldAnd(is_exterior_face_xi3_1, - fm.createFieldNot(aSeptumGroup.getFieldElementGroup(mesh2d)))) - epiGroup = findOrCreateAnnotationGroupForTerm(annotationGroups, region, get_heart_term("epicardium")) - epiGroup.getMeshGroup(mesh2d).addElementsConditional(is_a_epi) - laEndoGroup = findOrCreateAnnotationGroupForTerm(annotationGroups, region, get_heart_term("endocardium of left atrium")) - laEndoGroup.getMeshGroup(mesh2d).addElementsConditional(is_la_endo) - raEndoGroup = findOrCreateAnnotationGroupForTerm(annotationGroups, region, get_heart_term("endocardium of right atrium")) - raEndoGroup.getMeshGroup(mesh2d).addElementsConditional(is_ra_endo) - del is_exterior - del is_exterior_face_xi3_0 - del is_exterior_face_xi3_1 - del is_la - del is_ra - del is_la_endo - del is_ra_endo - del is_a_epi - - fm.endChange() return annotationGroups + @classmethod def refineMesh(cls, meshrefinement, options): """ @@ -3238,21 +3204,41 @@ def refineMesh(cls, meshrefinement, options): return # finish on last so can continue elsewhere element = meshrefinement._sourceElementiterator.next() + @classmethod - def generateMesh(cls, region, options): + def defineFaceAnnotations(cls, region, options, annotationGroups): """ - Generate base or refined mesh. - :param region: Zinc region to create mesh in. Must be empty. + Add face annotation groups from the highest dimension mesh. + Must have defined faces and added subelements for highest dimension groups. + :param region: Zinc region containing model. :param options: Dict containing options. See getDefaultOptions(). - :return: list of AnnotationGroup for mesh. + :param annotationGroups: List of annotation groups for top-level elements. + New face annotation groups are appended to this list. """ - if not options['Refine']: - return cls.generateBaseMesh(region, options) - baseRegion = region.createRegion() - baseAnnotationGroups = cls.generateBaseMesh(baseRegion, options) - meshrefinement = MeshRefinement(baseRegion, region, baseAnnotationGroups) - cls.refineMesh(meshrefinement, options) - return meshrefinement.getAnnotationGroups() + # create endocardium and epicardium groups + fm = region.getFieldmodule() + laGroup = getAnnotationGroupForTerm(annotationGroups, get_heart_term("left atrium myocardium")) + raGroup = getAnnotationGroupForTerm(annotationGroups, get_heart_term("right atrium myocardium")) + aSeptumGroup = getAnnotationGroupForTerm(annotationGroups, get_heart_term("interatrial septum")) + mesh2d = fm.findMeshByDimension(2) + is_exterior = fm.createFieldIsExterior() + is_exterior_face_xi3_0 = fm.createFieldAnd(is_exterior, fm.createFieldIsOnFace(Element.FACE_TYPE_XI3_0)) + is_exterior_face_xi3_1 = fm.createFieldAnd(is_exterior, fm.createFieldIsOnFace(Element.FACE_TYPE_XI3_1)) + is_la = laGroup.getFieldElementGroup(mesh2d) + is_ra = raGroup.getFieldElementGroup(mesh2d) + is_la_endo = fm.createFieldAnd(is_la, is_exterior_face_xi3_0) + is_ra_endo = fm.createFieldOr(fm.createFieldAnd(fm.createFieldAnd(is_ra, is_exterior_face_xi3_0), + fm.createFieldNot(is_la_endo)), + fm.createFieldAnd(aSeptumGroup.getFieldElementGroup(mesh2d), is_exterior_face_xi3_1)) + is_a_epi = fm.createFieldAnd(fm.createFieldOr(is_la, is_ra), + fm.createFieldAnd(is_exterior_face_xi3_1, + fm.createFieldNot(aSeptumGroup.getFieldElementGroup(mesh2d)))) + epiGroup = findOrCreateAnnotationGroupForTerm(annotationGroups, region, get_heart_term("epicardium")) + epiGroup.getMeshGroup(mesh2d).addElementsConditional(is_a_epi) + laEndoGroup = findOrCreateAnnotationGroupForTerm(annotationGroups, region, get_heart_term("endocardium of left atrium")) + laEndoGroup.getMeshGroup(mesh2d).addElementsConditional(is_la_endo) + raEndoGroup = findOrCreateAnnotationGroupForTerm(annotationGroups, region, get_heart_term("endocardium of right atrium")) + raEndoGroup.getMeshGroup(mesh2d).addElementsConditional(is_ra_endo) def getLeftAtriumPulmonaryVeinOstiaElementsCounts(elementsCountAroundLeftAtriumFreeWall, elementsCountOverAtria, commonLeftRightPvOstium): diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_heartventricles1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_heartventricles1.py index 9c76fd9c..493d8a28 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_heartventricles1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_heartventricles1.py @@ -10,7 +10,7 @@ from opencmiss.zinc.element import Element, Elementbasis, Elementfieldtemplate from opencmiss.zinc.field import Field from opencmiss.zinc.node import Node -from scaffoldmaker.annotation.annotationgroup import AnnotationGroup +from scaffoldmaker.annotation.annotationgroup import AnnotationGroup, findOrCreateAnnotationGroupForTerm, getAnnotationGroupForTerm from scaffoldmaker.annotation.heart_terms import get_heart_term from scaffoldmaker.meshtypes.scaffold_base import Scaffold_base from scaffoldmaker.utils import vector @@ -244,7 +244,6 @@ def generateBaseMesh(cls, region, options): useCrossDerivatives = options['Use cross derivatives'] fm = region.getFieldmodule() - fm.beginChange() coordinates = findOrCreateFieldCoordinates(fm) cache = fm.createFieldcache() @@ -1034,40 +1033,7 @@ def generateBaseMesh(cls, region, options): for meshGroup in meshGroups: meshGroup.addElement(element) - # create endocardium and epicardium groups - fm.defineAllFaces() - lvGroup.addSubelements() - rvGroup.addSubelements() - vSeptumGroup.addSubelements() - mesh2d = fm.findMeshByDimension(2) - is_exterior = fm.createFieldIsExterior() - is_exterior_face_xi3_0 = fm.createFieldAnd(is_exterior, fm.createFieldIsOnFace(Element.FACE_TYPE_XI3_0)) - is_exterior_face_xi3_1 = fm.createFieldAnd(is_exterior, fm.createFieldIsOnFace(Element.FACE_TYPE_XI3_1)) - is_lv = lvGroup.getFieldElementGroup(mesh2d) - is_rv = rvGroup.getFieldElementGroup(mesh2d) - is_lv_endo = fm.createFieldAnd(is_lv, is_exterior_face_xi3_0) - is_rv_endo = fm.createFieldOr(fm.createFieldAnd(fm.createFieldAnd(is_rv, is_exterior_face_xi3_0), - fm.createFieldNot(is_lv_endo)), - fm.createFieldAnd(vSeptumGroup.getFieldElementGroup(mesh2d), is_exterior_face_xi3_1)) - is_epi = fm.createFieldAnd(is_exterior_face_xi3_1, - fm.createFieldNot(vSeptumGroup.getFieldElementGroup(mesh2d))) - epiGroup = AnnotationGroup(region, get_heart_term("epicardium")) - epiGroup.getMeshGroup(mesh2d).addElementsConditional(is_epi) - lvEndoGroup = AnnotationGroup(region, get_heart_term("endocardium of left ventricle")) - lvEndoGroup.getMeshGroup(mesh2d).addElementsConditional(is_lv_endo) - rvEndoGroup = AnnotationGroup(region, get_heart_term("endocardium of right ventricle")) - rvEndoGroup.getMeshGroup(mesh2d).addElementsConditional(is_rv_endo) - del is_exterior - del is_exterior_face_xi3_0 - del is_exterior_face_xi3_1 - del is_lv - del is_rv - del is_lv_endo - del is_rv_endo - del is_epi - annotationGroups += [ lvEndoGroup, rvEndoGroup, epiGroup ] - - # apex annotation points + # apex annotation point element1 = mesh.findElementByIdentifier(1) markerPoint = markerPoints.createNode(nodeIdentifier, markerTemplateInternal) nodeIdentifier += 1 @@ -1075,7 +1041,6 @@ def generateBaseMesh(cls, region, options): markerName.assignString(cache, 'APEX') # interlex.org has 'apex of heart' markerLocation.assignMeshLocation(cache, element1, [ 0.0, 0.0, 1.0 ]) - fm.endChange() return annotationGroups @classmethod @@ -1126,21 +1091,41 @@ def refineMesh(cls, meshrefinement, options): return # finish on last so can continue in ventriclesbase element = meshrefinement._sourceElementiterator.next() + @classmethod - def generateMesh(cls, region, options): + def defineFaceAnnotations(cls, region, options, annotationGroups): """ - Generate base or refined mesh. - :param region: Zinc region to create mesh in. Must be empty. + Add face annotation groups from the highest dimension mesh. + Must have defined faces and added subelements for highest dimension groups. + :param region: Zinc region containing model. :param options: Dict containing options. See getDefaultOptions(). - :return: list of AnnotationGroup for mesh. + :param annotationGroups: List of annotation groups for top-level elements. + New face annotation groups are appended to this list. """ - if not options['Refine']: - return cls.generateBaseMesh(region, options) - baseRegion = region.createRegion() - baseAnnotationGroups = cls.generateBaseMesh(baseRegion, options) - meshrefinement = MeshRefinement(baseRegion, region, baseAnnotationGroups) - cls.refineMesh(meshrefinement, options) - return meshrefinement.getAnnotationGroups() + # create endocardium and epicardium groups + fm = region.getFieldmodule() + lvGroup = getAnnotationGroupForTerm(annotationGroups, get_heart_term("left ventricle myocardium")) + rvGroup = getAnnotationGroupForTerm(annotationGroups, get_heart_term("right ventricle myocardium")) + vSeptumGroup = getAnnotationGroupForTerm(annotationGroups, get_heart_term("interventricular septum")) + mesh2d = fm.findMeshByDimension(2) + is_exterior = fm.createFieldIsExterior() + is_exterior_face_xi3_0 = fm.createFieldAnd(is_exterior, fm.createFieldIsOnFace(Element.FACE_TYPE_XI3_0)) + is_exterior_face_xi3_1 = fm.createFieldAnd(is_exterior, fm.createFieldIsOnFace(Element.FACE_TYPE_XI3_1)) + is_lv = lvGroup.getFieldElementGroup(mesh2d) + is_rv = rvGroup.getFieldElementGroup(mesh2d) + is_lv_endo = fm.createFieldAnd(is_lv, is_exterior_face_xi3_0) + is_rv_endo = fm.createFieldOr(fm.createFieldAnd(fm.createFieldAnd(is_rv, is_exterior_face_xi3_0), + fm.createFieldNot(is_lv_endo)), + fm.createFieldAnd(vSeptumGroup.getFieldElementGroup(mesh2d), is_exterior_face_xi3_1)) + is_v_epi = fm.createFieldAnd(fm.createFieldOr(is_lv, is_rv), + fm.createFieldAnd(is_exterior_face_xi3_1, + fm.createFieldNot(vSeptumGroup.getFieldElementGroup(mesh2d)))) + epiGroup = findOrCreateAnnotationGroupForTerm(annotationGroups, region, get_heart_term("epicardium")) + epiGroup.getMeshGroup(mesh2d).addElementsConditional(is_v_epi) + lvEndoGroup = findOrCreateAnnotationGroupForTerm(annotationGroups, region, get_heart_term("endocardium of left ventricle")) + lvEndoGroup.getMeshGroup(mesh2d).addElementsConditional(is_lv_endo) + rvEndoGroup = findOrCreateAnnotationGroupForTerm(annotationGroups, region, get_heart_term("endocardium of right ventricle")) + rvEndoGroup.getMeshGroup(mesh2d).addElementsConditional(is_rv_endo) def getSeptumPoints(septumArcRadians, lvRadius, radialDisplacement, elementsCountAroundLVFreeWall, elementsCountAroundVSeptum, z, n3): diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_heartventriclesbase1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_heartventriclesbase1.py index beb39ca9..17ae0ce9 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_heartventriclesbase1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_heartventriclesbase1.py @@ -291,7 +291,6 @@ def generateBaseMesh(cls, region, options): useCrossDerivatives = False fm = region.getFieldmodule() - fm.beginChange() coordinates = findOrCreateFieldCoordinates(fm) cache = fm.createFieldcache() @@ -1338,39 +1337,6 @@ def generateBaseMesh(cls, region, options): for meshGroup in meshGroups: meshGroup.addElement(element) - # create endocardium and epicardium groups - fm.defineAllFaces() - lvGroup.addSubelements() - rvGroup.addSubelements() - vSeptumGroup.addSubelements() - mesh2d = fm.findMeshByDimension(2) - is_exterior = fm.createFieldIsExterior() - is_exterior_face_xi3_0 = fm.createFieldAnd(is_exterior, fm.createFieldIsOnFace(Element.FACE_TYPE_XI3_0)) - is_exterior_face_xi3_1 = fm.createFieldAnd(is_exterior, fm.createFieldIsOnFace(Element.FACE_TYPE_XI3_1)) - is_lv = lvGroup.getFieldElementGroup(mesh2d) - is_rv = rvGroup.getFieldElementGroup(mesh2d) - is_lv_endo = fm.createFieldAnd(is_lv, is_exterior_face_xi3_0) - is_rv_endo = fm.createFieldOr(fm.createFieldAnd(fm.createFieldAnd(is_rv, is_exterior_face_xi3_0), - fm.createFieldNot(is_lv_endo)), - fm.createFieldAnd(vSeptumGroup.getFieldElementGroup(mesh2d), is_exterior_face_xi3_1)) - is_epi = fm.createFieldAnd(is_exterior_face_xi3_1, - fm.createFieldNot(vSeptumGroup.getFieldElementGroup(mesh2d))) - epiGroup = getAnnotationGroupForTerm(annotationGroups, get_heart_term("epicardium")) - epiGroup.getMeshGroup(mesh2d).addElementsConditional(is_epi) - lvEndoGroup = getAnnotationGroupForTerm(annotationGroups, get_heart_term("endocardium of left ventricle")) - lvEndoGroup.getMeshGroup(mesh2d).addElementsConditional(is_lv_endo) - rvEndoGroup = getAnnotationGroupForTerm(annotationGroups, get_heart_term("endocardium of right ventricle")) - rvEndoGroup.getMeshGroup(mesh2d).addElementsConditional(is_rv_endo) - del is_exterior - del is_exterior_face_xi3_0 - del is_exterior_face_xi3_1 - del is_lv - del is_rv - del is_lv_endo - del is_rv_endo - del is_epi - - fm.endChange() return annotationGroups @classmethod @@ -1438,18 +1404,15 @@ def refineMesh(cls, meshrefinement, options): return # finish on last so can continue in full heart mesh element = meshrefinement._sourceElementiterator.next() + @classmethod - def generateMesh(cls, region, options): + def defineFaceAnnotations(cls, region, options, annotationGroups): """ - Generate base or refined mesh. - :param region: Zinc region to create mesh in. Must be empty. + Add face annotation groups from the highest dimension mesh. + Must have defined faces and added subelements for highest dimension groups. + :param region: Zinc region containing model. :param options: Dict containing options. See getDefaultOptions(). - :return: list of AnnotationGroup for mesh. + :param annotationGroups: List of annotation groups for top-level elements. + New face annotation groups are appended to this list. """ - if not options['Refine']: - return cls.generateBaseMesh(region, options) - baseRegion = region.createRegion() - baseAnnotationGroups = cls.generateBaseMesh(baseRegion, options) - meshrefinement = MeshRefinement(baseRegion, region, baseAnnotationGroups) - cls.refineMesh(meshrefinement, options) - return meshrefinement.getAnnotationGroups() + MeshType_3d_heartventricles1.defineFaceAnnotations(region, options, annotationGroups) diff --git a/src/scaffoldmaker/meshtypes/scaffold_base.py b/src/scaffoldmaker/meshtypes/scaffold_base.py index 33fd792c..a20e8d5c 100644 --- a/src/scaffoldmaker/meshtypes/scaffold_base.py +++ b/src/scaffoldmaker/meshtypes/scaffold_base.py @@ -2,6 +2,9 @@ Scaffold abstract base class. Describes methods each scaffold must or may override. """ +import copy +from opencmiss.utils.zinc.general import ChangeManager +from scaffoldmaker.utils.meshrefinement import MeshRefinement class Scaffold_base: ''' @@ -9,16 +12,16 @@ class Scaffold_base: Not intended to be instantiated. Most methods must be overridden by actual scaffolds. ''' - @staticmethod - def getName(): + @classmethod + def getName(cls): ''' Must override. :return: Unique type name for scaffold, for display in user interface. ''' return None - @staticmethod - def getParameterSetNames(): + @classmethod + def getParameterSetNames(cls): ''' Optionally override to return additional default parameter set names supported by getDefaultOptions(). Always have first set name 'Default'. Do not use name 'Custom' as clients may use internally. @@ -38,8 +41,8 @@ def getDefaultOptions(cls, parameterSetName='Default'): 'Boolean option' : True } - @staticmethod - def getOrderedOptionNames(): + @classmethod + def getOrderedOptionNames(cls): ''' Must override to get list of parameters in order for display and editing in user interface. Note can omit parameter names to remove from interface. @@ -87,18 +90,63 @@ def checkOptions(cls, options): options['Real option'] = 0.0 @classmethod - def generateMesh(cls, region, options): - ''' - Must override to generate scaffold mesh in region using Zinc API with options. - Scaffolds supporting refinement may switch to call other functions: - @classmethod - def generateBaseMesh(cls, region, options): - # generates base high order scaffold - @classmethod - def refineMesh(cls, meshrefinement, options): - # performs refinement + def generateBaseMesh(cls, region, options): + """ + Override to generate scaffold mesh in region using Zinc API with options. + Some older classes may do this in an override of generateMesh(). :param region: Zinc region to define model in. Must be empty. :param options: Dict containing options. See getDefaultOptions(). - :return: List of AnnotationGroup or None if not supported. - ''' - return None + :return: list of AnnotationGroup + """ + return [] + + @classmethod + def refineMesh(cls, meshrefinement, options): + """ + Override to refine/resample mesh if supported. + :param meshrefinement: MeshRefinement, which knows source and target region. + :param options: Dict containing options. See getDefaultOptions(). + """ + pass + + @classmethod + def defineFaceAnnotations(cls, region, options, annotationGroups): + """ + Override in classes with face annotation groups. + Add face annotation groups from the highest dimension mesh. + Must have defined faces and added subelements for highest dimension groups. + :param region: Zinc region containing model. + :param options: Dict containing options. See getDefaultOptions(). + :param annotationGroups: List of annotation groups for top-level elements. + New face annotation groups are appended to this list. + """ + pass + + @classmethod + def generateMesh(cls, region, options): + """ + Generate base or refined mesh. + Some classes may override to a simpler version just generating the base mesh. + :param region: Zinc region to create mesh in. Must be empty. + :param options: Dict containing options. See getDefaultOptions(). + :return: list of AnnotationGroup for mesh. + """ + fieldmodule = region.getFieldmodule() + with ChangeManager(fieldmodule): + if options.get('Refine'): + baseRegion = region.createRegion() + annotationGroups = cls.generateBaseMesh(baseRegion, options) + meshrefinement = MeshRefinement(baseRegion, region, annotationGroups) + cls.refineMesh(meshrefinement, options) + annotationGroups = meshrefinement.getAnnotationGroups() + else: + annotationGroups = cls.generateBaseMesh(region, options) + fieldmodule.defineAllFaces() + oldAnnotationGroups = copy.copy(annotationGroups) + for annotationGroup in annotationGroups: + annotationGroup.addSubelements() + cls.defineFaceAnnotations(region, options, annotationGroups) + for annotation in annotationGroups: + if annotation not in oldAnnotationGroups: + annotationGroup.addSubelements() + return annotationGroups From 9756afa97a1e627e363855a0693613489832ddb2 Mon Sep 17 00:00:00 2001 From: Richard Christie Date: Wed, 22 Apr 2020 11:07:43 +1200 Subject: [PATCH 08/15] Rename appendage to auricle, add endocardium --- src/scaffoldmaker/annotation/heart_terms.py | 6 ++++-- .../meshtypes/meshtype_3d_heartatria1.py | 14 +++++++++++--- 2 files changed, 15 insertions(+), 5 deletions(-) diff --git a/src/scaffoldmaker/annotation/heart_terms.py b/src/scaffoldmaker/annotation/heart_terms.py index fd49ed50..cd729941 100644 --- a/src/scaffoldmaker/annotation/heart_terms.py +++ b/src/scaffoldmaker/annotation/heart_terms.py @@ -24,8 +24,10 @@ ( "endocardium of right atrium", "FMA:7281", "UBERON:0009129" ), ( "interatrial septum", "FMA:7108" ), ( "fossa ovalis", "FMA:9246" ), - ( "left atrial appendage", "FMA:7219" ), # GRC rename auricle - ( "right atrial appendage", "FMA:7218" ), # GRC rename auricle + ( "left auricle", "FMA:7219", "UBERON:0006630" ), # uncertain if just the tissue like myocardium + ( "right auricle", "FMA:7218", "UBERON:0006631" ), # uncertain if just the tissue like myocardium + ( "endocardium of left auricle", "FMA:13236", "UBERON:0011006" ), + ( "endocardium of right auricle", "FMA:13235", "UBERON:0011007" ), ( "pulmonary vein", "FMA:66643", "UBERON:0002016" ), ( "left pulmonary vein", "UBERON:0009030" ), ( "left inferior pulmonary vein", "FMA:49913" ), diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_heartatria1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_heartatria1.py index 714a3672..29b583a1 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_heartatria1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_heartatria1.py @@ -740,8 +740,8 @@ def generateBaseMesh(cls, region, options): raGroup = AnnotationGroup(region, get_heart_term("right atrium myocardium")) aSeptumGroup = AnnotationGroup(region, get_heart_term("interatrial septum")) fossaGroup = AnnotationGroup(region, get_heart_term("fossa ovalis")) - laaGroup = AnnotationGroup(region, get_heart_term("left atrial appendage")) - raaGroup = AnnotationGroup(region, get_heart_term("right atrial appendage")) + laaGroup = AnnotationGroup(region, get_heart_term("left auricle")) + raaGroup = AnnotationGroup(region, get_heart_term("right auricle")) annotationGroups = [ laGroup, raGroup, aSeptumGroup, fossaGroup, laaGroup, raaGroup ] lpvOstiumSettings = lpvOstium.getScaffoldSettings() @@ -3173,7 +3173,7 @@ def refineMesh(cls, meshrefinement, options): isSeptumEdgeWedge = sourceFm.createFieldXor(sourceFm.createFieldAnd(laElementGroupField, raElementGroupField), aSeptumElementGroupField) # last atria element is last element in following group: - lastGroup = getAnnotationGroupForTerm(annotationGroups, get_heart_term("right atrial appendage")) + lastGroup = getAnnotationGroupForTerm(annotationGroups, get_heart_term("right auricle")) lastMeshGroup = lastGroup.getMeshGroup(meshrefinement._sourceMesh) lastElementIdentifier = -1 elementIter = lastMeshGroup.createElementiterator() @@ -3220,6 +3220,8 @@ def defineFaceAnnotations(cls, region, options, annotationGroups): laGroup = getAnnotationGroupForTerm(annotationGroups, get_heart_term("left atrium myocardium")) raGroup = getAnnotationGroupForTerm(annotationGroups, get_heart_term("right atrium myocardium")) aSeptumGroup = getAnnotationGroupForTerm(annotationGroups, get_heart_term("interatrial septum")) + laaGroup = getAnnotationGroupForTerm(annotationGroups, get_heart_term("left auricle")) + raaGroup = getAnnotationGroupForTerm(annotationGroups, get_heart_term("right auricle")) mesh2d = fm.findMeshByDimension(2) is_exterior = fm.createFieldIsExterior() is_exterior_face_xi3_0 = fm.createFieldAnd(is_exterior, fm.createFieldIsOnFace(Element.FACE_TYPE_XI3_0)) @@ -3233,12 +3235,18 @@ def defineFaceAnnotations(cls, region, options, annotationGroups): is_a_epi = fm.createFieldAnd(fm.createFieldOr(is_la, is_ra), fm.createFieldAnd(is_exterior_face_xi3_1, fm.createFieldNot(aSeptumGroup.getFieldElementGroup(mesh2d)))) + is_laa_endo = fm.createFieldAnd(laaGroup.getFieldElementGroup(mesh2d), is_exterior_face_xi3_0) + is_raa_endo = fm.createFieldAnd(raaGroup.getFieldElementGroup(mesh2d), is_exterior_face_xi3_0) epiGroup = findOrCreateAnnotationGroupForTerm(annotationGroups, region, get_heart_term("epicardium")) epiGroup.getMeshGroup(mesh2d).addElementsConditional(is_a_epi) laEndoGroup = findOrCreateAnnotationGroupForTerm(annotationGroups, region, get_heart_term("endocardium of left atrium")) laEndoGroup.getMeshGroup(mesh2d).addElementsConditional(is_la_endo) raEndoGroup = findOrCreateAnnotationGroupForTerm(annotationGroups, region, get_heart_term("endocardium of right atrium")) raEndoGroup.getMeshGroup(mesh2d).addElementsConditional(is_ra_endo) + laaEndoGroup = findOrCreateAnnotationGroupForTerm(annotationGroups, region, get_heart_term("endocardium of left auricle")) + laaEndoGroup.getMeshGroup(mesh2d).addElementsConditional(is_laa_endo) + raaEndoGroup = findOrCreateAnnotationGroupForTerm(annotationGroups, region, get_heart_term("endocardium of right auricle")) + raaEndoGroup.getMeshGroup(mesh2d).addElementsConditional(is_raa_endo) def getLeftAtriumPulmonaryVeinOstiaElementsCounts(elementsCountAroundLeftAtriumFreeWall, elementsCountOverAtria, commonLeftRightPvOstium): From b489446fdcd4a430b4498496f40fb8916d8c50b5 Mon Sep 17 00:00:00 2001 From: Richard Christie Date: Wed, 22 Apr 2020 13:21:05 +1200 Subject: [PATCH 09/15] Redefine right auricle to just the protrusion --- .../meshtypes/meshtype_3d_heartatria1.py | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_heartatria1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_heartatria1.py index 29b583a1..c0c3d459 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_heartatria1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_heartatria1.py @@ -2142,8 +2142,9 @@ def generateBaseMesh(cls, region, options): continue # right atrial appendage scalefactors = None meshGroups = [ raMeshGroup ] - if (e1 >= elementsCountAroundRightAtriumPosteriorVenous) and (e1 < elementsCountAroundRightAtriumFreeWall - elementsCountAroundRightAtriumAorta - 1): - meshGroups += [ raaMeshGroup ] + # Anderson definition of right atrial appendage starts at crista terminalis: + #if (e1 >= elementsCountAroundRightAtriumPosteriorVenous) and (e1 < elementsCountAroundRightAtriumFreeWall - elementsCountAroundRightAtriumAorta - 1): + # meshGroups += [ raaMeshGroup ] if e1 == -1: # crux/posterior interatrial groove straddles left and right atria, collapsed to 6 node wedge nids[0] = labNodeId[0][elementsCountAroundLeftAtriumFreeWall] @@ -2933,8 +2934,10 @@ def generateBaseMesh(cls, region, options): elementsCountRadial = elementsCountLaaRadial, meshGroups = [ laMeshGroup, laaMeshGroup ]) - # create right atrial appendage 'plain' elements - meshGroups = [ raMeshGroup, raaMeshGroup ] + # create right atrium plain elements + # Anderson considers these part of the right atrial appendage: + #meshGroups = [ raMeshGroup, raaMeshGroup ] + meshGroups = [ raMeshGroup ] for e2 in range(2): for e1 in range(elementsCountOverSideRightAtriumPouch): eft1 = eft From f966d9a4ef32c46bb8870a678e98fd70f5cbabb2 Mon Sep 17 00:00:00 2001 From: Richard Christie Date: Wed, 22 Apr 2020 13:58:32 +1200 Subject: [PATCH 10/15] Add heart marker SVC-RA --- .../meshtypes/meshtype_3d_heartatria1.py | 30 +++++++++++++++---- 1 file changed, 25 insertions(+), 5 deletions(-) diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_heartatria1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_heartatria1.py index c0c3d459..8d406320 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_heartatria1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_heartatria1.py @@ -6,7 +6,8 @@ from __future__ import division import copy import math -from opencmiss.utils.zinc.field import findOrCreateFieldCoordinates +from opencmiss.utils.zinc.field import findOrCreateFieldCoordinates, findOrCreateFieldGroup, \ + findOrCreateFieldNodeGroup, findOrCreateFieldStoredMeshLocation, findOrCreateFieldStoredString from opencmiss.utils.zinc.finiteelement import getMaximumElementIdentifier, getMaximumNodeIdentifier from opencmiss.zinc.element import Element, Elementbasis from opencmiss.zinc.field import Field @@ -733,6 +734,8 @@ def generateBaseMesh(cls, region, options): useCrossDerivatives = options['Use cross derivatives'] fm = region.getFieldmodule() + mesh = fm.findMeshByDimension(3) + nodes = fm.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES) coordinates = findOrCreateFieldCoordinates(fm) cache = fm.createFieldcache() @@ -777,12 +780,20 @@ def generateBaseMesh(cls, region, options): lFibrousRingGroup = AnnotationGroup(region, get_heart_term("left fibrous ring")) rFibrousRingGroup = AnnotationGroup(region, get_heart_term("right fibrous ring")) + # annotation fiducial points + markerGroup = findOrCreateFieldGroup(fm, "marker") + markerName = findOrCreateFieldStoredString(fm, name="marker_name") + markerLocation = findOrCreateFieldStoredMeshLocation(fm, mesh, name="marker_location") + + markerPoints = findOrCreateFieldNodeGroup(markerGroup, nodes).getNodesetGroup() + markerTemplateInternal = nodes.createNodetemplate() + markerTemplateInternal.defineField(markerName) + markerTemplateInternal.defineField(markerLocation) + ############## # Create nodes ############## - nodes = fm.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES) - nodetemplate = nodes.createNodetemplate() nodetemplate.defineField(coordinates) nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_VALUE, 1) @@ -798,8 +809,6 @@ def generateBaseMesh(cls, region, options): nodeIdentifier = max(1, getMaximumNodeIdentifier(nodes) + 1) - mesh = fm.findMeshByDimension(3) - elementIdentifier = max(1, getMaximumElementIdentifier(mesh) + 1) laMeshGroup = laGroup.getMeshGroup(mesh) @@ -2943,6 +2952,8 @@ def generateBaseMesh(cls, region, options): eft1 = eft elementtemplate1 = elementtemplate scalefactors = None + addMarker = None + nc = 2 + e1 if e2 == 0: nids = [ ractNodeId[0][nc], ractNodeId[0][nc + 1], raaqNodeId[0][e1], raaqNodeId[0][e1 + 1], @@ -2963,6 +2974,7 @@ def generateBaseMesh(cls, region, options): remapEftNodeValueLabel(eft1, [ 2, 6 ], Node.VALUE_LABEL_D_DS1, [ ( Node.VALUE_LABEL_D_DS1, [] ), ( Node.VALUE_LABEL_D_DS2, [] ) ]) remapEftNodeValueLabel(eft1, [ 3, 7 ], Node.VALUE_LABEL_D_DS1, [ ( Node.VALUE_LABEL_D_DS1, [1] ), ( Node.VALUE_LABEL_D_DS2, [] ) ]) remapEftNodeValueLabel(eft1, [ 3, 7 ], Node.VALUE_LABEL_D_DS2, [ ( Node.VALUE_LABEL_D_DS2, [1] ) ]) + addMarker = { "name" : "SVC-RA", "xi" : [ 0.0, 1.0, 0.0 ] } elif (e2 == 1) and (e1 < (elementsCountOverSideRightAtriumPouch - 1)): pass # regular elements else: @@ -2992,6 +3004,14 @@ def generateBaseMesh(cls, region, options): result3 = '-' #print('create element raa plain', element.isValid(), elementIdentifier, result2, result3, nids) elementIdentifier += 1 + + if addMarker: + markerPoint = markerPoints.createNode(nodeIdentifier, markerTemplateInternal) + nodeIdentifier += 1 + cache.setNode(markerPoint) + markerName.assignString(cache, addMarker["name"]) + markerLocation.assignMeshLocation(cache, element, addMarker["xi"]) + for meshGroup in meshGroups: meshGroup.addElement(element) From 45a51634ceb4cb998c68226d4451886c1394d8a9 Mon Sep 17 00:00:00 2001 From: Richard Christie Date: Thu, 23 Apr 2020 11:20:35 +1200 Subject: [PATCH 11/15] Fix linear map for hermite cross derivative terms --- src/scaffoldmaker/utils/annulusmesh.py | 6 +-- .../utils/eftfactory_tricubichermite.py | 51 +++++++++++++++++++ 2 files changed, 53 insertions(+), 4 deletions(-) diff --git a/src/scaffoldmaker/utils/annulusmesh.py b/src/scaffoldmaker/utils/annulusmesh.py index ca918d87..0daed94b 100644 --- a/src/scaffoldmaker/utils/annulusmesh.py +++ b/src/scaffoldmaker/utils/annulusmesh.py @@ -343,8 +343,7 @@ def createAnnulusMesh3d(nodes, mesh, nextNodeIdentifier, nextElementIdentifier, eft1 = eftFactory.createEftNoCrossDerivatives() setEftScaleFactorIds(eft1, [1], []) if mapStartLinearDerivativeXi3: - eftFactory.setEftLinearDerivative(eft1, [ 1, 5 ], Node.VALUE_LABEL_D_DS3, 1, 5, 1) - eftFactory.setEftLinearDerivative(eft1, [ 2, 6 ], Node.VALUE_LABEL_D_DS3, 2, 6, 1) + eftFactory.setEftLinearDerivative2(eft1, [ 1, 5, 2, 6 ], Node.VALUE_LABEL_D_DS3, [ Node.VALUE_LABEL_D2_DS1DS3 ]) if mapStartDerivatives: for i in range(2): lns = [ 1, 5 ] if (i == 0) else [ 2, 6 ] @@ -370,8 +369,7 @@ def createAnnulusMesh3d(nodes, mesh, nextNodeIdentifier, nextElementIdentifier, remapEftNodeValueLabel(eft1, ln, Node.VALUE_LABEL_D2_DS2DS3, \ derivativeSignsToExpressionTerms( ( Node.VALUE_LABEL_D_DS1, Node.VALUE_LABEL_D_DS2, Node.VALUE_LABEL_D_DS3 ), d3Map)) if mapEndLinearDerivativeXi3: - eftFactory.setEftLinearDerivative(eft1, [ 3, 7 ], Node.VALUE_LABEL_D_DS3, 3, 7, 1) - eftFactory.setEftLinearDerivative(eft1, [ 4, 8 ], Node.VALUE_LABEL_D_DS3, 4, 8, 1) + eftFactory.setEftLinearDerivative2(eft1, [ 3, 7, 4, 8 ], Node.VALUE_LABEL_D_DS3, [ Node.VALUE_LABEL_D2_DS1DS3 ]) if mapEndDerivatives: for i in range(2): lns = [ 3, 7 ] if (i == 0) else [ 4, 8 ] diff --git a/src/scaffoldmaker/utils/eftfactory_tricubichermite.py b/src/scaffoldmaker/utils/eftfactory_tricubichermite.py index be3b1f1e..0934420a 100644 --- a/src/scaffoldmaker/utils/eftfactory_tricubichermite.py +++ b/src/scaffoldmaker/utils/eftfactory_tricubichermite.py @@ -833,6 +833,57 @@ def setEftLinearDerivative(self, eft, basisNodes, derivative, localNode1, localN eft.setTermNodeParameter(f, 2, localNode1, Node.VALUE_LABEL_VALUE, 1) eft.setTermScaling(f, 2, [minus1scaleFactorIndex]) + def setEftLinearDerivative2(self, eft, basisNodes, derivative, crossDerivatives, localNodes=None, minus1scaleFactorIndex=1): + ''' + Makes the derivative at basisNodes equal to the difference in the value parameter at the corresponding + localNodes to make the basis equivalent to linear Lagrange. + Specified cross derivatives are also set to equal linear differences in the other derivative. + :basisNodes: List of pairs of basis nodes to set derivative at, numbering from 1. + e.g. [ 1, 5, 2, 6 ] will set derivative at 1 and 5 to value[5] - value[1], and 2 and 6 to value[6] - value[2] + :derivative: Derivative to set from difference in value e.g. Node.VALUE_LABEL_D_DS3. + :crossDerivatives: List of cross derivatives to set from difference in other derivative. + Valid constants are mixed second cross derivatives including the specified derivative above, + e.g. [ Node.VALUE_LABEL_D2_DS1DS3 ] would be set to difference of Node.VALUE_LABEL_D_DS1 for derivative Node.VALUE_LABEL_D_DS3 + :localNodes: List of pairs of local nodes to get values and derivatives for, numbering from 1. If None, use basisNodes. + :param minus1scaleFactorIndex: Local scale factor index for general value -1.0, default=1. + ''' + if derivative == Node.VALUE_LABEL_D_DS3: + allowedSourceDerivatives = [ Node.VALUE_LABEL_D_DS1, Node.VALUE_LABEL_D_DS2 ] + allowedCrossDerivatives = [ Node.VALUE_LABEL_D2_DS1DS3, Node.VALUE_LABEL_D2_DS2DS3 ] + elif derivative == Node.VALUE_LABEL_D_DS2: + allowedSourceDerivatives = [ Node.VALUE_LABEL_D_DS1, Node.VALUE_LABEL_D_DS3 ] + allowedCrossDerivatives = [ Node.VALUE_LABEL_D2_DS1DS2, Node.VALUE_LABEL_D2_DS2DS3 ] + elif derivative == Node.VALUE_LABEL_D_DS1: + allowedSourceDerivatives = [ Node.VALUE_LABEL_D_DS2, Node.VALUE_LABEL_D_DS3 ] + allowedCrossDerivatives = [ Node.VALUE_LABEL_D2_DS1DS2, Node.VALUE_LABEL_D2_DS1DS3 ] + else: + assert False, 'setEftLinearDerivative2. Invalid derivative' + sourceValueType = [ Node.VALUE_LABEL_VALUE ] + targetValueType = [ derivative ] + for crossDerivative in crossDerivatives: + for i in range(2): + if crossDerivative == allowedCrossDerivatives[i]: + sourceValueType.append(allowedSourceDerivatives[i]) + targetValueType.append(allowedCrossDerivatives[i]) + break + else: + assert False, 'setEftLinearDerivative2. Invalid cross derivative' + assert len(basisNodes) in [ 2, 4, 6, 8 ] + if localNodes: + assert len(localNodes) == len(basisNodes) + else: + localNodes = basisNodes + for i in range(len(sourceValueType)): + nodeFunctionIndex = targetValueType[i] - Node.VALUE_LABEL_VALUE + 1 + for bi in range(len(basisNodes)): + li = bi - (bi % 2) + f = (basisNodes[bi] - 1)*8 + nodeFunctionIndex + eft.setFunctionNumberOfTerms(f, 2) + eft.setTermNodeParameter(f, 1, localNodes[li + 1], sourceValueType[i], 1) + eft.setTermScaling(f, 1, []) + eft.setTermNodeParameter(f, 2, localNodes[li ], sourceValueType[i], 1) + eft.setTermScaling(f, 2, [minus1scaleFactorIndex]) + def setEftMidsideXi1HangingNode(self, eft, hangingBasisNode, otherBasisNode, localNode1, localNode2, scaleFactorIndexes): ''' Makes the functions for eft at basisNode work as a hanging node interpolating the parameters From 9000ba6d094bf2247a315badab6077fd48aaba25 Mon Sep 17 00:00:00 2001 From: Richard Christie Date: Thu, 23 Apr 2020 11:27:13 +1200 Subject: [PATCH 12/15] Export surface annotation groups to VTK --- .../annotation/annotationgroup.py | 16 +++++ src/scaffoldmaker/utils/exportvtk.py | 67 +++++++++++++++---- 2 files changed, 69 insertions(+), 14 deletions(-) diff --git a/src/scaffoldmaker/annotation/annotationgroup.py b/src/scaffoldmaker/annotation/annotationgroup.py index cfa48004..aab2ede3 100644 --- a/src/scaffoldmaker/annotation/annotationgroup.py +++ b/src/scaffoldmaker/annotation/annotationgroup.py @@ -75,6 +75,14 @@ def getMeshGroup(self, mesh): ''' return self.getFieldElementGroup(mesh).getMeshGroup() + def hasMeshGroup(self, mesh): + ''' + :param mesh: The Zinc mesh to query a sub group of. + :return: True if MeshGroup for mesh exists and is not empty, otherwise False. + ''' + elementGroup = self._group.getFieldElementGroup(mesh) + return elementGroup.isValid() and (elementGroup.getMeshGroup().getSize() > 0) + def getNodesetGroup(self, nodeset): ''' :param nodeset: The Zinc nodeset to manage a sub group of. @@ -82,6 +90,14 @@ def getNodesetGroup(self, nodeset): ''' return self.getFieldNodeGroup(nodeset).getNodesetGroup() + def hasNodesetGroup(self, nodeset): + ''' + :param nodeset: The Zinc nodeset to query a sub group of. + :return: True if NodesetGroup for nodeset exists and is not empty, otherwise False. + ''' + nodeGroup = self._group.getFieldNodeGroup(nodeset) + return nodeGroup.isValid() and (nodeGroup.getNodesetGroup().getSize() > 0) + def addSubelements(self): ''' Call after group is complete and faces have been defined to add faces and diff --git a/src/scaffoldmaker/utils/exportvtk.py b/src/scaffoldmaker/utils/exportvtk.py index a72d80f1..2d4cd02d 100644 --- a/src/scaffoldmaker/utils/exportvtk.py +++ b/src/scaffoldmaker/utils/exportvtk.py @@ -7,7 +7,7 @@ from opencmiss.zinc.field import Field from opencmiss.utils.zinc.field import findOrCreateFieldCoordinates from opencmiss.utils.zinc.finiteelement import getElementNodeIdentifiersBasisOrder - +from opencmiss.zinc.result import RESULT_OK class ExportVtk: ''' @@ -45,23 +45,37 @@ def _write(self, outstream): cache = self._fieldmodule.createFieldcache() nodes = self._fieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES) + # exclude marker point nodes from output + markerNodes = None + markerGroup = self._fieldmodule.findFieldByName("marker") + if markerGroup.isValid(): + markerGroup = markerGroup.castGroup() + markerNodeGroup = markerGroup.getFieldNodeGroup(nodes) + if markerNodeGroup.isValid(): + markerNodes = markerNodeGroup.getNodesetGroup() pointCount = nodes.getSize() + if markerNodes: + pointCount -= markerNodes.getSize() outstream.write('POINTS ' + str(pointCount) + ' double\n') nodeIter = nodes.createNodeiterator() node = nodeIter.next() index = 0 while node.isValid(): - nodeIdentifierToIndex[node.getIdentifier()] = index - cache.setNode(node) - result, x = self._coordinates.evaluateReal(cache, coordinatesCount) - first = True - for s in x: - if not first: - outstream.write(' ') - outstream.write(str(s)) - first = False - outstream.write('\n') - index += 1 + if not (markerNodes and markerNodes.containsNode(node)): + nodeIdentifierToIndex[node.getIdentifier()] = index + cache.setNode(node) + result, x = self._coordinates.evaluateReal(cache, coordinatesCount) + if result != RESULT_OK: + print("Coordinates not found for node", node.getIdentifier()) + x = [ 0.0, 0.0, 0.0 ] + first = True + for s in x: + if not first: + outstream.write(' ') + outstream.write(str(s)) + first = False + outstream.write('\n') + index += 1 node = nodeIter.next() # following assumes all hex (3-D) or all quad (2-D) elements @@ -93,9 +107,19 @@ def _write(self, outstream): outstream.write(cellTypeString + ' ') outstream.write(cellTypeString + '\n') - if self._annotationGroups: + # use cell data for annotation groups containing elements of mesh dimension + # use point data for lower dimensional annotation groups + pointAnnotationGroups = [] + cellAnnotationGroups = [] + for annotationGroup in self._annotationGroups: + if annotationGroup.hasMeshGroup(self._mesh): + cellAnnotationGroups.append(annotationGroup) + elif annotationGroup.hasNodesetGroup(nodes): + pointAnnotationGroups.append(annotationGroup) + + if cellAnnotationGroups: outstream.write('CELL_DATA ' + str(cellCount) + '\n') - for annotationGroup in self._annotationGroups: + for annotationGroup in cellAnnotationGroups: safeName = annotationGroup.getName().replace(' ', '_') outstream.write('SCALARS ' + safeName + ' int 1\n') outstream.write('LOOKUP_TABLE default\n') @@ -107,6 +131,21 @@ def _write(self, outstream): element = elementIter.next() outstream.write('\n') + if pointAnnotationGroups: + outstream.write('POINT_DATA ' + str(pointCount) + '\n') + for annotationGroup in pointAnnotationGroups: + safeName = annotationGroup.getName().replace(' ', '_') + outstream.write('SCALARS ' + safeName + ' int 1\n') + outstream.write('LOOKUP_TABLE default\n') + nodesetGroup = annotationGroup.getNodesetGroup(nodes) + nodeIter = nodes.createNodeiterator() + node = nodeIter.next() + while node.isValid(): + if not (markerNodes and markerNodes.containsNode(node)): + outstream.write('1 ' if nodesetGroup.containsNode(node) else '0 ') + node = nodeIter.next() + outstream.write('\n') + def writeFile(self, filename): ''' From d33c6a6c9e3e7da34050ae342bedc6e31629e2f0 Mon Sep 17 00:00:00 2001 From: Richard Christie Date: Thu, 23 Apr 2020 14:14:49 +1200 Subject: [PATCH 13/15] Write fiducial markers to csv file with VTK export --- src/scaffoldmaker/utils/exportvtk.py | 79 +++++++++++++++++----------- 1 file changed, 49 insertions(+), 30 deletions(-) diff --git a/src/scaffoldmaker/utils/exportvtk.py b/src/scaffoldmaker/utils/exportvtk.py index 2d4cd02d..e68e2607 100644 --- a/src/scaffoldmaker/utils/exportvtk.py +++ b/src/scaffoldmaker/utils/exportvtk.py @@ -3,10 +3,12 @@ ''' import io +import os from sys import version_info -from opencmiss.zinc.field import Field from opencmiss.utils.zinc.field import findOrCreateFieldCoordinates from opencmiss.utils.zinc.finiteelement import getElementNodeIdentifiersBasisOrder +from opencmiss.utils.zinc.general import ChangeManager +from opencmiss.zinc.field import Field from opencmiss.zinc.result import RESULT_OK class ExportVtk: @@ -27,9 +29,17 @@ def __init__(self, region, description, annotationGroups = None): self._mesh = self._fieldmodule.findMeshByDimension(dimension) if self._mesh.getSize() > 0: break + self._nodes = self._fieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES) self._coordinates = findOrCreateFieldCoordinates(self._fieldmodule) self._description = description self._annotationGroups = annotationGroups if annotationGroups else [] + self._markerNodes = None + markerGroup = self._fieldmodule.findFieldByName("marker") + if markerGroup.isValid(): + markerGroup = markerGroup.castGroup() + markerNodeGroup = markerGroup.getFieldNodeGroup(self._nodes) + if markerNodeGroup.isValid(): + self._markerNodes = markerNodeGroup.getNodesetGroup() def _write(self, outstream): @@ -40,41 +50,26 @@ def _write(self, outstream): outstream.write('ASCII\n') outstream.write('DATASET UNSTRUCTURED_GRID\n') nodeIdentifierToIndex = {} # map needed since vtk points are zero index based, i.e. have no identifier - coordinates = self._coordinates - coordinatesCount = coordinates.getNumberOfComponents() + coordinatesCount = self._coordinates.getNumberOfComponents() cache = self._fieldmodule.createFieldcache() - nodes = self._fieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES) - # exclude marker point nodes from output - markerNodes = None - markerGroup = self._fieldmodule.findFieldByName("marker") - if markerGroup.isValid(): - markerGroup = markerGroup.castGroup() - markerNodeGroup = markerGroup.getFieldNodeGroup(nodes) - if markerNodeGroup.isValid(): - markerNodes = markerNodeGroup.getNodesetGroup() - pointCount = nodes.getSize() - if markerNodes: - pointCount -= markerNodes.getSize() + # exclude marker nodes from output + pointCount = self._nodes.getSize() + if self._markerNodes: + pointCount -= self._markerNodes.getSize() outstream.write('POINTS ' + str(pointCount) + ' double\n') - nodeIter = nodes.createNodeiterator() + nodeIter = self._nodes.createNodeiterator() node = nodeIter.next() index = 0 while node.isValid(): - if not (markerNodes and markerNodes.containsNode(node)): + if not (self._markerNodes and self._markerNodes.containsNode(node)): nodeIdentifierToIndex[node.getIdentifier()] = index cache.setNode(node) result, x = self._coordinates.evaluateReal(cache, coordinatesCount) if result != RESULT_OK: print("Coordinates not found for node", node.getIdentifier()) - x = [ 0.0, 0.0, 0.0 ] - first = True - for s in x: - if not first: - outstream.write(' ') - outstream.write(str(s)) - first = False - outstream.write('\n') + x = [ 0.0 ]*coordinatesCount + outstream.write(" ".join(str(s) for s in x) + "\n") index += 1 node = nodeIter.next() @@ -94,7 +89,7 @@ def _write(self, outstream): elementIter = self._mesh.createElementiterator() element = elementIter.next() while element.isValid(): - eft = element.getElementfieldtemplate(coordinates, -1) # assumes all components same + eft = element.getElementfieldtemplate(self._coordinates, -1) # assumes all components same nodeIdentifiers = getElementNodeIdentifiersBasisOrder(element, eft) outstream.write(localNodeCountStr) for localIndex in vtkIndexing: @@ -114,7 +109,7 @@ def _write(self, outstream): for annotationGroup in self._annotationGroups: if annotationGroup.hasMeshGroup(self._mesh): cellAnnotationGroups.append(annotationGroup) - elif annotationGroup.hasNodesetGroup(nodes): + elif annotationGroup.hasNodesetGroup(self._nodes): pointAnnotationGroups.append(annotationGroup) if cellAnnotationGroups: @@ -137,19 +132,43 @@ def _write(self, outstream): safeName = annotationGroup.getName().replace(' ', '_') outstream.write('SCALARS ' + safeName + ' int 1\n') outstream.write('LOOKUP_TABLE default\n') - nodesetGroup = annotationGroup.getNodesetGroup(nodes) - nodeIter = nodes.createNodeiterator() + nodesetGroup = annotationGroup.getNodesetGroup(self._nodes) + nodeIter = self._nodes.createNodeiterator() node = nodeIter.next() while node.isValid(): - if not (markerNodes and markerNodes.containsNode(node)): + if not (self._markerNodes and self._markerNodes.containsNode(node)): outstream.write('1 ' if nodesetGroup.containsNode(node) else '0 ') node = nodeIter.next() outstream.write('\n') + def _writeMarkers(self, outstream): + coordinatesCount = self._coordinates.getNumberOfComponents() + cache = self._fieldmodule.createFieldcache() + markerLocation = self._fieldmodule.findFieldByName("marker_location") + markerName = self._fieldmodule.findFieldByName("marker_name") + if markerLocation.isValid() and markerName.isValid(): + with ChangeManager(self._fieldmodule): + markerCoordinates = self._fieldmodule.createFieldEmbedded(self._coordinates, markerLocation) + nodeIter = self._markerNodes.createNodeiterator() + node = nodeIter.next() + while node.isValid(): + cache.setNode(node) + result, x = markerCoordinates.evaluateReal(cache, coordinatesCount) + if result == RESULT_OK: + name = markerName.evaluateString(cache) + outstream.write(",".join(str(s) for s in x) + "," + name + "\n") + node = nodeIter.next() + del markerCoordinates + + def writeFile(self, filename): ''' Export to legacy vtk file/ ''' with open(filename, 'w') as outstream: self._write(outstream) + if self._markerNodes and (self._markerNodes.getSize() > 0): + markerFilename = os.path.splitext(filename)[0] + "_marker.csv" + with open(markerFilename, 'w') as outstream: + self._writeMarkers(outstream) From d2370c380ea452b3fc25e91247a63907b714f804 Mon Sep 17 00:00:00 2001 From: Richard Christie Date: Fri, 24 Apr 2020 17:03:26 +1200 Subject: [PATCH 14/15] Review fixes, add epicardium of auricles --- src/scaffoldmaker/annotation/annotationgroup.py | 2 +- src/scaffoldmaker/annotation/heart_terms.py | 2 ++ src/scaffoldmaker/meshtypes/meshtype_3d_heart2.py | 7 ++++--- src/scaffoldmaker/meshtypes/meshtype_3d_heartatria1.py | 6 ++++++ 4 files changed, 13 insertions(+), 4 deletions(-) diff --git a/src/scaffoldmaker/annotation/annotationgroup.py b/src/scaffoldmaker/annotation/annotationgroup.py index aab2ede3..f8e95247 100644 --- a/src/scaffoldmaker/annotation/annotationgroup.py +++ b/src/scaffoldmaker/annotation/annotationgroup.py @@ -129,7 +129,7 @@ def findAnnotationGroupByName(annotationGroups: list, name: str): def findOrCreateAnnotationGroupForTerm(annotationGroups: list, region, term) -> AnnotationGroup: ''' - Find existing annotation group for term, or create it for region if not gound. + Find existing annotation group for term, or create it for region if not found. If annotation group created here, append it to annotationGroups. :param annotationGroups: list(AnnotationGroup) :param region: Zinc region to create group for. diff --git a/src/scaffoldmaker/annotation/heart_terms.py b/src/scaffoldmaker/annotation/heart_terms.py index cd729941..f96194dd 100644 --- a/src/scaffoldmaker/annotation/heart_terms.py +++ b/src/scaffoldmaker/annotation/heart_terms.py @@ -28,6 +28,8 @@ ( "right auricle", "FMA:7218", "UBERON:0006631" ), # uncertain if just the tissue like myocardium ( "endocardium of left auricle", "FMA:13236", "UBERON:0011006" ), ( "endocardium of right auricle", "FMA:13235", "UBERON:0011007" ), + ( "epicardium of left auricle", "FMA:13233" ), + ( "epicardium of right auricle", "FMA:13232" ), ( "pulmonary vein", "FMA:66643", "UBERON:0002016" ), ( "left pulmonary vein", "UBERON:0009030" ), ( "left inferior pulmonary vein", "FMA:49913" ), diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_heart2.py b/src/scaffoldmaker/meshtypes/meshtype_3d_heart2.py index 7dfd9ae5..5923b01a 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_heart2.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_heart2.py @@ -5,7 +5,7 @@ from __future__ import division import math from opencmiss.utils.zinc.field import findOrCreateFieldCoordinates -from scaffoldmaker.annotation.annotationgroup import AnnotationGroup +from scaffoldmaker.annotation.annotationgroup import AnnotationGroup, mergeAnnotationGroups from scaffoldmaker.meshtypes.meshtype_3d_heartatria2 import MeshType_3d_heartatria2 from scaffoldmaker.meshtypes.meshtype_3d_heartventriclesbase2 import MeshType_3d_heartventriclesbase2 from scaffoldmaker.meshtypes.scaffold_base import Scaffold_base @@ -82,8 +82,9 @@ def generateBaseMesh(cls, region, options): cache = fm.createFieldcache() # generate heartventriclesbase2 model and put atria2 on it - annotationGroups = MeshType_3d_heartventriclesbase2.generateBaseMesh(region, options) - annotationGroups += MeshType_3d_heartatria2.generateBaseMesh(region, options) + ventriclesAnnotationGroups = MeshType_3d_heartventriclesbase2.generateBaseMesh(region, options) + atriaAnnotationGroups = MeshType_3d_heartatria2.generateBaseMesh(region, options) + annotationGroups = mergeAnnotationGroups(ventriclesAnnotationGroups, atriaAnnotationGroups) fm.endChange() return annotationGroups diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_heartatria1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_heartatria1.py index 8d406320..93c7e140 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_heartatria1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_heartatria1.py @@ -3260,6 +3260,8 @@ def defineFaceAnnotations(cls, region, options, annotationGroups): fm.createFieldNot(aSeptumGroup.getFieldElementGroup(mesh2d)))) is_laa_endo = fm.createFieldAnd(laaGroup.getFieldElementGroup(mesh2d), is_exterior_face_xi3_0) is_raa_endo = fm.createFieldAnd(raaGroup.getFieldElementGroup(mesh2d), is_exterior_face_xi3_0) + is_laa_epi = fm.createFieldAnd(laaGroup.getFieldElementGroup(mesh2d), is_exterior_face_xi3_1) + is_raa_epi = fm.createFieldAnd(raaGroup.getFieldElementGroup(mesh2d), is_exterior_face_xi3_1) epiGroup = findOrCreateAnnotationGroupForTerm(annotationGroups, region, get_heart_term("epicardium")) epiGroup.getMeshGroup(mesh2d).addElementsConditional(is_a_epi) laEndoGroup = findOrCreateAnnotationGroupForTerm(annotationGroups, region, get_heart_term("endocardium of left atrium")) @@ -3270,6 +3272,10 @@ def defineFaceAnnotations(cls, region, options, annotationGroups): laaEndoGroup.getMeshGroup(mesh2d).addElementsConditional(is_laa_endo) raaEndoGroup = findOrCreateAnnotationGroupForTerm(annotationGroups, region, get_heart_term("endocardium of right auricle")) raaEndoGroup.getMeshGroup(mesh2d).addElementsConditional(is_raa_endo) + laaEpiGroup = findOrCreateAnnotationGroupForTerm(annotationGroups, region, get_heart_term("epicardium of left auricle")) + laaEpiGroup.getMeshGroup(mesh2d).addElementsConditional(is_laa_epi) + raaEpiGroup = findOrCreateAnnotationGroupForTerm(annotationGroups, region, get_heart_term("epicardium of right auricle")) + raaEpiGroup.getMeshGroup(mesh2d).addElementsConditional(is_raa_epi) def getLeftAtriumPulmonaryVeinOstiaElementsCounts(elementsCountAroundLeftAtriumFreeWall, elementsCountOverAtria, commonLeftRightPvOstium): From 7bc9f347febf893e9566e37ec4aa0e962c12491d Mon Sep 17 00:00:00 2001 From: Richard Christie Date: Fri, 24 Apr 2020 20:04:53 +1200 Subject: [PATCH 15/15] Test heart refinement, surface groups and markers --- tests/test_heart.py | 121 +++++++++++++++++++++++++++++++++++++++----- 1 file changed, 109 insertions(+), 12 deletions(-) diff --git a/tests/test_heart.py b/tests/test_heart.py index b5d08df5..ec4de5d0 100644 --- a/tests/test_heart.py +++ b/tests/test_heart.py @@ -1,10 +1,13 @@ import unittest -from opencmiss.utils.zinc.finiteelement import evaluateFieldNodesetRange +import copy +from opencmiss.utils.zinc.finiteelement import evaluateFieldNodesetRange, findNodeWithName from opencmiss.zinc.context import Context from opencmiss.zinc.field import Field from opencmiss.zinc.result import RESULT_OK -from scaffoldmaker.annotation.annotationgroup import AnnotationGroup +from scaffoldmaker.annotation.annotationgroup import AnnotationGroup, getAnnotationGroupForTerm +from scaffoldmaker.annotation.heart_terms import get_heart_term from scaffoldmaker.meshtypes.meshtype_3d_heart1 import MeshType_3d_heart1 +from scaffoldmaker.utils.meshrefinement import MeshRefinement from testutils import assertAlmostEqualList class HeartScaffoldTestCase(unittest.TestCase): @@ -13,10 +16,11 @@ def test_heart1(self): """ Test creation of heart scaffold. """ - parameterSetNames = MeshType_3d_heart1.getParameterSetNames() + scaffold = MeshType_3d_heart1 + parameterSetNames = scaffold.getParameterSetNames() self.assertEqual(parameterSetNames, [ "Default", "Human 1", "Mouse 1", "Pig 1", "Rat 1", "Unit Human 1", "Unit Mouse 1", "Unit Pig 1", "Unit Rat 1" ]); - options = MeshType_3d_heart1.getDefaultOptions("Human 1") + options = scaffold.getDefaultOptions("Human 1") self.assertEqual(116, len(options)) self.assertEqual(0.9, options.get("LV outer height")) self.assertEqual(80.0, options.get("Unit scale")) @@ -25,14 +29,9 @@ def test_heart1(self): context = Context("Test") region = context.getDefaultRegion() self.assertTrue(region.isValid()) - annotationGroups = MeshType_3d_heart1.generateBaseMesh(region, options) - - self.assertEqual(24, len(annotationGroups)) + annotationGroups = scaffold.generateMesh(region, options) + self.assertEqual(25, len(annotationGroups)) fieldmodule = region.getFieldmodule() - self.assertEqual(RESULT_OK, fieldmodule.defineAllFaces()) - if annotationGroups is not None: - for annotationGroup in annotationGroups: - annotationGroup.addSubelements() mesh3d = fieldmodule.findMeshByDimension(3) self.assertEqual(289, mesh3d.getSize()) mesh2d = fieldmodule.findMeshByDimension(2) @@ -40,7 +39,7 @@ def test_heart1(self): mesh1d = fieldmodule.findMeshByDimension(1) self.assertEqual(1356, mesh1d.getSize()) nodes = fieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES) - self.assertEqual(530, nodes.getSize()) + self.assertEqual(528, nodes.getSize()) datapoints = fieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_DATAPOINTS) self.assertEqual(0, datapoints.getSize()) @@ -50,5 +49,103 @@ def test_heart1(self): assertAlmostEqualList(self, minimums, [ -50.7876375290527, -57.76590573823474, -91.6 ], 1.0E-6) assertAlmostEqualList(self, maximums, [ 43.81084359764995, 39.03925080604259, 40.71693637558552 ], 1.0E-6) + # check some annotationGroups: + expectedSizes3d = { + "left ventricle myocardium" : 94, + "right ventricle myocardium" : 79, + "interventricular septum" : 30, + "left atrium myocardium" : 88, + "right atrium myocardium" : 76, + "interatrial septum" : 17 + } + for name in expectedSizes3d: + group = getAnnotationGroupForTerm(annotationGroups, get_heart_term(name)) + size = group.getMeshGroup(mesh3d).getSize() + self.assertEqual(expectedSizes3d[name], size, name) + expectedSizes2d = { + "endocardium of left ventricle" : 74, + "endocardium of right ventricle" : 59, + "endocardium of left atrium" : 82, + "endocardium of right atrium" : 70, + "epicardium" : 211 + } + for name in expectedSizes2d: + group = getAnnotationGroupForTerm(annotationGroups, get_heart_term(name)) + size = group.getMeshGroup(mesh2d).getSize() + self.assertEqual(expectedSizes2d[name], size, name) + + # refine 2x2x2 and check result + # first remove any surface annotation groups as they are re-added by defineFaceAnnotations + removeAnnotationGroups = [] + for annotationGroup in annotationGroups: + if annotationGroup.getMeshGroup(mesh3d).getSize() == 0: + removeAnnotationGroups.append(annotationGroup) + for annotationGroup in removeAnnotationGroups: + annotationGroups.remove(annotationGroup) + self.assertEqual(16, len(annotationGroups)) + + refineRegion = region.createRegion() + refineFieldmodule = refineRegion.getFieldmodule() + options['Refine number of elements surface'] = 2 + options['Refine number of elements through LV wall'] = 2 + options['Refine number of elements through wall'] = 2 + 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(25, len(annotationGroups)) + + mesh3d = refineFieldmodule.findMeshByDimension(3) + self.assertEqual(2236, mesh3d.getSize()) + mesh2d = refineFieldmodule.findMeshByDimension(2) + self.assertEqual(7676, mesh2d.getSize()) + mesh1d = refineFieldmodule.findMeshByDimension(1) + self.assertEqual(8623, mesh1d.getSize()) + nodes = refineFieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES) + self.assertEqual(3183, 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_heart_term(name)) + size = group.getMeshGroup(mesh3d).getSize() + self.assertEqual(expectedSizes3d[name]*8, size, name) + for name in expectedSizes2d: + group = getAnnotationGroupForTerm(annotationGroups, get_heart_term(name)) + size = group.getMeshGroup(mesh2d).getSize() + if name == "epicardium": + # fibrous ring is only refined by 1 through its thickness + fibrousRingReduction = (options['Refine number of elements surface'] - 1)*options['Refine number of elements surface']*( + options['Number of elements around left atrium free wall'] + options['Number of elements around right atrium free wall']) + self.assertEqual(expectedSizes2d[name]*4 - fibrousRingReduction, size, name) + else: + self.assertEqual(expectedSizes2d[name]*4, 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(5, 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, "APEX") + self.assertTrue(node.isValid()) + cache.setNode(node) + element, xi = markerLocation.evaluateMeshLocation(cache, 3) + self.assertEqual(5, element.getIdentifier()) + assertAlmostEqualList(self, xi, [ 0.0, 0.0, 1.0 ], 1.0E-10) + if __name__ == "__main__": unittest.main()