Commit 24bc5f77 authored by Edward Andò's avatar Edward Andò
Browse files

start of integration of spam-pixelSearch into tests (warning commented tests in test_scripts.py)

parent 770b7f0f
Pipeline #56031 failed with stages
in 11 minutes and 59 seconds
......@@ -285,7 +285,7 @@ def getImagettes(nodePosition, Phi, searchRange, boundingBox, im1, im2, im1mask,
imagette1mask = None
else:
imagette1mask = spam.helpers.slicePadded(im1mask, boundingBox + numpy.array([0, 1, 0, 1, 0, 1])) != 0
maskVolumeCondition = imagette1mask.mean() > minMaskCoverage
maskVolumeCondition = imagette1mask.mean() >= minMaskCoverage
# Make sure imagette is not 0-dimensional in any dimension
......@@ -339,7 +339,12 @@ def getImagettes(nodePosition, Phi, searchRange, boundingBox, im1, im2, im1mask,
# Create pixelSearchCC vector
pixelSearchCC = numpy.zeros((numberOfNodes), dtype=float)
# Error compatible with register()
error = numpy.zeros((numberOfNodes), dtype=float)
returnStatus = numpy.ones((numberOfNodes), dtype=int)
deltaPhiNorm = numpy.ones((numberOfNodes), dtype=int)
iterations = numpy.ones((numberOfNodes), dtype=int)
# Add nodes to a queue -- mostly useful for MPI
q = queue.Queue()
......@@ -408,7 +413,8 @@ while finishedNodes != numberOfNodes:
if imagetteReturns['returnStatus'] == 1:
returns = spam.DIC.correlate.pixelSearch(imagetteReturns['imagette1'],
imagetteReturns['imagette2'],
imagette1mask = imagetteReturns['imagette1mask'])
imagette1mask = imagetteReturns['imagette1mask'],
returnError = True)
pixelSearchOffset = imagetteReturns['pixelSearchOffset']
writeReturns = True
......@@ -432,6 +438,7 @@ while finishedNodes != numberOfNodes:
PhiField[nodeNumber, 0:3, 3] += numpy.array(returns[0]) + pixelSearchOffset
pixelSearchCC[nodeNumber] = returns[1]
error[nodeNumber] = returns[2]
widgets[0] = progressbar.FormatLabel(" CC={:0>7.5f} ".format(pixelSearchCC[nodeNumber]))
pbar.update(finishedNodes)
......@@ -446,15 +453,19 @@ if args.TSV:
# F[0:3,0:2]
# Pixel-search CC
if args.LAB1 is not None:
TSVheader = "Label\tZpos\tYpos\tXpos\tFzz\tFzy\tFzx\tZdisp\tFyz\tFyy\tFyx\tYdisp\tFxz\tFxy\tFxx\tXdisp\tpixelSearchCC\treturnStatus"
TSVheader = "Label\tZpos\tYpos\tXpos\tFzz\tFzy\tFzx\tZdisp\tFyz\tFyy\tFyx\tYdisp\tFxz\tFxy\tFxx\tXdisp\tpixelSearchCC\treturnStatus\terror\tdeltaPhiNorm\titerations"
else:
TSVheader = "NodeNumber\tZpos\tYpos\tXpos\tFzz\tFzy\tFzx\tZdisp\tFyz\tFyy\tFyx\tYdisp\tFxz\tFxy\tFxx\tXdisp\tpixelSearchCC\treturnStatus"
TSVheader = "NodeNumber\tZpos\tYpos\tXpos\tFzz\tFzy\tFzx\tZdisp\tFyz\tFyy\tFyx\tYdisp\tFxz\tFxy\tFxx\tXdisp\tpixelSearchCC\treturnStatus\terror\tdeltaPhiNorm\titerations"
outMatrix = numpy.array([numpy.array(range(nodePositions.shape[0])),
nodePositions[:, 0], nodePositions[:, 1], nodePositions[:, 2],
nodePositions[:, 0], nodePositions[:, 1], nodePositions[:, 2],
PhiField[:, 0, 0], PhiField[:, 0, 1], PhiField[:, 0, 2], PhiField[:, 0, 3],
PhiField[:, 1, 0], PhiField[:, 1, 1], PhiField[:, 1, 2], PhiField[:, 1, 3],
PhiField[:, 2, 0], PhiField[:, 2, 1], PhiField[:, 2, 2], PhiField[:, 2, 3],
pixelSearchCC, returnStatus]).T
pixelSearchCC,
returnStatus,
error,
deltaPhiNorm,
iterations]).T
numpy.savetxt(args.OUT_DIR+"/"+args.PREFIX+".tsv",
outMatrix,
......
......@@ -20,8 +20,8 @@ import coverage
FNULL = open(os.devnull, 'w')
# We also check 2D slice so don't displace in Z or rotate around an axis other than Z
refTranslation = [0.0, 2.0, 0.0]
refRotation = [5.0, 0.0, 0.0]
refTranslation = numpy.array([0.0, 15.0, 0.0])
refRotation = numpy.array([5.0, 0.0, 0.0])
testFolder = "./"
try:
......@@ -42,6 +42,7 @@ class testAll(unittest.TestCase):
try:
rm(testFolder+"snow-ref.tif")
rm(testFolder+"snow-def.tif")
rm(testFolder+"snow-mask.tif")
rm(testFolder+"snow-def-cgs.tif")
rm(testFolder+"snow-ref-lab-displaced.tif")
rm(testFolder+"snow-ref-lab.tif")
......@@ -88,19 +89,20 @@ class testAll(unittest.TestCase):
pass
def test_ldic_and_regularStrain(self):
def test_gridDICtools(self):
# Make test directory
if not os.path.isdir(testFolder):
os.mkdir(testFolder)
# We'll do def -> ref to avoid interpolation errors
#######################################################
### Step 1 generate files
#######################################################
# Run on snow data, with bin2 registration, without output
# load 3D image
snowRef = spam.datasets.loadSnow()
# save it locally
tifffile.imsave(testFolder + "snow-ref.tif", snowRef)
os.chmod(testFolder + "snow-ref.tif", 0o666)
os.chmod( testFolder + "snow-ref.tif", 0o666)
# Generate deformed image
Phi = spam.deformation.computePhi({'t': refTranslation,
'r': refRotation})
......@@ -109,20 +111,53 @@ class testAll(unittest.TestCase):
tifffile.imsave(testFolder + "snow-def.tif", snowDef)
os.chmod(testFolder + "snow-def.tif", 0o666)
# make sure it runs the help without error
exitCode = subprocess.call(["spam-ldic", "--help"], stdout=FNULL, stderr=FNULL)
# Mask for pixel search
snowMask = numpy.ones_like(snowDef, dtype='<u1')
snowMask[0] = 0
snowMask[-1] = 0
tifffile.imsave(testFolder + "snow-mask.tif", snowDef)
os.chmod(testFolder + "snow-mask.tif", 0o666)
# generate a fake "eye reg" initial guess which is close
spam.helpers.writeRegistrationTSV(testFolder + "snow-ref-snow-def-ereg.tsv", (numpy.array( snowRef.shape )-1)/2.0,
{'Phi': spam.deformation.computePhi({'t': 0.8*refTranslation, 'r': 0.8*refRotation}),
'error': 0,
'iterations': 0,
'returnStatus': 2,
'deltaPhiNorm': 1})
#######################################################
### Step 2 check spam-reg functionality
#######################################################
## Just run a simple registration
exitCode = subprocess.call(["spam-reg", "-bb", "2", testFolder + "snow-ref.tif", testFolder + "snow-def.tif", "-od", testFolder + "", "-mar", "20", "-pf", testFolder + "snow-ref-snow-def-ereg.tsv", "-mf1", testFolder + "snow-mask.tif"])
self.assertEqual(exitCode, 0)
regResult = spam.helpers.readCorrelationTSV(testFolder + "snow-ref-snow-def-registration.tsv")
transformation = spam.deformation.decomposePhi(regResult['PhiField'][0])
print(transformation['t'])
print(transformation['r'])
self.assertTrue(numpy.allclose(refTranslation, transformation['t'], atol=0.01))
self.assertTrue(numpy.allclose(refRotation, transformation['r'], atol=0.01))
#######################################################
### Step 3 check spam-pixelSearch functionality
#######################################################
### Step 3.1 load initial guess and just search around it...
exitCode = subprocess.call(["spam-pixelSearch", "-pf", testFolder + "snow-ref-snow-def-registration.tsv", "-sr", "-2", "2", "-2", "2", "-2", "2", "-glt", "5000", "-hws", "10", testFolder + "snow-ref.tif", testFolder + "snow-def.tif", "-od", testFolder + "", "-tif", "-mf1", testFolder + "snow-mask.tif", "-mc", "1.0"])
self.assertEqual(exitCode, 0)
PSresult = spam.helpers.readCorrelationTSV(testFolder + "snow-ref-snow-def-pixelSearch.tsv", readPixelSearchCC=True)
print(PSresult['PhiField'][PSresult['returnStatus']==1,0,-1].mean())
print(PSresult['pixelSearchCC'][PSresult['returnStatus']==1].mean())
# Assert that the CC is nice and high
self.assertTrue(0.95 < PSresult['pixelSearchCC'][PSresult['returnStatus']==1].mean())
# And the z-displacement is low
self.assertTrue(numpy.isclose(0, PSresult['PhiField'][PSresult['returnStatus']==1,0,-1].mean(), atol=1.0))
## Just run a simple DVC with no outputs
# exitCode = subprocess.call(["spam-ldic", "-reg", "-regbb", "2", "-hws", "10", testFolder+"snow-ref.tif", testFolder+"snow-def.tif", "-od", testFolder+""], stdout=FNULL, stderr=FNULL)
exitCode = subprocess.call(["spam-ldic", "-reg", "-regbb", "2", "-hws", "10", testFolder + "snow-ref.tif", testFolder + "snow-def.tif", "-od", testFolder + ""])
#self.assertEqual(exitCode, 0)
# Check output results with TSV output, bin1 registration + TIFFs
# Decreasing node spacing in order to have more points for Geers strain calculation
# exitCode = subprocess.call(["spam-ldic", "-Ffile", testFolder+"snow-ref-snow-def-bin2-registration.tsv", "-Ffb", "2", "-hws", "15", testFolder+"snow-ref.tif", testFolder+"snow-def.tif", "-od", testFolder+"", "-tsv"])
exitCode = subprocess.call(["spam-ldic", "-ps", "off", "-reg", "-regbb", "2", "-regbe", "1", "-regm", "0.1", "-regu", "-gpi", "10", "-hws", "10", '-ns', '15', testFolder + "snow-ref.tif", testFolder + "snow-def.tif", "-od", testFolder + "", "-tsv", "-tif", '-vtk'])
self.assertEqual(exitCode, 0)
# Check registration
#######################################################
......@@ -214,7 +249,7 @@ class testAll(unittest.TestCase):
self.assertAlmostEqual(numpy.mean(TSV['PhiField'][:,1,-1][TSV['returnStatus']==2]), refTranslation[1], places=2)
self.assertAlmostEqual(numpy.mean(TSV['PhiField'][:,2,-1][TSV['returnStatus']==2]), refTranslation[2], places=2)
def test_ITKwatershed(self):
def _test_ITKwatershed(self):
# make sure it runs the help without error
exitCode = subprocess.call(["spam-ITKwatershed", "--help"],
stdout=FNULL, stderr=FNULL)
......@@ -236,12 +271,12 @@ class testAll(unittest.TestCase):
lab = tifffile.imread(testFolder + "snow-ref-watershed.tif")
self.assertEqual(lab.max()>10, True)
def test_gdic(self):
def _test_gdic(self):
# make sure it runs the help without error
exitCode = subprocess.call(["spam-gdic", "--help"], stdout=FNULL, stderr=FNULL)
self.assertEqual(exitCode, 0)
def test_deformImage(self):
def _test_deformImage(self):
# make sure it runs the help without error
exitCode = subprocess.call(["spam-deformImageFromField", "--help"], stdout=FNULL, stderr=FNULL)
self.assertEqual(exitCode, 0)
......@@ -366,7 +401,7 @@ class testAll(unittest.TestCase):
self.assertEqual(numpy.mean(imDefScript[10:-10, 10:-10, 10:-10]) > numpy.mean(imDefScriptCGS[10:-10, 10:-10, 10:-10]), True)
def test_ddic_and_grainMover(self):
def _test_ddic_and_grainMover(self):
#######################################################
### 1. Generate Data
#######################################################
......@@ -713,7 +748,7 @@ class testAll(unittest.TestCase):
for label in [1,2]:
self.assertAlmostEqual(numpy.abs(TSVextreme['PhiField'][label, axis, -1] - 4.0), 0.0, places=1)
def test_discreteStrain(self):
def _test_discreteStrain(self):
# make sure it runs the help without error
exitCode = subprocess.call(["spam-discreteStrain", "--help"],
stdout=FNULL,
......@@ -746,6 +781,18 @@ class testAll(unittest.TestCase):
#self.assertAlmostEqual(VTK[3]['yy'].mean(), (1.01 - 1.00), places=3)
#self.assertAlmostEqual(VTK[3]['zz'].mean(), (1.01 - 1.00), places=3)
def test_all_help(self):
exitCode = subprocess.call(["spam-reg", "--help"], stdout=FNULL, stderr=FNULL)
self.assertEqual(exitCode, 0)
exitCode = subprocess.call(["spam-pixelSearch", "--help"], stdout=FNULL, stderr=FNULL)
self.assertEqual(exitCode, 0)
exitCode = subprocess.call(["spam-ldic", "--help"], stdout=FNULL, stderr=FNULL)
self.assertEqual(exitCode, 0)
if __name__ == '__main__':
unittest.main()
......@@ -32,6 +32,9 @@ import progressbar
#numpy.set_printoptions(precision=3, suppress=True)
def _errorCalc(im1, im2):
return numpy.nansum( numpy.square(numpy.subtract(im1, im2)) ) / numpy.nansum(im1)
def 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, imShowProgress=False, imShowProgressNewFig=False):
"""
Perform subpixel image correlation between im1 and im2.
......@@ -257,20 +260,6 @@ def register(im1, im2, im1mask=None, PhiInit=None, PhiRigid=False, PhiInitBinRat
plt.axis([im1crop.shape[1], 0, im1crop.shape[0], 0])
plt.ion()
# Numerical value for normalising the error
# errorNormalisationTmp = im1[crop1]
# errorNormalisation = errorNormalisationTmp[numpy.isfinite(errorNormalisationTmp)].sum()
# del errorNormalisationTmp
# 2018-07-10 EA and OS: Removing im1 mask, so errorNormalisation is one-shot
# 2019-01-18 EA putting it back in so error is more meaningful
if im1mask is not None:
errorNormalisation = im1crop[numpy.where(im1mask[crop1] == True)].sum()
else:
errorNormalisation = im1crop.sum()
# 2019-09-11 EA: Avoid nans
if errorNormalisation == 0: errorNormalisation = 1
###########################################################
# Important -- since we're moving im2, initial Phis will be
# pointing the wrong way, they need to be inversed
......@@ -330,8 +319,7 @@ def register(im1, im2, im1mask=None, PhiInit=None, PhiRigid=False, PhiInitBinRat
returnStatus = 0
# Big value to start with to ensure the first iteration
deltaPhiNorm = 100.0
errorTmp = numpy.square(numpy.subtract(im1crop, im2def[crop2]))
error = errorTmp[numpy.isfinite(errorTmp)].sum() / errorNormalisation
error = _errorCalc(im1crop, im2def[crop2])
if verbose:
print("Start correlation with Error = {:0.2f}".format(error))
......@@ -435,8 +423,7 @@ def register(im1, im2, im1mask=None, PhiInit=None, PhiRigid=False, PhiInitBinRat
elif interpolator == 'python': im2def = spam.DIC.applyPhiPython(im2, Phi=PhiInv, interpolationOrder=interpolationOrder)
# Error calculation
errorTmp = numpy.square(numpy.subtract(im1crop, im2def[crop2]))
error = errorTmp[numpy.isfinite(errorTmp)].sum() / errorNormalisation
error = _errorCalc(im1crop, im2def[crop2])
# Keep interested people up to date with what's happening
# if verbose:
......@@ -762,7 +749,7 @@ def registerMultiscale(im1, im2, binStart, binStop=1, im1mask=None, PhiInit=None
return reg
def pixelSearch(imagette1, imagette2, imagette1mask=None, imagette2mask=None):
def pixelSearch(imagette1, imagette2, imagette1mask=None, imagette2mask=None, returnError=False):
"""
This function performs a pixel-by-pixel search in 3D for a small reference imagette1 within a larger imagette2.
......@@ -789,6 +776,11 @@ def pixelSearch(imagette1, imagette2, imagette1mask=None, imagette2mask=None):
(True = Correlate these pixels, False = Skip),
Default = no mask
returnError : bool, optional
Return a normalised sum of squared differences error compatible with
register() family of functions? If yes, it's the last returned variable
Default = False
Returns
-------
p : 3-component vector
......@@ -820,7 +812,14 @@ def pixelSearch(imagette1, imagette2, imagette1mask=None, imagette2mask=None):
imagette2.astype("<f4"),
returns)
return numpy.array(returns[0:3]), returns[3]
if returnError:
error = _errorCalc(imagette1, imagette2[int(returns[0]):int(returns[0])+imagette1.shape[0],
int(returns[1]):int(returns[1])+imagette1.shape[1],
int(returns[2]):int(returns[2])+imagette1.shape[2]])
return numpy.array(returns[0:3]), returns[3], error
else:
return numpy.array(returns[0:3]), returns[3]
def globalCorrelation(im1, im2, mesh, convergenceCriterion=0.01, debugFiles=False, maxIterations=20, initialDisplacements=None, boundaryConditions=None, prefix="./"):
......
......@@ -180,7 +180,7 @@ def writeStrainTSV(fileName, points, decomposedFfield, firstColumn="StrainPointN
header=header)
def readCorrelationTSV(fileName, fieldBinRatio=1.0, readOnlyDisplacements=False, readConvergence=True, readPSCC=False):
def readCorrelationTSV(fileName, fieldBinRatio=1.0, readOnlyDisplacements=False, readConvergence=True, readPixelSearchCC=False):
"""
This function reads a TSV file containing a field of deformation functions **Phi** at one or a number of points.
This is typically the output of the spam-ldic and spam-ddic scripts,
......@@ -205,8 +205,8 @@ def readCorrelationTSV(fileName, fieldBinRatio=1.0, readOnlyDisplacements=False,
Read "returnStatus", "deltaPhiNorm", "iterations" from file
Default = True
readPSCC : bool, optional
Read "PSCC" from file
readPixelSearchCC : bool, optional
Read "pixelSearchCC" from file
Default = False
Returns
......@@ -226,7 +226,7 @@ def readCorrelationTSV(fileName, fieldBinRatio=1.0, readOnlyDisplacements=False,
iterations: nx1 array of n points iterations from the correlation
PSCC: nx1 array of n points PSCC from the correlation
pixelSearchCC: nx1 array of n points pixelSearchCC from the correlation
"""
if not os.path.isfile(fileName):
......@@ -236,7 +236,7 @@ def readCorrelationTSV(fileName, fieldBinRatio=1.0, readOnlyDisplacements=False,
f = numpy.genfromtxt(fileName, delimiter="\t", names=True)
# RS = []
# deltaPhiNorm = []
# PSCC = []
# pixelSearchCC = []
# If this is a one-line TSV file (an initial registration for example)
if f.size == 1:
......@@ -303,7 +303,7 @@ def readCorrelationTSV(fileName, fieldBinRatio=1.0, readOnlyDisplacements=False,
deltaPhiNorm = f['SubPixDeltaPhiNorm']
iterations = f['SubPixIterations']
PSCC = 0
pixelSearchCC = 0
# there is more than one line in the TSV file -- a field -- typical case
else:
......@@ -354,10 +354,10 @@ def readCorrelationTSV(fileName, fieldBinRatio=1.0, readOnlyDisplacements=False,
iterations = f[:]['SubPixIterations']
# Return pixelSearchCC
if readPSCC:
PSCC = numpy.zeros(nPoints)
if readPixelSearchCC:
pixelSearchCC = numpy.zeros(nPoints)
try:
PSCC = f[:]['pixelSearchCC']
pixelSearchCC = f[:]['pixelSearchCC']
except ValueError:
pass
......@@ -409,8 +409,8 @@ def readCorrelationTSV(fileName, fieldBinRatio=1.0, readOnlyDisplacements=False,
output.update({"returnStatus": RS,
"deltaPhiNorm": deltaPhiNorm,
"iterations": iterations})
if readPSCC:
output.update({"PSCC": PSCC})
if readPixelSearchCC:
output.update({"pixelSearchCC": pixelSearchCC})
if readOnlyDisplacements:
output.update({"displacements": PhiField[:, 0:3, -1]})
......@@ -432,13 +432,13 @@ def readStrainTSV(fileName):
Returns
-------
Dictionary containing:
fieldDims: 1x3 array of the field dimensions (ZYX)
fieldCoords : nx3 array of the field coordinates (ZYX)
numberOfLabels: number of labels (for a discrete strain result)
vol: nx1 array of n points with volumetric strain computed under the hypotesis of large strains (if computed)
dev: nx1 array of n points with deviatoric strain computed under the hypotesis of large strains (if computed)
......@@ -478,7 +478,7 @@ def readStrainTSV(fileName):
fieldCoords[:, 1] = f['Ypos']
fieldCoords[:, 2] = f['Xpos']
output['fieldCoords'] = fieldCoords
#Check if we are working with a regular grid or discrete
grid = False
discrete = False
......@@ -486,7 +486,7 @@ def readStrainTSV(fileName):
grid = True
else:
discrete = True
if grid:
fieldDims = numpy.array([len(numpy.unique(f['Zpos'])), len(numpy.unique(f['Ypos'])), len(numpy.unique(f['Xpos']))])
output['fieldDims'] = fieldDims
......@@ -494,7 +494,7 @@ def readStrainTSV(fileName):
else:
output['fieldDims'] = [0, 0, 0]
output['numberOfLabels'] = nPoints
#Check for all the possible keys
if 'vol' in keys:
volStrain = numpy.zeros((nPoints, 1))
......
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