Transform raw data by moving nodes and renumbering data

Update to use opencmiss.utils.

setup.py

@@ -18,8 +18,8 @@ def readfile(filename, split=False):
 # into the 'requirements.txt' file.
 requires = [
     # minimal requirements listing
-    "opencmiss.utils @ https://api.github.com/repos/OpenCMISS-Bindings/opencmiss.utils/tarball/master"
-    #"opencmiss.zinc"  # not yet on pypi - need manual install
+    "opencmiss.utils >= 0.2",
+    "opencmiss.zinc"  # not yet on pypi - need manual install from opencmiss.org
 ]
 source_license = readfile("LICENSE")
 

src/scaffoldfitter/fitter.py

@@ -4,9 +4,9 @@
 
 import json
 from opencmiss.utils.zinc.field import assignFieldParameters, createFieldFiniteElementClone, getGroupList, getManagedFieldNames, \
-    getOrCreateFieldFiniteElement, getOrCreateFieldStoredMeshLocation, getUniqueFieldName, orphanFieldOfName
-from opencmiss.utils.zinc.finiteelement import evaluateNodesetMeanCoordinates, findNodeWithName
-from opencmiss.utils.zinc.general import ZincCacheChanges
+    findOrCreateFieldFiniteElement, findOrCreateFieldStoredMeshLocation, getUniqueFieldName, orphanFieldByName
+from opencmiss.utils.zinc.finiteelement import evaluateFieldNodesetMean, findNodeWithName, getMaximumNodeIdentifier
+from opencmiss.utils.zinc.general import ChangeManager
 from opencmiss.zinc.context import Context
 from opencmiss.zinc.field import Field, FieldFindMeshLocation, FieldGroup
 from opencmiss.zinc.result import RESULT_OK, RESULT_WARNING_PART_DONE
@@ -23,6 +23,7 @@ def __init__(self, zincModelFileName : str, zincDataFileName : str):
         self._zincDataFileName = zincDataFileName
         self._context = Context("Scaffoldfitter")
         self._region = None
+        self._rawDataRegion = None
         self._fieldmodule = None
         self._modelCoordinatesField = None
         self._modelCoordinatesFieldName = None
