spam.DIC package

spam.DIC package#

Submodules#

spam.DIC.ddic module#

spam.DIC.ddic.ddic(im1, im2, lab1, labelsToCorrelate=None, PhiField=None, boundingBoxes=None, centresOfMass=None, processes=None, labelDilate=1, margin=5, maskOthers=True, volThreshold=100, multiScaleBin=1, updateGrad=False, correlateRigid=True, maxIter=50, deltaPhiMin=0.001, interpolationOrder=1, debug=False, twoD=False)[source]#

spam.DIC.deform module#

spam.DIC.deform.applyPhi(im, Phi=None, PhiCentre=None, interpolationOrder=1)[source]#

Deform a 3D image using a deformation function “Phi”, applied using spam’s C++ interpolator. Only interpolation order = 1 is implemented.

Parameters:

im (3D numpy array) – 3D numpy array of grey levels to be deformed
Phi (4x4 array, optional) – “Phi” deformation function. Highly recommended additional argument (why are you calling this function otherwise?)
PhiCentre (3x1 array of floats, optional) – Centre of application of Phi. Default = (numpy.array(im1.shape)-1)/2.0 i.e., the centre of the image
interpolationOrder (int, optional) – Order of image interpolation to use, options are either 0 (strict nearest neighbour) or 1 (trilinear interpolation) Default = 1

Returns:

imDef – Deformed greyscales by Phi

Return type:

3D array

spam.DIC.deform.applyPhiPython(im, Phi=None, PhiCentre=None, interpolationOrder=3)[source]#

Deform a 3D image using a deformation function “Phi”, applied using scipy.ndimage.map_coordinates Can have orders > 1 but is hungry in memory.

Parameters:

im (2D/3D numpy array) – 2D/3D numpy array of grey levels to be deformed
Phi (4x4 array, optional) – “Phi” linear deformation function. Highly recommended additional argument (why are you calling this function otherwise?)
PhiCentre (3x1 array of floats, optional) – Centre of application of Phi. Default = (numpy.array(im1.shape)-1)/2.0 i.e., the centre of the image
interpolationOrder (int, optional) – Order of image interpolation to use. This value is passed directly to scipy.ndimage.map_coordinates as “order”. Default = 3

Returns:

imSub – Deformed greyscales by Phi

Return type:

3D array

spam.DIC.deform.applyPhiField(im, fieldCoordsRef, PhiField, imMaskDef=None, displacementMode='applyPhi', nNeighbours=8, interpolationOrder=1, nProcesses=32, verbose=False)[source]#

Deform a 3D image using a field of deformation functions “Phi” coming from a regularGrid, applied using scipy.ndimage.map_coordinates.

Parameters:

im (3D array) – 3D array of grey levels to be deformed
fieldCoordsRef (2D array, optional) – nx3 array of n points coordinates (ZYX) centre where each deformation function “Phi” has been calculated
PhiField (3D array, optional) – nx4x4 array of n points deformation functions
imMaskDef (3D array of bools, optional) – 3D array same size as im but DEFINED IN THE DEFORMED CONFIGURATION which should be True in the pixels to fill in in the deformed configuration. Default = None
displacementMode (string, optional) – How do you want to calculate displacements? With “interpolate” they are just interpolated from the PhiField With “applyPhi” each neighbour’s Phi function is applied to the pixel position and the resulting translations weighted and summed. Default = “applyPhi”
nNeighbours (int, optional) – Number of the nearest neighbours to consider #This OR neighbourRadius must be set. Default = 8
interpolationOrder (int, optional) – Order of image interpolation to use. This value is passed directly to scipy.ndimage.map_coordinates as “order”. Default = 1
nProcesses (integer, optional) – Number of processes for multiprocessing Default = number of CPUs in the system
verbose (boolean, optional) – Print progress? Default = True

Returns:

imDef – deformed greylevels by a field of deformation functions “Phi”

Return type:

3D array

spam.DIC.deform.binning(im, binning, returnCropAndCentre=False)[source]#

This function downscales images by averaging NxNxN voxels together in 3D and NxN pixels in 2D. This is useful for reducing data volumes, and denoising data (due to averaging procedure).

Parameters:

im (2D/3D numpy array) – Input measurement field
binning (int) – The number of pixels/voxels to average together
returnCropAndCentre (bool (optional)) – Return the position of the centre of the binned image in the coordinates of the original image, and the crop Default = False

Returns:

imBin (2/3D numpy array) – binning-binned array
(otherwise if returnCropAndCentre) (list containing:) – imBin, topCrop, bottomCrop centre of imBin in im coordinates (useful for re-stitching)

Notes

Here we will only bin pixels/voxels if they is a sufficient number of neighbours to perform the binning. This means that the number of pixels that will be rejected is the dimensions of the image, modulo the binning amount.

