Closes #162 and adds an assert for non-symmetric Us in computePhi

tests/test_deformation.py

@@ -36,6 +36,11 @@ class testAll(unittest.TestCase):
 
 
     def test_computePhi(self):
+        # U MUST be symmetric, make sure this is checked and sets off an assertion
+        badU = numpy.eye(3)
+        badU[2,1] = 5
+        self.assertRaises(AssertionError, spam.deformation.computePhi, {'U': badU})
+
         trans1 = {'t': [0.0, 3.0, 3.0]}
         trans2 = {'r': [-5.0, 0.0, 0.0]}
         trans3 = {'z': [2, 2, 2]}
@@ -336,8 +341,6 @@ class testAll(unittest.TestCase):
         self.assertAlmostEqual(decomposedF['vol'][1:-1,1:-1,1:-1].mean(), 3. * (shrink - 1.), places=3)
 
 
-
-
     def test_geers(self):
         #nodeSpacing = [150,150,150]
         nodeSpacing = [200,200,200]
@@ -506,6 +509,7 @@ class testAll(unittest.TestCase):
         Ffield = spam.deformation.FfieldRegularGeers(displacements, nodeSpacing, mask=True)
         self.assertEqual(numpy.isnan(Ffield[:,:,:,0,0]).sum(), 1)
 
+
     def test_bagiStrain(self):
         #######################################################
         # Check triangulation and strains
@@ -688,6 +692,7 @@ class testAll(unittest.TestCase):
         self.assertEqual(decomposedF['dev'][otherTets].tolist(), numpy.zeros(otherTets.sum()).tolist())
         self.assertTrue(numpy.sum(decomposedF['dev'][tetsWithPointToMove] != 0) == len(tetsWithPointToMove))
 
+
     def test_merge(self):
         #######################################################
         ### We're using the DDIC test from scripts here, lightly modified
@@ -841,7 +846,7 @@ class testAll(unittest.TestCase):
         interpCoods = numpy.array([[50.,50.,50.]])
 
         interpolatedPhi = spam.deformation.interpolatePhiField(fieldCoords, fieldValues, interpCoods)
-        
+
         decomposedInterpolatedPhi = spam.deformation.decomposePhi(interpolatedPhi[0])
 
         for key in transformation.keys():
@@ -849,6 +854,7 @@ class testAll(unittest.TestCase):
                 self.assertAlmostEqual(decomposedInterpolatedPhi[key][i],
                                        transformation[key][i]/2.0 ,places=3)
 
+
     def test_computeRigidPhi(self):
         transformation = {'t': [5., 0., 0.],
                           'r': [3., 0., 0.]}
@@ -869,7 +875,8 @@ class testAll(unittest.TestCase):
                     self.assertEqual(PhiR1[i,j] > PhiR2[i,j], True)
                 else:
                     self.assertAlmostEqual(PhiR1[i,j], PhiR2[i,j], places=2)
-                    
+
+
     def test_decomposeF(self):
         #Test that it runs without problems
         F = numpy.array(([5, 2, 3],[2,5,1],[3, 1, 5]))

tools/deformation/deformationField.py

@@ -128,7 +128,7 @@ def FfieldRegularQ8(displacementField, nodeSpacing):
                 bar.update(nodesDone)
                 nodesDone += 1
     bar.finish()
-    
+
     return Ffield
 
 
@@ -387,7 +387,7 @@ def decomposeFfield(Ffield, components, twoD=False):
 
     Parameters
     ----------
-        Ffield : multidimensional x 3 x 3 numpy array
+        Ffield : multidimensional x 3 x 3 numpy array of floats
             Spatial field of Fs
 
         components : list of strings
@@ -410,9 +410,9 @@ def decomposeFfield(Ffield, components, twoD=False):
     import progressbar
     pbar = progressbar.ProgressBar()
     # The last two are for the 3x3 F field
