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

1""" 

2Library of SPAM multimodal image correlation functions using a Gaussian Mixture model from Tudisco et al. 2017. 

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 

19# 2017-05-29 ER and EA 

20 

21import matplotlib as mpl 

22import matplotlib.pyplot as plt 

23import numpy 

24import scipy.ndimage 

25import spam.deformation 

26import spam.DIC 

27import spam.helpers 

28from spam.DIC.DICToolkit import computeDICoperatorsGM, computeGMresidualAndPhase 

29 

30numpy.set_printoptions(precision=3, suppress=True) 

31mpl.rc("font", size=6) 

32cmapPhases = "Set1_r" 

33 

34 

35def multimodalRegistration( 

36 im1, 

37 im2, 

38 phaseDiagram, 

39 gaussianParameters, 

40 im1mask=None, 

41 PhiInit=None, 

42 margin=None, 

43 maxIterations=50, 

44 deltaPhiMin=0.005, 

45 interpolationOrder=1, 

46 verbose=False, 

47 GRAPHS=False, 

48 INTERACTIVE=False, 

49 sliceAxis=0, 

50 suffix="", 

51 rootPath=".", 

52 BINS=64, 

53): 

54 """ 

55 Perform subpixel image correlation between im1 and im2 -- images of the same object acquired with different modalities. 

56 

57 This function will deform im2 until it best matches im1. 

58 The matching includes sub-pixel displacements, rotation, and linear straining of the whole image. 

59 The correlation of im1, im2 will give a deformationFunction F which maps im1 into im2. 

60 Phi(im1(x),im2(F.x)) = 0 

61 As per Equation (3) of Tudisco `et al.` 

62 

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

64 

65 Parameters 

66 ---------- 

67 im1 : 3D numpy array 

68 The greyscale image that will not move 

69 

70 im2 : 3D numpy array 

71 The greyscale image that will be deformed 

72 

73 phaseDiagram : Nbins x Nbins numpy array of ints 

74 Pre-labelled phase diagram, which is a joint histogram with the phases labelled 

75 

76 gaussianParameters : 2D numpy array Nx6 

77 Key parameter which is the result of a 2D gaussian fit of the first N peaks in the joint histograms of 

78 im1 and im2. 

79 The 6 parameters of the fit are: 

80 φ, μ(im1), μ(im2), and a, b, c, where {a,b,c} are the parameter that can be found for `two-dimensional elliptical Gaussian function` here: 

81 https://en.wikipedia.org/wiki/Gaussian_function, `i.e.`, coupled with im1², im1*im2 and im2² respectively 

82 

83 im1mask : 3D numpy array, float, optional 

84 A mask for the zone to correlate in im1 with NaNs in the zone to not correlate. Default = None, `i.e.`, correlate all of im1 minus the margin 

85 

86 PhiInit : 4x4 numpy array, optional 

87 Initial transformation to apply. Default = numpy.eye(4), `i.e.`, no transformation 

88 

89 margin : int, optional 

90 Margin, in pixels, to take in im1 and im2 to allow space for interpolation and movement. 

91 Default = None (`i.e.`, enough margin for a worse-case 45degree rotation with no displacement) 

92 

93 maxIterations : int, optional 

94 Maximum number of quasi-Newton interations to perform before stopping. Default = 25 

95 

96 deltaPhiMin : float, optional 

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

98 

99 interpolationOrder : int, optional 

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

101 

102 verbose : bool, optional 

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

104 

105 GRAPHS : bool, optional 

106 Pop up a window showing the image differences (im1-im2) as im1 is progressively deformed. 

107 

108 Returns 

109 -------- 

110 Dictionary: 

111 

112 'transformation' 

113 Dictionary containing: 

114 

115 't' : 3-component vector 

116 

117 Z, Y, X displacements in pixels 

118 

119 'r' : 3-component vector 

120 

121 Z, Y, X components of rotation vector 

122 

123 'Phi': 4 x 4 numpy array of floats 

124 Deformation function, Phi 

125 

126 'returnStatus': int 

127 Return status from the correlation: 

128 

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

130 

131 1 : Hit maximum number of iterations while iterating 

132 

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

134 

135 -2 : Singular matrix M cannot be inverted 

136 

137 -3 : Displacement > 5*margin 

138 

139 'iterations': int 

140 Number of iterations 

141 

142 'logLikelyhood' : float 

143 Number indicating quality of match 

144 

145 'deltaPhiNorm' : float 

146 Size of last Phi step 

147 