The returned volume is computed with only fully binned voxels, meaning that some voxels on the edges may be excluded. This means that the output volume size is the input volume size / binning or less (in fact the crop in the input volume is the input volume size % binning

spam.DIC.deform.applyMeshTransformation(im, points, connectivity, displacements, imTetLabel=None, nThreads=1)[source]#

This function deforms an image based on a tetrahedral mesh and nodal displacements (normally from Global DVC), using the mesh’s shape functions to interpolate.

Parameters:

im (3D numpy array of greyvalues) – Input image to be deformed
points (m x 3 numpy array) – M nodal coordinates in reference configuration
connectivity (n x 4 numpy array) – Tetrahedral connectivity generated by spam.mesh.triangulate() for example
displacements (m x 3 numpy array) – M displacements defined at the nodes
imTetLabel (3D numpy array of ints (optional)) – Pixels labelled with the tetrahedron (i.e., line number in connectivity matrix) they belong to. If this is not passed, it’s calculated in this function (can be slow). WARNING: This is in the deformed configuration.
nThreads (int (optional, default=1)) – The number of threads used for the cpp parallelisation.

Returns:

imDef – Deformed image

Return type:

3D numpy array of greyvalues

spam.DIC.globalDVC module#

Set of functions for performing mechanically regularised global DVC. The regularisation implementation is taken from [Mendoza2019].

Here is how it can be used:

import spam.mesh
import spam.DIC

# Mesh
points, connectivity = spam.mesh.createCylinder([-24, 27], 12.5, 97.3, 10, zOrigin=-4.2)

# Regularisation: parameters
myParameters = {
    "young": 10,
    "poisson": 0.25,
    "dirichlet": {
        "z_start": {"dof": [0]},
        "z_end": {"dof": [0]},
    },
}
p = spam.DIC.regularisationParameters(myParameters)

# Regularisation step 1: get labels
labels = spam.DIC.surfaceLabels(points, surfaces=p["surfaces"])

# Regularisation step 2: build regularisation matrix
regularisationMatrix, regularisationField = spam.DIC.regularisationMatrix(
    points,
    connectivity,
    young=p["young"],
    poisson=p["poisson"],
    xiBulk=p["xi"],
    dirichlet=p["dirichlet"],
    periods=p["periods"],
    labels=labels,
)

# Global DVC
globalCorrelation(
    im1,
    im2,
    points,
    connectivity,
    regularisationMatrix=regularisationMatrix,
    regularisationField=regularisationField,
)

[Mendoza2019] (1,2,3,4,5)

A. Mendoza, J. Neggers, F. Hild, S Roux (2019). Complete mechanical regularization applied to digital image and volume correlation. Computer Methods in Applied Mechanics and Engineering, Volume 355, 1 October 2019, Pages 27-43 https://doi.org/10.1016/j.cma.2019.06.005

Make a link to the script

spam.DIC.globalDVC.surfaceLabels(points, surfaces=[], connectivity=None)[source]#

Creates a label for each points based on a list of keywords that corresponds to surfaces (ie. ["z_start", "z_end"]). The label value is based on the position of the surface in the list.

Parameters:

points (Nx3 array) – List of coordinates of the mesh nodes.
surfaces (list of string) –
A list of keywords corresponding to surfaces.
- z_start: plane at z == min(points[:,0])
- z_end: plane at z == max(points[:,0])
- y_start: plane at y == min(points[:,1])
- y_end: plane at y == max(points[:,1])
- x_start: plane at x == min(points[:,2])
- x_end: plane at x == max(points[:,2])
- z_lateral: lateral surface of a cylinder oriented in the first direction.
- y_lateral: lateral surface of a cylinder oriented in the second direction.
- x_lateral: lateral surface of a cylinder oriented in the third direction.
connectivity (array (only for debug purposes)) – Connectivity matrix. If set, creates a VTK file with labels.

Returns:

Surface label for each points.

Return type:

N x 1 array

Note

Surface labels start at 1, 0 corresponding to bulk or not specified surfaces.
Points belong to a single surface. The first surface in surfaces prevails.

Example

>>> import spam.mesh
>>> import spam.DIC
>>>
>>> # create mesh
>>> points, connectivity = spam.mesh.createCylinder([-24, 27], 12.5, 97.3, 10, zOrigin=-4.2)
>>> # compute labels for bottom and top surface only
>>> labels = spam.DIC.surfaceLabels(points, surfaces=["z_start", "z_end"], connectivity=connectivity)

spam.DIC.globalDVC.projectionMatrices(points, connectivity, labels, dirichlet=[])[source]#

Compute binary diagonal matrices and Laplace-Beltrami operator. Details on the meaning of these matrices can be found in [Mendoza2019] eq. (12) to (14) and eq. (17) to (19).

Parameters:

points (Nx3 array) – List of coordinates of the mesh nodes.
connectivity (Mx4 numpay array) – Connectivity matrice of the mesh (tetrahedra).
labels (Nx1 array) – Surface labels for each points as defined in spam.DIC.surfaceLabels()
dirichlet (list of (int, list, )) –
Each element of the list defines a surface with Dirichlet boundary conditions by a tuple.
- The first element of the tuple is the surface label which should belong to labels.
- The second element of the tuple is a list of degrees of freedom directions to consider for the Dirichlet boundary conditions.
- Extra elements in the tuple are ignored.
```
dirichlet = [
    # surface 1, dof in x, y
    (
        1,
        [1, 2],
    ),
    # surface 3, dof in z
    (
        3,
        [0],
    ),
]
```

Returns:

3Nx3N array – $D_m$ The binary diagonal matrix corresponding to all dofs of the bulk and neumann surfaces nodes (ie the non dirichlet)
List of 3Nx3N array – $D_{S_i}$ Binary diagonal matrices corresponding to the dofs of each Dirichlet surfaces
List of 3Nx3N array – $L_{S_i}$ List of Laplace-Beltrami operators corresponding to the same Dirichlet surfaces.

spam.DIC.globalDVC.isochoricField(points, periods=3, connectivity=None)[source]#

Helper function to build the isochoric test function used to normalised the mechanical regularisation functionals. The function is a shear displacement field with its wave vector in the direction of the longest dimension of the mesh.

\[\textbf{v}(\textbf{x}) = \sin(2 \pi \textbf{k}\cdot \textbf{x} )\]

[Mendoza2019] eq. (25).

Parameters:

points (Nx3 array) – List of coordinates of the mesh nodes.
periods (int) – The number of periods of the shear wave.
connectivity (array (only for debug purposes)) – Connectivity matrix. If set, creates a VTK file with the field.

Returns:

float – The magnitude of the wave vector $\bf{k}$.
Nx3 array – The field $\bf{v}$ at each nodes.

spam.DIC.globalDVC.regularisationMatrix(points, connectivity, young=1, poisson=0.25, voxelSize=1, xiBulk=None, dirichlet=[], labels=[], periods=3)[source]#

Computes the mechanical regularisation matrix $M_{reg}$ for the global DVC: $$(M_{reg} + M_{c}) \delta\textbf{u} = \textbf{b} - M_{reg} \textbf{u} $$ where $$M_{reg} = M_m + \sum_{i} M_{S_i}$$ corresponds to both bulk/Neuman and Dirichlet surfaces contributions ([Mendoza2019] eq. (29) and (30)).

Parameters:

points (Nx3 array) – List of coordinates of the mesh nodes.
connectivity (Mx4 array) – Connectivity matrix of the mesh.
young (float (optional)) – The Young modulus used to build the stiffness matrix $K$ from which $M_m$ and $M_{S_i}$ are derived. Default = 1 (it is useless to change this value if you don’t impose forces with meaningful units on the Neumann surfaces)
poisson (float (optional)) – The Poisson ratio used to build the stiffness matrix $K$ from which $M_m$ and $M_{S_i}$ are derived. Default = 0.25
voxelSize (float (optional)) – The size of a voxel in the same unit used to set the Young modulus assuming that the mesh geometry is in voxels. This voxel size $v_s$ is used to convert the Young modulus in Newton per voxels squared. Default = 1 (it is useless to change this value if you don’t impose forces with meaningful units on the Neumann surfaces).
xiBulk (float (optional)) – The regularisation length $\xi_m$ of the bulk/Neumann contribution $M_{reg}$. It has to be compared to the characteristic length of the mesh. The bigger the regularisation length the more important the regularisation contribution. If None it is taken to be equal to the mesh characteristic length. Default = None
dirichlet (list of (int, list, float) (optional)) –
Each element of the list defines a surface with Dirichlet boundary conditions by a tuple.
- The first element of the tuple is the surface label which should belong to labels.
- The second element of the tuple is a list of degrees of freedom directions to consider for the Dirichlet boundary conditions.
- The third element of the tuple correspond to the regularisation length of the surface $\xi_{S_i}$
```
dirichlet = [
    # surface 1, dof in x, y, xi = 0.1
    [1, [1, 2], 0.1],
    # surface 3, dof in z, xi not set
    [3, [0], None],
]
```
labels (Nx1 array (optional)) – Surface labels for each points as defined in spam.DIC.surfaceLabels(). Mandatory only if dirichlet is set.
periods (float (optional)) – Number of periods of the isochoric function $v$ used to compute the normalized energy ([Mendoza2019] eq. (27)). $v$ is computed with spam.DIC.isochoricField() with a default number of periods of 3.

Returns:

3Nx3N array – $M_{req}$ the regularisation matrix.
Nx3 array – $v$ the isochoric field at each nodes of the mesh.

Note

The dirichlet argument is compatible with the one in spam.DIC.projectionMatrices()
As long as you don’t impose forces to the Neumann surfaces it is usless to set a specific Young modulus and a voxel size.
Imposing forces is not implemented yet :)

Warning

This function is not tested yet

spam.DIC.globalDVC.regularisationParameters(userParameters)[source]#

Convert a user friendly dictionary of parameters into a dictionary of variables compatible with the regularisation functions of this module.

Parameters:

userParameters ((str | dict)) –

The user parameters for the mechanical regularisation scheme. It can be:

if str: the path to a YAML file. A dummy example can be downloaded here: RegularisationParameter
if dict: a dictionary containing the following keys and values:

{
    # bulk / Neumann regularisation
    "young": 10,  # mandatory (the Young modulus)
    "poisson": 0.25,  # mandatory (the Poisson ratio)
    "xiBulk": 30,  # optional (the bulk/Neumann regularisation lentgh)
    "periods": 3,  # the number of periods of the isochoric function (optional)
    "voxelSize": 0.01,  # the voxel size (optional)
    # Information about the Dirichlet surface regularisation
    # Each surface is defined by a search keyword among
    # z_start, z_end, y_start, y_end, x_start and x_end for plane lookups
    # z_lateral, y_lateral, x_lateral for cylinder lateral surface lookups
    "dirichlet": {
        "z_start": {  # surface label 1
            "xi": 5,  # the regularisation length (optional)
            "dof": [0, 1, 2],  # mandatory
        },
        "z_end": {"dof": [0]},  # surface label 2  # mandatory
    },
}

Returns:

A dictionary with the variables needed for the regularisation functions. The dictionary keys are named to match the function’s signatures:

surface: the dirichlet surfaces for spam.DIC.surfaceLabels()
dirichlet: the dirichlet surfaces for spam.DIC.regularisationMatrix()
young: the Young modulus for spam.DIC.regularisationMatrix()
poisson: the Poisson ratio for spam.DIC.regularisationMatrix()
xiBulk: the bulk characteristic lentgh for spam.DIC.regularisationMatrix()
periods: the Poisson ratio for spam.DIC.regularisationMatrix()
voxelSize: the Poisson ratio for spam.DIC.regularisationMatrix()

Return type:

dict

spam.DIC.globalDVC.globalCorrelation(im1, im2, points, connectivity, regularisationMatrix=None, regularisationField=None, initialDisplacements=None, convergenceCriterion=0.01, maxIterations=20, medianFilterEachIteration=False, debugFiles=False, prefix='globalCorrelation', nThreads=None)[source]#

Global DVC (works only with 3D images).

Parameters:

im1 (3D array) – Reference image in which the mesh is defined
im2 (3D array) – Deformed image, should be same size as im1
points (M x 3 array) – M nodal coordinates in reference configuration
connectivity (N x 4 array) – connectivityhedral connectivity generated by spam.mesh.triangulate() for example
regularisationMatrix – Mechanical regularisation stiffness matrix. If None no mechanical regularisation is applied. First output of spam.DIC.regularisationMatrix

