Commit c0206920 authored by Olga Stamati's avatar Olga Stamati
Browse files

increase test coverage in grid

parent 429c8b3c
Pipeline #50500 passed with stages
in 23 minutes and 55 seconds
......@@ -407,7 +407,7 @@ class TestFunctionDVC(unittest.TestCase):
nodePositions = spam.DIC.grid.makeGrid(im.shape, 10)[0]
Phi = spam.deformation.computePhi({'t': [0, 2, 1]})
imDef = spam.DIC.applyPhi(im, Phi=Phi)
ps = spam.DIC.grid.pixelSearchOnGrid(im, imDef,
ps = spam.DIC.pixelSearchOnGrid(im, imDef,
nodePositions, 5,
searchRange={"zRange": [0, 0], "yRange": [-3, 3], "xRange": [-3, 3]})
......@@ -415,13 +415,23 @@ class TestFunctionDVC(unittest.TestCase):
self.assertEqual(ps["PhiField"][:, 1, -1].sum(), nodePositions.shape[0] * 2)
self.assertEqual(ps["PhiField"][:, 2, -1].sum(), nodePositions.shape[0])
# force to fail for greyThreshold condition
psb = spam.DIC.pixelSearchOnGrid(im, imDef,
nodePositions, 5,
searchRange={"zRange": [0, 0], "yRange": [-3, 3], "xRange": [-3, 3]},
greyThreshold=[numpy.inf, numpy.inf])
for node in range(nodePositions.shape[0]):
self.assertTrue(numpy.allclose(psb["PhiField"][node, 0:3, 0:3], numpy.eye(3)))
self.assertEqual(numpy.isfinite(psb["PhiField"][:, 0:3, -1]).sum(), 0)
# pass default optional parameters
ps2 = spam.DIC.grid.pixelSearchOnGrid(im, imDef,
nodePositions, 5,
searchRange={"zRange": [0, 0], "yRange": [-3, 3], "xRange": [-3, 3]},
im1mask=numpy.ones_like(im),
minMaskCoverage=0.5,
greyThreshold=[-numpy.inf, numpy.inf])
ps2 = spam.DIC.pixelSearchOnGrid(im, imDef,
nodePositions, 5,
searchRange={"zRange": [0, 0], "yRange": [-3, 3], "xRange": [-3, 3]},
im1mask=numpy.ones_like(im),
minMaskCoverage=0.5,
greyThreshold=[-numpy.inf, numpy.inf])
self.assertEqual(ps2["PhiField"][:, 0, -1].sum(), 0)
self.assertEqual(ps2["PhiField"][:, 1, -1].sum(), nodePositions.shape[0] * 2)
......@@ -441,13 +451,14 @@ class TestFunctionDVC(unittest.TestCase):
PhiField[nodeNumber, 0:3, 0:3] = Phi[0:3, 0:3]
# t
PhiField[nodeNumber, 0:3, -1 ] = spam.deformation.decomposePhi(Phi.copy(), PhiCentre=imCentre, PhiPoint=nodePosition)['t']
ps2 = spam.DIC.grid.pixelSearchOnGrid(im, imDef,
nodePositions, 5,
searchRange={"zRange": [-3, 3], "yRange": [-3, 3], "xRange": [-3, 3]},
PhiField=PhiField)
ps3 = spam.DIC.pixelSearchOnGrid(im, imDef,
nodePositions, 5,
searchRange={"zRange": [-3, 3], "yRange": [-3, 3], "xRange": [-3, 3]},
PhiField=PhiField)
for nodeNumber, nodePosition in enumerate(nodePositions):
self.assertAlmostEqual(numpy.sum(PhiField[nodeNumber] - ps2["PhiField"][nodeNumber]), 0, places=1)
self.assertAlmostEqual(numpy.sum(PhiField[nodeNumber] - ps3["PhiField"][nodeNumber]), 0, places=1)
def test_registerOnGrid(self):
im = spam.datasets.loadSnow()
......@@ -459,35 +470,65 @@ class TestFunctionDVC(unittest.TestCase):
PhiField[:, 1, -1] = 2
PhiField[:, 2, -1] = 1
subPix = spam.DIC.grid.registerOnGrid(im, imDef,
nodePositions, 10,
PhiField=PhiField,
margin=1,
maxIterations=100,
deltaPhiMin=0.001,
interpolationOrder=1,
minMaskCoverage=0.5, im1mask=numpy.ones_like(im),
greyThreshold=[-numpy.inf, numpy.inf])
subPix = spam.DIC.registerOnGrid(im, imDef,
nodePositions, 10,
PhiField=PhiField,
margin=1,
maxIterations=100,
deltaPhiMin=0.001,
interpolationOrder=1,
minMaskCoverage=0.5, im1mask=numpy.ones_like(im),
greyThreshold=[-numpy.inf, numpy.inf])
self.assertEqual(subPix["PhiField"][:, 0, -1].sum(), 0)
self.assertEqual(subPix["PhiField"][:, 1, -1].sum(), nodePositions.shape[0] * 2)
self.assertEqual(subPix["PhiField"][:, 2, -1].sum(), nodePositions.shape[0])
self.assertEqual(subPix["returnStatus"].sum(), nodePositions.shape[0] * 2)
# check 2d + threshold + python interpolator
nodePositions[:, 0] = 1
subPix2 = spam.DIC.registerOnGrid(im[25], imDef[25],
nodePositions, 10,
interpolationOrder=3,
greyThreshold=[numpy.inf, numpy.inf])
self.assertEqual(numpy.isfinite(subPix2["PhiField"][:, 0:3, -1]).sum(), 0)
for node in range(nodePositions.shape[0]):
self.assertTrue(numpy.allclose(subPix2["PhiField"][node, 0:3, 0:3], numpy.eye(3)))
self.assertEqual(subPix2["returnStatus"][node], -5)
# check nans in PhiField
PhiField[:, 0:3, -1] *= numpy.nan
subPix3 = spam.DIC.registerOnGrid(im, imDef,
nodePositions, 10,
PhiField=PhiField,
margin=1,
maxIterations=100,
deltaPhiMin=0.001,
interpolationOrder=1,
minMaskCoverage=0.5, im1mask=numpy.ones_like(im),
greyThreshold=[-numpy.inf, numpy.inf])
self.assertEqual(numpy.isfinite(subPix3["PhiField"][:, 0:3, -1]).sum(), 0)
for node in range(nodePositions.shape[0]):
self.assertTrue(numpy.allclose(subPix3["PhiField"][node, 0:3, 0:3], numpy.eye(3)))
self.assertEqual(subPix3["returnStatus"][node], -6)
def test_grid(self):
im = spam.datasets.loadSnow()
imageSize = im.shape
# Make sure it fails with 2D image
grid = spam.DIC.makeGrid(im[im.shape[0] // 2], 10)
self.assertTrue(grid is None)
grid1 = spam.DIC.makeGrid(imageSize[1:], 10)
self.assertTrue(grid1 is None)
# Make sure it fails with wrong number of node spacings
grid = spam.DIC.makeGrid(im, [10, 10])
self.assertTrue(grid is None)
# Make sure it fails with wrong number of node spacings
grid = spam.DIC.makeGrid(im, [10, 10, 10, 10])
self.assertTrue(grid is None)
# Make sure it fails with wrong number of node spacing
grid2 = spam.DIC.makeGrid(imageSize, [10, 10])
self.assertTrue(grid2 is None)
# Make sure it fails with wrong number of node spacing
grid3 = spam.DIC.makeGrid(imageSize, [10, 10, 10, 10])
self.assertTrue(grid3 is None)
if __name__ == '__main__':
unittest.main()
......@@ -349,7 +349,7 @@ def pixelSearchOnGrid(im1, im2, nodePositions, halfWindowSize, searchRange, PhiF
'pixelSearchCC': pixelSearchCC}
def registerOnGrid(im1, im2, nodePositions, halfWindowSize, PhiField=None, margin=None, maxIterations=None, deltaPhiMin=None, interpolationOrder=None, minMaskCoverage=None, im1mask=None, greyThreshold=[-numpy.inf, numpy.inf], mpi=False, killWorkersWhenDone=True):
def registerOnGrid(im1, im2, nodePositions, halfWindowSize, PhiField=None, margin=None, maxIterations=None, deltaPhiMin=None, interpolationOrder=None, minMaskCoverage=0.5, im1mask=None, greyThreshold=[-numpy.inf, numpy.inf], mpi=False, killWorkersWhenDone=True):
"""
This function handles grid-based local correlation, performing a "register" subpixel refinement.
Here we minimise a residual which is the difference between im1 and im2.
......
......@@ -394,120 +394,120 @@ def projectTetFieldToGrains(points, connectivity, tetField):
return(grainField)
def projectBagiStrainToGrid(points, connectivity, tet_strain, nx=8, ny=8, nz=10):
"""
Projects/coarse-grains strains from tetrahedra onto a grid. The first step is
constructing a Cartesian grid based on the extremal positions of all tetrahedra
nodes. The next step is finding all tetrahedra having a node within each grid cell.
the final step is volume-averaging over all of the corresponding tetrahedra.
Parameters
----------
points: m x 3 numpy array
M Particles' coordinates in deformed configuration
connectivity: n x 4 numpy array
Delaunay triangulation connectivity generated by spam.mesh.delaunayFromVoronoi for example
tet_strain: n x 3 x 3 numpy array
Tetrahedra strail from spam.mesh.bagiStrain. This can be
either tetStrainsSym or tetStrainAsym. The resulting projected/coarse-grained
will be either the symmetric or skew-symmetric strain, respectively
Returns
-------
gridStrain:
Ng array containing (3x3) arrays of strain. (Ng = number of grains)
Example
-------
grid_bounds,num_in_grid,gridStrain
= spam.mesh.projectBagiStrainToGrid(connectivity,tet_strain,points,nx,ny,nz)
Returns strains for each grid zone and boundaries of grid zone.
"""
# gridStrainVoigt: Ng array containing (6x1) arrays of strain in Voigt notation.
# RCH Oct 2018.
# Import modules
import numpy
#import scipy.spatial as sps
from progressbar import ProgressBar
pbar = ProgressBar()
tetVolumes2 = spam.mesh.tetVolumes(points, connectivity)
#def _projectBagiStrainToGrid(points, connectivity, tet_strain, nx=8, ny=8, nz=10):
#"""
#Projects/coarse-grains strains from tetrahedra onto a grid. The first step is
#constructing a Cartesian grid based on the extremal positions of all tetrahedra
#nodes. The next step is finding all tetrahedra having a node within each grid cell.
#the final step is volume-averaging over all of the corresponding tetrahedra.
#Parameters
#----------
#points: m x 3 numpy array
#M Particles' coordinates in deformed configuration
#connectivity: n x 4 numpy array
#Delaunay triangulation connectivity generated by spam.mesh.delaunayFromVoronoi for example
#tet_strain: n x 3 x 3 numpy array
#Tetrahedra strail from spam.mesh.bagiStrain. This can be
#either tetStrainsSym or tetStrainAsym. The resulting projected/coarse-grained
#will be either the symmetric or skew-symmetric strain, respectively
#Returns
#-------
#gridStrain:
#Ng array containing (3x3) arrays of strain. (Ng = number of grains)
#Example
#-------
#grid_bounds,num_in_grid,gridStrain
#= spam.mesh.projectBagiStrainToGrid(connectivity,tet_strain,points,nx,ny,nz)
#Returns strains for each grid zone and boundaries of grid zone.
#"""
## gridStrainVoigt: Ng array containing (6x1) arrays of strain in Voigt notation.
## RCH Oct 2018.
## Import modules
#import numpy
##import scipy.spatial as sps
#from progressbar import ProgressBar
#pbar = ProgressBar()
# Artificial edge buffer to make sure we grab the outer-most particles
edgeBuf = 0.01
# Get sample extents for grid
edges = numpy.array([[points[:, 0].min() - edgeBuf, points[:, 0].max() + edgeBuf],
[points[:, 1].min() - edgeBuf, points[:, 1].max() + edgeBuf],
[points[:, 2].min() - edgeBuf, points[:, 2].max() + edgeBuf]])
# Construct grid with convention xmin,xmax,ymin,ymax,zmin,zmax
grid_bounds = []
xspace = (edges[0, 1] - edges[0, 0]) / (nx)
yspace = (edges[1, 1] - edges[1, 0]) / (ny)
zspace = (edges[2, 1] - edges[2, 0]) / (nz)
for ii in range(nz):
for jj in range(nx):
for kk in range(ny):
grid_bounds.append(numpy.array([edges[0, 0] + jj * xspace, edges[0, 0] + (jj + 1) * xspace,
edges[1, 0] + kk * yspace, edges[1, 0] + (kk + 1) * yspace,
edges[2, 0] + ii * zspace, edges[2, 0] + (ii + 1) * zspace]))
grid_bounds = numpy.array(grid_bounds)
# Initialize list of tet strains
gridStrain = []
# gridStrainVoigt = []
# First, find all grains in each grid:
print("spam.mesh.projectBagiStrainToGrid(): Projecting tetrahedal field to grid...\n")
for ii in pbar(range(grid_bounds.shape[0])):
xcond = numpy.logical_and(points[:, 0] > grid_bounds[ii, 0], points[:, 0] <= grid_bounds[ii, 1])
ycond = numpy.logical_and(points[:, 1] > grid_bounds[ii, 2], points[:, 1] <= grid_bounds[ii, 3])
zcond = numpy.logical_and(points[:, 2] > grid_bounds[ii, 4], points[:, 2] <= grid_bounds[ii, 5])
indgrain = numpy.where(numpy.logical_and(xcond, numpy.logical_and(ycond, zcond)))
# Find all tets for which these grains act as a node:
indtet = []
for jj in range(indgrain[0].shape[0]):
indtettmp, _ = numpy.where(connectivity == indgrain[0][jj])
for kk in range(indtettmp.shape[0]):
indtet.append(indtettmp[kk])
# Find unique values of these tets
indtet = numpy.unique(indtet, return_index=False)
# Loop over unique tets to get volume-averaged strains
volTot = 0.0
strainTot = numpy.zeros((3, 3))
num_in_grid = []
num_in_grid_tmp = 0
if (indtet.shape[0] > 0):
for jj in range(indtet.shape[0]):
# Get volume
tet_ids = connectivity[indtet[jj], :]
tetCoords = points[tet_ids, :]
#cvxh = sps.ConvexHull(tetCoords)
vol = tetVolumes2[indtet[jj]]
# Add volume and strain
strainTot = strainTot + tet_strain[indtet[jj]] * vol
volTot += vol
# Divide strainTot by volTot
strainTot = strainTot / volTot
num_in_grid_tmp = num_in_grid_tmp + indtet.shape[0]
# Store in grid
num_in_grid.append(num_in_grid_tmp)
gridStrain.append(numpy.array(strainTot))
# gridStrainVoigt.append(numpy.array([strainTot[0, 0], strainTot[1, 1], strainTot[2, 2],
# strainTot[1, 2], strainTot[0, 2], strainTot[0, 1]]))
# return(grid_bounds, num_in_grid, gridStrain, gridStrainVoigt)
return(grid_bounds, num_in_grid, numpy.array(gridStrain).reshape(nz, ny, nx))
#tetVolumes2 = spam.mesh.tetVolumes(points, connectivity)
## Artificial edge buffer to make sure we grab the outer-most particles
#edgeBuf = 0.01
## Get sample extents for grid
#edges = numpy.array([[points[:, 0].min() - edgeBuf, points[:, 0].max() + edgeBuf],
#[points[:, 1].min() - edgeBuf, points[:, 1].max() + edgeBuf],
#[points[:, 2].min() - edgeBuf, points[:, 2].max() + edgeBuf]])
## Construct grid with convention xmin,xmax,ymin,ymax,zmin,zmax
#grid_bounds = []
#xspace = (edges[0, 1] - edges[0, 0]) / (nx)
#yspace = (edges[1, 1] - edges[1, 0]) / (ny)
#zspace = (edges[2, 1] - edges[2, 0]) / (nz)
#for ii in range(nz):
#for jj in range(nx):
#for kk in range(ny):
#grid_bounds.append(numpy.array([edges[0, 0] + jj * xspace, edges[0, 0] + (jj + 1) * xspace,
#edges[1, 0] + kk * yspace, edges[1, 0] + (kk + 1) * yspace,
#edges[2, 0] + ii * zspace, edges[2, 0] + (ii + 1) * zspace]))
#grid_bounds = numpy.array(grid_bounds)
## Initialize list of tet strains
#gridStrain = []
## gridStrainVoigt = []
## First, find all grains in each grid:
#print("spam.mesh.projectBagiStrainToGrid(): Projecting tetrahedal field to grid...\n")
#for ii in pbar(range(grid_bounds.shape[0])):
#xcond = numpy.logical_and(points[:, 0] > grid_bounds[ii, 0], points[:, 0] <= grid_bounds[ii, 1])
#ycond = numpy.logical_and(points[:, 1] > grid_bounds[ii, 2], points[:, 1] <= grid_bounds[ii, 3])
#zcond = numpy.logical_and(points[:, 2] > grid_bounds[ii, 4], points[:, 2] <= grid_bounds[ii, 5])
#indgrain = numpy.where(numpy.logical_and(xcond, numpy.logical_and(ycond, zcond)))
## Find all tets for which these grains act as a node:
#indtet = []
#for jj in range(indgrain[0].shape[0]):
#indtettmp, _ = numpy.where(connectivity == indgrain[0][jj])
#for kk in range(indtettmp.shape[0]):
#indtet.append(indtettmp[kk])
## Find unique values of these tets
#indtet = numpy.unique(indtet, return_index=False)
## Loop over unique tets to get volume-averaged strains
#volTot = 0.0
#strainTot = numpy.zeros((3, 3))
#num_in_grid = []
#num_in_grid_tmp = 0
#if (indtet.shape[0] > 0):
#for jj in range(indtet.shape[0]):
## Get volume
#tet_ids = connectivity[indtet[jj], :]
#tetCoords = points[tet_ids, :]
##cvxh = sps.ConvexHull(tetCoords)
#vol = tetVolumes2[indtet[jj]]
## Add volume and strain
#strainTot = strainTot + tet_strain[indtet[jj]] * vol
#volTot += vol
## Divide strainTot by volTot
#strainTot = strainTot / volTot
#num_in_grid_tmp = num_in_grid_tmp + indtet.shape[0]
## Store in grid
#num_in_grid.append(num_in_grid_tmp)
#gridStrain.append(numpy.array(strainTot))
## gridStrainVoigt.append(numpy.array([strainTot[0, 0], strainTot[1, 1], strainTot[2, 2],
## strainTot[1, 2], strainTot[0, 2], strainTot[0, 1]]))
## return(grid_bounds, num_in_grid, gridStrain, gridStrainVoigt)
#return(grid_bounds, num_in_grid, numpy.array(gridStrain).reshape(nz, ny, nx))
def BCFieldFromDVCField(points, dvcField, mask=None, pixelSize=1.0, meshType="cube", centre=[0, 0], radius=1.0, topBottom=False, tol=1e-6, neighbours=4):
......
Supports Markdown
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