Commit 69f8675e authored by Edward Andò's avatar Edward Andò
Browse files

new ddic script: lc options instead of gp; gradient update; debug mode;...

new ddic script: lc options instead of gp; gradient update; debug mode; applyPhi nicely for pixel search
parent f94ea529
Pipeline #47832 failed with stages
in 11 minutes and 45 seconds
...@@ -345,7 +345,7 @@ if mpiRank == boss or not mpi: ...@@ -345,7 +345,7 @@ if mpiRank == boss or not mpi:
extractCube=False, extractCube=False,
boundingBoxes=boundingBoxes, boundingBoxes=boundingBoxes,
centresOfMass=centresOfMass, centresOfMass=centresOfMass,
margin=labelDilateCurrent + args.GRID_POINT_MARGIN, margin=labelDilateCurrent + args.LABEL_CORRELATE_MARGIN,
maskOtherLabels=args.MASK, maskOtherLabels=args.MASK,
labelDilate=labelDilateCurrent) labelDilate=labelDilateCurrent)
...@@ -382,12 +382,12 @@ if mpiRank == boss or not mpi: ...@@ -382,12 +382,12 @@ if mpiRank == boss or not mpi:
(searchRangeForThisLabel['xRange'][0] + searchRangeForThisLabel['xRange'][1])/2] ).astype(int) (searchRangeForThisLabel['xRange'][0] + searchRangeForThisLabel['xRange'][1])/2] ).astype(int)
# Slice for image 2 # Slice for image 2
subVolSlice2 = (slice(max(int(boundingBoxes[label][0] - max(labelDilateCurrent, 0) - args.GRID_POINT_MARGIN + searchRangeForThisLabel['zRange'][0]), 0), 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.GRID_POINT_MARGIN + searchRangeForThisLabel['zRange'][1] + 1), im1.shape[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.GRID_POINT_MARGIN + searchRangeForThisLabel['yRange'][0]), 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.GRID_POINT_MARGIN + searchRangeForThisLabel['yRange'][1] + 1), im1.shape[1])), 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.GRID_POINT_MARGIN + searchRangeForThisLabel['xRange'][0]), 0), 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.GRID_POINT_MARGIN + searchRangeForThisLabel['xRange'][1] + 1), im1.shape[2]))) min(int(boundingBoxes[label][5] + max(labelDilateCurrent, 0) + args.LABEL_CORRELATE_MARGIN + searchRangeForThisLabel['xRange'][1] + 1), im1.shape[2])))
# Extract it... # Extract it...
imagette2 = im2[subVolSlice2] imagette2 = im2[subVolSlice2]
...@@ -403,12 +403,12 @@ if mpiRank == boss or not mpi: ...@@ -403,12 +403,12 @@ if mpiRank == boss or not mpi:
# Deform imagette1, using relative COM as point of application # Deform imagette1, using relative COM as point of application
imagette1def = spam.DIC.applyPhi(imagette1, PhiNoDisp, PhiPoint=getLabel['centreOfMassREL']) imagette1def = spam.DIC.applyPhi(imagette1, PhiNoDisp, PhiPoint=getLabel['centreOfMassREL'])
imagette1defCrop = imagette1def[args.GRID_POINT_MARGIN:-args.GRID_POINT_MARGIN, #imagette1defCrop = imagette1def[args.LABEL_CORRELATE_MARGIN:-args.LABEL_CORRELATE_MARGIN,
args.GRID_POINT_MARGIN:-args.GRID_POINT_MARGIN, # This is to remove edge artefacts of applyPhi #args.LABEL_CORRELATE_MARGIN:-args.LABEL_CORRELATE_MARGIN, # This is to remove edge artefacts of applyPhi
args.GRID_POINT_MARGIN:-args.GRID_POINT_MARGIN] #args.LABEL_CORRELATE_MARGIN:-args.LABEL_CORRELATE_MARGIN]
imagette1crop = imagette1[args.GRID_POINT_MARGIN:-args.GRID_POINT_MARGIN, #imagette1crop = imagette1[args.LABEL_CORRELATE_MARGIN:-args.LABEL_CORRELATE_MARGIN,
args.GRID_POINT_MARGIN:-args.GRID_POINT_MARGIN, # This is to remove edge artefacts of applyPhi #args.LABEL_CORRELATE_MARGIN:-args.LABEL_CORRELATE_MARGIN, # This is to remove edge artefacts of applyPhi
args.GRID_POINT_MARGIN:-args.GRID_POINT_MARGIN] #args.LABEL_CORRELATE_MARGIN:-args.LABEL_CORRELATE_MARGIN]
if args.MASK: if args.MASK:
imagette1toCorrelate = imagette1def.copy() imagette1toCorrelate = imagette1def.copy()
...@@ -456,15 +456,15 @@ if mpiRank == boss or not mpi: ...@@ -456,15 +456,15 @@ if mpiRank == boss or not mpi:
print("Phi after pixel search:\n", PhiField[label]) print("Phi after pixel search:\n", PhiField[label])
if args.GRID_POINT: if args.LABEL_CORRELATE:
labelDisplacementInt = numpy.round(PhiField[label][0:3, -1]) labelDisplacementInt = numpy.round(PhiField[label][0:3, -1])
slicette2 = (slice(max(int(boundingBoxes[label][0] - args.GRID_POINT_MARGIN - max(labelDilateCurrent, 0) + labelDisplacementInt[0] ), 0 ), slicette2 = (slice(max(int(boundingBoxes[label][0] - args.LABEL_CORRELATE_MARGIN - max(labelDilateCurrent, 0) + labelDisplacementInt[0] ), 0 ),
min(int(boundingBoxes[label][1] + args.GRID_POINT_MARGIN + max(labelDilateCurrent, 0) + labelDisplacementInt[0] + 1), im1.shape[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.GRID_POINT_MARGIN - max(labelDilateCurrent, 0) + labelDisplacementInt[1] ), 0 ), slice(max(int(boundingBoxes[label][2] - args.LABEL_CORRELATE_MARGIN - max(labelDilateCurrent, 0) + labelDisplacementInt[1] ), 0 ),
min(int(boundingBoxes[label][3] + args.GRID_POINT_MARGIN + max(labelDilateCurrent, 0) + labelDisplacementInt[1] + 1), im1.shape[1])), 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.GRID_POINT_MARGIN - max(labelDilateCurrent, 0) + labelDisplacementInt[2] ), 0 ), slice(max(int(boundingBoxes[label][4] - args.LABEL_CORRELATE_MARGIN - max(labelDilateCurrent, 0) + labelDisplacementInt[2] ), 0 ),
min(int(boundingBoxes[label][5] + args.GRID_POINT_MARGIN + max(labelDilateCurrent, 0) + labelDisplacementInt[2] + 1), im1.shape[2]))) min(int(boundingBoxes[label][5] + args.LABEL_CORRELATE_MARGIN + max(labelDilateCurrent, 0) + labelDisplacementInt[2] + 1), im1.shape[2])))
#print(slicette1) #print(slicette1)
#print(slicette2) #print(slicette2)
...@@ -485,9 +485,11 @@ if mpiRank == boss or not mpi: ...@@ -485,9 +485,11 @@ if mpiRank == boss or not mpi:
im1mask=maskette1, im1mask=maskette1,
margin=1, margin=1,
PhiInit=PhiTemp, PhiInit=PhiTemp,
maxIterations=args.GRID_POINT_MAX_ITERATIONS, PhiRigid=args.LABEL_CORRELATE_RIGID,
deltaPhiMin=args.GRID_POINT_MIN_PHI_CHANGE, updateGradient=args.LABEL_CORRELATE_UPDATE_GRADIENT,
interpolationOrder=args.GRID_POINT_INTERPOLATION_ORDER, maxIterations=args.LABEL_CORRELATE_MAX_ITERATIONS,
deltaPhiMin=args.LABEL_CORRELATE_MIN_PHI_CHANGE,
interpolationOrder=args.LABEL_CORRELATE_INTERPOLATION_ORDER,
verbose=args.DEBUG, verbose=args.DEBUG,
imShowProgress='Z' if args.DEBUG else None) imShowProgress='Z' if args.DEBUG else None)
writeReturns = True writeReturns = True
...@@ -499,10 +501,12 @@ if mpiRank == boss or not mpi: ...@@ -499,10 +501,12 @@ if mpiRank == boss or not mpi:
'im2': im2[slicette2], 'im2': im2[slicette2],
'im1mask': maskette1, 'im1mask': maskette1,
'PhiInit': PhiTemp, 'PhiInit': PhiTemp,
'PhiRigid': args.LABEL_CORRELATE_RIGID,
'updateGradient': args.LABEL_CORRELATE_UPDATE_GRADIENT,
'margin': 1, # see top of this file for compensation 'margin': 1, # see top of this file for compensation
'maxIterations': args.GRID_POINT_MAX_ITERATIONS, 'maxIterations': args.LABEL_CORRELATE_MAX_ITERATIONS,
'deltaPhiMin': args.GRID_POINT_MIN_PHI_CHANGE, 'deltaPhiMin': args.LABEL_CORRELATE_MIN_PHI_CHANGE,
'interpolationOrder': args.GRID_POINT_INTERPOLATION_ORDER, 'interpolationOrder': args.LABEL_CORRELATE_INTERPOLATION_ORDER,
'labelDisplacementInt': labelDisplacementInt, 'labelDisplacementInt': labelDisplacementInt,
'writeReturns': writeReturns, 'writeReturns': writeReturns,
'labelDilate': labelDilateCurrent 'labelDilate': labelDilateCurrent
...@@ -519,7 +523,7 @@ if mpiRank == boss or not mpi: ...@@ -519,7 +523,7 @@ if mpiRank == boss or not mpi:
print("\t\tBoss: Fail imDiff", imagette2imagette1sizeDiff) print("\t\tBoss: Fail imDiff", imagette2imagette1sizeDiff)
grainOK = True grainOK = True
else: # No args.GRID_POINT else: # No args.LABEL_CORRELATE
finishedLabels += 1 finishedLabels += 1
grainOK = True grainOK = True
...@@ -665,6 +669,8 @@ elif mpi: # We are not the mpi boss, so we are a lukasKanade worker ...@@ -665,6 +669,8 @@ elif mpi: # We are not the mpi boss, so we are a lukasKanade worker
m['im2'], m['im2'],
im1mask=m['im1mask'], im1mask=m['im1mask'],
PhiInit=m['PhiInit'], PhiInit=m['PhiInit'],
PhiRigid=m['PhiRigid'],
updateGradient=m['updateGradient'],
margin=m['margin'], margin=m['margin'],
maxIterations=m['maxIterations'], maxIterations=m['maxIterations'],
deltaPhiMin=m['deltaPhiMin'], deltaPhiMin=m['deltaPhiMin'],
......
...@@ -262,6 +262,17 @@ class TestFunctionDVC(unittest.TestCase): ...@@ -262,6 +262,17 @@ class TestFunctionDVC(unittest.TestCase):
self.assertAlmostEqual(numpy.array(returns7['transformation']["t"][i]), 0, places=1) self.assertAlmostEqual(numpy.array(returns7['transformation']["t"][i]), 0, places=1)
self.assertAlmostEqual(numpy.array(returns7['transformation']["r"][i]) - [rot, 0, 0][i], 0, places=1) self.assertAlmostEqual(numpy.array(returns7['transformation']["r"][i]) - [rot, 0, 0][i], 0, places=1)
# CASE 7b: lasrge rot and only rigid
returns7b = spam.DIC.register(im, imDef, margin=10, PhiInit=PhiGuess, PhiRigid=True)
self.assertEqual(returns7b['returnStatus'], 2)
returns7b['transformation'] = spam.deformation.decomposePhi(returns7b['Phi'])
for i in range(3):
self.assertAlmostEqual(numpy.array(returns7b['transformation']["t"][i]), 0, places=1)
self.assertAlmostEqual(numpy.array(returns7b['transformation']["r"][i]) - [rot, 0, 0][i], 0, places=1)
# Check that it's rigid!!
self.assertAlmostEqual(numpy.linalg.norm(numpy.array(returns7b['transformation']["U"])-numpy.eye(3)), 0, places=4)
# CASE 8: 3D case with large initial rotation (say 120 deg) with gradient update, should need fewer iterations # CASE 8: 3D case with large initial rotation (say 120 deg) with gradient update, should need fewer iterations
rot = 120 rot = 120
Phi = spam.deformation.computePhi({'r': [rot, 0, 0]}) Phi = spam.deformation.computePhi({'r': [rot, 0, 0]})
...@@ -274,12 +285,12 @@ class TestFunctionDVC(unittest.TestCase): ...@@ -274,12 +285,12 @@ class TestFunctionDVC(unittest.TestCase):
for i in range(3): for i in range(3):
self.assertAlmostEqual(numpy.array(returns8['transformation']["t"][i]), 0, places=1) self.assertAlmostEqual(numpy.array(returns8['transformation']["t"][i]), 0, places=1)
self.assertAlmostEqual(numpy.array(returns8['transformation']["r"][i]) - [rot, 0, 0][i], 0, places=1) self.assertAlmostEqual(numpy.array(returns8['transformation']["r"][i]) - [rot, 0, 0][i], 0, places=1)
#Test for a singular PhiInit #Test for a singular PhiInit
PhiGuess = numpy.array([[1., 2, 3, 0],[2, 4, 6, 0],[3, 6, 9, 0]]) PhiGuess = numpy.array([[1., 2, 3, 0],[2, 4, 6, 0],[3, 6, 9, 0]])
returns8 = spam.DIC.register(im, imDef, margin=10, PhiInit=PhiGuess, verbose=True, updateGradient=True) returns8 = spam.DIC.register(im, imDef, margin=10, PhiInit=PhiGuess, verbose=True, updateGradient=True)
self.assertEqual(returns8['returnStatus'], -2) self.assertEqual(returns8['returnStatus'], -2)
def test_registerMultiScale(self): def test_registerMultiScale(self):
t = {'t': [0, 0, 0], t = {'t': [0, 0, 0],
...@@ -298,7 +309,9 @@ class TestFunctionDVC(unittest.TestCase): ...@@ -298,7 +309,9 @@ class TestFunctionDVC(unittest.TestCase):
for i in range(3): for i in range(3):
self.assertAlmostEqual(t2[c][i] - t[c][i], 0, places=1) self.assertAlmostEqual(t2[c][i] - t[c][i], 0, places=1)
# Check transaltion example with translation guess
# Check translation example with translation guess
t = {'t': [4.0, 0, 0], t = {'t': [4.0, 0, 0],
'r': [0.0, 0, 0]} 'r': [0.0, 0, 0]}
Phi = spam.deformation.computePhi(t) Phi = spam.deformation.computePhi(t)
...@@ -315,13 +328,13 @@ class TestFunctionDVC(unittest.TestCase): ...@@ -315,13 +328,13 @@ class TestFunctionDVC(unittest.TestCase):
for i in range(3): for i in range(3):
self.assertAlmostEqual(t4[c][i] - t[c][i], 0, places=1) self.assertAlmostEqual(t4[c][i] - t[c][i], 0, places=1)
# check graident update # check gradient update
reg5 = spam.DIC.registerMultiscale(im, imDef, 2, margin=8, PhiInit=spam.deformation.computePhi({'t': [5,0,0]}), PhiInitBinRatio=2, updateGradient=True) reg5 = spam.DIC.registerMultiscale(im, imDef, 2, margin=8, PhiInit=spam.deformation.computePhi({'t': [5,0,0]}), PhiInitBinRatio=2, updateGradient=True)
t5 = spam.deformation.decomposePhi(reg4['Phi']) t5 = spam.deformation.decomposePhi(reg4['Phi'])
for c in ['t', 'r']: for c in ['t', 'r']:
for i in range(3): for i in range(3):
self.assertAlmostEqual(t5[c][i] - t[c][i], 0, places=1) self.assertAlmostEqual(t5[c][i] - t[c][i], 0, places=1)
# Test for margin = None # Test for margin = None
reg6 = spam.DIC.registerMultiscale(im, imDef, 2, margin=None, PhiInit=spam.deformation.computePhi({'t': [5,0,0]}), PhiInitBinRatio=2, updateGradient=True) reg6 = spam.DIC.registerMultiscale(im, imDef, 2, margin=None, PhiInit=spam.deformation.computePhi({'t': [5,0,0]}), PhiInitBinRatio=2, updateGradient=True)
t6 = spam.deformation.decomposePhi(reg6['Phi']) t6 = spam.deformation.decomposePhi(reg6['Phi'])
...@@ -364,7 +377,7 @@ class TestFunctionDVC(unittest.TestCase): ...@@ -364,7 +377,7 @@ class TestFunctionDVC(unittest.TestCase):
# self.assertTrue(ps2["cc"]>0.9) # self.assertTrue(ps2["cc"]>0.9)
# for i in range(3): # for i in range(3):
# self.assertAlmostEqual(tVector[i], ps2["transformation"]["t"][i], places=4) # self.assertAlmostEqual(tVector[i], ps2["transformation"]["t"][i], places=4)
# Test for imagette2.shape >= imagette1.shape # Test for imagette2.shape >= imagette1.shape
ps2 = spam.DIC.pixelSearch(imDef, im[subVolSlice1].copy(), ps2 = spam.DIC.pixelSearch(imDef, im[subVolSlice1].copy(),
searchCentre=p, searchCentre=p,
......
...@@ -461,11 +461,11 @@ def ddicParser(parser): ...@@ -461,11 +461,11 @@ def ddicParser(parser):
dest='PSR', dest='PSR',
help='Z- Z+ Y- Y+ X- X+ ranges (in pixels) for the pxiel search. Requires pixel search to be activated. Default = +-3px') 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('-nogp', parser.add_argument('-nolc',
'--no-grid-point', '--no-label-correlation',
action="store_false", action="store_false",
dest='GRID_POINT', dest='LABEL_CORRELATE',
help='Disactivate grid-point registration?') help='Disactivate label registration?')
parser.add_argument('-ld', parser.add_argument('-ld',
'--label-dilate', '--label-dilate',
...@@ -534,7 +534,6 @@ def ddicParser(parser): ...@@ -534,7 +534,6 @@ def ddicParser(parser):
dest='PHIFILE_BIN_RATIO', dest='PHIFILE_BIN_RATIO',
help="Ratio of binning level between loaded F file and current calculation. If the input F file has been obtained on a 500x500x500 image and now the calculation is on 1000x1000x1000, this should be 2. Default = 1") help="Ratio of binning level between loaded F file and current calculation. If the input F file has been obtained on a 500x500x500 image and now the calculation is on 1000x1000x1000, this should be 2. Default = 1")
parser.add_argument('-pfd', parser.add_argument('-pfd',
'--phiFile-direct', '--phiFile-direct',
action="store_true", action="store_true",
...@@ -561,34 +560,46 @@ def ddicParser(parser): ...@@ -561,34 +560,46 @@ def ddicParser(parser):
dest='MASK', dest='MASK',
help='Don\'t mask each label\'s image') help='Don\'t mask each label\'s image')
parser.add_argument('-gpm', parser.add_argument('-lcm',
'--grid-point-margin', '--label-correlate-margin',
type=numpy.uint, type=numpy.uint,
default=2, default=5,
dest='GRID_POINT_MARGIN', dest='LABEL_CORRELATE_MARGIN',
help="Margin in pixels for grid-point registration to be added to label dilate. Default = 2px") help="Margin in pixels for grid-point registration to be added to label dilate. Default = 2px")
parser.add_argument('-gpi', parser.add_argument('-lci',
'--grid-point-max-iterations', '--label-correlate-max-iterations',
type=numpy.uint, type=numpy.uint,
default=50, default=50,
dest='GRID_POINT_MAX_ITERATIONS', dest='LABEL_CORRELATE_MAX_ITERATIONS',
help="Maximum iterations for grid-point registration. Default = 50") help="Maximum iterations for grid-point registration. Default = 50")
parser.add_argument('-gpp', parser.add_argument('-lcp',
'--grid-point-min-phi-change', '--label-correlate-min-phi-change',
type=numpy.float, type=numpy.float,
default=0.001, default=0.001,
dest='GRID_POINT_MIN_PHI_CHANGE', dest='LABEL_CORRELATE_MIN_PHI_CHANGE',
help="Minimum change in Phi to consider grid-point registration as converged. Default = 0.001") help="Minimum change in Phi to consider grid-point registration as converged. Default = 0.001")
parser.add_argument('-gpo', parser.add_argument('-lco',
'--grid-point-interpolation-order', '--label-correlate-interpolation-order',
type=numpy.uint, type=numpy.uint,
default=1, default=1,
dest='GRID_POINT_INTERPOLATION_ORDER', dest='LABEL_CORRELATE_INTERPOLATION_ORDER',
help="Interpolation order for grid-point registration. Default = 3") help="Interpolation order for grid-point registration. Default = 3")
parser.add_argument('-lnr',
'--label-correlate-non-rigid',
action="store_false",
dest='LABEL_CORRELATE_RIGID',
help='Activate non-rigid registration for each label')
parser.add_argument('-lcug',
'--label-correlate-update-gradient',
action="store_true",
dest='LABEL_CORRELATE_UPDATE_GRADIENT',
help='Update gradient in label registration? More computation time but more robust and possibly fewer iterations. Default = False')
parser.add_argument('-od', parser.add_argument('-od',
'--out-dir', '--out-dir',
type=str, type=str,
......
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