From 3cd66ed4a3902a1ea0cdf74c90cd07a51da8a609 Mon Sep 17 00:00:00 2001 From: Hugh Sorby Date: Sun, 9 Feb 2025 11:16:28 +1300 Subject: [PATCH 1/8] Add new API for getting Jacobian information about model. --- src/scaffoldfitter/fitter.py | 40 +++++++++++++++++++++++++++++++----- 1 file changed, 35 insertions(+), 5 deletions(-) diff --git a/src/scaffoldfitter/fitter.py b/src/scaffoldfitter/fitter.py index 8ae0d81..b2e0afe 100644 --- a/src/scaffoldfitter/fitter.py +++ b/src/scaffoldfitter/fitter.py @@ -7,8 +7,9 @@ from cmlibs.maths.vectorops import add, mult, sub from cmlibs.utils.zinc.field import assignFieldParameters, createFieldFiniteElementClone, getGroupList, \ findOrCreateFieldFiniteElement, findOrCreateFieldStoredMeshLocation, getUniqueFieldName, orphanFieldByName -from cmlibs.utils.zinc.finiteelement import evaluate_field_nodeset_range, findNodeWithName, getMaximumNodeIdentifier +from cmlibs.utils.zinc.finiteelement import evaluate_field_nodeset_range, findNodeWithName from cmlibs.utils.zinc.general import ChangeManager +from cmlibs.utils.zinc.mesh import calculate_jacobian, report_on_lowest_value from cmlibs.utils.zinc.region import write_to_buffer, read_from_buffer from cmlibs.zinc.context import Context from cmlibs.zinc.element import Elementbasis, Elementfieldtemplate @@ -186,7 +187,7 @@ def getInheritFitterStep(self, refFitterStep: FitterStep): """ refType = type(refFitterStep) for index in range(self._fitterSteps.index(refFitterStep) - 1, -1, -1): - if type(self._fitterSteps[index]) == refType: + if type(self._fitterSteps[index]) is refType: return self._fitterSteps[index] return None @@ -841,7 +842,7 @@ def getDataRMSAndMaximumProjectionErrorForGroup(self, groupName): with ChangeManager(self._fieldmodule): group = self._fieldmodule.findFieldByName(groupName).castGroup() if group.isValid(): - calculation_field = self._fieldmodule.createFieldAnd(group, self._activeDataGroupField) + calculation_field = self._fieldmodule.createFieldAnd(group, self._activeDataGroupField) nodeset = self._fieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_DATAPOINTS) temp_dataset_group = self._fieldmodule.createFieldGroup().createNodesetGroup(nodeset) temp_dataset_group.addNodesConditional(calculation_field) @@ -852,6 +853,37 @@ def getDataRMSAndMaximumProjectionErrorForGroup(self, groupName): return None, None + def getModelWorstElementJacobianInfo(self, mesh_group=None): + """ + Get the information on the 3D element with the worst jacobian value (most negative). + Optional mesh group parameter allows the user to make the calculation over a different group + from the whole mesh. + :param mesh_group: Optional parameter to specify a particular mesh group to make the calculation over. + :return: Element identifier, minimum jacobian value. Values are -1, inf if there is no data or bad fields. + """ + with ChangeManager(self._fieldmodule): + jacobian = calculate_jacobian(self._modelCoordinatesField) + report = report_on_lowest_value(jacobian, mesh_group if mesh_group else None) + + return report + + def getModelWorstElementJacobianInfoForGroup(self, group_name): + """ + Get the information on the 3D element with the worst jacobian + value (most negative) of the 3D mesh group with the given name. + If the group_name is not a valid group name then None, None is returned. + + :param group_name: Name of group to make calculation over. + :return: Element identifier, minimum jacobian value. + """ + with ChangeManager(self._fieldmodule): + group = self._fieldmodule.findFieldByName(group_name).castGroup() + if group.isValid(): + result = self.getMeshWorstElementJacobianInfo(group) + return result + + return None, None + def getMarkerDataFields(self): """ Only call if markerGroup exists. @@ -1096,7 +1128,6 @@ def _discoverModelFitGroup(self): Discover modelFitGroup from set name. """ self._modelFitGroup = None - field = None if self._modelFitGroupName: field = self._fieldmodule.findFieldByName(self._modelFitGroupName) if field.isValid(): @@ -1161,7 +1192,6 @@ def _discoverFlattenGroup(self): Discover flattenGroup from set name. """ self._flattenGroup = None - field = None if self._flattenGroupName: field = self._fieldmodule.findFieldByName(self._flattenGroupName) if field.isValid(): From 6f9991547d1fc10ba4f95c4e3f72a13ea208ac7b Mon Sep 17 00:00:00 2001 From: Hugh Sorby Date: Mon, 10 Feb 2025 14:05:25 +1300 Subject: [PATCH 2/8] Improve element jacobian calculation. --- src/scaffoldfitter/fitter.py | 24 ++++++++++++++---------- 1 file changed, 14 insertions(+), 10 deletions(-) diff --git a/src/scaffoldfitter/fitter.py b/src/scaffoldfitter/fitter.py index b2e0afe..0742582 100644 --- a/src/scaffoldfitter/fitter.py +++ b/src/scaffoldfitter/fitter.py @@ -853,9 +853,11 @@ def getDataRMSAndMaximumProjectionErrorForGroup(self, groupName): return None, None - def getModelWorstElementJacobianInfo(self, mesh_group=None): + def getLowestElementJacobian(self, mesh_group=None): """ Get the information on the 3D element with the worst jacobian value (most negative). + For a right-handed element values <= 0.0 are bad, the opposite holds for left-handed + elements. Optional mesh group parameter allows the user to make the calculation over a different group from the whole mesh. :param mesh_group: Optional parameter to specify a particular mesh group to make the calculation over. @@ -863,24 +865,26 @@ def getModelWorstElementJacobianInfo(self, mesh_group=None): """ with ChangeManager(self._fieldmodule): jacobian = calculate_jacobian(self._modelCoordinatesField) - report = report_on_lowest_value(jacobian, mesh_group if mesh_group else None) + report = report_on_lowest_value(jacobian, mesh_group) + del jacobian return report - def getModelWorstElementJacobianInfoForGroup(self, group_name): + def getLowestElementJacobianForGroup(self, group_name): """ - Get the information on the 3D element with the worst jacobian + Get the information on the 3D element with the worst Jacobian value (most negative) of the 3D mesh group with the given name. If the group_name is not a valid group name then None, None is returned. + If the fields for the calculation of the Jacobian are invalid then + -1, infinity is returned. :param group_name: Name of group to make calculation over. - :return: Element identifier, minimum jacobian value. + :return: Element identifier, minimum Jacobian value. """ - with ChangeManager(self._fieldmodule): - group = self._fieldmodule.findFieldByName(group_name).castGroup() - if group.isValid(): - result = self.getMeshWorstElementJacobianInfo(group) - return result + group = self._fieldmodule.findFieldByName(group_name).castGroup() + if group.isValid(): + result = self.getLowestElementJacobian(group) + return result return None, None From cab3318c0b0b01b36ea1b56644e7185cb532a252 Mon Sep 17 00:00:00 2001 From: Hugh Sorby Date: Mon, 10 Feb 2025 15:03:45 +1300 Subject: [PATCH 3/8] Tidy up getLowestElementJacobianForGroup. --- src/scaffoldfitter/fitter.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/scaffoldfitter/fitter.py b/src/scaffoldfitter/fitter.py index 0742582..341420d 100644 --- a/src/scaffoldfitter/fitter.py +++ b/src/scaffoldfitter/fitter.py @@ -883,8 +883,7 @@ def getLowestElementJacobianForGroup(self, group_name): """ group = self._fieldmodule.findFieldByName(group_name).castGroup() if group.isValid(): - result = self.getLowestElementJacobian(group) - return result + return self.getLowestElementJacobian(group) return None, None From 9548b23f2943c5e9692e3b0c9291bcd5601f4a81 Mon Sep 17 00:00:00 2001 From: Hugh Sorby Date: Tue, 11 Feb 2025 14:40:41 +1300 Subject: [PATCH 4/8] Update to changed API for Jacobian fields. --- src/scaffoldfitter/fitter.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/scaffoldfitter/fitter.py b/src/scaffoldfitter/fitter.py index 341420d..937cf05 100644 --- a/src/scaffoldfitter/fitter.py +++ b/src/scaffoldfitter/fitter.py @@ -9,7 +9,7 @@ findOrCreateFieldFiniteElement, findOrCreateFieldStoredMeshLocation, getUniqueFieldName, orphanFieldByName from cmlibs.utils.zinc.finiteelement import evaluate_field_nodeset_range, findNodeWithName from cmlibs.utils.zinc.general import ChangeManager -from cmlibs.utils.zinc.mesh import calculate_jacobian, report_on_lowest_value +from cmlibs.utils.zinc.field import get_element_jacobian_field, get_scalar_field_minimum_in_mesh from cmlibs.utils.zinc.region import write_to_buffer, read_from_buffer from cmlibs.zinc.context import Context from cmlibs.zinc.element import Elementbasis, Elementfieldtemplate @@ -864,11 +864,11 @@ def getLowestElementJacobian(self, mesh_group=None): :return: Element identifier, minimum jacobian value. Values are -1, inf if there is no data or bad fields. """ with ChangeManager(self._fieldmodule): - jacobian = calculate_jacobian(self._modelCoordinatesField) - report = report_on_lowest_value(jacobian, mesh_group) + jacobian = get_element_jacobian_field(self._modelCoordinatesField) + result = get_scalar_field_minimum_in_mesh(jacobian, mesh_group) del jacobian - return report + return result def getLowestElementJacobianForGroup(self, group_name): """ From 47b852e72ec126663d24a28e3f6b3e8ceb18b68e Mon Sep 17 00:00:00 2001 From: Hugh Sorby Date: Wed, 12 Feb 2025 11:26:26 +1300 Subject: [PATCH 5/8] Update to changes in API to zinc.utils. --- src/scaffoldfitter/fitter.py | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/src/scaffoldfitter/fitter.py b/src/scaffoldfitter/fitter.py index 937cf05..da7863c 100644 --- a/src/scaffoldfitter/fitter.py +++ b/src/scaffoldfitter/fitter.py @@ -6,10 +6,9 @@ from cmlibs.maths.vectorops import add, mult, sub from cmlibs.utils.zinc.field import assignFieldParameters, createFieldFiniteElementClone, getGroupList, \ - findOrCreateFieldFiniteElement, findOrCreateFieldStoredMeshLocation, getUniqueFieldName, orphanFieldByName -from cmlibs.utils.zinc.finiteelement import evaluate_field_nodeset_range, findNodeWithName + findOrCreateFieldFiniteElement, findOrCreateFieldStoredMeshLocation, getUniqueFieldName, orphanFieldByName, create_jacobian_determinant_field +from cmlibs.utils.zinc.finiteelement import evaluate_field_nodeset_range, findNodeWithName, get_scalar_field_minimum_in_mesh from cmlibs.utils.zinc.general import ChangeManager -from cmlibs.utils.zinc.field import get_element_jacobian_field, get_scalar_field_minimum_in_mesh from cmlibs.utils.zinc.region import write_to_buffer, read_from_buffer from cmlibs.zinc.context import Context from cmlibs.zinc.element import Elementbasis, Elementfieldtemplate @@ -864,7 +863,7 @@ def getLowestElementJacobian(self, mesh_group=None): :return: Element identifier, minimum jacobian value. Values are -1, inf if there is no data or bad fields. """ with ChangeManager(self._fieldmodule): - jacobian = get_element_jacobian_field(self._modelCoordinatesField) + jacobian = create_jacobian_determinant_field(self._modelCoordinatesField) result = get_scalar_field_minimum_in_mesh(jacobian, mesh_group) del jacobian From dd2f33741a9d18cd8cdeedbb26406fc05efa9e2a Mon Sep 17 00:00:00 2001 From: Hugh Sorby Date: Wed, 12 Feb 2025 11:30:37 +1300 Subject: [PATCH 6/8] Get mesh group from group field for current mesh. --- src/scaffoldfitter/fitter.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/scaffoldfitter/fitter.py b/src/scaffoldfitter/fitter.py index da7863c..6126d74 100644 --- a/src/scaffoldfitter/fitter.py +++ b/src/scaffoldfitter/fitter.py @@ -882,7 +882,8 @@ def getLowestElementJacobianForGroup(self, group_name): """ group = self._fieldmodule.findFieldByName(group_name).castGroup() if group.isValid(): - return self.getLowestElementJacobian(group) + mesh_group = group.getMeshGroup(self._mesh) + return self.getLowestElementJacobian(mesh_group) return None, None From 9556aad5a713fdb5c34ebfb5320ec99f68472d92 Mon Sep 17 00:00:00 2001 From: Hugh Sorby Date: Thu, 13 Feb 2025 15:52:44 +1300 Subject: [PATCH 7/8] Use create_xi_reference_jacobian_determinant_field in getLowestElementJacobian. --- src/scaffoldfitter/fitter.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/scaffoldfitter/fitter.py b/src/scaffoldfitter/fitter.py index 6126d74..94227ff 100644 --- a/src/scaffoldfitter/fitter.py +++ b/src/scaffoldfitter/fitter.py @@ -6,7 +6,8 @@ from cmlibs.maths.vectorops import add, mult, sub from cmlibs.utils.zinc.field import assignFieldParameters, createFieldFiniteElementClone, getGroupList, \ - findOrCreateFieldFiniteElement, findOrCreateFieldStoredMeshLocation, getUniqueFieldName, orphanFieldByName, create_jacobian_determinant_field + findOrCreateFieldFiniteElement, findOrCreateFieldStoredMeshLocation, getUniqueFieldName, orphanFieldByName, \ + create_xi_reference_jacobian_determinant_field from cmlibs.utils.zinc.finiteelement import evaluate_field_nodeset_range, findNodeWithName, get_scalar_field_minimum_in_mesh from cmlibs.utils.zinc.general import ChangeManager from cmlibs.utils.zinc.region import write_to_buffer, read_from_buffer @@ -863,7 +864,7 @@ def getLowestElementJacobian(self, mesh_group=None): :return: Element identifier, minimum jacobian value. Values are -1, inf if there is no data or bad fields. """ with ChangeManager(self._fieldmodule): - jacobian = create_jacobian_determinant_field(self._modelCoordinatesField) + jacobian = create_xi_reference_jacobian_determinant_field(self._modelCoordinatesField) result = get_scalar_field_minimum_in_mesh(jacobian, mesh_group) del jacobian From 4e3d1885b5605d122ef84f71b2b0037d6db439aa Mon Sep 17 00:00:00 2001 From: Hugh Sorby Date: Thu, 13 Feb 2025 15:57:34 +1300 Subject: [PATCH 8/8] Set minimum cmlibs.utils required to 0.11.2. --- pyproject.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyproject.toml b/pyproject.toml index a4e1501..af9c92a 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -16,7 +16,7 @@ authors = [ ] dependencies = [ "cmlibs.maths >= 0.6.2", - "cmlibs.utils >= 0.10.0", + "cmlibs.utils >= 0.11.2", "cmlibs.zinc >= 4.1" ] description = "Scaffold/model geometric fitting library using Zinc."