@@ -89,16 +90,72 @@ def load(self):
         """
         self._region = self._context.createRegion()
         self._fieldmodule = self._region.getFieldmodule()
+        self._rawDataRegion = self._region.createChild("raw_data")
+        self._loadModel()
+        self._loadData()
+        self._defineDataProjectionFields()
+        self.calculateDataProjections()
+
+    def _loadModel(self):
         result = self._region.readFile(self._zincModelFileName)
         assert result == RESULT_OK, "Failed to load model file" + str(self._zincModelFileName)
-        result = self._region.readFile(self._zincDataFileName)
-        assert result == RESULT_OK, "Failed to load data file" + str(self._zincDataFileName)
         self._mesh = [ self._fieldmodule.findMeshByDimension(d + 1) for d in range(3) ]
         self._discoverModelCoordinatesField()
+
+    def _loadData(self):
+        """
+        Load zinc data file into self._rawDataRegion.
+        Transfer data points (and converted nodes) into self._region.
+        """
+        result = self._rawDataRegion.readFile(self._zincDataFileName)
+        assert result == RESULT_OK, "Failed to load data file " + str(self._zincDataFileName)
+        # if there both nodes and datapoints, offset datapoint identifiers to ensure different
+        fieldmodule = self._rawDataRegion.getFieldmodule()
+        nodes = fieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES)
+        if nodes.getSize() > 0:
+            datapoints = fieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_DATAPOINTS)
+            if datapoints.getSize() > 0:
+                maximumDatapointIdentifier = max(0, getMaximumNodeIdentifier(datapoints))
+                maximumNodeIdentifier = max(0, getMaximumNodeIdentifier(nodes))
+                # this assumes identifiers are in low ranges and can be improved if there is a problem:
+                identifierOffset = 100000
+                while (maximumDatapointIdentifier > identifierOffset) or (maximumNodeIdentifier > identifierOffset):
+                    assert identifierOffset < 1000000000, "Invalid node and datapoint identifier ranges"
+                    identifierOffset *= 10
+                with ChangeManager(fieldmodule):
+                    while True:
+                        # logic relies on datapoints being in identifier order
+                        datapoint = datapoints.createNodeiterator().next()
+                        identifier = datapoint.getIdentifier()
+                        if identifier >= identifierOffset:
+                            break;
+                        result = datapoint.setIdentifier(identifier + identifierOffset)
+                        assert result == RESULT_OK, "Failed to offset datapoint identifier"
+            # transfer nodes as datapoints to self._region
+            sir = self._rawDataRegion.createStreaminformationRegion()
+            srm = sir.createStreamresourceMemory()
+            sir.setResourceDomainTypes(srm, Field.DOMAIN_TYPE_NODES)
+            self._rawDataRegion.write(sir)
+            result, buffer = srm.getBuffer()
+            assert result == RESULT_OK, "Failed to write nodes"
+            buffer = buffer.replace(bytes("!#nodeset nodes", "utf-8"), bytes("!#nodeset datapoints", "utf-8"))
+            sir = self._region.createStreaminformationRegion()
+            srm = sir.createStreamresourceMemoryBuffer(buffer)
+            result = self._region.read(sir)
+            assert result == RESULT_OK, "Failed to load nodes as datapoints"
+        # transfer datapoints to self._region
+        sir = self._rawDataRegion.createStreaminformationRegion()
+        srm = sir.createStreamresourceMemory()
+        sir.setResourceDomainTypes(srm, Field.DOMAIN_TYPE_DATAPOINTS)
+        self._rawDataRegion.write(sir)
+        result, buffer = srm.getBuffer()
+        assert result == RESULT_OK, "Failed to write datapoints"
+        sir = self._region.createStreaminformationRegion()
+        srm = sir.createStreamresourceMemoryBuffer(buffer)
+        result = self._region.read(sir)
+        assert result == RESULT_OK, "Failed to load datapoints"
         self._discoverDataCoordinatesField()
         self._discoverMarkerGroup()
-        self._defineDataProjectionFields()
-        self.calculateDataProjections()
 
     def getDataCoordinatesField(self):
         return self._dataCoordinatesField
@@ -248,8 +305,8 @@ def _calculateMarkerDataLocations(self):
         datapoints = self._fieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_DATAPOINTS)
         meshDimension = mesh.getDimension()
         fieldcache = self._fieldmodule.createFieldcache()
-        with ZincCacheChanges(self._fieldmodule):
-            self._markerDataLocationField = getOrCreateFieldStoredMeshLocation(self._fieldmodule, mesh, name=markerPrefix + "_data_location_" + mesh.getName())
+        with ChangeManager(self._fieldmodule):
+            self._markerDataLocationField = findOrCreateFieldStoredMeshLocation(self._fieldmodule, mesh, name=markerPrefix + "_data_location_" + mesh.getName(), managed=False)
             self._markerDataLocationGroupField = self._fieldmodule.createFieldNodeGroup(datapoints)
             self._markerDataLocationGroupField.setName(getUniqueFieldName(self._fieldmodule, markerPrefix + "_data_location_group"))
             self._markerDataLocationGroup = self._markerDataLocationGroupField.getNodesetGroup()
@@ -296,7 +353,7 @@ def _discoverMarkerGroup(self):
 
     def _updateMarkerCoordinatesField(self):
         if self._modelCoordinatesField and self._markerLocationField:
-            with ZincCacheChanges(self._fieldmodule):
+            with ChangeManager(self._fieldmodule):
                 markerPrefix = self._markerGroup.getName()
                 self._markerCoordinatesField = self._fieldmodule.createFieldEmbedded(self._modelCoordinatesField, self._markerLocationField)
                 self._markerCoordinatesField.setName(getUniqueFieldName(self._fieldmodule, markerPrefix + "_coordinates"))
@@ -305,7 +362,7 @@ def _updateMarkerCoordinatesField(self):
 
     def _updateMarkerDataLocationCoordinatesField(self):
         if self._modelCoordinatesField and self._markerDataLocationField:
-            with ZincCacheChanges(self._fieldmodule):
+            with ChangeManager(self._fieldmodule):
                 markerPrefix = self._markerGroup.getName()
                 self._markerDataLocationCoordinatesField = self._fieldmodule.createFieldEmbedded(self._modelCoordinatesField, self._markerDataLocationField)
                 self._markerDataLocationCoordinatesField.setName(getUniqueFieldName(self._fieldmodule, markerPrefix + "_data_location_coordinates"))
@@ -326,7 +383,7 @@ def setModelCoordinatesField(self, modelCoordinatesField : Field):
         self._modelCoordinatesField = finiteElementField
         self._modelCoordinatesFieldName = modelCoordinatesField.getName()
         modelReferenceCoordinatesFieldName = "reference_" + self._modelCoordinatesField.getName();
-        orphanFieldOfName(self._fieldmodule, modelReferenceCoordinatesFieldName);
+        orphanFieldByName(self._fieldmodule, modelReferenceCoordinatesFieldName);
         self._modelReferenceCoordinatesField = createFieldFiniteElementClone(self._modelCoordinatesField, modelReferenceCoordinatesFieldName)
         self._updateMarkerCoordinatesField()
         self._updateMarkerDataLocationCoordinatesField()
@@ -374,11 +431,11 @@ def _defineDataProjectionFields(self):
         self._dataProjectionErrorFields = []
         self._dataProjectionNodeGroupFields = []
         self._dataProjectionNodesetGroups = []
-        with ZincCacheChanges(self._fieldmodule):
+        with ChangeManager(self._fieldmodule):
             datapoints = self._fieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_DATAPOINTS)
             for d in range(2):
                 mesh = self._mesh[d]
-                dataProjectionLocation = getOrCreateFieldStoredMeshLocation(self._fieldmodule, mesh, name="data_projection_location_" + mesh.getName())
+                dataProjectionLocation = findOrCreateFieldStoredMeshLocation(self._fieldmodule, mesh, name="data_projection_location_" + mesh.getName(), managed=False)
                 self._dataProjectionLocationFields.append(dataProjectionLocation)
                 dataProjectionCoordinates = self._fieldmodule.createFieldEmbedded(self._modelCoordinatesField, dataProjectionLocation)
                 dataProjectionCoordinates.setName(getUniqueFieldName(self._fieldmodule, "data_projection_coordinates_" + mesh.getName()))
@@ -393,8 +450,8 @@ def _defineDataProjectionFields(self):
                 field.setName(getUniqueFieldName(self._fieldmodule, "data_projection_group_" + mesh.getName()))
                 self._dataProjectionNodeGroupFields.append(field)
                 self._dataProjectionNodesetGroups.append(field.getNodesetGroup())
-            self._dataProjectionDirectionField = getOrCreateFieldFiniteElement(self._fieldmodule, "data_projection_direction",
-                componentsCount = 3, componentNames = [ "x", "y", "z" ])
+            self._dataProjectionDirectionField = findOrCreateFieldFiniteElement(self._fieldmodule, "data_projection_direction",
+                components_count = 3, component_names = [ "x", "y", "z" ])
 
     def calculateDataProjections(self):
         """
