Skip to content
GitLab
Commit f965e6cc authored by Edward Andò's avatar Edward Andò
Browse files

[skip-ci] 2D pixel search tests and big improvement of rounding errors

parent a731c949
Pipeline #56332 skipped
Hide whitespace changes
Inline Side-by-side
scripts/spam-pixelSearch
View file @ f965e6cc
......@@ -237,14 +237,14 @@ def pixelSearchOneNode(nodeNumber):
imagetteReturns['returnStatus'] = 0
else:
imagetteReturns = spam.DIC.getImagettes(nodePositions[nodeNumber], PhiField[nodeNumber].copy(), searchRange.copy(), boundingBoxes[nodeNumber], im1, im2, im1mask, args.MASK_COVERAGE, greyThreshold)
imagetteReturns = spam.DIC.getImagettes(nodePositions[nodeNumber], PhiField[nodeNumber].copy(), searchRange.copy(), boundingBoxes[nodeNumber], im1, im2, im1mask, args.MASK_COVERAGE, greyThreshold, twoD)
# If getImagettes was successful (size check and mask coverage check)
if imagetteReturns['returnStatus'] == 1:
PSreturns = spam.DIC.correlate.pixelSearch(imagetteReturns['imagette1'],
imagetteReturns['imagette2'],
imagette1mask = imagetteReturns['imagette1mask'],
returnError = True)
imagetteReturns['imagette2'],
imagette1mask = imagetteReturns['imagette1mask'],
returnError = True)
pixelSearchOffset = imagetteReturns['pixelSearchOffset']
writeReturns = True
......
tests/test_scripts.py
View file @ f965e6cc
......@@ -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 = numpy.array([0.0, 15.0, 0.0])
refRotation = numpy.array([5.0, 0.0, 0.0])
refTranslation = numpy.array([0.0, 15, 0.0])
refRotation = numpy.array([3.0, 0.0, 0.0])
testFolder = "./"
try:
......@@ -40,18 +40,21 @@ class testAll(unittest.TestCase):
def rm(fn):
if os.path.isfile(fn): os.remove(fn)
try:
rm(testFolder+"snow-ref.tif")
rm(testFolder+"snow-def.tif")
#rm(testFolder+"snow-ref.tif")
#rm(testFolder+"snow-ref-2D.tif")
#rm(testFolder+"snow-def.tif")
#rm(testFolder+"snow-def-2D.tif")
rm(testFolder+"snow-def-onlyDisp.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")
rm(testFolder+"snow-ref-snow-def-ereg-onlyRot.tsv")
rm(testFolder + "snow-ref-snow-def-ereg-onlyRot.tsv")
rm(testFolder+"snow-ref-snow-def-ereg.tsv")
rm(testFolder+"snow-ref-snow-def-pixelSearch-CC.tif")
rm(testFolder+"snow-ref-snow-def-pixelSearch-returnStatus.tif")
rm(testFolder+"snow-ref-snow-def-pixelSearch-Xdisp.tif")
#rm(testFolder+"snow-ref-snow-def-pixelSearch-CC.tif")
#rm(testFolder+"snow-ref-snow-def-pixelSearch-returnStatus.tif")
#rm(testFolder+"snow-ref-snow-def-pixelSearch-Xdisp.tif")
rm(testFolder+"snow-ref-snow-def-pixelSearch-Ydisp.tif")
rm(testFolder+"snow-ref-snow-def-pixelSearch-Zdisp.tif")
rm(testFolder+"snow-ref-snow-def-PSCC.tif")
......@@ -112,9 +115,13 @@ class testAll(unittest.TestCase):
# Run on snow data, with bin2 registration, without output
# load 3D image
snowRef = spam.datasets.loadSnow()
snowRef2D = snowRef[50]
# save it locally
tifffile.imsave(testFolder + "snow-ref.tif", snowRef)
os.chmod( testFolder + "snow-ref.tif", 0o666)
tifffile.imsave(testFolder + "snow-ref-2D.tif", snowRef2D)
os.chmod(testFolder + "snow-ref-2D.tif", 0o666)
# Generate deformed image
Phi = spam.deformation.computePhi({'t': refTranslation,
'r': refRotation})
......@@ -122,6 +129,9 @@ class testAll(unittest.TestCase):
# save it locally
tifffile.imsave(testFolder + "snow-def.tif", snowDef)
os.chmod(testFolder + "snow-def.tif", 0o666)
# Plus 2D version
tifffile.imsave(testFolder + "snow-def-2D.tif", spam.DIC.applyPhiPython(snowRef2D[numpy.newaxis, ...], Phi=Phi))
os.chmod(testFolder + "snow-def-2D.tif", 0o666)
snowDefOnlyDisp = spam.DIC.applyPhi(snowRef, Phi=spam.deformation.computePhi({'t': refTranslation}))
# save it locally
......@@ -132,8 +142,10 @@ class testAll(unittest.TestCase):
snowMask = numpy.ones_like(snowDef, dtype='<u1')
snowMask[0] = 0
snowMask[-1] = 0
tifffile.imsave(testFolder + "snow-mask.tif", snowDef)
tifffile.imsave(testFolder + "snow-mask.tif", snowMask)
os.chmod(testFolder + "snow-mask.tif", 0o666)
#tifffile.imsave(testFolder + "snow-mask-2D.tif", snowMask[50])
#os.chmod(testFolder + "snow-mask-2D.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,
......@@ -151,6 +163,14 @@ class testAll(unittest.TestCase):
'returnStatus': 2,
'deltaPhiNorm': 1})
# generate a fake "eye reg" initial guess which is close in 2D
spam.helpers.writeRegistrationTSV(testFolder + "snow-ref-2D-snow-def-2D-ereg-onlyRot.tsv", [0, snowRef.shape[1]-1/2, snowRef.shape[2]-1/2],
{'Phi': spam.deformation.computePhi({'r': refRotation}),
'error': 0,
'iterations': 0,
'returnStatus': 2,
'deltaPhiNorm': 1})
#######################################################
### Step 2 check spam-reg functionality
#######################################################
......@@ -174,15 +194,16 @@ class testAll(unittest.TestCase):
"-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"])
"-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())
self.assertTrue(0.95 < numpy.median(PSresult['pixelSearchCC'][PSresult['returnStatus']==1]))
# And the z-displacement is low
self.assertTrue(numpy.isclose(0, PSresult['PhiField'][PSresult['returnStatus']==1,0,-1].mean(), atol=1.0))
self.assertTrue(numpy.isclose(0, numpy.median(PSresult['PhiField'][PSresult['returnStatus']==1,0,-1]), atol=1.0))
### Step 3.2 load initial ONLY ROTATION guess and do a big search around the applied displacement
exitCode = subprocess.call(["spam-pixelSearch", "-pf", testFolder + "snow-ref-snow-def-ereg-onlyRot.tsv",
......@@ -192,17 +213,37 @@ class testAll(unittest.TestCase):
"{}".format(int(refTranslation[2]-2)), "{}".format(int(refTranslation[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"])
"-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())
self.assertTrue(0.95 < numpy.median(PSresult['pixelSearchCC'][PSresult['returnStatus']==1]))
# And the z-displacement is low
self.assertTrue(numpy.isclose(0, PSresult['PhiField'][PSresult['returnStatus']==1,0,-1].mean(), atol=1.0))
self.assertTrue(numpy.isclose(0, numpy.median(PSresult['PhiField'][PSresult['returnStatus']==1,0,-1]), atol=1.0))
### Step 3.3 NO ROTATION IMAGE big search around the applied displacement
### Step 3.3 check 2D mode
exitCode = subprocess.call(["spam-pixelSearch", "-pf", testFolder + "snow-ref-2D-snow-def-2D-ereg-onlyRot.tsv",
"-sr",
"{}".format(int(refTranslation[0]-2)), "{}".format(int(refTranslation[0]+2)),
"{}".format(int(refTranslation[1]-2)), "{}".format(int(refTranslation[1]+2)),
"{}".format(int(refTranslation[2]-2)), "{}".format(int(refTranslation[2]+2)),
"-glt", "5000", "-hws", "10", "-ns", "5",
testFolder + "snow-ref-2D.tif", testFolder + "snow-def-2D.tif",
"-od", testFolder + "", "-tif",
#"-mf1", testFolder + "snow-mask.tif", "-mc", "1.0"
])
self.assertEqual(exitCode, 0)
PSresult = spam.helpers.readCorrelationTSV(testFolder + "snow-ref-2D-snow-def-2D-pixelSearch.tsv", readPixelSearchCC=True)
# Assert that the CC is nice and high
self.assertTrue(0.90 < numpy.median(PSresult['pixelSearchCC'][PSresult['returnStatus']==1]))
# And the z-displacement is low
self.assertTrue(numpy.isclose(0, numpy.median(PSresult['PhiField'][PSresult['returnStatus']==1,0,-1]), atol=1.0))
### Step 3.4 NO ROTATION IMAGE big search around the applied displacement
exitCode = subprocess.call(["spam-pixelSearch",
"-sr",
"{}".format(int(refTranslation[0]-2)), "{}".format(int(refTranslation[0]+2)),
......@@ -214,13 +255,14 @@ class testAll(unittest.TestCase):
self.assertEqual(exitCode, 0)
PSresult = spam.helpers.readCorrelationTSV(testFolder + "snow-ref-snow-def-onlyDisp-pixelSearch.tsv", readPixelSearchCC=True)
# 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(refTranslation[0], PSresult['PhiField'][PSresult['returnStatus']==1,0,-1].mean(), atol=1.0))
self.assertTrue(numpy.isclose(refTranslation[1], PSresult['PhiField'][PSresult['returnStatus']==1,1,-1].mean(), atol=1.0))
self.assertTrue(numpy.isclose(refTranslation[2], PSresult['PhiField'][PSresult['returnStatus']==1,2,-1].mean(), atol=1.0))
self.assertTrue(0.95 < numpy.median(PSresult['pixelSearchCC'][PSresult['returnStatus']==1]))
# Check correct values measured
self.assertTrue(numpy.isclose(refTranslation[0], numpy.median(PSresult['PhiField'][PSresult['returnStatus']==1,0,-1]), atol=0.1))
self.assertTrue(numpy.isclose(refTranslation[1], numpy.median(PSresult['PhiField'][PSresult['returnStatus']==1,1,-1]), atol=0.1))
self.assertTrue(numpy.isclose(refTranslation[2], numpy.median(PSresult['PhiField'][PSresult['returnStatus']==1,2,-1]), atol=0.1))
exit()
# Check output results with TSV output, bin1 registration + TIFFs
......
tools/DIC/grid.py
View file @ f965e6cc
......@@ -128,6 +128,9 @@ Parameters
If the mean is below the first value or above the second value, the grid point is not correlated.
Default = [ -inf, inf ]
twoD : bool, optional
Are passed imagette1 and imagette2 3D images?
Returns
-------
Dictionary containing:
......@@ -141,28 +144,33 @@ Returns
"""
# WARNING ------------------------- VVVVVVVVVVV--- gets easily overwritten, pass a .copy()!
def getImagettes(nodePosition, Phi, searchRange, boundingBox, im1, im2, im1mask, minMaskCoverage, greyThreshold):
def getImagettes(nodePosition, Phi, searchRange, boundingBox, im1, im2, im1mask, minMaskCoverage, greyThreshold, twoD=False):
returnStatus = 1
imagette1mask = None
initialDisplacement = Phi[0:3, 3].astype(int)
initialDisplacement = numpy.round(Phi[0:3, 3]).astype(int)
if twoD: m = 0
else: m = 0
# Catch bad bounding boxes:
if numpy.all((boundingBox[1::2] - boundingBox[0::2]) > [0, 0, 0]) and numpy.all(boundingBox[1::2] - boundingBox[0::2] != [0, 0, 0]):
if numpy.all((boundingBox[1::2] - boundingBox[0::2]) > 0) or ( numpy.all(boundingBox[3::2] - boundingBox[2::2] > 0) and twoD):
# 2020-09-25 OS and EA: Prepare startStop array for imagette 1 to be extracted with new slicePadded
PhiNoDisp = Phi.copy()
PhiNoDisp[0:3,-1] = 0.0
PhiNoDisp[0:3,-1] -= initialDisplacement
#PhiNoDisp[0:3,-1] = 0.0
# If F is not the identity, create a pad to be able to apply F to imagette 1
if numpy.allclose(PhiNoDisp, numpy.eye(4)):
applyPhiPad = 0
applyPhiPad = (0, 0, 0)
else:
# 2020-10-06 OS and EA: Add a pad to each dimension of 25% of max(halfWindowSize) to allow space to apply F (no displacement) to imagette1
applyPhiPad = int(0.5*numpy.ceil(max(boundingBox[1::2]-boundingBox[0::2])))
startStopIm1 = [int(boundingBox[0] - applyPhiPad), int(boundingBox[1] + applyPhiPad + 1),
int(boundingBox[2] - applyPhiPad), int(boundingBox[3] + applyPhiPad + 1),
int(boundingBox[4] - applyPhiPad), int(boundingBox[5] + applyPhiPad + 1)]
if twoD: applyPhiPad = ( 0, applyPhiPad, applyPhiPad)
else: applyPhiPad = (applyPhiPad, applyPhiPad, applyPhiPad)
startStopIm1 = [int(boundingBox[0] - applyPhiPad[0]), int(boundingBox[1] + applyPhiPad[0] + 1),
int(boundingBox[2] - applyPhiPad[1]), int(boundingBox[3] + applyPhiPad[1] + 1),
int(boundingBox[4] - applyPhiPad[2]), int(boundingBox[5] + applyPhiPad[2] + 1)]
# In either case, extract imagette1, now guaranteed to be the right size
imagette1padded = spam.helpers.slicePadded(im1, startStopIm1)
......@@ -173,11 +181,18 @@ def getImagettes(nodePosition, Phi, searchRange, boundingBox, im1, im2, im1mask,
imagette1def = imagette1padded
else:
# apply F to imagette 1 padded
imagette1paddedDef = spam.DIC.applyPhi(imagette1padded, PhiNoDisp)
if twoD: imagette1paddedDef = spam.DIC.applyPhiPython(imagette1padded, PhiNoDisp)
else: imagette1paddedDef = spam.DIC.applyPhi(imagette1padded, PhiNoDisp)
# undo padding
imagette1def = imagette1paddedDef[applyPhiPad:-applyPhiPad,
applyPhiPad:-applyPhiPad,
applyPhiPad:-applyPhiPad]
if twoD:
imagette1def = imagette1paddedDef[ : ,
applyPhiPad[1]:-applyPhiPad[1],
applyPhiPad[2]:-applyPhiPad[2]]
else:
imagette1def = imagette1paddedDef[applyPhiPad[0]:-applyPhiPad[0],
applyPhiPad[1]:-applyPhiPad[1],
applyPhiPad[2]:-applyPhiPad[2]]
### Check mask
if im1mask is None:
......@@ -191,7 +206,7 @@ def getImagettes(nodePosition, Phi, searchRange, boundingBox, im1, im2, im1mask,
# Make sure imagette is not 0-dimensional in any dimension
# Check minMaskVolume
if numpy.all(numpy.array(imagette1def.shape) > 0):
if numpy.all(numpy.array(imagette1def.shape) > 0) or (twoD and numpy.all(numpy.array(imagette1def.shape[1:3]) > 0)):
# ------------ Grey threshold low --------------- and -------------- Grey threshold high -----------
if numpy.nanmean(imagette1def) > greyThreshold[0] and numpy.nanmean(imagette1def) < greyThreshold[1]:
if maskVolumeCondition:
......
tools/helpers/optionsParser.py
View file @ f965e6cc
......@@ -18,6 +18,7 @@ this program. If not, see <http://www.gnu.org/licenses/>.
import argparse
import numpy
import tifffile
import os
# Nice str2bool suggestion from Maxim (https://stackoverflow.com/questions/15008758/parsing-boolean-values-with-argparse)
......@@ -229,9 +230,8 @@ def ldicParser(parser):
# print("DICregularGrid: Input TIFF files need to be writeable in order to guess their dimensionality")
# exit()
# 2019-03-21 EA: better check for dimensions, doesn't depend on writability of files
import tifffile
tiff = tifffile.TiffFile(args.inFiles[0].name)
imagejSingleSlice = True
#imagejSingleSlice = True
# if tiff.imagej_metadata is not None:
# if 'slices' in tiff.imagej_metadata:
# if tiff.imagej_metadata['slices'] > 1:
......@@ -1974,13 +1974,6 @@ def pixelSearch(parser):
dest='SEARCH_RANGE',
help='Z- Z+ Y- Y+ X- X+ ranges (in pixels) for the pxiel search. Requires pixel search to be activated. Default = +-3px')
#parser.add_argument('-psf',
#'--pixel-search-filter',
#type=int,
#default=0,
#dest='PS_FILTER',
#help='Median filter pixel search results. Default = 0')
# Default: node spacing equal in all three directions
parser.add_argument('-ns',
'--node-spacing',
......@@ -2017,20 +2010,6 @@ def pixelSearch(parser):
dest='HWS',
help="Half correlation window size, measured each side of the node pixel (different in 3 directions). Default = 10, 10, 10px")
#parser.add_argument('-bb',
#'--binning-begin',
#type=int,
#default=4,
#dest='BIN_BEGIN',
#help='Initial binning to apply to input images for initial registration. Default = 4')
#parser.add_argument('-be',
#'--binning-end',
#type=int,
#default=1,
#dest='BIN_END',
#help='Binning level to stop at for initial registration. Default = 1')
parser.add_argument('-glt',
'--grey-low-threshold',
type=float,
......@@ -2059,12 +2038,12 @@ def pixelSearch(parser):
dest='PREFIX',
help='Prefix for output files (without extension). Default is basename of mesh file')
parser.add_argument('-def',
'--save-deformed-image1',
action="store_true",
default=False,
dest='DEF',
help="Activate the saving of a deformed image 1 (as <im1>-reg-def.tif)")
#parser.add_argument('-def',
#'--save-deformed-image1',
#action="store_true",
#default=False,
#dest='DEF',
#help="Activate the saving of a deformed image 1 (as <im1>-reg-def.tif)")
parser.add_argument('-vtk',
'--VTKout',
......@@ -2088,6 +2067,14 @@ def pixelSearch(parser):
args = parser.parse_args()
# 2019-04-05 EA: 2D image detection approved by Christophe Golke, update for shape 2019-08-29
tiff = tifffile.TiffFile(args.im1.name)
if len(tiff.pages) == 1 and len(tiff.series[0].shape) == 2:
twoD = True
else:
twoD = False
tiff.close()
# If we have no out dir specified, deliver on our default promise -- this can't be done inline before since parser.parse_args() has not been run at that stage.
if args.OUT_DIR is None:
args.OUT_DIR = os.path.dirname(args.im1.name)
......@@ -2136,6 +2123,16 @@ def pixelSearch(parser):
args.NS = None
args.MASK1 = None
args.MASK_COVERAGE = 0
# Catch and overwrite 2D options
if twoD:
args.SEARCH_RANGE[0] = 0
args.SEARCH_RANGE[1] = 0
# Catch and overwrite 2D options
elif twoD:
args.NS[0] = 1
args.HWS[0] = 0
args.SEARCH_RANGE[0] = 0
args.SEARCH_RANGE[1] = 0
# Output file name prefix
if args.PREFIX is None:
......
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