Coverage for /usr/local/lib/python3.10/site-packages/spam/DIC/kinematics.py: 97%

1""" 

2Library of SPAM functions for post processing a deformation field. 

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 multiprocessing 

20 

21import numpy 

22import scipy.spatial 

23import spam.deformation 

24 

25try: 

26 multiprocessing.set_start_method("fork") 

27except RuntimeError: 

28 pass 

29import progressbar 

30 

31nProcessesDefault = multiprocessing.cpu_count() 

32 

33 

34def estimateLocalQuadraticCoherency(points, displacements, neighbourRadius=None, nNeighbours=None, epsilon=0.1, nProcesses=nProcessesDefault, verbose=False): 

35 """ 

36 This function computes the local quadratic coherency (LQC) of a set of displacement vectors as per Masullo and Theunissen 2016. 

37 LQC is the average residual between the point's displacement and a second-order (parabolic) surface Phi. 

38 The quadratic surface Phi is fitted to the point's closest N neighbours and evaluated at the point's position. 

39 Neighbours are selected based on: radius (default option) or number (activated if nNeighbours is not None). 

40 A point with a LQC value smaller than a threshold (0.1 in Masullo and Theunissen 2016) is classified as coherent 

41 

42 Parameters 

43 ---------- 

44 points : n x 3 numpy array of floats 

45 Coordinates of the points Z, Y, X 

46 

47 displacements : n x 3 numpy array of floats 

48 Displacements of the points 

49 

50 neighbourRadius: float, optional 

51 Distance in pixels around the point to extract neighbours. 

52 This OR nNeighbours must be set. 

53 Default = None 

54 

55 nNeighbours : int, optional 

56 Number of the nearest neighbours to consider 

57 This OR neighbourRadius must be set. 

58 Default = None 

59 

60 epsilon: float, optional 

61 Background error as per (Westerweel and Scarano 2005) 

62 Default = 0.1 

63 

64 nProcesses : integer, optional 

65 Number of processes for multiprocessing 

66 Default = number of CPUs in the system 

67 

68 verbose : boolean, optional 

69 Print progress? 

70 Default = True 

71 

72 Returns 

73 ------- 

74 LQC: n x 1 array of floats 

75 The local quadratic coherency for each point 

76 

77 Note 

78 ----- 

79 Based on: https://gricad-gitlab.univ-grenoble-alpes.fr/DVC/pt4d 