148 'residualField': 3D numpy array of floats 

149 Same size as input image, residual field 

150 

151 'phaseField': phaseField 

152 Same size as input image, labelled phases 

153 

154 Note 

155 ---- 

156 This correlation is what is proposed in Tudisco et al. "An extension of Digital Image Correlation for intermodality image registration", section 4.3. 

157 """ 

158 

159 # if verbose: 

160 # print("Enter registration") 

161 

162 # for interactive graphs 

163 if INTERACTIVE: 

164 GRAPHS = True 

165 plt.ion() 

166 

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

168 if margin is None: 

169 # sqrt 

170 margin = int((3**0.5 - 1.0) * max(im1.shape) * 0.5) 

171 else: 

172 # Make sure margin is an int 

173 margin = int(margin) 

174 

175 # Exit clause for singular matrices 

176 singular = False 

177 

178 crop = ( 

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

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

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

182 ) 

183 

184 # Figure out coordinates on which the correlation should happen, i.e., the non NaN ones 

185 # NOTE: these coordinates are within the cropped margin for interpolation 

186 # cropImCoordsFG = numpy.where( numpy.isfinite( im1[crop] ) ) 

187 if im1mask is not None: 

188 # TODO: This could just be directly = mask and equals inv of mask 

189 # i.e., arry a boolean mask and not a series of coord array 

190 # cropImCoordsFG = numpy.where(im1mask[crop] is True) 

191 # cropImCoordsBG = numpy.where(im1mask[crop] is False) 

192 pass 

193 else: 

194 # Everything regular, except that we might have NaNs in the CW... 

195 # cropImCoordsFG = numpy.ones_like( im1[crop], dtype='bool' ) 

196 # cropImCoordsBG = numpy.zeros_like( im1[crop], dtype='bool' ) 

197 # cropImCoordsFG = numpy.isfinite(im1[crop], dtype="bool") 

198 # cropImCoordsBG = numpy.isnan(im1[crop], dtype="bool") 

199 # HACK 

200 # print cropImCoordsFG 

201 # print cropImCoordsBG 

202 # HACK: set nans in im2 to zero 

203 im2[numpy.isnan(im2)] = 0.0 

204 # if numpy.isnan( im1[crop] ).sum() > 0 numpy.isnan( im2[crop] ).sum() > 0: 

205 

206 # compute image centre 

207 # im1centre = (numpy.array(im1.shape)-1)/2.0 

208 

209 # If there is no initial Phi, initalise it to zero. 

210 if PhiInit is None: 

211 Phi = numpy.eye(4) 

212 im2def = im2.copy() 

213 

214 # If there is an initial guess... 

215 else: 

216 Phi = PhiInit.copy() 

217 

218 # invert PhiInit to apply it to im2 

219 try: 

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

221 except numpy.linalg.linalg.LinAlgError: 

222 PhiInv = numpy.eye(4) 

223 

224 im2def = spam.DIC.applyPhi(im2, Phi=PhiInv, interpolationOrder=interpolationOrder).astype(im2.dtype) 

225 

226 # At this stage we've computed gradients which we are not going to update, im2 and it's gradients will be set equal to 

227 # their cropped versions: 

228 # im2def = im2def[crop] 

229 # im2defGradZ = im2defGradZ[crop] 

230 # im2defGradY = im2defGradY[crop] 

231 # im2defGradX = im2defGradX[crop] 

232 

233 # Mask nans in gradient (but not before or there are jumps when the grey goes to zero 

234 # beacuse of this a label dilate of at least 1 is recommended) 

235 # im2defGradZ[ numpy.where( numpy.isnan( im1gradZ ) ) ] = 0 

236 # im2defGradY[ numpy.where( numpy.isnan( im1gradY ) ) ] = 0 

237 # im2defGradX[ numpy.where( numpy.isnan( im1gradX ) ) ] = 0 

238 

239 iterations = 0 

240 returnStatus = 0 

241 deltaPhiNorm = 0.0 

242 

243 # compute initial Log likelyhood 

244 # p, _, _ = numpy.histogram2d(im1[crop].ravel(), im2def[crop].ravel(), bins=BINS, range=[[0.0, BINS], [0.0, BINS]], normed=False, weights=None) 

245 # LL = 0.0 

246 # for v in p.ravel(): 

247 # if v: 

248 # LL += numpy.log(v) 

249 # LL = numpy.prod(p[numpy.where(p > 0)]) 

250 

251 p, _, _ = numpy.histogram2d( 

252 im1[crop].ravel(), 

253 im2def[crop].ravel(), 

254 bins=BINS, 

255 range=[[0.0, BINS], [0.0, BINS]], 

256 # normed=False, 

257 weights=None, 

258 ) 

259 LL = numpy.sum(numpy.log(p[numpy.where(p > 0)])) 

260 

261 if verbose: 

262 print("\tInitial state LL = {:0.2f}".format(LL)) 

263 print("\tIteration Number {:03d} ".format(iterations), end="") 

264 print("LL = {:0.2f} ".format(LL), end="") 

265 print("dPhi = {:0.4f} ".format(deltaPhiNorm), end="") 

266 # print("\nPhi", Phi ) 

267 # currentTransformation = spam.deformation.decomposePhi(Phi, PhiPoint=im1centre) 

268 currentTransformation = spam.deformation.decomposePhi(Phi) 

269 print( 

270 "Tr = {: .3f}, {: .3f}, {: .3f} ".format(*currentTransformation["t"]), 

271 end="", 

272 ) 

273 print( 

274 "Ro = {: .3f}, {: .3f}, {: .3f} ".format(*currentTransformation["r"]), 

275 end="", 

276 ) 

277 print("Zo = {: .3f}, {: .3f}, {: .3f} ".format(*currentTransformation["z"])) 

278 

279 # Use gradient of image 1 which does not move 

280 # im1GradZ, im1GradY, im1GradX = [g.astype('<f4') for g in numpy.gradient(im1)] 

281 # im2GradZ, im2GradY, im2GradX = [g.astype('<f4') for g in numpy.gradient(im2def)] 

282 

283 while (iterations <= maxIterations and deltaPhiNorm > deltaPhiMin) or iterations == 0: 

284 previousLL = LL 

285 

286 if verbose: 

287 print("\tIteration Number {:03d} ".format(iterations + 1), end="") 

288 

289 im2defGradZ, im2defGradY, im2defGradX = (g.astype("<f4") for g in numpy.gradient(im2def)) 

290 

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

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

293 # im2 updated 

294 computeDICoperatorsGM( 

295 im1[crop].astype("<f4"), 

296 im2def[crop].astype("<f4"), 

297 im2defGradZ[crop].astype("<f4"), 

298 im2defGradY[crop].astype("<f4"), 

299 im2defGradX[crop].astype("<f4"), 

300 phaseDiagram.astype("<u1"), 

301 gaussianParameters.astype("<f8"), 

302 M, 

303 A, 

304 ) 

305 # im2 with guess 

306 # computeDICoperatorsGM(im1[crop].astype("<f4"), 

307 # im2def[crop].astype("<f4"), 

308 # im2GradZ[crop].astype("<f4"), 

309 # im2GradY[crop].astype("<f4"), 

310 # im2GradX[crop].astype("<f4"), 

311 # phaseDiagram.astype("<u1"), 

312 # gaussianParameters.astype("<f8"), 

313 # M, A) 

314 # im1 

315 # computeDICoperatorsGM(im1[crop].astype("<f4"), 

316 # im2def[crop].astype("<f4"), 

317 # im1GradZ[crop].astype("<f4"), 

318 # im1GradY[crop].astype("<f4"), 

319 # im1GradX[crop].astype("<f4"), 

320 # phaseDiagram.astype("<u1"), 

321 # gaussianParameters.astype("<f8"), 

322 # M, A) 

323 

324 try: 

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

326 

327 except numpy.linalg.linalg.LinAlgError: # no cover 

328 singular = True 

329 # TODO: Calculate error for clean break. 

330 print("\tsingular M matrix") 

331 print("exiting") 

332 exit() 

333 # break 

334 

335 deltaPhiNorm = numpy.linalg.norm(deltaPhi) 

336 

337 # Add padding zeros 

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

339 

340 # In Roux X-N paper equation number 11 

341 # Phi = numpy.dot((numpy.eye(4) - deltaPhi), Phi) 

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

343 

344 # currentTransformation = spam.deformation.decomposePhi(Phi, PhiPoint=im1centre) 

345 currentTransformation = spam.deformation.decomposePhi(Phi) 

346 

347 # Solve for delta Phi 

348 try: 

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

350 except numpy.linalg.linalg.LinAlgError: 

351 singular = True 

352 break 

353 

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

355 

356 residualField = numpy.zeros_like(im2[crop], dtype="<f4") 

357 phaseField = numpy.zeros_like(im2[crop], dtype="<u1") 

358 

359 computeGMresidualAndPhase( 

360 im1[crop].astype("<f4"), 

361 im2def[crop].astype("<f4"), 

362 phaseDiagram.astype("<u1"), 

363 gaussianParameters.astype("<f8"), 

364 residualField, 

365 phaseField, 

366 ) 

367 

368 # computeGMresidualAndPhase(im1[crop].astype("<f4"), 

369 # im2def[crop].astype("<f4"), 

370 # phaseDiagram.copy().astype("<u1"), 

371 # gaussianParameters.copy().astype("<f8"), 

372 # residualField, 

373 # phaseField) 

374 

375 # compute current log likelyhood 

376 # p, _, _ = numpy.histogram2d(im1[crop].ravel(), im2def[crop].ravel(), bins=BINS, range=[[0.0, BINS], [0.0, BINS]], normed=False, weights=None) 

377 # LL = 0.0 

378 # for v in p.ravel(): 

379 # if v: 

380 # LL += numpy.log(v) 

381 p, _, _ = numpy.histogram2d( 

382 im1[crop].ravel(), 

383 im2def[crop].ravel(), 

384 bins=BINS, 

385 range=[[0.0, BINS], [0.0, BINS]], 

386 # normed=False, 

387 weights=None, 

388 ) 

389 LL = numpy.sum(numpy.log(p[numpy.where(p > 0)])) 

390 

391 if verbose: 

392 print("LL = {:0.2f} ".format(LL), end="") 

393 print("dPhi = {:0.4f} ".format(deltaPhiNorm), end="") 

394 # print("\nPhi", Phi ) 

395 print( 

396 "Tr = {: .3f}, {: .3f}, {: .3f} ".format(*currentTransformation["t"]), 

397 end="", 

398 ) 

399 print( 

400 "Ro = {: .3f}, {: .3f}, {: .3f} ".format(*currentTransformation["r"]), 

401 end="", 

402 ) 

403 print("Zo = {: .3f}, {: .3f}, {: .3f} ".format(*currentTransformation["z"])) 

404 

405 if previousLL < LL * 0.6: 

406 # undo this bad Phi which has increased the LL: 

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

408 returnStatus = -1 

409 print("Log-likelyhood increasing, divergence condition detected, exiting.") 

410 # break 

411 print("...no actually continuing...") 

412 

413 if GRAPHS: 

414 NPHASES = gaussianParameters.shape[0] 

415 grid = plt.GridSpec(2, 4) 

416 plt.clf() 

417 plt.suptitle( 

418 "Gaussian Mixture {} iteration number {} " 

419 "|deltaPhi|={:.5f} \nT = [{: 2.4f} {: 2.4f} {:.4f}]\n" 

420 "R = [{: 2.4f} {: 2.4f} {: 2.4f}]\n" 

421 "Z = [{: 2.4f} {: 2.4f} {: 2.4f}]".format( 

422 suffix, 

423 iterations, 

424 deltaPhiNorm, 

425 currentTransformation["t"][0], 

426 currentTransformation["t"][1], 

427 currentTransformation["t"][2], 

428 currentTransformation["r"][0], 

429 currentTransformation["r"][1], 

430 currentTransformation["r"][2], 

431 currentTransformation["z"][0], 

432 currentTransformation["z"][1], 

433 currentTransformation["z"][2], 

434 ) 

435 ) 

436 

437 plt.subplot(grid[0, 0]) 

438 plt.axis("off") 

439 plt.title("Residual field") 

440 if sliceAxis == 0: 

441 plt.imshow(residualField[residualField.shape[0] // 2, :, :]) 

442 elif sliceAxis == 1: 

443 plt.imshow(residualField[:, residualField.shape[1] // 2, :]) 

444 elif sliceAxis == 2: 

445 plt.imshow(residualField[:, :, residualField.shape[2] // 2]) 

446 # plt.colorbar() 

447 

448 plt.subplot(grid[0, 1]) 

449 plt.axis("off") 

450 plt.title("Phase field") 

451 if sliceAxis == 0: 

452 plt.imshow( 

453 phaseField[phaseField.shape[0] // 2, :, :], 

454 vmin=-0.5, 

455 vmax=NPHASES + 0.5, 

456 cmap=mpl.cm.get_cmap(cmapPhases, NPHASES + 1), 

457 ) 

458 elif sliceAxis == 1: 

459 plt.imshow( 

460 phaseField[:, phaseField.shape[1] // 2, :], 

461 vmin=-0.5, 

462 vmax=NPHASES + 0.5, 

463 cmap=mpl.cm.get_cmap(cmapPhases, NPHASES + 1), 

464 ) 

465 elif sliceAxis == 2: 

466 plt.imshow( 

467 phaseField[:, :, phaseField.shape[2] // 2], 

468 vmin=-0.5, 

469 vmax=NPHASES + 0.5, 

470 cmap=mpl.cm.get_cmap(cmapPhases, NPHASES + 1), 

471 ) 

472 # plt.colorbar(ticks=numpy.arange(0, NPHASES+1)) 

473 

474 plt.subplot(grid[1, 0]) 

475 plt.axis("off") 

476 plt.title("im1") 

477 if sliceAxis == 0: 

478 plt.imshow(im1[crop][im1[crop].shape[0] // 2, :, :]) 

479 elif sliceAxis == 1: 

480 plt.imshow(im1[crop][:, im1[crop].shape[1] // 2, :]) 

481 elif sliceAxis == 2: 

482 plt.imshow(im1[crop][:, :, im1[crop].shape[2] // 2]) 

483 # plt.colorbar() 

484 

485 plt.subplot(grid[1, 1]) 

486 plt.axis("off") 

487 plt.title("im2def") 

488 if sliceAxis == 0: 

489 plt.imshow(im2def[crop][im2def[crop].shape[0] // 2, :, :]) 

490 elif sliceAxis == 1: 

491 plt.imshow(im2def[crop][:, im2def[crop].shape[1] // 2, :]) 

492 elif sliceAxis == 2: 

493 plt.imshow(im2def[crop][:, :, im2def[crop].shape[2] // 2]) 

494 # plt.colorbar() 

495 

496 plt.subplot(grid[:, 2:]) 

497 plt.axis("off") 

498 plt.title("Checker Board") 

499 if sliceAxis == 0: 

500 plt.imshow( 

501 spam.helpers.checkerBoard( 

502 im1[crop][im2def[crop].shape[0] // 2, :, :], 

503 im2def[crop][im2def[crop].shape[0] // 2, :, :], 

504 ) 

505 ) 

506 elif sliceAxis == 1: 

507 plt.imshow( 

508 spam.helpers.checkerBoard( 

509 im1[crop][:, im2def[crop].shape[1] // 2, :], 

510 im2def[crop][:, im2def[crop].shape[1] // 2, :], 

511 ) 

512 ) 

513 elif sliceAxis == 2: 

514 plt.imshow( 

515 spam.helpers.checkerBoard( 

516 im1[crop][:, :, im2def[crop].shape[2] // 2], 

517 im2def[crop][:, :, im2def[crop].shape[2] // 2], 

518 ) 

519 ) 

520 

521 if INTERACTIVE: 

522 plt.show() 

523 plt.pause(1.0) 

524 else: 

525 plt.savefig( 

526 "{}/GaussianMixture_Iteration-{}-{}.png".format(rootPath, iterations, suffix), 

527 dpi=600, 

528 ) 

529 

530 iterations += 1 

531 

532 if INTERACTIVE: 

533 plt.ioff() 

534 plt.close() 

535 

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

537 if iterations >= maxIterations: 

538 returnStatus = 1 

539 if deltaPhiNorm <= deltaPhiMin: 

540 returnStatus = 2 

541 if singular: 

542 returnStatus = -2 

543 

544 if verbose: 

545 print() 

546 # pbar.finish() 

547 if iterations > maxIterations: 

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

549 if deltaPhiNorm <= deltaPhiMin: 

550 print("\t -> Converged") 

551 if singular: 

552 print("\t -> Singular") 

553 

554 return { 

555 "transformation": currentTransformation, 

556 "Phi": Phi, 

557 "returnStatus": returnStatus, 

558 "iterations": iterations, 

559 "residualField": residualField, 

560 "phaseField": phaseField, 

561 "logLikelyhood": LL, 

562 "deltaPhiNorm": deltaPhiNorm, 

563 } 

564 

565 

566################################################################ 

567# helper functions for correlation with different modalities 

568################################################################ 

569 

570 

571def gaussianMixtureParameters( 

572 im1, 

573 im2, 

574 BINS=64, 

575 NPHASES=2, 

576 im1threshold=0, 

577 im2threshold=0, 

578 distanceMaxima=None, 

579 excludeBorder=False, 

580 fitDistance=None, 

581 GRAPHS=False, 

582 INTERACTIVE=False, 

583 sliceAxis=0, 

584 rootPath=".", 

585 suffix="", 

586): 

587 """ 

