Coverage for /usr/local/lib/python3.8/site-packages/spam/DIC/registration.py: 83%

1""" 

2Library of SPAM image correlation functions. 

3Copyright (C) 2020 SPAM Contributors 

4 

5This program is free software: you can redistribute it and/or modify it 

6under the terms of the GNU General Public License as published by the Free 

7Software Foundation, either version 3 of the License, or (at your option) 

8any later version. 

9 

10This program is distributed in the hope that it will be useful, but WITHOUT 

11ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or 

12FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for 

13more details. 

14 

15You should have received a copy of the GNU General Public License along with 

16this program. If not, see <http://www.gnu.org/licenses/>. 

17""" 

18 

19import numpy 

20import progressbar 

21import spam.deformation 

22import spam.DIC 

23import spam.label # for im1mask 

24 

25# 2017-05-29 ER and EA 

26# This is spam's C++ DIC toolkit, but since we're in the tools/ directory we can import it directly 

27from spam.DIC.DICToolkit import computeDICjacobian, computeDICoperators 

28 

29# numpy.set_printoptions(precision=3, suppress=True) 

30 

31 

32def _errorCalc(im1, im2): 

33 return numpy.nansum(numpy.square(numpy.subtract(im1, im2))) / numpy.nansum(im1) 

34 

35 

36def register( 

37 im1, 

38 im2, 

39 im1mask=None, 

40 PhiInit=None, 

41 PhiRigid=False, 

42 PhiInitBinRatio=1.0, 

43 margin=None, 

44 maxIterations=25, 

45 deltaPhiMin=0.001, 

46 updateGradient=False, 

47 interpolationOrder=1, 

48 interpolator="python", 

49 verbose=False, 

50 imShowProgress=False, 

51 imShowProgressNewFig=False, 

52): 

53 r""" 

54 Perform subpixel image correlation between im1 and im2. 

55 

56 The result of register(im1, im2) will give a deformation function :math:`\Phi` which maps im1 into im2. 

57 The Phi function used here allows the measurement of sub-pixel displacements, rotation, and linear straining of the whole image. 

58 However, this function will numerically deform im2 until it best matches im1. 

59 

60 :math:`im1(x) = im2(\Phi.x)` 

61 

62 If im1 and im2 follow each other in time, then the resulting Phi is im1 -> im2 which makes sense in most cases. 

63 "Discrete correlation" can be performed by masking im1. 

64 

65 im1 and im2 do not necessarily have to be the same size (`i.e.`, im2 can be bigger) -- this is good since there 

66 is a zone to accommodate movement. In the case of a bigger im2, im1 and im2 are assumed to be centred with respect to each other. 

67 

68 Parameters 

69 ---------- 

70 im1 : 3D numpy array 

71 The greyscale image that will not move -- must not contain NaNs 

72 

73 im2 : 3D numpy array 

74 The greyscale image that will be deformed -- must not contain NaNs 

75 

76 im1mask : 3D boolean numpy array, optional 

77 A mask for the zone to correlate in im1 with `False` in the zone to not correlate. 

78 Default = None, `i.e.`, correlate all of im1 minus the margin. 

79 If this is defined, the Phi returned is in the centre of mass of the mask 

80 

81 PhiInit : 4x4 numpy array, optional 

82 Initial deformation to apply to im1. 

83 Default = numpy.eye(4), `i.e.`, no transformation 

84 

85 PhiRigid : bool, optional 

86 Run a rigid correlation? Only the rigid part of your PhiInit will be kept. 

87 Default = False 

88 

89 PhiInitBinRatio : float, optional 

90 Change translations in PhiInit, if it's been calculated on a differently-binned image. Default = 1 

91 

92 margin : int, optional 

93 Margin, in pixels, to take in im1. 

94 Can also be a N-component list of ints, representing the margin in ND. 

95 If im2 has the same size as im1 this is strictly necessary to allow space for interpolation and movement 

96 Default = None (`i.e.`, 10% of max dimension of im1) 

97 

98 maxIterations : int, optional 

99 Maximum number of quasi-Newton iterations to perform before stopping. Default = 25 

100 

101 deltaPhiMin : float, optional 

102 Smallest change in the norm of Phi (the transformation operator) before stopping. Default = 0.001 

103 

104 updateGradient : bool, optional 

105 Should the gradient of the image be computed (and updated) on the deforming im2? 

106 Default = False (it is computed once on im1) 

107 

108 interpolationOrder : int, optional 

109 Order of the greylevel interpolation for applying Phi to im1 when correlating. Recommended value is 3, but you can get away with 1 for faster calculations. Default = 3 

110 

111 interpolator : string, optional 

112 Which interpolation function to use from `spam`. 

113 Default = 'python'. 'C' is also an option 

114 

115 verbose : bool, optional 

116 Get to know what the function is really thinking, recommended for debugging only. 

117 Default = False 

118 

119 imShowProgress : bool, optional 

120 Pop up a window showing a ``imShowProgress`` slice of the image differences (im1-im2) as im1 is progressively deformed. 

121 Default = False 

122 

123 imShowProgressNewFig : bool, optional (defaul = False) 

124 Make a new plt.figure for each iteration, useful for examples gallery 

125 

126 Returns 

127 ------- 

128 Dictionary : 

129 

130 'Phi' : 4x4 float array 

131 Deformation function defined at the centre of the image 

132 

133 'returnStatus' : signed int 

134 Return status from the correlation: 

135 

136 2 : Achieved desired precision in the norm of delta Phi 

137 

138 1 : Hit maximum number of iterations while iterating 

139 

140 -1 : Error is more than 80% of previous error, we're probably diverging 

141 

142 -2 : Singular matrix M cannot be inverted 

143 

144 -3 : Displacement > 5*margin 

145 

146 'error' : float 

147 Error float describing mismatch between images, it's the sum of the squared difference divided by the sum of im1 

148 

149 'iterations' : int 

150 Number of iterations 

151 

152 'deltaPhiNorm' : float 

153 Norm of deltaPhi 

154 

155 Note 

156 ---- 

157 This correlation was written in the style of S. Roux (especially "An extension of Digital Image Correlation for intermodality image registration") 

158 especially equations 12 and 13. 

159 """ 

