From a81a7dc42c1c63ac8a342621f4c023c17b32d257 Mon Sep 17 00:00:00 2001 From: Richard Christie Date: Sat, 19 Mar 2022 17:27:06 +1300 Subject: [PATCH 1/3] Add user marker point annotation groups Use simpler marker annotation API in brainstem1 and heartventricles1. --- .../annotation/annotationgroup.py | 302 +++++++++++++++++- .../meshtypes/meshtype_3d_brainstem.py | 152 ++++----- .../meshtypes/meshtype_3d_heartventricles1.py | 20 +- src/scaffoldmaker/scaffoldpackage.py | 32 +- src/scaffoldmaker/utils/zinc_utils.py | 35 +- tests/test_general.py | 160 +++++++++- 6 files changed, 572 insertions(+), 129 deletions(-) diff --git a/src/scaffoldmaker/annotation/annotationgroup.py b/src/scaffoldmaker/annotation/annotationgroup.py index 1fa8ab75..ce25f9fe 100644 --- a/src/scaffoldmaker/annotation/annotationgroup.py +++ b/src/scaffoldmaker/annotation/annotationgroup.py @@ -1,12 +1,17 @@ """ Describes subdomains of a scaffold with attached names and terms. """ +import sys from opencmiss.utils.zinc.general import ChangeManager -from opencmiss.zinc.field import Field, FieldGroup +from opencmiss.utils.zinc.field import find_or_create_field_coordinates, find_or_create_field_group, \ + find_or_create_field_stored_mesh_location, find_or_create_field_stored_string +from opencmiss.zinc.element import Element +from opencmiss.zinc.field import Field, FieldFiniteElement, FieldGroup +from opencmiss.zinc.node import Node from opencmiss.zinc.result import RESULT_OK -from scaffoldmaker.utils.zinc_utils import group_get_highest_dimension, \ - identifier_ranges_from_string, identifier_ranges_to_string, \ +from scaffoldmaker.utils.zinc_utils import get_highest_dimension_mesh, get_next_unused_node_identifier, \ + group_get_highest_dimension, identifier_ranges_from_string, identifier_ranges_to_string, \ mesh_group_add_identifier_ranges, mesh_group_to_identifier_ranges, \ nodeset_group_add_identifier_ranges, nodeset_group_to_identifier_ranges @@ -34,6 +39,24 @@ def __init__(self, region, term): self._group = fieldmodule.createFieldGroup() self._group.setName(self._name) self._group.setManaged(True) + self._isMarker = False + self._markerMaterialCoordinatesField = None + + def clear(self): + """ + Manual clean up of user defined annotation group. See use in ScaffoldPackage. + If this is a marker, destroys the marker node. + Finally clears group + """ + fieldmodule = self._group.getFieldmodule() + with ChangeManager(fieldmodule): + if self._isMarker: + nodes = fieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES) + markerNode = self.getNodesetGroup(nodes).createNodeiterator().next() + nodes.destroyNode(markerNode) + self._isMarker = False + self._markerMaterialCoordinatesField = None + self._group.clear() def toDict(self): ''' @@ -62,8 +85,28 @@ def toDict(self): 'dimension' : dimension, 'identifierRanges' : identifier_ranges_to_string(identifierRanges) } + if self._isMarker: + dct['marker'] = self._markerToDict() return dct + def _markerToDict(self): + """ + If the node is in the marker group, return a dict mapping "marker_location" to [elementId, xi] and if another + entry is present it is a map from the materialCoordinatesField name to the coordinates value. Eg. + { + "heart coordinates": [11.5, 16.2, 7.6], + "marker_location": [5, [0.0, 0.0, 1.0]] + } + :return: A python dict as above or None if not a marker point. + """ + assert self._isMarker + markerDct = {} + if self._markerMaterialCoordinatesField: + markerDct[self._markerMaterialCoordinatesField.getName()] = self.getMarkerMaterialCoordinates()[1] + element, xi = self.getMarkerLocation() + markerDct['marker_location'] = [element.getIdentifier(), xi] + return markerDct + @classmethod def fromDict(cls, dct, region): ''' @@ -76,6 +119,7 @@ def fromDict(cls, dct, region): ontId = dct['ontId'] dimension = dct['dimension'] identifierRangesString = dct['identifierRanges'] + markerDct = dct.get('marker') identifierRanges = identifier_ranges_from_string(identifierRangesString) fieldmodule = region.getFieldmodule() with ChangeManager(fieldmodule): @@ -85,10 +129,245 @@ def fromDict(cls, dct, region): meshGroup = annotationGroup.getMeshGroup(fieldmodule.findMeshByDimension(dimension)) mesh_group_add_identifier_ranges(meshGroup, identifierRanges) else: - nodesetGroup = annotationGroup.getNodesetGroup(fieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES)) + if markerDct: + nodeIdentifier = int(identifierRanges[0][0]) + annotationGroup._markerFromDict(nodeIdentifier, markerDct) + nodesetGroup = annotationGroup.getNodesetGroup( + fieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES)) nodeset_group_add_identifier_ranges(nodesetGroup, identifierRanges) return annotationGroup + def _markerFromDict(self, nodeIdentifier, markerDct: dict): + """ + Define a new marker point node with nodeIdentifier and fields defined as in markerDct. + Warn but do nothing if node with that identifier exists. + :param nodeIdentifier: + :param markerDct: A dict mapping names of fields to values, e.g. + { + "heart coordinates": [11.5, 16.2, 7.6], + "marker_location": [5, [0.0, 0.0, 1.0]] + } + """ + assert not self._isMarker + fieldmodule = self._group.getFieldmodule() + nodes = fieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES) + node = nodes.findNodeByIdentifier(nodeIdentifier) + if node.isValid(): + print("Error: Read annotation group " + self._name + ". "\ + "Marker point node " + str(nodeIdentifier) + " already exists", file=sys.stderr) + return + mesh = get_highest_dimension_mesh(fieldmodule) + if not mesh: + print("Error: Read annotation group " + self._name + ". Marker cannot be made for empty mesh", + file=sys.stderr) + return + with ChangeManager(fieldmodule): + materialCoordinatesField = None + materialCoordinates = None + element = None + xi = None + for key in markerDct.keys(): + if key != "marker_location": + fieldName = key + materialCoordinatesField = fieldmodule.findFieldByName(fieldName).castFiniteElement() + materialCoordinates = markerDct[fieldName] + if materialCoordinatesField.isValid() and isinstance(materialCoordinates, list): + break; + print("Error: Read annotation group " + self._name + ". " \ + "Invalid marker material coordinates field " + fieldName + " or value.", file=sys.stderr) + materialCoordinatesField = None + materialCoordinates = None + break + if not materialCoordinatesField: + element_xi = markerDct.get("marker_location") + if not element_xi: + print("Error: Read annotation group " + self._name + ". " \ + "Marker missing marker_location entry", file=sys.stderr) + return + elementIdentifier, xi = element_xi + element = mesh.findElementByIdentifier(elementIdentifier) + if not element.isValid(): + print("Error: Read annotation group " + self._name + ". " \ + "Marker element " + str(elementIdentifier) + " not found", file=sys.stderr) + return + self.createMarkerNode(nodeIdentifier, materialCoordinatesField, materialCoordinates, element, xi) + + def isMarker(self): + """ + Query if this annotation group is a created marker node. + """ + return self._isMarker + + def createMarkerNode(self, startNodeIdentifier=1, + materialCoordinatesField: FieldFiniteElement=None, materialCoordinates=None, + element=None, xi=[0.0, 0.0, 0.0]): + """ + Convert annotation group into a marker point annotation. + Important: annotation group must currently be empty, and elements must exist. + The marker node is added to the marker group in addition to this group. + Raises an exception if marker creation cannot be achieved. + :param startNodeIdentifier: First unused node identifier >= this is uses for marker node. + :param materialCoordinatesField: Material coordinates field to define location of marker point in. + Must be a finite element type field for which isTypeCoordinates() is True, with up to 3 components, + and at least as many components as the highest mesh dimension. + Only one of materialCoordinatesField or element may be specified. + :param materialCoordinates: The coordinates to assign to the materialCoordinatesField, if field supplied. + The element:xi coordinates are computed from it. + :param element: Optional element to set initial location from; must be in highest dimension mesh. + Only one of materialCoordinatesField or element may be specified. + If neither is specified the first element in the highest dimension mesh is used. + :param xi: If element is supplied or first is assumed, the xi coordinates to embed at. + :return: Zinc Node representing marker point, with either default location in first element or that supplied. + """ + assert not self._isMarker + assert self._group.isEmpty(), "Can only create marker in empty annotation group" + fieldmodule = self._group.getFieldmodule() + mesh = get_highest_dimension_mesh(fieldmodule) + assert mesh, "Can only create marker point if there is a mesh with elements" + assert not (element and materialCoordinatesField) + if materialCoordinatesField: + coordinatesCount = materialCoordinatesField.getNumberOfComponents() + assert materialCoordinates and (len(materialCoordinates) >= coordinatesCount) and \ + materialCoordinatesField.castFiniteElement().isValid() and \ + materialCoordinatesField.isTypeCoordinate() and (mesh.getDimension() <= coordinatesCount <= 3) + nodes = fieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES) + nodeIdentifier = get_next_unused_node_identifier(nodes, startNodeIdentifier) + with ChangeManager(fieldmodule): + markerGroup = find_or_create_field_group(fieldmodule, "marker") + markerNodeGroup = markerGroup.getFieldNodeGroup(nodes) + if not markerNodeGroup.isValid(): + markerNodeGroup = markerGroup.createFieldNodeGroup(nodes) + markerNodesetGroup = markerNodeGroup.getNodesetGroup() + markerLocation = find_or_create_field_stored_mesh_location(fieldmodule, mesh, name="marker_location") + markerName = find_or_create_field_stored_string(fieldmodule, name="marker_name") + nodetemplate = nodes.createNodetemplate() + if materialCoordinatesField: + assert RESULT_OK == nodetemplate.defineField(materialCoordinatesField) + assert RESULT_OK == nodetemplate.defineField(markerLocation) + assert RESULT_OK == nodetemplate.defineField(markerName) + markerNode = nodes.createNode(nodeIdentifier, nodetemplate) + assert RESULT_OK == self.getNodesetGroup(nodes).addNode(markerNode) + assert RESULT_OK == markerNodesetGroup.addNode(markerNode) + self._isMarker = True + # assign marker name to be same as this group's name. This needs to be maintained. See setName() + fieldcache = fieldmodule.createFieldcache() + fieldcache.setNode(markerNode) + markerName.assignString(fieldcache, self._name) + if materialCoordinatesField: + self.setMarkerMaterialCoordinates(materialCoordinatesField, materialCoordinates) + else: + if not element: + element = mesh.createElementiterator().next() + self.setMarkerLocation(element, xi) + return markerNode + + def getMarkerLocation(self): + """ + If the annotation is a created marker point, get its element:xi location. + :return: Zinc Element, xi (list of float) + """ + assert self._isMarker + fieldmodule = self._group.getFieldmodule() + mesh = get_highest_dimension_mesh(fieldmodule) + markerLocation = find_or_create_field_stored_mesh_location(fieldmodule, mesh, name="marker_location") + fieldcache = fieldmodule.createFieldcache() + nodes = fieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES) + markerNode = self.getNodesetGroup(nodes).createNodeiterator().next() + fieldcache.setNode(markerNode) + element, xi = markerLocation.evaluateMeshLocation(fieldcache, mesh.getDimension()) + return element, xi + + def setMarkerLocation(self, element: Element, xi: list): + """ + If the annotation group is a created marker point, set its element:xi location. + If the marker also has a materialCoordinatesField, it is updated to the value + at the supplied element:xi location. + :param element: Element to set location in + :param xi: xi coordinates (list of float) + """ + assert self._isMarker + fieldmodule = self._group.getFieldmodule() + mesh = get_highest_dimension_mesh(fieldmodule) + assert mesh.containsElement(element), "Invalid element, not in highest dimension mesh" + assert len(xi) >= mesh.getDimension() + with ChangeManager(fieldmodule): + markerLocation = find_or_create_field_stored_mesh_location(fieldmodule, mesh, name="marker_location") + fieldcache = fieldmodule.createFieldcache() + nodes = fieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES) + markerNode = self.getNodesetGroup(nodes).createNodeiterator().next() + fieldcache.setNode(markerNode) + markerLocation.assignMeshLocation(fieldcache, element, xi) + if self._markerMaterialCoordinatesField: + fieldcache.setMeshLocation(element, xi) + result, materialCoordinates = self._markerMaterialCoordinatesField.evaluateReal( + fieldcache, self._markerMaterialCoordinatesField.getNumberOfComponents()) + fieldcache.setNode(markerNode) + self._markerMaterialCoordinatesField.assignReal(fieldcache, materialCoordinates) + + def getMarkerMaterialCoordinates(self): + """ + If the annotation is a created marker point, get its material coordinates field and location, if set. + :return: materialCoordinatesField, materialCoordinates or None, None + """ + assert self._isMarker + if not self._markerMaterialCoordinatesField: + return None, None + fieldmodule = self._group.getFieldmodule() + fieldcache = fieldmodule.createFieldcache() + nodes = fieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES) + markerNode = self.getNodesetGroup(nodes).createNodeiterator().next() + fieldcache.setNode(markerNode) + result, materialCoordinates = self._markerMaterialCoordinatesField.evaluateReal( + fieldcache, self._markerMaterialCoordinatesField.getNumberOfComponents()) + return self._markerMaterialCoordinatesField, materialCoordinates + + def setMarkerMaterialCoordinates(self, materialCoordinatesField, materialCoordinates=None): + """ + Also updates the marker location when this is assigned, forcing it to be within the mesh. + The material coordinates are then recalculated to be in the mesh. + Some approximations may occur if the point is outside the mesh - user beware. + :param materialCoordinatesField: Material coordinates field to set or change to, or None to remove + material coordinates field so marker point only has an element:xi location. Must be defined on elements + of highest dimension mesh. + :param materialCoordinates: If None, evaluate materialCoordinatesField at current marker location. + """ + assert self._isMarker + if not (self._markerMaterialCoordinatesField or materialCoordinatesField): + return + fieldmodule = self._group.getFieldmodule() + mesh = get_highest_dimension_mesh(fieldmodule) + nodes = fieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES) + markerNode = self.getNodesetGroup(nodes).createNodeiterator().next() + with ChangeManager(fieldmodule): + oldMaterialCoordinatesField = self._markerMaterialCoordinatesField + if materialCoordinatesField != oldMaterialCoordinatesField: + nodetemplate = nodes.createNodetemplate() + if self._markerMaterialCoordinatesField: + assert RESULT_OK == nodetemplate.undefineField(self._markerMaterialCoordinatesField) + if materialCoordinatesField: + assert RESULT_OK == nodetemplate.defineField(materialCoordinatesField) + assert RESULT_OK == markerNode.merge(nodetemplate) + self._markerMaterialCoordinatesField = materialCoordinatesField + if materialCoordinatesField: + if materialCoordinates is not None: + # find nearest location in mesh + constCoordinates = fieldmodule.createFieldConstant(materialCoordinates) + findMeshLocation = fieldmodule.createFieldFindMeshLocation( + constCoordinates, materialCoordinatesField, mesh) + fieldcache = fieldmodule.createFieldcache() + fieldcache.setNode(markerNode) + element, xi = findMeshLocation.evaluateMeshLocation(fieldcache, mesh.getDimension()) + if not element.isValid(): + self.setMarkerMaterialCoordinates(oldMaterialCoordinatesField) + print("AnnotationGroup setMarkerMaterialCoordinates. Field is not defined on mesh. Reverting.") + return + del findMeshLocation + del constCoordinates + else: + element, xi = self.getMarkerLocation() + # use this function to reassign material coordinates to be within mesh + self.setMarkerLocation(element, xi) + def getName(self): return self._name @@ -99,6 +378,8 @@ def setName(self, name): as the name is already in use. :return: True on success, otherwise False ''' + if self._name == name: + return True fieldmodule = self._group.getFieldmodule() # use ChangeManager so multiple name changes are atomic with ChangeManager(fieldmodule): @@ -114,6 +395,13 @@ def setName(self, name): if nodeGroup.isValid(): nodeGroup.setName(name + '.' + nodes.getName()) self._name = name + if self._isMarker: + # assign marker name to be same as this group's name. This needs to be maintained + markerName = find_or_create_field_stored_string(fieldmodule, name="marker_name") + fieldcache = fieldmodule.createFieldcache() + markerNode = self.getNodesetGroup(nodes).createNodeiterator().next() + fieldcache.setNode(markerNode) + markerName.assignString(fieldcache, self._name) return True return False @@ -239,7 +527,8 @@ def findOrCreateAnnotationGroupForTerm(annotationGroups: list, region, term) -> 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 + "'" + 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) @@ -256,7 +545,8 @@ def getAnnotationGroupForTerm(annotationGroups: list, term) -> 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 + "'" + 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/meshtypes/meshtype_3d_brainstem.py b/src/scaffoldmaker/meshtypes/meshtype_3d_brainstem.py index 1fe045a4..eb04c297 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_brainstem.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_brainstem.py @@ -11,6 +11,7 @@ from opencmiss.zinc.field import FieldFindMeshLocation from opencmiss.utils.zinc.field import Field, findOrCreateFieldCoordinates, findOrCreateFieldGroup, findOrCreateFieldNodeGroup, findOrCreateFieldStoredMeshLocation, \ findOrCreateFieldStoredString +from opencmiss.utils.zinc.general import ChangeManager from opencmiss.utils.zinc.finiteelement import getMaximumNodeIdentifier from scaffoldmaker.annotation.annotationgroup import AnnotationGroup, findOrCreateAnnotationGroupForTerm from scaffoldmaker.annotation.brainstem_terms import get_brainstem_term @@ -37,7 +38,7 @@ class MeshType_3d_brainstem1(Scaffold_base): 'D2 derivatives': True, 'D3 derivatives': True, 'Length': 3.0, - 'Number of elements': 6 + 'Number of elements': 8 }, 'meshEdits': exnodeStringFromNodeValues( # dimensional. [Node.VALUE_LABEL_VALUE, Node.VALUE_LABEL_D_DS1, Node.VALUE_LABEL_D_DS2, Node.VALUE_LABEL_D2_DS1DS2, @@ -56,21 +57,21 @@ class MeshType_3d_brainstem1(Scaffold_base): { '_AnnotationGroup': True, 'dimension': 1, - 'identifierRanges': '1', + 'identifierRanges': '1-3', 'name': get_brainstem_term('medulla oblongata')[0], 'ontId': get_brainstem_term('medulla oblongata')[1] }, { '_AnnotationGroup': True, 'dimension': 1, - 'identifierRanges': '2', + 'identifierRanges': '4-6', 'name': get_brainstem_term('pons')[0], 'ontId': get_brainstem_term('pons')[1] }, { '_AnnotationGroup': True, 'dimension': 1, - 'identifierRanges': '3', + 'identifierRanges': '7-8', 'name': get_brainstem_term('midbrain')[0], 'ontId': get_brainstem_term('midbrain')[1] }] @@ -103,21 +104,21 @@ class MeshType_3d_brainstem1(Scaffold_base): { '_AnnotationGroup': True, 'dimension': 1, - 'identifierRanges': '1', + 'identifierRanges': '1-2', 'name': get_brainstem_term('medulla oblongata')[0], 'ontId': get_brainstem_term('medulla oblongata')[1] }, { '_AnnotationGroup': True, 'dimension': 1, - 'identifierRanges': '2', + 'identifierRanges': '3-4', 'name': get_brainstem_term('pons')[0], 'ontId': get_brainstem_term('pons')[1] }, { '_AnnotationGroup': True, 'dimension': 1, - 'identifierRanges': '3', + 'identifierRanges': '5-6', 'name': get_brainstem_term('midbrain')[0], 'ontId': get_brainstem_term('midbrain')[1] }] @@ -501,37 +502,18 @@ def generateBaseMesh(cls, region, options): annotationGroupAlong = [[brainstemGroup, midbrainGroup], [brainstemGroup, ponsGroup], [brainstemGroup, medullaGroup]] - - # point markers # centralCanal = findOrCreateAnnotationGroupForTerm(annotationGroups, region, # get_brainstem_term('central canal of spinal cord')) # cerebralAqueduct = findOrCreateAnnotationGroupForTerm(annotationGroups, region, # get_brainstem_term('cerebral aqueduct')) # foramenCaecum = findOrCreateAnnotationGroupForTerm(annotationGroups, region, - # get_brainstem_term('foramen caecum of medulla oblongata')) - dorsalMidCaudalGroup = findOrCreateAnnotationGroupForTerm(annotationGroups, region, - get_brainstem_term('brainstem dorsal midline caudal point')) - ventralMidCaudalGroup = findOrCreateAnnotationGroupForTerm(annotationGroups, region, - get_brainstem_term('brainstem ventral midline caudal point')) - dorsalMidCranGroup = findOrCreateAnnotationGroupForTerm(annotationGroups, region, - get_brainstem_term('brainstem dorsal midline cranial point')) - ventralMidCranGroup = findOrCreateAnnotationGroupForTerm(annotationGroups, region, - get_brainstem_term('brainstem ventral midline cranial point')) - dorsalMidMedullaPonsJunction = findOrCreateAnnotationGroupForTerm(annotationGroups, region, - get_brainstem_term('brainstem dorsal midline pons-medulla junction')) - ventralMidMedullaPonsJunction = findOrCreateAnnotationGroupForTerm(annotationGroups, region, - get_brainstem_term('brainstem ventral midline pons-medulla junction')) - dorsalMidMidbrainPonsJunction = findOrCreateAnnotationGroupForTerm(annotationGroups, region, - get_brainstem_term('brainstem dorsal midline midbrain-pons junction')) - ventralMidMidbrainPonsJunction = findOrCreateAnnotationGroupForTerm(annotationGroups, region, - get_brainstem_term('brainstem ventral midline midbrain-pons junction')) ####################### # CREATE MAIN BODY MESH ####################### cylinderShape = CylinderShape.CYLINDER_SHAPE_FULL if not halfBrainStem else CylinderShape.CYLINDER_SHAPE_LOWER_HALF - # Body coordinates + # brainstem coordinates cylinderCentralPath = CylinderCentralPath(region, centralPath, elementsCountAlong) base = CylinderEnds(elementsCountAcrossMajor, elementsCountAcrossMinor, centre=[0.0, 0.0, 0.0], @@ -544,40 +526,43 @@ def generateBaseMesh(cls, region, options): cylinderCentralPath=cylinderCentralPath, useCrossDerivatives=False) + # workaround for old Zinc field wrapper bug: must create brainstem coordinates field before reading file brainstem_coordinates = findOrCreateFieldCoordinates(fm, name="brainstem coordinates") - # Brain coordinates + # generate brainstem coordinates field in temporary region tmp_region = region.createRegion() tmp_fm = tmp_region.getFieldmodule() - tmp_brainstem_coordinates = findOrCreateFieldCoordinates(tmp_fm, name="brainstem coordinates") - - cylinderCentralPath1 = CylinderCentralPath(tmp_region, brainstemPath, elementsCountAlong) - - base1 = CylinderEnds(elementsCountAcrossMajor, elementsCountAcrossMinor, - centre=[0.0, 0.0, 0.0], - alongAxis=cylinderCentralPath1.alongAxis[0], - majorAxis=cylinderCentralPath1.majorAxis[0], - minorRadius=cylinderCentralPath1.minorRadii[0]) - - cylinder2 = CylinderMesh(tmp_fm, tmp_brainstem_coordinates, elementsCountAlong, base1, - cylinderShape=cylinderShape, - cylinderCentralPath=cylinderCentralPath1, - useCrossDerivatives=False) - - # Write two coordinates - sir = tmp_region.createStreaminformationRegion() - srm = sir.createStreamresourceMemory() - tmp_region.write(sir) - result, buffer = srm.getBuffer() - - sir = region.createStreaminformationRegion() - srm = sir.createStreamresourceMemoryBuffer(buffer) - region.read(sir) - - del srm - del sir + with ChangeManager(tmp_fm): + tmp_brainstem_coordinates = findOrCreateFieldCoordinates(tmp_fm, name="brainstem coordinates") + + cylinderCentralPath1 = CylinderCentralPath(tmp_region, brainstemPath, elementsCountAlong) + + base1 = CylinderEnds(elementsCountAcrossMajor, elementsCountAcrossMinor, + centre=[0.0, 0.0, 0.0], + alongAxis=cylinderCentralPath1.alongAxis[0], + majorAxis=cylinderCentralPath1.majorAxis[0], + minorRadius=cylinderCentralPath1.minorRadii[0]) + + cylinder2 = CylinderMesh(tmp_fm, tmp_brainstem_coordinates, elementsCountAlong, base1, + cylinderShape=cylinderShape, + cylinderCentralPath=cylinderCentralPath1, + useCrossDerivatives=False) + + # write to memory buffer + sir = tmp_region.createStreaminformationRegion() + srm = sir.createStreamresourceMemory() + tmp_region.write(sir) + result, buffer = srm.getBuffer() + + # read into main region + sir = region.createStreaminformationRegion() + srm = sir.createStreamresourceMemoryBuffer(buffer) + region.read(sir) + + del srm + del sir + del tmp_brainstem_coordinates del tmp_fm - del tmp_brainstem_coordinates del tmp_region # Annotating groups @@ -596,49 +581,22 @@ def generateBaseMesh(cls, region, options): ################ # point markers ################ - pointMarkers = [ - {"group": dorsalMidCaudalGroup, "marker_brainstem_coordinates": [0.0, 1.0, 0.0]}, - {"group": ventralMidCaudalGroup, "marker_brainstem_coordinates": [0.0, -1.0, 0.0]}, - {"group": dorsalMidCranGroup, "marker_brainstem_coordinates": [0.0, 1.0, 8.0]}, - {"group": ventralMidCranGroup, "marker_brainstem_coordinates": [0.0, -1.0, 8.0]}, - {"group": dorsalMidMedullaPonsJunction, "marker_brainstem_coordinates": [0.0, 1.0, 3.0]}, - {"group": ventralMidMedullaPonsJunction, "marker_brainstem_coordinates": [0.0, -1.0, 3.0]}, - {"group": dorsalMidMidbrainPonsJunction, "marker_brainstem_coordinates": [0.0, 1.0, 6.0]}, - {"group": ventralMidMidbrainPonsJunction, "marker_brainstem_coordinates": [0.0, -1.0, 6.0]}, - ] - - markerGroup = findOrCreateFieldGroup(fm, "marker") - markerName = findOrCreateFieldStoredString(fm, name="marker_name") - markerLocation = findOrCreateFieldStoredMeshLocation(fm, mesh, name="marker_location") + markerTermNameBrainstemCoordinatesMap = { + 'brainstem dorsal midline caudal point': [0.0, 1.0, 0.0], + 'brainstem ventral midline caudal point': [0.0, -1.0, 0.0], + 'brainstem dorsal midline cranial point': [0.0, 1.0, 8.0], + 'brainstem ventral midline cranial point': [0.0, -1.0, 8.0], + 'brainstem dorsal midline pons-medulla junction': [0.0, 1.0, 3.0], + 'brainstem ventral midline pons-medulla junction': [0.0, -1.0, 3.0], + 'brainstem dorsal midline midbrain-pons junction': [0.0, 1.0, 6.0], + 'brainstem ventral midline midbrain-pons junction': [0.0, -1.0, 6.0] + } nodes = fm.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES) - markerPoints = findOrCreateFieldNodeGroup(markerGroup, nodes).getNodesetGroup() - markerBrainstemCoordinates = findOrCreateFieldCoordinates(fm, name="marker_body_coordinates") - markerTemplateInternal = nodes.createNodetemplate() - markerTemplateInternal.defineField(markerName) - markerTemplateInternal.defineField(markerLocation) - markerTemplateInternal.defineField(markerBrainstemCoordinates) - - cache = fm.createFieldcache() - - brainstemNodesetGroup = brainstemGroup.getNodesetGroup(nodes) - nodeIdentifier = max(1, getMaximumNodeIdentifier(nodes) + 1) - findMarkerLocation = fm.createFieldFindMeshLocation(markerBrainstemCoordinates, brainstem_coordinates, mesh) - findMarkerLocation.setSearchMode(FieldFindMeshLocation.SEARCH_MODE_EXACT) - - for pointMarker in pointMarkers: - group = pointMarker["group"] - markerPoint = markerPoints.createNode(nodeIdentifier, markerTemplateInternal) - cache.setNode(markerPoint) - - markerBrainstemCoordinates.assignReal(cache, pointMarker["marker_brainstem_coordinates"]) - markerName.assignString(cache, group.getName()) - - element, xi = findMarkerLocation.evaluateMeshLocation(cache, 3) - markerLocation.assignMeshLocation(cache, element, xi) - group.getNodesetGroup(nodes).addNode(markerPoint) - brainstemNodesetGroup.addNode(markerPoint) + for termName, brainstemCoordinatesValues in markerTermNameBrainstemCoordinatesMap.items(): + annotationGroup = findOrCreateAnnotationGroupForTerm(annotationGroups, region, get_brainstem_term(termName)) + annotationGroup.createMarkerNode(nodeIdentifier, brainstem_coordinates, brainstemCoordinatesValues) nodeIdentifier += 1 return annotationGroups @@ -701,7 +659,7 @@ def createCranialNerveEmergentMarkers(region, mesh, coordinatesName): # return element xi # use findMeshLocation to find the elementxi in an arbitrary mesh of given number of elements. - if coordinatesName == "bodyCoordinates": + if coordinatesName == "brainstem coordinates": # brainstem_coordinates: the left-side nerves nerveDict = {'OCULOMOTOR_left': [-0.13912257342955267, -0.5345161733750351, -0.7374762051676923], 'TROCHLEAR_left': [-0.13148279719950992, 0.4218745504359067, -0.7375838988856348], diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_heartventricles1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_heartventricles1.py index ad812ac8..0690edd3 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_heartventricles1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_heartventricles1.py @@ -255,16 +255,7 @@ def generateBaseMesh(cls, region, options): vSeptumGroup = AnnotationGroup(region, get_heart_term("interventricular septum")) annotationGroups = [ heartGroup, apexGroup, lvGroup, rvGroup, vSeptumGroup ] - # 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 @@ -1025,13 +1016,10 @@ def generateBaseMesh(cls, region, options): # apex annotation point element1 = mesh.findElementByIdentifier(1) - markerPoint = markerPoints.createNode(nodeIdentifier, markerTemplateInternal) - nodeIdentifier += 1 - cache.setNode(markerPoint) - markerName.assignString(cache, apexGroup.getName()) - markerLocation.assignMeshLocation(cache, element1, [ 0.0, 0.0, 1.0 ]) - for group in [ heartGroup, lvGroup, apexGroup ]: - group.getNodesetGroup(nodes).addNode(markerPoint) + markerNode = apexGroup.createMarkerNode(nodeIdentifier, element=element1, xi=[0.0, 0.0, 1.0]) + nodeIdentifier = markerNode.getIdentifier() + 1 + for group in [ heartGroup, lvGroup ]: + group.getNodesetGroup(nodes).addNode(markerNode) return annotationGroups diff --git a/src/scaffoldmaker/scaffoldpackage.py b/src/scaffoldmaker/scaffoldpackage.py index 918ed833..66f0267e 100644 --- a/src/scaffoldmaker/scaffoldpackage.py +++ b/src/scaffoldmaker/scaffoldpackage.py @@ -8,7 +8,9 @@ from opencmiss.maths.vectorops import euler_to_rotation_matrix from opencmiss.utils.zinc.field import createFieldEulerAnglesRotationMatrix +from opencmiss.utils.zinc.finiteelement import get_maximum_node_identifier from opencmiss.utils.zinc.general import ChangeManager +from opencmiss.zinc.field import Field from scaffoldmaker.annotation.annotationgroup import AnnotationGroup, findAnnotationGroupByName from scaffoldmaker.meshtypes.scaffold_base import Scaffold_base @@ -51,7 +53,7 @@ def __init__(self, scaffoldType, dct={}, defaultParameterSetName='Default'): meshEdits = copy.deepcopy(meshEdits) self._meshEdits = meshEdits self._autoAnnotationGroups = [] - # read user AnnotationGroups dict: + # read user AnnotationGroups list in dict form: userAnnotationGroupsDict = dct.get('userAnnotationGroups') # serialised form of user annotation groups, read from serialisation before generate(), updated before writing self._userAnnotationGroupsDict = copy.deepcopy(userAnnotationGroupsDict) if userAnnotationGroupsDict else [] @@ -59,6 +61,7 @@ def __init__(self, scaffoldType, dct={}, defaultParameterSetName='Default'): self._userAnnotationGroups = [] # region is set in generate(); can only instantiate user AnnotationGroups then self._region = None + self._nextNodeIdentifier = 1 def __eq__(self, other): ''' @@ -232,6 +235,9 @@ def generate(self, region, applyTransformation=True): self._region = region with ChangeManager(region.getFieldmodule()): self._autoAnnotationGroups = self._scaffoldType.generateMesh(region, self._scaffoldSettings) + # need next node identifier for creating user-defined marker points + nodes = region.getFieldmodule().findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES) + self._nextNodeIdentifier = get_maximum_node_identifier(nodes) + 1 if self._meshEdits: # apply mesh edits, a Zinc-readable model file containing node edits # Note: these are untransformed coordinates @@ -243,6 +249,12 @@ def generate(self, region, applyTransformation=True): if applyTransformation: self.applyTransformation() + def getNextNodeIdentifier(self): + """ + :return: First node identifier to try using after generating. Used for user-defined marker points. + """ + return self._nextNodeIdentifier + def getAnnotationGroups(self): ''' Empty until after call to generate(). @@ -261,14 +273,19 @@ def createUserAnnotationGroup(self, term=None): ''' Create a new, empty user annotation group. Only call after generate(). - :param term: Identifier for anatomical term, currently a tuple of name, id. - e.g. ('heart', 'FMA:7088'). Or None to generate a unique name. Name must be - unique if supplied; id should be unique but may be None. + :param term: Identifier for anatomical term, a tuple of (name, id), or None + to generate a unique name 'group#' with ID "None". + If a term is supplied, both its name and id must be strings, the name + must be unique i.e. unused in the list of annotation groups. + The id should be a unique string, or use "None" if unknown. + e.g. ('heart', 'FMA:7088') or None. :return: New AnnotationGroup. ''' assert self._region - if term: - assert not self.findAnnotationGroupByName(term[0]) + if term is not None: + assert isinstance(term, tuple) and (len(term) == 2) and all(isinstance(s, str) for s in term),\ + "Invalid annotation term " + str(term) + assert not self.findAnnotationGroupByName(term[0]), "Annotation term " + str(term) + " name is in use" useTerm = term else: number = 1 @@ -277,7 +294,7 @@ def createUserAnnotationGroup(self, term=None): if not self.findAnnotationGroupByName(name): break number += 1 - useTerm = (name, None) + useTerm = (name, "None") annotationGroup = AnnotationGroup(self._region, useTerm) self._userAnnotationGroups.append(annotationGroup) return annotationGroup @@ -288,6 +305,7 @@ def deleteAnnotationGroup(self, annotationGroup): :return: True on success, otherwise False ''' if annotationGroup and self.isUserAnnotationGroup(annotationGroup): + annotationGroup.clear() self._userAnnotationGroups.remove(annotationGroup) return True return False diff --git a/src/scaffoldmaker/utils/zinc_utils.py b/src/scaffoldmaker/utils/zinc_utils.py index 5413d7cd..95605ea6 100644 --- a/src/scaffoldmaker/utils/zinc_utils.py +++ b/src/scaffoldmaker/utils/zinc_utils.py @@ -5,10 +5,10 @@ from opencmiss.utils.zinc.field import findOrCreateFieldCoordinates from opencmiss.utils.zinc.general import ChangeManager from opencmiss.zinc.context import Context -from opencmiss.zinc.element import MeshGroup +from opencmiss.zinc.element import Mesh, MeshGroup from opencmiss.zinc.field import Field, FieldGroup from opencmiss.zinc.fieldmodule import Fieldmodule -from opencmiss.zinc.node import Node +from opencmiss.zinc.node import Node, Nodeset from opencmiss.zinc.result import RESULT_OK from scaffoldmaker.utils import interpolation as interp from scaffoldmaker.utils import vector @@ -146,6 +146,31 @@ def createFaceMeshGroupExteriorOnFace(fieldmodule : Fieldmodule, elementFaceType del isOnFace return faceMeshGroup + +def get_highest_dimension_mesh(fieldmodule : Fieldmodule) -> Mesh: + ''' + Get highest dimension non-empty mesh. + :return: Zinc Mesh or None if all are empty. + ''' + for dimension in range(3, 0, -1): + mesh = fieldmodule.findMeshByDimension(dimension) + if mesh.getSize() > 0: + return mesh + return None + + +def get_next_unused_node_identifier(nodeset: Nodeset, start_identifier=1) -> int: + """ + :return: Unused node identifier >= start_identifier. + """ + identifier = start_identifier + node = nodeset.findNodeByIdentifier(identifier) + while node.isValid(): + identifier += 1 + node = nodeset.findNodeByIdentifier(identifier) + return identifier + + def group_add_group_elements(group : FieldGroup, other_group : FieldGroup, only_dimension=None): ''' Add to group elements and/or nodes from other_group. @@ -169,6 +194,7 @@ def group_add_group_elements(group : FieldGroup, other_group : FieldGroup, only_ nodeset_group = node_group.getNodesetGroup() nodeset_group.addNodesConditional(other_group.getFieldNodeGroup(nodeset)) + def group_get_highest_dimension(group : FieldGroup): ''' Get highest dimension of elements or nodes in group. @@ -186,6 +212,7 @@ def group_get_highest_dimension(group : FieldGroup): return 0 return -1 + def identifier_ranges_fix(identifier_ranges): ''' Sort from lowest to highest identifier and merge adjacent and overlapping @@ -202,11 +229,13 @@ def identifier_ranges_fix(identifier_ranges): else: i += 1 + def identifier_ranges_from_string(identifier_ranges_string): ''' Parse string containing identifiers and identifier ranges. Function is suitable for processing manual input with whitespace, trailing non-digits. Ranges are sorted so strictly increasing. Overlapping ranges are merged. + Future: migrate to use .. as separator for compatibility with EX file groups and cmgui. :param identifier_ranges_string: Identifier ranges as a string e.g. '1-30,55,66-70'. '30-1, 55,66-70s' also produces the same result. :return: Ordered list of identifier ranges e.g. [[1,30],[55,55],[66,70]] @@ -242,6 +271,7 @@ def identifier_ranges_from_string(identifier_ranges_string): def identifier_ranges_to_string(identifier_ranges): ''' Convert ranges to a string, contracting single object ranges. + Future: migrate to use .. as separator for compatibility with EX file groups and cmgui. :param identifier_ranges: Ordered list of identifier ranges e.g. [[1,30],[55,55],[66,70]] :return: Identifier ranges as a string e.g. '1-30,55,66-70' ''' @@ -328,6 +358,7 @@ def nodeset_group_to_identifier_ranges(nodeset_group): ''' return domain_iterator_to_identifier_ranges(nodeset_group.createNodeiterator()) + def mesh_destroy_elements_and_nodes_by_identifiers(mesh, element_identifiers): ''' Deletes elements and related nodes using element identifiers. diff --git a/tests/test_general.py b/tests/test_general.py index 1743837d..4825bfff 100644 --- a/tests/test_general.py +++ b/tests/test_general.py @@ -8,6 +8,7 @@ from opencmiss.zinc.result import RESULT_OK from scaffoldmaker.annotation.annotationgroup import AnnotationGroup from scaffoldmaker.meshtypes.meshtype_3d_box1 import MeshType_3d_box1 +from scaffoldmaker.meshtypes.meshtype_3d_brainstem import MeshType_3d_brainstem1 from scaffoldmaker.meshtypes.meshtype_3d_heartatria1 import MeshType_3d_heartatria1 from scaffoldmaker.scaffoldpackage import ScaffoldPackage from scaffoldmaker.scaffolds import Scaffolds @@ -115,7 +116,7 @@ def test_user_annotation_groups(self): annotationGroup1 = scaffoldPackage.createUserAnnotationGroup() self.assertEqual('group1', annotationGroup1.getName()) # default name - self.assertIsNone(annotationGroup1.getId()) + self.assertEqual('None', annotationGroup1.getId()) self.assertTrue(scaffoldPackage.isUserAnnotationGroup(annotationGroup1)) self.assertEqual(-1, annotationGroup1.getDimension()) # -1 = empty group = annotationGroup1.getGroup() @@ -193,6 +194,163 @@ def test_user_annotation_groups(self): identifier_ranges_string = identifier_ranges_to_string(nodeset_group_to_identifier_ranges(nodesetGroup2)) self.assertEqual('1,3-5,7', identifier_ranges_string) + def test_user_marker_points(self): + """ + Test user marker point on brainstem1 scaffold which defined "brainstem coordinates". + """ + scaffoldPackage = ScaffoldPackage(MeshType_3d_brainstem1) + + context = Context("Test") + region = context.getDefaultRegion() + self.assertTrue(region.isValid()) + fieldmodule = region.getFieldmodule() + mesh = fieldmodule.findMeshByDimension(3) + nodes = fieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES) + + scaffoldPackage.generate(region) + # 1 higher than last node in scaffold. Make marker nodes from this number + nextNodeIdentifier = scaffoldPackage.getNextNodeIdentifier() + + brainstemCoordinatesField = fieldmodule.findFieldByName("brainstem coordinates") + self.assertTrue(brainstemCoordinatesField.isValid()) + + annotationGroups = scaffoldPackage.getAnnotationGroups() + self.assertEqual(18, len(annotationGroups)) + TOL = 1.0E-6 # coarse to handle find xi tolerances + + # check a built-in non-marker annotation group + ponsGroup = scaffoldPackage.findAnnotationGroupByName('pons') + self.assertFalse(ponsGroup.isMarker()) + self.assertRaises(AssertionError, lambda: ponsGroup.getMarkerMaterialCoordinates()) + self.assertRaises(AssertionError, lambda: ponsGroup.getMarkerLocation()) + self.assertRaises(AssertionError, lambda: ponsGroup.createMarkerNode(nextNodeIdentifier)) + + # check a built-in marker annotation group + brainstemVentralCranialPointGroup = \ + scaffoldPackage.findAnnotationGroupByName('brainstem ventral midline cranial point') + self.assertIsNotNone(brainstemVentralCranialPointGroup) + self.assertTrue(brainstemVentralCranialPointGroup.isMarker()) + brainstemCoordinatesFieldOut, brainstemCoordinatesValueOut = \ + brainstemVentralCranialPointGroup.getMarkerMaterialCoordinates() + self.assertEqual(brainstemCoordinatesFieldOut, brainstemCoordinatesField) + assertAlmostEqualList(self, [0.0, -1.0, 8.0], brainstemCoordinatesValueOut, delta=TOL) + elementOut, xiOut = brainstemVentralCranialPointGroup.getMarkerLocation() + self.assertEqual(235, elementOut.getIdentifier()) + assertAlmostEqualList(self, [1.0, 1.0, 0.0], xiOut, delta=TOL) + self.assertRaises(AssertionError, lambda: brainstemVentralCranialPointGroup.createMarkerNode(nextNodeIdentifier)) + + # check a non-existant annotation group + bobGroup = scaffoldPackage.findAnnotationGroupByName("bob") + self.assertIsNone(bobGroup) + + # now make a marker annotation named "bob" at the default location + bobGroup = scaffoldPackage.createUserAnnotationGroup(('bob', 'BOB:1')) + self.assertTrue(scaffoldPackage.isUserAnnotationGroup(bobGroup)) + self.assertFalse(bobGroup.isMarker()) + node = bobGroup.createMarkerNode(nextNodeIdentifier) + bobNodeIdentifier = node.getIdentifier() + self.assertEqual(nextNodeIdentifier, bobNodeIdentifier) + brainstemCoordinatesFieldOut, brainstemCoordinatesValueOut = bobGroup.getMarkerMaterialCoordinates() + self.assertIsNone(brainstemCoordinatesFieldOut) + self.assertIsNone(brainstemCoordinatesValueOut) + elementOut, xiOut = bobGroup.getMarkerLocation() + self.assertEqual(1, elementOut.getIdentifier()) + assertAlmostEqualList(self, [0.0, 0.0, 0.0], xiOut, delta=TOL) + # now request brainstem coordinates and let the annotation group determine its values from element:xi + bobGroup.setMarkerMaterialCoordinates(brainstemCoordinatesField) + brainstemCoordinatesFieldOut, brainstemCoordinatesValueOut = bobGroup.getMarkerMaterialCoordinates() + self.assertEqual(brainstemCoordinatesFieldOut, brainstemCoordinatesField) + # these should be precisely cos(45) but are not due to ellipse approximations + assertAlmostEqualList(self, [0.707016, -0.707198, 0], brainstemCoordinatesValueOut, delta=TOL) + # set element:xi location and check brainstem coordinates change + bobGroup.setMarkerLocation(mesh.findElementByIdentifier(33), [0.0, 0.5, 0.0]) + brainstemCoordinatesFieldOut, brainstemCoordinatesValueOut = bobGroup.getMarkerMaterialCoordinates() + self.assertEqual(brainstemCoordinatesFieldOut, brainstemCoordinatesField) + assertAlmostEqualList(self, [0.707016, -0.707198, 1.5], brainstemCoordinatesValueOut, delta=TOL) + # assign brainstem coordinates and check element:xi has moved + bobGroup.setMarkerMaterialCoordinates(brainstemCoordinatesField, [-0.1, -0.5, 2.2]) + brainstemCoordinatesFieldOut, brainstemCoordinatesValueOut = bobGroup.getMarkerMaterialCoordinates() + self.assertEqual(brainstemCoordinatesFieldOut, brainstemCoordinatesField) + assertAlmostEqualList(self, [-0.1, -0.5, 2.2], brainstemCoordinatesValueOut, delta=TOL) + elementOut, xiOut = bobGroup.getMarkerLocation() + self.assertEqual(82, elementOut.getIdentifier()) + assertAlmostEqualList(self, [0.305595, 0.2, 0.485941], xiOut, delta=TOL) + + # now make a marker annotation named "fred" with brainstem coordinates from the start + fredGroup = scaffoldPackage.createUserAnnotationGroup(('fred', 'FRED:1')) + # AnnotationGroup.createMarkerNode increments nextNodeIdentifier to one not used by existing node + node = fredGroup.createMarkerNode(nextNodeIdentifier, brainstemCoordinatesField, [0.5, 0.5, 4]) + fredNodeIdentifier = node.getIdentifier() + self.assertEqual(nextNodeIdentifier + 1, fredNodeIdentifier) + del node + brainstemCoordinatesFieldOut, brainstemCoordinatesValueOut = fredGroup.getMarkerMaterialCoordinates() + self.assertEqual(brainstemCoordinatesFieldOut, brainstemCoordinatesField) + assertAlmostEqualList(self, [0.5, 0.5, 4], brainstemCoordinatesValueOut, delta=TOL) + elementOut, xiOut = fredGroup.getMarkerLocation() + self.assertEqual(105, elementOut.getIdentifier()) + assertAlmostEqualList(self, [0.346095, 1, 0.66399], xiOut, delta=TOL) + + annotationGroups = scaffoldPackage.getAnnotationGroups() + self.assertEqual(20, len(annotationGroups)) + + # test deleting a marker annotation group + scaffoldPackage.deleteAnnotationGroup(fredGroup) + + annotationGroups = scaffoldPackage.getAnnotationGroups() + self.assertEqual(19, len(annotationGroups)) + + # check node fred has gone + node = nodes.findNodeByIdentifier(fredNodeIdentifier) + self.assertFalse(node.isValid()) + + # re-recreate fred with just element:xi location + fredGroup = scaffoldPackage.createUserAnnotationGroup(('fred', 'FRED:1')) + element = mesh.findElementByIdentifier(105) + node = fredGroup.createMarkerNode(nextNodeIdentifier, element=element, xi=[0.346095, 1, 0.66399]) + self.assertEqual(fredNodeIdentifier, node.getIdentifier()) + del node + elementOut, xiOut = fredGroup.getMarkerLocation() + self.assertEqual(105, elementOut.getIdentifier()) + assertAlmostEqualList(self, [0.346095, 1, 0.66399], xiOut, delta=TOL) + + # check the total number of groups + annotationGroups = scaffoldPackage.getAnnotationGroups() + self.assertEqual(20, len(annotationGroups)) + + # test serialisation + dct = scaffoldPackage.toDict() + scaffoldType = Scaffolds().findScaffoldTypeByName(dct['scaffoldTypeName']) + + scaffoldPackage2 = ScaffoldPackage(scaffoldType, dct) + region2 = context.createRegion() + fieldmodule2 = region2.getFieldmodule() + + scaffoldPackage2.generate(region2) + + brainstemCoordinatesField2 = fieldmodule2.findFieldByName("brainstem coordinates") + self.assertTrue(brainstemCoordinatesField2.isValid()) + + annotationGroups2 = scaffoldPackage2.getAnnotationGroups() + self.assertEqual(20, len(annotationGroups2)) + + # check user markers have been defined correctly for scaffoldPackage2 + + bobGroup2 = scaffoldPackage2.findAnnotationGroupByName('bob') + brainstemCoordinatesFieldOut, brainstemCoordinatesValueOut = bobGroup2.getMarkerMaterialCoordinates() + self.assertEqual(brainstemCoordinatesFieldOut, brainstemCoordinatesField2) + assertAlmostEqualList(self, [-0.1, -0.5, 2.2], brainstemCoordinatesValueOut, delta=TOL) + elementOut, xiOut = bobGroup2.getMarkerLocation() + self.assertEqual(82, elementOut.getIdentifier()) + assertAlmostEqualList(self, [0.305595, 0.2, 0.485941], xiOut, delta=TOL) + + fredGroup2 = scaffoldPackage2.findAnnotationGroupByName('fred') + brainstemCoordinatesFieldOut, brainstemCoordinatesValueOut = fredGroup2.getMarkerMaterialCoordinates() + self.assertIsNone(brainstemCoordinatesFieldOut) + self.assertIsNone(brainstemCoordinatesValueOut) + elementOut, xiOut = fredGroup2.getMarkerLocation() + self.assertEqual(105, elementOut.getIdentifier()) + assertAlmostEqualList(self, [0.346095, 1, 0.66399], xiOut, delta=TOL) + if __name__ == "__main__": unittest.main() From d3cc1ef9623290673cd80cceddbab8b52dd476a2 Mon Sep 17 00:00:00 2001 From: Richard Christie Date: Fri, 25 Mar 2022 12:47:13 +1300 Subject: [PATCH 2/3] Fix annotation marker points Fix marker points on 1D elements. Tidy handling of marker group and fields. --- .../annotation/annotationgroup.py | 72 +++++++++++++------ src/scaffoldmaker/scaffoldpackage.py | 10 ++- 2 files changed, 59 insertions(+), 23 deletions(-) diff --git a/src/scaffoldmaker/annotation/annotationgroup.py b/src/scaffoldmaker/annotation/annotationgroup.py index ce25f9fe..534f2524 100644 --- a/src/scaffoldmaker/annotation/annotationgroup.py +++ b/src/scaffoldmaker/annotation/annotationgroup.py @@ -6,8 +6,9 @@ from opencmiss.utils.zinc.general import ChangeManager from opencmiss.utils.zinc.field import find_or_create_field_coordinates, find_or_create_field_group, \ find_or_create_field_stored_mesh_location, find_or_create_field_stored_string -from opencmiss.zinc.element import Element -from opencmiss.zinc.field import Field, FieldFiniteElement, FieldGroup +from opencmiss.zinc.element import Element, Mesh +from opencmiss.zinc.field import Field, FieldFiniteElement, FieldGroup, FieldStoredMeshLocation, FieldStoredString +from opencmiss.zinc.fieldmodule import Fieldmodule from opencmiss.zinc.node import Node from opencmiss.zinc.result import RESULT_OK from scaffoldmaker.utils.zinc_utils import get_highest_dimension_mesh, get_next_unused_node_identifier, \ @@ -29,17 +30,11 @@ def __init__(self, region, term): ''' self._name = term[0] self._id = term[1] - fieldmodule = region.getFieldmodule() - field = fieldmodule.findFieldByName(self._name) - if field.isValid(): - self._group = field.castGroup() - assert self._group.isValid(), 'AnnotationGroup found existing non-group field called ' + self._name - else: - with ChangeManager(fieldmodule): - self._group = fieldmodule.createFieldGroup() - self._group.setName(self._name) - self._group.setManaged(True) + self._group = find_or_create_field_group(region.getFieldmodule(), self._name, managed=False) + assert self._group.getName() == self._name, \ + 'AnnotationGroup found existing non-group field called ' + self._name self._isMarker = False + self._markerGroup = None # held while marker point exists self._markerMaterialCoordinatesField = None def clear(self): @@ -55,6 +50,7 @@ def clear(self): markerNode = self.getNodesetGroup(nodes).createNodeiterator().next() nodes.destroyNode(markerNode) self._isMarker = False + self._markerGroup = None self._markerMaterialCoordinatesField = None self._group.clear() @@ -206,7 +202,8 @@ def createMarkerNode(self, startNodeIdentifier=1, Important: annotation group must currently be empty, and elements must exist. The marker node is added to the marker group in addition to this group. Raises an exception if marker creation cannot be achieved. - :param startNodeIdentifier: First unused node identifier >= this is uses for marker node. + :param startNodeIdentifier: First unused node identifier >= this may use for marker node. Incremented until + an unused node identifier is found for the marker node. :param materialCoordinatesField: Material coordinates field to define location of marker point in. Must be a finite element type field for which isTypeCoordinates() is True, with up to 3 components, and at least as many components as the highest mesh dimension. @@ -233,13 +230,13 @@ def createMarkerNode(self, startNodeIdentifier=1, nodes = fieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES) nodeIdentifier = get_next_unused_node_identifier(nodes, startNodeIdentifier) with ChangeManager(fieldmodule): - markerGroup = find_or_create_field_group(fieldmodule, "marker") + markerGroup = getAnnotationMarkerGroup(fieldmodule) markerNodeGroup = markerGroup.getFieldNodeGroup(nodes) if not markerNodeGroup.isValid(): markerNodeGroup = markerGroup.createFieldNodeGroup(nodes) markerNodesetGroup = markerNodeGroup.getNodesetGroup() - markerLocation = find_or_create_field_stored_mesh_location(fieldmodule, mesh, name="marker_location") - markerName = find_or_create_field_stored_string(fieldmodule, name="marker_name") + markerLocation = getAnnotationMarkerLocationField(fieldmodule, mesh) + markerName = getAnnotationMarkerNameField(fieldmodule) nodetemplate = nodes.createNodetemplate() if materialCoordinatesField: assert RESULT_OK == nodetemplate.defineField(materialCoordinatesField) @@ -249,6 +246,7 @@ def createMarkerNode(self, startNodeIdentifier=1, assert RESULT_OK == self.getNodesetGroup(nodes).addNode(markerNode) assert RESULT_OK == markerNodesetGroup.addNode(markerNode) self._isMarker = True + self._markerGroup = markerGroup # maintain reference so group is not destroyed # assign marker name to be same as this group's name. This needs to be maintained. See setName() fieldcache = fieldmodule.createFieldcache() fieldcache.setNode(markerNode) @@ -269,12 +267,14 @@ def getMarkerLocation(self): assert self._isMarker fieldmodule = self._group.getFieldmodule() mesh = get_highest_dimension_mesh(fieldmodule) - markerLocation = find_or_create_field_stored_mesh_location(fieldmodule, mesh, name="marker_location") + markerLocation = getAnnotationMarkerLocationField(fieldmodule, mesh) fieldcache = fieldmodule.createFieldcache() nodes = fieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES) markerNode = self.getNodesetGroup(nodes).createNodeiterator().next() fieldcache.setNode(markerNode) element, xi = markerLocation.evaluateMeshLocation(fieldcache, mesh.getDimension()) + if not isinstance(xi, list): + xi = [xi] # workaround for Zinc 1-D xi being a plain float return element, xi def setMarkerLocation(self, element: Element, xi: list): @@ -289,9 +289,9 @@ def setMarkerLocation(self, element: Element, xi: list): fieldmodule = self._group.getFieldmodule() mesh = get_highest_dimension_mesh(fieldmodule) assert mesh.containsElement(element), "Invalid element, not in highest dimension mesh" - assert len(xi) >= mesh.getDimension() + assert isinstance(xi, list) and (len(xi) >= mesh.getDimension()) with ChangeManager(fieldmodule): - markerLocation = find_or_create_field_stored_mesh_location(fieldmodule, mesh, name="marker_location") + markerLocation = getAnnotationMarkerLocationField(fieldmodule, mesh) fieldcache = fieldmodule.createFieldcache() nodes = fieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES) markerNode = self.getNodesetGroup(nodes).createNodeiterator().next() @@ -361,6 +361,8 @@ def setMarkerMaterialCoordinates(self, materialCoordinatesField, materialCoordin self.setMarkerMaterialCoordinates(oldMaterialCoordinatesField) print("AnnotationGroup setMarkerMaterialCoordinates. Field is not defined on mesh. Reverting.") return + if not isinstance(xi, list): + xi = [xi] # workaround for Zinc 1-D xi being a plain float del findMeshLocation del constCoordinates else: @@ -397,7 +399,7 @@ def setName(self, name): self._name = name if self._isMarker: # assign marker name to be same as this group's name. This needs to be maintained - markerName = find_or_create_field_stored_string(fieldmodule, name="marker_name") + markerName = getAnnotationMarkerNameField(fieldmodule) fieldcache = fieldmodule.createFieldcache() markerNode = self.getNodesetGroup(nodes).createNodeiterator().next() fieldcache.setNode(markerNode) @@ -565,3 +567,33 @@ def mergeAnnotationGroups(*annotationGroupsIn): if not findAnnotationGroupByName(annotationGroups, agroup._name): annotationGroups.append(agroup) return annotationGroups + +def getAnnotationMarkerGroup(fieldmodule: Fieldmodule) -> FieldGroup: + """ + Find or create the standard Zinc Group which marker points are created in. + Clients should use this method rather than finding by standard name "marker". + :param fieldmodule: Zinc Fieldmodule to find/create group in. + :return: FieldGroup + """ + return find_or_create_field_group(fieldmodule, "marker", managed=False) + +def getAnnotationMarkerLocationField(fieldmodule: Fieldmodule, mesh: Mesh) -> FieldStoredMeshLocation: + """ + Find or create the standard Zinc Field used to store marker mesh locations in the highest dimension mesh. + Clients should use this method rather than finding by standard name "marker_location". + :param fieldmodule: Zinc Fieldmodule to find/create field in. + :param mesh: The highest dimension mesh in region. + :return: FieldStoredMeshLocation + """ + return find_or_create_field_stored_mesh_location(fieldmodule, mesh, name="marker_location", managed=False) + +def getAnnotationMarkerNameField(fieldmodule: Fieldmodule) -> FieldStoredString: + """ + Find or create the standard Zinc Field used to store marker names. + Clients should use this method rather than finding by standard name "marker_name". + :param fieldmodule: Zinc Fieldmodule to find/create group in. + :return: FieldStoredString + """ + return find_or_create_field_stored_string(fieldmodule, name="marker_name", managed=False) + + diff --git a/src/scaffoldmaker/scaffoldpackage.py b/src/scaffoldmaker/scaffoldpackage.py index 66f0267e..0c680544 100644 --- a/src/scaffoldmaker/scaffoldpackage.py +++ b/src/scaffoldmaker/scaffoldpackage.py @@ -52,6 +52,8 @@ def __init__(self, scaffoldType, dct={}, defaultParameterSetName='Default'): else: meshEdits = copy.deepcopy(meshEdits) self._meshEdits = meshEdits + self._isGenerated = False # set to True when generate() is called + # annotation groups automatically created by scaffold script = set in generate() self._autoAnnotationGroups = [] # read user AnnotationGroups list in dict form: userAnnotationGroupsDict = dct.get('userAnnotationGroups') @@ -99,8 +101,7 @@ def toDict(self): } if self._meshEdits: dct['meshEdits'] = self._meshEdits - if self._userAnnotationGroups: - self.updateUserAnnotationGroups() + self.updateUserAnnotationGroups() if self._userAnnotationGroupsDict: dct['userAnnotationGroups'] = self._userAnnotationGroupsDict return dct @@ -108,8 +109,10 @@ def toDict(self): def updateUserAnnotationGroups(self): ''' Ensure user annotation groups are present in serialised form (dict). + Only done if scaffold has been generated. + :param force: If not True, ''' - if self._userAnnotationGroups: + if self._isGenerated: self._userAnnotationGroupsDict = [ annotationGroup.toDict() for annotationGroup in self._userAnnotationGroups ] def getMeshEdits(self): @@ -246,6 +249,7 @@ def generate(self, region, applyTransformation=True): region.read(sir) # define user AnnotationGroups from serialised Dict self._userAnnotationGroups = [ AnnotationGroup.fromDict(dct, self._region) for dct in self._userAnnotationGroupsDict ] + self._isGenerated = True if applyTransformation: self.applyTransformation() From eaab9988144aea533cbe69bbfb72d9216dce118c Mon Sep 17 00:00:00 2001 From: Richard Christie Date: Fri, 25 Mar 2022 20:38:09 +1300 Subject: [PATCH 3/3] Fix heart test Add fibrous ring surface groups to atria and ventricles with base. Catch exceptions in vtkExport as unsupported meshes cause it to crash. --- .../meshtypes/meshtype_3d_heartatria1.py | 39 +++++++++++- .../meshtype_3d_heartventriclesbase1.py | 62 ++++++++++++++++++- src/scaffoldmaker/utils/exportvtk.py | 19 +++--- tests/test_general.py | 8 +-- tests/test_heart.py | 4 +- 5 files changed, 113 insertions(+), 19 deletions(-) diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_heartatria1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_heartatria1.py index 5cfbc0ce..baf638d7 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_heartatria1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_heartatria1.py @@ -12,7 +12,7 @@ findOrCreateFieldNodeGroup, findOrCreateFieldStoredMeshLocation, findOrCreateFieldStoredString from opencmiss.utils.zinc.finiteelement import getMaximumElementIdentifier, getMaximumNodeIdentifier from opencmiss.zinc.element import Element, Elementbasis -from opencmiss.zinc.field import Field +from opencmiss.zinc.field import Field, FieldGroup from opencmiss.zinc.node import Node from scaffoldmaker.annotation.annotationgroup import AnnotationGroup, findOrCreateAnnotationGroupForTerm, getAnnotationGroupForTerm from scaffoldmaker.annotation.heart_terms import get_heart_term @@ -837,10 +837,11 @@ def generateBaseMesh(cls, region, options): svcGroup = AnnotationGroup(region, get_heart_term("superior vena cava")) laeVenousMidpointGroup = AnnotationGroup(region, get_heart_term("left atrium epicardium venous midpoint")) raeVenousMidpointGroup = AnnotationGroup(region, get_heart_term("right atrium epicardium venous midpoint")) - annotationGroups += [ laeVenousMidpointGroup, ivcGroup, ivcInletGroup, raeVenousMidpointGroup, svcGroup, svcInletGroup ] - # av boundary nodes are put in left and right fibrous ring groups only so they can be found by heart1 + # av boundary nodes are put in left and right fibrous ring groups 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 += [laeVenousMidpointGroup, ivcGroup, ivcInletGroup, raeVenousMidpointGroup, + svcGroup, svcInletGroup, lFibrousRingGroup, rFibrousRingGroup] # annotation fiducial points markerGroup = findOrCreateFieldGroup(fm, "marker") @@ -3312,6 +3313,38 @@ def defineFaceAnnotations(cls, region, options, annotationGroups): raaEpiGroup = findOrCreateAnnotationGroupForTerm(annotationGroups, region, get_heart_term("epicardium of right auricle")) raaEpiGroup.getMeshGroup(mesh2d).addElementsConditional(is_raa_epi) + 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_la, is_fibrous_ring) + lFibrousRingGroup.getMeshGroup(mesh2d).addElementsConditional(is_left_fibrous_ring) + is_right_fibrous_ring = fm.createFieldAnd(is_ra, is_fibrous_ring) + rFibrousRingGroup.getMeshGroup(mesh2d).addElementsConditional(is_right_fibrous_ring) def getLeftAtriumPulmonaryVeinOstiaElementsCounts(elementsCountAroundLeftAtriumFreeWall, elementsCountOverAtria, commonLeftRightPvOstium): ''' diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_heartventriclesbase1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_heartventriclesbase1.py index 1defc75d..c17c8d15 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_heartventriclesbase1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_heartventriclesbase1.py @@ -13,7 +13,7 @@ findOrCreateFieldNodeGroup, findOrCreateFieldStoredMeshLocation, findOrCreateFieldStoredString from opencmiss.utils.zinc.finiteelement import getMaximumElementIdentifier, getMaximumNodeIdentifier from opencmiss.zinc.element import Element -from opencmiss.zinc.field import Field +from opencmiss.zinc.field import Field, FieldGroup from opencmiss.zinc.node import Node from scaffoldmaker.annotation.annotationgroup import AnnotationGroup, getAnnotationGroupForTerm from scaffoldmaker.annotation.heart_terms import get_heart_term @@ -307,10 +307,10 @@ def generateBaseMesh(cls, region, options): 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")) + annotationGroups += [conusArteriosusGroup, lFibrousRingGroup, rFibrousRingGroup] # annotation fiducial points markerGroup = findOrCreateFieldGroup(fm, "marker") @@ -1405,3 +1405,61 @@ 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() + 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() + 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) diff --git a/src/scaffoldmaker/utils/exportvtk.py b/src/scaffoldmaker/utils/exportvtk.py index 3203a6dd..3ec632d9 100644 --- a/src/scaffoldmaker/utils/exportvtk.py +++ b/src/scaffoldmaker/utils/exportvtk.py @@ -4,6 +4,7 @@ import io import os +import sys from sys import version_info from opencmiss.utils.zinc.finiteelement import getElementNodeIdentifiersBasisOrder @@ -42,7 +43,6 @@ def __init__(self, region, description, annotationGroups = None): if markerNodeGroup.isValid(): self._markerNodes = markerNodeGroup.getNodesetGroup() - def _write(self, outstream): if version_info.major > 2: assert isinstance(outstream, io.TextIOBase), 'ExportVtk.write: Invalid outstream argument' @@ -168,11 +168,14 @@ def _writeMarkers(self, outstream): def writeFile(self, filename): ''' - Export to legacy vtk file/ + 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) + try: + 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) + except Exception as e: + print("Failed to write VTK file", filename, file=sys.stderr); diff --git a/tests/test_general.py b/tests/test_general.py index 4825bfff..17f0e57a 100644 --- a/tests/test_general.py +++ b/tests/test_general.py @@ -107,7 +107,7 @@ def test_user_annotation_groups(self): scaffoldPackage.generate(region) annotationGroups = scaffoldPackage.getAnnotationGroups() - self.assertEqual(22, len(annotationGroups)) + self.assertEqual(24, len(annotationGroups)) endocardium_of_la = scaffoldPackage.findAnnotationGroupByName('endocardium of left atrium') self.assertTrue(isinstance(endocardium_of_la, AnnotationGroup)) @@ -146,7 +146,7 @@ def test_user_annotation_groups(self): self.assertEqual('group2', annotationGroup3.getName()) # default name self.assertTrue(scaffoldPackage.isUserAnnotationGroup(annotationGroup3)) annotationGroups = scaffoldPackage.getAnnotationGroups() - self.assertEqual(25, len(annotationGroups)) + self.assertEqual(27, len(annotationGroups)) # rename group1 to fred self.assertTrue(annotationGroup1.setName('fred')) @@ -156,7 +156,7 @@ def test_user_annotation_groups(self): self.assertTrue(scaffoldPackage.deleteAnnotationGroup(annotationGroup3)) annotationGroups = scaffoldPackage.getAnnotationGroups() - self.assertEqual(24, len(annotationGroups)) + self.assertEqual(26, len(annotationGroups)) # test serialisation dct = scaffoldPackage.toDict() @@ -170,7 +170,7 @@ def test_user_annotation_groups(self): scaffoldPackage2.generate(region2) annotationGroups2 = scaffoldPackage2.getAnnotationGroups() - self.assertEqual(24, len(annotationGroups2)) + self.assertEqual(26, len(annotationGroups2)) annotationGroup1 = scaffoldPackage2.findAnnotationGroupByName('fred') self.assertEqual('fred', annotationGroup1.getName()) diff --git a/tests/test_heart.py b/tests/test_heart.py index 52c1f30a..742c158b 100644 --- a/tests/test_heart.py +++ b/tests/test_heart.py @@ -72,10 +72,10 @@ def test_heart1(self): fieldcache = fieldmodule.createFieldcache() result, surfaceArea = surfaceAreaField.evaluateReal(fieldcache, 1) self.assertEqual(result, RESULT_OK) - self.assertAlmostEqual(surfaceArea, 36508.985114306226, delta=1.0E-2) + self.assertAlmostEqual(surfaceArea, 36497.828754581264, delta=1.0E-2) result, volume = volumeField.evaluateReal(fieldcache, 1) self.assertEqual(result, RESULT_OK) - self.assertAlmostEqual(volume, 221285.56778831664, delta=1.0E-2) + self.assertAlmostEqual(volume, 221275.23964745103, delta=1.0E-2) # check some annotationGroups: expectedSizes3d = {