Vous avez reçu un message "Your GitLab account has been locked ..." ? Pas d'inquiétude : lisez cet article https://docs.gricad-pages.univ-grenoble-alpes.fr/help/unlock/

Commit a8f2db17 authored by Edward Andò's avatar Edward Andò
Browse files

Merge branch 'noBorders'

parents d361cba8 96ef61a9
......@@ -120,7 +120,7 @@ if mpiRank == boss or not mpi:
###############################################################
numberOfLabels = (lab1.max() + 1).astype('u4')
print("spam-ddic: Number of labels = {}\n".format(numberOfLabels))
print("spam-ddic: Number of labels = {}\n".format(numberOfLabels-1))
print("spam-ddic: Calculating Bounding Boxes and Centres of Mass of all labels.")
boundingBoxes = ltk.boundingBoxes(lab1)
......@@ -364,8 +364,11 @@ if mpiRank == boss or not mpi:
else:
maskette1 = None
maskette1vol = numpy.inf # must pass "if" below
slicette1 = getLabel['slice']
imagette1 = im1[slicette1].copy()
# Use new padded slicer, to remain aligned with getLabel['subvol']
# + add 1 on the "max" side for bounding box -> slice
imagette1 = spam.helpers.slicePadded(im1, getLabel['boundingBox'] + numpy.array([0,1,0,1,0,1]))
#imagette1 = im1[getLabel['slice']].copy()
if maskette1vol > args.VOLUME_THRESHOLD:
if args.PS == 'on' or (not registrationSuccessful and args.PS == 'auto'):
......@@ -381,16 +384,17 @@ if mpiRank == boss or not mpi:
(searchRangeForThisLabel['yRange'][0] + searchRangeForThisLabel['yRange'][1])/2,
(searchRangeForThisLabel['xRange'][0] + searchRangeForThisLabel['xRange'][1])/2] ).astype(int)
# Slice for image 2
subVolSlice2 = (slice(max(int(boundingBoxes[label][0] - max(labelDilateCurrent, 0) - args.LABEL_CORRELATE_MARGIN + searchRangeForThisLabel['zRange'][0]), 0),
min(int(boundingBoxes[label][1] + max(labelDilateCurrent, 0) + args.LABEL_CORRELATE_MARGIN + searchRangeForThisLabel['zRange'][1] + 1), im1.shape[0])),
slice(max(int(boundingBoxes[label][2] - max(labelDilateCurrent, 0) - args.LABEL_CORRELATE_MARGIN + searchRangeForThisLabel['yRange'][0]), 0),
min(int(boundingBoxes[label][3] + max(labelDilateCurrent, 0) + args.LABEL_CORRELATE_MARGIN + searchRangeForThisLabel['yRange'][1] + 1), im1.shape[1])),
slice(max(int(boundingBoxes[label][4] - max(labelDilateCurrent, 0) - args.LABEL_CORRELATE_MARGIN + searchRangeForThisLabel['xRange'][0]), 0),
min(int(boundingBoxes[label][5] + max(labelDilateCurrent, 0) + args.LABEL_CORRELATE_MARGIN + searchRangeForThisLabel['xRange'][1] + 1), im1.shape[2])))
# 2020-09-25 OS and EA: Prepare startStop array for imagette 2 to be extracted with new slicePadded
startStopIm2 = [int(int(boundingBoxes[label][0]) - max(labelDilateCurrent, 0) - args.LABEL_CORRELATE_MARGIN + searchRangeForThisLabel['zRange'][0] ),
int(int(boundingBoxes[label][1]) + max(labelDilateCurrent, 0) + args.LABEL_CORRELATE_MARGIN + searchRangeForThisLabel['zRange'][1] + 1),
int(int(boundingBoxes[label][2]) - max(labelDilateCurrent, 0) - args.LABEL_CORRELATE_MARGIN + searchRangeForThisLabel['yRange'][0] ),
int(int(boundingBoxes[label][3]) + max(labelDilateCurrent, 0) + args.LABEL_CORRELATE_MARGIN + searchRangeForThisLabel['yRange'][1] + 1),
int(int(boundingBoxes[label][4]) - max(labelDilateCurrent, 0) - args.LABEL_CORRELATE_MARGIN + searchRangeForThisLabel['xRange'][0] ),
int(int(boundingBoxes[label][5]) + max(labelDilateCurrent, 0) + args.LABEL_CORRELATE_MARGIN + searchRangeForThisLabel['xRange'][1] + 1)]
imagette2 = spam.helpers.slicePadded(im2, startStopIm2)
# Extract it...
imagette2 = im2[subVolSlice2]
# Catch
# point in im2 that we are searching around
searchCentre = (numpy.array(imagette2.shape, dtype='<f4') - 1) / 2.0 + labelDisplacementInt - middleOfSearchRange
......@@ -414,6 +418,7 @@ if mpiRank == boss or not mpi:
imagette1toCorrelate = imagette1def.copy()
maskette1def = spam.DIC.applyPhi(maskette1, PhiNoDisp, PhiPoint=getLabel['centreOfMassREL'], interpolationOrder=0)
imagette1toCorrelate[maskette1def == 0] = numpy.nan
else:
......@@ -460,18 +465,21 @@ if mpiRank == boss or not mpi:
if args.LABEL_CORRELATE:
labelDisplacementInt = numpy.round(PhiField[label][0:3, -1])
slicette2 = (slice(max(int(boundingBoxes[label][0] - args.LABEL_CORRELATE_MARGIN - max(labelDilateCurrent, 0) + labelDisplacementInt[0] ), 0 ),
min(int(boundingBoxes[label][1] + args.LABEL_CORRELATE_MARGIN + max(labelDilateCurrent, 0) + labelDisplacementInt[0] + 1), im1.shape[0])),
slice(max(int(boundingBoxes[label][2] - args.LABEL_CORRELATE_MARGIN - max(labelDilateCurrent, 0) + labelDisplacementInt[1] ), 0 ),
min(int(boundingBoxes[label][3] + args.LABEL_CORRELATE_MARGIN + max(labelDilateCurrent, 0) + labelDisplacementInt[1] + 1), im1.shape[1])),
slice(max(int(boundingBoxes[label][4] - args.LABEL_CORRELATE_MARGIN - max(labelDilateCurrent, 0) + labelDisplacementInt[2] ), 0 ),
min(int(boundingBoxes[label][5] + args.LABEL_CORRELATE_MARGIN + max(labelDilateCurrent, 0) + labelDisplacementInt[2] + 1), im1.shape[2])))
# 2020-09-25 OS and EA: Prepare startStop array for imagette 2 to be extracted with new
# slicePadded, this should solved "Boss: failed imDiff" and RS=-5 forever
startStopIm2 = [int(int(boundingBoxes[label][0]) - args.LABEL_CORRELATE_MARGIN - max(labelDilateCurrent, 0) + labelDisplacementInt[0] ),
int(int(boundingBoxes[label][1]) + args.LABEL_CORRELATE_MARGIN + max(labelDilateCurrent, 0) + labelDisplacementInt[0] + 1),
int(int(boundingBoxes[label][2]) - args.LABEL_CORRELATE_MARGIN - max(labelDilateCurrent, 0) + labelDisplacementInt[1] ),
int(int(boundingBoxes[label][3]) + args.LABEL_CORRELATE_MARGIN + max(labelDilateCurrent, 0) + labelDisplacementInt[1] + 1),
int(int(boundingBoxes[label][4]) - args.LABEL_CORRELATE_MARGIN - max(labelDilateCurrent, 0) + labelDisplacementInt[2] ),
int(int(boundingBoxes[label][5]) + args.LABEL_CORRELATE_MARGIN + max(labelDilateCurrent, 0) + labelDisplacementInt[2] + 1)]
imagette2imagette1sizeDiff = numpy.array(im2[slicette2].shape) - numpy.array(imagette1.shape)
imagette2 = spam.helpers.slicePadded(im2, startStopIm2)
# Catch register DimProblem
if not (imagette2imagette1sizeDiff < 0).any():
imagette2imagette1sizeDiff = numpy.array(imagette2.shape) - numpy.array(imagette1.shape)
# If all of imagette2 is nans it fell outside im2 (or in any case it's going to be difficult to correlate)
if not numpy.all(numpy.isnan(imagette2)):
# Remove int() part of displacement since it's already used to extract imagette2
PhiTemp = PhiField[label].copy()
PhiTemp[0:3, -1] -= labelDisplacementInt
......@@ -481,7 +489,7 @@ if mpiRank == boss or not mpi:
if not mpi:
returns = spam.DIC.correlate.registerMultiscale(imagette1,
im2[slicette2],
imagette2,
args.LABEL_CORRELATE_MULTISCALE_BINNING,
im1mask=maskette1,
margin=1,
......@@ -499,7 +507,7 @@ if mpiRank == boss or not mpi:
# build message for lukasKanade worker
m = {'label': label,
'im1': imagette1,
'im2': im2[slicette2],
'im2': imagette2,
'binStart': args.LABEL_CORRELATE_MULTISCALE_BINNING,
'im1mask': maskette1,
'PhiInit': PhiTemp,
......@@ -521,8 +529,8 @@ if mpiRank == boss or not mpi:
workersActive[worker] = True
else: # Regardless of MPI or single proc failed imDiff condition
finishedLabels += 1
returnStatus[label] = -5
print("\t\tBoss: Fail imDiff", imagette2imagette1sizeDiff)
returnStatus[label] = -4
print("\t\tBoss: Empty imagette 2 with initial displacement", labelDisplacementInt)
grainOK = True
else: # No args.LABEL_CORRELATE
......
......@@ -160,22 +160,17 @@ if mpiRank == boss or not mpi:
registrationSuccessful = False
if args.REG:
print("\tldic: Starting registration")
if twoD:
regMargin = int(args.REG_MARGIN*min(im1.shape[1:3]))
else:
regMargin = int(args.REG_MARGIN*min(im1.shape))
# Run registration
regReturns = spam.DIC.correlate.registerMultiscale(im1, im2,
args.REG_BIN_BEGIN, binStop=args.REG_BIN_END,
margin=int(args.REG_MARGIN * min(im1.shape)),
im1mask=im1mask,
interpolationOrder=1,
maxIterations=500,
deltaPhiMin=0.0001,
updateGradient=args.REG_UPDATE,
interpolator='C',
verbose=True)
args.REG_BIN_BEGIN, binStop=args.REG_BIN_END,
margin=int(args.REG_MARGIN * min(im1.shape)),
im1mask=im1mask,
interpolationOrder=1,
maxIterations=500,
deltaPhiMin=0.0001,
updateGradient=args.REG_UPDATE,
interpolator='C',
verbose=True)
if regReturns['returnStatus'] != 2:
print("spam-ddic: Registration did not converge, try increasing the registration margin?")
......
......@@ -145,7 +145,14 @@ class testAll(unittest.TestCase):
self.assertEqual(mask[10:110, 10:110, 10:110].sum(), 100**3)
## case 6: not touching im
startStop = numpy.array([100, 120, 100, 120, 100, 120])
startStop = numpy.array([100, 120, 40, 60, 40, 60])
imSliced, mask = spam.helpers.slicePadded(im, startStop, createMask = True)
self.assertEqual(list(imSliced.shape), [20, 20, 20])
self.assertEqual(numpy.isfinite(imSliced).sum(), 0)
self.assertEqual(mask.sum(), 0)
## case 7: not touching im
startStop = numpy.array([-120, -100, 40, 60, 40, 60])
imSliced, mask = spam.helpers.slicePadded(im, startStop, createMask = True)
self.assertEqual(list(imSliced.shape), [20, 20, 20])
self.assertEqual(numpy.isfinite(imSliced).sum(), 0)
......
......@@ -178,51 +178,66 @@ class TestFunctionLabel(unittest.TestCase):
def test_getLabel(self):
# No-option get label
gl = spam.label.getLabel(threeCubedLabelVol, 1)
self.assertEqual(numpy.array([[[1]]], dtype=bool), gl['subvol'].tolist())
self.assertEqual(gl['slice'][0].start, 1)
self.assertEqual(gl['slice'][0].stop, 2)
self.assertEqual(gl['slice'][1].start, 1)
self.assertEqual(gl['slice'][1].stop, 2)
self.assertEqual(gl['slice'][2].start, 1)
self.assertEqual(gl['slice'][2].stop, 2)
gottenLabel = spam.label.getLabel(threeCubedLabelVol, 1)
self.assertEqual(numpy.array([[[1]]], dtype=bool), gottenLabel['subvol'].tolist())
self.assertEqual(gottenLabel['slice'][0].start, 1)
self.assertEqual(gottenLabel['slice'][0].stop, 2)
self.assertEqual(gottenLabel['slice'][1].start, 1)
self.assertEqual(gottenLabel['slice'][1].stop, 2)
self.assertEqual(gottenLabel['slice'][2].start, 1)
self.assertEqual(gottenLabel['slice'][2].stop, 2)
self.assertEqual(gottenLabel['boundingBox'][0], 1)
self.assertEqual(gottenLabel['boundingBox'][1], 1)
self.assertEqual(gottenLabel['boundingBox'][2], 1)
self.assertEqual(gottenLabel['boundingBox'][3], 1)
self.assertEqual(gottenLabel['boundingBox'][4], 1)
self.assertEqual(gottenLabel['boundingBox'][5], 1)
COM = spam.label.centresOfMass(threeCubedLabelVol)
self.assertEqual(COM[1].tolist(), gl['centreOfMassABS'].tolist())
self.assertEqual((COM[1] - numpy.array((1, 1, 1))).tolist(), gl['centreOfMassREL'].tolist())
self.assertEqual(gl['volumeInitial'], 1)
self.assertEqual(COM[1].tolist(), gottenLabel['centreOfMassABS'].tolist())
self.assertEqual((COM[1] - numpy.array((1, 1, 1))).tolist(), gottenLabel['centreOfMassREL'].tolist())
self.assertEqual(gottenLabel['volumeInitial'], 1)
BB = spam.label.boundingBoxes(threeCubedLabelVol)
# With extract cube
gl = spam.label.getLabel(threeCubedLabelVol, 1, boundingBoxes=BB, centresOfMass=COM, extractCube=True)
self.assertEqual(gl['sliceCube'][0].start, 1)
self.assertEqual(gl['sliceCube'][0].stop, 2)
self.assertEqual(gl['sliceCube'][1].start, 1)
self.assertEqual(gl['sliceCube'][1].stop, 2)
self.assertEqual(gl['sliceCube'][2].start, 1)
self.assertEqual(gl['sliceCube'][2].stop, 2)
# print gl['slice']
gottenLabel = spam.label.getLabel(threeCubedLabelVol, 1, boundingBoxes=BB, centresOfMass=COM, extractCube=True)
self.assertEqual(gottenLabel['sliceCube'][0].start, 1)
self.assertEqual(gottenLabel['sliceCube'][0].stop, 2)
self.assertEqual(gottenLabel['sliceCube'][1].start, 1)
self.assertEqual(gottenLabel['sliceCube'][1].stop, 2)
self.assertEqual(gottenLabel['sliceCube'][2].start, 1)
self.assertEqual(gottenLabel['boundingBoxCube'][0], 1)
self.assertEqual(gottenLabel['boundingBoxCube'][1], 1)
self.assertEqual(gottenLabel['boundingBoxCube'][2], 1)
self.assertEqual(gottenLabel['boundingBoxCube'][3], 1)
self.assertEqual(gottenLabel['boundingBoxCube'][4], 1)
self.assertEqual(gottenLabel['boundingBoxCube'][5], 1)
# print gottenLabel['slice']
COM = spam.label.centresOfMass(threeCubedLabelVol)
self.assertEqual(COM[1].tolist(), gl['centreOfMassABS'].tolist())
self.assertEqual((COM[1] - numpy.array((1, 1, 1))).tolist(), gl['centreOfMassREL'].tolist())
self.assertEqual(gl['volumeInitial'], 1)
self.assertEqual(COM[1].tolist(), gottenLabel['centreOfMassABS'].tolist())
self.assertEqual((COM[1] - numpy.array((1, 1, 1))).tolist(), gottenLabel['centreOfMassREL'].tolist())
self.assertEqual(gottenLabel['volumeInitial'], 1)
# Asking for a label that is not there
gl = spam.label.getLabel(threeCubedLabelVol, 100, boundingBoxes=BB, centresOfMass=COM, extractCube=True)
self.assertEqual(gl, None)
gottenLabel = spam.label.getLabel(threeCubedLabelVol, 100, boundingBoxes=BB, centresOfMass=COM, extractCube=True)
self.assertEqual(gottenLabel, None)
# Just run the remaining cases
gl = spam.label.getLabel(threeCubedLabelVol, 1, boundingBoxes=BB, centresOfMass=COM, extractCube=True, extractCubeSize=10)
gl = spam.label.getLabel(threeCubedLabelVol, 1, boundingBoxes=BB, centresOfMass=COM, extractCube=True, extractCubeSize=0)
gl = spam.label.getLabel(threeCubedLabelVol, 1, boundingBoxes=BB, centresOfMass=COM, margin=1, labelDilate=-2)
gl = spam.label.getLabel(threeCubedLabelVol, 1, boundingBoxes=BB, centresOfMass=COM, margin=2, labelDilate=2)
volRef = gl['volumeInitial']
gottenLabel = spam.label.getLabel(threeCubedLabelVol, 1, boundingBoxes=BB, centresOfMass=COM, extractCube=True, extractCubeSize=10)
gottenLabel = spam.label.getLabel(threeCubedLabelVol, 1, boundingBoxes=BB, centresOfMass=COM, extractCube=True, extractCubeSize=0)
gottenLabel = spam.label.getLabel(threeCubedLabelVol, 1, boundingBoxes=BB, centresOfMass=COM, margin=1, labelDilate=-2)
gottenLabel = spam.label.getLabel(threeCubedLabelVol, 1, boundingBoxes=BB, centresOfMass=COM, margin=2, labelDilate=2)
volRef = gottenLabel['volumeInitial']
# New case for "labelDilateMaskOtherLabels"
threeCubedLabelVol2 = threeCubedLabelVol.copy()
threeCubedLabelVol2[0,1,1] = 2
gl = spam.label.getLabel(threeCubedLabelVol2, 1, boundingBoxes=BB, centresOfMass=COM, margin=2, labelDilate=2, labelDilateMaskOtherLabels=True)
#print(gl['volume'], volRef)
self.assertEqual(gl['volumeDilated'] > volRef, True)
gottenLabel = spam.label.getLabel(threeCubedLabelVol2, 1, boundingBoxes=BB, centresOfMass=COM, margin=2, labelDilate=2, labelDilateMaskOtherLabels=True)
#print(gottenLabel['volume'], volRef)
self.assertEqual(gottenLabel['volumeDilated'] > volRef, True)
def test_labelsOnEdges(self):
labelled = numpy.zeros((5, 5, 5), dtype="<u2")
......
......@@ -44,22 +44,22 @@ class testAll(unittest.TestCase):
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-bin1-registration.tsv")
rm(testFolder+"snow-ref-snow-def-bin2-registration.tsv")
rm(testFolder+"snow-ref-snow-def-discreteDVC.tsv")
rm(testFolder+"snow-ref-snow-def-discreteDVC.vtk")
rm(testFolder+"snow-ref-snow-def.tsv")
rm(testFolder+"snow-ref-snow-def.vtk")
rm(testFolder+"snow-ref-snow-def-deltaPhiNorm.tif")
rm(testFolder+"snow-ref-snow-def-error.tif")
rm(testFolder+"snow-ref-snow-def-iterations.tif")
rm(testFolder+"snow-ref-snow-def-returnStatus.tif")
rm(testFolder+"snow-ref-snow-def-Xdisp.tif")
rm(testFolder+"snow-ref-snow-def-Ydisp.tif")
rm(testFolder+"snow-ref-snow-def-Zdisp.tif")
rm(testFolder+"snow-ref-snow-def-strain-Geers.tsv")
rm(testFolder+"snow-ref-snow-def-strain-Q8.tsv")
rm(testFolder+"snow-ref-snow-def-strain-Q8.vtk")
rm(testFolder+"snow-def-snow-ref-bin1-registration.tsv")
rm(testFolder+"snow-def-snow-ref-bin2-registration.tsv")
rm(testFolder+"snow-def-snow-ref-discreteDVC.tsv")
rm(testFolder+"snow-def-snow-ref-discreteDVC.vtk")
rm(testFolder+"snow-def-snow-ref.tsv")
rm(testFolder+"snow-def-snow-ref.vtk")
rm(testFolder+"snow-def-snow-ref-deltaPhiNorm.tif")
rm(testFolder+"snow-def-snow-ref-error.tif")
rm(testFolder+"snow-def-snow-ref-iterations.tif")
rm(testFolder+"snow-def-snow-ref-returnStatus.tif")
rm(testFolder+"snow-def-snow-ref-Xdisp.tif")
rm(testFolder+"snow-def-snow-ref-Ydisp.tif")
rm(testFolder+"snow-def-snow-ref-Zdisp.tif")
rm(testFolder+"snow-def-snow-ref-strain-Geers.tsv")
rm(testFolder+"snow-def-snow-ref-strain-Q8.tsv")
rm(testFolder+"snow-def-snow-ref-strain-Q8.vtk")
rm(testFolder+"test-strain.vtk")
rm(testFolder+"test.vtk")
rm(testFolder+"Step0-Step1-discreteDVC.tsv")
......@@ -68,7 +68,7 @@ class testAll(unittest.TestCase):
rm(testFolder+"Step1.tif")
rm(testFolder+"Step2.tif")
rm(testFolder+"Lab0-displaced.tif")
rm(testFolder+"Step0-Step1-bin2-registration.tsv")
rm(testFolder+"Step0-Step1-bin1-registration.tsv")
rm(testFolder+"Step0-Step1-discreteDVC.vtk")
rm(testFolder+"Step0-Step2-discreteDVC.tsv")
rm(testFolder+"Step0-Step2-discreteDVC.vtk")
......@@ -76,6 +76,9 @@ class testAll(unittest.TestCase):
rm(testFolder+"zoom.tsv")
rm(testFolder+"snow-ref-bin.tif")
rm(testFolder+"snow-ref-watershed.tif")
rm(testFolder+"extreme-im1.tif")
rm(testFolder+"extreme-im1lab.tif")
rm(testFolder+"extreme-im2.tif")
pass
except OSError:
pass
......@@ -87,6 +90,8 @@ class testAll(unittest.TestCase):
if not os.path.isdir(testFolder):
os.mkdir(testFolder)
# We'll do def -> ref to avoid interpolation errors
# Run on snow data, with bin2 registration, without output
# load 3D image
snowRef = spam.datasets.loadSnow()
......@@ -105,15 +110,15 @@ class testAll(unittest.TestCase):
exitCode = subprocess.call(["spam-ldic", "--help"], stdout=FNULL, stderr=FNULL)
self.assertEqual(exitCode, 0)
# Just run a simple DVC with no outputs
## 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)
#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", "-reg", "-regbb", "1", "-hws", "10", '-ns', '15', testFolder + "snow-ref.tif", testFolder + "snow-def.tif", "-od", testFolder + "", "-tsv", "-tif", '-vtk'])
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
......@@ -121,8 +126,8 @@ class testAll(unittest.TestCase):
TSV = spam.helpers.readCorrelationTSV(testFolder + "snow-ref-snow-def-bin1-registration.tsv")
transformation = spam.deformation.decomposePhi(TSV['PhiField'][0])
for i in range(3):
self.assertAlmostEqual(refTranslation[i], transformation['t'][i], places=2)
self.assertAlmostEqual(refRotation[i], transformation['r'][i], places=2)
self.assertAlmostEqual(refTranslation[i], transformation['t'][i], places=1)
self.assertAlmostEqual(refRotation[i], transformation['r'][i], places=1)
# Check that saved TIFF corresponds to the TSV
TSV = spam.helpers.readCorrelationTSV(testFolder + "snow-ref-snow-def.tsv")
......@@ -140,11 +145,11 @@ class testAll(unittest.TestCase):
# There's Z-rotation so translations need to be decomposed properly... just check Z
self.assertAlmostEqual(refTranslation[0], numpy.mean(Phis[converged,0,-1]), places=1)
# Check only rotations, with a rotation there will be a non-trivial displacement field
for i in range(3):
#self.assertAlmostEqual(refTranslation[i], numpy.mean(Phis[converged,i,-1]), places=1)
self.assertAlmostEqual(refRotation[i], numpy.mean(decomposedPhisConverged['r'][converged, i]), places=1)
#######################################################
## OK onto the regularStrain
#######################################################
......@@ -160,7 +165,7 @@ class testAll(unittest.TestCase):
# Now read output
strain = spam.helpers.readStrainTSV("snow-ref-snow-def-strain-Geers.tsv")
dims = [strain['fieldDims'][0], strain['fieldDims'][1], strain['fieldDims'][2], 3]
# Check each rotation component while having a 2-pixel crop of the boundaries
for i in range(3):
self.assertAlmostEqual(numpy.mean(strain['r'].reshape(dims)[2:-2, 2:-2, 2:-2, i]),
......@@ -188,7 +193,7 @@ class testAll(unittest.TestCase):
VTK = spam.helpers.readStructuredVTK("snow-ref-snow-def-strain-Q8.vtk")
for component in ["U", "e", "vol", "volss", "dev", "devss"]:
self.assertAlmostEqual(numpy.nanmean(strain[component]), numpy.nanmean(VTK[component]), places=3)
self.assertAlmostEqual(numpy.nanmean(strain[component]), numpy.nanmean(VTK[component]), places=2)
def test_ITKwatershed(self):
# make sure it runs the help without error
......@@ -212,7 +217,6 @@ class testAll(unittest.TestCase):
lab = tifffile.imread(testFolder + "snow-ref-watershed.tif")
self.assertEqual(lab.max()>10, True)
def test_gdic(self):
# make sure it runs the help without error
exitCode = subprocess.call(["spam-gdic", "--help"], stdout=FNULL, stderr=FNULL)
......@@ -257,11 +261,11 @@ class testAll(unittest.TestCase):
newline='\n',
comments='',
header="Label\tZpos\tYpos\tXpos\t" +
"Zdisp\tYdisp\tXdisp\t" +
"Fzz\tFzy\tFzx\t" +
"Fyz\tFyy\tFyx\t" +
"Fxz\tFxy\tFxx\t" +
"PSCC\terror\titerations\treturnStatus\tdeltaPhiNorm")
"Zdisp\tYdisp\tXdisp\t" +
"Fzz\tFzy\tFzx\t" +
"Fyz\tFyy\tFyx\t" +
"Fxz\tFxy\tFxx\t" +
"PSCC\terror\titerations\treturnStatus\tdeltaPhiNorm")
### Run the script
exitCode = subprocess.call(["spam-deformImageFromField", "-pf", testFolder+"rotation.tsv", testFolder+"snow-ref.tif",
......@@ -664,6 +668,33 @@ class testAll(unittest.TestCase):
for u in range(3):
self.assertAlmostEqual(TSV['PhiField'][i][u,-1], TSV2['PhiField'][i][u,-1], places=1)
#######################################################
### 2020-09-05 EA and OS: New "extreme" test with particles touching the boundaries
#######################################################
im1 = spam.kalisphera.makeBlurryNoisySphere([100,100,100], [[50,50,50], [50,10,10]], [10, 10], 0.5, 0.05, background=0.0, foreground=1.0)
im2 = spam.kalisphera.makeBlurryNoisySphere([100,100,100], [[54,54,54], [54,14,14]], [10, 10], 0.5, 0.05, background=0.0, foreground=1.0)
im1lab = scipy.ndimage.label(im1>0.5)[0]
tifffile.imsave(testFolder + "extreme-im1.tif", im1)
tifffile.imsave(testFolder + "extreme-im1lab.tif", im1lab)
tifffile.imsave(testFolder + "extreme-im2.tif", im2)
# Just run a simple DVC
exitCode = subprocess.call(["spam-ddic",
testFolder + "extreme-im1.tif",
testFolder + "extreme-im1lab.tif",
testFolder + "extreme-im2.tif",
"-regbb", "4", "-regbe", "2",
"-lcm", "10", "-ld", "2"])
self.assertEqual(exitCode, 0)
TSVextreme = spam.helpers.readCorrelationTSV(testFolder + "extreme-im1-extreme-im2-discreteDVC.tsv")
#self.assertAlmostEqual(TSV['PhiField'][i][u,-1], TSV2['PhiField'][i][u,-1], places=1)
self.assertEqual(list(TSVextreme['returnStatus'][1:]), [2, 2])
for axis in range(3):
for label in [1,2]:
self.assertAlmostEqual(numpy.abs(TSVextreme['PhiField'][label, axis, -1] - 4.0), 0.0, places=1)
def test_discreteStrain(self):
# make sure it runs the help without error
exitCode = subprocess.call(["spam-discreteStrain", "--help"],
......
......@@ -508,7 +508,7 @@ def register(im1, im2, im1mask=None, PhiInit=None, PhiRigid=False, PhiInitBinRat
plt.subplot(3,3,9)
plt.title('im1-im2def X-slice')
plt.imshow(numpy.subtract(im1crop, im2def[crop2])[:, :, im1crop.shape[2]//2], cmap='coolwarm', vmin=vmin, vmax=vmax)
plt.pause(0.1)
plt.pause(0.01)
iterations += 1
......
......@@ -160,31 +160,38 @@ def pixelSearchOnGrid(im1, im2, nodePositions, halfWindowSize, searchRange, PhiF
halfWindowSize[1] - searchRangeForThisNode['yRange'][0],
halfWindowSize[2] - searchRangeForThisNode['xRange'][0]]
# Prepare slice for im1 -- this is always in the same place
subVolSlice1 = (slice(int(nodePositions[nodeNumber, 0] - halfWindowSize[0]),
int(nodePositions[nodeNumber, 0] + halfWindowSize[0] + 1)),
slice(int(nodePositions[nodeNumber, 1] - halfWindowSize[1]),
int(nodePositions[nodeNumber, 1] + halfWindowSize[1] + 1)),
slice(int(nodePositions[nodeNumber, 2] - halfWindowSize[2]),
int(nodePositions[nodeNumber, 2] + halfWindowSize[2] + 1)))
# Extract it...
imagette1 = im1[subVolSlice1]
# 2020-09-25 OS and EA: Prepare startStop array for imagette 1 to be extracted with new slicePadded
## Prepare slice for im1 -- this is always in the same place
#subVolSlice1 = (slice(int(nodePositions[nodeNumber, 0] - halfWindowSize[0]), int(nodePositions[nodeNumber, 0] + halfWindowSize[0] + 1)),
#slice(int(nodePositions[nodeNumber, 1] - halfWindowSize[1]), int(nodePositions[nodeNumber, 1] + halfWindowSize[1] + 1)),
#slice(int(nodePositions[nodeNumber, 2] - halfWindowSize[2]), int(nodePositions[nodeNumber, 2] + halfWindowSize[2] + 1)))
## Extract it...
#imagette1 = im1[subVolSlice1]
startStopIm1 = [int(nodePositions[nodeNumber, 0] - halfWindowSize[0]), int(nodePositions[nodeNumber, 0] + halfWindowSize[0] + 1),
int(nodePositions[nodeNumber, 1] - halfWindowSize[1]), int(nodePositions[nodeNumber, 1] + halfWindowSize[1] + 1),
int(nodePositions[nodeNumber, 2] - halfWindowSize[2]), int(nodePositions[nodeNumber, 2] + halfWindowSize[2] + 1)]
imagette1 = spam.helpers.slicePadded(im1, startStopIm1)
# Make sure imagette is not 0-dimensional in any dimension
if numpy.all(numpy.array(imagette1.shape) > 0):
if numpy.nanmean(imagette1) > greyThreshold[0] and numpy.nanmean(imagette1) < greyThreshold[1] and len(imagette1.ravel()) > minMaskVolume:
# Slice for image 2
# TODO: Check we extraced what we asked for -- could put displacements off at the boundaries
# could even give a bad return status for this condition...
subVolSlice2 = (slice(int(nodePositions[nodeNumber, 0] - halfWindowSize[0] + searchRangeForThisNode['zRange'][0]),
int(nodePositions[nodeNumber, 0] + halfWindowSize[0] + searchRangeForThisNode['zRange'][1] + 1)),
slice(int(nodePositions[nodeNumber, 1] - halfWindowSize[1] + searchRangeForThisNode['yRange'][0]),
int(nodePositions[nodeNumber, 1] + halfWindowSize[1] + searchRangeForThisNode['yRange'][1] + 1)),
slice(int(nodePositions[nodeNumber, 2] - halfWindowSize[2] + searchRangeForThisNode['xRange'][0]),
int(nodePositions[nodeNumber, 2] + halfWindowSize[2] + searchRangeForThisNode['xRange'][1] + 1)))
# Extract it...
imagette2 = im2[subVolSlice2]
## 2020-09-25 OS and EA: Prepare startStop array for imagette 2 to be extracted with new slicePadded
#subVolSlice2 = (slice(int(nodePositions[nodeNumber, 0] - halfWindowSize[0] + searchRangeForThisNode['zRange'][0]),
#int(nodePositions[nodeNumber, 0] + halfWindowSize[0] + searchRangeForThisNode['zRange'][1] + 1)),
#slice(int(nodePositions[nodeNumber, 1] - halfWindowSize[1] + searchRangeForThisNode['yRange'][0]),
#int(nodePositions[nodeNumber, 1] + halfWindowSize[1] + searchRangeForThisNode['yRange'][1] + 1)),
#slice(int(nodePositions[nodeNumber, 2] - halfWindowSize[2] + searchRangeForThisNode['xRange'][0]),
#int(nodePositions[nodeNumber, 2] + halfWindowSize[2] + searchRangeForThisNode['xRange'][1] + 1)))
## Extract it...
#imagette2 = im2[subVolSlice2]
startStopIm2 = [int(nodePositions[nodeNumber, 0] - halfWindowSize[0] + searchRangeForThisNode['zRange'][0]),
int(nodePositions[nodeNumber, 0] + halfWindowSize[0] + searchRangeForThisNode['zRange'][1] + 1),
int(nodePositions[nodeNumber, 1] - halfWindowSize[1] + searchRangeForThisNode['yRange'][0]),
int(nodePositions[nodeNumber, 1] + halfWindowSize[1] + searchRangeForThisNode['yRange'][1] + 1),
int(nodePositions[nodeNumber, 2] - halfWindowSize[2] + searchRangeForThisNode['xRange'][0]),
int(nodePositions[nodeNumber, 2] + halfWindowSize[2] + searchRangeForThisNode['xRange'][1] + 1)]
imagette2 = spam.helpers.slicePadded(im2, startStopIm2)
# Failed minMaskVolume or greylevel condition
else:
......@@ -198,7 +205,8 @@ def pixelSearchOnGrid(im1, im2, nodePositions, halfWindowSize, searchRange, PhiF
imagette1 = None
imagette2 = None
return {'imagette1': imagette1, 'imagette2': imagette2,
return {'imagette1': imagette1,
'imagette2': imagette2,
'returnStatus': returnStatus,
'initialDisplacement': initialDisplacement,
'searchRangeForThisNode': searchRangeForThisNode,
......@@ -436,18 +444,24 @@ def registerOnGrid(im1, im2, nodePositions, halfWindowSize, PhiField=None, margi
nodeDisplacement = numpy.round(PhiField[nodeNumber][0:3, -1])
# prepare slice for imagette 1 -- this is fixed +1 for the margin which is always 1
subVolSlice1 = (slice(int(nodePositions[nodeNumber, 0] - halfWindowSize[0] - 1), int(nodePositions[nodeNumber, 0] + halfWindowSize[0] + 1 + 1)),
slice(int(nodePositions[nodeNumber, 1] - halfWindowSize[1] - 1), int(nodePositions[nodeNumber, 1] + halfWindowSize[1] + 1 + 1)),
slice(int(nodePositions[nodeNumber, 2] - halfWindowSize[2] - 1), int(nodePositions[nodeNumber, 2] + halfWindowSize[2] + 1 + 1)))
imagette1 = im1[subVolSlice1].copy()
# 2020-09-25 OS and EA: Prepare startStop array for imagette 1 to be extracted with new slicePadded
#subVolSlice1 = (slice(int(nodePositions[nodeNumber, 0] - halfWindowSize[0] - 1), int(nodePositions[nodeNumber, 0] + halfWindowSize[0] + 1 + 1)),
#slice(int(nodePositions[nodeNumber, 1] - halfWindowSize[1] - 1), int(nodePositions[nodeNumber, 1] + halfWindowSize[1] + 1 + 1)),
#slice(int(nodePositions[nodeNumber, 2] - halfWindowSize[2] - 1), int(nodePositions[nodeNumber, 2] + halfWindowSize[2] + 1 + 1)))
#imagette1 = im1[subVolSlice1].copy()
startStopIm1 = [int(nodePositions[nodeNumber, 0] - halfWindowSize[0] - 1), int(nodePositions[nodeNumber, 0] + halfWindowSize[0] + 1 + 1),
int(nodePositions[nodeNumber, 1] - halfWindowSize[1] - 1), int(nodePositions[nodeNumber, 1] + halfWindowSize[1] + 1 + 1),
int(nodePositions[nodeNumber, 2] - halfWindowSize[2] - 1), int(nodePositions[nodeNumber, 2] + halfWindowSize[2] + 1 + 1)]
imagette1, imagette1mask = spam.helpers.slicePadded(im1, startStopIm1, createMask=True)
# If there is a mask defined for im1, compute the coverage of this correlation window
# Check Mask volume condition
if im1mask is not None:
imagette1mask = im1mask[subVolSlice1]
maskVolCondition = imagette1mask.sum() > minMaskVolume
else:
maskVolCondition = True
# AND masks together, the one coming from the im1mask -- i.e., allowed correlation zones
# and the imagette1sliceMask indicating whether we've gone outside the volume
imagette1mask = numpy.logical_and(imagette1mask, spam.helpers.slicePadded(im1mask, startStopIm1))
maskVolCondition = imagette1mask.sum() > minMaskVolume
# Make sure imagette is not 0-dimensional in any dimension
if numpy.all(numpy.array(imagette1.shape) > 0):
......@@ -457,11 +471,18 @@ def registerOnGrid(im1, im2, nodePositions, halfWindowSize, PhiField=None, margi
# Prepare im2 imagette
# start from nodePosition in im1 and move it according to the nodeDisplacement