-    fieldDimensions = Ffield.shape[0:-2]
+    fieldDimensions  = Ffield.shape[0:-2]
     fieldRavelLength = numpy.prod(numpy.array(fieldDimensions))
-    FfieldFlat = Ffield.reshape(fieldRavelLength, 3, 3)
+    FfieldFlat       = Ffield.reshape(fieldRavelLength, 3, 3)
 
     output = {}
     for component in components:
@@ -439,6 +439,66 @@ def decomposeFfield(Ffield, components, twoD=False):
     return output
 
 
+def decomposePhiField(PhiField, components, twoD=False):
+    """
+    This function takes in a Phi field (from readCorrelationTSV?) and
+    returns fields of desired transformation components.
+
+    Parameters
+    ----------
+        PhiField : multidimensional x 4 x 4 numpy array of floats
+            Spatial field of Phis
+
+        components : list of strings
+            These indicate the desired components consistent with spam.deformation.decomposePhi
+
+        twoD : bool, optional
+            Is the PhiField in 2D? This changes the strain calculation.
+            Default = False
+
+    Returns
+    -------
+        Dictionary containing appropriately reshaped fields of the transformation components requested.
+
+        Keys:
+            vol, dev, volss, devss are scalars
+            t, r, and z are 3-component vectors
+            e and U are 3x3 tensors
+    """
+    import spam.deformation
+    import progressbar
+    pbar = progressbar.ProgressBar()
+
+    # The last two are for the 4x4 Phi field
+    fieldDimensions  = PhiField.shape[0:-2]
+    fieldRavelLength = numpy.prod(numpy.array(fieldDimensions))
+    PhiFieldFlat     = PhiField.reshape(fieldRavelLength, 4, 4)
+
+    output = {}
+    for component in components:
+        output[component] = []
+
+    # Iterate through flat field of Fs
+    for n in pbar(range(PhiFieldFlat.shape[0])):
+        Phi = PhiFieldFlat[n]
+        decomposedPhi = spam.deformation.decomposePhi(Phi, twoD=twoD)
+        for component in components:
+            output[component].append(decomposedPhi[component])
+
+    # Reshape on the output
+    for component in components:
+        if component == 'vol' or component == 'dev' or component == 'volss' or component == 'devss':
+            output[component] = numpy.array(output[component]).reshape(*PhiField.shape[0:-2])
+
+        if component == 't' or component == 'r' or component == 'z':
+            output[component] = numpy.array(output[component]).reshape(*PhiField.shape[0:-2], 3)
+
+        if component == 'U' or component == 'e':
+            output[component] = numpy.array(output[component]).reshape(*PhiField.shape[0:-2], 3, 3)
+
+    return output
+
+
 def correctPhiField(fileName=None, fieldCoords=None, fieldValues=None, fieldRS=None, fieldDPhi=None, fieldPSCC=None, fieldIT=None, fieldBinRatio=1.0, ignoreBadPoints=False, ignoreBackGround=False, correctBadPoints=False, deltaPhiNormMin=0.001, pixelSearchCCmin=0.98, nNeighbours=12, filterPoints=False, filterPointsRadius=3, verbose=False, saveFile=False, saveFileName=None):
     """
     This function corrects a field of deformation functions **Phi** calculated at a number of points.

tools/deformation/deformationFunction.py

@@ -71,6 +71,8 @@ def computePhi(transformation, PhiCentre=[0.0, 0.0, 0.0], PhiPoint=[0.0, 0.0, 0.
 
     # Stretch tensor overrides 'z' and 's', hence elif next
     if 'U' in transformation:
+        # symmetric check from https://stackoverflow.com/questions/42908334/checking-if-a-matrix-is-symmetric-in-numpy
+        assert numpy.allclose(transformation['U'], transformation['U'].T), 'U not symmetric'
         tmp = numpy.eye(4)
         tmp[:3, :3] = transformation['U']
         Phi = numpy.dot(tmp, Phi)