160 # Explicitly set input images to floats 

161 im1 = im1.astype("<f4") 

162 im2 = im2.astype("<f4") 

163 

164 # initialise exit clause for singular "M" matrices 

165 singular = False 

166 

167 # 2022-06-03 GP: Setting deltaPhiMin to only positive values 

168 deltaPhiMin = numpy.abs(deltaPhiMin) 

169 

170 # Detect unpadded 2D image first: 

171 if im1.ndim == 2: 

172 # pad them 

173 im1 = im1[numpy.newaxis, ...] 

174 im2 = im2[numpy.newaxis, ...] 

175 if im1mask is not None: 

176 im1mask = im1mask[numpy.newaxis, ...] 

177 

178 # Detect 2D images 

179 if im1.shape[0] == 1: 

180 twoD = True 

181 

182 # Override interpolator for python in 2D 

183 interpolator = "python" 

184 

185 # Define masks for M and A in 2D since we'll ignore the Z components 

186 # Components of M and A which don't include Z 

187 twoDmaskA = numpy.zeros((12), dtype=bool) 

188 for i in [5, 6, 7, 9, 10, 11]: 

189 twoDmaskA[i] = True 

190 

191 twoDmaskM = numpy.zeros((12, 12), dtype=bool) 

192 for y in range(12): 

193 for x in range(12): 

194 if twoDmaskA[y] and twoDmaskA[x]: 

195 twoDmaskM[y, x] = True 

196 

197 else: 

198 twoD = False 

199 

200 if interpolationOrder > 1: 

201 # Override interpolator for python for higher than linear 

202 interpolator = "python" 

203 

204 # Automatically calculate margin if none is passed 

205 # Detect default case and calculate maring necessary for a 45deg rotation with no displacement 

206 if margin is None: 

207 if twoD: 

208 # z-margin will be overwritten below 

209 margin = [1 + int(0.1 * min(im1.shape[1:]))] * 3 

210 else: 

211 margin = [1 + int(0.1 * min(im1.shape))] * 3 

212 elif type(margin) == list: 

213 pass 

214 else: 

215 # Make sure margin is an int 

216 margin = int(margin) 

217 margin = [margin] * 3 

218 

219 # Make sure im2 is bigger than im1 and check difference in size 

220 # Get difference in image sizes. This should be positive, since we must always have enough data for im2 interpolation 

221 im1im2sizeDiff = numpy.array(im2.shape) - numpy.array(im1.shape) 

222 

223 # Check im2 is bigger or same size 

224 if (im1im2sizeDiff < 0).any(): 

225 print("\tcorrelate.register(): im2 is smaller than im1 in at least one dimension: im2.shape: {}, im1.shape: {}".format(im2.shape, im1.shape)) 

226 raise ValueError("correlate.register():DimProblem") 

227 

228 # Make sure margin is at least 1 for the gradient calculation 

229 if twoD: 

230 margin[0] = 0 

231 elif min(margin) < 1 and min(im1im2sizeDiff) == 0: 