80 """ 

81 

82 # initialise the coherency matrix 

83 LQC = numpy.ones((points.shape[0]), dtype=float) 

84 

85 # build KD-tree for quick neighbour identification 

86 treeCoord = scipy.spatial.KDTree(points) 

87 

88 # check if neighbours are selected based on radius or based on number, default is radius 

89 if (nNeighbours is None) == (neighbourRadius is None): 

90 print("spam.DIC.estimateLocalQuadraticCoherency(): One and only one of nNeighbours and neighbourRadius must be passed") 

91 if nNeighbours is not None: 

92 ball = False 

93 elif neighbourRadius is not None: 

94 ball = True 

95 

96 # calculate coherency for each point 

97 global _multiprocessingCoherencyOnePoint 

98 

99 def _multiprocessingCoherencyOnePoint(point): 

100 # select neighbours based on radius 

101 if ball: 

102 radius = neighbourRadius 

103 indices = numpy.array(treeCoord.query_ball_point(points[point], radius)) 

104 # make sure that at least 27 neighbours are selected 

105 while len(indices) <= 27: 

106 radius *= 2 

107 indices = numpy.array(treeCoord.query_ball_point(points[point], radius)) 

108 N = len(indices) 

109 # select neighbours based on number 

110 else: 

111 _, indices = treeCoord.query(points[point], k=nNeighbours) 

112 N = nNeighbours 

113 

114 # fill in point+neighbours positions for the parabolic surface coefficients 

115 X = numpy.zeros((N, 10), dtype=float) 

116 for i, neighbour in enumerate(indices): 

117 pos = points[neighbour] 

118 X[i, 0] = 1 

119 X[i, 1] = pos[0] 

120 X[i, 2] = pos[1] 

121 X[i, 3] = pos[2] 

122 X[i, 4] = pos[0] * pos[1] 

123 X[i, 5] = pos[0] * pos[2] 

124 X[i, 6] = pos[1] * pos[2] 

125 X[i, 7] = pos[0] * pos[0] 

126 X[i, 8] = pos[1] * pos[1] 

127 X[i, 9] = pos[2] * pos[2] 

128 

129 # keep point's index 

130 i0 = numpy.where(indices == point)[0][0] 

131 

132 # fill in disp 

133 u = displacements[indices, 0] 

134 v = displacements[indices, 1] 

135 w = displacements[indices, 2] 

136 UnormMedian = numpy.median(numpy.linalg.norm(displacements[indices], axis=1)) 

137 

138 # deviation of each disp vector from local median 

139 sigma2 = (u - numpy.median(u)) ** 2 + (v - numpy.median(v)) ** 2 + (w - numpy.median(w)) ** 2 

140 

141 # coefficient for gaussian weighting 

142 K = (numpy.sqrt(sigma2).sum()) / N 

143 K += epsilon 

144 

145 # fill in gaussian weighting diag components 

146 Wg = numpy.exp(-0.5 * sigma2 * K ** (-0.5)) 

147 # create the diag matrix 

148 Wdiag = numpy.diag(Wg) 

149 

150 # create matrices to solve with least-squares 

151 XtWX = numpy.dot(X.T, numpy.dot(Wdiag, X)) 

152 XtWXInv = numpy.linalg.inv(XtWX) # TODO: check for singular matrix 

153 XtWUu = numpy.dot(X.T, numpy.dot(Wdiag, u)) 

154 XtWUv = numpy.dot(X.T, numpy.dot(Wdiag, v)) 

155 XtWUw = numpy.dot(X.T, numpy.dot(Wdiag, w)) 

156 

157 # solve least-squares to compute the coefficients of the parabolic surface 

158 au = numpy.dot(XtWXInv, XtWUu) 

159 av = numpy.dot(XtWXInv, XtWUv) 

160 aw = numpy.dot(XtWXInv, XtWUw) 

161 

162 # evaluate parabolic surface at point's position 

163 phiu = numpy.dot(au, X[i0, :]) 

164 phiv = numpy.dot(av, X[i0, :]) 

165 phiw = numpy.dot(aw, X[i0, :]) 

166 

167 # compute normalised residuals 

168 Cu = (phiu - u[i0]) ** 2 / (UnormMedian + epsilon) ** 2 

169 Cv = (phiv - v[i0]) ** 2 / (UnormMedian + epsilon) ** 2 

170 Cw = (phiw - w[i0]) ** 2 / (UnormMedian + epsilon) ** 2 

171 

172 # return coherency as the average normalised residual 

173 return point, (Cu + Cv + Cw) / 3 

174 

175 if verbose: 

176 pbar = progressbar.ProgressBar(maxval=points.shape[0]).start() 

177 finishedPoints = 0 

178 

179 with multiprocessing.Pool(processes=nProcesses) as pool: 

180 for returns in pool.imap_unordered(_multiprocessingCoherencyOnePoint, range(points.shape[0])): 

181 if verbose: 

182 finishedPoints += 1 

183 pbar.update(finishedPoints) 

184 LQC[returns[0]] = returns[1] 

185 pool.close() 

186 pool.join() 

187 

188 if verbose: 

189 pbar.finish() 

190 

191 return LQC 

192 

193 

194def estimateDisplacementFromQuadraticFit(fieldCoords, displacements, pointsToEstimate, neighbourRadius=None, nNeighbours=None, epsilon=0.1, nProcesses=nProcessesDefault, verbose=False): 

195 """ 

196 This function estimates the displacement of an incoherent point based on a local quadratic fit 

197 of the displacements of N coherent neighbours, as per Masullo and Theunissen 2016. 

198 A quadratic surface Phi is fitted to the point's closest coherent neighbours. 

199 Neighbours are selected based on: radius (default option) or number (activated if nNeighbours is not None) 

200 

201 Parameters 

202 ---------- 

203 fieldCoords : n x 3 numpy array of floats 

204 Coordinates of the points Z, Y, X where displacement is defined 

205 

206 displacements : n x 3 numpy array of floats 

207 Displacements of the points 

208 

209 pointsToEstimate : m x 3 numpy array of floats 

210 Coordinates of the points Z, Y, X where displacement should be estimated 

211 

212 neighbourRadius: float, optional 

213 Distance in pixels around the point to extract neighbours. 

214 This OR nNeighbours must be set. 

215 Default = None 

216 

217 nNeighbours : int, optional 

218 Number of the nearest neighbours to consider 

219 This OR neighbourRadius must be set. 

220 Default = None 

221 

222 epsilon: float, optional 

223 Background error as per (Westerweel and Scarano 2005) 

224 Default = 0.1 

225 

226 nProcesses : integer, optional 

227 Number of processes for multiprocessing 

228 Default = number of CPUs in the system 

229 

230 verbose : boolean, optional 

231 Print progress? 

232 Default = True 

233 

234 Returns 

235 ------- 

236 displacements: m x 3 array of floats 

237 The estimated displacements at the requested positions. 

238 

239 Note 

240 ----- 

241 Based on https://gricad-gitlab.univ-grenoble-alpes.fr/DVC/pt4d 