588 This function, given two different modality images, builds the joint histogram in BINS bins, 

589 and fits NPHASES peaks with bivariate Gaussian distributions. 

590 

591 Parameters 

592 ---------- 

593 im1 : 3D numpy array of floats 

594 One image, values should be rescaled to 0:BIN-1 

595 

596 im2 : 3D numpy array of floats 

597 Another image, values should be rescaled to 0:BIN-1 

598 

599 BINS : int, optional 

600 Number of bins for the joint histogram, recommend 2^x 

601 Default = 64 

602 

603 NPHASES : int, optional 

604 Number of phases to automatically fit 

605 Default = 2 

606 

607 im1threshold : float, optional 

608 

609 im2threshold : float, optional 

610 

611 distanceMaxima : float, optional 

612 

613 excludeBorder : False, optional 

614 Exclude edges in peak finding? This is passed to skimage.feature.peak_local_max() 

615 

616 fitDistance : float, optional 

617 

618 GRAPHS : bool, optional 

619 

620 INTERACTIVE : bool, optional 

621 

622 sliceAxis=0 

623 

624 rootPath="." 

625 

626 suffix="" 

627 

628 Returns 

629 ------- 

630 gaussianParameters : list of lists 

631 List of NPHASES components, containing the parameters of the bivariate Gaussian fit for each phase: 