Returns:

displacements – (converged?) Nodal displacements

Return type:

N x 3 array of floats

Example

>>> import spam.DIC
>>> spam.DIC.globalCorrelation(
    imRef,
    imDef
)

spam.DIC.grid module#

This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.

This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.

You should have received a copy of the GNU General Public License along with this program. If not, see <http://www.gnu.org/licenses/>.

spam.DIC.grid.makeGrid(imageSize, nodeSpacing)[source]#

Define a grid of correlation points.

Parameters:

imageSize (3-item list) – Size of volume to spread the grid inside
nodeSpacing (3-item list or int) – Spacing between nodes

Returns:

nodePositions – Array containing Z, Y, X positions of each point in the grid

Return type:

nPointsx3 numpy.array

spam.DIC.grid.getImagettes(im1, nodePosition, halfWindowSize, Phi, im2, searchRange, im1mask=None, im2mask=None, minMaskCoverage=0.0, greyThreshold=[-inf, inf], applyF='all', twoD=False)[source]#

This function is responsible for extracting correlation windows (“imagettes”) from two larger images (im1 and im2). Both spam.correlate.pixelSearch and spam.correlate.register[Multiscale] want a fixed, smaller imagette1 and a larger imagette 2 in which to search/interpolate.

Parameters:

im1 (3D numpy array) – This is the large input reference image
nodePosition (3-component numpy array of ints) – This defines the centre of the window to extract from im1. Note: for 2D Z = 0
halfWindowSize (3-component numpy array of ints) – This defines the half-size of the correlation window, i.e., how many pixels to extract in Z, Y, X either side of the centre. Note: for 2D Z = 0
Phi (4x4 numpy array of floats) – Phi matrix representing the movement of imagette1, if not equal to I, imagette1 is deformed by the non-translation parts of Phi (F) and the displacement is added to the search range (see below)
im2 (3D numpy array) – This is the large input deformed image
searchRange (6-component numpy array of ints) – This defines where imagette2 should be extracted with respect to imagette1’s position in im1. The 6 components correspond to [ Zbot Ztop Ybot Ytop Xbot Xtop ]. If Z, Y and X values are the same, then imagette2 will be displaced and the same size as imagette1. If ‘bot’ is lower than ‘top’, imagette2 will be larger in that dimension
im1mask (3D numpy array, optional) – This needs to be same size as im1, but can be None if no mask is wanted. This defines a mask for zones to correlate in im1, 0 means zone not to correlate Default = None
im2mask (3D numpy array, optional) – This needs to be same size as im2, but can be None if no mask is wanted. This defines a mask for zones to correlate in im2, 0 means zone not to correlate Default = None
minMaskCoverage (float, optional) – Threshold for imagette1 non-mask coverage, i.e. how much of imagette1 can be full of mask before it is rejected with returnStatus = -5? Default = 0
greyThreshold (two-component list of floats, optional) – Bottom and top threshold values for mean value of imagette1 to reject it with returnStatus = -5 Default = no threshold
applyF (string, optional) – If a non-identity Phi is passed, should the F be applied to the returned imagette1? Options are: ‘all’, ‘rigid’, ‘no’ Default = ‘all’ Note: as of January 2021, it seems to make more sense to have this as ‘all’ for pixelSearch, and ‘no’ for local DIC
twoD (bool, optional) – Are the images two-dimensional?

Returns:

‘imagette1’ : 3D numpy array,

’imagette1mask’: 3D numpy array of same size as imagette1 or None,

’imagette2’: 3D numpy array, bigger or equal size to imagette1

’imagette2mask’: 3D numpy array of same size as imagette2 or None,

’returnStatus’: int,: Describes success in extracting imagette1 and imagette2. If == 1 success, otherwise negative means failure.
’pixelSearchOffset’: 3-component list of ints: Coordinates of the top of the pixelSearch range in im1, i.e., the displacement that needs to be added to the raw pixelSearch output to make it a im1 -> im2 displacement

Return type:

Dictionary

spam.DIC.kinematics module#

This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.

This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.

You should have received a copy of the GNU General Public License along with this program. If not, see <http://www.gnu.org/licenses/>.

spam.DIC.kinematics.estimateLocalQuadraticCoherency(points, displacements, neighbourRadius=None, nNeighbours=None, epsilon=0.1, nProcesses=32, verbose=False)[source]#

This function computes the local quadratic coherency (LQC) of a set of displacement vectors as per Masullo and Theunissen 2016. LQC is the average residual between the point’s displacement and a second-order (parabolic) surface Phi. The quadratic surface Phi is fitted to the point’s closest N neighbours and evaluated at the point’s position. Neighbours are selected based on: radius (default option) or number (activated if nNeighbours is not None). A point with a LQC value smaller than a threshold (0.1 in Masullo and Theunissen 2016) is classified as coherent

Parameters:

points (n x 3 numpy array of floats) – Coordinates of the points Z, Y, X
displacements (n x 3 numpy array of floats) – Displacements of the points
neighbourRadius (float, optional) – Distance in pixels around the point to extract neighbours. This OR nNeighbours must be set. Default = None
nNeighbours (int, optional) – Number of the nearest neighbours to consider This OR neighbourRadius must be set. Default = None
epsilon (float, optional) – Background error as per (Westerweel and Scarano 2005) Default = 0.1
nProcesses (integer, optional) – Number of processes for multiprocessing Default = number of CPUs in the system
verbose (boolean, optional) – Print progress? Default = True

Returns:

LQC – The local quadratic coherency for each point

Return type:

n x 1 array of floats

Note

Based on: https://gricad-gitlab.univ-grenoble-alpes.fr/DVC/pt4d