242 """ 

243 estimatedDisplacements = numpy.zeros_like(pointsToEstimate) 

244 

245 # build KD-tree of coherent points for quick neighbour identification 

246 treeCoord = scipy.spatial.KDTree(fieldCoords) 

247 

248 # check if neighbours are selected based on radius or based on number, default is radius 

249 if (nNeighbours is None) == (neighbourRadius is None): 

250 print("spam.DIC.estimateDisplacementFromQuadraticFit(): One and only one of nNeighbours and neighbourRadius must be passed") 

251 

252 ball = None 

253 if nNeighbours is not None: 

254 ball = False 

255 elif neighbourRadius is not None: 

256 ball = True 

257 

258 # estimate disp for each incoherent point 

259 global _multiprocessingDispOnePoint 

260 

261 def _multiprocessingDispOnePoint(pointToEstimate): 

262 # select neighbours based on radius 

263 if ball: 

264 radius = neighbourRadius 

265 indices = numpy.array(treeCoord.query_ball_point(pointsToEstimate[pointToEstimate], radius)) 

266 # make sure that at least 27 neighbours are selected 

267 while len(indices) <= 27: 

268 radius *= 2 

269 indices = numpy.array(treeCoord.query_ball_point(pointsToEstimate[pointToEstimate], radius)) 

270 N = len(indices) 

271 # select neighbours based on number 

272 else: 

273 _, indices = treeCoord.query(pointsToEstimate[pointToEstimate], k=nNeighbours) 

274 N = nNeighbours 

275 

276 # fill in neighbours positions for the parabolic surface coefficients 

277 X = numpy.zeros((N, 10), dtype=float) 

278 for i, neighbour in enumerate(indices): 

279 pos = fieldCoords[neighbour] 

280 X[i, 0] = 1 

281 X[i, 1] = pos[0] 

282 X[i, 2] = pos[1] 

283 X[i, 3] = pos[2] 

284 X[i, 4] = pos[0] * pos[1] 

285 X[i, 5] = pos[0] * pos[2] 

286 X[i, 6] = pos[1] * pos[2] 

287 X[i, 7] = pos[0] * pos[0] 

288 X[i, 8] = pos[1] * pos[1] 

289 X[i, 9] = pos[2] * pos[2] 

290 

291 # fill in point's position for the evaluation of the parabolic surface 

292 pos0 = pointsToEstimate[pointToEstimate] 

293 X0 = numpy.zeros((10), dtype=float) 

294 X0[0] = 1 

295 X0[1] = pos0[0] 

296 X0[2] = pos0[1] 

297 X0[3] = pos0[2] 

298 X0[4] = pos0[0] * pos0[1] 

299 X0[5] = pos0[0] * pos0[2] 

300 X0[6] = pos0[1] * pos0[2] 

301 X0[7] = pos0[0] * pos0[0] 

302 X0[8] = pos0[1] * pos0[1] 

303 X0[9] = pos0[2] * pos0[2] 

304 

305 # fill in disp of neighbours 

306 u = displacements[indices, 0] 

307 v = displacements[indices, 1] 

308 w = displacements[indices, 2] 

309 # UnormMedian = numpy.median(numpy.linalg.norm(displacements[indices], axis=1)) 

310 

311 # deviation of each disp vector from local median 

312 sigma2 = (u - numpy.median(u)) ** 2 + (v - numpy.median(v)) ** 2 + (w - numpy.median(w)) ** 2 

313 

314 # coefficient for gaussian weighting 

315 K = (numpy.sqrt(sigma2).sum()) / N 

316 K += epsilon 

317 

318 # fill in gaussian weighting diag components 

319 Wg = numpy.exp(-0.5 * sigma2 * K ** (-0.5)) # careful I think the first 0.5 was missing 

320 # create the diag matrix 

321 Wdiag = numpy.diag(Wg) 

322 

323 # create matrices to solve with least-squares 

324 XtWX = numpy.dot(X.T, numpy.dot(Wdiag, X)) 

325 XtWXInv = numpy.linalg.inv(XtWX) # TODO: check for singular matrix 

326 XtWUu = numpy.dot(X.T, numpy.dot(Wdiag, u)) 

327 XtWUv = numpy.dot(X.T, numpy.dot(Wdiag, v)) 

328 XtWUw = numpy.dot(X.T, numpy.dot(Wdiag, w)) 

329 

330 # solve least-squares to compute the coefficients of the parabolic surface 

331 au = numpy.dot(XtWXInv, XtWUu) 

332 av = numpy.dot(XtWXInv, XtWUv) 

333 aw = numpy.dot(XtWXInv, XtWUw) 

334 

335 # evaluate parabolic surface at incoherent point's position 

336 phiu = numpy.dot(au, X0) 

337 phiv = numpy.dot(av, X0) 

338 phiw = numpy.dot(aw, X0) 

339 

340 return pointToEstimate, [phiu, phiv, phiw] 

341 

342 # Iterate through flat field of Fs 

343 if verbose: 

344 pbar = progressbar.ProgressBar(maxval=pointsToEstimate.shape[0]).start() 

345 finishedPoints = 0 

346 

347 # Run multiprocessing filling in FfieldFlatGood, which will then update FfieldFlat 

348 with multiprocessing.Pool(processes=nProcesses) as pool: 

349 for returns in pool.imap_unordered(_multiprocessingDispOnePoint, range(pointsToEstimate.shape[0])): 

350 if verbose: 

351 finishedPoints += 1 

352 pbar.update(finishedPoints) 

353 estimatedDisplacements[returns[0]] = returns[1] 

354 pool.close() 

355 pool.join() 

356 

357 if verbose: 

358 pbar.finish() 

359 

360 # overwrite bad points displacements 

361 return estimatedDisplacements 

362 

363 

364def interpolatePhiField( 

365 fieldCoords, 

366 PhiField, 

367 pointsToInterpolate, 

368 nNeighbours=None, 

369 neighbourRadius=None, 

370 interpolateF="all", 

371 neighbourDistanceWeight="inverse", 

372 checkPointSurrounded=False, 

373 nProcesses=nProcessesDefault, 

374 verbose=False, 

375): 

376 """ 