632 [0] = GLOBALphi -- Normlised "Height" of the fit 

633 [1] = GLOBALmu1 -- Mean in F of bivairate Gaussian 

634 [2] = GLOBALmu2 -- Mean in G of bivairate Gaussian 

635 [3] = a -- 

636 [4] = b -- 

637 [5] = c -- 

638 

639 

640 p : BINSxBINS 2D numpy array of floats 

641 The normalised joint histogram, the sum of which will be =1 if all of your image values are 0:BIN-1 

642 """ 

643 

644 # To fit 2D likelyhood with gaussian ellipsoids 

645 # To find maxima in likelyhood which correspond to phases 

646 import skimage.feature 

647 from scipy.optimize import curve_fit 

648 

649 # for interactive graphs 

650 if INTERACTIVE: 

651 GRAPHS = True 

652 plt.ion() 

653 

654 # DEFINE the global variables needed for curve fitting 

655 GLOBALmu1 = 0.0 # mean of the f image 

656 GLOBALmu2 = 0.0 # mean of the g image 

657 GLOBALphi = 0.0 # number of hits (value of the maxima) 

658 

659 # DEFINE fitting functions 

660 # https://en.wikipedia.org/wiki/Gaussian_function#Two-dimensional_Gaussian_function 

661 # def gaussian2Delliptical( XY, GLOBALphi, GLOBALmu2, GLOBALmu1, a, b, c ): 

662 # Thsi needs to be in here in order to pass GLOBALphi as a global variable. 

663 # Perhaps it should be optimised? 

664 def computeLambda(a, b, c, x, GLOBALmu1, y, GLOBALmu2): 

665 return numpy.longfloat(0.5 * (a * (x - GLOBALmu1) ** 2 + 2.0 * b * (x - GLOBALmu1) * (y - GLOBALmu2) + c * (y - GLOBALmu2) ** 2)) 

666 

667 def gaussian2Delliptical(XY, a, b, c): 

668 # invert x and y on purpose to be consistent with H 

669 grid = numpy.array(numpy.meshgrid(XY[1], XY[0])) 

670 field = numpy.zeros(grid.shape[1:3]) 

671 for ny in range(grid.shape[2]): 

672 y = grid[0, 0, ny] 

673 for nx in range(grid.shape[1]): 

674 x = grid[1, nx, 0] 

675 field[nx, ny] = float(GLOBALphi) * numpy.exp(-computeLambda(a, b, c, x, GLOBALmu1, y, GLOBALmu2)) 

676 return field.ravel() 

677 

678 # START function 

679 

680 im1min = im1.min() 

681 im1max = im1.max() 

682 im2min = im2.min() 

683 im2max = im2.max() 

684 

685 # f,g discretisation 

686 X = numpy.linspace(0, BINS - 1, BINS) 

687 Y = numpy.linspace(0, BINS - 1, BINS) 

688 

689 print("\tim1 from {:.2f} to {:.2f}".format(im1min, im1max)) 

690 print("\tim2 from {:.2f} to {:.2f}".format(im2min, im2max)) 

691 

692 # Calculate probability distribution p(f,g) 

693 p, _, _ = numpy.histogram2d( 

694 im1.ravel(), 

695 im2.ravel(), 

696 bins=BINS, 

697 range=[[0.0, BINS], [0.0, BINS]], 

698 # normed=False, 

699 weights=None, 

700 ) 

701 p /= float(len(im1.ravel())) 

702 p_sum = p.sum() 

703 

704 print(f"\tp normalisation: {p_sum:.2f}") 

705 

706 # write joint histogram in a dat file for tikz 

707 # if DATA: 

708 # tmp = p.copy() 

709 # with open("GaussianMixture_jointHistogram-{}-{}.dat".format(0, suffix), "w") as f: 

710 # string = "{}\t {}\t {}\n".format("x", "y", "c") 

711 # f.write(string) 

712 # for xbin in range(tmp.shape[0]): 

713 # x = (xbin+0.5)/tmp.shape[0] 

714 # for ybin in range(tmp.shape[1]): 

715 # y = (ybin+0.5)/tmp.shape[1] 

716 # if tmp[xbin, ybin]: 

717 # string = "{}\t {}\t {}\n".format(x, y, tmp[xbin, ybin]) 

718 # f.write(string) 

719 

720 # Get disordered maxima of the joint histogram 

721 print("\tFind maxima") 

722 if distanceMaxima is None: 

723 distanceMaxima = int(BINS / 25.0) 

724 distanceMaxima = max(1, int(distanceMaxima)) 

725 print(f"\t\tMin distance between maxima: {distanceMaxima:.0f}") 

726 maxima = skimage.feature.peak_local_max(p, min_distance=int(distanceMaxima), exclude_border=excludeBorder) 

727 nMax = len(maxima) 

728 

729 if nMax < NPHASES: 

730 print(f"\t\t{NPHASES} phases asked but only {nMax} maxima...") 

731 NPHASES = min(NPHASES, nMax) 

732 

733 # print "\t\t Number of maxima: {:2}".format(nMax) 

734 if nMax == 0: # no cover 

735 print("In gaussian fit: no maxium found... Stoping...") 

736 exit() 

737 

738 # Organise maxima into muF, muG, p(f,g) the sort at take the maximum 

739 maxTmp = numpy.zeros((nMax, 3)) 

740 for i, (f, g) in enumerate(maxima): 

741 if f > im1threshold or g > im2threshold: 

742 maxTmp[i] = [f, g, p[f, g]] 

743 # print("\t\t max {:2} --- f: {:.2f} g: {:.2f} hits: {}".format(i+1,*maxTmp[i])) 

744 

745 nMax = 0 

746 percentage = 0.1 

747 while nMax < NPHASES: 

748 maxSorted = [m for m in maxTmp[maxTmp[:, 2].argsort()[::-1]] if float(m[2]) > percentage * float(p_sum)] 

749 nMax = len(maxSorted) 

750 print(f"\t\t{nMax:02d} maxima found over the {NPHASES} needed with {percentage:.2e} times of the total count") 

751 percentage /= 10.0 

752 

753 for i, (GLOBALmu1, GLOBALmu2, GLOBALphi) in enumerate(maxSorted): 

754 print(f"\t\tMaximum {i + 1}:\t mu1={GLOBALmu1:.2f}\t mu2={GLOBALmu2:.2f}\t Phi={GLOBALphi:.2e}") 

755 print("") 

756 

757 # output of the function: list of fitting gaussian parameters 

758 # size NPHASES x 6, the 6 parameters being GLOBALlogPhi, GLOBALmu1, GLOBALmu2, a, b, c 

759 gaussianParameters = numpy.zeros((NPHASES, 6)).astype("<f4") 

760 

761 # loop over phases 

762 maxEllipsoid = numpy.zeros_like(p) 

763 for iPhase in range(NPHASES): 

764 GLOBALmu1, GLOBALmu2, GLOBALphi = maxSorted[iPhase] 

765 print(f"\tPhase {iPhase + 1:2}:\t\t mu1={GLOBALmu1:.2f}\t mu2={GLOBALmu2:.2f}\t Phi={GLOBALphi:.2e}") 

766 if fitDistance is None: 

767 fitDistance = BINS / 2.0 

768 

769 # fit the probability distribution p(f,g) with an gaussian ellipsoids 

770 pFit = p.copy() 

771 

772 for nf in range(pFit.shape[0]): 

773 for ng in range(pFit.shape[1]): 

774 posF = nf 

775 posG = ng 

776 dist = numpy.sqrt((posF - GLOBALmu1) ** 2.0 + (posG - GLOBALmu2) ** 2.0) # cicrle 

777 # dist = abs(posF-GLOBALmu1)+abs(posG-GLOBALmu2) # square 

778 if dist > fitDistance: 

779 pFit[nf, ng] = 0.0 

780 

781 (a, b, c), _ = curve_fit(gaussian2Delliptical, (X, Y), pFit.ravel(), p0=(1, 1, 1)) 

782 print(f"\t\tFit:\t\t a={a:.2f}\t b={b:.2f}\t c={c:.2f}\t Hessian: {a * c - b**2:.2f}") 

783 while a * c - b**2 < 0: 

784 print("\t\t\tWarning: Hessian < 0") 

785 print("\t\t\t The gaussian doesn't have a local maximum.") 

786 fitDistance /= 2.0 

787 print(f"\t\t\t Try with fit distance = {fitDistance} ") 

788 

789 for nf in range(pFit.shape[0]): 

790 for ng in range(pFit.shape[1]): 

791 posF = nf / float(pFit.shape[0] - 1) 

792 posG = ng / float(pFit.shape[1] - 1) 

793 dist = numpy.sqrt((posF - GLOBALmu1) ** 2.0 + (posG - GLOBALmu2) ** 2.0) # cicrle 

794 # dist = abs(posF-GLOBALmu1)+abs(posG-GLOBALmu2) # square 

795 if dist > fitDistance: 

796 pFit[nf, ng] = 0.0 

797 

798 (a, b, c), _ = curve_fit(gaussian2Delliptical, (X, Y), pFit.ravel(), p0=(1, 1, 1)) 

799 print(f"\t\tFit:\t\t a={a:.2f}\t b={b:.2f}\t c={c:.2f}\t Hessian: {a * c - b**2:.2f}") 

800 

801 # compute covariance function 

802 cov = scipy.linalg.inv(numpy.array([[a, b], [b, c]])) 

803 print(f"\t\tCov:\t\t 1,1={cov[0, 0]:.4f}\t 1,2={cov[1, 0]:.4f}\t 2,2={cov[1, 1]:.4f}") 

804 

805 gaussianParameters[iPhase, 0] = GLOBALphi 

806 gaussianParameters[iPhase, 1] = GLOBALmu1 

807 gaussianParameters[iPhase, 2] = GLOBALmu2 

808 gaussianParameters[iPhase, 3] = a 

809 gaussianParameters[iPhase, 4] = b 

810 gaussianParameters[iPhase, 5] = c 

811 

812 currentEllipsoid = gaussian2Delliptical((X, Y), a, b, c).reshape((len(X), len(Y))) 

813 

814 maxEllipsoid = numpy.maximum(maxEllipsoid, currentEllipsoid) 

815 

816 if GRAPHS: 

817 plt.clf() 

818 plt.suptitle(f"Gaussian Mixture fitting Phase number {iPhase + 1}") 

819 plt.subplot(221) 

820 plt.title(f"im1 (from {im1.min():.2f} to {im1.max():.2f})") 

821 plt.axis("off") 

822 if sliceAxis == 0: 

823 plt.imshow(im1[im1.shape[0] // 2, :, :], vmin=0.0, vmax=BINS) 

824 elif sliceAxis == 1: 

825 plt.imshow(im1[:, im1.shape[1] // 2, :], vmin=0.0, vmax=BINS) 

826 elif sliceAxis == 2: 

827 plt.imshow(im1[:, :, im1.shape[2] // 2], vmin=0.0, vmax=BINS) 

828 plt.colorbar() 

829 

830 plt.subplot(222) 

831 plt.title(f"im2 (from {im2.min():.2f} to {im2.max():.2f})") 

832 plt.axis("off") 

833 if sliceAxis == 0: 

834 plt.imshow(im2[im2.shape[0] // 2, :, :], vmin=0.0, vmax=BINS) 

835 elif sliceAxis == 1: 

836 plt.imshow(im2[:, im2.shape[1] // 2, :], vmin=0.0, vmax=BINS) 

837 elif sliceAxis == 2: 

838 plt.imshow(im2[:, :, im2.shape[2] // 2], vmin=0.0, vmax=BINS) 

839 plt.colorbar() 

840 

841 plt.subplot(223) 

842 tmp = p.copy() 

843 tmp[p <= 0] = numpy.nan 

844 tmp = numpy.log(tmp) 

845 plt.title("Log Probability log(p(im1,im2))") 

846 plt.imshow(tmp.T, origin="lower", extent=[0.0, BINS, 0.0, BINS])