Skip to content

GitLab

Commit 03a90685 authored by Edward Andò's avatar Edward Andò
Browse files

decomposeFfield and decomposePhiField parallelised

parent ccb4a434
Hide whitespace changes
Inline Side-by-side
tools/deformation/deformationField.py
View file @ 03a90685
......@@ -17,11 +17,16 @@ this program. If not, see <http://www.gnu.org/licenses/>.
"""
# 2020-02-24 Olga Stamati and Edward Ando
from __future__ import print_function
import spam.deformation
import numpy
import progressbar
import spam.label
def FfieldRegularQ8(displacementField, nodeSpacing):
import multiprocessing
# Global number of processes
nProcessesDefault = multiprocessing.cpu_count()
def FfieldRegularQ8(displacementField, nodeSpacing, nProcesses=nProcessesDefault, verbose=True):
"""
This function computes the transformation gradient field F from a given displacement field.
Please note: the transformation gradient tensor: F = I + du/dx.
......@@ -37,6 +42,14 @@ def FfieldRegularQ8(displacementField, nodeSpacing):
nodeSpacing : 3-component list of floats
Length between two nodes in every direction (*i.e.,* size of a cell)
nProcesses : integer, optional
Number of processes for multiprocessing
Default = number of CPUs in the system
verbose : boolean, optional
Print progress?
Default = True
Returns
-------
F : (nz-1) x (ny-1) x (nx-1) x 3x3 array of n cells
......@@ -47,6 +60,8 @@ def FfieldRegularQ8(displacementField, nodeSpacing):
# Define dimensions
dims = list(displacementField.shape[0:3])
# Q8 has 1 element fewer than the number of displacement points
nCells = [n - 1 for n in dims]
# Check if a 2D field is passed
......@@ -89,10 +104,37 @@ def FfieldRegularQ8(displacementField, nodeSpacing):
jacY = 2.0 / float(nodeSpacing[1])
jacX = 2.0 / float(nodeSpacing[2])
bar = progressbar.ProgressBar(maxval=nCells[0]*nCells[1]*nCells[2]).start()
if verbose: bar = progressbar.ProgressBar(maxval=nCells[0]*nCells[1]*nCells[2]).start()
nodesDone = 0
# Loop over the cells
## Loop over the cells
#global computeConvexVolume
#def computeConvexVolume(label):
#labelI = spam.label.getLabel(lab, label, boundingBoxes = boundingBoxes, centresOfMass = centresOfMass)
#subvol = labelI['subvol']
#points = numpy.transpose(numpy.where(subvol))
#try:
#hull = scipy.spatial.ConvexHull(points)
#deln = scipy.spatial.Delaunay(points[hull.vertices])
#idx = numpy.stack(numpy.indices(subvol.shape), axis = -1)
#out_idx = numpy.nonzero(deln.find_simplex(idx) + 1)
#hullIm = numpy.zeros(subvol.shape)
#hullIm[out_idx] = 1
#hullVol = spam.label.volumes(hullIm)
#return label, hullVol[-1]
#except:
#return label, 0
## Run multiprocessing
#with multiprocessing.Pool(processes=nProcesses) as pool:
#for returns in pool.imap_unordered(computeConvexVolume, range(1, nLabels+1)):
#if verbose:
#finishedNodes += 1
#widgets[0] = progressbar.FormatLabel(" vol={:0>7.5f} ".format(returns[1]))
#pbar.update(finishedNodes)
#convexVolume[returns[0]] = returns[1]
for kCell in range(nCells[0]):
for jCell in range(nCells[1]):
for iCell in range(nCells[2]):
......@@ -172,7 +214,6 @@ def FfieldRegularGeers(displacementField, nodeSpacing, neighbourRadius=1.5, mask
Taken from the implementation in "TomoWarp2: A local digital volume correlation code", Tudisco et al., 2017
"""
import numpy
import progressbar
import scipy.spatial
pbar = progressbar.ProgressBar()
......@@ -397,7 +438,7 @@ def FfieldBagi(points, connectivity, displacements):
return Ffield
def decomposeFfield(Ffield, components, twoD=False):
def decomposeFfield(Ffield, components, twoD=False, nProcesses=nProcessesDefault, verbose=True):
"""
This function takes in an F field (from either FfieldRegularQ8, FfieldRegularGeers, FfieldBagi) and
returns fields of desired transformation components.
......@@ -414,6 +455,14 @@ def decomposeFfield(Ffield, components, twoD=False):
Is the Ffield in 2D? This changes the strain calculation.
Default = False
nProcesses : integer, optional
Number of processes for multiprocessing
Default = number of CPUs in the system
verbose : boolean, optional
Print progress?
Default = True
Returns
-------
Dictionary containing appropriately reshaped fields of the transformation components requested.
......@@ -423,9 +472,6 @@ def decomposeFfield(Ffield, components, twoD=False):
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 3x3 F field
fieldDimensions = Ffield.shape[0:-2]
fieldRavelLength = numpy.prod(numpy.array(fieldDimensions))
......@@ -440,12 +486,38 @@ def decomposeFfield(Ffield, components, twoD=False):
if component == 'U' or component == 'e':
output[component] = numpy.zeros((fieldRavelLength, 3, 3))
# Iterate through flat field of Fs
for n in pbar(range(FfieldFlat.shape[0])):
# Function for parallel mode
global decomposeOneF
def decomposeOneF(n):
F = FfieldFlat[n]
decomposedF = spam.deformation.decomposeF(F, twoD=twoD)
for component in components:
output[component][n] = decomposedF[component]
if numpy.isfinite(F).sum() == 9:
decomposedF = spam.deformation.decomposeF(F, twoD=twoD)
return n, decomposedF
else:
return n, {'r': numpy.array([numpy.nan] * 3),
'z': numpy.array([numpy.nan] * 3),
'U': numpy.eye(3) * numpy.nan,
'e': numpy.eye(3) * numpy.nan,
'vol': numpy.nan,
'dev': numpy.nan,
'volss': numpy.nan,
'devss': numpy.nan}
# Iterate through flat field of Fs
if verbose: pbar = progressbar.ProgressBar(maxval=fieldRavelLength).start()
finishedPoints = 0
# Run multiprocessing
with multiprocessing.Pool(processes=nProcesses) as pool:
for returns in pool.imap_unordered(decomposeOneF, range(fieldRavelLength)):
if verbose:
finishedPoints += 1
pbar.update(finishedPoints)
for component in components:
output[component][returns[0]] = returns[1][component]
if verbose: pbar.finish()
# Reshape on the output
for component in components:
......@@ -461,7 +533,7 @@ def decomposeFfield(Ffield, components, twoD=False):
return output
def decomposePhiField(PhiField, components, twoD=False):
def decomposePhiField(PhiField, components, twoD=False, nProcesses=nProcessesDefault, verbose=True):
"""
This function takes in a Phi field (from readCorrelationTSV?) and
returns fields of desired transformation components.
......@@ -478,6 +550,14 @@ def decomposePhiField(PhiField, components, twoD=False):
Is the PhiField in 2D? This changes the strain calculation.
Default = False
nProcesses : integer, optional
Number of processes for multiprocessing
Default = number of CPUs in the system
verbose : boolean, optional
Print progress?
Default = True
Returns
-------
Dictionary containing appropriately reshaped fields of the transformation components requested.
......@@ -487,10 +567,6 @@ def decomposePhiField(PhiField, components, twoD=False):
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))
......@@ -505,12 +581,39 @@ def decomposePhiField(PhiField, components, twoD=False):
if component == 'U' or component == 'e':
output[component] = numpy.zeros((fieldRavelLength, 3, 3))
# Iterate through flat field of Phis
for n in pbar(range(PhiFieldFlat.shape[0])):
# Function for parallel mode
global decomposeOnePhi
def decomposeOnePhi(n):
Phi = PhiFieldFlat[n]
decomposedPhi = spam.deformation.decomposePhi(Phi, twoD=twoD)
for component in components:
output[component][n] = decomposedPhi[component]
if numpy.isfinite(Phi).sum() == 16:
decomposedPhi = spam.deformation.decomposePhi(Phi, twoD=twoD)
return n, decomposedPhi
else:
return n, {'t': numpy.array([numpy.nan] * 3),
'r': numpy.array([numpy.nan] * 3),
'z': numpy.array([numpy.nan] * 3),
'U': numpy.eye(3) * numpy.nan,
'e': numpy.eye(3) * numpy.nan,
'vol': numpy.nan,
'dev': numpy.nan,
'volss': numpy.nan,
'devss': numpy.nan}
# Iterate through flat field of Fs
if verbose: pbar = progressbar.ProgressBar(maxval=fieldRavelLength).start()
finishedPoints = 0
# Run multiprocessing
with multiprocessing.Pool(processes=nProcesses) as pool:
for returns in pool.imap_unordered(decomposeOnePhi, range(fieldRavelLength)):
if verbose:
finishedPoints += 1
pbar.update(finishedPoints)
for component in components:
output[component][returns[0]] = returns[1][component]
if verbose: pbar.finish()
# Reshape on the output
for component in components:
......@@ -1042,30 +1145,30 @@ def interpolatePhiField(fieldCoords, fieldValues, interpCoords, fieldInterpBinRa
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()` with `readConvergence=True, readPSCC=True, readLabelDilate=True`
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.
......@@ -1077,7 +1180,7 @@ def getDisplacementFromNeighbours(labIm, DVC, fileName, method = 'getLabel', cen
Dictionary
TSV file with the same columns as the input
"""
# Compute centreOfMass if needed
if centresOfMass == None:
centresOfMass = spam.label.centresOfMass(labIm)
......@@ -1093,7 +1196,7 @@ def getDisplacementFromNeighbours(labIm, DVC, fileName, method = 'getLabel', cen
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
......@@ -1173,34 +1276,34 @@ def getDisplacementFromNeighbours(labIm, DVC, fileName, method = 'getLabel', cen
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')
def mergeRegistrationAndDiscreteFields(regTSV, discreteTSV, fileName, regTSVBinRatio = 1):
"""
This function merges a registration TSV with a discrete TSV.
Can be used to create the first guess for `spam-ddic`, using the registration over the whole file, and a previous result from `spam-ddic`.
Parameters
-----------
regTSV : dictionary
Dictionary with deformation field, obtained from a registrion, usually from the whole sample, and read using `spam.helpers.tsvio.readCorrelationTSV()`
discreteTSV : 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
regTSVBinRatio : float, optional
Change translations from regTSV, if it's been calculated on a differently-binned image. Default = 1
Returns
--------
Dictionary
TSV file with the same columns as the input
"""
# Create a first guess
phiGuess = discreteTSV['PhiField'].copy()
# Main loop
......@@ -1215,7 +1318,7 @@ def mergeRegistrationAndDiscreteFields(regTSV, discreteTSV, fileName, regTSVBinR
PhiPoint = deformPos)['t']
# Add the extra disp to the phi guess
phiGuess[lab][:-1,-1] += extraDisp*regTSVBinRatio
# Save
outMatrix = numpy.array([numpy.array(range(discreteTSV['numberOfLabels'])),
discreteTSV['fieldCoords'][:, 0],
......
tools/deformation/deformationFunction.py
View file @ 03a90685
......@@ -193,7 +193,7 @@ def decomposeF(F, twoD=False):
if numpy.isinf(F).sum() > 0:
print("deformationFunction.decomposeF(): Inf value in F. Exiting")
return transformation
# Check for singular F if yes quit
try:
numpy.linalg.inv(F)
......
tools/tests/test_deformation.py
View file @ 03a90685
......@@ -35,7 +35,7 @@ class testAll(unittest.TestCase):
os.remove(testFolder+'merged.tsv')
os.remove(testFolder+'TSV_getDisplacementFromNeighbours.tsv')
os.remove(testFolder+'TSV_mergeRegistrationAndDiscreteFields.tsv')
except OSError:
pass
......@@ -953,7 +953,7 @@ 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
......@@ -1038,8 +1038,8 @@ class testAll(unittest.TestCase):
"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',
......@@ -1049,7 +1049,7 @@ class testAll(unittest.TestCase):
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
......@@ -1066,7 +1066,7 @@ class testAll(unittest.TestCase):
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
......@@ -1083,7 +1083,7 @@ class testAll(unittest.TestCase):
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())
# Case 3: Run the function with an incomplete TSV file
TSV = spam.helpers.readCorrelationTSV('Step0-Step1-discreteDVC.tsv',
readConvergence=False,
......@@ -1092,11 +1092,11 @@ class testAll(unittest.TestCase):
readPSCC=False)
res = spam.deformation.deformationField.getDisplacementFromNeighbours(imLab, TSV, './TSV_getDisplacementFromNeighbours.tsv', previousDVC = TSV)
self.assertIsNone(res)
def test_mergeRegistrationAndDiscreteFields(self):
#Create an initial Lab
#######################################################
### We're using the DDIC test from scripts here, lightly modified
#######################################################
......@@ -1132,9 +1132,9 @@ class testAll(unittest.TestCase):
#Save images
tifffile.imsave(testFolder + "Step0.tif", box.astype('<f4'))
tifffile.imsave(testFolder + "Lab0.tif", labIm0.astype(spam.label.labelType))
# Create the first step
#Create Phi and Apply (25 px displacement on Y-axis, and 0 degree rotation along Z axis)
translationStep1 = [5, 2, 0]
rotationStep1 = [0, 0, 0]
......@@ -1156,9 +1156,9 @@ class testAll(unittest.TestCase):
boxDeformed = numpy.random.normal(boxDeformed, scale=noiseSTD)
#Save images
tifffile.imsave(testFolder + "Step1.tif", boxDeformed.astype('<f4'))
# Create the second step
#Create Phi and Apply (25 px displacement on Y-axis, and 5 degree rotation along Z axis)
translationStep2 = [10, 4, 0]
rotationStep2 = [0, 0, 0]
......@@ -1180,7 +1180,7 @@ class testAll(unittest.TestCase):
boxDeformed = numpy.random.normal(boxDeformed, scale=noiseSTD)
#Save images
tifffile.imsave(testFolder + "Step2.tif", boxDeformed.astype('<f4'))
# Create the TSV from macro reg between 0 and 1
exitCode = subprocess.call(["spam-ldic",
"-glt", "0.5",
......@@ -1192,16 +1192,16 @@ class testAll(unittest.TestCase):
"Step0.tif", "Step1.tif",
"-od", testFolder+""])
self.assertEqual(exitCode, 0)
# Create the TSV from ddic between 0 and 1
#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)
# Create the TSV from macro reg between 1 and 2
exitCode = subprocess.call(["spam-ldic",
"-glt", "0.5",
......@@ -1213,14 +1213,14 @@ class testAll(unittest.TestCase):
"Step1.tif", "Step2.tif",
"-od", testFolder+""])
self.assertEqual(exitCode, 0)
# Read ddic TSV
discreteTSV = spam.helpers.readCorrelationTSV('Step0-Step1-discreteDVC.tsv',
readConvergence=True,
readError=True,
readLabelDilate=True,
readPSCC=True)
# Read macro TSV
macroTSV = spam.helpers.readCorrelationTSV('Step1-Step2-bin2-registration.tsv',
fieldBinRatio=2)
......
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