Commit 22c35a60 authored by Olga Stamati's avatar Olga Stamati
Browse files

spam-pixelSearchPropagate fixed and tested for grid and label modes

parent d28ed01f
Pipeline #63086 failed with stages
in 13 minutes and 12 seconds
......@@ -102,7 +102,8 @@ elif args.LAB1 is not None:
centresOfMass = spam.label.centresOfMass(lab1, boundingBoxes=boundingBoxes)
im1mask = None
im2mask = None
gp = numpy.vstack([startPoint, centresOfMass])
gp = centresOfMass.copy()
gp[0] = startPoint
else:
gp, nodesDim = spam.DIC.makeGrid(im1.shape, args.NS)
......@@ -110,16 +111,15 @@ else:
print("\n\tRanking points")
guidingPoints = spam.mesh.rankPoints(gp, neighbourRadius=args.RADIUS)
guidingPoints, rowNumbers = spam.mesh.rankPoints(gp, neighbourRadius=args.RADIUS)
numberOfPoints = guidingPoints.shape[0]
print(guidingPoints[0:10])
# Initialise arrays
PhiField = numpy.zeros((numberOfPoints, 4, 4))
for point in range(numberOfPoints):
PhiField[point] = numpy.eye(4)
PhiField[0, 0:3, -1] += startPointDisplacement
pixelSearchCC = numpy.zeros((numberOfPoints), dtype=float)
# Returns compatible with register()
......@@ -133,16 +133,17 @@ widgets = [progressbar.FormatLabel(''), ' ', progressbar.Bar(), ' ', progressbar
pbar = progressbar.ProgressBar(widgets=widgets, maxval=numberOfPoints)
pbar.start()
# Step 1: simple PS for first point
if args.LAB1 is not None:
imagetteReturnsTop = spam.label.getImagettesLabelled(lab1, nodeNumber,
imagetteReturnsTop = spam.label.getImagettesLabelled(lab1, lab1[startPoint[0], startPoint[1], startPoint[2]],
PhiField[0].copy(),
im1, im2,
searchRange.copy(),
boundingBoxes, nodePositions,
boundingBoxes, centresOfMass,
margin=args.LABEL_DILATE,
labelDilate=args.LABEL_DILATE,
applyF=args.APPLY_F,
applyF="no",
volumeThreshold=3**3)
imagetteReturnsTop['imagette2mask'] = None
else:
......@@ -156,12 +157,13 @@ if imagetteReturnsTop['returnStatus'] == 1:
imagette1mask = imagetteReturnsTop['imagette1mask'],
imagette2mask = imagetteReturnsTop['imagette2mask'],
returnError = True)
PhiField[0, 0:3, -1] += imagetteReturnsTop['pixelSearchOffset'] + PSreturnsTop[0]
PhiField[0, 0:3, -1] = PSreturnsTop[0] + imagetteReturnsTop['pixelSearchOffset']
pixelSearchCC[0] = PSreturnsTop[1]
error[0] = PSreturnsTop[2]
# Failed to extract imagettes or something
else:
print("Failed to extact correlation window for starting point, exiting")
print("Failed to extract correlation window for starting point, exiting")
exit()
if pixelSearchCC[0] < args.CC_MIN:
print("CC obtained for starting point is less than threshold, not continuing")
......@@ -188,8 +190,28 @@ for point in range(1, numberOfPoints):
initialDisplacement = numpy.sum(PhiField[indices, 0:3, -1]*weights[:, numpy.newaxis], axis=0)/weights.sum()
# 2.4: Call PS around the estimated position
PhiField[point, 0:3, -1] += initialDisplacement
imagetteReturns = spam.DIC.getImagettes(im1, guidingPoints[point], args.HWS, PhiField[point].copy(), im2, searchRange.copy(), applyF='no')
PhiField[point, 0:3, -1] = initialDisplacement
if args.LAB1 is not None:
imagetteReturns = spam.label.getImagettesLabelled(lab1, rowNumbers[point],
PhiField[0].copy(),
im1, im2,
searchRange.copy(),
boundingBoxes, centresOfMass,
margin=args.LABEL_DILATE,
labelDilate=args.LABEL_DILATE,
applyF='no',
volumeThreshold=3**3)
imagetteReturns['imagette2mask'] = None
else:
imagetteReturns = spam.DIC.getImagettes(im1,
guidingPoints[point],
args.HWS,
PhiField[point].copy(),
im2, searchRange.copy(),
greyThreshold=[args.GREY_LOW_THRESH, args.GREY_HIGH_THRESH],
applyF='no')
if imagetteReturns['returnStatus'] == 1:
PSreturns = spam.DIC.pixelSearch(imagetteReturns['imagette1'],
......@@ -197,7 +219,8 @@ for point in range(1, numberOfPoints):
imagette1mask = imagetteReturns['imagette1mask'],
imagette2mask = imagetteReturns['imagette2mask'],
returnError = True)
PhiField[point, 0:3, -1] += imagetteReturns['pixelSearchOffset'] + PSreturns[0]
PhiField[point, 0:3, -1] = PSreturns[0] + imagetteReturns['pixelSearchOffset']
pixelSearchCC[point] = PSreturns[1]
error[point] = PSreturns[2]
returnStatus[point] = imagetteReturns['returnStatus']
......@@ -209,7 +232,26 @@ for point in range(1, numberOfPoints):
error[point] = numpy.inf
returnStatus[point] = imagetteReturns['returnStatus']
print("")
# Detect regular grid mode
if args.GUIDING_POINTS_FILE is None and args.LAB1 is None:
rowNumbers = rowNumbers[1:]
guidingPoints = guidingPoints[rowNumbers]
PhiField = PhiField[rowNumbers]
error = error[rowNumbers]
returnStatus = returnStatus[rowNumbers]
pixelSearchCC = pixelSearchCC[rowNumbers]
deltaPhiNorm = deltaPhiNorm[rowNumbers]
iterations = iterations[rowNumbers]
if args.GUIDING_POINTS_FILE is None and args.LAB1 is None:
if args.TIFF:
tifffile.imsave(args.OUT_DIR+"/"+args.PREFIX+"-Zdisp.tif", PhiField[:, 0, -1].astype('<f4').reshape(nodesDim))
tifffile.imsave(args.OUT_DIR+"/"+args.PREFIX+"-Ydisp.tif", PhiField[:, 1, -1].astype('<f4').reshape(nodesDim))
tifffile.imsave(args.OUT_DIR+"/"+args.PREFIX+"-Xdisp.tif", PhiField[:, 2, -1].astype('<f4').reshape(nodesDim))
tifffile.imsave(args.OUT_DIR+"/"+args.PREFIX+"-CC.tif", pixelSearchCC.astype('<f4').reshape(nodesDim))
tifffile.imsave(args.OUT_DIR+"/"+args.PREFIX+"-returnStatus.tif", returnStatus.astype('<f4').reshape(nodesDim))
if args.TSV:
# Make one big array for writing:
# First the node number,
......@@ -217,7 +259,7 @@ if args.TSV:
# F[0:3,0:3]
# Pixel-search CC
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(numberOfPoints)),
outMatrix = numpy.array([numpy.array(range(guidingPoints.shape[0])),
guidingPoints[:, 0], guidingPoints[:, 1], guidingPoints[:, 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],
......@@ -229,24 +271,39 @@ if args.TSV:
iterations]).T
numpy.savetxt(args.OUT_DIR+"/"+args.PREFIX+".tsv",
outMatrix,
fmt='%.7f',
delimiter='\t',
newline='\n',
comments='',
header=TSVheader)
outMatrix,
fmt='%.7f',
delimiter='\t',
newline='\n',
comments='',
header=TSVheader)
# Collect data into VTK output
if args.VTK:
# boring nans overwriting
disp = PhiField[:, 0:3, -1]
disp[numpy.logical_not(numpy.isfinite(disp))] = 0
if args.GUIDING_POINTS_FILE is None and args.LAB1 is None:
cellData = {}
cellData['displacements'] = PhiField[:, :-1, 3].reshape((nodesDim[0], nodesDim[1], nodesDim[2], 3))
cellData['pixelSearchCC'] = pixelSearchCC.reshape(nodesDim)
# Overwrite nans and infs with 0, rubbish I know
cellData['displacements'][numpy.logical_not(numpy.isfinite(cellData['displacements']))] = 0
# This is perfect in the case where NS = 2xHWS, these cells will all be in the right place
# In the case of overlapping of under use of data, it should be approximately correct
# If you insist on overlapping, then perhaps it's better to save each point as a cube glyph
# and actually *have* overlapping
spam.helpers.writeStructuredVTK(origin=nodePositions[0]-args.HWS, aspectRatio=args.NS, cellData=cellData, fileName=args.OUT_DIR+"/"+args.PREFIX+".vtk")
else:
# boring nans overwriting
disp = PhiField[:, 0:3, -1]
disp[numpy.logical_not(numpy.isfinite(disp))] = 0
magDisp = numpy.linalg.norm(disp, axis=1)
magDisp[numpy.logical_not(numpy.isfinite(magDisp))] = 0
magDisp = numpy.linalg.norm(disp, axis=1)
VTKglyphDict = {'displacements': disp,
'mag(displacements)': magDisp ,
'pixelSearchCC': pixelSearchCC}
VTKglyphDict = {'displacements': disp,
'mag(displacements)': magDisp ,
'pixelSearchCC': pixelSearchCC
}
spam.helpers.writeGlyphsVTK(guidingPoints, VTKglyphDict, fileName=args.OUT_DIR + "/" + args.PREFIX + ".vtk")
spam.helpers.writeGlyphsVTK(guidingPoints, VTKglyphDict, fileName=args.OUT_DIR + "/" + args.PREFIX + ".vtk")
\ No newline at end of file
......@@ -47,6 +47,9 @@ class testAll(unittest.TestCase):
rm(testFolder+"balls-5a-pixelSearch.tsv")
rm(testFolder+"balls-5b-ddic.tsv")
rm(testFolder+"balls-5b-ddic.vtk")
rm(testFolder+"balls-5c-pixelSearchPropagate.tsv")
rm(testFolder+"balls-5d-ddic.tsv")
rm(testFolder+"balls-5d-ddic.vtk")
rm(testFolder+"balls-6a-pixelSearch.tsv")
rm(testFolder+"extreme-im1-extreme-im2-ddic.tsv")
rm(testFolder+"extreme-im1-extreme-im2-ddic.vtk")
......@@ -79,6 +82,7 @@ class testAll(unittest.TestCase):
rm(testFolder+"snow-grid-3b-pixelSearch.tsv")
rm(testFolder+"snow-grid-3c-pixelSearch.tsv")
rm(testFolder+"snow-grid-3e-pixelSearch.tsv")
rm(testFolder+"snow-grid-4a-pixelSearchPropagate.tsv")
rm(testFolder+"snow-grid-5a-ldic-deltaPhiNorm.tif")
rm(testFolder+"snow-grid-5a-ldic-error.tif")
rm(testFolder+"snow-grid-5a-ldic-iterations.tif")
......@@ -88,6 +92,7 @@ class testAll(unittest.TestCase):
rm(testFolder+"snow-grid-5a-ldic-Ydisp.tif")
rm(testFolder+"snow-grid-5a-ldic-Zdisp.tif")
rm(testFolder+"snow-grid-5b-ldic.tsv")
rm(testFolder+"snow-grid-5d-ldic.tsv")
rm(testFolder+"snow-grid-6a-rx-Geers.tif")
rm(testFolder+"snow-grid-6a-ry-Geers.tif")
rm(testFolder+"snow-grid-6a-rz-Geers.tif")
......@@ -254,7 +259,6 @@ class testAll(unittest.TestCase):
# And the z-displacement is low
self.assertTrue(numpy.isclose(0, numpy.median(PSresult3a['PhiField'][PSresult3a['returnStatus']==1,0,-1]), atol=1.0))
### Step 3b load initial ONLY ROTATION guess and do a big search around the applied displacement
exitCode = subprocess.call(["spam-pixelSearch", "-pf", testFolder + "snow-grid-ereg-onlyRot.tsv",
"-sr",
......@@ -338,7 +342,26 @@ class testAll(unittest.TestCase):
#######################################################
### Step 4 check spam-pixelSearchPropagate
#######################################################
# ???
### check only with dip
exitCode = subprocess.call(["spam-pixelSearchPropagate",
"-sp", "{}".format(snowRef.shape[0]//2), "{}".format(snowRef.shape[1]//2), "{}".format(snowRef.shape[2]//2),
"{}".format(int(refTranslation[0])), "{}".format(int(refTranslation[1])), "{}".format(int(refTranslation[2])),
"-sr", "-2", "2", "-2", "2", "-2", "2",
"-glt", "5000", "-hws", "15",
testFolder + "snow-ref.tif", testFolder + "snow-def-onlyDisp.tif",
"-od", testFolder + "",
"-pre", "snow-grid-4a"
])
self.assertEqual(exitCode, 0)
PSresult4a = spam.helpers.readCorrelationTSV(testFolder + "snow-grid-4a-pixelSearchPropagate.tsv", readPixelSearchCC=True)
# Assert that the CC is nice and high
self.assertTrue(0.95 < numpy.median(PSresult4a['pixelSearchCC'][PSresult4a['returnStatus']==1]))
# Check displacement vector
self.assertTrue(numpy.isclose(refTranslation[0], numpy.median(PSresult4a['PhiField'][PSresult4a['returnStatus']==1,0,-1]), atol=0.1))
self.assertTrue(numpy.isclose(refTranslation[1], numpy.median(PSresult4a['PhiField'][PSresult4a['returnStatus']==1,1,-1]), atol=0.1))
self.assertTrue(numpy.isclose(refTranslation[2], numpy.median(PSresult4a['PhiField'][PSresult4a['returnStatus']==1,2,-1]), atol=0.1))
#######################################################
### Step 5 check spam-ldic loading previous results
......@@ -387,12 +410,11 @@ class testAll(unittest.TestCase):
self.assertTrue(0.8 < numpy.sum(LDICresult5b['returnStatus']==2)/len(LDICresult5b['returnStatus']))
# median Rotation within 0.5 deg?
self.assertTrue(numpy.isclose(refRotation[0], numpy.median(spam.deformation.decomposePhiField(LDICresult5a['PhiField'][LDICresult5a['returnStatus']>0], components='r')['r'][:,0]), atol=0.5))
self.assertTrue(numpy.isclose(refRotation[0], numpy.median(spam.deformation.decomposePhiField(LDICresult5b['PhiField'][LDICresult5b['returnStatus']>0], components='r')['r'][:,0]), atol=0.5))
# And the z-displacement is low
self.assertTrue(numpy.isclose(0, numpy.median(LDICresult5b['PhiField'][LDICresult5b['returnStatus']==2,0,-1]), atol=1.0))
### Step "5c" check 2D mode
exitCode = subprocess.call(["spam-ldic",
"-pf", testFolder + "snow-grid-2D-registration.tsv",
......@@ -414,8 +436,21 @@ class testAll(unittest.TestCase):
# And the z-displacement is low
self.assertTrue(numpy.isclose(0, numpy.median(LDICresult5c['PhiField'][LDICresult5c['returnStatus']==2,0,-1]), atol=1.0))
### Step 5d: load pixelSearchPropagate results with -pfd (same grid) + ug
exitCode = subprocess.call(["spam-ldic",
"-pf", testFolder + "snow-grid-4a-pixelSearchPropagate.tsv", "-pfd",
"-glt", "5000", "-hws", "15",
testFolder + "snow-ref.tif", testFolder + "snow-def-onlyDisp.tif",
"-od", testFolder + "",
"-pre", "snow-grid-5d",
])
self.assertEqual(exitCode, 0)
LDICresult5d = spam.helpers.readCorrelationTSV(testFolder + "snow-grid-5d-ldic.tsv")
# Check displacement vector
self.assertTrue(numpy.isclose(refTranslation[0], numpy.median(LDICresult5d['PhiField'][LDICresult5d['returnStatus']==2,0,-1]), atol=0.1))
self.assertTrue(numpy.isclose(refTranslation[1], numpy.median(LDICresult5d['PhiField'][LDICresult5d['returnStatus']==2,1,-1]), atol=0.1))
self.assertTrue(numpy.isclose(refTranslation[2], numpy.median(LDICresult5d['PhiField'][LDICresult5d['returnStatus']==2,2,-1]), atol=0.1))
#######################################################
## Step 6 onto the regularStrain
......@@ -807,6 +842,44 @@ class testAll(unittest.TestCase):
self.assertTrue(numpy.isclose(numpy.median(DDICresult5b['PhiField'][1:,2,-1]), translationStep2[2], atol=0.3))
self.assertTrue(numpy.median(DDICresult5b['returnStatus']) == 2)
######################################################
### 5c Run ddic -- Step0 -> Step2 with pixel search propagate
#######################################################
# Just run a simple DVC with no outputs
exitCode = subprocess.call(["spam-pixelSearchPropagate",
"-sr", "-2", "2", "-2", "2", "-2", "2",
"-sp", "{}".format(int(COMref[1, 0])), "{}".format(int(COMref[1, 1])), "{}".format(int(COMref[1, 2])),
"{}".format(int(translationStep2[0])), "{}".format(int(translationStep2[1])), "{}".format(int(translationStep2[2])),
testFolder + "Step0.tif",
testFolder + "Step2.tif",
"-lab1", testFolder + "Lab0.tif",
"-od", testFolder + "",
"-pre", "balls-5c"
])
self.assertEqual(exitCode, 0)
PSresult5c = spam.helpers.readCorrelationTSV(testFolder + "balls-5c-pixelSearchPropagate.tsv")
self.assertTrue(numpy.sum(PSresult5c['returnStatus']) == len(radii))
# Check displacement vector
self.assertTrue(numpy.isclose(numpy.median(PSresult5c['PhiField'][1:,0,-1]), translationStep2[0], atol=0.3))
self.assertTrue(numpy.isclose(numpy.median(PSresult5c['PhiField'][1:,1,-1]), translationStep2[1], atol=0.3))
self.assertTrue(numpy.isclose(numpy.median(PSresult5c['PhiField'][1:,2,-1]), translationStep2[2], atol=0.3))
### 5d ddic
# Just run a simple DVC with no outputs except TSV
exitCode = subprocess.call(["spam-ddic",
"-pf", testFolder + "balls-5c-pixelSearchPropagate.tsv", "-pfd",
testFolder + "Step0.tif",
testFolder + "Lab0.tif",
testFolder + "Step2.tif",
"-od", testFolder + "",
"-pre", "balls-5d"])
self.assertEqual(exitCode, 0)
DDICresult5d = spam.helpers.readCorrelationTSV(testFolder + "balls-5d-ddic.tsv")
self.assertTrue(numpy.isclose(numpy.median(DDICresult5d['PhiField'][1:,0,-1]), translationStep2[0], atol=0.3))
self.assertTrue(numpy.isclose(numpy.median(DDICresult5d['PhiField'][1:,1,-1]), translationStep2[1], atol=0.3))
self.assertTrue(numpy.isclose(numpy.median(DDICresult5d['PhiField'][1:,2,-1]), translationStep2[2], atol=0.3))
self.assertTrue(numpy.median(DDICresult5d['returnStatus']) == 2)
#######################################################
### 6. Rerun ddic -- Step0 -> Step2 with prev result
......@@ -860,7 +933,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,
......@@ -900,6 +973,9 @@ class testAll(unittest.TestCase):
exitCode = subprocess.call(["spam-pixelSearch", "--help"], stdout=FNULL, stderr=FNULL)
self.assertEqual(exitCode, 0)
exitCode = subprocess.call(["spam-pixelSearchPropagate", "--help"], stdout=FNULL, stderr=FNULL)
self.assertEqual(exitCode, 0)
exitCode = subprocess.call(["spam-ldic", "--help"], stdout=FNULL, stderr=FNULL)
self.assertEqual(exitCode, 0)
......
......@@ -2222,12 +2222,12 @@ def pixelSearchPropagate(parser):
type=argparse.FileType('r'),
help="Path to tiff file containing a labelled image 1 that defines zones to correlate. Disactivates -hws and -ns options")
#parser.add_argument('-ld',
#'--label-dilate',
#type=int,
#default=1,
#dest='LABEL_DILATE',
#help="Only if -lab1 is defined: Number of times to dilate labels. Default = 1")
parser.add_argument('-ld',
'--label-dilate',
type=int,
default=1,
dest='LABEL_DILATE',
help="Only if -lab1 is defined: Number of times to dilate labels. Default = 1")
parser.add_argument('-mf1',
'--maskFile1',
......@@ -2354,6 +2354,14 @@ def pixelSearchPropagate(parser):
dest='TSV',
help='Disactivate TSV output format. Default = False')
parser.add_argument('-tif',
'-tiff',
'--TIFFout',
'--TIFout',
action="store_true",
dest='TIFF',
help='Activate TIFFoutput format. Default = False')
args = parser.parse_args()
# 2019-04-05 EA: 2D image detection approved by Christophe Golke, update for shape 2019-08-29
......@@ -2364,10 +2372,11 @@ def pixelSearchPropagate(parser):
# twoD = False
#tiff.close()
# You really need a start point...
if args.START_POINT_DISP[0:3] == [-1, -1, -1]:
print("You need to input a starting point from which to propagate!\n(even if displacement is zero)")
exit()
if args.GUIDING_POINTS_FILE is None:
# You really need a start point...
if args.START_POINT_DISP[0:3] == [-1, -1, -1]:
print("You need to input a starting point from which to propagate!\n(even if displacement is zero)")
exit()
# 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:
......
......@@ -1089,6 +1089,10 @@ def rankPoints(points, neighbourRadius=20):
from scipy.spatial import KDTree
import progressbar
rowNumbers = numpy.zeros((points.shape[0]), dtype=int)
points = numpy.array([points[:,0], points[:,1], points[:,2], numpy.arange(points.shape[0])]).T
# counters
p = pR = 0
......@@ -1106,17 +1110,19 @@ def rankPoints(points, neighbourRadius=20):
while points.shape[0] >= 1:
# Try to see how many points can be found around the current point based on the given radius
treeCoord = KDTree(points)
indRadius = numpy.array(treeCoord.query_ball_point(pointsRanked[pR], neighbourRadius))
treeCoord = KDTree(points[:, 0:3])
indRadius = numpy.array(treeCoord.query_ball_point(pointsRanked[pR, 0:3], neighbourRadius))
indRadius = indRadius[numpy.argsort(indRadius)]
# if no points inside the given radius, just find the closest point
if len(indRadius)<1:
distance, ind = treeCoord.query(pointsRanked[pR], k=1)
distance, ind = treeCoord.query(pointsRanked[pR, 0:3], k=1)
indRadius = numpy.array([ind])
# fill in the ranked array with the point(s) found
pointsRanked[p+1:p+1+len(indRadius)] = points[indRadius]
for qn, q in enumerate (range(p+1, p+1+len(indRadius))):
rowNumbers[int(points[indRadius[qn], -1])] = q
# remove ranked point(s) from input array
points = numpy.delete(points, indRadius, 0)
......@@ -1124,8 +1130,8 @@ def rankPoints(points, neighbourRadius=20):
p += len(indRadius)
pR += 1
widgets[0] = progressbar.FormatLabel("{:.1f}%%".format((pR/numberOfPoints)*100))
widgets[0] = progressbar.FormatLabel("{:.1f}%%".format((p/numberOfPoints)*100))
pbar.update(pR)
return pointsRanked
return pointsRanked[:,0:3], rowNumbers
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment