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

merging ddic-debug -- there is now a debug mode for ddic, its correlation is...

merging ddic-debug -- there is now a debug mode for ddic, its correlation is by default rigid, and there is a new spam-ereg-discrete that helps align label manually
parents 384f67ff 5764bd18
Pipeline #49337 passed with stages
in 23 minutes and 32 seconds
-r requirements.txt
coverage
sphinx>=2.0
sphinx-gallery>=0.1.13
......
numpy==1.18.5
scipy==1.4.1
opencv-python==3.4.10.35
SimpleITK==1.2.4
mpi4py==3.0.3
......
This diff is collapsed.
#!/usr/bin/env python
"""
This python facilitates eye-alignment with a graphical QT interface
for Discrete Digital Image Correlation using SPAM functions
Copyright (C) 2020 SPAM Contributors
This program is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the Free
Software Foundation, either version 3 of the License, or (at your option)
any later version.
This program is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for
more details.
You should have received a copy of the GNU General Public License along with
this program. If not, see <http://www.gnu.org/licenses/>.
"""
import os
import sys
import subprocess
import spam.helpers
import spam.label
import numpy
import tifffile
import spam.visual.visualClass as visual
import argparse
from PyQt5.QtWidgets import QApplication, QWidget, QFileDialog, QGridLayout, QPushButton
# Define argument parser object
parser = argparse.ArgumentParser(description="spam-ereg-discrete "+spam.helpers.optionsParser.GLPv3descriptionHeader +\
"This script facilitates eye-alignment for Discrete Digital Image Correlation two 3D greyscale images"+\
" (reference and deformed configurations) and requires the input of a labelled image for the reference configuration",
formatter_class=argparse.RawTextHelpFormatter)
# Parse arguments with external helper function
args = spam.helpers.optionsParser.eregDiscreteParser(parser)
print("spam-ereg-discrete: Current Settings:")
argsDict = vars(args)
for key in sorted(argsDict):
print("\t{}: {}".format(key, argsDict[key]))
outFile = args.OUT_DIR+"/"+args.PREFIX+".tsv"
REFlab = tifffile.imread(args.lab1.name)
REFlabBB = spam.label.boundingBoxes(REFlab)
REFlabCOM = spam.label.centresOfMass(REFlab)
REFgrey = tifffile.imread(args.im1.name)
DEFgrey = tifffile.imread(args.im2.name)
# Some variable for nice filenames
REFstr = os.path.basename(args.im1.name)
DEFstr = os.path.basename(args.im2.name)
if args.PHIFILE is not None:
DDIC = spam.helpers.readCorrelationTSV(args.PHIFILE.name)
DDIC['error'] = numpy.genfromtxt(args.PHIFILE.name, delimiter='\t', names=True)['error']
DDIC['LabelDilate'] = numpy.genfromtxt(args.PHIFILE.name, delimiter='\t', names=True)['LabelDilate']
DDIC['PSCC'] = numpy.genfromtxt(args.PHIFILE.name, delimiter='\t', names=True)['PSCC']
else:
labMax = REFlab.max()
DDIC = {}
DDIC['fieldCoords'] = REFlabCOM
DDIC['PhiField'] = numpy.zeros([labMax, 4, 4])
for label in range(labMax):
DDIC['PhiField'][label] = numpy.eye(4)
DDIC['returnStatus'] = numpy.zeros(labMax, dtype=int)
DDIC['deltaPhiNorm'] = numpy.zeros(labMax, dtype=int)
DDIC['iterations'] = numpy.zeros(labMax, dtype=int)
DDIC['error'] = numpy.zeros(labMax, dtype=int)
DDIC['LabelDilate'] = numpy.zeros(labMax, dtype=int)
DDIC['PSCC'] = numpy.zeros(labMax, dtype=int)
class MainWindow(QWidget):
def __init__(self):
QWidget.__init__(self)
self.Phi = numpy.eye(4)
self.mainWindowGrid = QGridLayout(self)
self.nonConvergedGrains = numpy.where(DDIC['returnStatus'] < args.RETURN_STAT_THRESHOLD)[0][1:]
self.N = 0 # Number of the current nonConvergedGrain that's being studied
print("Going to work on these labels:\n", self.nonConvergedGrains)
if len(self.nonConvergedGrains) > 0:
self.alignOneLabel()
self.labAndPhi = []
else:
print("No labels to work on")
exit()
def alignOneLabel(self):
nonConvergedGrain = self.nonConvergedGrains[self.N]
print("\tGrain {}".format(nonConvergedGrain))
Phi = DDIC['PhiField'][nonConvergedGrain]
displacement = Phi[0:3,-1]
displacementInt = displacement.astype(int)
self.diplacementInt = displacementInt
# Remove the int part of displacement
Phi[0:3,-1] -= displacementInt
print("\t\tSubtracted this displacement:", displacementInt)
#REFsubvol = REFgrey[REFlabBB[nonConvergedGrain,0]-args.margin : REFlabBB[nonConvergedGrain,1]+1+args.margin,
#REFlabBB[nonConvergedGrain,2]-args.margin : REFlabBB[nonConvergedGrain,3]+1+args.margin,
#REFlabBB[nonConvergedGrain,4]-args.margin : REFlabBB[nonConvergedGrain,5]+1+args.margin]
REFgl = spam.label.getLabel(REFlab, nonConvergedGrain,
boundingBoxes=REFlabBB, centresOfMass=REFlabCOM,
labelDilate=args.LABEL_DILATE, margin=args.margin,
maskOtherLabels=args.MASK)
REFsubvol = REFgrey[REFgl['slice']]
if args.MASK:
# If mask asked, also flatten greylevels
REFsubvol[REFgl['subvol'] == 0] = 0
DEFsubvol = DEFgrey[REFlabBB[nonConvergedGrain,0]+displacementInt[0]-args.margin : REFlabBB[nonConvergedGrain,1]+1+displacementInt[0]+args.margin,
REFlabBB[nonConvergedGrain,2]+displacementInt[1]-args.margin : REFlabBB[nonConvergedGrain,3]+1+displacementInt[1]+args.margin,
REFlabBB[nonConvergedGrain,4]+displacementInt[2]-args.margin : REFlabBB[nonConvergedGrain,5]+1+displacementInt[2]+args.margin]
self.eregWidget = visual.ereg( [REFsubvol, DEFsubvol],
Phi, 1,
["{} - label {}".format(REFstr, nonConvergedGrain), "{} - label {}".format(DEFstr, nonConvergedGrain)], 0)
self.mainWindowGrid.addWidget(self.eregWidget, 1, 1)
self.nextLabelButton = QPushButton("Accept and move on to next grain", self)
self.nextLabelButton.clicked.connect(self.nextLabel)
self.mainWindowGrid.addWidget(self.nextLabelButton, 2, 1)
def nextLabel(self):
#print("Phi:", self.eregWidget.output())
self.eregWidget.close()
# Get Phi output from graphical
PhiTmp = self.eregWidget.output()
# Add back in int displacement
PhiTmp[0:3, -1] += self.diplacementInt
# nonConvergedGrain label number, eye-Phi
self.labAndPhi.append([self.nonConvergedGrains[self.N], PhiTmp])
# Move onto next grain, otherwise write and quit
self.N += 1
if self.N < len(self.nonConvergedGrains):
self.alignOneLabel()
else:
self.nextLabelButton.close()
self.eregWidget.close()
#print(self.labAndPhi)
print("Updating output...")
for nonConvergedGrain, Phi in self.labAndPhi:
DDIC['PhiField'][nonConvergedGrain] = Phi
print("Writing output to {}...".format(outFile), end='')
outMatrix = numpy.array([numpy.array(range(DDIC['numberOfLabels'])),
DDIC['fieldCoords'][:, 0], DDIC['fieldCoords'][:, 1], DDIC['fieldCoords'][:, 2],
DDIC['PhiField'][:, 0, 3], DDIC['PhiField'][:, 1, 3], DDIC['PhiField'][:, 2, 3],
DDIC['PhiField'][:, 0, 0], DDIC['PhiField'][:, 0, 1], DDIC['PhiField'][:, 0, 2],
DDIC['PhiField'][:, 1, 0], DDIC['PhiField'][:, 1, 1], DDIC['PhiField'][:, 1, 2],
DDIC['PhiField'][:, 2, 0], DDIC['PhiField'][:, 2, 1], DDIC['PhiField'][:, 2, 2],
DDIC['error'], DDIC['iterations'],
DDIC['returnStatus'], DDIC['deltaPhiNorm'],
DDIC['LabelDilate'], DDIC['PSCC']]).T
numpy.savetxt(outFile,
outMatrix,
fmt='%.7f',
delimiter='\t',
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" +
"error\titerations\treturnStatus\tdeltaPhiNorm\tLabelDilate\tPSCC")
print("...done")
self.close()
#self.mainWindowGrid.close()
app = QApplication(["Label Registration"])
window = MainWindow()
window.show()
app.exec_()
......@@ -131,6 +131,7 @@ scripts = ['scripts/spam-mmr',
"scripts/spam-ddic",
"scripts/spam-ldic",
"scripts/spam-ereg",
"scripts/spam-ereg-discrete",
"scripts/spam-deformImageFromField",
"scripts/spam-mmr-graphical",
"scripts/spam-moveGrains",
......
......@@ -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']["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
rot = 120
Phi = spam.deformation.computePhi({'r': [rot, 0, 0]})
......@@ -292,18 +303,20 @@ class TestFunctionDVC(unittest.TestCase):
for i in range(3):
self.assertAlmostEqual(t1[c][i] - t[c][i], 0, places=1)
reg2 = spam.DIC.registerMultiscale(im, imDef, 1.1, PhiInit=numpy.eye(4), margin=[8, 8, 8])
reg2 = spam.DIC.registerMultiscale(im, imDef, 1.1, PhiInit=numpy.eye(4), margin=8)
t2 = spam.deformation.decomposePhi(reg2['Phi'])
for c in ['t', 'r']:
for i in range(3):
self.assertAlmostEqual(t2[c][i] - t[c][i], 0, places=1)
# Check translation example with translation guess
t = {'t': [4.0, 0, 0],
'r': [0.0, 0, 0]}
Phi = spam.deformation.computePhi(t)
imDef = spam.DIC.applyPhi(im, Phi=Phi)
reg3 = spam.DIC.registerMultiscale(im, imDef, 2, deltaPhiMin=0.00001, PhiInit=Phi)
reg3 = spam.DIC.registerMultiscale(im, imDef, 2, deltaPhiMin=0.00001, margin=4, PhiInit=Phi)
t3 = spam.deformation.decomposePhi(reg3['Phi'])
for c in ['t', 'r']:
for i in range(3):
......@@ -315,7 +328,7 @@ class TestFunctionDVC(unittest.TestCase):
for i in range(3):
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)
t5 = spam.deformation.decomposePhi(reg4['Phi'])
for c in ['t', 'r']:
......
......@@ -74,8 +74,8 @@ class testAll(unittest.TestCase):
self.assertEqual(numpy.add(stack, -snow).sum(), 0)
# test all cases
spam.helpers.stackToArray("tmp_", stack=range(100), createMask=True, verbose=True)
spam.helpers.stackToArray("tmp_", stack=range(100), createMask=True, erosion=True)
spam.helpers.stackToArray("tmp_", stack=range(100), verbose=True)
spam.helpers.stackToArray("tmp_", stack=range(100), erosion=True)
def test_shift(self):
# create slices from snow.tif, restack
......
......@@ -213,10 +213,16 @@ class TestFunctionLabel(unittest.TestCase):
# 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=1, labelDilate=-2)
gl = spam.label.getLabel(threeCubedLabelVol, 1, boundingBoxes=BB, centresOfMass=COM, margin=2, labelDilate=2)
volRef = gl['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)
def test_labelsOnEdges(self):
labelled = numpy.zeros((5, 5, 5), dtype="<u2")
......
......@@ -356,10 +356,10 @@ class testAll(unittest.TestCase):
# put 0 in the middle
centres -= numpy.mean(centres, axis=0)
rMax = numpy.amax(radii)
# pad box size
boxSizeDEM = boxSizeDEM + 5 * rMax
# add half box size to centres
centres += numpy.array(boxSizeDEM)/2.0
boxSizePx = (boxSizeDEM / pixelSize).astype(int)
......@@ -379,7 +379,7 @@ class testAll(unittest.TestCase):
tifffile.imsave(testFolder + "Lab0.tif", labIm0.astype(spam.label.labelType))
#test of rigid translation and rotation
#Create Phi and Apply (25 px displacement on Y-axis, and 5 degree rotation along Z axis)
#Create Phi and Apply (2 px displacement on X-axis, and 5 degree rotation along Z axis)
translationStep1 = [0, 0, 2]
rotationStep1 = [5, 0, 0]
transformation = {'t': translationStep1, 'r': rotationStep1}
......@@ -401,7 +401,6 @@ class testAll(unittest.TestCase):
tifffile.imsave(testFolder + "Step1.tif", boxDeformed.astype('<f4'))
#test of rigid translation and rotation
#Create Phi and Apply (25 px displacement on Y-axis, and 5 degree rotation along Z axis)
translationStep2 = [0, 0, 18]
rotationStep2 = [0, 0, 0]
transformation = {'t': translationStep2, 'r': rotationStep2}
......@@ -471,7 +470,8 @@ class testAll(unittest.TestCase):
testFolder + "Lab0.tif",
testFolder + "Step1.tif",
"-od", testFolder + "",
"-pf", testFolder + "Step0-Step1-discreteDVC.tsv"])
"-pf", testFolder + "Step0-Step1-discreteDVC.tsv",
"-pfd"])
self.assertEqual(exitCode, 0)
#Check number of converged grains
......@@ -579,7 +579,8 @@ class testAll(unittest.TestCase):
testFolder + "Lab0.tif",
testFolder + "Step2.tif",
"-od", testFolder + "",
"-pf", testFolder + "Step0-Step2-discreteDVC.tsv"])
"-pf", testFolder + "Step0-Step2-discreteDVC.tsv",
"-pfd"])
self.assertEqual(exitCode, 0)
# Check TSV for local results
......@@ -606,7 +607,7 @@ class testAll(unittest.TestCase):
self.assertAlmostEqual(rotationStep2[2], averageRot[2], places=0)
#######################################################
### 4. Check -skp
### 4. Check -skp
#######################################################
# Manually modify the status of the last label
TSV['PhiField'][-1][:] = 0
......@@ -622,7 +623,8 @@ class testAll(unittest.TestCase):
TSV2 = spam.helpers.readCorrelationTSV(testFolder + "Step0-Step2-discreteDVC.tsv")
# Compare "deltaPhiNorm" for all but the last one (they should be the same since they were skipped)
for i in range(1, TSV2['PhiField'].shape[0]-1):
self.assertEqual(TSV['deltaPhiNorm'][i], TSV2['deltaPhiNorm'][i])
for u in range(3):
self.assertAlmostEqual(TSV['PhiField'][i][u,-1], TSV2['PhiField'][i][u,-1], places=1)
def test_discreteStrain(self):
# make sure it runs the help without error
......
......@@ -32,7 +32,7 @@ import progressbar
#numpy.set_printoptions(precision=3, suppress=True)
def register(im1, im2, im1mask=None, PhiInit=None, PhiInitBinRatio=1.0, margin=None, maxIterations=25, deltaPhiMin=0.001, updateGradient=False, interpolationOrder=1, interpolator='python', verbose=False, imShowProgress=None, imShowProgressNewFig=False):
def register(im1, im2, im1mask=None, PhiInit=None, PhiRigid=False, PhiInitBinRatio=1.0, margin=None, maxIterations=25, deltaPhiMin=0.001, updateGradient=False, interpolationOrder=1, interpolator='python', verbose=False, imShowProgress=False, imShowProgressNewFig=False):
"""
Perform subpixel image correlation between im1 and im2.
......@@ -65,6 +65,10 @@ def register(im1, im2, im1mask=None, PhiInit=None, PhiInitBinRatio=1.0, margin=N
Initial deformation to apply to im1.
Default = numpy.eye(4), `i.e.`, no transformation
PhiRigid : bool, optional
Run a rigid correlation? Only the rigid part of your PhiInit will be kept.
Default = False
PhiInitBinRatio : float, optional
Change translations in PhiInit, if it's been calculated on a differently-binned image. Default = 1
......@@ -92,11 +96,12 @@ def register(im1, im2, im1mask=None, PhiInit=None, PhiInitBinRatio=1.0, margin=N
Default = 'python'. 'C' is also an option
verbose : bool, optional
Get to know what the function is really thinking, recommended for debugging only. Default = False
Get to know what the function is really thinking, recommended for debugging only.
Default = False
imShowProgress : String, optional (default = None)
imShowProgress : Bool, optional
Pop up a window showing a ``imShowProgress`` slice of the image differences (im1-im2) as im1 is progressively deformed.
Accepted options are "Z", "Y" and "X" -- the slicing direction.
Default = False
imShowProgressNewFig : bool, optional (defaul = False)
Make a new plt.figure for each iteration, useful for examples gallery
......@@ -138,10 +143,6 @@ def register(im1, im2, im1mask=None, PhiInit=None, PhiInitBinRatio=1.0, margin=N
Since we're using classical DIC (without multimodal registration) there is no need to update the gradient image,
so this is calculated once at the beginning.
"""
# History:
# 2018-03-20 EA: Trying to make margin a list of 3 values as z, y, x margins in order to try
# clean-ish 2D compatilibility
# Explicitly set input images to floats
im1 = im1.astype('<f4')
im2 = im2.astype('<f4')
......@@ -227,18 +228,33 @@ def register(im1, im2, im1mask=None, PhiInit=None, PhiInitBinRatio=1.0, margin=N
# print( "\tcorrelate.register(): realMargin is:", realMargin)
# If live plot is asked for, initialise canvas
if imShowProgress is not None:
if imShowProgress:
import matplotlib.pyplot as plt
# Plot ranges for signed residual
vmin = -im1crop.max() / 2
vmax = im1crop.max() / 2
vmin = -im1crop.max()
vmax = im1crop.max()
if not imShowProgressNewFig:
if imShowProgress == "Z" or imShowProgress == "z":
plt.axis([im1crop.shape[2], 0, im1crop.shape[1], 0])
if imShowProgress == "Y" or imShowProgress == "y":
plt.axis([im1crop.shape[2], 0, im1crop.shape[0], 0])
if imShowProgress == "X" or imShowProgress == "x":
plt.axis([im1crop.shape[1], 0, im1crop.shape[0], 0])
#if imShowProgress == "Z" or imShowProgress == "z":
plt.subplot(3,3,1)
plt.axis([im1crop.shape[2], 0, im1crop.shape[1], 0])
plt.subplot(3,3,2)
plt.axis([im1crop.shape[2], 0, im1crop.shape[1], 0])
plt.subplot(3,3,3)
plt.axis([im1crop.shape[2], 0, im1crop.shape[1], 0])
#if imShowProgress == "Y" or imShowProgress == "y":
plt.subplot(3,3,4)
plt.axis([im1crop.shape[2], 0, im1crop.shape[0], 0])
plt.subplot(3,3,5)
plt.axis([im1crop.shape[2], 0, im1crop.shape[0], 0])
plt.subplot(3,3,6)
plt.axis([im1crop.shape[2], 0, im1crop.shape[0], 0])
#if imShowProgress == "X" or imShowProgress == "x":
plt.subplot(3,3,7)
plt.axis([im1crop.shape[1], 0, im1crop.shape[0], 0])
plt.subplot(3,3,8)
plt.axis([im1crop.shape[1], 0, im1crop.shape[0], 0])
plt.subplot(3,3,9)
plt.axis([im1crop.shape[1], 0, im1crop.shape[0], 0])
plt.ion()
# Numerical value for normalising the error
......@@ -269,6 +285,10 @@ def register(im1, im2, im1mask=None, PhiInit=None, PhiInitBinRatio=1.0, margin=N
# Copying into different variable
# Apply binning on displacement
Phi = PhiInit.copy()
# If we're in rigid mode, keep only translations and rotations for this guess
# If you don't do this it goes mad (i.e., rigid updates to non-rigid guess don't seem to work)
if PhiRigid: Phi = spam.deformation.computeRigidPhi(Phi.copy())
Phi[0:3, -1] *= PhiInitBinRatio
# invert PhiInit to apply it to im2
......@@ -301,21 +321,6 @@ def register(im1, im2, im1mask=None, PhiInit=None, PhiInitBinRatio=1.0, margin=N
#if verbose: print("done ")
return imGradZ, imGradY, imGradX
# Calculate gradient of the non-moving im1
# 2017-11-06 EA: This is going in a try, since it's possible to have an error like this:
# "Shape of array too small to calculate a numerical gradient, "
# ValueError: Shape of array too small to calculate a numerical gradient, at least (edge_order + 1) elements are required.
#try:
#if updateGradient:
#imGradZ, imGradY, imGradX = computeGradient(im2def, twoD)
#else:
#imGradZ, imGradY, imGradX = computeGradient(im1, twoD)
## If gradient calculation failed, set singular to true, means early exit
#except ValueError:
## Override max iteration and also set singular
#maxIterations = 0
#singular = True
# Apply stationary im1 mask
if im1mask is not None:
im1crop[im1mask[crop1] == False] = numpy.nan
......@@ -337,6 +342,7 @@ def register(im1, im2, im1mask=None, PhiInit=None, PhiInitBinRatio=1.0, margin=N
# pbar = progressbar.ProgressBar(widgets=widgets, maxval=numberOfNodes)
pbar.start()
# --- Start Iterations ---
while iterations < maxIterations and deltaPhiNorm > deltaPhiMin:
errorPrev = error
......@@ -394,23 +400,39 @@ def register(im1, im2, im1mask=None, PhiInit=None, PhiInitBinRatio=1.0, margin=N
singular = True
break
deltaPhiNorm = numpy.linalg.norm(deltaPhi)
# If we're doing a rigid registration...
if PhiRigid:
# Add padding zeros
deltaPhi = numpy.hstack([deltaPhi, numpy.zeros(4)]).reshape((4, 4))
# Add padding zeros
deltaPhi = numpy.hstack([deltaPhi, numpy.zeros(4)]).reshape((4, 4))
deltaPhiPlusI = (numpy.eye(4) + deltaPhi)
# Keep only rigid part of deltaPhi
deltaPhiPlusIrigid = spam.deformation.computeRigidPhi(deltaPhiPlusI.copy())
# Apply Delta Phi correction to Phi In Roux X-N paper equation number 11
#Phi = numpy.dot((numpy.eye(4) + deltaPhi), Phi)
Phi = numpy.dot(Phi, (numpy.eye(4) + deltaPhi))
# Subtract I from the rigid dPhi+1, and compute norm only on first 3 rows
# ...basically recompute deltaPhiNorm only on rigid part
deltaPhiNorm = numpy.linalg.norm((deltaPhiPlusIrigid-numpy.eye(4))[0:3].ravel())
# Apply Delta Phi correction to Phi In Roux X-N paper equation number 11
Phi = numpy.dot(Phi, deltaPhiPlusIrigid)
else:
# The general, non-rigid case
deltaPhiNorm = numpy.linalg.norm(deltaPhi)
# Add padding zeros
deltaPhi = numpy.hstack([deltaPhi, numpy.zeros(4)]).reshape((4, 4))
# Update Phi
Phi = numpy.dot(Phi, (numpy.eye(4) + deltaPhi))
# Solve for delta Phi
try: PhiInv = numpy.linalg.inv(Phi.copy())
except numpy.linalg.linalg.LinAlgError:
singular = True
break
# reset im1def as emtpy matrix for deformed image
if interpolator == 'C': im2def = spam.DIC.applyPhi(im2, Phi=PhiInv, interpolationOrder=interpolationOrder)
if interpolator == 'C': im2def = spam.DIC.applyPhi(im2, Phi=PhiInv, interpolationOrder=interpolationOrder)
elif interpolator == 'python': im2def = spam.DIC.applyPhiPython(im2, Phi=PhiInv, interpolationOrder=interpolationOrder)
# Error calculation
......@@ -432,10 +454,10 @@ def register(im1, im2, im1mask=None, PhiInit=None, PhiInitBinRatio=1.0, margin=N
# Second divergence criterion on displacement (Issue #62)
# If any displcement is bigger than 5* the margin...
if (numpy.abs(spam.deformation.decomposePhi(Phi.copy())['t']) > 5 * realMargin).any():
if verbose: print("\t -> diverging on displacement condition")
returnStatus = -3
break
#if (numpy.abs(spam.deformation.decomposePhi(Phi.copy())['t']) > 5 * realMargin).any():
#if verbose: print("\t -> diverging on displacement condition")
#returnStatus = -3
#break
# 2018-10-02 - EA: Add divergence condition on U
trans = spam.deformation.decomposePhi(Phi.copy())
......@@ -449,32 +471,46 @@ def register(im1, im2, im1mask=None, PhiInit=None, PhiInitBinRatio=1.0, margin=N
returnStatus = -3
break
if imShowProgress is not None:
if imShowProgress == "Z" or imShowProgress == "z":
if imShowProgressNewFig: plt.figure()
else: plt.clf()
plt.imshow(numpy.subtract(im1crop, im2def[crop2])[im1crop.shape[0] // 2, :, :], cmap='coolwarm', vmin=vmin, vmax=vmax)
if imShowProgress == "Y" or imShowProgress == "y":
if imShowProgressNewFig: plt.figure()
else: plt.clf()
plt.imshow(numpy.subtract(im1crop, im2def[crop2])[:, im1crop.shape[1] // 2, :], cmap='coolwarm', vmin=vmin, vmax=vmax)
if imShowProgress == "X" or imShowProgress == "x":
if imShowProgressNewFig: plt.figure()
else: plt.clf()
plt.imshow(numpy.subtract(im1crop, im2def[crop2])[:, :, im1crop.shape[2] // 2], cmap='coolwarm', vmin=vmin, vmax=vmax)
if imShowProgress:
#if imShowProgress == "Z" or imShowProgress == "z":
if imShowProgressNewFig: plt.figure()
else: plt.clf()
plt.subplot(3,3,1)
plt.imshow(im1crop[im1crop.shape[0]//2, :, :], cmap='Greys_r', vmin=0, vmax=vmax)
plt.subplot(3,3,2)
plt.imshow(im2def[crop2][im1crop.shape[0]//2, :, :], cmap='Greys_r', vmin=0, vmax=vmax)
plt.subplot(3,3,3)
plt.imshow(numpy.subtract(im1crop, im2def[crop2])[im1crop.shape[0]//2, :, :], cmap='coolwarm', vmin=vmin, vmax=vmax)
#if imShowProgress == "Y" or imShowProgress == "y":
#if imShowProgressNewFig: plt.figure()
#else: plt.clf()
plt.subplot(3,3,4)
plt.imshow(im1crop[:, im1crop.shape[1]//2, :], cmap='Greys_r', vmin=0, vmax=vmax