232 margin = [1] * 3 

233 

234 # Calculate crops -- margin for im2 and more for im1 if it is bigger 

235 # Margin + half the difference in size for im2 -- im1 will start in the middle. 

236 crop2 = ( 

237 slice( 

238 int(im1im2sizeDiff[0] / 2 + margin[0]), 

239 int(im1im2sizeDiff[0] / 2 + im1.shape[0] - margin[0]), 

240 ), 

241 slice( 

242 int(im1im2sizeDiff[1] / 2 + margin[1]), 

243 int(im1im2sizeDiff[1] / 2 + im1.shape[1] - margin[1]), 

244 ), 

245 slice( 

246 int(im1im2sizeDiff[2] / 2 + margin[2]), 

247 int(im1im2sizeDiff[2] / 2 + im1.shape[2] - margin[2]), 

248 ), 

249 ) 

250 

251 # Get subvolume crops from both images -- just the margin for im1 

252 crop1 = ( 

253 slice(int(margin[0]), int(im1.shape[0] - margin[0])), 

254 slice(int(margin[1]), int(im1.shape[1] - margin[1])), 

255 slice(int(margin[2]), int(im1.shape[2] - margin[2])), 

256 ) 

257 

258 # Create im1 crop to shift less data 

259 im1crop = im1[crop1].copy() 

260 

261 # Calculate effective margin 

262 # to calculate displacement divergence 

263 # using max for the margin -- subjective choice 

264 max(margin) + min(im1im2sizeDiff) / 2 

265 # print( "\tcorrelate.register(): realMargin is:", realMargin) 

266 

267 # If live plot is asked for, initialise canvas 

268 if imShowProgress: 

269 import matplotlib.pyplot as plt 

270 

271 # Plot ranges for signed residual 

272 vmin = -im1crop.max() 

273 vmax = im1crop.max() 

274 if not imShowProgressNewFig: 

275 # if imShowProgress == "Z" or imShowProgress == "z": 

276 plt.subplot(3, 3, 1) 

277 plt.axis([im1crop.shape[2], 0, im1crop.shape[1], 0]) 

278 plt.subplot(3, 3, 2) 

279 plt.axis([im1crop.shape[2], 0, im1crop.shape[1], 0]) 

280 plt.subplot(3, 3, 3) 

281 plt.axis([im1crop.shape[2], 0, im1crop.shape[1], 0]) 

282 # if imShowProgress == "Y" or imShowProgress == "y": 

283 plt.subplot(3, 3, 4) 

284 plt.axis([im1crop.shape[2], 0, im1crop.shape[0], 0]) 

285 plt.subplot(3, 3, 5) 

286 plt.axis([im1crop.shape[2], 0, im1crop.shape[0], 0]) 

287 plt.subplot(3, 3, 6) 

288 plt.axis([im1crop.shape[2], 0, im1crop.shape[0], 0]) 

289 # if imShowProgress == "X" or imShowProgress == "x": 

290 plt.subplot(3, 3, 7) 

291 plt.axis([im1crop.shape[1], 0, im1crop.shape[0], 0]) 

292 plt.subplot(3, 3, 8) 

293 plt.axis([im1crop.shape[1], 0, im1crop.shape[0], 0]) 

294 plt.subplot(3, 3, 9) 

295 plt.axis([im1crop.shape[1], 0, im1crop.shape[0], 0]) 

296 plt.ion() 

297 

298 ########################################################### 

299 # Important -- since we're moving im2, initial Phis will be 

300 # pointing the wrong way, they need to be inversed 

301 ########################################################### 

302 # If there is no initial Phi, initalise it and im1defCrop to zero. 

303 if PhiInit is None: 

304 Phi = numpy.eye(4) 

305 im2def = im2.copy() 

306 

307 else: 

308 # 2020-03-17 in isolation from COVID-19 EA and OS: Apparently this changes the PhiInit outside this function, 

309 # Copying into different variable 

310 # Apply binning on displacement 

311 Phi = PhiInit.copy() 

312 

313 # If we're in rigid mode, keep only translations and rotations for this guess 

314 # If you don't do this it goes mad (i.e., rigid updates to non-rigid guess don't seem to work) 

315 if PhiRigid: 

316 Phi = spam.deformation.computeRigidPhi(Phi.copy()) 

317 Phi[0:3, -1] *= PhiInitBinRatio 

318 

319 # invert PhiInit to apply it to im2 

320 try: 

321 PhiInv = numpy.linalg.inv(Phi.copy()) 

322 except numpy.linalg.linalg.LinAlgError: 

323 PhiInv = numpy.eye(4) 