377 This function interpolates components of a Phi field at a given number of points, using scipy's KD-tree to find neighbours. 

378 

379 Parameters 

380 ---------- 

381 fieldCoords : 2D array 

382 nx3 array of n points coordinates (ZYX) 

383 centre where each deformation function Phi has been measured 

384 

385 PhiField : 3D array 

386 nx4x4 array of n points deformation functions 

387 

388 pointsToInterpolate : 2D array 

389 mx3 array of m points coordinates (ZYX) 

390 Points where the deformation function Phi should be interpolated 

391 

392 nNeighbours : int, optional 

393 Number of the nearest neighbours to consider 

394 This OR neighbourRadius must be set. 

395 Default = None 

396 

397 neighbourRadius: float, optional 

398 Distance in pixels around the point to extract neighbours. 

399 This OR nNeighbours must be set. 

400 Default = None 

401 

402 interpolateF : string, optional 

403 Interpolate the whole Phi, just the rigid part, or just the displacement? 

404 Corresponding options are 'all', 'rigid', 'no' 

405 Default = "all" 

406 

407 neighbourDistanceWeight : string, optional 

408 How to weight neigbouring points? 

409 Possible approaches: inverse of distance, gaussian weighting, straight average, median 

410 Corresponding options: 'inverse', 'gaussian', 'mean', 'median' 

411 

412 checkPointSurrounded : bool, optional 

413 Only interpolate points whose neighbours surround them in Z, Y, X directions 

414 (or who fall exactly on a give point)? 

415 Default = False 

416 

417 nProcesses : integer, optional 

418 Number of processes for multiprocessing 

419 Default = number of CPUs in the system 

420 

421 verbose : bool, optional 

422 follow the progress of the function. 

423 Default = False 

424 

425 Returns 

426 ------- 

427 PhiField : mx4x4 array 

428 Interpolated **Phi** functions at the requested positions 

429 """ 

430 tol = 1e-6 # OS is responsible for the validitidy of this magic number 

431 

432 numberOfPointsToInterpolate = pointsToInterpolate.shape[0] 

433 # create the k-d tree of the coordinates of good points, we need this to search for the k nearest neighbours easily 

434 # for details see: https://en.wikipedia.org/wiki/K-d_tree & 

435 # https://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.spatial.KDTree.query.html 

436 treeCoord = scipy.spatial.KDTree(fieldCoords) 

437 

438 # extract the Phi matrices of the bad points 

439 outputPhiField = numpy.zeros((numberOfPointsToInterpolate, 4, 4), dtype=PhiField.dtype) 

440 

441 assert interpolateF in ["all", "rigid", "no"], "spam.DIC.interpolatePhiField(): interpolateF argument should either be 'all', 'rigid', or 'no'" 

442 assert neighbourDistanceWeight in ["inverse", "gaussian", "mean", "median"], "spam.DIC.interpolatePhiField(): neighbourDistanceWeight argument should be 'inverse', 'gaussian', 'mean', or 'median'" 

443 # check if neighbours are selected based on radius or based on number, default is radius 

444 assert (nNeighbours is None) != (neighbourRadius is None), "spam.DIC.interpolatePhiField(): One and only one of nNeighbours and neighbourRadius must be passed" 

445 

446 if nNeighbours is not None: 

447 ball = False 

448 elif neighbourRadius is not None: 

449 ball = True 

450 

451 global _multiprocessingInterpolateOnePoint 

452 

453 def _multiprocessingInterpolateOnePoint(pointNumber): 

454 pointToInterpolate = pointsToInterpolate[pointNumber] 

455 outputPhi = numpy.zeros((4, 4), dtype=PhiField.dtype) 

456 outputPhi[-1] = [0, 0, 0, 1] 

457 if interpolateF == "no": 

458 outputPhi[0:3, 0:3] = numpy.eye(3) 

459 

460 ####################################################### 

461 # Find neighbours 

462 ####################################################### 

463 if ball: 

464 # Ball lookup 

465 indices = treeCoord.query_ball_point(pointToInterpolate, neighbourRadius) 

466 if len(indices) == 0: 

467 # No point! 

468 return pointNumber, numpy.eye(4) * numpy.nan 

469 else: 

470 distances = numpy.linalg.norm(pointToInterpolate - fieldCoords[indices], axis=1) 

471 else: 

472 # Number of Neighbour lookup 

473 distances, indices = treeCoord.query(pointToInterpolate, k=nNeighbours) 

474 indices = numpy.array(indices) 

475 distances = numpy.array(distances) 

476 

477 ####################################################### 

478 # Check if there is only one neighbour 

479 ####################################################### 

480 if indices.size == 1: 

481 if checkPointSurrounded: 

482 # unless they're the same point, can't be surrounded 

483 if not numpy.allclose(fieldCoords[indices], pointToInterpolate): 

484 return pointNumber, numpy.eye(4) * numpy.nan 

485 

486 if interpolateF in ["all", "rigid"]: # We need to interpolate all 12 components of Phi 

487 outputPhi = PhiField[indices].copy() 

488 if interpolateF == "rigid": 

489 outputPhi = spam.deformation.computeRigidPhi(outputPhi) 

490 else: # interpolate only displacements 

491 outputPhi[0:3, -1] = PhiField[indices, 0:3, -1].copy() 

492 

493 return pointNumber, outputPhi 

494 

495 ####################################################### 

496 # If > 1 neighbour, interpolate Phi or displacements 

497 ####################################################### 

498 else: 

499 if neighbourDistanceWeight == "inverse": 

500 weights = 1 / (distances + tol) 

501 elif neighbourDistanceWeight == "gaussian": 

502 # This could be an input variable VVVVVVVVVVVVVVVVVVVVVV--- the gaussian weighting distance 

503 weights = numpy.exp(-(distances**2) / numpy.max(distances / 2) ** 2) 

504 elif neighbourDistanceWeight == "mean": 

505 weights = numpy.ones_like(distances) 

506 elif neighbourDistanceWeight == "median": 

507 # is this the equivalent kernel to a median, we think so... 

508 weights = numpy.zeros_like(distances) 

509 weights[len(distances) // 2] = 1 

510 

511 if checkPointSurrounded: 

512 posMax = numpy.array([fieldCoords[indices, i].max() for i in range(3)]) 

513 posMin = numpy.array([fieldCoords[indices, i].min() for i in range(3)]) 

514 if not numpy.logical_and(numpy.all(pointToInterpolate >= posMin), numpy.all(pointToInterpolate <= posMax)): 

515 return pointNumber, numpy.eye(4) * numpy.nan 

516 

517 if interpolateF == "no": 

518 outputPhi[0:3, -1] = numpy.sum(PhiField[indices, 0:3, -1] * weights[:, numpy.newaxis], axis=0) / weights.sum() 

519 else: 

520 outputPhi[:-1] = numpy.sum(PhiField[indices, :-1] * weights[:, numpy.newaxis, numpy.newaxis], axis=0) / weights.sum() 

521 

522 if interpolateF == "rigid": 

523 outputPhi = spam.deformation.computeRigidPhi(outputPhi) 

524 

525 return pointNumber, outputPhi 

526 

527 if verbose: 

528 print("\nStarting Phi field interpolation (with {} process{})".format(nProcesses, "es" if nProcesses > 1 else "")) 

529 widgets = [progressbar.Bar(), " ", progressbar.AdaptiveETA()] 

530 pbar = progressbar.ProgressBar(widgets=widgets, maxval=numberOfPointsToInterpolate) 

531 pbar.start() 

532 finishedNodes = 0 

533 

534 with multiprocessing.Pool(processes=nProcesses) as pool: 

535 for returns in pool.imap_unordered(_multiprocessingInterpolateOnePoint, range(numberOfPointsToInterpolate)): 

536 # Update progres bar if point is not skipped 

537 if verbose: 

538 pbar.update(finishedNodes) 

539 finishedNodes += 1 

540 

541 outputPhiField[returns[0]] = returns[1] 

542 pool.close() 

543 pool.join() 

544 

545 if verbose: 

546 pbar.finish() 

547 

548 return outputPhiField 

549 

550 

551def mergeRegularGridAndDiscrete(regularGrid=None, discrete=None, labelledImage=None, binningLabelled=1, alwaysLabel=True): 

552 """ 