spam.DIC.kinematics.estimateDisplacementFromQuadraticFit(fieldCoords, displacements, pointsToEstimate, neighbourRadius=None, nNeighbours=None, epsilon=0.1, nProcesses=32, verbose=False)[source]#

This function estimates the displacement of an incoherent point based on a local quadratic fit of the displacements of N coherent neighbours, as per Masullo and Theunissen 2016. A quadratic surface Phi is fitted to the point’s closest coherent neighbours. Neighbours are selected based on: radius (default option) or number (activated if nNeighbours is not None)

Parameters:

fieldCoords (n x 3 numpy array of floats) – Coordinates of the points Z, Y, X where displacement is defined
displacements (n x 3 numpy array of floats) – Displacements of the points
pointsToEstimate (m x 3 numpy array of floats) – Coordinates of the points Z, Y, X where displacement should be estimated
neighbourRadius (float, optional) – Distance in pixels around the point to extract neighbours. This OR nNeighbours must be set. Default = None
nNeighbours (int, optional) – Number of the nearest neighbours to consider This OR neighbourRadius must be set. Default = None
epsilon (float, optional) – Background error as per (Westerweel and Scarano 2005) Default = 0.1
nProcesses (integer, optional) – Number of processes for multiprocessing Default = number of CPUs in the system
verbose (boolean, optional) – Print progress? Default = True

Returns:

displacements – The estimated displacements at the requested positions.

Return type:

m x 3 array of floats

Note

Based on https://gricad-gitlab.univ-grenoble-alpes.fr/DVC/pt4d

spam.DIC.kinematics.interpolatePhiField(fieldCoords, PhiField, pointsToInterpolate, nNeighbours=None, neighbourRadius=None, interpolateF='all', neighbourDistanceWeight='inverse', checkPointSurrounded=False, nProcesses=32, verbose=False)[source]#

This function interpolates components of a Phi field at a given number of points, using scipy’s KD-tree to find neighbours.

Parameters:

fieldCoords (2D array) – nx3 array of n points coordinates (ZYX) centre where each deformation function Phi has been measured
PhiField (3D array) – nx4x4 array of n points deformation functions
pointsToInterpolate (2D array) – mx3 array of m points coordinates (ZYX) Points where the deformation function Phi should be interpolated
nNeighbours (int, optional) – Number of the nearest neighbours to consider This OR neighbourRadius must be set. Default = None
neighbourRadius (float, optional) – Distance in pixels around the point to extract neighbours. This OR nNeighbours must be set. Default = None
interpolateF (string, optional) – Interpolate the whole Phi, just the rigid part, or just the displacement? Corresponding options are ‘all’, ‘rigid’, ‘no’ 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’
checkPointSurrounded (bool, optional) – Only interpolate points whose neighbours surround them in Z, Y, X directions (or who fall exactly on a give point)? Default = False
nProcesses (integer, optional) – Number of processes for multiprocessing Default = number of CPUs in the system
verbose (bool, optional) – follow the progress of the function. Default = False

Returns:

PhiField – Interpolated Phi functions at the requested positions

Return type:

mx4x4 array

spam.DIC.kinematics.mergeRegularGridAndDiscrete(regularGrid=None, discrete=None, labelledImage=None, binningLabelled=1, alwaysLabel=True)[source]#

This function corrects a Phi field from the spam-ldic script (measured on a regular grid) by looking into the results from one or more spam-ddic results (measured on individual labels) by applying the discrete measurements to the grid points.

This can be useful where there are large flat zones in the image that cannot be correlated with small correlation windows, but can be identified and tracked with a spam-ddic computation (concrete is a good example).

Parameters:

regularGrid (dictionary) – Dictionary containing readCorrelationTSV of regular grid correlation script, spam-ldic. Default = None
discrete (dictionary or list of dictionaries) – Dictionary (or list thereof) containing readCorrelationTSV of discrete correlation script, spam-ddic. File name of TSV from DICdiscrete client, or list of filenames Default = None
labelledImage (3D numpy array of ints, or list of numpy arrays) – Labelled volume used for discrete computation Default = None
binningLabelled (int) – Are the labelled images and their PhiField at a different bin level than the regular field? Default = 1
alwaysLabel (bool) – If regularGrid point falls inside the label, should we use the label displacement automatically? Otherwise if the regularGrid point has converged should we use that? Default = True (always use Label displacement)

Returns:

Output matrix, with number of rows equal to spam-ldic (the node spacing of the regular grid) and with columns: “NodeNumber”, “Zpos”, “Ypos”, “Xpos”, “Zdisp”, “Ydisp”, “Xdisp”, “deltaPhiNorm”, “returnStatus”, “iterations”

Return type:

Either dictionary or TSV file

spam.DIC.kinematics.getDisplacementFromNeighbours(labIm, DVC, fileName, method='getLabel', neighboursRange=None, centresOfMass=None, previousDVC=None)[source]#

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’
neighboursRange (int) – Parameter controlling the search range to detect neighbours for each method. For ‘getLabel’, it correspond to the size of the subset. Default = meanRadii For ‘mesh’, it correspond to the size of the alpha shape used for carving the mesh. Default = 5*meanDiameter.
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

Returns:

TSV file with the same columns as the input

Return type:

Dictionary

spam.DIC.kinematics.applyRegistrationToPoints(Phi, PhiCentre, points, applyF='all', nProcesses=32, verbose=False)[source]#

This function takes a whole-image registration and applies it to a set of points

Parameters:

