Commit e9238d3e authored by Gustavo Pinzon's avatar Gustavo Pinzon
Browse files

New function spam.deformation.deformationField.getDisplacementFromNeighbours() with test

parent d30c6577
Pipeline #52843 passed with stages
in 13 minutes and 15 seconds
......@@ -19,7 +19,8 @@ this program. If not, see <http://www.gnu.org/licenses/>.
# 2020-02-24 Olga Stamati and Edward Ando
from __future__ import print_function
import numpy
import progressbar
import spam.label
def FfieldRegularQ8(displacementField, nodeSpacing):
"""
This function computes the transformation gradient field F from a given displacement field.
......@@ -41,8 +42,6 @@ def FfieldRegularQ8(displacementField, nodeSpacing):
F : (nz-1) x (ny-1) x (nx-1) x 3x3 array of n cells
The field of the transformation gradient tensors
"""
import numpy
import progressbar
# import spam.DIC.deformationFunction
#import spam.mesh.strain
......@@ -1020,3 +1019,143 @@ def interpolatePhiField(fieldCoords, fieldValues, interpCoords, fieldInterpBinRa
output[:, 0:3, 3] *= fieldInterpBinRatio
return output
def getDisplacementFromNeighbours(labIm, DVC, fileName, method = 'getLabel', centresOfMass = None, previousDVC = None):
"""
This function computes the displacement as the mean displacement from the neighbours, for non-converged grains using a TSV file obtained from `spam-ddic` script. Returns a new TSV file with the new Phi (composed only by the displacement part).
The generated TSV can be used as an input for `spam-ddic`.
Parameters
-----------
lab : 3D array of integers
Labelled volume, with lab.max() labels
DVC : dictionary
Dictionary with deformation field, obtained from `spam-ddic` script, and read using `spam.helpers.tsvio.readCorrelationTSV()`
fileName : string
FileName including full path and .tsv at the end to write
method : string, optional
Method to compute the neighbours using `spam.label.getNeighbours()`.
'getLabel' : The neighbours are the labels inside the subset obtained through spam.getLabel()
'mesh' : The neighbours are computed using a tetrahedral connectivity matrix
Default = 'getLabel'
centresOfMass : lab.max()x3 array of floats, optional
Centres of mass in format returned by ``centresOfMass``.
If not defined (Default = None), it is recomputed by running ``centresOfMass``
previousDVC : dictionary, optional
Dictionary with deformation field, obtained from `spam-ddic` script, and read using `spam.helpers.tsvio.readCorrelationTSV()` for the previous step.
This allows the to compute only the displacement increment from the neighbours, while using the F tensor from a previous (converged) step.
If `previousDVS = None`, then the resulting Phi would be composed only by the displacement of the neighbours.
Default = None
fileName : string
Output filename, if None return dictionary as from spam.helpers.readCorrelationTSV()
Default = None
Returns
--------
Dictionary
Output dictionary with the same columns as the input
"""
# Compute centreOfMass if needed
if centresOfMass == None:
centresOfMass = spam.label.centresOfMass(labIm)
# Get number of labels
numberOfLabels = (labIm.max() + 1).astype('u4')
# Create Phi field
PhiField = numpy.zeros((numberOfLabels, 4, 4), dtype='<f4')
# Rest of arrays
try:
iterations = DVC['iterations']
returnStatus = DVC['returnStatus']
deltaPhiNorm = DVC['deltaPhiNorm']
PSCC = DVC['PSCC']
labelDilateList = DVC['LabelDilate']
error = DVC['error']
# Get the problematic labels
probLab = numpy.where(DVC['returnStatus']!= 2)[0]
# Remove the 0 from the wrongLab list
probLab = probLab[~numpy.in1d(probLab, 0)]
# Get neighbours
neighbours = spam.label.getNeighbours(labIm, probLab, method = method)
# Solve first the converged particles - make a copy of the PhiField
for i in range(numberOfLabels):
PhiField[i] = DVC['PhiField'][i]
# Work on the problematic labels
widgets = [progressbar.FormatLabel(''), ' ', progressbar.Bar(), ' ', progressbar.AdaptiveETA()]
pbar = progressbar.ProgressBar(widgets=widgets, maxval=len(probLab))
pbar.start()
for i in range(0, len(probLab), 1):
wrongLab = probLab[i]
neighboursLabel = neighbours[i]
t = []
for n in neighboursLabel:
# Check the return status of each neighbour
if DVC['returnStatus'][n] == 2:
# We dont have a previous DVC file loaded
if previousDVC is None:
# Get translation, rotation and zoom from Phi
nPhi = spam.deformation.deformationFunction.decomposePhi(DVC['PhiField'][n])
# Append the results
t.append(nPhi['t'])
# We have a previous DVC file loaded
else:
# Get translation, rotation and zoom from Phi at t=step
nPhi = spam.deformation.deformationFunction.decomposePhi(DVC['PhiField'][n])
# Get translation, rotation and zoom from Phi at t=step-1
nPhiP = spam.deformation.deformationFunction.decomposePhi(previousDVC['PhiField'][n])
# Append the incremental results
t.append(nPhi['t']-nPhiP['t'])
# Transform list to array
if not t:
# This is a non-working label, take care of it
Phi = spam.deformation.computePhi({'t': [0,0,0]})
PhiField[wrongLab] = Phi
else:
t = numpy.asarray(t)
# Compute mean
meanT = numpy.mean(t, axis = 0)
# Reconstruct
transformation = {'t': meanT}
Phi = spam.deformation.computePhi(transformation)
# Save
if previousDVC is None:
PhiField[wrongLab] = Phi
else:
PhiField[wrongLab] = previousDVC['PhiField'][wrongLab]
# Add the incremental displacement
PhiField[wrongLab][:-1,-1] += Phi[:-1,-1]
# Save
outMatrix = numpy.array([numpy.array(range(numberOfLabels)),
centresOfMass[:, 0], centresOfMass[:, 1], centresOfMass[:, 2],
PhiField[:, 0, 3], PhiField[:, 1, 3], PhiField[:, 2, 3],
PhiField[:, 0, 0], PhiField[:, 0, 1], PhiField[:, 0, 2],
PhiField[:, 1, 0], PhiField[:, 1, 1], PhiField[:, 1, 2],
PhiField[:, 2, 0], PhiField[:, 2, 1], PhiField[:, 2, 2],
PSCC, error, iterations, returnStatus, deltaPhiNorm, labelDilateList]).T
numpy.savetxt(fileName,
outMatrix,
fmt='%.7f',
delimiter='\t',
newline='\n',
comments='',
header="Label\tZpos\tYpos\tXpos\t" +
"Zdisp\tYdisp\tXdisp\t" +
"Fzz\tFzy\tFzx\t" +
"Fyz\tFyy\tFyx\t" +
"Fxz\tFxy\tFxx\t" +
"PSCC\terror\titerations\treturnStatus\tdeltaPhiNorm\tLabelDilate")
except:
print('spam.deformation.deformationField.getDisplacementFromNeighbours(): Missing information in the input TSV. Make sure you are reading iterations, returnStatus, deltaPhiNorm, PSCC, LabelDilate, and error.')
print('spam.deformation.deformationField.getDisplacementFromNeighbours(): Aborting')
......@@ -30,6 +30,7 @@ class testAll(unittest.TestCase):
os.remove(testFolder+'Step0-Step1-discreteDVC.tsv')
os.remove(testFolder+'Step0-Step1-discreteDVC.vtk')
os.remove(testFolder+'merged.tsv')
os.remove(testFolder+'TSV_getDisplacementFromNeighbours.tsv')
except OSError:
pass
......@@ -909,6 +910,136 @@ class testAll(unittest.TestCase):
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))
def test_getDisplacementFromNeighbours(self):
# Create a TSV and Lab file just as in test_merge
#######################################################
### 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="<f8")
spam.kalisphera.makeSphere(box, centresPx, radiiPx)
box[numpy.where(box > 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('<f4'))
tifffile.imsave(testFolder + "Lab0.tif", labIm0.astype(spam.label.labelType))
#test of rigid translation and rotation
#Create Phi and Apply (25 px displacement on Y-axis, and 5 degree rotation along Z axis)
translationStep1 = [5, 2, 0]
rotationStep1 = [0, 0, 0]
transformation = {'t': translationStep1,
'r': rotationStep1}
Phi = spam.deformation.computePhi(transformation)
# transform centres around the centres of the box
centresPxDeformed = numpy.zeros_like(centresPx)
for i, centrePx in enumerate(centresPx):
centresPxDeformed[i] = centrePx + spam.deformation.decomposePhi(Phi, PhiPoint=centrePx, PhiCentre=numpy.array(boxSizePx)/2.0)['t']
boxDeformed = numpy.zeros(boxSizePx, dtype="<f8")
spam.kalisphera.makeSphere(boxDeformed, centresPxDeformed, radiiPx)
boxDeformed[numpy.where(boxDeformed > 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('<f4'))
#######################################################
### Now use the ddic and ldic
#######################################################
exitCode = subprocess.call(["spam-ldic",
"-glt", "0.5",
"-hws", "10",
"-ns", "10",
"-reg", "-regbb", "2", "-regbe", "2",
"-tsv",
"-gpi", "5",
"Step0.tif", "Step1.tif",
"-od", testFolder+""])
self.assertEqual(exitCode, 0)
exitCode = subprocess.call(["spam-ddic",
"-pf", "Step0-Step1-bin2-registration.tsv", "-pfb", "2",
"-ld", "2",
"Step0.tif", "Lab0.tif", "Step1.tif",
"-od", testFolder+""])
self.assertEqual(exitCode, 0)
# Read lab file and TSV from previous code
imLab = tifffile.imread('Lab0.tif')
TSV = spam.helpers.readCorrelationTSV('Step0-Step1-discreteDVC.tsv',
readConvergence=True,
readError=True,
readLabelDilate=True,
readPSCC=True)
# Manually change the RS of one grain
TSV['returnStatus'][10] = -1
# Case 1: Run the function normally and check the displacement and F
exitCode = spam.deformation.deformationField.getDisplacementFromNeighbours(imLab, TSV, './TSV_getDisplacementFromNeighbours.tsv')
# Read the resulting TSV
tsvRes = spam.helpers.readCorrelationTSV('TSV_getDisplacementFromNeighbours.tsv',
readConvergence=True,
readError=True,
readLabelDilate=True,
readPSCC=True)
# Get the mean translation
meanT = spam.deformation.deformationFunction.decomposePhi(tsvRes['PhiField'][10])['t']
# Check the values (The values are linked to the displacement imposed at the beginning of the code on test_merge)
self.assertAlmostEqual(meanT[0], 5, places=2)
self.assertAlmostEqual(meanT[1], 2, places=2)
self.assertAlmostEqual(meanT[2], 0, places=2)
# Check that the F matrix is equal to numpy.eye(3)
self.assertTrue((tsvRes['PhiField'][10][:-1, :-1] == numpy.eye(3)).all())
# Case 2: Run the function using the same TSV as the previous step to preserve F
spam.deformation.deformationField.getDisplacementFromNeighbours(imLab, TSV, './TSV_getDisplacementFromNeighbours.tsv', previousDVC = TSV)
# Read the resulting TSV
tsvRes = spam.helpers.readCorrelationTSV('TSV_getDisplacementFromNeighbours.tsv',
readConvergence=True,
readError=True,
readLabelDilate=True,
readPSCC=True)
# Get the mean translation
meanT = spam.deformation.deformationFunction.decomposePhi(tsvRes['PhiField'][10])['t']
# Check the values (The values are linked to the displacement imposed at the beginning of the code on test_merge)
self.assertAlmostEqual(meanT[0], 5, places=2)
self.assertAlmostEqual(meanT[1], 2, places=2)
self.assertAlmostEqual(meanT[2], 0, places=2)
# Check that the F matrix is equal to the initial
self.assertTrue((tsvRes['PhiField'][10][:-1, :-1] == TSV['PhiField'][10][:-1, :-1]).all())
if __name__ == '__main__':
unittest.main()
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment