Transformations in ScaffoldPackage

src/scaffoldmaker/meshtypes/meshtype_3d_heartatria1.py

@@ -1677,7 +1677,8 @@ def generateBaseMesh(cls, region, options):
             ractProportions[nc][0], ractProportions[nc][1],
             ravtProportions[ran1raap][0], ravtProportions[ran1raap][1],
             elementsCount = elementsCountOverSideRightAtriumPouch,
-            derivativeStart = [ (ractd2[1][nc][c] - ractd1[1][nc][c]) for c in range(3) ],
+            #derivativeStart = [ (ractd2[1][nc][c] - ractd1[1][nc][c]) for c in range(3) ],
+            derivativeStart = [ -d for d in ractd1[1][nc] ],
             derivativeEnd = [ -d for d in ravtd2[1][ran1raap] ])
         # get inner points
         raapx  = [ [ None ], raapx  ]
@@ -1708,7 +1709,8 @@ def generateBaseMesh(cls, region, options):
             ractProportions[nc][0], ractProportions[nc][1],
             ravtProportions[ran1raaq][0], ravtProportions[ran1raaq][1],
             elementsCount = elementsCountOverAtria//2 + elementsCountOverCristaTerminalisAnterior - 3,
-            derivativeStart = [ (ractd2[1][nc][c] - 0.5*ractd1[1][nc][c]) for c in range(3) ],
+            #derivativeStart = [ (ractd2[1][nc][c] - 0.5*ractd1[1][nc][c]) for c in range(3) ],
+            derivativeStart = [ (ractd2[1][nc][c] -ractd1[1][nc][c]) for c in range(3) ],
             derivativeEnd = [ -d for d in ravtd2[1][ran1raaq] ])
         # get inner points
         raaqx  = [ [ None ], raaqx  ]
@@ -2972,7 +2974,8 @@ def generateBaseMesh(cls, region, options):
                     remapEftNodeValueLabel(eft1, [ 1, 5 ], Node.VALUE_LABEL_D_DS1, [ ( Node.VALUE_LABEL_D_DS1, [1] ) ])
                     remapEftNodeValueLabel(eft1, [ 1, 5 ], Node.VALUE_LABEL_D_DS2, [ ( Node.VALUE_LABEL_D_DS2, [1] ) ])
                     remapEftNodeValueLabel(eft1, [ 2, 6 ], Node.VALUE_LABEL_D_DS1, [ ( Node.VALUE_LABEL_D_DS1, [] ), ( Node.VALUE_LABEL_D_DS2, [] ) ])
-                    remapEftNodeValueLabel(eft1, [ 3, 7 ], Node.VALUE_LABEL_D_DS1, [ ( Node.VALUE_LABEL_D_DS1, [1] ), ( Node.VALUE_LABEL_D_DS2, [] ) ])
+                    #remapEftNodeValueLabel(eft1, [ 3, 7 ], Node.VALUE_LABEL_D_DS1, [ ( Node.VALUE_LABEL_D_DS1, [1] ), ( Node.VALUE_LABEL_D_DS2, [] ) ])
+                    remapEftNodeValueLabel(eft1, [ 3, 7 ], Node.VALUE_LABEL_D_DS1, [ ( Node.VALUE_LABEL_D_DS1, [1] ) ])
                     remapEftNodeValueLabel(eft1, [ 3, 7 ], Node.VALUE_LABEL_D_DS2, [ ( Node.VALUE_LABEL_D_DS2, [1] ) ])
                     addMarker = { "name" : "SVC-RA", "xi" : [ 0.0, 1.0, 0.0 ] }
                 elif (e2 == 1) and (e1 < (elementsCountOverSideRightAtriumPouch - 1)):
@@ -3072,7 +3075,8 @@ def generateBaseMesh(cls, region, options):
                 if n1 == 0:
                     raasDerivativesMap[n3][ix] = ( None, (-1, 0, 0), None, (0, 1, 0) )  # compare ( None, (-1, 1, 0), None, (0, 1, 0) )
                 else: # n1 == 1:
