diff --git a/src/scaffoldmaker/annotation/annotationgroup.py b/src/scaffoldmaker/annotation/annotationgroup.py index 1fa8ab75..534f2524 100644 --- a/src/scaffoldmaker/annotation/annotationgroup.py +++ b/src/scaffoldmaker/annotation/annotationgroup.py @@ -1,12 +1,18 @@ """ 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, 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 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 @@ -24,16 +30,29 @@ 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): + """ + 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._markerGroup = None + self._markerMaterialCoordinatesField = None + self._group.clear() def toDict(self): ''' @@ -62,8 +81,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 +115,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 +125,251 @@ 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 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. + 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 = getAnnotationMarkerGroup(fieldmodule) + markerNodeGroup = markerGroup.getFieldNodeGroup(nodes) + if not markerNodeGroup.isValid(): + markerNodeGroup = markerGroup.createFieldNodeGroup(nodes) + markerNodesetGroup = markerNodeGroup.getNodesetGroup() + markerLocation = getAnnotationMarkerLocationField(fieldmodule, mesh) + markerName = getAnnotationMarkerNameField(fieldmodule) + 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 + 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) + 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 = 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): + """ + 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 isinstance(xi, list) and (len(xi) >= mesh.getDimension()) + with ChangeManager(fieldmodule): + markerLocation = getAnnotationMarkerLocationField(fieldmodule, mesh) + 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 + if not isinstance(xi, list): + xi = [xi] # workaround for Zinc 1-D xi being a plain float + 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 +380,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 +397,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 = getAnnotationMarkerNameField(fieldmodule) + fieldcache = fieldmodule.createFieldcache() + markerNode = self.getNodesetGroup(nodes).createNodeiterator().next() + fieldcache.setNode(markerNode) + markerName.assignString(fieldcache, self._name) return True return False @@ -239,7 +529,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 +547,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.") @@ -275,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/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_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_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/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/scaffoldpackage.py b/src/scaffoldmaker/scaffoldpackage.py index 918ed833..0c680544 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 @@ -50,8 +52,10 @@ 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 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 +63,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): ''' @@ -96,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 @@ -105,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): @@ -232,6 +238,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 @@ -240,9 +249,16 @@ 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() + 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 +277,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 +298,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 +309,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/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/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..17f0e57a 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 @@ -106,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)) @@ -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() @@ -145,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')) @@ -155,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() @@ -169,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()) @@ -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() 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 = {