# -*- coding: utf-8 -*- from __future__ import print_function import unittest import numpy import os import scipy.ndimage import tifffile import subprocess import spam.helpers import spam.kalisphera import spam.datasets import spam.mesh import spam.label import spam.DIC import spam.deformation testFolder = './' class testAll(unittest.TestCase): def tearDown(self): try: os.remove(testFolder+'Step0.tif') os.remove(testFolder+'Lab0.tif') os.remove(testFolder+'Step1.tif') os.remove(testFolder+'Step0-Step1-bin2-registration.tsv') os.remove(testFolder+'Step0-Step1.tsv') os.remove(testFolder+'Step0-Step1-discreteDVC.tsv') os.remove(testFolder+'Step0-Step1-discreteDVC.vtk') os.remove(testFolder+'merged.tsv') except OSError: pass 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]} trans4 = {'s': [0.9, 0.8, 0.7]} Phi1 = spam.deformation.computePhi(trans1) self.assertEqual(numpy.sum([Phi1[0, -1], Phi1[1, -1], Phi1[2, -1]]), 6) Phi2 = spam.deformation.computePhi(trans2) self.assertEqual(numpy.sum([Phi2[0, -1], Phi2[1, -1], Phi2[2, -1]]), 0) Phi3 = spam.deformation.computePhi(trans2, PhiCentre=[50.0, 50.0, 50.0], PhiPoint=[50.0, 16.0, 84.0]) self.assertAlmostEqual(numpy.sum([Phi3[0, -1], Phi3[1, -1], Phi3[2, -1]]), 5.926, places=2) Phi4 = spam.deformation.computePhi(trans3) self.assertEqual([Phi4[0, 0], Phi4[1, 1], Phi4[2, 2]], [2., 2., 2.]) Phi5 = spam.deformation.computePhi(trans4) self.assertEqual(Phi5[0, 1], Phi5[1, 0], 0.9) self.assertEqual(Phi5[0, 2], Phi5[2, 0], 0.8) self.assertEqual(Phi5[1, 3], Phi5[3, 1], 0.7) Phi6 = spam.deformation.computePhi({'U':numpy.eye(3)}) self.assertEqual(Phi6.tolist(), numpy.eye(4).tolist()) U = numpy.eye(3) U[0,0] = 0.8 Phi7 = spam.deformation.computePhi({'U':U}) self.assertEqual(Phi7[0,0], 0.8) def test_decomposePhi(self): # CASE 1: singular Phi Phi = numpy.zeros((4,4)) self.assertEqual(numpy.isnan(spam.deformation.decomposePhi(Phi)["t"]).sum(), 3) self.assertEqual(numpy.isnan(spam.deformation.decomposePhi(Phi)["r"]).sum(), 3) # CASE 2: Phi contains nans Phin = numpy.eye(4) Phin[0, 3] = numpy.nan self.assertEqual(numpy.isnan(spam.deformation.decomposePhi(Phin)["t"]).sum(), 3) self.assertEqual(numpy.isnan(spam.deformation.decomposePhi(Phin)["r"]).sum(), 3) # CASE 3: Phi contains inf Phii = numpy.eye(4) Phii[0, 3] = numpy.inf self.assertEqual(numpy.isnan(spam.deformation.decomposePhi(Phii)["t"]).sum(), 3) self.assertEqual(numpy.isnan(spam.deformation.decomposePhi(Phii)["r"]).sum(), 3) # CASE 4: negative eigenvalues in transpose(Phi).Phi Phib = numpy.array([[0.25383295079467, -0.028793281703941, 0.05805266499330, 10.722113444672508], [-0.03543852689277, 0.004282141971071, -0.008072359472034, 18.14869851125678], [0.041196242920521, -0.004961945218488, 0.009385858417623, 17.144152186875889], [0, 0, 0, 1]]) self.assertAlmostEqual(spam.deformation.decomposePhi(Phib)["t"].tolist(), Phib[0:3, -1].tolist(), places=5) # CASE 5: error in inverting U to create the rotation matrix # it doesn't seem to pass the test on gricad... # Phib = numpy.array([[0.000000001, 0.000000001, 0.000000001, 1], # [0.000000001, 0.000000001, 0.000000001, 1], # [0.000000001, 0.0, 0.000000001, 1], # [0, 0, 0, 1]]) # self.assertEqual(spam.deformation.decomposePhi(Phib)["t"], [0,0,0]) # CASE 6: back calculation of transformation transIn = {'t': [4.6, 11.2, 0.3], 'r': [-5.0, 2.0, 8.0]} a = numpy.random.uniform(1, 100, 3) b = numpy.random.uniform(1, 100, 3) Phi = spam.deformation.computePhi(transIn, PhiCentre=a, PhiPoint=b) transBack = spam.deformation.decomposePhi(Phi, PhiCentre=b, PhiPoint=a) self.assertTrue(numpy.linalg.norm(numpy.subtract(transBack["t"], transIn["t"])) < 1e-5) self.assertTrue(numpy.linalg.norm(numpy.subtract(transBack["r"], transIn["r"])) < 1e-5) def test_Q8(self): # case 0: 2D field with nans under large strains dispField0 = numpy.zeros((4, 3)) dispField0[:, 1] = [20, 20, -20, numpy.nan] # dipsY dispField0[:, 2] = [20, -20, 20, numpy.nan] # dipsX F0 = spam.deformation.FfieldRegularQ8(dispField0.reshape(1, 2, 2, 3), nodeSpacing=[0, 100, 100]) decomposedF = spam.deformation.decomposeF(F0) self.assertTrue(numpy.isnan((decomposedF['U']).sum())) # return nan ## case 0b: 2D field with nans under small strains self.assertTrue(numpy.isnan((decomposedF['e']).sum())) # return nan # case 1: 2D field of 1 square with isotropic compression under large strains dispField1 = numpy.zeros((4, 3)) dispField1[:, 1] = [20, 20, -20, -20] # dipsY dispField1[:, 2] = [20, -20, 20, -20] # dipsX Vo1 = 100 * 100 # initial volume Vf1 = 60 * 60 # final volume F1 = spam.deformation.FfieldRegularQ8(dispField1.reshape(1, 2, 2, 3), nodeSpacing=[0, 100, 100]) decomposedF1 = spam.deformation.decomposeF(F1[0, 0, 0], twoD = True) self.assertAlmostEqual(decomposedF1["vol"], (Vf1 - Vo1) / float(Vo1), places=3) # volStrain should be the volume change self.assertAlmostEqual(decomposedF1["dev"], 0, places=6) # devStrain must be 0 # case 1b: 2D field of 1 rectangle with isotropic compression under small strains dispField1b = numpy.zeros((4, 3)) dispField1b[:, 1] = [1, 1, -1, -1] # dipsY dispField1b[:, 2] = [1, -1, 1, -1] # dipsX Vo1b = 100 * 100 # initial volume Vf1b = 98 * 98 # final volume F1b = spam.deformation.FfieldRegularQ8(dispField1b.reshape(1, 2, 2, 3), nodeSpacing=[0, 100, 100]) decomposedF1b = spam.deformation.decomposeF(F1b[0, 0, 0], twoD = True) self.assertAlmostEqual(decomposedF1b["volss"], (Vf1b - Vo1b) / float(Vo1b), places=3) # volStrain should be the volume change self.assertEqual(decomposedF1b["devss"], 0) # devStrain must be 0 # case 2: 3D field of 1 cell (2x2x2nodes) with isotropic compression under large strains dispField2 = numpy.zeros((8, 3)) dispField2[:, 0] = [20, 20, 20, 20, -20, -20, -20, -20] # dipsZ dispField2[:, 1] = [20, 20, -20, -20, 20, 20, -20, -20] # dipsY dispField2[:, 2] = [20, -20, 20, -20, 20, -20, 20, -20] # dipsX Vo2 = 100 * 100 * 100 # initial volume Vf2 = 60 * 60 * 60 # final volume F2 = spam.deformation.FfieldRegularQ8(dispField2.reshape(2, 2, 2, 3), nodeSpacing=[100, 100, 100]) decomposedF2 = spam.deformation.decomposeF(F2[0, 0, 0]) self.assertAlmostEqual(decomposedF2["vol"], (Vf2 - Vo2) / float(Vo2), places=3) # volStrain should be the volume change self.assertAlmostEqual(decomposedF2["dev"], 0, places=6) # devStrain must be 0 # case 2b: 3D field of 1 cell (2x2x2nodes) with isotropic compression under small strains dispField2b = numpy.zeros((8, 3)) dispField2b[:, 0] = [1, 1, 1, 1, -1, -1, -1, -1] # dipsZ dispField2b[:, 1] = [1, 1, -1, -1, 1, 1, -1, -1] # dipsY dispField2b[:, 2] = [1, -1, 1, -1, 1, -1, 1, -1] # dipsX Vo2b = 100 * 100 * 100 # initial volume Vf2b = 98 * 98 * 98 # final volume F2b = spam.deformation.FfieldRegularQ8(dispField2b.reshape(2, 2, 2, 3), nodeSpacing=[100, 100, 100]) decomposedF2b = spam.deformation.decomposeF(F2b[0, 0, 0]) self.assertAlmostEqual(decomposedF2b["volss"], (Vf2b - Vo2b) / float(Vo2b), places=2) self.assertEqual(decomposedF2b["devss"], 0) # case 3: 2D field of 1 rectangle with shear under large strains dispField3 = numpy.zeros((4, 3)) dispField3[:, 2] = [0, 0, 20, 20] # dipsX F3 = spam.deformation.FfieldRegularQ8(dispField3.reshape(1, 2, 2, 3), nodeSpacing=[0, 100, 100]) decomposedF3 = spam.deformation.decomposeF(F3[0, 0, 0], twoD = True) self.assertEqual(decomposedF3["vol"], 0) # volStrain must be 0 self.assertAlmostEqual(decomposedF3["U"][1, 2], decomposedF3["U"][2, 1], places=4) # diag strain matrix self.assertAlmostEqual(abs(decomposedF3["U"][1, 2] - ((20 / 100.0) * 0.5)), 0, places=3) # case 3b: 2D field of 1 rectangle with shear under small strains dispField3b = numpy.zeros((4, 3)) dispField3b[:, 2] = [0, 0, 2, 2] # dipsX F3b = spam.deformation.FfieldRegularQ8(dispField3b.reshape(1, 2, 2, 3), nodeSpacing=[0, 100, 100]) decomposedF3b = spam.deformation.decomposeF(F3b[0, 0, 0], twoD = True) self.assertEqual(decomposedF3b["volss"], 0) # volStrain must be 0 self.assertEqual(decomposedF3b["e"][1, 2], decomposedF3b["e"][2, 1]) # diag strain matrix self.assertEqual(decomposedF3b["e"][1, 2], (2 / 100.0) * 0.5) # case 4: 3D field with shear under large strains dispField4 = numpy.zeros((8, 3)) dispField4[:, 2] = [0, 0, 20, 20, 0, 0, 20, 20] # dipsX F4 = spam.deformation.FfieldRegularQ8(dispField4.reshape(2, 2, 2, 3), nodeSpacing=[100, 100, 100]) decomposedF4 = spam.deformation.decomposeF(F4[0, 0, 0]) self.assertEqual(decomposedF4["vol"], 0) # volStrain must be 0 self.assertAlmostEqual(decomposedF4["U"][1, 2], decomposedF4["U"][2, 1], places=4) # diag strain matrix self.assertAlmostEqual(abs(decomposedF4["U"][1, 2] - (20 / 100.0) * 0.5), 0, places=3) # case 4b: 3D field with shear under small strains dispField4b = numpy.zeros((8, 3)) dispField4b[:, 2] = [0, 0, 2, 2, 0, 0, 2, 2] # dipsX F4b = spam.deformation.FfieldRegularQ8(dispField4b.reshape(2, 2, 2, 3), nodeSpacing=[100, 100, 100]) decomposedF4b = spam.deformation.decomposeF(F4b[0, 0, 0]) self.assertEqual(decomposedF4b["volss"], 0) # volStrain must be 0 self.assertEqual(decomposedF4b["e"][1, 2], decomposedF4b["e"][2, 1]) # diag strain matrix self.assertEqual(decomposedF4b["e"][1, 2], (2 / 100.0) * 0.5) #nodeSpacing = [150,150,150] nodeSpacing = [200,200,200] imSize = [2000,1000,1000] # Random points pointsRef, dims = spam.DIC.grid.makeGrid(imSize, nodeSpacing=nodeSpacing) ####################################################### # Dilation of points for grain strain calc ####################################################### # Test with 1% dilation dilation = 1.01 # make sure average volumetric strain is the same as the projected one pointsDef = pointsRef * dilation displacements = (pointsDef - pointsRef).reshape(dims[0], dims[1], dims[2], 3) Ffield = spam.deformation.FfieldRegularQ8(displacements, nodeSpacing) ev = spam.deformation.decomposeFfield(Ffield, ["vol"])["vol"] # print("\tU trace mean new:",numpy.array(U)[:,[0,1,2],[0,1,2]].mean()) # print("\tev",numpy.array(ev)) # print("\tev",numpy.array(ev).mean()) self.assertAlmostEqual(ev.mean(), 3. * (dilation - 1), places=3) #self.assertAlmostEqual(ev.mean(), 3. * (dilation - 1), places=3) ######################################################## ## Z-stretch of points for strain calc ######################################################## zStretch = 1.02 print("\n\nSpreading points in Z by by 1.2") # Deform the points by spreading them apart -- this is dilating so positive volumetric strain displacements = (pointsRef * [zStretch, 1.0, 1.0] - pointsRef).reshape(dims[0], dims[1], dims[2], 3) Ffield = spam.deformation.FfieldRegularQ8(displacements, nodeSpacing) # only compute on internal nodes, we know there's a problem at the edges self.assertAlmostEqual(Ffield[1:-1,1:-1,1:-1,0,0].mean(), zStretch, places=3) self.assertAlmostEqual(Ffield[1:-1,1:-1,1:-1,1,1].mean(), 1.0, places=3) self.assertAlmostEqual(Ffield[1:-1,1:-1,1:-1,2,2].mean(), 1.0, places=3) ######################################################## ## Y-stretch of points for strain calc ######################################################## yStretch = 1.02 print("\n\nSpreading points in Y by by 1.2") # Deform the points by spreading them apart -- this is dilating so positive volumetric strain displacements = (pointsRef * [1.0, yStretch, 1.0] - pointsRef).reshape(dims[0], dims[1], dims[2], 3) Ffield = spam.deformation.FfieldRegularQ8(displacements, nodeSpacing) # only compute on internal nodes, we know there's a problem at the edges self.assertAlmostEqual(Ffield[1:-1,1:-1,1:-1,0,0].mean(), 1.0, places=3) self.assertAlmostEqual(Ffield[1:-1,1:-1,1:-1,1,1].mean(), yStretch, places=3) self.assertAlmostEqual(Ffield[1:-1,1:-1,1:-1,2,2].mean(), 1.0, places=3) ####################################################### # Rotation of points for grain strain calc ####################################################### pointsRotated = pointsRef.copy() rotAngle = 3.0 for n, point in enumerate(pointsRef): Phi = spam.deformation.computePhi({'r': [rotAngle, 0.0, 0.0]}, PhiCentre = numpy.array(((numpy.array(imSize)-1)/2)-1), PhiPoint = point) pointsRotated[n] += Phi[0:3, -1] # print( point, pointsRotated[n], '\n\n' ) displacements = (pointsRotated - pointsRef).reshape(dims[0], dims[1], dims[2], 3) Ffield = spam.deformation.FfieldRegularQ8(displacements, nodeSpacing) # Compute strains decomposedF = spam.deformation.decomposeFfield(Ffield, ["r", "vol", "dev"]) self.assertAlmostEqual(decomposedF['r'][1:-1,1:-1,1:-1, 0].mean(), rotAngle, places=3) self.assertAlmostEqual(decomposedF['r'][1:-1,1:-1,1:-1, 1].mean(), 0.0, places=3) self.assertAlmostEqual(decomposedF['r'][1:-1,1:-1,1:-1, 2].mean(), 0.0, places=3) self.assertAlmostEqual(decomposedF['vol'][1:-1,1:-1,1:-1].mean(), 0, places=3) self.assertAlmostEqual(decomposedF['dev'][1:-1,1:-1,1:-1].mean(), 0, places=3) ####################################################### # Shear of points for grain strain calc: # N.B.: this is a shear which results in a SYMMETRIC F ####################################################### pointsRotated = pointsRef.copy() shearVal = 0.1 for n, point in enumerate(pointsRef): Phi = spam.deformation.computePhi({'s': [shearVal, 0.0, 0.0]}, PhiCentre = numpy.array(((numpy.array(imSize)-1)/2)-1), PhiPoint = point) pointsRotated[n] += Phi[0:3, -1] # print( point, pointsRotated[n], '\n\n' ) displacements = (pointsRotated - pointsRef).reshape(dims[0], dims[1], dims[2], 3) Ffield = spam.deformation.FfieldRegularQ8(displacements, nodeSpacing) # Compute strains decomposedF = spam.deformation.decomposeFfield(Ffield, ["U"]) UmeanWithoutBoundaries = numpy.mean(decomposedF['U'][1:-1,1:-1,1:-1], axis=(0, 1, 2)) for i in range(3): for j in range(i,3): self.assertAlmostEqual(UmeanWithoutBoundaries[i,j], UmeanWithoutBoundaries[j,i], places=3) ####################################################### # Homogeneous shrinking together ####################################################### shrink = 0.99 # Deform the points by pushing them closer together -- this is compressing so negative volumetric strain displacements = (pointsRef * shrink - pointsRef).reshape(dims[0], dims[1], dims[2], 3) Ffield = spam.deformation.FfieldRegularQ8(displacements, nodeSpacing) decomposedF = spam.deformation.decomposeFfield(Ffield, ["vol"]) self.assertAlmostEqual(decomposedF['vol'][1:-1,1:-1,1:-1].mean(), 3. * (shrink - 1.), places=3) ### while we're at it, the brute force distance too, why not here Ffield = spam.deformation.FfieldRegularQ8(displacements, nodeSpacing) decomposedF = spam.deformation.decomposeFfield(Ffield, ["vol"]) 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] imSize = [2000,1000,1000] # Random points pointsRef, dims = spam.DIC.grid.makeGrid(imSize, nodeSpacing=nodeSpacing) ####################################################### # Dilation of points for grain strain calc ####################################################### # Test with 1% dilation dilation = 1.01 # make sure average volumetric strain is the same as the projected one pointsDef = pointsRef * dilation displacements = (pointsDef - pointsRef).reshape(dims[0], dims[1], dims[2], 3) Ffield = spam.deformation.FfieldRegularGeers(displacements, nodeSpacing) ev = spam.deformation.decomposeFfield(Ffield, ["vol"])["vol"] # print("\tU trace mean new:",numpy.array(U)[:,[0,1,2],[0,1,2]].mean()) # print("\tev",numpy.array(ev)) # print("\tev",numpy.array(ev).mean()) self.assertAlmostEqual(ev.mean(), 3. * (dilation - 1), places=3) #self.assertAlmostEqual(ev.mean(), 3. * (dilation - 1), places=3) ######################################################## ## Z-stretch of points for strain calc ######################################################## zStretch = 1.02 print("\n\nSpreading points in Z by by 1.2") # Deform the points by spreading them apart -- this is dilating so positive volumetric strain displacements = (pointsRef * [zStretch, 1.0, 1.0] - pointsRef).reshape(dims[0], dims[1], dims[2], 3) Ffield = spam.deformation.FfieldRegularGeers(displacements, nodeSpacing) # only compute on internal nodes, we know there's a problem at the edges self.assertAlmostEqual(Ffield[1:-1,1:-1,1:-1,0,0].mean(), zStretch, places=3) self.assertAlmostEqual(Ffield[1:-1,1:-1,1:-1,1,1].mean(), 1.0, places=3) self.assertAlmostEqual(Ffield[1:-1,1:-1,1:-1,2,2].mean(), 1.0, places=3) ######################################################## ## Y-stretch of points for strain calc ######################################################## yStretch = 1.02 print("\n\nSpreading points in Y by by 1.2") # Deform the points by spreading them apart -- this is dilating so positive volumetric strain displacements = (pointsRef * [1.0, yStretch, 1.0] - pointsRef).reshape(dims[0], dims[1], dims[2], 3) Ffield = spam.deformation.FfieldRegularGeers(displacements, nodeSpacing) # only compute on internal nodes, we know there's a problem at the edges self.assertAlmostEqual(Ffield[1:-1,1:-1,1:-1,0,0].mean(), 1.0, places=3) self.assertAlmostEqual(Ffield[1:-1,1:-1,1:-1,1,1].mean(), yStretch, places=3) self.assertAlmostEqual(Ffield[1:-1,1:-1,1:-1,2,2].mean(), 1.0, places=3) ####################################################### # Rotation of points for grain strain calc ####################################################### pointsRotated = pointsRef.copy() rotAngle = 3.0 for n, point in enumerate(pointsRef): Phi = spam.deformation.computePhi({'r': [rotAngle, 0.0, 0.0]}, PhiCentre = numpy.array(((numpy.array(imSize)-1)/2)-1), PhiPoint = point) pointsRotated[n] += Phi[0:3, -1] # print( point, pointsRotated[n], '\n\n' ) displacements = (pointsRotated - pointsRef).reshape(dims[0], dims[1], dims[2], 3) Ffield = spam.deformation.FfieldRegularGeers(displacements, nodeSpacing) # Compute strains decomposedF = spam.deformation.decomposeFfield(Ffield, ["r", "vol", "dev"]) self.assertAlmostEqual(decomposedF['r'][1:-1,1:-1,1:-1, 0].mean(), rotAngle, places=3) self.assertAlmostEqual(decomposedF['r'][1:-1,1:-1,1:-1, 1].mean(), 0.0, places=3) self.assertAlmostEqual(decomposedF['r'][1:-1,1:-1,1:-1, 2].mean(), 0.0, places=3) self.assertAlmostEqual(decomposedF['vol'][1:-1,1:-1,1:-1].mean(), 0, places=3) self.assertAlmostEqual(decomposedF['dev'][1:-1,1:-1,1:-1].mean(), 0, places=3) ####################################################### # Shear of points for grain strain calc: # N.B.: this is a shear which results in a SYMMETRIC F ####################################################### pointsRotated = pointsRef.copy() shearVal = 0.1 for n, point in enumerate(pointsRef): Phi = spam.deformation.computePhi({'s': [shearVal, 0.0, 0.0]}, PhiCentre = numpy.array(((numpy.array(imSize)-1)/2)-1), PhiPoint = point) pointsRotated[n] += Phi[0:3, -1] # print( point, pointsRotated[n], '\n\n' ) displacements = (pointsRotated - pointsRef).reshape(dims[0], dims[1], dims[2], 3) Ffield = spam.deformation.FfieldRegularGeers(displacements, nodeSpacing) # Compute strains decomposedF = spam.deformation.decomposeFfield(Ffield, ["U"]) UmeanWithoutBoundaries = numpy.mean(decomposedF['U'][1:-1,1:-1,1:-1], axis=(0, 1, 2)) for i in range(3): for j in range(i,3): self.assertAlmostEqual(UmeanWithoutBoundaries[i,j], UmeanWithoutBoundaries[j,i], places=3) ####################################################### # Homogeneous shrinking together ####################################################### shrink = 0.99 # Deform the points by pushing them closer together -- this is compressing so negative volumetric strain displacements = (pointsRef * shrink - pointsRef).reshape(dims[0], dims[1], dims[2], 3) Ffield = spam.deformation.FfieldRegularGeers(displacements, nodeSpacing) decomposedF = spam.deformation.decomposeFfield(Ffield, ["vol"]) self.assertAlmostEqual(decomposedF['vol'][1:-1,1:-1,1:-1].mean(), 3. * (shrink - 1.), places=3) ### while we're at it, the brute force distance too, why not here Ffield = spam.deformation.FfieldRegularGeers(displacements, nodeSpacing, bruteForceDistance=True) decomposedF = spam.deformation.decomposeFfield(Ffield, ["vol"]) self.assertAlmostEqual(decomposedF['vol'][1:-1,1:-1,1:-1].mean(), 3. * (shrink - 1.), places=3) ####################################################### # Increasing neighbourhood radius should decrease noise ####################################################### #Add noise to displacements displacementAmplitude = displacements.max()-displacements.min() displacements += numpy.random.random(displacements.shape)*0.1*displacementAmplitude FfieldR1 = spam.deformation.FfieldRegularGeers(displacements, nodeSpacing, neighbourRadius=1) decomposedFR1 = spam.deformation.decomposeFfield(FfieldR1, ["vol"]) FfieldR1p5 = spam.deformation.FfieldRegularGeers(displacements, nodeSpacing, neighbourRadius=1.5) decomposedFR1p5 = spam.deformation.decomposeFfield(FfieldR1p5, ["vol"]) FfieldR2 = spam.deformation.FfieldRegularGeers(displacements, nodeSpacing, neighbourRadius=2) decomposedFR2 = spam.deformation.decomposeFfield(FfieldR2, ["vol"]) #FfieldR3 = spam.deformation.FfieldRegularGeers(displacements, nodeSpacing, neighbourRadius=3) #decomposedFR3 = spam.deformation.decomposeFfield(FfieldR3, ["vol"]) std1 = decomposedFR1[ 'vol'][1:-1,1:-1,1:-1].std() std1p5 = decomposedFR1p5['vol'][1:-1,1:-1,1:-1].std() std2 = decomposedFR2[ 'vol'][1:-1,1:-1,1:-1].std() #compute stdev of vol strain (most sensitive!!) and try for three different neighbourhoodRadii self.assertEqual(std1 > std1p5, True) self.assertEqual(std1p5 > std2, True) ####################################################### # Application of mask ####################################################### zStretch = 1.02 print("\n\nSpreading points in Z by by 1.2") # Deform the points by spreading them apart -- this is dilating so positive volumetric strain displacements = (pointsRef * [zStretch, 1.0, 1.0] - pointsRef).reshape(dims[0], dims[1], dims[2], 3) # manually add a NaN in the middle of the displacement field displacements[dims[0]//2, dims[1]//2, dims[2]//2] = [numpy.nan, numpy.nan, numpy.nan] Ffield = spam.deformation.FfieldRegularGeers(displacements, nodeSpacing) # only compute on internal nodes, we know there's a problem at the edges self.assertEqual(numpy.isnan(Ffield[:,:,:,0,0]).sum(), 1) 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 ####################################################### # New example with larger distance # pointsRef = numpy.array([[ 0, 0, 0], # [100, 0, 0], # [ 0, 100, 0], # [100, 100, 0], # [ 0, 0, 100], # [100, 0, 100], # [ 0, 100, 100], # [100, 100, 100]], dtype='= 6) self.assertEqual(connectivity.shape[1], 4) ####################################################### # Dilation of points for grain strain calc ####################################################### # Test with 1% dilation dilation = 1.01 # make sure average volumetric strain is the same as the projected one pointsDef = pointsRef * dilation tetVolumesRef = spam.mesh.tetVolumes(pointsRef, connectivity) tetVolumesDef = spam.mesh.tetVolumes(pointsDef, connectivity) # print("volStrains@tet:", (tetVolumesDef-tetVolumesRef)/tetVolumesRef) self.assertAlmostEqual(numpy.mean((tetVolumesDef - tetVolumesRef) / tetVolumesRef), 3. * (dilation - 1), places=3) #### # Check spam.deformation.FfieldBagi ### # Deform the points by spreading them apart -- this is dilating so positive volumetric strain # U, ev, ed = spam.deformation.FfieldBagi(points, connectivity, points - points) Ffield = spam.deformation.FfieldBagi(pointsRef, connectivity, pointsDef - pointsRef) ev = spam.deformation.decomposeFfield(Ffield, ["vol"])["vol"] # print("\tU trace mean new:",numpy.array(U)[:,[0,1,2],[0,1,2]].mean()) # print("\tev",numpy.array(ev)) # print("\tev",numpy.array(ev).mean()) self.assertAlmostEqual(ev.mean(), 3. * (dilation - 1), places=3) #print(spam.deformation.decomposeFfield(Ffield, ['vol'])) #self.assertAlmostEqual(ev.mean(), 3. * (dilation - 1), places=3) ######################################################## ## Z-stretch of points for grain strain calc ######################################################## zStretch = 1.02 print("\n\nSpreading points in Z by by 1.2") # Deform the points by spreading them apart -- this is dilating so positive volumetric strain Ffield = spam.deformation.FfieldBagi(pointsRef, connectivity, pointsRef * [zStretch, 1.0, 1.0] - pointsRef) self.assertAlmostEqual(Ffield[:,0,0].mean(), zStretch, places=3) # Make sure that a single NaN does not pollute strains and is safely ignored pointsRefOneBadPoint = pointsRef.copy() badPointCoord = pointsRef.shape[0]//2 pointsRefOneBadPoint[badPointCoord] = numpy.nan Ffield = spam.deformation.FfieldBagi(pointsRefOneBadPoint, connectivity, pointsRef * [zStretch, 1.0, 1.0] - pointsRef) # make sure that the nans in the strain are where there is our bad particle self.assertEqual(list(numpy.isnan(Ffield[:,0,0])), list(numpy.sum(connectivity==badPointCoord, axis=1)>0)) # check the mean strain without the bad ones self.assertAlmostEqual(numpy.nanmean(Ffield[:,0,0]), zStretch, places=3) # Same as above with all strain components (F and R) to pass them to grain projector Ffield = spam.deformation.FfieldBagi(pointsRef, connectivity, pointsRef * [zStretch, 1.0, 1.0] - pointsRef) self.assertAlmostEqual(Ffield[:,0,0].mean(), zStretch, places=3) # Project F back to grains grainF = spam.mesh.projectTetFieldToGrains(pointsRef * [zStretch, 1.0, 1.0], connectivity, Ffield) self.assertAlmostEqual(numpy.mean(grainF[:,0,0]), zStretch, places=3) self.assertAlmostEqual(numpy.mean(grainF[:,1,1]), 1.00, places=3) self.assertAlmostEqual(numpy.mean(grainF[:,2,2]), 1.00, places=3) ####################################################### # Rotation of points for grain strain calc ####################################################### pointsRotated = pointsRef.copy() rotAngle = 3.0 for n, point in enumerate(pointsRef): pointPad = numpy.ones(4) pointPad[0:3] = point # pointsRotated[n] = point + spam.DIC.computeTransformationOperator( {'r': [10.0,0.0,0.0]}, Phi = spam.deformation.computePhi({'r': [rotAngle, 0.0, 0.0]}, PhiCentre=[0.0, 50.0, 50.0], PhiPoint=point) # pointsRotated[n] = numpy.dot(F,pointPad)[0:3] pointsRotated[n] += Phi[0:3, -1] # print( point, pointsRotated[n], '\n\n' ) # Compute strains Ffield = spam.deformation.FfieldBagi(pointsRef, connectivity, pointsRotated - pointsRef) decomposedF = spam.deformation.decomposeFfield(Ffield, ["r", "vol", "dev"]) self.assertAlmostEqual(decomposedF['r'][:, 0].mean(), rotAngle, places=3) self.assertAlmostEqual(decomposedF['r'][:, 1].mean(), 0.0, places=3) self.assertAlmostEqual(decomposedF['r'][:, 2].mean(), 0.0, places=3) self.assertAlmostEqual(decomposedF['vol'].mean(), 0, places=3) self.assertAlmostEqual(decomposedF['dev'].mean(), 0, places=3) # Project to grains grain_r = spam.mesh.projectTetFieldToGrains(pointsRotated, connectivity, decomposedF['r']) self.assertAlmostEqual(numpy.mean(grain_r[:,0]), rotAngle, places=3) self.assertAlmostEqual(numpy.mean(grain_r[:,1]), 0.00, places=3) self.assertAlmostEqual(numpy.mean(grain_r[:,2]), 0.00, places=3) ####################################################### # Homogeneous shrinking together ####################################################### shrink = 0.99 # Deform the points by pushing them closer together -- this is compressing so negative volumetric strain Ffield = spam.deformation.FfieldBagi(pointsRef, connectivity, pointsRef * shrink - pointsRef) decomposedF = spam.deformation.decomposeFfield(Ffield, ["vol"]) self.assertAlmostEqual(decomposedF['vol'].mean(), 3. * (shrink - 1.), places=3) Ffield = spam.deformation.FfieldBagi(pointsRef, connectivity, pointsRef * shrink - pointsRef) # Project F back to grains grainF = spam.mesh.projectTetFieldToGrains(pointsRef * shrink, connectivity, Ffield) self.assertAlmostEqual(numpy.mean(grainF[:,0,0]), shrink, places=3) self.assertAlmostEqual(numpy.mean(grainF[:,1,1]), shrink, places=3) self.assertAlmostEqual(numpy.mean(grainF[:,2,2]), shrink, places=3) ### Attempt a grid projection #grid_bounds, num_in_grid, gridStrain = spam.mesh.projectBagiStrainToGrid(pointsRef * shrink, connectivity, F, nx=3, ny=3, nz=3) #print(grid_bounds) #print(num_in_grid) #print(gridStrain) ####################################################### # Move one point and see if the right elements strain ####################################################### # make sure that the connectivity is well respected in the output -- move only one point, # and make sure only the elements concerned are moved. # Move point two: pointToMove = 2 newPoints = pointsRef.copy() # newPoints[pointToMove] += numpy.random.randint( -10, 10, (3)) newPoints[pointToMove] += [10, 10, 10] # list of numbers tetsWithPointToMove = numpy.where(connectivity == pointToMove)[0] Ffield = spam.deformation.FfieldBagi(pointsRef, connectivity, newPoints - pointsRef) decomposedF = spam.deformation.decomposeFfield(Ffield, ["dev"]) # full length indexing bool array otherTets = numpy.ones(connectivity.shape[0], dtype=bool) otherTets[tetsWithPointToMove] = 0 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 ####################################################### #First we need to create some data using DEM dataset pixelSize = 0.0001 blurSTD = 0.8 noiseSTD = 0.01 boxSizeDEM, centres, radii = spam.datasets.loadUniformDEMboxsizeCentreRadius() # put 0 in the middle centres -= numpy.mean(centres, axis=0) rMax = numpy.amax(radii) # pad box size boxSizeDEM = boxSizeDEM + 5 * rMax # add half box size to centres centres += numpy.array(boxSizeDEM)/2.0 boxSizePx = (boxSizeDEM / pixelSize).astype(int) centresPx = centres / pixelSize radiiPx = radii / pixelSize box = numpy.zeros(boxSizePx, dtype=" 1.0)] = 1.0 box[numpy.where(box < 0.0)] = 0.0 box = box * 0.5 box = box + 0.25 box = scipy.ndimage.filters.gaussian_filter(box, sigma=blurSTD) box = numpy.random.normal(box, scale=noiseSTD) binIm0 = box >= 0.5 #Run watershed labIm0 = spam.label.ITKwatershed.watershed(binIm0) #Save images tifffile.imsave(testFolder + "Step0.tif", box.astype(' 1.0)] = 1.0 boxDeformed[numpy.where(boxDeformed < 0.0)] = 0.0 boxDeformed = boxDeformed * 0.5 boxDeformed = boxDeformed + 0.25 boxDeformed = scipy.ndimage.filters.gaussian_filter(boxDeformed, sigma=blurSTD) boxDeformed = numpy.random.normal(boxDeformed, scale=noiseSTD) #Save images tifffile.imsave(testFolder + "Step1.tif", boxDeformed.astype(' rsTwoLdic, True) # Check that there some mergeSource = 1 mergeSource = numpy.unique(numpy.genfromtxt(testFolder+'merged.tsv', names=True)['mergeSource']) self.assertEqual(mergeSource[0], 0) self.assertEqual(mergeSource[1], 1) def test_interpolatePhiField(self): # Fake Phi field -- this will be defined at 0,0,0 and +100 in all directions # z = 0 --> transformation below: transformation = {'t': [5., 0., 0.], 'r': [0., 0., 0.]} # z = 100 --> nothing fieldCoords = numpy.array([[ 0, 0, 0], [ 0, 0,100], [ 0,100, 0], [ 0,100,100], [100, 0, 0], [100, 0,100], [100,100, 0], [100,100,100]]) fieldValues = numpy.zeros([8,4,4]) fieldValues[0:4] = spam.deformation.computePhi(transformation) fieldValues[4:8] = numpy.eye(4) interpCoods = numpy.array([[50.,50.,50.]]) interpolatedPhi = spam.deformation.interpolatePhiField(fieldCoords, fieldValues, interpCoods) decomposedInterpolatedPhi = spam.deformation.decomposePhi(interpolatedPhi[0]) for key in transformation.keys(): for i in range(3): 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.]} PhiR1 = spam.deformation.computePhi(transformation) PhiR2 = spam.deformation.computeRigidPhi(PhiR1) for i in range(4): for j in range(4): self.assertAlmostEqual(PhiR1[i,j], PhiR2[i,j], places=2) transformation = {'t': [5., 0., 0.], 'r': [3., 0., 0.], 'z': [1.01, 1.01, 1.01]} PhiR1 = spam.deformation.computePhi(transformation) PhiR2 = spam.deformation.computeRigidPhi(PhiR1) for i in range(4): for j in range(4): if i==j and i < 3: 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])) res = spam.deformation.deformationFunction.decomposeF(F) self.assertIsNotNone(res) #Test for singular F F = numpy.array(([1, 2, 3],[2, 4, 6],[3, 6, 9])) res = spam.deformation.deformationFunction.decomposeF(F) self.assertTrue(numpy.isnan(res['r'][0])) #Test for nan F = numpy.array(([1, 2, 3],[2, 4, 6],[3, numpy.nan, 9])) res = spam.deformation.deformationFunction.decomposeF(F) self.assertTrue(numpy.isnan(res['r'][0])) #Test for inf F = numpy.array(([1, 2, 3],[2, 4, 6],[3, numpy.inf, 9])) res = spam.deformation.deformationFunction.decomposeF(F) self.assertTrue(numpy.isnan(res['r'][0])) def test_decomposePhiField(self): PhiField = numpy.array([numpy.eye(4) for point in range(3)]) decomposedPhiField = spam.deformation.decomposePhiField(PhiField, ["vol", "dev", "volss", "devss", "t", "r", "z", "U", "e"]) for point in range(PhiField.shape[0]): self.assertTrue(numpy.allclose(decomposedPhiField["t"][point], numpy.zeros(3))) self.assertTrue(numpy.allclose(decomposedPhiField["r"][point], numpy.zeros(3))) self.assertTrue(numpy.allclose(decomposedPhiField["U"][point], numpy.eye(3))) self.assertTrue(numpy.allclose(decomposedPhiField["vol"][point], 0.0)) self.assertTrue(numpy.allclose(decomposedPhiField["dev"][point], 0.0)) self.assertTrue(numpy.allclose(decomposedPhiField["e"][point], numpy.zeros((3,3)))) self.assertTrue(numpy.allclose(decomposedPhiField["volss"][point], 0.0)) self.assertTrue(numpy.allclose(decomposedPhiField["devss"][point], 0.0)) if __name__ == '__main__': unittest.main()