553 This function corrects a Phi field from the spam-ldic script (measured on a regular grid) 

554 by looking into the results from one or more spam-ddic results (measured on individual labels) 

555 by applying the discrete measurements to the grid points. 

556 

557 This can be useful where there are large flat zones in the image that cannot 

558 be correlated with small correlation windows, but can be identified and 

559 tracked with a spam-ddic computation (concrete is a good example). 

560 

561 Parameters 

562 ----------- 

563 regularGrid : dictionary 

564 Dictionary containing readCorrelationTSV of regular grid correlation script, `spam-ldic`. 

565 Default = None 

566 

567 discrete : dictionary or list of dictionaries 

568 Dictionary (or list thereof) containing readCorrelationTSV of discrete correlation script, `spam-ddic`. 

569 File name of TSV from DICdiscrete client, or list of filenames 

570 Default = None 

571 

572 labelledImage : 3D numpy array of ints, or list of numpy arrays 

573 Labelled volume used for discrete computation 

574 Default = None 

575 

576 binningLabelled : int 

577 Are the labelled images and their PhiField at a different bin level than 

578 the regular field? 

579 Default = 1 

580 

581 alwaysLabel : bool 

582 If regularGrid point falls inside the label, should we use the 

583 label displacement automatically? 

584 Otherwise if the regularGrid point has converged should we use that? 

585 Default = True (always use Label displacement) 

586 

587 Returns 

588 -------- 

589 Either dictionary or TSV file 

590 Output matrix, with number of rows equal to spam-ldic (the node spacing of the regular grid) and with columns: 

591 "NodeNumber", "Zpos", "Ypos", "Xpos", "Zdisp", "Ydisp", "Xdisp", "deltaPhiNorm", "returnStatus", "iterations" 

592 """ 

593 import spam.helpers 

594 

595 # If we have a list of input discrete files, we also need a list of labelled images 

596 if type(discrete) == list: 

597 if type(labelledImage) != list: 

598 print("spam.deformation.deformationFunction.mergeRegularGridAndDiscrete(): if you pass a list of discreteTSV you must also pass a list of labelled images") 

599 return 

600 if len(discrete) != len(labelledImage): 

601 print("spam.deformation.deformationFunction.mergeRegularGridAndDiscrete(): len(discrete) must be equal to len(labelledImage)") 

602 return 

603 nDiscrete = len(discrete) 

604 

605 # We have only one TSV and labelled image, it should be a number array 

606 else: 

607 if type(labelledImage) != numpy.ndarray: 

608 print("spam.deformation.deformationFunction.mergeRegularGridAndDiscrete(): with a single discrete TSV file, labelledImage must be a numpy array") 

609 return 

610 discrete = [discrete] 

611 labelledImage = [labelledImage] 

612 nDiscrete = 1 

613 

614 output = {} 

615 

616 # Regular grid is the master, and so we copy dimensions and positions 

617 output["fieldDims"] = regularGrid["fieldDims"] 

618 output["fieldCoords"] = regularGrid["fieldCoords"] 

619 

620 output["PhiField"] = numpy.zeros_like(regularGrid["PhiField"]) 

621 output["iterations"] = numpy.zeros_like(regularGrid["iterations"]) 

622 output["deltaPhiNorm"] = numpy.zeros_like(regularGrid["deltaPhiNorm"]) 

623 output["returnStatus"] = numpy.zeros_like(regularGrid["returnStatus"]) 

624 output["pixelSearchCC"] = numpy.zeros_like(regularGrid["returnStatus"]) 

625 output["error"] = numpy.zeros_like(regularGrid["returnStatus"]) 

626 output["mergeSource"] = numpy.zeros_like(regularGrid["iterations"]) 

627 

628 # from progressbar import ProgressBar 

629 # pbar = ProgressBar() 

630 

631 # For each point on the regular grid... 

632 # for n, gridPoint in pbar(enumerate(regularGrid['fieldCoords'].astype(int))): 

633 for n, gridPoint in enumerate(regularGrid["fieldCoords"].astype(int)): 

634 # Find labels corresponding to this grid position for the labelledImage images 

635 labels = [] 

636 for m in range(nDiscrete): 

637 labels.append(int(labelledImage[m][int(gridPoint[0] / float(binningLabelled)), int(gridPoint[1] / float(binningLabelled)), int(gridPoint[2] / float(binningLabelled))])) 

638 labels = numpy.array(labels) 

639 

640 # Is the point inside a discrete label? 

641 if (labels == 0).all() or (not alwaysLabel and regularGrid["returnStatus"][n] == 2): 

642 # Use the REGULAR GRID MEASUREMENT 

643 # If we're not in a label, copy the results from DICregularGrid 

644 output["PhiField"][n] = regularGrid["PhiField"][n] 

645 output["deltaPhiNorm"][n] = regularGrid["deltaPhiNorm"][n] 

646 output["returnStatus"][n] = regularGrid["returnStatus"][n] 

647 output["iterations"][n] = regularGrid["iterations"][n] 

648 # output['error'][n] = regularGrid['error'][n] 

649 # output['pixelSearchCC'][n] = regularGrid['pixelSearchCC'][n] 

650 else: 

651 # Use the DISCRETE MEASUREMENT 

652 # Give precedence to earliest non-zero-labelled discrete field, conflicts not handled 

653 m = numpy.where(labels != 0)[0][0] 

654 label = labels[m] 

655 # print("m,label = ", m, label) 

656 tmp = discrete[m]["PhiField"][label].copy() 

657 tmp[0:3, -1] *= float(binningLabelled) 

658 translation = spam.deformation.decomposePhi(tmp, PhiCentre=discrete[m]["fieldCoords"][label] * float(binningLabelled), PhiPoint=gridPoint)["t"] 

659 # This is the Phi we will save for this point -- take the F part of the labelled's Phi 

660 phi = discrete[m]["PhiField"][label].copy() 

661 # ...and add the computed displacement as applied to the grid point 

662 phi[0:3, -1] = translation 

663 

664 output["PhiField"][n] = phi 

665 output["deltaPhiNorm"][n] = discrete[m]["deltaPhiNorm"][label] 

666 output["returnStatus"][n] = discrete[m]["returnStatus"][label] 

667 output["iterations"][n] = discrete[m]["iterations"][label] 

668 # output['error'][n] = discrete[m]['error'][label] 

669 # output['pixelSearchCC'][n] = discrete[m]['pixelSearchCC'][label] 

670 output["mergeSource"][n] = m + 1 

671 

672 # if fileName is not None: 

673 # TSVheader = "NodeNumber\tZpos\tYpos\tXpos\tFzz\tFzy\tFzx\tZdisp\tFyz\tFyy\tFyx\tYdisp\tFxz\tFxy\tFxx\tXdisp" 

674 # outMatrix = numpy.array([numpy.array(range(output['fieldCoords'].shape[0])), 

675 # output['fieldCoords'][:, 0], output['fieldCoords'][:, 1], output['fieldCoords'][:, 2], 

676 # output['PhiField'][:, 0, 0], output['PhiField'][:, 0, 1], output['PhiField'][:, 0, 2], output['PhiField'][:, 0, 3], 

677 # output['PhiField'][:, 1, 0], output['PhiField'][:, 1, 1], output['PhiField'][:, 1, 2], output['PhiField'][:, 1, 3], 

678 # output['PhiField'][:, 2, 0], output['PhiField'][:, 2, 1], output['PhiField'][:, 2, 2], output['PhiField'][:, 2, 3]]).T 

679 

680 # outMatrix = numpy.hstack([outMatrix, numpy.array([output['iterations'], 

681 # output['returnStatus'], 

682 # output['deltaPhiNorm'], 

683 # output['mergeSource']]).T]) 

684 # TSVheader = TSVheader+"\titerations\treturnStatus\tdeltaPhiNorm\tmergeSource" 

685 

686 # numpy.savetxt(fileName, 

687 # outMatrix, 

688 # fmt='%.7f', 

689 # delimiter='\t', 

690 # newline='\n', 

691 # comments='', 

692 # header=TSVheader) 

693 # else: 

694 return output 

695 

696 

697def getDisplacementFromNeighbours(labIm, DVC, fileName, method="getLabel", neighboursRange=None, centresOfMass=None, previousDVC=None): 

698 """ 

699 This function computes the displacement as the mean displacement from the neighbours, 

700 for non-converged grains using a TSV file obtained from `spam-ddic` script. 

701 Returns a new TSV file with the new Phi (composed only by the displacement part). 

702 

703 The generated TSV can be used as an input for `spam-ddic`. 

704 

705 Parameters 

706 ----------- 

707 lab : 3D array of integers 

708 Labelled volume, with lab.max() labels 

709 

710 DVC : dictionary 

711 Dictionary with deformation field, obtained from `spam-ddic` script, 

712 and read using `spam.helpers.tsvio.readCorrelationTSV()` with `readConvergence=True, readPSCC=True, readLabelDilate=True` 

713 

714 fileName : string 

715 FileName including full path and .tsv at the end to write 

716 

717 method : string, optional 

718 Method to compute the neighbours using `spam.label.getNeighbours()`. 

719 'getLabel' : The neighbours are the labels inside the subset obtained through spam.getLabel() 