324 

325 # Since we are now using Fcentred for iterations, do nothing 

326 # call decomposePhi to apply PhiInit (calculated on the centre of the image) to the origin (0,0,0) 

327 if interpolator == "C": 

328 im2def = spam.DIC.applyPhi(im2, Phi=PhiInv, interpolationOrder=interpolationOrder) 

329 

330 elif interpolator == "python": 

331 im2def = spam.DIC.applyPhiPython(im2, Phi=PhiInv, interpolationOrder=interpolationOrder) 

332 

333 def computeGradient(im, twoD): 

334 # Function to compute gradients 

335 if twoD: 

336 # If 2D image we have no gradients in the 1st direction 

337 # if verbose: print("Calculating gradients...", end="") 

338 imGradY, imGradX = numpy.gradient(im[0]) 

339 imGradX = imGradX[numpy.newaxis, ...] 

340 imGradY = imGradY[numpy.newaxis, ...] 

341 imGradZ = numpy.zeros_like(imGradX) 

342 # if verbose: print("done") 

343 else: 

344 # if verbose: print("Calculating gradients...", end="") 

345 imGradZ, imGradY, imGradX = numpy.gradient(im) 

346 # if verbose: print("done ") 

347 return imGradZ, imGradY, imGradX 

348 

349 # Apply stationary im1 mask 

350 if im1mask is not None: 

351 im1crop[im1mask[crop1] is False] = numpy.nan 

352 

353 # Initialise iteration variables 

354 iterations = 0 

355 returnStatus = 0 

356 # Big value to start with to ensure the first iteration 

357 deltaPhiNorm = 100.0 

358 error = _errorCalc(im1crop, im2def[crop2]) 

359 

360 if verbose: 

361 print("Start correlation with Error = {:0.2f}".format(error)) 

362 

363 widgets = [ 

364 " Iteration Number:", 

365 progressbar.Counter(), 

366 " ", 

367 progressbar.FormatLabel(""), 

368 " (", 

369 progressbar.Timer(), 

370 ")", 

371 ] 

372 pbar = progressbar.ProgressBar(widgets=widgets, maxval=maxIterations) 

373 # widgets = [progressbar.FormatLabel(''), ' ', progressbar.Bar(), ' ', progressbar.AdaptiveETA()] 

374 # pbar = progressbar.ProgressBar(widgets=widgets, maxval=numberOfNodes) 

375 pbar.start() 

376 

377 # --- Start Iterations --- 

378 while iterations < maxIterations and deltaPhiNorm > deltaPhiMin: 

379 errorPrev = error 

380 

381 # On first iteration, compute hessian and jacobian in any case 

382 # ...or if we've been asked to update Gradient 

383 if iterations == 0 or updateGradient: 

384 # If we've been asked to update gradient, compute it on im2, which is moving 

385 if updateGradient: 

386 imGradZ, imGradY, imGradX = computeGradient(im2def, twoD) 

387 crop = crop2 

388 # Otherwise compute it once and for all on the non-moving im1 

389 else: 

390 imGradZ, imGradY, imGradX = computeGradient(im1, twoD) 

391 crop = crop1 

392 

393 M = numpy.zeros((12, 12), dtype="<f8") 

394 A = numpy.zeros((12), dtype="<f8") 

395 

396 # Compute both DIC operators A and M with C library 

397 computeDICoperators( 

398 im1crop.astype("<f4"), 

399 im2def[crop2].astype("<f4"), 

400 imGradZ[crop].astype("<f4"), 

401 imGradY[crop].astype("<f4"), 

402 imGradX[crop].astype("<f4"), 

403 M, 

404 A, 

405 ) 

406 else: 

407 # just update jacobian A 

408 A = numpy.zeros((12), dtype="<f8") 

409 computeDICjacobian( 

410 im1crop.astype("<f4"), 

411 im2def[crop2].copy().astype("<f4"), 

412 imGradZ[crop].copy().astype("<f4"), 

413 imGradY[crop].copy().astype("<f4"), 

414 imGradX[crop].copy().astype("<f4"), 

415 A, 

416 ) 

417 

418 # Solve for delta Phi 

419 if twoD: 

420 # If a twoD image, cut out the bits of the M and A matrices that interest us 

421 # This is necessary since the rest is super singular 

422 # Solve for delta Phi 

423 try: 

424 deltaPhi = numpy.dot(numpy.linalg.inv(M[twoDmaskM].reshape(6, 6)), A[twoDmaskA]) 

425 except numpy.linalg.linalg.LinAlgError: 

426 singular = True 

427 break 

