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 3713ab91 authored by Edward Andò's avatar Edward Andò
Browse files

slicePadded in spam-ddic + tests

parent 56144098
......@@ -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,7 +364,11 @@ if mpiRank == boss or not mpi:
else:
maskette1 = None
maskette1vol = numpy.inf # must pass "if" below
imagette1 = im1[getLabel['slice']].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'):
......@@ -380,16 +384,15 @@ 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)]
# Extract it...
imagette2 = im2[subVolSlice2]
imagette2 = spam.helpers.slicePadded(im2, startStopIm2)
# point in im2 that we are searching around
searchCentre = (numpy.array(imagette2.shape, dtype='<f4') - 1) / 2.0 + labelDisplacementInt - middleOfSearchRange
......@@ -413,6 +416,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:
......@@ -461,12 +465,12 @@ if mpiRank == boss or not mpi:
# 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(boundingBoxes[label][0] - args.LABEL_CORRELATE_MARGIN - max(labelDilateCurrent, 0) + labelDisplacementInt[0] ),
int(boundingBoxes[label][1] + args.LABEL_CORRELATE_MARGIN + max(labelDilateCurrent, 0) + labelDisplacementInt[0] + 1),
int(boundingBoxes[label][2] - args.LABEL_CORRELATE_MARGIN - max(labelDilateCurrent, 0) + labelDisplacementInt[1] ),
int(boundingBoxes[label][3] + args.LABEL_CORRELATE_MARGIN + max(labelDilateCurrent, 0) + labelDisplacementInt[1] + 1),
int(boundingBoxes[label][4] - args.LABEL_CORRELATE_MARGIN - max(labelDilateCurrent, 0) + labelDisplacementInt[2] ),
int(boundingBoxes[label][5] + args.LABEL_CORRELATE_MARGIN + max(labelDilateCurrent, 0) + labelDisplacementInt[2] + 1)]
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)]
imagette2 = spam.helpers.slicePadded(im2, startStopIm2)
......
......@@ -82,7 +82,7 @@ class testAll(unittest.TestCase):
pass
def test_ldic_and_regularStrain(self):
def _test_ldic_and_regularStrain(self):
# Make test directory
if not os.path.isdir(testFolder):
os.mkdir(testFolder)
......@@ -160,7 +160,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]),
......@@ -190,7 +190,7 @@ class testAll(unittest.TestCase):
for component in ["U", "e", "vol", "volss", "dev", "devss"]:
self.assertAlmostEqual(numpy.nanmean(strain[component]), numpy.nanmean(VTK[component]), places=3)
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)
......@@ -212,13 +212,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)
......@@ -664,7 +663,34 @@ class testAll(unittest.TestCase):
for u in range(3):
self.assertAlmostEqual(TSV['PhiField'][i][u,-1], TSV2['PhiField'][i][u,-1], places=1)
def test_discreteStrain(self):
#######################################################
### 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"],
stdout=FNULL,
......
......@@ -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.startStopIm1(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,16 @@ def registerOnGrid(im1, im2, nodePositions, halfWindowSize, PhiField=None, margi
# Prepare im2 imagette
# start from nodePosition in im1 and move it according to the nodeDisplacement
# and move it + 1 for the margin which is always 1
subVolSlice2 = (slice(int(nodePositions[nodeNumber, 0] - halfWindowSize[0] + nodeDisplacement[0] - margin[0] - 1), int(nodePositions[nodeNumber, 0] + halfWindowSize[0] + nodeDisplacement[0] + margin[0] + 1 + 1)),
slice(int(nodePositions[nodeNumber, 1] - halfWindowSize[1] + nodeDisplacement[1] - margin[1] - 1), int(nodePositions[nodeNumber, 1] + halfWindowSize[1] + nodeDisplacement[1] + margin[1] + 1 + 1)),
slice(int(nodePositions[nodeNumber, 2] - halfWindowSize[2] + nodeDisplacement[2] - margin[2] - 1), int(nodePositions[nodeNumber, 2] + halfWindowSize[2] + nodeDisplacement[2] + margin[2] + 1 + 1)))
# Extract it
imagette2 = im2[subVolSlice2].copy()
# 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] + nodeDisplacement[0] - margin[0] - 1), int(nodePositions[nodeNumber, 0] + halfWindowSize[0] + nodeDisplacement[0] + margin[0] + 1 + 1)),
#slice(int(nodePositions[nodeNumber, 1] - halfWindowSize[1] + nodeDisplacement[1] - margin[1] - 1), int(nodePositions[nodeNumber, 1] + halfWindowSize[1] + nodeDisplacement[1] + margin[1] + 1 + 1)),
#slice(int(nodePositions[nodeNumber, 2] - halfWindowSize[2] + nodeDisplacement[2] - margin[2] - 1), int(nodePositions[nodeNumber, 2] + halfWindowSize[2] + nodeDisplacement[2] + margin[2] + 1 + 1)))
## Extract it
#imagette2 = im2[subVolSlice2].copy()
startStopIm2 = [int(nodePositions[nodeNumber, 0] - halfWindowSize[0] + nodeDisplacement[0] - margin[0] - 1), int(nodePositions[nodeNumber, 0] + halfWindowSize[0] + nodeDisplacement[0] + margin[0] + 1 + 1),
int(nodePositions[nodeNumber, 1] - halfWindowSize[1] + nodeDisplacement[1] - margin[1] - 1), int(nodePositions[nodeNumber, 1] + halfWindowSize[1] + nodeDisplacement[1] + margin[1] + 1 + 1),
int(nodePositions[nodeNumber, 2] - halfWindowSize[2] + nodeDisplacement[2] - margin[2] - 1), int(nodePositions[nodeNumber, 2] + halfWindowSize[2] + nodeDisplacement[2] + margin[2] + 1 + 1)]
imagette2 = spam.helpers.slicePadded(im2, startStopIm2)
# Extract initial F for correlation, remove int() part of displacement since it's already used to extract imagette2
PhiInit = PhiField[nodeNumber].copy()
......
......@@ -481,11 +481,19 @@ def slicePadded(im, startStop, createMask = False):
mask : 3D numpy array of bools
The 3D mask, only returned if createMask is True
"""
startStop = numpy.array(startStop).astype(numpy.int)
assert (startStop[1]>startStop[0]), "spam.helpers.slicePadded(): Zmax should be bigger than Zmin"
assert (startStop[3]>startStop[2]), "spam.helpers.slicePadded(): Ymax should be bigger than Ymin"
assert (startStop[5]>startStop[4]), "spam.helpers.slicePadded(): Xmax should be bigger than Xmin"
imSliced = numpy.zeros((startStop[1]-startStop[0],
startStop[3]-startStop[2],
startStop[5]-startStop[4]), dtype=im.dtype)
start = numpy.array([startStop[0], startStop[2], startStop[4]])
stop = numpy.array([startStop[1], startStop[3], startStop[5]])
start = numpy.array([startStop[0], startStop[2], startStop[4]])
stop = numpy.array([startStop[1], startStop[3], startStop[5]])
startOffset = numpy.array([max(0, start[0]), max(0, start[1]), max(0, start[2])])
stopOffset = numpy.array([min(im.shape[0], stop[0]), min(im.shape[1], stop[1]), min(im.shape[2], stop[2])])
startLocal = startOffset - start
......
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