720 'mesh' : The neighbours are computed using a tetrahedral connectivity matrix 

721 Default = 'getLabel' 

722 

723 neighboursRange : int 

724 Parameter controlling the search range to detect neighbours for each method. 

725 For 'getLabel', it correspond to the size of the subset. Default = meanRadii 

726 For 'mesh', it correspond to the size of the alpha shape used for carving the mesh. Default = 5*meanDiameter. 

727 

728 centresOfMass : lab.max()x3 array of floats, optional 

729 Centres of mass in format returned by ``centresOfMass``. 

730 If not defined (Default = None), it is recomputed by running ``centresOfMass`` 

731 

732 previousDVC : dictionary, optional 

733 Dictionary with deformation field, obtained from `spam-ddic` script, and read using `spam.helpers.tsvio.readCorrelationTSV()` for the previous step. 

734 This allows the to compute only the displacement increment from the neighbours, while using the F tensor from a previous (converged) step. 

735 If `previousDVS = None`, then the resulting Phi would be composed only by the displacement of the neighbours. 

736 Default = None 

737 

738 Returns 

739 -------- 

740 Dictionary 

741 TSV file with the same columns as the input 

742 """ 

743 import spam.label 

744 

745 # Compute centreOfMass if needed 

746 if centresOfMass is None: 

747 centresOfMass = spam.label.centresOfMass(labIm) 

748 # Get number of labels 

749 numberOfLabels = (labIm.max() + 1).astype("u4") 

750 # Create Phi field 

751 PhiField = numpy.zeros((numberOfLabels, 4, 4), dtype="<f4") 

752 # Rest of arrays 

753 try: 

754 iterations = DVC["iterations"] 

755 returnStatus = DVC["returnStatus"] 

756 deltaPhiNorm = DVC["deltaPhiNorm"] 

757 PSCC = DVC["pixelSearchCC"] 

758 labelDilateList = DVC["LabelDilate"] 

759 error = DVC["error"] 

760 

761 # Get the problematic labels 

762 probLab = numpy.where(DVC["returnStatus"] != 2)[0] 

763 # Remove the 0 from the wrongLab list 

764 probLab = probLab[~numpy.isin(probLab, 0)] 

765 # Get neighbours 

766 neighbours = spam.label.getNeighbours(labIm, probLab, method=method, neighboursRange=neighboursRange) 

767 # Solve first the converged particles - make a copy of the PhiField 

768 for i in range(numberOfLabels): 

769 PhiField[i] = DVC["PhiField"][i] 

770 # Work on the problematic labels 

771 widgets = [progressbar.FormatLabel(""), " ", progressbar.Bar(), " ", progressbar.AdaptiveETA()] 

772 pbar = progressbar.ProgressBar(widgets=widgets, maxval=len(probLab)) 

773 pbar.start() 

774 for i in range(0, len(probLab), 1): 

775 wrongLab = probLab[i] 

776 neighboursLabel = neighbours[i] 

777 t = [] 

778 for n in neighboursLabel: 

779 # Check the return status of each neighbour 

780 if DVC["returnStatus"][n] == 2: 

781 # We dont have a previous DVC file loaded 

782 if previousDVC is None: 

783 # Get translation, rotation and zoom from Phi 

784 nPhi = spam.deformation.deformationFunction.decomposePhi(DVC["PhiField"][n]) 

785 # Append the results 

786 t.append(nPhi["t"]) 

787 # We have a previous DVC file loaded 

788 else: 

789 # Get translation, rotation and zoom from Phi at t=step 

790 nPhi = spam.deformation.deformationFunction.decomposePhi(DVC["PhiField"][n]) 

791 # Get translation, rotation and zoom from Phi at t=step-1 

792 nPhiP = spam.deformation.deformationFunction.decomposePhi(previousDVC["PhiField"][n]) 

793 # Append the incremental results 

794 t.append(nPhi["t"] - nPhiP["t"]) 

795 # Transform list to array 

796 if not t: 

797 # This is a non-working label, take care of it 

798 Phi = spam.deformation.computePhi({"t": [0, 0, 0]}) 

799 PhiField[wrongLab] = Phi 

800 else: 

801 t = numpy.asarray(t) 

802 # Compute mean 

803 meanT = numpy.mean(t, axis=0) 

804 # Reconstruct 

805 transformation = {"t": meanT} 

806 Phi = spam.deformation.computePhi(transformation) 

807 # Save 

808 if previousDVC is None: 

809 PhiField[wrongLab] = Phi 

810 else: 

811 PhiField[wrongLab] = previousDVC["PhiField"][wrongLab] 

812 # Add the incremental displacement 

813 PhiField[wrongLab][:-1, -1] += Phi[:-1, -1] 

814 # Update the progressbar