Phi (4x4 numpy array of floats) – Measured Phi function to apply to points
PhiCentre (3-component list of floats) – Origin where the Phi is measured (normally the middle of the image unless masked)
points (nx3 numpy array of floats) – Points to apply the Phi to
applyF (string, optional) – The whole Phi is always applied to the positions of the points to get their displacement. This mode only controls what is copied into the F for each point, everything, only rigid, or only displacements? Corresponding options are ‘all’, ‘rigid’, ‘no’ Default = “all”
nProcesses (integer, optional) – Number of processes for multiprocessing Default = number of CPUs in the system
verbose (boolean, optional) – Print progress? Default = True

Returns:

PhiField – Output Phi field

Return type:

nx4x4 numpy array of floats

spam.DIC.kinematics.mergeRegistrationAndDiscreteFields(regTSV, discreteTSV, fileName, regTSVBinRatio=1)[source]#

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 registration, 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:

TSV file with the same columns as the input

Return type:

Dictionary

spam.DIC.ldic module#

spam.DIC.ldic.ldic(im1, im2, nodePositions, hws, skipNodesMask=False, processes=None, im1mask=None, PhiField=None, margin=[-1, 1, -1, 1, -1, 1], maskCoverage=0.5, greyThreshold=[-inf, inf], applyF='no', maxIterations=50, deltaPhiMin=0.001, PhiRigid=True, updateGradient=False, interpolationOrder=1)[source]#

Function to perform local DIC/DVC (i.e., running the spam.DIC.register function) correlating many 2D/3D boxes spread out typically on a regular grid of “nodes”.

Parameters:

im1 (-) – Reference image in which the nodes are defined, and from which the output Phi field is measured
im2 (-) – Deformed image
nodePositions (-) – Nx2 or Nx3 matrix defining centres of boxes. This can be generated with nodePositions, nodesDim = spam.DIC.makeGrid(im1.shape, nodeSpacing)
hws (-) – This defines the half-size of the correlation window, i.e., how many pixels to extract in Z, Y, X either side of the centre. Note: for 2D Z = 0
skipNodes (-) – Vector of N bools which are true when nodes should be skipped. They will simply not be correlated, so ignore the outputs related to these nodes. If you have guesses for them, remember to merge them with the outputs for this function
processes (-) – Number of processes to run the ldic on, by default it’s the number of detected threads on your machine
spam.DIC.register() (- for all other parameters see)

spam.DIC.multimodal module#

This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.

This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.

You should have received a copy of the GNU General Public License along with this program. If not, see <http://www.gnu.org/licenses/>.

spam.DIC.multimodal.multimodalRegistration(im1, im2, phaseDiagram, gaussianParameters, im1mask=None, PhiInit=None, margin=None, maxIterations=50, deltaPhiMin=0.005, interpolationOrder=1, verbose=False, GRAPHS=False, INTERACTIVE=False, sliceAxis=0, suffix='', rootPath='.', BINS=64)[source]#

Perform subpixel image correlation between im1 and im2 – images of the same object acquired with different modalities.

This function will deform im2 until it best matches im1. The matching includes sub-pixel displacements, rotation, and linear straining of the whole image. The correlation of im1, im2 will give a deformationFunction F which maps im1 into im2. Phi(im1(x),im2(F.x)) = 0 As per Equation (3) of Tudisco et al.

“Discrete correlation” can be performed by masking im1.

Parameters:

im1 (3D numpy array) – The greyscale image that will not move
im2 (3D numpy array) – The greyscale image that will be deformed
phaseDiagram (Nbins x Nbins numpy array of ints) – Pre-labelled phase diagram, which is a joint histogram with the phases labelled
gaussianParameters (2D numpy array Nx6) – Key parameter which is the result of a 2D gaussian fit of the first N peaks in the joint histograms of im1 and im2. The 6 parameters of the fit are: φ, μ(im1), μ(im2), and a, b, c, where {a,b,c} are the parameter that can be found for two-dimensional elliptical Gaussian function here: https://en.wikipedia.org/wiki/Gaussian_function, i.e., coupled with im1², im1*im2 and im2² respectively
im1mask (3D numpy array, float, optional) – A mask for the zone to correlate in im1 with NaNs in the zone to not correlate. Default = None, i.e., correlate all of im1 minus the margin
PhiInit (4x4 numpy array, optional) – Initial transformation to apply. Default = numpy.eye(4), i.e., no transformation
margin (int, optional) – Margin, in pixels, to take in im1 and im2 to allow space for interpolation and movement. Default = None (i.e., enough margin for a worse-case 45degree rotation with no displacement)
maxIterations (int, optional) – Maximum number of quasi-Newton interations to perform before stopping. Default = 25
deltaPhiMin (float, optional) – Smallest change in the norm of F (the transformation operator) before stopping. Default = 0.001
interpolationOrder (int, optional) – Order of the greylevel interpolation for applying F to im1 when correlating. Reccommended value is 3, but you can get away with 1 for faster calculations. Default = 3
verbose (bool, optional) – Get to know what the function is really thinking, recommended for debugging only. Default = False
GRAPHS (bool, optional) – Pop up a window showing the image differences (im1-im2) as im1 is progressively deformed.

Returns:

‘transformation’

Dictionary containing:

’t’ : 3-component vector

Z, Y, X displacements in pixels

’r’ : 3-component vector

Z, Y, X components of rotation vector

’Phi’: 4 x 4 numpy array of floats

Deformation function, Phi

’returnStatus’: int

Return status from the correlation:

2 : Achieved desired precision in the norm of delta Phi

1 : Hit maximum number of iterations while iterating

-1 : Error is more than 80% of previous error, we’re probably diverging

-2 : Singular matrix M cannot be inverted

-3 : Displacement > 5*margin

’iterations’: int

Number of iterations

’logLikelyhood’float

Number indicating quality of match

’deltaPhiNorm’float

Size of last Phi step

’residualField’: 3D numpy array of floats

Same size as input image, residual field

’phaseField’: phaseField

Same size as input image, labelled phases

Return type:

Dictionary

Note

This correlation is what is proposed in Tudisco et al. “An extension of Digital Image Correlation for intermodality image registration”, section 4.3.

spam.DIC.multimodal.gaussianMixtureParameters(im1, im2, BINS=64, NPHASES=2, im1threshold=0, im2threshold=0, distanceMaxima=None, excludeBorder=False, fitDistance=None, GRAPHS=False, INTERACTIVE=False, sliceAxis=0, rootPath='.', suffix='')[source]#

This function, given two different modality images, builds the joint histogram in BINS bins, and fits NPHASES peaks with bivariate Gaussian distributions.

Parameters:

im1 (3D numpy array of floats) – One image, values should be rescaled to 0:BIN-1
im2 (3D numpy array of floats) – Another image, values should be rescaled to 0:BIN-1
BINS (int, optional) – Number of bins for the joint histogram, recommend 2^x Default = 64
NPHASES (int, optional) – Number of phases to automatically fit Default = 2
im1threshold (float, optional)
im2threshold (float, optional)
distanceMaxima (float, optional)
excludeBorder (False, optional) – Exclude edges in peak finding? This is passed to skimage.feature.peak_local_max()
fitDistance (float, optional)
GRAPHS (bool, optional)
INTERACTIVE (bool, optional)
sliceAxis=0
rootPath="."
suffix=""

Returns:

gaussianParameters – List of NPHASES components, containing the parameters of the bivariate Gaussian fit for each phase: [0] = GLOBALphi – Normlised “Height” of the fit [1] = GLOBALmu1 – Mean in F of bivairate Gaussian [2] = GLOBALmu2 – Mean in G of bivairate Gaussian [3] = a – [4] = b – [5] = c –

Return type:

list of lists

pBINSxBINS 2D numpy array of floats: The normalised joint histogram, the sum of which will be =1 if all of your image values are 0:BIN-1

spam.DIC.multimodal.phaseDiagram(gaussianParameters, jointHistogram, voxelCoverage=None, sigmaMax=9, BINS=64, GRAPHS=False, INTERACTIVE=False, suffix='', rootPath='.')[source]#: To be commented too

spam.DIC.pixelSearch module#

This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.

This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.

You should have received a copy of the GNU General Public License along with this program. If not, see <http://www.gnu.org/licenses/>.

spam.DIC.pixelSearch.pixelSearchDiscrete(lab1, im1, im2, PhiField, searchRange, boundingBoxes, numberOfNodes, nodePositions, applyF, labelDilate, volThreshold, numProc)[source]#

spam.DIC.pixelSearch.pixelSearchLocal(im1, im2, PhiField, searchRange, hws, im1mask, im2mask, numberOfNodes, nodePositions, applyF, maskCoverage, greyLowThresh, greyHighThresh, twoD, numProc)[source]#

spam.DIC.registration module#

spam.DIC.registration.register(im1, im2, im1mask=None, PhiInit=None, PhiRigid=False, PhiInitBinRatio=1.0, margin=None, maxIterations=25, deltaPhiMin=0.001, updateGradient=False, interpolationOrder=1, interpolator='python', verbose=False, returnPhiMaskCentre=True, imShowProgress=False, imShowProgressNewFig=False, imShowProgressPause=0.01)[source]#

Perform subpixel image correlation between im1 and im2.

The result of register(im1, im2) will give a deformation function $\Phi$ which maps im1 into im2. The Phi function used here allows the measurement of sub-pixel displacements, rotation, and linear straining of the whole image. However, this function will numerically deform im2 until it best matches im1.

$im1(x) = im2(\Phi.x)$

If im1 and im2 follow each other in time, then the resulting Phi is im1 -> im2 which makes sense in most cases. “Discrete correlation” can be performed by masking im1.

im1 and im2 do not necessarily have to be the same size (i.e., im2 can be bigger) – this is good since there is a zone to accommodate movement. In the case of a bigger im2, im1 and im2 are assumed to be centred with respect to each other.

Parameters:

im1 (3D numpy array) – The greyscale image that will not move – must not contain NaNs
im2 (3D numpy array) – The greyscale image that will be deformed – must not contain NaNs
im1mask (3D boolean numpy array, optional) – A mask for the zone to correlate in im1 with False in the zone to not correlate. Default = None, i.e., correlate all of im1 minus the margin. If this is defined, the Phi returned is in the centre of mass of the mask
PhiInit (4x4 numpy array, optional) – Initial deformation to apply to im1. Default = numpy.eye(4), i.e., no transformation
PhiRigid (bool, optional) – Run a rigid correlation? Only the rigid part of your PhiInit will be kept. Default = False
PhiInitBinRatio (float, optional) – Change translations in PhiInit, if it’s been calculated on a differently-binned image. Default = 1
margin (int, optional) – Margin, in pixels, to take in im1. Can also be a N-component list of ints, representing the margin in ND. If im2 has the same size as im1 this is strictly necessary to allow space for interpolation and movement Default = None (i.e., 10% of max dimension of im1)
maxIterations (int, optional) – Maximum number of quasi-Newton iterations to perform before stopping. Default = 25
deltaPhiMin (float, optional) – Smallest change in the norm of Phi (the transformation operator) before stopping. Default = 0.001
updateGradient (bool, optional) – Should the gradient of the image be computed (and updated) on the deforming im2? Default = False (it is computed once on im1)
interpolationOrder (int, optional) – Order of the greylevel interpolation for applying Phi to im1 when correlating. Recommended value is 3, but you can get away with 1 for faster calculations. Default = 3
interpolator (string, optional) – Which interpolation function to use from spam. Default = ‘python’. ‘C’ is also an option
verbose (bool, optional) – Get to know what the function is really thinking, recommended for debugging only. Default = False
returnPhiMaskCentre (bool, optional) – In case a mask is passed, should the Phi returned by the function be returned at the centre of mass of the mask? If False, it’s returned (as per no mask) in the centre of the image. Default = True
imShowProgress (bool, optional) – Pop up a window showing a imShowProgress slice of the image differences (im1-im2) as im1 is progressively deformed. Default = False
imShowProgressNewFig (bool, optional (defaul = False)) – Make a new plt.figure for each iteration, useful for examples gallery
imShowProgressPause (float, optional (defaul = 0.01)) – Seconds to pause between imShowProgress updates

