diff --git a/docs/_images/scaffoldmaker_human_heart.jpg b/docs/_images/scaffoldmaker_human_heart.jpg new file mode 100644 index 00000000..6619c139 Binary files /dev/null and b/docs/_images/scaffoldmaker_human_heart.jpg differ diff --git a/docs/conf.py b/docs/conf.py index db8b0960..1e5cb046 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -2,7 +2,7 @@ # -- Project information -project = 'Scaffold Maker' +project = 'scaffoldmaker' copyright = '2022, University of Auckland' author = 'University of Auckland' @@ -33,3 +33,6 @@ # -- Options for EPUB output epub_show_urls = 'footnote' + +# If true, enable figure numbering +numfig = True diff --git a/docs/index.rst b/docs/index.rst index 7d83f837..ee6659e1 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -1,14 +1,13 @@ -Scaffold Maker -============== +*scaffoldmaker library* +======================= -The **Scaffold Maker** is part of the software that is used in the collection of tools used for mapping data to scaffolds. +The *scaffoldmaker library* contains scripts for programmatically generating anatomical scaffolds for a range of organs and other structures, and includes utility code for building common geometric shapes, annotation, and refinement. The scripts generate the scaffold model in the data structures of the underlying *OpenCMISS-Zinc library*, and through the Zinc API clients can interrogate the model or export it for subsequent use. -.. note:: - - This project is under active development. +Most users will make scaffolds using the ABI Mapping Tools' **Scaffold Creator** user interface for this library, and its documentation gives a good introduction to anatomical scaffolds and some common features of them. +This documentation is intended to be a resource describing details about using individual scaffolds, as well as developing new ones. .. toctree:: install - + scaffolds/heart diff --git a/docs/install.rst b/docs/install.rst index 7b241847..dee6e930 100644 --- a/docs/install.rst +++ b/docs/install.rst @@ -1,6 +1,5 @@ - -Installing -========== +Installation +============ Install with:: diff --git a/docs/scaffolds/heart.rst b/docs/scaffolds/heart.rst new file mode 100644 index 00000000..8dca1ee1 --- /dev/null +++ b/docs/scaffolds/heart.rst @@ -0,0 +1,83 @@ +Heart Scaffold +============== + +The current recommended heart scaffold is ``3D Heart 1`` built from ``class MeshType_3d_heart1``; the human variant is shown in :numref:`fig-scaffoldmaker-human-heart`. + +.. _fig-scaffoldmaker-human-heart: + +.. figure:: ../_images/scaffoldmaker_human_heart.jpg + :align: center + + Human heart scaffold. + +The heart scaffold is a 3-D volumetric model of the heart representing solid walls of all four chambers. It has pulmonary vein and vena cava inlets, simple representations of auricles (atrial appendages), but does not have representations of the heart valves (only their orifices, or boundaries with the ventricles). Only the myocardium layer is fully represented (but see below). + +.. note:: + + Separate scaffolds are provided for the heart ventricles (below base/valve plane), ventricles with base, and atria. Beware that the heart atria scaffold is not compatible with the full heart scaffold because it uses different numbering for elements, faces and nodes; if working on a subset of the heart it is recommended to create the whole heart scaffold and delete the unwanted elements, which allows the resulting subset to be merged with the whole heart. Further to this, all scaffolds used must be built with the same parameters for numbers of elements, inlets, optional layers etc. to be able to merge together. + +Variants +-------- + +The heart scaffold is provided with parameter sets for the following four species, which differ in shape, and in particular have different numbers of pulmonary veins: + +* Human (2 left, 2 right pulmonary veins) +* Pig (1 left, 1 right pulmonary veins) +* Rat (3 pulmonary veins: left, middle, right) +* Mouse (3 pulmonary veins: left, middle, right) + +These variants' geometry and annotations are best viewed in the **Scaffold Creator** tool in the ABI Mapping Tools. On the web, the latest published generic heart scaffold variants can be viewed on the `SPARC Portal `_ by searching for ``heart``, filtering for models, selecting a variant and viewing the scaffold in its Gallery tab. + +The heart scaffold script generates the scaffold mesh and geometry from ellipsoid and cubic functions with many parameters controlling the shape. The parameters were carefully tuned for each species, and it is not recommended that these be edited. + +An advanced optional feature is to check *Define epicardium layer* (set parameter to ``true``) which adds a layer of 3-D elements outside the myocardium to represent the thick epicardium layer consisting of epicardial fat and other tissue. This is currently only implemented over the atria, excluding the auricles. + +Coordinates +----------- + +The heart scaffold only defines geometric coordinates (field ``coordinates``) which give geometry at approximately unit scale. Note that the scaffold has a *Unit scale* parameter (default value ``1.0``) which scales the entire scaffold efficiently. + +A material coordinates field is not provided, so to perform embedding at this time, it is necessary to use the generic ``coordinates`` field as material coordinates, built with the standard, *unmodified* parameter set (incl. *Unit scale* ``1.0``) for the species. + +The heart scaffold supports limited refinement/resampling by checking *Refine* (set parameter to ``true``) with chosen *Refine number of elements~* parameters. Be aware that only the ``coordinates`` field is currently defined on the refined mesh (but annotations are transferred), and the refined whole heart is only conformant if all *Refine number of elements~* parameters have the same value. + +Annotations +----------- + +Important anatomical regions of the heart are defined by groups of elements (or faces, edges and nodes/points) and annotated with standard term names and identifiers from a controlled vocabulary. + +Annotated 3-dimensional volume regions are defined by groups of 3-D elements including (using only one of the items separated by slash /): + +* left/right atrium/ventricle myocardium +* left/middle/right pulmonary vein +* inferior/superior vena cava +* left/right auricle + +**Terms for volume regions such as the above are not to be used for digitized contours!** They are used for applying different material properties in models and the strain/curvature penalty (stiffness) parameters in fitting. + +Annotated 2-dimensional surface regions are defined for matching annotated contours digitized from medical images including (where ``luminal`` means bordering the lumen or cavity of a tubular structure, ``outer`` is the outside boundary of a layer on the wall relative to the cavity): + +* luminal surface of left/right atrium/ventricle +* luminal surface of left/middle/right pulmonary vein +* luminal surface of inferior/superior vena cava +* outer surface of myocardium +* outer surface of myocardium of left/right atrium/ventricle (if need to distinguish) +* outer surface of epicardium (only if *Define epicardium layer* is checked) + +The following are proposed for future heart scaffolds including great vessels (the root variants being the subset of the luminal surface within the arterial valve up to the sinotubular junction) so should be used for annotating their contours: + +* luminal surface of [root of] aorta +* luminal surface of [root of] pulmonary trunk + +Note that luminal surfaces of atria structures are defined with an overlap in the scaffold as it is difficult to determine the exact boundary between pulmonary veins/vena cavae with the myocardium for the left/right atrium, which gives leeway for variation between digitized datasets. + +Note that at organ scale additional terms such as ``endocardium of left ventricle`` which represent microscopically thin layers are defined identically to ``luminal surface of left ventricle`` etc. surface groups. However, this is not the case on the epicardium which can be a layer of finite thickness at whole organ scale. + +The heart scaffold currently has no annotated 1-dimensional line regions. + +Several fiducial marker points are defined on the heart scaffold, of which the following two are potentially usable when digitizing: + +* apex of heart +* crux cordis + +At present these are both defined on the outer surface of myocardium, but when a volumetric epicardium layer is defined over the whole heart these will either be defined on the outer surface of epicardium, or separate points defined to distinguish distinct points on the two surfaces. diff --git a/src/scaffoldmaker/annotation/heart_terms.py b/src/scaffoldmaker/annotation/heart_terms.py index a88aa4d5..04961f42 100644 --- a/src/scaffoldmaker/annotation/heart_terms.py +++ b/src/scaffoldmaker/annotation/heart_terms.py @@ -4,60 +4,93 @@ # convention: preferred name, preferred id, followed by any other ids and alternative names heart_terms = [ - ( "heart", "UBERON:0000948", "FMA:7088" ), - # 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" ), - ( "epicardial fat", "UBERON:0015129"), - ( "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 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" ), - ( "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" ), - ( "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" ), + # heart - volume terms + ("heart", "UBERON:0000948", "FMA:7088"), # group of the entire heart + ("epicardial fat", "UBERON:0015129"), # not used + ("epicardium", "UBERON:0002348", "FMA:9461"), # volumetric layer outside myocardium to pericardial cavity + # heart - surface terms + ("outer surface of epicardium", "ILX:0793548"), + ("outer surface of myocardium", "ILX:0793530"), + + # ventricles - volume terms + ("conus arteriosus", "UBERON:0003983"), + ("endocardium of left ventricle", "UBERON:0009713", "FMA:9559"), + ("endocardium of right ventricle", "UBERON:0009712", "FMA:9536"), + ("heart left ventricle", "UBERON:0002084"), # the whole left ventricle + ("heart right ventricle", "UBERON:0002080"), # the whole right ventricle + ("interventricular septum", "UBERON:0002094", "FMA:7133"), + ("left fibrous ring", "FMA:77124"), + ("left ventricle myocardium", "UBERON:0006566", "FMA:9558"), + ("right fibrous ring", "FMA:77125"), + ("right ventricle myocardium", "UBERON:0006567", "FMA:9535"), + # ventricles - surface terms + ("luminal surface of left ventricle", "ILX:0793537"), + ("luminal surface of right ventricle", "ILX:0793538"), + ("outer surface of myocardium of left ventricle", "ILX:0793533"), + ("outer surface of myocardium of right ventricle", "ILX:0793534"), + + # atria - volume terms + ("left atrium endocardium", "UBERON:0034903", "FMA:7286"), + ("endocardium of left auricle", "UBERON:0011006", "FMA:13236"), + ("right atrium endocardium", "UBERON:0009129", "FMA:7281"), + ("endocardium of right auricle", "UBERON:0011007", "FMA:13235"), + ("epicardium of left auricle", "FMA:13233"), + ("epicardium of right auricle", "FMA:13232"), + ("fossa ovalis", "UBERON:0003369", "FMA:9246"), + ("inferior vena cava", "UBERON:0001072", "FMA:10951", "posterior vena cava"), + ("inferior vena cava inlet", "ILX:0738358"), + ("interatrial septum", "UBERON:0002085", "FMA:7108"), + ("left atrium myocardium", "FMA:7285"), + ("left auricle", "UBERON:0006630", "FMA:7219"), + ("left cardiac atrium", "UBERON:0002079"), # the whole left atrium + ("left inferior pulmonary vein", "FMA:49913"), + ("left pulmonary vein", "UBERON:0009030"), + ("left superior pulmonary vein", "FMA:49916"), + ("middle pulmonary vein", "ILX:0739222"), # in mouse, rat, rabbit + ("pulmonary vein", "UBERON:0002016", "FMA:66643"), + ("right atrium myocardium", "FMA:7282"), + ("right auricle", "UBERON:0006631", "FMA:7218"), + ("right cardiac atrium", "UBERON:0002078"), # the whole right atrium + ("right inferior pulmonary vein", "FMA:49911"), + ("right pulmonary vein", "UBERON:0009032"), + ("right superior pulmonary vein", "FMA:49914"), + ("superior vena cava", "UBERON:0001585", "FMA:4720", "anterior vena cava"), + ("superior vena cava inlet", "ILX:0738367"), + # atria - surface terms + ("luminal surface of left atrium", "ILX:0793535"), + ("luminal surface of right atrium", "ILX:0793536"), + ("luminal surface of inferior vena cava", "ILX:0793542"), + ("luminal surface of left pulmonary vein", "ILX:0793539"), + ("luminal surface of middle pulmonary vein", "ILX:0793540"), + ("luminal surface of right pulmonary vein", "ILX:0793541"), + ("luminal surface of superior vena cava", "ILX:0793543"), + ("outer surface of myocardium of left atrium", "ILX:0793531"), + ("outer surface of myocardium of right atrium", "ILX:0793532"), + + # arterial valves and great vessels - volume terms + ("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"), + # arterial valves and great vessels - surface terms + ("luminal surface of aorta", "ILX:0793544"), + ("luminal surface of pulmonary trunk", "ILX:0793546"), + ("luminal surface of root of aorta", "ILX:0793545"), + ("luminal surface of root of pulmonary trunk", "ILX:0793547"), + # fiducial markers - ( "apex of heart", "UBERON:0002098", "FMA:7164"), - ( "crux of heart", "ILX:0777104", "FMA:7220" ), # point on posterior surface where the four chambers meet on interventricular, atrio-ventricular and interatrial sulci - ( "left atrium epicardium venous midpoint", "ILX:0778116"), # point at centre of pulmonary vein ostia on left atrium epicardium, on anterior/ventral side for rodents - ( "right atrium epicardium venous midpoint", "ILX:0778117") # point at centre of inferior & superior vena cavae on right atrium epicardium - ] + ("apex of heart", "UBERON:0002098", "FMA:7164"), + # point on posterior surface where the four chambers meet on interventricular, A-V and interatrial sulci + ("crux cordis", "ILX:0777104", "FMA:7220"), + # point at centre of pulmonary vein ostia on left atrium epicardium, on anterior/ventral side for rodents + ("left atrium epicardium venous midpoint", "ILX:0778116"), + # point at centre of inferior & superior vena cavae on right atrium epicardium + ("right atrium epicardium venous midpoint", "ILX:0778117") +] def get_heart_term(name : str): """ diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_heart1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_heart1.py index f5cb1568..5beb4b5c 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_heart1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_heart1.py @@ -36,11 +36,7 @@ def getParameterSetNames(): 'Human 1', 'Mouse 1', 'Pig 1', - 'Rat 1', - 'Unit Human 1', - 'Unit Mouse 1', - 'Unit Pig 1', - 'Unit Rat 1'] + 'Rat 1'] @staticmethod def getDefaultOptions(parameterSetName='Default'): @@ -83,7 +79,7 @@ def getOrderedOptionNames(): 'Refine number of elements surface', 'Refine number of elements through LV wall', 'Refine number of elements through wall', - 'Refine number of elements through epicardial fat layer']: + 'Refine number of elements through epicardium layer']: optionNames.remove(optionName) optionNames.append(optionName) return optionNames @@ -133,32 +129,20 @@ def generateBaseMesh(cls, region, options): coordinates = findOrCreateFieldCoordinates(fm) cache = fm.createFieldcache() - mesh = fm.findMeshByDimension(3) - # generate heartventriclesbase1 model and put atria1 on it ventriclesAnnotationGroups = MeshType_3d_heartventriclesbase1.generateBaseMesh(region, options) atriaAnnotationGroups = MeshType_3d_heartatria1.generateBaseMesh(region, options) annotationGroups = mergeAnnotationGroups(ventriclesAnnotationGroups, atriaAnnotationGroups) heartGroup = findOrCreateAnnotationGroupForTerm(annotationGroups, region, get_heart_term("heart")) - cruxGroup = findOrCreateAnnotationGroupForTerm(annotationGroups, region, get_heart_term("crux of heart")) + cruxGroup = findOrCreateAnnotationGroupForTerm(annotationGroups, region, get_heart_term("crux cordis")) 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") - 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(markerName) - markerTemplateInternal.defineField(markerLocation) - ############## # Create nodes ############## + nodes = fm.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES) nodeIdentifier = max(1, getMaximumNodeIdentifier(nodes) + 1) # discover left and right fibrous ring nodes from ventricles and atria @@ -219,6 +203,7 @@ def generateBaseMesh(cls, region, options): # Create elements ################# + mesh = fm.findMeshByDimension(3) heartMeshGroup = heartGroup.getMeshGroup(mesh) lFibrousRingMeshGroup = lFibrousRingGroup.getMeshGroup(mesh) rFibrousRingMeshGroup = rFibrousRingGroup.getMeshGroup(mesh) @@ -397,16 +382,13 @@ def generateBaseMesh(cls, region, options): for meshGroup in meshGroups: meshGroup.addElement(element) - # annotation fiducial points + # crux cordis annotation point cruxElement = mesh.findElementByIdentifier(cruxElementId) cruxXi = [ 0.5, 0.5, 1.0 ] - markerPoint = markerPoints.createNode(nodeIdentifier, markerTemplateInternal) - nodeIdentifier += 1 - cache.setNode(markerPoint) - markerName.assignString(cache, cruxGroup.getName()) - markerLocation.assignMeshLocation(cache, cruxElement, cruxXi) - for group in [ heartGroup, cruxGroup ]: - group.getNodesetGroup(nodes).addNode(markerPoint) + markerNode = cruxGroup.createMarkerNode(nodeIdentifier, element=cruxElement, xi=cruxXi) + nodeIdentifier = markerNode.getIdentifier() + 1 + for group in [ heartGroup ]: + group.getNodesetGroup(nodes).addNode(markerNode) return annotationGroups diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_heartatria1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_heartatria1.py index a587e481..0610cb84 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_heartatria1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_heartatria1.py @@ -16,8 +16,8 @@ from opencmiss.zinc.field import Field, FieldGroup from opencmiss.zinc.node import Node from opencmiss.zinc.result import RESULT_OK -from scaffoldmaker.annotation.annotationgroup import AnnotationGroup, findOrCreateAnnotationGroupForTerm, \ - getAnnotationGroupForTerm +from scaffoldmaker.annotation.annotationgroup import AnnotationGroup, findAnnotationGroupByName, \ + findOrCreateAnnotationGroupForTerm, getAnnotationGroupForTerm from scaffoldmaker.annotation.heart_terms import get_heart_term from scaffoldmaker.meshtypes.meshtype_3d_ostium1 import MeshType_3d_ostium1, generateOstiumMesh from scaffoldmaker.meshtypes.scaffold_base import Scaffold_base @@ -203,11 +203,7 @@ def getParameterSetNames(): 'Human 1', 'Mouse 1', 'Pig 1', - 'Rat 1', - 'Unit Human 1', - 'Unit Mouse 1', - 'Unit Pig 1', - 'Unit Rat 1'] + 'Rat 1'] @classmethod def getDefaultOptions(cls, parameterSetName='Default'): @@ -215,7 +211,6 @@ def getDefaultOptions(cls, parameterSetName='Default'): isMouse = 'Mouse' in parameterSetName isPig = 'Pig' in parameterSetName isRat = 'Rat' in parameterSetName - notUnitScale = 'Unit' not in parameterSetName if isPig: lpvOstium = cls.lpvOstiumDefaultScaffoldPackages['LPV Pig 1'] rpvOstium = cls.rpvOstiumDefaultScaffoldPackages['RPV Pig 1'] @@ -314,20 +309,15 @@ def getDefaultOptions(cls, parameterSetName='Default'): options['Superior vena cava inlet length'] = 0.1 options['Superior vena cava inlet inner diameter'] = 0.18 options['Superior vena cava inlet wall thickness'] = 0.015 - options['Define epicardial fat layer'] = False - options['Epicardial fat minimum thickness'] = 0.01 + options['Define epicardium layer'] = False + options['Epicardium layer minimum thickness'] = 0.01 options['Refine'] = False options['Refine number of elements surface'] = 4 options['Refine number of elements through wall'] = 1 - options['Refine number of elements through epicardial fat layer'] = 1 + options['Refine number of elements through epicardium layer'] = 1 options['Use cross derivatives'] = False - if isHuman: - if notUnitScale: - options['Unit scale'] = 80.0 - elif isMouse or isRat: - if notUnitScale: - options['Unit scale'] = 6.0 if isMouse else 12.0 + if isMouse or isRat: options['Atria base inner major axis length'] = 0.32 options['Atria base inner minor axis length'] = 0.26 options['Atria major axis rotation degrees'] = 30.0 @@ -394,9 +384,7 @@ def getDefaultOptions(cls, parameterSetName='Default'): options['Superior vena cava inlet position over'] = 0.62 options['Superior vena cava inlet position right'] = 0.25 options['Superior vena cava inlet wall thickness'] = 0.012 - elif 'Pig' in parameterSetName: - if notUnitScale: - options['Unit scale'] = 80.0 + elif isPig: options['Atrial base side incline degrees'] = 0.0 options['Left atrial appendage angle axial degrees'] = -10.0 options['Left atrial appendage angle left degrees'] = 20.0 @@ -517,12 +505,12 @@ def getOrderedOptionNames(): 'Superior vena cava inlet length', 'Superior vena cava inlet inner diameter', 'Superior vena cava inlet wall thickness', - 'Define epicardial fat layer', - 'Epicardial fat minimum thickness', + 'Define epicardium layer', + 'Epicardium layer minimum thickness', 'Refine', 'Refine number of elements surface', 'Refine number of elements through wall', - 'Refine number of elements through epicardial fat layer' + 'Refine number of elements through epicardium layer' #,'Use cross derivatives' ] @@ -622,7 +610,7 @@ def checkOptions(cls, options): 'Superior vena cava inlet length', 'Superior vena cava inlet inner diameter', 'Superior vena cava inlet wall thickness', - 'Epicardial fat minimum thickness']: + 'Epicardium layer minimum thickness']: if options[key] < 0.0: options[key] = 0.0 for key in [ @@ -805,13 +793,11 @@ def generateBaseMesh(cls, region, options): svcLength = unitScale*options['Superior vena cava inlet length'] svcInnerRadius = unitScale*0.5*options['Superior vena cava inlet inner diameter'] svcWallThickness = unitScale*options['Superior vena cava inlet wall thickness'] - defineEpicardialFatLayer = options['Define epicardial fat layer'] - epicardialFatMinimumThickness = unitScale*options['Epicardial fat minimum thickness'] + defineEpicardiumLayer = options['Define epicardium layer'] + epicardiumLayerMinimumThickness = unitScale*options['Epicardium layer minimum thickness'] 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() @@ -823,9 +809,9 @@ def generateBaseMesh(cls, region, options): laaGroup = AnnotationGroup(region, get_heart_term("left auricle")) raaGroup = AnnotationGroup(region, get_heart_term("right auricle")) annotationGroups = [heartGroup, lamGroup, ramGroup, aSeptumGroup, fossaGroup, laaGroup, raaGroup] - if defineEpicardialFatLayer: - epicardialFatGroup = AnnotationGroup(region, get_heart_term("epicardial fat")) - annotationGroups.append(epicardialFatGroup) + if defineEpicardiumLayer: + epicardiumGroup = AnnotationGroup(region, get_heart_term("epicardium")) + annotationGroups.append(epicardiumGroup) lpvOstiumSettings = lpvOstium.getScaffoldSettings() lpvCount = lpvOstiumSettings['Number of vessels'] @@ -865,20 +851,11 @@ def generateBaseMesh(cls, region, options): annotationGroups += [laeVenousMidpointGroup, ivcGroup, ivcInletGroup, raeVenousMidpointGroup, svcGroup, svcInletGroup, lFibrousRingGroup, rFibrousRingGroup] - # 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) @@ -894,6 +871,7 @@ def generateBaseMesh(cls, region, options): nodeIdentifier = max(1, getMaximumNodeIdentifier(nodes) + 1) + mesh = fm.findMeshByDimension(3) elementIdentifier = max(1, getMaximumElementIdentifier(mesh) + 1) heartMeshGroup = heartGroup.getMeshGroup(mesh) @@ -2052,27 +2030,25 @@ def generateBaseMesh(cls, region, options): coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, ltBaseOuterd2[n1]) nodeIdentifier += 1 - if defineEpicardialFatLayer: + if defineEpicardiumLayer: # epicardial fat pad in RAGP region -- track surface bridging interatrial groove fpx = [] fpd1 = [] fpd2 = [] - fpd3 = [] - proportionAcross = 0.1 + proportionAcross = 0.2 for i in range(2): if i == 0: - nx, nd2, nd1, nd3 = laTrackSurface.createHermiteCurvePoints( - 0.0, proportionAcross, 1.0, proportionAcross, elementsCountAcrossTrackSurface)[0:4] + nx, nd2, nd1, nd3, nProportions = laTrackSurface.createHermiteCurvePoints( + 0.0, proportionAcross, 1.0, proportionAcross, elementsCountAcrossTrackSurface) else: nx, nd2, nd1, nd3 = raTrackSurface.createHermiteCurvePoints( 1.0, proportionAcross, 0.0, proportionAcross, elementsCountAcrossTrackSurface)[0:4] if nx: - nx = [add(nx[i], mult(nd3[i], epicardialFatMinimumThickness)) for i in range(len(nx))] + nx = [add(nx[i], mult(nd3[i], epicardiumLayerMinimumThickness)) for i in range(len(nx))] nd1 = [[-c for c in d] for d in nd1] fpx.append(nx) fpd1.append(nd1) fpd2.append(nd2) - fpd3.append(nd3) # put into single arrays cycling left to right fastest, smoothing each d1 row nx = [] @@ -2082,7 +2058,8 @@ def generateBaseMesh(cls, region, options): scale = interp.computeCubicHermiteDerivativeScaling(fpx[0][n], fpd1[0][n], fpx[1][n], fpd1[1][n]) for i in range(2): nx.append(fpx[i][n]) - nd1.append(mult(fpd1[i][n], scale)) + fpd1[i][n] = mult(fpd1[i][n], scale) + nd1.append(fpd1[i][n]) nd2.append(fpd2[i][n]) fpTrackSurface = TrackSurface(1, elementsCountAcrossTrackSurface, nx, nd1, nd2) @@ -2644,13 +2621,10 @@ def generateBaseMesh(cls, region, options): laevmElementId = elementIdentifier - elementsCountAroundLpvOstium laevmXi = [ 0.0, 1.0, 1.0 ] laevmElement = mesh.findElementByIdentifier(laevmElementId) - markerPoint = markerPoints.createNode(nodeIdentifier, markerTemplateInternal) - nodeIdentifier += 1 - cache.setNode(markerPoint) - markerName.assignString(cache, laeVenousMidpointGroup.getName()) - markerLocation.assignMeshLocation(cache, laevmElement, laevmXi) - for group in [ heartGroup, lamGroup, laeVenousMidpointGroup ]: - group.getNodesetGroup(nodes).addNode(markerPoint) + markerNode = laeVenousMidpointGroup.createMarkerNode(nodeIdentifier, element=laevmElement, xi=laevmXi) + nodeIdentifier = markerNode.getIdentifier() + 1 + for group in [ heartGroup, lamGroup ]: + group.getNodesetGroup(nodes).addNode(markerNode) if not commonLeftRightPvOstium: # create right pulmonary vein annulus @@ -2952,13 +2926,10 @@ def generateBaseMesh(cls, region, options): raevmElementId = elementIdentifier - elementsCountAroundVC + elementsCountOverRightAtriumVenous//2 + elementsCountOverSideRightAtriumVC//2 raevmXi = [ 0.5 if (elementsCountOverSideRightAtriumVC % 2) else 0.0, 1.0, 1.0 ] raevmElement = mesh.findElementByIdentifier(raevmElementId) - markerPoint = markerPoints.createNode(nodeIdentifier, markerTemplateInternal) - nodeIdentifier += 1 - cache.setNode(markerPoint) - markerName.assignString(cache, raeVenousMidpointGroup.getName()) - markerLocation.assignMeshLocation(cache, raevmElement, raevmXi) - for group in [ heartGroup, ramGroup, raeVenousMidpointGroup ]: - group.getNodesetGroup(nodes).addNode(markerPoint) + markerNode = raeVenousMidpointGroup.createMarkerNode(nodeIdentifier, element=raevmElement, xi=raevmXi) + nodeIdentifier = markerNode.getIdentifier() + 1 + for group in [heartGroup, ramGroup]: + group.getNodesetGroup(nodes).addNode(markerNode) # create left atrial appendage position = laTrackSurface.createPositionProportion(laaMidpointOver, laaMidpointLeft) @@ -3085,7 +3056,6 @@ def generateBaseMesh(cls, region, options): eft1 = eft elementtemplate1 = elementtemplate scalefactors = None - addMarker = None nc = 2 + e1 if e2 == 0: @@ -3108,7 +3078,6 @@ def generateBaseMesh(cls, region, options): #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_DS1, [ ( Node.VALUE_LABEL_D_DS1, [1] ) ]) 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: @@ -3139,13 +3108,6 @@ def generateBaseMesh(cls, region, options): #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) @@ -3269,7 +3231,7 @@ def generateBaseMesh(cls, region, options): elementsCountRadial = elementsCountAlongAtrialAppendages, meshGroups = [ heartMeshGroup, ramMeshGroup, raaMeshGroup ]) - if defineEpicardialFatLayer: + if defineEpicardiumLayer: # project epicardial points over atria to build fat pad epiGroup = fm.createFieldGroup() epiMesh = epiGroup.createFieldElementGroup(mesh).getMeshGroup() @@ -3315,6 +3277,20 @@ def generateBaseMesh(cls, region, options): bridgeNodes = bridgeGroup.createFieldNodeGroup(nodes).getNodesetGroup() bridgeNodeTangents = {} + # parameters for shifting centre of fatpad to spread elements around PV, VC inlets + # have ivcPositionOver, svcPositionOver + rpvPositionOver = lpvOstiumPositionOver if commonLeftRightPvOstium else rpvOstiumPositionOver + # GRC fudge factors + svcSpread = ivcSpread = ivcPositionOver + rpvSpread = ivcSpread + ivcMag = 0.15 + svcMag = 0.0 + rpvMag = 0.15 + + for n in range(elementsCountAcrossTrackSurface + 1): + xi_ivc = max(1.0, math.fabs(ivcPositionOver - nProportions[n][1])) + xi_svc = max(1.0, math.fabs(ivcPositionOver - nProportions[n][1])) + for epiNodeIdentifier in epiFatPadNodeIdentifiersMap.keys(): epiNode = nodes.findNodeByIdentifier(epiNodeIdentifier) cache.setNode(epiNode) @@ -3324,16 +3300,43 @@ def generateBaseMesh(cls, region, options): result, epid3 = coordinates.getNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS3, 1, 3) if result != RESULT_OK: epid3 = vector.crossproduct3(epid1, epid2) - fatx = add(epix, vector.setMagnitude(epid3, epicardialFatMinimumThickness)) + fatx = add(epix, vector.setMagnitude(epid3, epicardiumLayerMinimumThickness)) epifx = None epifPosition = fpTrackSurface.findNearestPosition(epix, startPosition=None) - if 0.0001 < epifPosition.xi1 < 0.9999: - epifx, epifd1, epifd2 = fpTrackSurface.evaluateCoordinates(epifPosition, derivatives=True) + if not (((epifPosition.e1 == 0) and (epifPosition.xi1 < 0.0001)) or + ((epifPosition.e1 == (fpTrackSurface.elementsCount1 - 1)) and (0.9999 < epifPosition.xi1))): + # shift proportion 1 around inlet + xi_centre = 1.0 - (2.0 * math.fabs(0.5 - epifPosition.xi1)) + proportion1, proportion2 = fpTrackSurface.getProportion(epifPosition) + # scale to spread out centre + # proportion1Scaled = interp.interpolateCubicHermite([0.0], [0.5], [0.5], [0.6], xi_centre)[0] + # if proportion1 < 0.5: + # proportion1 = proportion1Scaled + # else: + # proportion1 = 1.0 - proportion1Scaled + proportion1Shift = 0.0 + xi_ivc = math.fabs((ivcPositionOver - proportion2) / ivcSpread) + if xi_ivc < 1.0: + proportion1Shift -= interp.interpolateCubicHermite([ivcMag], [0.0], [0.0], [0.0], xi_ivc)[0] + xi_svc = math.fabs((svcPositionOver - proportion2) / svcSpread) + if xi_svc < 1.0: + proportion1Shift -= interp.interpolateCubicHermite([svcMag], [0.0], [0.0], [0.0], xi_svc)[0] + xi_rpv = math.fabs((rpvPositionOver - proportion2) / rpvSpread) + if xi_rpv < 1.0: + proportion1Shift += interp.interpolateCubicHermite([rpvMag], [0.0], [0.0], [0.0], xi_rpv)[0] + proportion1Shift *= xi_centre + # proportion1Shift = interp.interpolateCubicHermite( + # [0.0], [0.0], [proportion1Shift], [0.0], xi_centre)[0] + proportion1 += proportion1Shift + + # convert shifted proportions to shifted position on fpTrackSurface + epifPosition2 = fpTrackSurface.createPositionProportion(proportion1, proportion2) + epifx, epifd1, epifd2 = fpTrackSurface.evaluateCoordinates(epifPosition2, derivatives=True) delta_epi = sub(epifx, epix) # epifx must be above the epicardium surface - # and at least epicardialFatMinimumThickness away from epix - if (dot(delta_epi, epid3) > 0.0) and (magnitude(delta_epi) >= epicardialFatMinimumThickness): + # and at least epicardiumLayerMinimumThickness away from epix + if (dot(delta_epi, epid3) > 0.0) and (magnitude(delta_epi) >= epicardiumLayerMinimumThickness): epifNormal = normalize(cross(epifd1, epifd2)) # epix must be under the fatpad if dot(delta_epi, epifNormal) > 0.0: @@ -3361,8 +3364,8 @@ def generateBaseMesh(cls, region, options): elementtemplateX = mesh.createElementtemplate() elementtemplateX.setElementShapeType(Element.SHAPE_TYPE_CUBE) - epicardialFatMeshGroup = epicardialFatGroup.getMeshGroup(mesh) - meshGroups = [heartMeshGroup, epicardialFatMeshGroup] + epicardiumMeshGroup = epicardiumGroup.getMeshGroup(mesh) + meshGroups = [heartMeshGroup, epicardiumMeshGroup] elementtemplate = mesh.createElementtemplate() # iterate over list of identifiers since can't iterate over mesh while modifying it for epiElementIdentifier in epiElementIdentifiers: @@ -3476,9 +3479,9 @@ def refineMesh(cls, meshrefinement, options): elementsCountAroundRightAtriumFreeWall = options['Number of elements around right atrium free wall'] refineElementsCountSurface = options['Refine number of elements surface'] refineElementsCountThroughWall = options['Refine number of elements through wall'] - refineElementsCountThroughEpicardialFatLayer =\ - options['Refine number of elements through epicardial fat layer'] - defineEpicardialFatLayer = options['Define epicardial fat layer'] + refineElementsCountThroughEpicardiumLayer =\ + options['Refine number of elements through epicardium layer'] + defineEpicardiumLayer = options['Define epicardium layer'] sourceFm = meshrefinement._sourceFm annotationGroups = meshrefinement._sourceAnnotationGroups @@ -3488,15 +3491,15 @@ def refineMesh(cls, meshrefinement, options): raElementGroupField = raGroup.getFieldElementGroup(meshrefinement._sourceMesh) aSeptumGroup = getAnnotationGroupForTerm(annotationGroups, get_heart_term("interatrial septum")) aSeptumMeshGroup = aSeptumGroup.getMeshGroup(meshrefinement._sourceMesh) - epicardialFatGroup = None - epicardialFatMeshGroup = None - if defineEpicardialFatLayer: - epicardialFatGroup = getAnnotationGroupForTerm(annotationGroups, get_heart_term("epicardial fat")) - epicardialFatMeshGroup = epicardialFatGroup.getMeshGroup(meshrefinement._sourceMesh) + epicardiumGroup = None + epicardiumMeshGroup = None + if defineEpicardiumLayer: + epicardiumGroup = getAnnotationGroupForTerm(annotationGroups, get_heart_term("epicardium")) + epicardiumMeshGroup = epicardiumGroup.getMeshGroup(meshrefinement._sourceMesh) coordinates = findOrCreateFieldCoordinates(meshrefinement._sourceFm) # last atria element is last element in following group: - lastGroup = epicardialFatGroup + lastGroup = epicardiumGroup if not lastGroup: lastGroup = getAnnotationGroupForTerm(annotationGroups, get_heart_term("right auricle")) lastMeshGroup = lastGroup.getMeshGroup(meshrefinement._sourceMesh) @@ -3525,8 +3528,8 @@ def refineMesh(cls, meshrefinement, options): refineElements1 = refineElementsCountThroughWall else: refineElements2 = refineElementsCountThroughWall - elif epicardialFatGroup and epicardialFatMeshGroup.containsElement(element): - refineElements3 = refineElementsCountThroughEpicardialFatLayer + elif epicardiumGroup and epicardiumMeshGroup.containsElement(element): + refineElements3 = refineElementsCountThroughEpicardiumLayer meshrefinement.refineElementCubeStandard3d(element, refineElements1, refineElements2, refineElements3) if elementIdentifier == lastElementIdentifier: return # finish on last so can continue elsewhere @@ -3544,15 +3547,21 @@ def defineFaceAnnotations(cls, region, options, annotationGroups): New face annotation groups are appended to this list. """ # create endocardium and epicardium groups - defineEpicardialFatLayer = options['Define epicardial fat layer'] + defineEpicardiumLayer = options['Define epicardium layer'] fm = region.getFieldmodule() lamGroup = getAnnotationGroupForTerm(annotationGroups, get_heart_term("left atrium myocardium")) ramGroup = 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")) - epicardialFatGroup = getAnnotationGroupForTerm(annotationGroups, get_heart_term("epicardial fat")) \ - if defineEpicardialFatLayer else None + lpvGroup = getAnnotationGroupForTerm(annotationGroups, get_heart_term("left pulmonary vein")) + # middle pulmonary vein is only present in rodents: + mpvGroup = findAnnotationGroupByName(annotationGroups, "middle pulmonary vein") + rpvGroup = getAnnotationGroupForTerm(annotationGroups, get_heart_term("right pulmonary vein")) + ivcGroup = getAnnotationGroupForTerm(annotationGroups, get_heart_term("inferior vena cava")) + svcGroup = getAnnotationGroupForTerm(annotationGroups, get_heart_term("superior vena cava")) + # following will already be defined if defineEpicardiumLayer is true + epiGroup = findOrCreateAnnotationGroupForTerm(annotationGroups, region, get_heart_term("epicardium")) mesh2d = fm.findMeshByDimension(2) is_exterior = fm.createFieldIsExterior() @@ -3563,68 +3572,114 @@ def defineFaceAnnotations(cls, region, options, annotationGroups): is_lam_endo = fm.createFieldAnd(is_lam, is_exterior_face_xi3_0) is_ram_endo = fm.createFieldOr(fm.createFieldAnd(fm.createFieldAnd(is_ram, is_exterior_face_xi3_0), fm.createFieldNot(is_lam_endo)), - fm.createFieldAnd(aSeptumGroup.getFieldElementGroup(mesh2d), is_exterior_face_xi3_1)) + fm.createFieldAnd(aSeptumGroup.getFieldElementGroup(mesh2d), + is_exterior_face_xi3_1)) is_laa = laaGroup.getFieldElementGroup(mesh2d) is_raa = raaGroup.getFieldElementGroup(mesh2d) is_laa_endo = fm.createFieldAnd(is_laa, is_exterior_face_xi3_0) is_raa_endo = fm.createFieldAnd(is_raa, 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) - is_myocardium = fm.createFieldOr(is_lam, is_ram) - is_a_epi = fm.createFieldAnd(fm.createFieldOr(is_myocardium, fm.createFieldOr(is_laa, is_raa)), - 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) - if defineEpicardialFatLayer: - # add non-exterior surfaces between myocardium and epicardial fat to epicardium - is_epicardial_fat = epicardialFatGroup.getFieldElementGroup(mesh2d) - epiGroup.getMeshGroup(mesh2d).addElementsConditional(fm.createFieldAnd(is_myocardium, is_epicardial_fat)) - laEndoGroup = findOrCreateAnnotationGroupForTerm(annotationGroups, region, get_heart_term("endocardium of left atrium")) - laEndoGroup.getMeshGroup(mesh2d).addElementsConditional(is_lam_endo) - raEndoGroup = findOrCreateAnnotationGroupForTerm(annotationGroups, region, get_heart_term("endocardium of right atrium")) - raEndoGroup.getMeshGroup(mesh2d).addElementsConditional(is_ram_endo) - laaEndoGroup = findOrCreateAnnotationGroupForTerm(annotationGroups, region, get_heart_term("endocardium of left auricle")) + # is_myocardium = fm.createFieldOr(is_lam, is_ram) + is_ext_xi3_1_and_not_septum = fm.createFieldAnd( + is_exterior_face_xi3_1, fm.createFieldNot(aSeptumGroup.getFieldElementGroup(mesh2d))) + is_os_lam = fm.createFieldAnd(is_lam, is_ext_xi3_1_and_not_septum) + is_os_ram = fm.createFieldAnd(is_ram, is_ext_xi3_1_and_not_septum) + is_epi = epiGroup.getFieldElementGroup(mesh2d) + + # luminal surfaces of endocardium of left/right atrium + lslaEndoGroup = findOrCreateAnnotationGroupForTerm( + annotationGroups, region, get_heart_term("luminal surface of left atrium")) + lslaEndoGroup.getMeshGroup(mesh2d).addElementsConditional(is_lam_endo) + lsraEndoGroup = findOrCreateAnnotationGroupForTerm( + annotationGroups, region, get_heart_term("luminal surface of right atrium")) + lsraEndoGroup.getMeshGroup(mesh2d).addElementsConditional(is_ram_endo) + # endocardium groups are defined identically to luminal surfaces at scaffold scale + laEndoGroup = findOrCreateAnnotationGroupForTerm( + annotationGroups, region, get_heart_term("left atrium endocardium")) + laEndoGroup.getMeshGroup(mesh2d).addElementsConditional(lslaEndoGroup.getFieldElementGroup(mesh2d)) + raEndoGroup = findOrCreateAnnotationGroupForTerm( + annotationGroups, region, get_heart_term("right atrium endocardium")) + raEndoGroup.getMeshGroup(mesh2d).addElementsConditional(lsraEndoGroup.getFieldElementGroup(mesh2d)) + + 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 = 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 = 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 = findOrCreateAnnotationGroupForTerm( + annotationGroups, region, get_heart_term("epicardium of right auricle")) raaEpiGroup.getMeshGroup(mesh2d).addElementsConditional(is_raa_epi) + oslamGroup = findOrCreateAnnotationGroupForTerm( + annotationGroups, region, get_heart_term("outer surface of myocardium of left atrium")) + oslamGroup.getMeshGroup(mesh2d).addElementsConditional(is_os_lam) + osramGroup = findOrCreateAnnotationGroupForTerm( + annotationGroups, region, get_heart_term("outer surface of myocardium of right atrium")) + osramGroup.getMeshGroup(mesh2d).addElementsConditional(is_os_ram) + if defineEpicardiumLayer: + oslamGroup.getMeshGroup(mesh2d).addElementsConditional(fm.createFieldAnd(is_lam, is_epi)) + osramGroup.getMeshGroup(mesh2d).addElementsConditional(fm.createFieldAnd(is_ram, is_epi)) + osmGroup = findOrCreateAnnotationGroupForTerm( + annotationGroups, region, get_heart_term("outer surface of myocardium")) + osmGroup.getMeshGroup(mesh2d).addElementsConditional(oslamGroup.getFieldElementGroup(mesh2d)) + osmGroup.getMeshGroup(mesh2d).addElementsConditional(osramGroup.getFieldElementGroup(mesh2d)) + if defineEpicardiumLayer: + # future: limit to atria once ventricles have epicardium layer + is_os_epi = fm.createFieldAnd(is_epi, is_exterior_face_xi3_1) + osEpiGroup = findOrCreateAnnotationGroupForTerm( + annotationGroups, region,get_heart_term("outer surface of epicardium")) + osEpiGroup.getMeshGroup(mesh2d).addElementsConditional(is_os_epi) + # note: epiGroup only contains 3-D elements in this case + else: + # if no volumetric epicardium group, add outer surface of atrial myocardium + epiGroup.getMeshGroup(mesh2d).addElementsConditional(oslamGroup.getFieldElementGroup(mesh2d)) + epiGroup.getMeshGroup(mesh2d).addElementsConditional(osramGroup.getFieldElementGroup(mesh2d)) + + lslpvGroup = findOrCreateAnnotationGroupForTerm( + annotationGroups, region, get_heart_term("luminal surface of left pulmonary vein")) + lsrpvGroup = findOrCreateAnnotationGroupForTerm( + annotationGroups, region, get_heart_term("luminal surface of right pulmonary vein")) + is_lpv = lpvGroup.getFieldElementGroup(mesh2d) + is_mpv = mpvGroup.getFieldElementGroup(mesh2d) if mpvGroup else None + is_rpv = rpvGroup.getFieldElementGroup(mesh2d) + lslpvGroup.getMeshGroup(mesh2d).addElementsConditional( + fm.createFieldAnd(is_exterior_face_xi3_0, is_lpv)) + if mpvGroup: + lsmpvGroup = findOrCreateAnnotationGroupForTerm( + annotationGroups, region, get_heart_term("luminal surface of middle pulmonary vein")) + lsmpvGroup.getMeshGroup(mesh2d).addElementsConditional( + fm.createFieldAnd(is_exterior_face_xi3_0, is_mpv)) + lsrpvGroup.getMeshGroup(mesh2d).addElementsConditional( + fm.createFieldAnd(is_exterior_face_xi3_0, is_rpv)) + + lsivcGroup = findOrCreateAnnotationGroupForTerm( + annotationGroups, region, get_heart_term("luminal surface of inferior vena cava")) + lssvcGroup = findOrCreateAnnotationGroupForTerm( + annotationGroups, region, get_heart_term("luminal surface of superior vena cava")) + is_ivc = ivcGroup.getFieldElementGroup(mesh2d) + is_svc = svcGroup.getFieldElementGroup(mesh2d) + lsivcGroup.getMeshGroup(mesh2d).addElementsConditional( + fm.createFieldAnd(is_exterior_face_xi3_0, is_ivc)) + lssvcGroup.getMeshGroup(mesh2d).addElementsConditional( + fm.createFieldAnd(is_exterior_face_xi3_0, is_svc)) + lFibrousRingGroup = getAnnotationGroupForTerm(annotationGroups, get_heart_term("left fibrous ring")) rFibrousRingGroup = getAnnotationGroupForTerm(annotationGroups, get_heart_term("right fibrous ring")) - if (lFibrousRingGroup.getDimension() == 0) or (lFibrousRingGroup.getDimension() == 0): - # not already added by full heart scaffold - nodes = fm.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES) - lFibrousRingNodeGroup = lFibrousRingGroup.getNodesetGroup(nodes) - rFibrousRingNodeGroup = rFibrousRingGroup.getNodesetGroup(nodes) - # make temp group containing all elements, faces etc. if have any nodes in fibrous ring groups - tmpGroup = fm.createFieldGroup() - tmpGroup.setSubelementHandlingMode(FieldGroup.SUBELEMENT_HANDLING_MODE_FULL) - mesh3d = fm.findMeshByDimension(3) - tmp3dMeshGroup = tmpGroup.createFieldElementGroup(mesh3d).getMeshGroup() - coordinates = fm.findFieldByName("coordinates").castFiniteElement() - elementIter = mesh3d.createElementiterator() - element = elementIter.next() - while element.isValid(): - eft = element.getElementfieldtemplate(coordinates, -1) - if eft.isValid(): - for n in range(eft.getNumberOfLocalNodes()): - node = element.getNode(eft, n + 1) - if lFibrousRingNodeGroup.containsNode(node) or rFibrousRingNodeGroup.containsNode(node): - tmp3dMeshGroup.addElement(element) - break - element = elementIter.next() - tmp2dElementGroup = tmpGroup.getFieldElementGroup(mesh2d) - if tmp2dElementGroup.isValid(): - is_exterior_face_xi2_0 = fm.createFieldAnd(is_exterior, fm.createFieldIsOnFace(Element.FACE_TYPE_XI2_0)) - is_fibrous_ring = fm.createFieldAnd(tmp2dElementGroup, is_exterior_face_xi2_0) - is_left_fibrous_ring = fm.createFieldAnd(is_lam, is_fibrous_ring) - lFibrousRingGroup.getMeshGroup(mesh2d).addElementsConditional(is_left_fibrous_ring) - is_right_fibrous_ring = fm.createFieldAnd(is_ram, is_fibrous_ring) - rFibrousRingGroup.getMeshGroup(mesh2d).addElementsConditional(is_right_fibrous_ring) + if (lFibrousRingGroup.getDimension() <= 0) or (rFibrousRingGroup.getDimension() <= 0): + is_exterior_face_xi2_0 = fm.createFieldAnd(is_exterior, fm.createFieldIsOnFace(Element.FACE_TYPE_XI2_0)) + is_pv = fm.createFieldOr(is_lpv, is_rpv) + if mpvGroup: + is_pv = fm.createFieldOr(is_pv, is_mpv) + lFibrousRingGroup.getMeshGroup(mesh2d).addElementsConditional( + fm.createFieldAnd(is_lam, fm.createFieldAnd(is_exterior_face_xi2_0, fm.createFieldNot(is_pv)))) + is_vc = fm.createFieldOr(is_ivc, is_svc) + rFibrousRingGroup.getMeshGroup(mesh2d).addElementsConditional( + fm.createFieldAnd(is_ram, fm.createFieldAnd(is_exterior_face_xi2_0, fm.createFieldNot(is_vc)))) def getLeftAtriumPulmonaryVeinOstiaElementsCounts(elementsCountAroundLeftAtriumFreeWall, elementsCountOverAtria, commonLeftRightPvOstium): ''' diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_heartventricles1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_heartventricles1.py index 125633e2..3c1ad2bb 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_heartventricles1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_heartventricles1.py @@ -39,11 +39,7 @@ def getParameterSetNames(): 'Human 1', 'Mouse 1', 'Pig 1', - 'Rat 1', - 'Unit Human 1', - 'Unit Mouse 1', - 'Unit Pig 1', - 'Unit Rat 1'] + 'Rat 1'] @staticmethod def getDefaultOptions(parameterSetName='Default'): @@ -51,7 +47,6 @@ def getDefaultOptions(parameterSetName='Default'): isMouse = 'Mouse' in parameterSetName isPig = 'Pig' in parameterSetName isRat = 'Rat' in parameterSetName - notUnitScale = 'Unit' not in parameterSetName options = {} options['Number of elements around LV free wall'] = 5 options['Number of elements around RV free wall'] = 7 @@ -79,12 +74,7 @@ def getDefaultOptions(parameterSetName='Default'): options['Refine number of elements surface'] = 4 options['Refine number of elements through LV wall'] = 1 options['Refine number of elements through wall'] = 1 - if isHuman: - if 'Unit' not in parameterSetName: - options['Unit scale'] = 80.0 - elif isMouse or isRat: - if notUnitScale: - options['Unit scale'] = 5.0 if isMouse else 12.0 + if isMouse or isRat: options['Interventricular sulcus derivative factor'] = 0.8 options['LV outer height'] = 0.9 options['LV outer diameter'] = 0.85 @@ -103,8 +93,6 @@ def getDefaultOptions(parameterSetName='Default'): elif isPig: options['Number of elements up LV apex'] = 1 options['Number of elements up RV'] = 3 - if 'Unit' not in parameterSetName: - options['Unit scale'] = 80.0 options['LV outer height'] = 0.9 options['LV free wall thickness'] = 0.17 options['LV apex thickness'] = 0.07 @@ -115,12 +103,6 @@ def getDefaultOptions(parameterSetName='Default'): options['RV side extension'] = 0.1 options['RV side extension growth factor'] = 0.4 options['Ventricular septum thickness'] = 0.13 - elif 'Rat' in parameterSetName: - if 'Unit' not in parameterSetName: - options['Unit scale'] = 12.0 - options['LV outer height'] = 0.9 - options['LV apex thickness'] = 0.08 - options['RV width'] = 0.35 return options @staticmethod @@ -1239,28 +1221,53 @@ def defineFaceAnnotations(cls, region, options, annotationGroups): """ # 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")) + lvmGroup = getAnnotationGroupForTerm(annotationGroups, get_heart_term("left ventricle myocardium")) + rvmGroup = 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)))) + is_lvm = lvmGroup.getFieldElementGroup(mesh2d) + is_rvm = rvmGroup.getFieldElementGroup(mesh2d) + is_lvm_endo = fm.createFieldAnd(is_lvm, is_exterior_face_xi3_0) + is_rvm_endo = fm.createFieldOr(fm.createFieldAnd(fm.createFieldAnd(is_rvm, is_exterior_face_xi3_0), + fm.createFieldNot(is_lvm_endo)), + fm.createFieldAnd(vSeptumGroup.getFieldElementGroup(mesh2d), + is_exterior_face_xi3_1)) + is_ext_xi3_1_and_not_septum = fm.createFieldAnd( + is_exterior_face_xi3_1, fm.createFieldNot(vSeptumGroup.getFieldElementGroup(mesh2d))) + is_os_lvm = fm.createFieldAnd(is_lvm, is_ext_xi3_1_and_not_septum) + is_os_rvm = fm.createFieldAnd(is_rvm, is_ext_xi3_1_and_not_septum) + + # luminal surfaces of endocardium of left/right ventricle + lslvEndoGroup = findOrCreateAnnotationGroupForTerm( + annotationGroups, region, get_heart_term("luminal surface of left ventricle")) + lslvEndoGroup.getMeshGroup(mesh2d).addElementsConditional(is_lvm_endo) + lsrvEndoGroup = findOrCreateAnnotationGroupForTerm( + annotationGroups, region, get_heart_term("luminal surface of right ventricle")) + lsrvEndoGroup.getMeshGroup(mesh2d).addElementsConditional(is_rvm_endo) + # endocardium groups are defined identically to luminal surfaces at scaffold scale + lvEndoGroup = findOrCreateAnnotationGroupForTerm( + annotationGroups, region, get_heart_term("endocardium of left ventricle")) + lvEndoGroup.getMeshGroup(mesh2d).addElementsConditional(lslvEndoGroup.getFieldElementGroup(mesh2d)) + rvEndoGroup = findOrCreateAnnotationGroupForTerm( + annotationGroups, region, get_heart_term("endocardium of right ventricle")) + rvEndoGroup.getMeshGroup(mesh2d).addElementsConditional(lsrvEndoGroup.getFieldElementGroup(mesh2d)) + + oslvmGroup = findOrCreateAnnotationGroupForTerm( + annotationGroups, region, get_heart_term("outer surface of myocardium of left ventricle")) + oslvmGroup.getMeshGroup(mesh2d).addElementsConditional(is_os_lvm) + osrvmGroup = findOrCreateAnnotationGroupForTerm( + annotationGroups, region, get_heart_term("outer surface of myocardium of right ventricle")) + osrvmGroup.getMeshGroup(mesh2d).addElementsConditional(is_os_rvm) + osmGroup = findOrCreateAnnotationGroupForTerm( + annotationGroups, region, get_heart_term("outer surface of myocardium")) + osmGroup.getMeshGroup(mesh2d).addElementsConditional(oslvmGroup.getFieldElementGroup(mesh2d)) + osmGroup.getMeshGroup(mesh2d).addElementsConditional(osrvmGroup.getFieldElementGroup(mesh2d)) + # if no volumetric epicardium group, add outer surface of myocardium 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) + epiGroup.getMeshGroup(mesh2d).addElementsConditional(osmGroup.getFieldElementGroup(mesh2d)) 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 440009fc..4d5d7928 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_heartventriclesbase1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_heartventriclesbase1.py @@ -46,11 +46,7 @@ def getParameterSetNames(): 'Human 1', 'Mouse 1', 'Pig 1', - 'Rat 1', - 'Unit Human 1', - 'Unit Mouse 1', - 'Unit Pig 1', - 'Unit Rat 1'] + 'Rat 1'] @staticmethod def getDefaultOptions(parameterSetName='Default'): @@ -60,7 +56,6 @@ def getDefaultOptions(parameterSetName='Default'): isMouse = 'Mouse' in parameterSetName isPig = 'Pig' in parameterSetName isRat = 'Rat' in parameterSetName - notUnitScale = 'Unit' not in parameterSetName # only works with particular numbers of elements around options['Number of elements around LV free wall'] = 7 options['Number of elements around RV free wall'] = 7 @@ -316,8 +311,6 @@ def generateBaseMesh(cls, region, options): coordinates = findOrCreateFieldCoordinates(fm) cache = fm.createFieldcache() - mesh = fm.findMeshByDimension(3) - # generate heartventricles1 model to add base plane to annotationGroups = MeshType_3d_heartventricles1.generateBaseMesh(region, options) @@ -330,23 +323,21 @@ def generateBaseMesh(cls, region, options): # 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")) - annotationGroups += [conusArteriosusGroup, lFibrousRingGroup, rFibrousRingGroup] + # temporary groups for making face annotations + lvOutletGroup = AnnotationGroup(region, ("LV outlet", "None")) + mitralAorticCurtainGroup = AnnotationGroup(region, ("Mitral aortic curtain", "None")) + supraventricularCrestGroup = AnnotationGroup(region, ("Supraventricular crest", "None")) + annotationGroups += [conusArteriosusGroup, lFibrousRingGroup, rFibrousRingGroup, lvOutletGroup, + mitralAorticCurtainGroup, supraventricularCrestGroup] # annotation fiducial points markerGroup = findOrCreateFieldGroup(fm, "marker") - 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(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) @@ -764,8 +755,12 @@ def generateBaseMesh(cls, region, options): # Create elements ################# + mesh = fm.findMeshByDimension(3) heartMeshGroup = heartGroup.getMeshGroup(mesh) lvMeshGroup = lvGroup.getMeshGroup(mesh) + lvOutletMeshGroup = lvOutletGroup.getMeshGroup(mesh) + mitralAorticCurtainMeshGroup = mitralAorticCurtainGroup.getMeshGroup(mesh) + supraventricularCrestMeshGroup = supraventricularCrestGroup.getMeshGroup(mesh) rvMeshGroup = rvGroup.getMeshGroup(mesh) vSeptumMeshGroup = vSeptumGroup.getMeshGroup(mesh) conusArteriosusMeshGroup = conusArteriosusGroup.getMeshGroup(mesh) @@ -850,7 +845,6 @@ def generateBaseMesh(cls, region, options): nids = None scalefactors = None meshGroups = [ heartMeshGroup, rvMeshGroup ] - addMarker = None noa = e niv = e @@ -970,7 +964,6 @@ 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 ] @@ -994,13 +987,6 @@ 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) @@ -1112,8 +1098,7 @@ def generateBaseMesh(cls, region, options): eft1 = eft nids = None scalefactors = None - meshGroups = [ heartMeshGroup, lvMeshGroup ] - addMarker = None + meshGroups = [ heartMeshGroup, lvMeshGroup, lvOutletMeshGroup ] eft1 = tricubichermite.createEftNoCrossDerivatives() setEftScaleFactorIds(eft1, [1], []) @@ -1145,7 +1130,7 @@ def generateBaseMesh(cls, region, options): ln_map = [ 1, 2, 3, 3, 4, 5, 6, 6 ] remapEftLocalNodes(eft1, 6, ln_map) elif e == 2: - # 7 node collapsed element where lv outlet ring expands into + # 7 node collapsed element where lv outlet ring expands into LV freewall nids = [ lavNodeId[0][0][2], lavNodeId[0][0][1], lvOutletNodeId[0][4], lvOutletNodeId[0][5], lavNodeId[1][0][2], lavNodeId[1][0][1], lvOutletNodeId[1][4] ] remapEftNodeValueLabel(eft1, [ 1, 2, 5 ], Node.VALUE_LABEL_D_DS1, [ ( Node.VALUE_LABEL_D_DS1, [1] ) ]) @@ -1159,6 +1144,7 @@ def generateBaseMesh(cls, region, options): remapEftNodeValueLabel(eft1, [ 6 ], Node.VALUE_LABEL_D_DS3, [ ( Node.VALUE_LABEL_D_DS3, [1] ) ]) ln_map = [ 1, 2, 3, 4, 5, 6, 7, 6 ] remapEftLocalNodes(eft1, 7, ln_map) + meshGroups.append(mitralAorticCurtainMeshGroup) elif e == 3: # 6 node wedge element bridge/curtain between mitral and aortic valve orifices no = e - 4 @@ -1175,6 +1161,7 @@ def generateBaseMesh(cls, region, options): remapEftNodeValueLabel(eft1, [ 5, 6, 7, 8 ], Node.VALUE_LABEL_D_DS2, []) ln_map = [ 1, 2, 3, 4, 5, 6, 5, 6 ] remapEftLocalNodes(eft1, 6, ln_map) + meshGroups.append(mitralAorticCurtainMeshGroup) elif e == 4: # tetrahedral cfb shim-bridge connector element ni = elementsCountAroundLVFreeWall + elementsCountAroundAtrialSeptum @@ -1197,8 +1184,6 @@ 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, [] ) ]) @@ -1219,13 +1204,6 @@ 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) @@ -1280,6 +1258,7 @@ def generateBaseMesh(cls, region, options): remapEftNodeValueLabel(eft1, [ 8 ], Node.VALUE_LABEL_D_DS3, [ ( Node.VALUE_LABEL_D_DS2, [] ) ]) ln_map = [ 1, 2, 3, 4, 5, 6, 7, 7 ] remapEftLocalNodes(eft1, 7, ln_map) + meshGroups.append(supraventricularCrestMeshGroup) elif e == 2: # 8-node rv crest row 2 element 1 rvin1 = -elementsCountAroundAtrialSeptum - 2 @@ -1300,6 +1279,7 @@ def generateBaseMesh(cls, region, options): remapEftNodeValueLabel(eft1, [ 8 ], Node.VALUE_LABEL_D_DS1, [ ( Node.VALUE_LABEL_D_DS1, [] ), ( Node.VALUE_LABEL_D_DS3, [1] ) ]) remapEftNodeValueLabel(eft1, [ 8 ], Node.VALUE_LABEL_D_DS2, [ ( Node.VALUE_LABEL_D_DS3, [1] ) ]) remapEftNodeValueLabel(eft1, [ 8 ], Node.VALUE_LABEL_D_DS3, [ ( Node.VALUE_LABEL_D_DS2, [] ) ]) + meshGroups.append(supraventricularCrestMeshGroup) elif e == 3: # 8-node wedge rv crest row 2 element 2 rvin1 = -elementsCountAroundAtrialSeptum - 2 @@ -1318,7 +1298,7 @@ def generateBaseMesh(cls, region, options): remapEftNodeValueLabel(eft1, [ 7 ], Node.VALUE_LABEL_D_DS1, [ ( Node.VALUE_LABEL_D_DS1, [] ), ( Node.VALUE_LABEL_D_DS3, [] ) ]) remapEftNodeValueLabel(eft1, [ 7 ], Node.VALUE_LABEL_D_DS2, [ ( Node.VALUE_LABEL_D_DS3, [1] ) ]) remapEftNodeValueLabel(eft1, [ 7 ], Node.VALUE_LABEL_D_DS3, [ ( Node.VALUE_LABEL_D_DS2, [] ) ]) - meshGroups += [ conusArteriosusMeshGroup ] + meshGroups += [conusArteriosusMeshGroup, supraventricularCrestMeshGroup] elif e == 4: # 8-node rv crest inner 4 by rv outlet rvin1 = -elementsCountAroundAtrialSeptum - 2 @@ -1347,7 +1327,7 @@ def generateBaseMesh(cls, region, options): remapEftNodeValueLabel(eft1, [ 7 ], Node.VALUE_LABEL_D2_DS1DS2, [ ( Node.VALUE_LABEL_D_DS2, [] ) ]) remapEftNodeValueLabel(eft1, [ 8 ], Node.VALUE_LABEL_D_DS2, [ ( Node.VALUE_LABEL_D_DS1, [1] ) ]) remapEftNodeValueLabel(eft1, [ 8 ], Node.VALUE_LABEL_D_DS1, [ ( Node.VALUE_LABEL_D_DS2, [] ) ]) - meshGroups += [ conusArteriosusMeshGroup ] + meshGroups += [conusArteriosusMeshGroup, supraventricularCrestMeshGroup] else: continue @@ -1442,61 +1422,44 @@ def defineFaceAnnotations(cls, region, options, annotationGroups): New face annotation groups are appended to this list. """ MeshType_3d_heartventricles1.defineFaceAnnotations(region, options, annotationGroups) - - fm = region.getFieldmodule() + conusArteriosusGroup = getAnnotationGroupForTerm(annotationGroups, get_heart_term("conus arteriosus")) lFibrousRingGroup = getAnnotationGroupForTerm(annotationGroups, get_heart_term("left fibrous ring")) rFibrousRingGroup = getAnnotationGroupForTerm(annotationGroups, get_heart_term("right fibrous ring")) - if (lFibrousRingGroup.getDimension() == 0) or (lFibrousRingGroup.getDimension() == 0): - # not already added by full heart scaffold - nodes = fm.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES) - lFibrousRingNodeGroup = lFibrousRingGroup.getNodesetGroup(nodes) - rFibrousRingNodeGroup = rFibrousRingGroup.getNodesetGroup(nodes) - # make temp group containing all elements, faces etc. if have any nodes in fibrous ring groups - tmpGroup = fm.createFieldGroup() - tmpGroup.setSubelementHandlingMode(FieldGroup.SUBELEMENT_HANDLING_MODE_FULL) - mesh3d = fm.findMeshByDimension(3) - tmp3dMeshGroup = tmpGroup.createFieldElementGroup(mesh3d).getMeshGroup() - coordinates = fm.findFieldByName("coordinates").castFiniteElement() - elementIter = mesh3d.createElementiterator() - element = elementIter.next() - while element.isValid(): - eft = element.getElementfieldtemplate(coordinates, -1) - if eft.isValid(): - for n in range(eft.getNumberOfLocalNodes()): - node = element.getNode(eft, n + 1) - if lFibrousRingNodeGroup.containsNode(node) or rFibrousRingNodeGroup.containsNode(node): - tmp3dMeshGroup.addElement(element) - break - element = elementIter.next() + # temporary groups + lvOutletGroup = getAnnotationGroupForTerm(annotationGroups, ("LV outlet", "None")) + mitralAorticCurtainGroup = getAnnotationGroupForTerm(annotationGroups, ("Mitral aortic curtain", "None")) + supraventricularCrestGroup = getAnnotationGroupForTerm(annotationGroups, ("Supraventricular crest", "None")) + + fm = region.getFieldmodule() + if (lFibrousRingGroup.getDimension() <= 0) or (rFibrousRingGroup.getDimension() <= 0): mesh2d = fm.findMeshByDimension(2) - tmp2dElementGroup = tmpGroup.getFieldElementGroup(mesh2d) - if tmp2dElementGroup.isValid(): - lvGroup = getAnnotationGroupForTerm(annotationGroups, get_heart_term("left ventricle myocardium")) - rvGroup = getAnnotationGroupForTerm(annotationGroups, get_heart_term("right ventricle myocardium")) - is_lv = lvGroup.getFieldElementGroup(mesh2d) - is_rv = rvGroup.getFieldElementGroup(mesh2d) - is_exterior = fm.createFieldIsExterior() - is_face_xi1_0 = fm.createFieldIsOnFace(Element.FACE_TYPE_XI1_0) - is_face_xi2_0 = fm.createFieldIsOnFace(Element.FACE_TYPE_XI2_0) - is_face_xi2_1 = fm.createFieldIsOnFace(Element.FACE_TYPE_XI2_1) - is_face_xi1_0_or_xi2_1 = fm.createFieldOr(is_face_xi1_0, is_face_xi2_1) - is_exterior_face_xi1_0_or_xi2_1 = fm.createFieldAnd(is_exterior, is_face_xi1_0_or_xi2_1) - is_face_xi2_0_or_xi2_1 = fm.createFieldOr(is_face_xi2_0, is_face_xi2_1) - is_exterior_face_xi2_0_or_xi2_1 = fm.createFieldAnd(is_exterior, is_face_xi2_0_or_xi2_1) - is_left_fibrous_ring = fm.createFieldAnd( - is_lv, fm.createFieldAnd(tmp2dElementGroup, is_exterior_face_xi2_0_or_xi2_1)) - lFibrousRingGroup.getMeshGroup(mesh2d).addElementsConditional(is_left_fibrous_ring) - is_right_fibrous_ring = fm.createFieldAnd( - is_rv, fm.createFieldAnd(tmp2dElementGroup, is_exterior_face_xi1_0_or_xi2_1)) - rFibrousRingGroup.getMeshGroup(mesh2d).addElementsConditional(is_right_fibrous_ring) - # hacky correction around LV outflow - size = tmp3dMeshGroup.getSize() - elementIter = tmp3dMeshGroup.createElementiterator() - identifiers = [] - for i in range(size - 8): - element = elementIter.next() - identifiers.append(element.getIdentifier()) - for identifier in identifiers: - tmp3dMeshGroup.removeElement(mesh3d.findElementByIdentifier(identifier)) - is_lv_outflow = fm.createFieldAnd(tmp2dElementGroup, is_face_xi2_1) - lFibrousRingGroup.getMeshGroup(mesh2d).removeElementsConditional(is_lv_outflow) + is_exterior = fm.createFieldIsExterior() + is_face_xi1_0 = fm.createFieldIsOnFace(Element.FACE_TYPE_XI1_0) + is_face_xi2_0 = fm.createFieldIsOnFace(Element.FACE_TYPE_XI2_0) + is_face_xi2_1 = fm.createFieldIsOnFace(Element.FACE_TYPE_XI2_1) + lvGroup = getAnnotationGroupForTerm(annotationGroups, get_heart_term("left ventricle myocardium")) + rvGroup = getAnnotationGroupForTerm(annotationGroups, get_heart_term("right ventricle myocardium")) + is_lv = lvGroup.getFieldElementGroup(mesh2d) + is_rv = rvGroup.getFieldElementGroup(mesh2d) + is_lvo = lvOutletGroup.getFieldElementGroup(mesh2d) + is_mac = mitralAorticCurtainGroup.getFieldElementGroup(mesh2d) + is_lfr = fm.createFieldAnd( + is_exterior, + fm.createFieldOr( + fm.createFieldAnd(fm.createFieldAnd(is_lv, fm.createFieldNot(is_lvo)), is_face_xi2_1), + fm.createFieldAnd(is_mac, is_face_xi2_0))) + lFibrousRingGroup.getMeshGroup(mesh2d).addElementsConditional(is_lfr) + is_ca = conusArteriosusGroup.getFieldElementGroup(mesh2d) + is_svc = supraventricularCrestGroup.getFieldElementGroup(mesh2d) + is_rfr = fm.createFieldAnd( + is_exterior, + fm.createFieldAnd( + is_rv, + fm.createFieldOr( + fm.createFieldAnd(is_face_xi2_1, fm.createFieldNot(is_ca)), + fm.createFieldAnd(is_face_xi1_0, is_svc)))) + rFibrousRingGroup.getMeshGroup(mesh2d).addElementsConditional(is_rfr) + + annotationGroups.remove(lvOutletGroup) + annotationGroups.remove(mitralAorticCurtainGroup) + annotationGroups.remove(supraventricularCrestGroup) diff --git a/tests/test_general.py b/tests/test_general.py index 17f0e57a..b690c4c2 100644 --- a/tests/test_general.py +++ b/tests/test_general.py @@ -109,7 +109,7 @@ def test_user_annotation_groups(self): annotationGroups = scaffoldPackage.getAnnotationGroups() self.assertEqual(24, len(annotationGroups)) - endocardium_of_la = scaffoldPackage.findAnnotationGroupByName('endocardium of left atrium') + endocardium_of_la = scaffoldPackage.findAnnotationGroupByName('left atrium endocardium') self.assertTrue(isinstance(endocardium_of_la, AnnotationGroup)) self.assertFalse(scaffoldPackage.isUserAnnotationGroup(endocardium_of_la)) self.assertFalse(scaffoldPackage.deleteAnnotationGroup(endocardium_of_la)) # can't delete auto annotation groups diff --git a/tests/test_heart.py b/tests/test_heart.py index fb99cac4..430dff84 100644 --- a/tests/test_heart.py +++ b/tests/test_heart.py @@ -22,12 +22,12 @@ def test_heart1(self): """ 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" ]); + self.assertEqual(parameterSetNames, ["Default", "Human 1", "Mouse 1", "Pig 1", "Rat 1"]); options = scaffold.getDefaultOptions("Human 1") self.assertEqual(123, len(options)) self.assertEqual(0.9, options.get("LV outer height")) - self.assertEqual(80.0, options.get("Unit scale")) + self.assertEqual(1.0, options.get("Unit scale")) + options["Unit scale"] = 80.0 self.assertEqual(7, options.get("Number of elements around LV free wall")) self.assertEqual(7, options.get("Number of elements around RV free wall")) # simplify atria @@ -39,9 +39,22 @@ def test_heart1(self): context = Context("Test") region = context.getDefaultRegion() self.assertTrue(region.isValid()) - annotationGroups = scaffold.generateMesh(region, options) - self.assertEqual(32, len(annotationGroups)) fieldmodule = region.getFieldmodule() + + # Need to do the following manually to save originalAnnotationGroups which has some temporary groups + # annotationGroups = scaffold.generateMesh(region, options) + with ChangeManager(fieldmodule): + annotationGroups = scaffold.generateBaseMesh(region, options) + fieldmodule.defineAllFaces() + originalAnnotationGroups = copy.copy(annotationGroups) + for annotationGroup in annotationGroups: + annotationGroup.addSubelements() + scaffold.defineFaceAnnotations(region, options, annotationGroups) + for annotationGroup in annotationGroups: + if annotationGroup not in originalAnnotationGroups: + annotationGroup.addSubelements() + + self.assertEqual(45, len(annotationGroups)) mesh3d = fieldmodule.findMeshByDimension(3) self.assertEqual(332, mesh3d.getSize()) mesh2d = fieldmodule.findMeshByDimension(2) @@ -49,7 +62,7 @@ def test_heart1(self): mesh1d = fieldmodule.findMeshByDimension(1) self.assertEqual(1567, mesh1d.getSize()) nodes = fieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES) - self.assertEqual(614, nodes.getSize()) + self.assertEqual(611, nodes.getSize()) datapoints = fieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_DATAPOINTS) self.assertEqual(0, datapoints.getSize()) @@ -93,8 +106,8 @@ def test_heart1(self): expectedSizes2d = { "endocardium of left ventricle" : 88, "endocardium of right ventricle" : 73, - "endocardium of left atrium" : 82, - "endocardium of right atrium" : 74, + "left atrium endocardium" : 82, + "right atrium endocardium" : 74, "epicardium" : 229 } for name in expectedSizes2d: @@ -105,7 +118,7 @@ def test_heart1(self): # test finding a marker in scaffold markerGroup = fieldmodule.findFieldByName("marker").castGroup() markerNodes = markerGroup.getFieldNodeGroup(nodes).getNodesetGroup() - self.assertEqual(7, markerNodes.getSize()) + self.assertEqual(4, markerNodes.getSize()) markerName = fieldmodule.findFieldByName("marker_name") self.assertTrue(markerName.isValid()) markerLocation = fieldmodule.findFieldByName("marker_location") @@ -122,14 +135,9 @@ def test_heart1(self): self.assertTrue(apexGroup.getNodesetGroup(nodes).containsNode(node)) # refine 2x2x2 and check result - # first remove any face (but not point) annotation groups as they are re-added by defineFaceAnnotations - removeAnnotationGroups = [] - for annotationGroup in annotationGroups: - if (not annotationGroup.hasMeshGroup(mesh3d)) and (annotationGroup.hasMeshGroup(mesh2d) or annotationGroup.hasMeshGroup(mesh1d)): - removeAnnotationGroups.append(annotationGroup) - for annotationGroup in removeAnnotationGroups: - annotationGroups.remove(annotationGroup) - self.assertEqual(23, len(annotationGroups)) + # need to use original annotation groups to get temporaries + annotationGroups = originalAnnotationGroups + self.assertEqual(26, len(annotationGroups)) # 23 + 3 temporary groups refineRegion = region.createRegion() refineFieldmodule = refineRegion.getFieldmodule() @@ -148,7 +156,7 @@ def test_heart1(self): for annotation in annotationGroups: if annotation not in oldAnnotationGroups: annotationGroup.addSubelements() - self.assertEqual(32, len(annotationGroups)) + self.assertEqual(45, len(annotationGroups)) mesh3d = refineFieldmodule.findMeshByDimension(3) self.assertEqual(2580, mesh3d.getSize()) @@ -157,7 +165,7 @@ def test_heart1(self): mesh1d = refineFieldmodule.findMeshByDimension(1) self.assertEqual(9983, mesh1d.getSize()) nodes = refineFieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES) - self.assertEqual(3693, nodes.getSize()) + self.assertEqual(3690, nodes.getSize()) datapoints = refineFieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_DATAPOINTS) self.assertEqual(0, datapoints.getSize()) @@ -181,7 +189,7 @@ def test_heart1(self): markerGroup = refineFieldmodule.findFieldByName("marker").castGroup() refinedNodes = refineFieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES) markerNodes = markerGroup.getFieldNodeGroup(refinedNodes).getNodesetGroup() - self.assertEqual(7, markerNodes.getSize()) + self.assertEqual(4, markerNodes.getSize()) markerName = refineFieldmodule.findFieldByName("marker_name") self.assertTrue(markerName.isValid()) markerLocation = refineFieldmodule.findFieldByName("marker_location")