428 # ...and now put deltaPhi components back in place for a 3D deltaPhi 

429 deltaPhinew = numpy.zeros((12), dtype=float) 

430 deltaPhinew[twoDmaskA] = deltaPhi 

431 del deltaPhi 

432 deltaPhi = deltaPhinew 

433 else: 

434 # Solve for delta Phi 

435 try: 

436 deltaPhi = numpy.dot(numpy.linalg.inv(M), A) 

437 except numpy.linalg.linalg.LinAlgError: 

438 singular = True 

439 break 

440 

441 # If we're doing a rigid registration... 

442 if PhiRigid: 

443 # Add padding zeros 

444 deltaPhi = numpy.hstack([deltaPhi, numpy.zeros(4)]).reshape((4, 4)) 

445 

446 deltaPhiPlusI = numpy.eye(4) + deltaPhi 

447 # Keep only rigid part of deltaPhi 

448 deltaPhiPlusIrigid = spam.deformation.computeRigidPhi(deltaPhiPlusI.copy()) 

449 

450 # Subtract I from the rigid dPhi+1, and compute norm only on first 3 rows 

451 # ...basically recompute deltaPhiNorm only on rigid part 

452 deltaPhiNorm = numpy.linalg.norm((deltaPhiPlusIrigid - numpy.eye(4))[0:3].ravel()) 

453 

454 # Apply Delta Phi correction to Phi In Roux X-N paper equation number 11 

455 Phi = numpy.dot(Phi, deltaPhiPlusIrigid) 

456 

457 else: 

458 # The general, non-rigid case 

459 deltaPhiNorm = numpy.linalg.norm(deltaPhi) 

460 

461 # Add padding zeros 

462 deltaPhi = numpy.hstack([deltaPhi, numpy.zeros(4)]).reshape((4, 4)) 

463 

464 # Update Phi 

465 Phi = numpy.dot(Phi, (numpy.eye(4) + deltaPhi)) 

466 

467 try: 

468 PhiInv = numpy.linalg.inv(Phi.copy()) 

469 except numpy.linalg.linalg.LinAlgError: 

470 singular = True 

471 break 

472 

473 # reset im1def as emtpy matrix for deformed image 

474 if interpolator == "C": 

475 im2def = spam.DIC.applyPhi(im2, Phi=PhiInv, interpolationOrder=interpolationOrder) 

476 elif interpolator == "python": 

477 im2def = spam.DIC.applyPhiPython(im2, Phi=PhiInv, interpolationOrder=interpolationOrder) 

478 

479 # Error calculation 

480 error = _errorCalc(im1crop, im2def[crop2]) 

481 

482 # Keep interested people up to date with what's happening 

483 # if verbose: 

484 # print("Error = {:0.2f}".format(error)), 

485 # print("deltaPhiNorm = {:0.4f}".format(deltaPhiNorm)) 

486 

487 # Catch divergence condition after half of the max iterations 

488 if errorPrev < error * 0.8 and iterations > maxIterations / 2: 

489 # undo this bad Phi which has increased the error: 

490 # Phi = numpy.dot((numpy.eye(4) + deltaPhi), Phi) 

491 returnStatus = -1 

492 if verbose: 

493 print("\t -> diverging on error condition") 

494 break 

495 

496 # Second divergence criterion on displacement (Issue #62) 

497 # If any displcement is bigger than 5* the margin... 

498 # if (numpy.abs(spam.deformation.decomposePhi(Phi.copy())['t']) > 5 * realMargin).any(): 

499 # if verbose: print("\t -> diverging on displacement condition") 

500 # returnStatus = -3 

501 # break 

502 

503 # 2018-10-02 - EA: Add divergence condition on U 

504 trans = spam.deformation.decomposePhi(Phi.copy()) 

505 try: 

506 volumeChange = numpy.linalg.det(trans["U"]) 

507 if volumeChange > 3 or volumeChange < 0.2: 

508 if verbose: 

509 print("\t -> diverging on volumetric change condition") 

510 returnStatus = -3 

511 break 

512 except Exception: 

513 returnStatus = -3 

514 break 

515 

516 if imShowProgress: 

517 # if imShowProgress == "Z" or imShowProgress == "z": 

518 if imShowProgressNewFig: 

519 plt.figure() 

520 else: 

521 plt.clf() 

522 plt.suptitle("Iteration Number = {}".format(iterations), fontsize=10) 

523 plt.subplot(3, 3, 1) 

524 plt.title("im1 Z-slice") 