-                    raasDerivativesMap[n3][ix] = ( (0, 1, 0), (-1, -1, 0), None, (-1, 1, 0) )  # compare ( (0, 1, 0), (-1, 0, 0), None, (-1, 1, 0) )
+                    #raasDerivativesMap[n3][ix] = ( (0, 1, 0), (-1, -1, 0), None, (-1, 1, 0) )  # compare ( (0, 1, 0), (-1, 0, 0), None, (-1, 1, 0) )
+                    raasDerivativesMap[n3][ix] = ( (0, 1, 0), (-1, -1, 0), None, (-1, 0, 0) )
             ix += 1
         # down right atrial appendage pouch limit, raap
         n1Last = len(raapx[1]) - 1

src/scaffoldmaker/scaffoldpackage.py

@@ -4,6 +4,10 @@
 """
 
 import copy
+import math
+from opencmiss.utils.zinc.field import createFieldEulerAnglesRotationMatrix
+from opencmiss.utils.zinc.general import ChangeManager
+from opencmiss.utils.maths.vectorops import euler_to_rotation_matrix
 from scaffoldmaker.meshtypes.scaffold_base import Scaffold_base
 
 class ScaffoldPackage:
@@ -28,6 +32,16 @@ def __init__(self, scaffoldType, dct={}, defaultParameterSetName='Default'):
         if scaffoldSettings:
             # remove obsolete options? If so, deepcopy first?
             self._scaffoldSettings.update(scaffoldSettings)
+        # note rotation is stored in degrees
+        self._rotation = dct.get('rotation')
+        if not self._rotation:
+            self._rotation = [ 0.0, 0.0, 0.0 ]
+        self._scale = dct.get('scale')
+        if not self._scale:
+            self._scale = [ 1.0, 1.0, 1.0 ]
+        self._translation = dct.get('translation')
+        if not self._translation:
+            self._translation = [ 0.0, 0.0, 0.0 ]
         self._meshEdits = copy.deepcopy(dct.get('meshEdits'))
 
     def deepcopy(self, other):
@@ -56,8 +70,11 @@ def toDict(self):
         :return: Dictionary containing object encoding.
         '''
         dct = {
+            'rotation' : self._rotation,
             'scaffoldTypeName' : self._scaffoldType.getName(),
-            'scaffoldSettings' : self._scaffoldSettings
+            'scaffoldSettings' : self._scaffoldSettings,
+            'scale' : self._scale,
+            'translation' : self._translation
             }
         if self._meshEdits:
             dct['meshEdits'] = self._meshEdits
@@ -75,12 +92,116 @@ def getScaffoldSettings(self):
     def getScaffoldType(self):
         return self._scaffoldType
 
-    def generate(self, region):
+    def getRotation(self):
+        return self._rotation
+
+    def setRotation(self, rotation):
+        '''
+        :param rotation: list of 3 rotation angles about z, y', x'', in degrees.
+        :return: True if value changed, otherwise False.
+        '''
+        assert len(rotation) == 3
+        if self._rotation != rotation:
+            self._rotation = rotation
+            return True
+        return False
+
+    def getScale(self):
+        return self._scale
+
+    def setScale(self, scale):
+        '''
+        :param scale: list of 3 scales for x, y, z.
+        :return: True if value changed, otherwise False.
+        '''
+        assert len(scale) == 3
+        if self._scale != scale:
+            self._scale = scale
+            return True
+        return False
+
+    def getTranslation(self):
+        return self._translation
+
+    def setTranslation(self, translation):
+        '''
+        :param translation: list of 3 translation values for x, y, z.
+        :return: True if value changed, otherwise False.
+        '''
+        assert len(translation) == 3
+        if self._translation != translation:
+            self._translation = translation
+            return True
+        return False
+
+    def getTransformationMatrix(self):
+        '''
+        :return: 4x4 row-major transformation matrix with first index down rows, second across columns,
+        suitable for multiplication p' = Mp where p = [ x, y, z, h ], or None if identity.
+        '''
+        # apply transformation in order: scale then rotation then translation
+        if not all((v == 0.0) for v in self._rotation):
+            rotationMatrix = euler_to_rotation_matrix([ deg*math.pi/180.0 for deg in self._rotation ])
+            return [
+                [ rotationMatrix[0][0]*self._scale[0], rotationMatrix[0][1]*self._scale[1], rotationMatrix[0][2]*self._scale[2], self._translation[0] ],
+                [ rotationMatrix[1][0]*self._scale[0], rotationMatrix[1][1]*self._scale[1], rotationMatrix[1][2]*self._scale[2], self._translation[1] ],
+                [ rotationMatrix[2][0]*self._scale[0], rotationMatrix[2][1]*self._scale[1], rotationMatrix[2][2]*self._scale[2], self._translation[2] ],
+                [ 0.0, 0.0, 0.0, 1.0 ] ]
+        if not (all((v == 1.0) for v in self._scale) and all((v == 0.0) for v in self._translation)):
+            return [
+                [ self._scale[0], 0.0, 0.0, self._translation[0] ],
+                [ 0.0, self._scale[1], 0.0, self._translation[1] ],
+                [ 0.0, 0.0, self._scale[2], self._translation[2] ],
+                [ 0.0, 0.0, 0.0, 1.0 ] ]
+        return None
+
+    def applyTransformation(self, region):
+        '''
+        If rotation, scale or transformation are set, transform node coordinates.
+        '''
+        fieldmodule = region.getFieldmodule()
+        coordinates = fieldmodule.findFieldByName('coordinates').castFiniteElement()
+        if not coordinates.isValid():
+            print('Warning: ScaffoldPackage.applyTransformation: Missing coordinates field')
+            return
+        with ChangeManager(fieldmodule):
+            componentsCount = coordinates.getNumberOfComponents()
+            if componentsCount < 3:
+                # pad with zeros
+                coordinates = fieldmodule.createFieldConcatenate([ coordinates ] + [ fieldmodule.createFieldConstant([ 0.0 ]*(3 - componentsCount)) ])
+            newCoordinates = coordinates
+            # apply scale first so variable scaling in x, y, z doesn't interplay with rotation
+            if not all((v == 1.0) for v in self._scale):
+                #print("applyTransformation: apply scale", self._scale)
+                newCoordinates = newCoordinates*fieldmodule.createFieldConstant(self._scale)
+            if not all((v == 0.0) for v in self._rotation):
+                #print("applyTransformation: apply rotation", self._rotation)
+                newCoordinates = fieldmodule.createFieldMatrixMultiply(3,
+                    createFieldEulerAnglesRotationMatrix(fieldmodule, fieldmodule.createFieldConstant([ deg*math.pi/180.0 for deg in self._rotation ])), newCoordinates)
+            if not all((v == 0.0) for v in self._translation):
+                #print("applyTransformation: apply translation", self._translation)
+                newCoordinates = newCoordinates + fieldmodule.createFieldConstant(self._translation)
+            # be sure to delete temporary fields and fieldassignment to reduce messages
+            if newCoordinates is not coordinates:
+                fieldassignment = coordinates.createFieldassignment(newCoordinates)
+                fieldassignment.assign()
+                del fieldassignment
+            del newCoordinates
+            del coordinates
+
+    def generate(self, region, applyTransformation=True):
+        '''
+        :param applyTransformation: If True (default) apply scale, rotation and translation to
+        node coordinates. Specify False if client will transform, e.g. with graphics transformations.
+        '''
         #print('\nScaffoldPackage.generate: ', self.toDict())
         annotationGroups = self._scaffoldType.generateMesh(region, self._scaffoldSettings)
         if self._meshEdits:
             # apply mesh edits, a Zinc-readable model file containing node edits
+            # Note: these are untransformed coordinates
             sir = region.createStreaminformationRegion()
             srm = sir.createStreamresourceMemoryBuffer(self._meshEdits)
             region.read(sir)
+        if applyTransformation:
+            self.applyTransformation(region)
         return annotationGroups

tests/test_general.py

@@ -0,0 +1,82 @@
+import unittest
+from opencmiss.utils.maths.vectorops import magnitude
+from opencmiss.utils.zinc.finiteelement import evaluateFieldNodesetRange, findNodeWithName
+from opencmiss.zinc.context import Context
+from opencmiss.zinc.field import Field
+from opencmiss.zinc.node import Node
+from opencmiss.zinc.result import RESULT_OK
+from scaffoldmaker.meshtypes.meshtype_3d_box1 import MeshType_3d_box1
+from scaffoldmaker.scaffoldpackage import ScaffoldPackage
+from testutils import assertAlmostEqualList
+
+class GeneralScaffoldTestCase(unittest.TestCase):
+
+    def test_transformation(self):
+        """
+        Test transformation of a box scaffold with scaffold package.
+        """
+        scaffoldPackage = ScaffoldPackage(MeshType_3d_box1)
+
+        tmpScale = scaffoldPackage.getScale()
+        TOL = 1.0E-7
+        FINETOL = 1.0E-12
+        assertAlmostEqualList(self, tmpScale, [ 1.0, 1.0, 1.0 ], delta=FINETOL)
+        newScale = [ 2.0, 1.5, 0.5 ]
+        scaffoldPackage.setScale(newScale)
+        tmpScale = scaffoldPackage.getScale()
+        assertAlmostEqualList(self, tmpScale, newScale, delta=FINETOL)
+
+        tmpRotation = scaffoldPackage.getRotation()
+        assertAlmostEqualList(self, tmpRotation, [ 0.0, 0.0, 0.0 ], delta=FINETOL)
+        newRotation = [ 30.0, -10.0, 90.0 ]
+        scaffoldPackage.setRotation(newRotation)
+        tmpRotation = scaffoldPackage.getRotation()
+        assertAlmostEqualList(self, tmpRotation, newRotation, delta=FINETOL)
+
+        tmpTranslation = scaffoldPackage.getTranslation()
+        assertAlmostEqualList(self, tmpTranslation, [ 0.0, 0.0, 0.0 ], delta=FINETOL)
+        newTranslation = [ 0.5, 1.2, -0.1 ]
+        scaffoldPackage.setTranslation(newTranslation)
+        tmpTranslation = scaffoldPackage.getTranslation()
+        assertAlmostEqualList(self, tmpTranslation, newTranslation, delta=FINETOL)
+
+        context = Context("Test")
+        region = context.getDefaultRegion()
+        self.assertTrue(region.isValid())
+
+        scaffoldPackage.generate(region)
+
+        fieldmodule = region.getFieldmodule()
+        nodes = fieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES)
+        self.assertEqual(8, nodes.getSize())
+
+        coordinates = fieldmodule.findFieldByName("coordinates").castFiniteElement()
+        self.assertTrue(coordinates.isValid())
+        minimums, maximums = evaluateFieldNodesetRange(coordinates, nodes)
+        assertAlmostEqualList(self, minimums, [  2.744244002293470e-01,  6.367511648575830e-01, -1.000000000000000e-01 ], delta=TOL)
+        assertAlmostEqualList(self, maximums, [  2.455737063904887e+00,  2.184807753012208e+00,  1.724507984852172e+00 ], delta=TOL)
+
+        node = nodes.findNodeByIdentifier(8)
+        self.assertTrue(node.isValid())
+        fieldcache = fieldmodule.createFieldcache()
+        fieldcache.setNode(node)
+        result, x = coordinates.getNodeParameters(fieldcache, -1, Node.VALUE_LABEL_VALUE, 1, 3)
+        self.assertEqual(RESULT_OK, result)
+        assertAlmostEqualList(self, x , [  2.230161464134234e+00,  1.621558917869791e+00,  1.724507984852172e+00 ], delta=TOL)
+        # derivative magnitudes must also equal scale
+        result, d1 = coordinates.getNodeParameters(fieldcache, -1, Node.VALUE_LABEL_D_DS1, 1, 3)
+        self.assertEqual(RESULT_OK, result)
+        assertAlmostEqualList(self, d1, [  1.705737064039425e+00,  9.848077530127952e-01,  3.472963553408093e-01 ], delta=TOL)
+        self.assertAlmostEqual(newScale[0], magnitude(d1), delta=TOL)
+        result, d2 = coordinates.getNodeParameters(fieldcache, -1, Node.VALUE_LABEL_D_DS2, 1, 3)
+        self.assertEqual(RESULT_OK, result)
+        assertAlmostEqualList(self, d2, [ -2.255755995328457e-01, -1.302361332111701e-01,  1.477211629352659e+00 ], delta=TOL)
+        self.assertAlmostEqual(newScale[1], magnitude(d2), delta=TOL)
+        result, d3 = coordinates.getNodeParameters(fieldcache, -1, Node.VALUE_LABEL_D_DS3, 1, 3)
+        self.assertEqual(RESULT_OK, result)
+        assertAlmostEqualList(self, d3, [  2.499999998128999e-01, -4.330127019169794e-01,  0.000000000000000e+00 ], delta=TOL)
+        self.assertAlmostEqual(newScale[2], magnitude(d3), delta=TOL)
+
+
+if __name__ == "__main__":
+    unittest.main()