Returns:

‘Phi’4x4 float array

Deformation function defined at the centre of the image

’returnStatus’signed int

Return status from the correlation:

2 : Achieved desired precision in the norm of delta Phi

1 : Hit maximum number of iterations while iterating

-1 : Error is more than 80% of previous error, we’re probably diverging

-2 : Singular matrix M (most probably image texture problem)

-3 : Displacement > 5*margin

-4 : Singular Phi matrix (most probably due to divergence)

-5 : Not used in this funciton but reserved for the script (correlation skipped)

’error’float

Error float describing mismatch between images, it’s the sum of the squared difference divided by the sum of im1

’iterations’int

Number of iterations

’deltaPhiNorm’float

Norm of deltaPhi

Return type:

Dictionary

Note

This correlation was written in the style of S. Roux (especially “An extension of Digital Image Correlation for intermodality image registration”) especially equations 12 and 13.

spam.DIC.registration.registerMultiscale(im1, im2, binStart, binStop=1, im1mask=None, PhiInit=None, PhiRigid=False, PhiInitBinRatio=1.0, margin=None, maxIterations=100, deltaPhiMin=0.0001, updateGradient=False, interpolationOrder=1, interpolator='C', verbose=False, returnPhiMaskCentre=True, imShowProgress=False, forceChangeScale=False)[source]#

Perform multiscale subpixel image correlation between im1 and im2.

This means applying a downscale (binning) to the images, performing a Lucas and Kanade at that level, and then improving it on a 2* less downscaled image, all the way back to the full scale image.

If your input images have multiple scales of texture, this should save significant time.

Please see the documentation for register for the rest of the documentation.

Parameters:

im1 (3D numpy array) – The greyscale image that will not move – must not contain NaNs
im2 (3D numpy array) – The greyscale image that will be deformed – must not contain NaNs
binStart (int) – Maximum amount of binning to apply, please input a number which is 2^int
binStop (int, optional) – Which binning level to stop upscaling at. The value of 1 (full image resolution) is almost always recommended (unless memory/time problems). Default = 1
im1mask (3D boolean numpy array, optional) – A mask for the zone to correlate in im1 with False in the zone to not correlate. Default = None, i.e., correlate all of im1 minus the margin. If this is defined, the Phi returned is in the centre of mass of the mask
PhiInit (4x4 numpy array, optional) – Initial deformation to apply to im1, by default at bin1 scale Default = numpy.eye(4), i.e., no transformation
PhiRigid (bool, optional) – Run a rigid correlation? Only the rigid part of your PhiInit will be kept. Default = False
PhiInitBinRatio (float, optional) – Change translations in PhiInit, if it’s been calculated on a differently-binned image. Default = 1
margin (int, optional) – Margin, in pixels, to take in im1. Can also be a N-component list of ints, representing the margin in ND. If im2 has the same size as im1 this is strictly necessary to allow space for interpolation and movement Default = 0 (i.e., 10% of max dimension of im1)
maxIterations (int, optional) – Maximum number of quasi-Newton iterations to perform before stopping. Default = 25
deltaPhiMin (float, optional) – Smallest change in the norm of Phi (the transformation operator) before stopping. Default = 0.001
updateGradient (bool, optional) – Should the gradient of the image be computed (and updated) on the deforming im2? Default = False (it is computed once on im1)
interpolationOrder (int, optional) – Order of the greylevel interpolation for applying Phi to im1 when correlating. Recommended value is 3, but you can get away with 1 for faster calculations. Default = 3
interpolator (string, optional) – Which interpolation function to use from spam. Default = ‘python’. ‘C’ is also an option
verbose (bool, optional) – Get to know what the function is really thinking, recommended for debugging only. Default = False
returnPhiMaskCentre (bool, optional) – In case a mask is passed, should the Phi returned by the function be returned at the centre of mass of the mask? If False, it’s returned (as per no mask) in the centre of the image. Default = True
imShowProgress (bool, optional) – Pop up a window showing a imShowProgress slice of the image differences (im1-im2) as im1 is progressively deformed. Default = False
forceChangeScale (bool, optional) – Change up a scale even if not converged? Default = False

Returns:

‘Phi’: 4x4 float array

Deformation function defined at the centre of the image

’returnStatus’: signed int

Return status from the correlation:

2 : Achieved desired precision in the norm of delta Phi

1 : Hit maximum number of iterations while iterating

-1 : Error is more than 80% of previous error, we’re probably diverging

-2 : Singular matrix M (most probably image texture problem)

-3 : Displacement > 5*margin

-4 : Singular Phi matrix (most probably due to divergence)

-5 : Not used in this funciton but reserved for the script (correlation skipped)

’error’: float

Error float describing mismatch between images, it’s the sum of the squared difference divided by the sum of im1

’iterations’: int

Number of iterations

Return type:

Dictionary

spam.DIC package

Contents

spam.DIC package#

Submodules#

spam.DIC.ddic module#

spam.DIC.deform module#

spam.DIC.globalDVC module#

spam.DIC.grid module#

spam.DIC.kinematics module#

spam.DIC.ldic module#

spam.DIC.multimodal module#

spam.DIC.pixelSearch module#

spam.DIC.registration module#

Module contents#