Commit 17206f39 authored by Edward Andò's avatar Edward Andò
Browse files

[skip-ci] final changes to splitScripts

parent 57244cd8
......@@ -4,13 +4,6 @@
Scripts in spam
*******************
.. graphviz:: graphs/DIC.dot
|
|
|
|
|
**Spam** provides scripts that tie together a number of the functions available in python.
These are expected to be called from the command line (after activating your virtual environment):
......@@ -37,44 +30,51 @@ All scripts will report all their options with a short explanation if they are s
This can be done by opening it as a stack in ImageJ/Fiji and saving directly,
or loading each slice into a 3D numpy array and saving it with `tifffile.imsave`
Here's a list of existing **spam** scripts organised into categories:
- Registration:
- **spam-ereg** — Eye-alignment tool for pre-registration alignment
- **spam-reg** — Registration between two images, measuring a **single** Φ (homogeneous deformation function) that maps im1 into im2, see: :ref:`reg`
- **spam-mmr** — Multi-modal registration for images acquired in different modalities
- **spam-mmr-graphical** — Graphical version of the above
- Regular grid correlation:
- **spam-ldic** — "Standard" local image correlation with independent subvolumes distributed on a regular grid, see: :ref:`ldic`
- **spam-regularStrain** — Takes in displacements from a TSV file from above and computes strains :ref:`regularStrainScript`
.. Here's a list of existing **spam** scripts organised into categories:
..
.. - Registration:
..
.. - **spam-ereg** — Eye-alignment tool for pre-registration alignment
.. - **spam-reg** — Registration between two images, measuring a **single** Φ (homogeneous deformation function) that maps im1 into im2, see: :ref:`reg`
..
.. - **spam-mmr** — Multi-modal registration for images acquired in different modalities
.. - **spam-mmr-graphical** — Graphical version of the above
..
.. - Regular grid correlation:
..
.. - **spam-ldic** — "Standard" local image correlation with independent subvolumes distributed on a regular grid, see: :ref:`ldic`
.. - **spam-regularStrain** — Takes in displacements from a TSV file from above and computes strains :ref:`regularStrainScript`
..
.. - Global correlation:
..
.. - **spam-gdic** — Global image correlation in a tetrahedral mesh, see: :ref:`globalDICscript`
.. - **spam-deformImageFromField** — Deforms an image from a registration TSV file or a displacement field interpolated with a tetrahedral mesh
..
.. - Labelled/discrete operations:
..
.. - **spam-ITKwatershed** — Morphological watershed from ITK, requires binary image as input
.. - **spam-ddic** — Discrete image correlation with discrete objects defined with labelled image see: :ref:`discreteDICscript`
.. - **spam-moveGrains** — Deform a labelled image with displacements defined for each label
.. - **spam-discreteStrain** — Computes strains on a tetrahedral mesh linking grain centres
.. - **spam-ereg-discrete** — A graphical tool to pre-align grains for correlation
- Global correlation:
|
|
- **spam-gdic** — Global image correlation in a tetrahedral mesh, see: :ref:`globalDICscript`
- **spam-deformImageFromField** — Deforms an image from a registration TSV file or a displacement field interpolated with a tetrahedral mesh
Description of image correlation scripts
#########################################
|
- Labelled/discrete operations:
Here's the general flowchat for how we recommend to use the different image correlation scripts in spam:
- **spam-ITKwatershed** — Morphological watershed from ITK, requires binary image as input
- **spam-ddic** — Discrete image correlation with discrete objects defined with labelled image see: :ref:`discreteDICscript`
- **spam-moveGrains** — Deform a labelled image with displacements defined for each label
- **spam-discreteStrain** — Computes strains on a tetrahedral mesh linking grain centres
- **spam-ereg-discrete** — A graphical tool to pre-align grains for correlation
.. graphviz:: graphs/DIC.dot
|
|
|
|
|
Description of image correlation scripts
#########################################
|
|
.. _ereg:
Eye-registration script (spam-ereg)
......
......@@ -4,32 +4,32 @@
Tutorial: Image correlation -- Practice
****************************************
This is a general tutorial for doing image correlation: if your objective is to measure a displacement field with regularly-spacing points you're in the right place!
If you have discrete grains and you want to do particle tracking, read the "alignment" part of the tutorial, then go here: :ref:`discreteImageCorrelationTutorial`
Binning with illustration of how far to go
REMEMBER FOR register() to work you need to get your points close to the right transformation to be able to optimise
< flow chart >
When do you need alignment by eye? When it looks like this -> ROB
When do you need registration? When you have a clear and smoothly varing deformation like this -> NOT A SHEAR BAND!!
---------
OK let's define how to do the local part (copy from scripts page):
When do you need pixel search? When the initial guess from above is too rough introduce GLT GHT (copy from scripts page)
LDIC tricks, return status, etc...
********************
OLD VERSION
********************
.. This is a general tutorial for doing image correlation: if your objective is to measure a displacement field with regularly-spacing points you're in the right place!
..
.. If you have discrete grains and you want to do particle tracking, read the "alignment" part of the tutorial, then go here: :ref:`discreteImageCorrelationTutorial`
..
.. Binning with illustration of how far to go
..
.. REMEMBER FOR register() to work you need to get your points close to the right transformation to be able to optimise
..
.. < flow chart >
..
.. When do you need alignment by eye? When it looks like this -> ROB
..
.. When do you need registration? When you have a clear and smoothly varing deformation like this -> NOT A SHEAR BAND!!
..
.. ---------
..
.. OK let's define how to do the local part (copy from scripts page):
..
.. When do you need pixel search? When the initial guess from above is too rough introduce GLT GHT (copy from scripts page)
..
.. LDIC tricks, return status, etc...
..
..
.. ********************
.. OLD VERSION
.. ********************
Manipulating Φ
......
......@@ -264,13 +264,13 @@ else:
# TODO: Last case with DDIC in and DDIC out could be with NNEIGHBOURS
outputPhiField = spam.DIC.interpolatePhiField(goodInputNodePositions,
goodInputPhiField,
outputNodePositions,
nNeighbours=args.NUMBER_OF_NEIGHBOURS,
neighbourRadius=args.NEIGHBOUR_RADIUS,
mode=args.MODE,
nProcesses=args.PROCESSES,
verbose=True)
goodInputPhiField,
outputNodePositions,
nNeighbours=args.NUMBER_OF_NEIGHBOURS,
neighbourRadius=args.NEIGHBOUR_RADIUS,
mode=args.MODE,
nProcesses=args.PROCESSES,
verbose=True)
#- apply a registration (single Phi) to a series of points:
#- defined on a grid with NS
......
......@@ -113,7 +113,8 @@ if regReturns['returnStatus'] == 1:
print("\n\nRegistration hit max iterations, OK... saving")
if regReturns['returnStatus'] > 0 and args.DEF:
tifffile.imsave(args.OUT_DIR+"/"+os.path.splitext(os.path.basename(args.im1.name))[0]+'-reg-def.tif', spam.DIC.applyPhi(im1, Phi=regReturns['Phi']))
tifffile.imsave(args.OUT_DIR+"/"+os.path.splitext(os.path.basename(args.im1.name))[0]+'-reg-def.tif',
spam.DIC.applyPhi(im1, Phi=regReturns['Phi']).astype(im1.dtype))
if regReturns['returnStatus'] < 0:
print("\n\nWe're saving this registration but we don't trust it at all")
......
......@@ -72,44 +72,19 @@ if dims[0] == 1:
else:
twoD = False
dispFlat = f["PhiField"][:,:3,-1]
# Check if a mask of points (based on return status of the correlation) is asked
if args.MASK:
mask = numpy.zeros(fieldCoords[:,0].shape).reshape(dims)
ignoreBackGround = True
else:
ignoreBackGround = False
# Check if the correction of the input field is asked
if args.CORRECT_FIELD or args.CORRECT_MEDIAN_FILTER:
print("\nspam-regularStrain: Correcting/filtering field...")
fieldValues = spam.DIC.correctPhiField(fileName = args.inFile.name,
correctBadPoints = args.CORRECT_FIELD,
ignoreBackGround = ignoreBackGround,
deltaPhiNormMin = args.CORRECT_DELTA_PHI_NORM,
pixelSearchCCmin = args.CORRECT_PIXEL_SEARCH_CC,
nNeighbours = args.CORRECT_NEIGHBOURS,
filterPoints = args.CORRECT_MEDIAN_FILTER,
filterPointsRadius = args.CORRECT_MEDIAN_FILTER_RADIUS)
# Extract the corrected displacement field and reshape it
disp = fieldValues[:,:3,-1].reshape(dims[0], dims[1], dims[2], 3)
else:
# Directly extract the displacement field without correction/filtering
disp = f["PhiField"][:,:3,-1].reshape(dims[0], dims[1], dims[2], 3)
# Mask background if asked
if ignoreBackGround:
mask[numpy.where(f["returnStatus"].reshape(dims)<-4)] = numpy.nan
disp[:,:,:,0] += mask
disp[:,:,:,1] += mask
disp[:,:,:,2] += mask
mask = f["returnStatus"] < -4
dispFlat[mask] = numpy.nan
disp = dispFlat.reshape(dims[0], dims[1], dims[2], 3)
print("\nspam-regularStrain: Computing F=I+du/dx")
if not args.Q8:
Ffield = spam.deformation.FfieldRegularGeers(disp,
nodeSpacing=nodeSpacing,
mask=args.MASK,
neighbourRadius=args.STRAIN_NEIGHBOUR_RADIUS,
nProcesses=args.PROCESSES,
verbose=True)
......
......@@ -351,7 +351,7 @@ def estimateDisplacementFromQuadraticFit(fieldCoords, displacements, pointsToEst
return estimatedDisplacements
def interpolatePhiField(fieldCoords, PhiField, pointsToInterpolate, nNeighbours=None, neighbourRadius=None, mode="all", nProcesses=nProcessesDefault, verbose=False):
def interpolatePhiField(fieldCoords, PhiField, pointsToInterpolate, nNeighbours=None, neighbourRadius=None, mode="all", neighbourDistanceWeight='inverse', nProcesses=nProcessesDefault, verbose=False):
"""
This function interpolates components of a Phi field at a given number of points, using scipy's KD-tree to find neighbours.
......@@ -380,9 +380,14 @@ def interpolatePhiField(fieldCoords, PhiField, pointsToInterpolate, nNeighbours=
mode : string, optional
Interpolate the whole Phi, just the rigid part, or just the displacement?
Corresponding options are "all", "rigid", "disp"
Corresponding options are 'all', 'rigid', 'disp'
Default = "all"
neighbourDistanceWeight : string, optional
How to weight neigbouring points?
Possible approaches: inverse of distance, gaussian weighting, straight average, median
Corresponding options: 'inverse', 'gaussian', 'mean', 'median'
nProcesses : integer, optional
Number of processes for multiprocessing
Default = number of CPUs in the system
......@@ -407,12 +412,11 @@ def interpolatePhiField(fieldCoords, PhiField, pointsToInterpolate, nNeighbours=
# extract the Phi matrices of the bad points
outputPhiField = numpy.zeros((numberOfPointsToInterpolate, 4, 4), dtype=PhiField.dtype)
if mode not in ["all", "rigid", "disp"]:
print("spam.DIC.interpolatePhiField(): mode argument should either be 'all', 'rigid' or 'disp'")
assert mode in ["all", "rigid", "disp"], "spam.DIC.interpolatePhiField(): mode argument should either be 'all', 'rigid', or 'disp'"
assert neighbourDistanceWeight in ["inverse", "gaussian", "mean", "median"], "spam.DIC.interpolatePhiField(): neighbourDistanceWeight argument should be 'inverse', 'gaussian', 'mean', or 'median'"
# check if neighbours are selected based on radius or based on number, default is radius
if (nNeighbours is None) == (neighbourRadius is None):
print("spam.DIC.interpolatePhiField(): One and only one of nNeighbours and neighbourRadius must be passed")
assert (nNeighbours is None) != (neighbourRadius is None), "spam.DIC.interpolatePhiField(): One and only one of nNeighbours and neighbourRadius must be passed"
if nNeighbours is not None:
ball = False
elif neighbourRadius is not None:
......@@ -423,7 +427,12 @@ def interpolatePhiField(fieldCoords, PhiField, pointsToInterpolate, nNeighbours=
pointToInterpolate = pointsToInterpolate[pointNumber]
outputPhi = numpy.zeros((4, 4), dtype=PhiField.dtype)
outputPhi[-1] = [0, 0, 0, 1]
if mode == 'disp':
outputPhi[0:3, 0:3] = numpy.eye(3)
#######################################################
### Find neighbours
#######################################################
if ball:
# Ball lookup
indices = treeCoord.query_ball_point(pointToInterpolate, neighbourRadius)
......@@ -438,34 +447,42 @@ def interpolatePhiField(fieldCoords, PhiField, pointsToInterpolate, nNeighbours=
indices = numpy.array(indices)
distances = numpy.array(distances)
#print(numpy.max(indices))
# Check if there is only one neighbour
#######################################################
### Check if there is only one neighbour
#######################################################
if indices.size == 1:
if mode in ["all", "rigid"]: # We need to interpolate all 12 components of Phi
outputPhi = PhiField[indices].copy()
if mode == "rigid":
outputPhi = spam.deformation.computeRigidPhi(outputPhi)
else: # interpolate only displacements
outputPhi[0:3, 0:3] = numpy.eye(3)
outputPhi[0:3, -1] = PhiField[indices, 0:3, -1].copy()
return pointNumber, outputPhi
# If > 1 neighbour, interpolate based on distance
#######################################################
### If > 1 neighbour, interpolate Phi or displacements
#######################################################
else:
weights = (1/(distances+tol))
weightsTotal = float(weights.sum())
if mode in ["all", "rigid"]: # We need to interpolate all 12 components of Phi
for neighbourNumber, neighbourIndex in enumerate(indices):
weightInv = weights[neighbourNumber] / weightsTotal
outputPhi[:-1] += PhiField[neighbourIndex][:-1]*weightInv
if neighbourDistanceWeight == 'inverse':
weights = (1/(distances+tol))
elif neighbourDistanceWeight == 'gaussian':
# This could be an input variable VVVVVVVVVVVVVVVVVVVVVV--- the gaussian weighting distance
weights = numpy.exp(-distances**2/numpy.max(distances/2)**2)
elif neighbourDistanceWeight == 'mean':
weights = numpy.ones_like(distances)
elif neighbourDistanceWeight == 'median':
# is this the equivalent kernel to a median, we think so...
weights = numpy.zeros_like(distances)
weights[len(distances)//2] = 1
if mode == 'disp':
outputPhi[0:3,-1] = numpy.sum(PhiField[indices, 0:3,-1] * weights[:, numpy.newaxis], axis=0) / weights.sum()
else:
outputPhi[:-1] = numpy.sum(PhiField[indices, :-1] * weights[:, numpy.newaxis, numpy.newaxis], axis=0) / weights.sum()
if mode == "rigid":
outputPhi = spam.deformation.computeRigidPhi(outputPhi)
else: # interpolate only displacements
outputPhi[0:3, 0:3] = numpy.eye(3)
for neighbourNumber, neighbourIndex in enumerate(indices):
weightInv = weights[neighbourNumber] / weightsTotal
outputPhi[0:3,-1] += PhiField[neighbourIndex][0:3,-1]*weightInv
return pointNumber, outputPhi
......
......@@ -161,7 +161,7 @@ def FfieldRegularQ8(displacementField, nodeSpacing, nProcesses=nProcessesDefault
return Ffield
def FfieldRegularGeers(displacementField, nodeSpacing, neighbourRadius=1.5, mask=True, bruteForceDistance=False, nProcesses=nProcessesDefault, verbose=True):
def FfieldRegularGeers(displacementField, nodeSpacing, neighbourRadius=1.5, 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.
......@@ -187,11 +187,6 @@ def FfieldRegularGeers(displacementField, nodeSpacing, neighbourRadius=1.5, mask
Avoid non-correlated NaN points in the displacement field?
Default = True
bruteForceDistance : bool, optional
Use scipy.spatial.KDtree.query_ball_point to obtain neighbours (bruteForceDistance=False),
or explicitly compute distance to each point (bruteForceDistance=True).
Default = False
nProcesses : integer, optional
Number of processes for multiprocessing
Default = number of CPUs in the system
......@@ -234,23 +229,16 @@ def FfieldRegularGeers(displacementField, nodeSpacing, neighbourRadius=1.5, mask
nodeSpacing[2]:dims[2]*nodeSpacing[2]+nodeSpacing[2]:nodeSpacing[2]].reshape(3, nNodes).T
#Get non-nan displacements
if mask:
goodPointsMask = numpy.isfinite(displacementField[:,:,:,0].reshape(nNodes))
badPointsMask = numpy.isnan( displacementField[:,:,:,0].reshape(nNodes))
#Flattened variables
fieldCoordsFlatGood = fieldCoordsFlat[goodPointsMask]
displacementFieldFlatGood = displacementFieldFlat[goodPointsMask]
#set bad points to nan
FfieldFlat[badPointsMask] = numpy.eye(3) * numpy.nan
else:
#Flattened variables
goodPointsMask = numpy.ones(nNodes, dtype=bool)
fieldCoordsFlatGood = fieldCoordsFlat
displacementFieldFlatGood = displacementFieldFlat
goodPointsMask = numpy.isfinite(displacementField[:,:,:,0].reshape(nNodes))
badPointsMask = numpy.isnan( displacementField[:,:,:,0].reshape(nNodes))
#Flattened variables
fieldCoordsFlatGood = fieldCoordsFlat[goodPointsMask]
displacementFieldFlatGood = displacementFieldFlat[goodPointsMask]
#set bad points to nan
FfieldFlat[badPointsMask] = numpy.eye(3) * numpy.nan
#build KD-tree for neighbour identification
if not bruteForceDistance:
treeCoord = scipy.spatial.KDTree(fieldCoordsFlatGood)
treeCoord = scipy.spatial.KDTree(fieldCoordsFlatGood)
# Output array for good points
FfieldFlatGood = numpy.zeros_like(FfieldFlat[goodPointsMask])
......@@ -266,16 +254,9 @@ def FfieldRegularGeers(displacementField, nodeSpacing, neighbourRadius=1.5, mask
m0 = numpy.zeros((3))
mt = numpy.zeros((3))
if bruteForceDistance:
## Option 1: Manual calculation of distance
## superslow alternative to KDtree
distances = numpy.sqrt(numpy.sum(numpy.square(fieldCoordsFlatGood - centralNodePosition), axis=1))
ind = numpy.where(distances < neighbourRadius*max(nodeSpacing))[0]
else:
## Option 2: KDTree on distance
# KD-tree will always give the current point as zero-distance
ind = treeCoord.query_ball_point(centralNodePosition, neighbourRadius*max(nodeSpacing))
## Option 2: KDTree on distance
# KD-tree will always give the current point as zero-distance
ind = treeCoord.query_ball_point(centralNodePosition, neighbourRadius*max(nodeSpacing))
# We know that the current point will also be included, so remove it from the index list.
ind = numpy.array(ind)
......
......@@ -1044,46 +1044,6 @@ def regularStrainParser(parser):
dest='Q8',
help='Use Q8 element interpolation? More noisy and strain values not centred on displacement points')
parser.add_argument('-cif',
'--correct-input-field',
action="store_true",
dest='CORRECT_FIELD',
help='Activates correction of the input displacement field')
parser.add_argument('-cni',
'--correct-neighbours-for-field-interpolation',
type=int,
default=12,
dest='CORRECT_NEIGHBOURS',
help="Neighbours for field interpolation. Default = 12")
parser.add_argument('-cmf',
'--correct-median-filter',
action="store_true",
dest='CORRECT_MEDIAN_FILTER',
help="Activates an overall median filter on the input displacement field")
parser.add_argument('-cmfr',
'--correct-median-filter-radius',
type=int,
default=2,
dest='CORRECT_MEDIAN_FILTER_RADIUS',
help="Radius of median filter for correction of input displacement field. Default = 2")
parser.add_argument('-cdp',
'--correct-delta-phi-norm',
type=numpy.float,
default=0.001,
dest='CORRECT_DELTA_PHI_NORM',
help="Delta Phi norm for a return status = 1 correlation window to consider the point good. Default = 0.001")
parser.add_argument('-cpscc',
'--correct-pixel-search-cc',
type=numpy.float,
default=0,
dest='CORRECT_PIXEL_SEARCH_CC',
help="Pixel search correlation coefficient to consider the point good. Default = 0")
parser.add_argument('-od',
'--out-dir',
type=str,
......
......@@ -337,10 +337,6 @@ class testAll(unittest.TestCase):
decomposedF = spam.deformation.decomposeFfield(Ffield, ["vol"], verbose=False)
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, verbose=False)
decomposedF = spam.deformation.decomposeFfield(Ffield, ["vol"], verbose=False)
self.assertAlmostEqual(decomposedF['vol'][1:-1,1:-1,1:-1].mean(), 3. * (shrink - 1.), places=3)
def test_geers(self):
......@@ -499,11 +495,6 @@ class testAll(unittest.TestCase):
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
########################################################
......@@ -546,8 +537,8 @@ class testAll(unittest.TestCase):
self.assertEqual(numpy.isnan(Ffield[:,:,:,2,2]).sum(), 1)
self.assertEqual(numpy.isnan(Ffield2D[:,:,:,2,2]).sum(), 1)
Ffield = spam.deformation.FfieldRegularGeers(displacements, nodeSpacing, mask=True)
Ffield2D = spam.deformation.FfieldRegularGeers(displacements2D, nodeSpacing2D, mask=True)
Ffield = spam.deformation.FfieldRegularGeers(displacements, nodeSpacing)
Ffield2D = spam.deformation.FfieldRegularGeers(displacements2D, nodeSpacing2D)
self.assertEqual(numpy.isnan(Ffield[:,:,:,2,2]).sum(), 1)
self.assertEqual(numpy.isnan(Ffield2D[:,:,:,2,2]).sum(), 1)
......@@ -795,6 +786,18 @@ class testAll(unittest.TestCase):
self.assertTrue(numpy.allclose(interpolatedPhi[0, 0:3, -1], transformation['t'], atol=0.01))
self.assertTrue(numpy.allclose(interpolatedPhi[0, 0:3, 0:3], numpy.eye(3), atol=0.01))
# Test all new interplation options
interpolatedPhi = spam.DIC.interpolatePhiField(fieldCoords, fieldValues, interpCoords, nNeighbours=8, neighbourDistanceWeight='gaussian')
self.assertTrue(numpy.allclose(interpolatedPhi[0, 0:3, -1], transformation['t'], atol=0.01))
self.assertTrue(numpy.allclose(interpolatedPhi[0, 0:3, 0:3], Phi[0:3, 0:3], atol=0.01))
interpolatedPhi = spam.DIC.interpolatePhiField(fieldCoords, fieldValues, interpCoords, nNeighbours=8, neighbourDistanceWeight='mean')
self.assertTrue(numpy.allclose(interpolatedPhi[0, 0:3, -1], transformation['t'], atol=0.01))
self.assertTrue(numpy.allclose(interpolatedPhi[0, 0:3, 0:3], Phi[0:3, 0:3], atol=0.01))
interpolatedPhi = spam.DIC.interpolatePhiField(fieldCoords, fieldValues, interpCoords, nNeighbours=8, neighbourDistanceWeight='median')
self.assertTrue(numpy.allclose(interpolatedPhi[0, 0:3, -1], transformation['t'], atol=0.01))
self.assertTrue(numpy.allclose(interpolatedPhi[0, 0:3, 0:3], Phi[0:3, 0:3], atol=0.01))
def test_computeRigidPhi(self):
transformation = {'t': [5., 0., 0.],
......
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