@@ -403,7 +460,7 @@ def calculateDataProjections(self):
         Calculate and store projection direction unit vector.
         """
         assert self._dataCoordinatesField and self._modelCoordinatesField
-        with ZincCacheChanges(self._fieldmodule):
+        with ChangeManager(self._fieldmodule):
             findMeshLocation = None
             datapoints = self._fieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_DATAPOINTS)
             fieldcache = self._fieldmodule.createFieldcache()
@@ -485,7 +542,7 @@ def calculateDataProjections(self):
                     print("Warning: " + str(unprojected) + " data points with data coordinates have not been projected")
                 del unprojectedDatapoints
 
-            # remove temporary objects before clean up ZincCacheChanges
+            # remove temporary objects before ChangeManager exits
             del findMeshLocation
             del fieldcache
 
@@ -588,7 +645,7 @@ def evaluateNodeGroupMeanCoordinates(self, groupName, coordinatesFieldName, isDa
         nodesetGroup = group.getFieldNodeGroup(nodeset).getNodesetGroup()
         assert nodesetGroup.isValid()
         coordinates = self._fieldmodule.findFieldByName(coordinatesFieldName)
-        return evaluateNodesetMeanCoordinates(coordinates, nodesetGroup)
+        return evaluateFieldNodesetMean(coordinates, nodesetGroup)
 
     def getDiagnosticLevel(self):
         return self._diagnosticLevel

src/scaffoldfitter/fitterstepalign.py

@@ -2,9 +2,9 @@
 Fit step for gross alignment and scale.
 """
 