525 plt.imshow(im1crop[im1crop.shape[0] // 2, :, :], cmap="Greys_r", vmin=0, vmax=vmax) 

526 plt.subplot(3, 3, 2) 

527 plt.title("im2def Z-slice") 

528 plt.imshow( 

529 im2def[crop2][im1crop.shape[0] // 2, :, :], 

530 cmap="Greys_r", 

531 vmin=0, 

532 vmax=vmax, 

533 ) 

534 plt.subplot(3, 3, 3) 

535 plt.title("im1-im2def Z-slice") 

536 plt.imshow( 

537 numpy.subtract(im1crop, im2def[crop2])[im1crop.shape[0] // 2, :, :], 

538 cmap="coolwarm", 

539 vmin=vmin, 

540 vmax=vmax, 

541 ) 

542 # if imShowProgress == "Y" or imShowProgress == "y": 

543 # if imShowProgressNewFig: plt.figure() 

544 # else: plt.clf() 

545 plt.subplot(3, 3, 4) 

546 plt.title("im1 Y-slice") 

547 plt.imshow(im1crop[:, im1crop.shape[1] // 2, :], cmap="Greys_r", vmin=0, vmax=vmax) 

548 plt.subplot(3, 3, 5) 

549 plt.title("im2def Y-slice") 

550 plt.imshow( 

551 im2def[crop2][:, im1crop.shape[1] // 2, :], 

552 cmap="Greys_r", 

553 vmin=0, 

554 vmax=vmax, 

555 ) 

556 plt.subplot(3, 3, 6) 

557 plt.title("im1-im2def Y-slice") 

558 plt.imshow( 

559 numpy.subtract(im1crop, im2def[crop2])[:, im1crop.shape[1] // 2, :], 

560 cmap="coolwarm", 

561 vmin=vmin, 

562 vmax=vmax, 

563 ) 

564 # if imShowProgress == "X" or imShowProgress == "x": 

565 # if imShowProgressNewFig: plt.figure() 

566 # else: plt.clf() 

567 plt.subplot(3, 3, 7) 

568 plt.title("im1 X-slice") 

569 plt.imshow(im1crop[:, :, im1crop.shape[2] // 2], cmap="Greys_r", vmin=0, vmax=vmax) 

570 plt.subplot(3, 3, 8) 

571 plt.title("im2def X-slice") 

572 plt.imshow( 

573 im2def[crop2][:, :, im1crop.shape[2] // 2], 

574 cmap="Greys_r", 

575 vmin=0, 

576 vmax=vmax, 

577 ) 

578 plt.subplot(3, 3, 9) 

579 plt.title("im1-im2def X-slice") 

580 plt.imshow( 

581 numpy.subtract(im1crop, im2def[crop2])[:, :, im1crop.shape[2] // 2], 

582 cmap="coolwarm", 

583 vmin=vmin, 

584 vmax=vmax, 

585 ) 

586 plt.pause(0.01) 

587 

588 iterations += 1 

589 

590 if verbose: 

591 decomposedPhi = spam.deformation.decomposePhi(Phi.copy()) 

592 widgets[3] = progressbar.FormatLabel( 

593 " dPhiNorm={:0>7.5f} error={:0>4.2f} t=[{:0>5.3f} {:0>5.3f} {:0>5.3f}] r=[{:0>5.3f} {:0>5.3f} {:0>5.3f}] z=[{:0>5.3f} {:0>5.3f} {:0>5.3f}]".format( 

594 deltaPhiNorm, 

595 error, 

596 decomposedPhi["t"][0], 

597 decomposedPhi["t"][1], 

598 decomposedPhi["t"][2], 

599 decomposedPhi["r"][0], 

600 decomposedPhi["r"][1], 

601 decomposedPhi["r"][2], 

602 decomposedPhi["z"][0], 

603 decomposedPhi["z"][1], 

604 decomposedPhi["z"][2], 

605 ) 

606 ) 

607 pbar.update(iterations) 

608 

609 # Positive return status is a healthy end of while loop: 

610 if iterations >= maxIterations: 

611 returnStatus = 1 

612 if deltaPhiNorm <= deltaPhiMin: 

613 returnStatus = 2 

614 if singular: 

615 returnStatus = -2 

616 

617 if verbose: 

618 print() 

619 # pbar.finish() 

620 if iterations > maxIterations: 

621 print("\t -> No convergence before max iterations") 

622 if deltaPhiNorm <= deltaPhiMin: 

623 print("\t -> Converged") 

624 if singular: 

625 print("\t -> Singular") 

626 

627 if im1mask is not None: 

628 # If a mask on im1 is defined, return an Phi at the centre of the mass 

629 maskCOM = spam.label.centresOfMass(im1mask[crop1])[-1] 

630 # print("Mask COM", maskCOM) 

631 # print( "\nNormal Phi:\n", Phi) 

632 Phi[0:3, -1] = spam.deformation.decomposePhi( 

633 Phi.copy(), 

634 PhiCentre=(numpy.array(im1crop.shape) - 1) / 2.0, 

635 PhiPoint=maskCOM, 

636 )["t"] 

637 # print( "\nF in mask:\n", F) 

638 

639 return { 

640 "error": error, 

641 "Phi": Phi, 

642 "returnStatus": returnStatus, 

643 "iterations": iterations, 

644 "deltaPhiNorm": deltaPhiNorm, 

645 } 

646 

647 

648def registerMultiscale( 

649 im1, 

650 im2, 

651 binStart, 

652 binStop=1, 

653 im1mask=None, 

654 PhiInit=None, 

655 PhiRigid=False, 

656 PhiInitBinRatio=1.0, 

657 margin=None, 

658 maxIterations=100, 

659 deltaPhiMin=0.0001, 

660 updateGradient=False, 

661 interpolationOrder=1, 

662 interpolator="C", 

663 verbose=False, 

664 imShowProgress=False, 

665 forceChangeScale=False, 

666): 

667 """ 

668 Perform multiscale subpixel image correlation between im1 and im2. 

669 

670 This means applying a downscale (binning) to the images, performing a Lucas and Kanade at that level, 

671 and then improving it on a 2* less downscaled image, all the way back to the full scale image. 

672 

673 If your input images have multiple scales of texture, this should save significant time. 

674 

675 Please see the documentation for `register` for the rest of the documentation. 

676 

677 Parameters 

678 ---------- 

679 im1 : 3D numpy array 

680 The greyscale image that will not move -- must not contain NaNs 

681 

682 im2 : 3D numpy array 

683 The greyscale image that will be deformed -- must not contain NaNs 

684 

685 binStart : int 

686 Maximum amount of binning to apply, please input a number which is 2^int 

687 

688 binStop : int, optional 

689 Which binning level to stop upscaling at. 

690 The value of 1 (full image resolution) is almost always recommended (unless memory/time problems). 

691 Default = 1 

692 

693 im1mask : 3D boolean numpy array, optional 

694 A mask for the zone to correlate in im1 with `False` in the zone to not correlate. 

695 Default = None, `i.e.`, correlate all of im1 minus the margin. 

696 If this is defined, the Phi returned is in the centre of mass of the mask 

697 

698 PhiInit : 4x4 numpy array, optional 

699 Initial deformation to apply to im1, by default at bin1 scale 

700 Default = numpy.eye(4), `i.e.`, no transformation 

701 

702 PhiRigid : bool, optional 

703 Run a rigid correlation? Only the rigid part of your PhiInit will be kept. 

704 Default = False 

705 

706 PhiInitBinRatio : float, optional 

707 Change translations in PhiInit, if it's been calculated on a differently-binned image. Default = 1 

708 

709 margin : int, optional 

710 Margin, in pixels, to take in im1. 

711 Can also be a N-component list of ints, representing the margin in ND. 

712 If im2 has the same size as im1 this is strictly necessary to allow space for interpolation and movement 

713 Default = 0 (`i.e.`, 10% of max dimension of im1) 

714 

715 maxIterations : int, optional 

716 Maximum number of quasi-Newton iterations to perform before stopping. Default = 25 

717 

718 deltaPhiMin : float, optional 

719 Smallest change in the norm of Phi (the transformation operator) before stopping. Default = 0.001 

720 

721 updateGradient : bool, optional 

722 Should the gradient of the image be computed (and updated) on the deforming im2? 

723 Default = False (it is computed once on im1) 

724 

725 interpolationOrder : int, optional 

726 Order of the greylevel interpolation for applying Phi to im1 when correlating. Recommended value is 3, but you can get away with 1 for faster calculations. Default = 3 

727 

728 interpolator : string, optional 

729 Which interpolation function to use from `spam`. 

730 Default = 'python'. 'C' is also an option 

731 

732 verbose : bool, optional 

733 Get to know what the function is really thinking, recommended for debugging only. Default = False 

734 

735 imShowProgress : bool, optional 

736 Pop up a window showing a ``imShowProgress`` slice of the image differences (im1-im2) as im1 is progressively deformed. 

737 Default = False 

738 

739 forceChangeScale : bool, optional 

740 Change up a scale even if not converged? 

741 Default = False 

742 

743 Returns 

744 ------- 

745 Dictionary: 

746 

747 'Phi': 4x4 float array 

748 Deformation function defined at the centre of the image 

749 

750 'returnStatus': signed int 

751 Return status from the correlation: 

752 

753 2 : Achieved desired precision in the norm of delta Phi 

754 

755 1 : Hit maximum number of iterations while iterating 

756 

757 -1 : Error is more than 80% of previous error, we're probably diverging 

758 

759 -2 : Singular matrix M cannot be inverted 

760 

761 -3 : Displacement > 5*margin 

762 

763 'error': float 

764 Error float describing mismatch between images, it's the sum of the squared difference divided by the sum of im1 

765 

766 'iterations': int 

767 Number of iterations 

768 """ 

769 # Detect unpadded 2D image first: 

770 if len(im1.shape) == 2: 

771 # pad them 

772 im1 = im1[numpy.newaxis, ...] 

773 im2 = im2[numpy.newaxis, ...] 

774 if im1mask is not None: 

775 im1mask = im1mask[numpy.newaxis, ...] 

776 

777 # Detect 2D images 

778 if im1.shape[0] == 1: 

779 twoD = True 

780 else: 

781 twoD = False 

782 

783 import math 

784 

785 logbinstart = math.log(binStart, 2) 

786 if not logbinstart.is_integer(): 

787 print( 

788 "spam.DIC.registerMultiscale(): You asked for an initial binning of", 

789 binStart, 

790 ",rounding it to ", 

791 end="", 

792 ) 

793 binStart = 2 ** numpy.round(logbinstart) 

794 print(binStart) 

795 

796 logbinstop = math.log(binStop, 2) 

797 if not logbinstop.is_integer(): 

798 print( 

799 "spam.DIC.registerMultiscale(): You asked for a final binning of", 

800 binStop, 

801 ",rounding it to ", 

802 end="", 

803 ) 

804 binStop = 2 ** numpy.round(logbinstop) 

805 print(binStop) 

806 

807 # If there is no initial Phi, initalise it and im1defCrop to zero. 

808 if PhiInit is None: 

809 reg = {"Phi": numpy.eye(4)} 

810 else: 

811 # Apply binning on displacement -- the /2 is to be able to *2 it in the LK call 

812 tmp = PhiInit.copy() 

813 tmp[0:3, -1] *= PhiInitBinRatio / 2.0 / float(binStart) 

814 reg = {"Phi": tmp} 

815 

816 if im1mask is not None: 

817 # Multiply up to 100 so we can apply a threshold below on binning in % 

818 im1mask = im1mask.astype("<u1") * 100 

819 

820 # Sorry... This generates a list of binning levels, if binStart=8 and binStop=2 this will be [8, 4 ,2] 

821 binLevels = 2 ** numpy.arange(math.log(binStart, 2), math.log(binStop, 2) - 1, -1).astype(int) 

822 for binLevel in binLevels: 

823 if verbose: 

824 print( 

825 "spam.DIC.registerMultiscale(): working on binning: ", 

826 binLevel, 

827 ) 

828 if binLevel > 1: 

829 if twoD: 

830 import scipy.ndimage 

831 

832 im1b = scipy.ndimage.zoom(im1[0], 1 / binLevel, order=1) 

833 im2b = scipy.ndimage.zoom(im2[0], 1 / binLevel, order=1) 

834 # repad them 

835 im1b = im1b[numpy.newaxis, ...] 

836 im2b = im2b[numpy.newaxis, ...] 

837 if im1mask is not None: 

838 im1maskb = scipy.ndimage.zoom(im1mask[0], 1 / binLevel, order=1) 

839 im1maskb = im1maskb[numpy.newaxis, ...] 

840 else: 

841 im1maskb = None 

842 else: 

843 im1b = spam.DIC.binning(im1, binLevel) 

844 im2b = spam.DIC.binning(im2, binLevel) 

845 if im1mask is not None: 

846 im1maskb = spam.DIC.binning(im1mask, binLevel) > 0 

847 else: 

848 im1maskb = None 

849 else: 

850 im1b = im1 

851 im2b = im2 

852 if im1mask is not None: 

853 im1maskb = im1mask > 0 

854 else: 

855 im1maskb = None 

856 

857 # Automatically calculate margin if none is passed 

858 # Detect default case and calculate margin necessary for a 45deg rotation with no displacement 

859 if margin is None: 

860 if twoD: 

861 # z-margin will be overwritten below 

862 marginB = [1 + int(0.1 * min(im1b.shape[1:]))] * 3 

863 else: 

864 marginB = [1 + int(0.1 * min(im1b.shape))] * 3 

865 

866 elif type(margin) == list: 

867 marginB = (numpy.array(margin) // binLevel).tolist() 

868