-from opencmiss.utils.zinc.field import assignFieldParameters, createTransformationFields
+from opencmiss.utils.zinc.field import assignFieldParameters, createFieldsTransformations
 from opencmiss.utils.zinc.finiteelement import getNodeNameCentres
-from opencmiss.utils.zinc.general import ZincCacheChanges
+from opencmiss.utils.zinc.general import ChangeManager
 from opencmiss.zinc.field import Field
 from opencmiss.zinc.optimisation import Optimisation
 from opencmiss.zinc.result import RESULT_OK, RESULT_WARNING_PART_DONE
@@ -103,9 +103,9 @@ def run(self):
         if self._alignMarkers:
             self._doAlignMarkers()
         fieldmodule = self._fitter._fieldmodule
-        with ZincCacheChanges(fieldmodule):
+        with ChangeManager(fieldmodule):
             # rotate, scale and translate model
-            modelCoordinatesTransformed = createTransformationFields(
+            modelCoordinatesTransformed = createFieldsTransformations(
                 modelCoordinates, self._rotation, self._scale, self._translation)[0]
             fieldassignment = self._fitter._modelCoordinatesField.createFieldassignment(modelCoordinatesTransformed)
             result = fieldassignment.assign()
@@ -165,7 +165,7 @@ def _optimiseAlignment(self, markerMap):
         assert len(markerMap) >= 3, "Align:  Only " + str(len(markerMap)) + " markers - need at least 3"
         region = self._fitter._context.createRegion()
         fieldmodule = region.getFieldmodule()
-        with ZincCacheChanges(fieldmodule):
+        with ChangeManager(fieldmodule):
             modelCoordinates = fieldmodule.createFieldFiniteElement(3)
             dataCoordinates = fieldmodule.createFieldFiniteElement(3)
             nodes = fieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES)
@@ -182,7 +182,7 @@ def _optimiseAlignment(self, markerMap):
                 result2 = dataCoordinates.assignReal(fieldcache, positions[1])
                 assert (result1 == RESULT_OK) and (result2 == RESULT_OK), "Align:  Failed to set up data for alignment to markers optimisation"
             del fieldcache
-            modelCoordinatesTransformed, rotation, scale, translation = createTransformationFields(modelCoordinates)
+            modelCoordinatesTransformed, rotation, scale, translation = createFieldsTransformations(modelCoordinates)
             # create objective = sum of squares of vector from modelCoordinatesTransformed to dataCoordinates
             markerDiff = fieldmodule.createFieldSubtract(dataCoordinates, modelCoordinatesTransformed)
             objective = fieldmodule.createFieldNodesetSumSquares(markerDiff, nodes)

src/scaffoldfitter/fitterstepfit.py

@@ -2,8 +2,8 @@
 Fit step for gross alignment and scale.
 """
 
-from opencmiss.utils.zinc.field import assignFieldParameters, createDisplacementGradientFields
-from opencmiss.utils.zinc.general import ZincCacheChanges
+from opencmiss.utils.zinc.field import assignFieldParameters, createFieldsDisplacementGradients
+from opencmiss.utils.zinc.general import ChangeManager
 from opencmiss.zinc.field import Field, FieldFindMeshLocation
 from opencmiss.zinc.optimisation import Optimisation
 from opencmiss.zinc.result import RESULT_OK
@@ -166,7 +166,7 @@ def createDataProjectionObjectiveField(self, dimension):
         :return: Zinc FieldNodesetSumSquares.
         """
         fieldmodule = self._fitter.getFieldmodule()
-        with ZincCacheChanges(fieldmodule):
+        with ChangeManager(fieldmodule):
             dataProjectionDelta = self._fitter.getDataProjectionDeltaField(dimension)
             #dataProjectionInDirection = fieldmodule.createFieldDotProduct(dataProjectionDelta, self._fitter.getDataProjectionDirectionField())
             #dataProjectionInDirection = fieldmodule.createFieldMagnitude(dataProjectionDelta)
@@ -182,7 +182,7 @@ def createMarkerObjectiveField(self, weight):
         :return: Zinc FieldNodesetSumSquares.
         """
         fieldmodule = self._fitter.getFieldmodule()
-        with ZincCacheChanges(fieldmodule):
+        with ChangeManager(fieldmodule):
             markerDataLocation, markerDataLocationCoordinates, markerDataDelta = self._fitter.getMarkerDataLocationFields()
             markerDataWeightedDelta = markerDataDelta*fieldmodule.createFieldConstant([ weight ]*markerDataDelta.getNumberOfComponents())
             markerDataObjective = fieldmodule.createFieldNodesetSumSquares(markerDataWeightedDelta, self._fitter.getMarkerDataLocationNodesetGroup())
@@ -196,8 +196,8 @@ def createDeformationPenaltyObjectiveField(self):
         numberOfGaussPoints = 3
         fieldmodule = self._fitter.getFieldmodule()
         mesh = self._fitter.getHighestDimensionMesh()
-        with ZincCacheChanges(fieldmodule):
-            displacementGradient1, displacementGradient2 = createDisplacementGradientFields(self._fitter.getModelCoordinatesField(), self._fitter.getModelReferenceCoordinatesField(), mesh)
+        with ChangeManager(fieldmodule):
+            displacementGradient1, displacementGradient2 = createFieldsDisplacementGradients(self._fitter.getModelCoordinatesField(), self._fitter.getModelReferenceCoordinatesField(), mesh)
             if self._strainPenaltyWeight > 0.0:
                 weightedDisplacementGradient1 = displacementGradient1*fieldmodule.createFieldConstant([ self._strainPenaltyWeight ]*displacementGradient1.getNumberOfComponents())
             else:
@@ -229,7 +229,7 @@ def createEdgeDiscontinuityPenaltyObjectiveField(self):
         numberOfGaussPoints = 3
         fieldmodule = self._fitter.getFieldmodule()
         lineMesh = fieldmodule.findMeshByDimension(1)
-        with ZincCacheChanges(fieldmodule):
+        with ChangeManager(fieldmodule):
             edgeDiscontinuity = fieldmodule.createFieldEdgeDiscontinuity(self._fitter.getModelCoordinatesField())
             weightedEdgeDiscontinuity = edgeDiscontinuity*fieldmodule.createFieldConstant(self._edgeDiscontinuityPenaltyWeight)
             edgeDiscontinuityPenaltyObjective = fieldmodule.createFieldMeshIntegralSquares(weightedEdgeDiscontinuity, self._fitter.getModelReferenceCoordinatesField(), lineMesh)

tests/resources/cube_to_sphere_data_regular.exf

@@ -1,6 +1,6 @@
 EX Version: 2
 Region: /
-!#nodeset datapoints
+!#nodeset nodes
 Shape. Dimension=0
 #Fields=1
 1) data_coordinates, coordinate, rectangular cartesian, real, #Components=3
@@ -1159,35 +1159,36 @@ Node: 288
   6.481555104255676e-02
  -9.071839042007923e-03
   4.956691563129425e-01
+!#nodeset datapoints
 #Fields=2
 1) marker_data_coordinates, coordinate, rectangular cartesian, real, #Components=3
  x. #Values=1 (value)
  y. #Values=1 (value)
  z. #Values=1 (value)
 2) marker_data_name, field, rectangular cartesian, string, #Components=1
  1. #Values=1 (value)
-Node: 1000001
+Node: 1
   0.000000000000000e+00
   0.000000000000000e+00
  -5.000000000000000e-01
  bottom
-Node: 1000002
+Node: 2
   0.000000000000000e+00
   0.000000000000000e+00
   5.000000000000000e-01
  top
-Node: 1000003
+Node: 3
  -5.000000000000000e-01
   0.000000000000000e+00
   0.000000000000000e+00
  left
-Node: 1000004
+Node: 4
   5.000000000000000e-01
   0.000000000000000e+00
   0.000000000000000e+00
  right
  Group name: bottom
-!#nodeset datapoints
+!#nodeset nodes
 Shape. Dimension=0
 #Fields=0
 Node: 1
@@ -1266,12 +1267,12 @@ Node: 94
 !#nodeset datapoints
 Shape. Dimension=0
 #Fields=0
-Node: 1000001
-Node: 1000002
-Node: 1000003
-Node: 1000004
+Node: 1
+Node: 2
+Node: 3
+Node: 4
  Group name: sides
-!#nodeset datapoints
+!#nodeset nodes
 Shape. Dimension=0
 #Fields=0
 Node: 51
@@ -1419,7 +1420,7 @@ Node: 234
 Node: 237
 Node: 238
  Group name: top
-!#nodeset datapoints
+!#nodeset nodes
 Shape. Dimension=0
 #Fields=0
 Node: 195