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

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 spam.deformation 

20import numpy 

21import tifffile 

22import scipy.spatial 

23import multiprocessing 

24import progressbar 

25 

26nProcessesDefault = multiprocessing.cpu_count() 

27 

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

29 ''' 

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

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

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

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

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

35 

36 Parameters 

37 ---------- 

38 points : n x 3 numpy array of floats 

39 Coordinates of the points Z, Y, X 

40 

41 displacements : n x 3 numpy array of floats 

42 Displacements of the points 

43 

44 neighbourRadius: float, optional 

45 Distance in pixels around the point to extract neighbours. 

46 This OR nNeighbours must be set. 

47 Default = None 

48 

49 nNeighbours : int, optional 

50 Number of the nearest neighbours to consider 

51 This OR neighbourRadius must be set. 

52 Default = None 

53 

54 epsilon: float, optional 

55 Background error as per (Westerweel and Scarano 2005) 

56 Default = 0.1 

57 

58 nProcesses : integer, optional 

59 Number of processes for multiprocessing 

60 Default = number of CPUs in the system 

61 

62 verbose : boolean, optional 

63 Print progress? 

64 Default = True 

65 

66 Returns 

67 ------- 

68 LQC: n x 1 array of floats 

69 The local quadratic coherency for each point 

70 

71 Note 

72 ----- 

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

74 ''' 

75 

76 # initialise the coherency matrix 

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

78 

79 # build KD-tree for quick neighbour identification 

80 treeCoord = scipy.spatial.KDTree(points) 

81 

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

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

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

85 if nNeighbours is not None: 

86 ball = False 

87 elif neighbourRadius is not None: 

88 ball = True 

89 

90 # calculate coherency for each point 

91 global coherencyOnePoint 

92 def coherencyOnePoint(point): 

93 # select neighbours based on radius 

94 if ball: 

95 radius = neighbourRadius 

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

97 # make sure that at least 27 neighbours are selected 

98 while len(indices) <= 27: 

99 radius *= 2 

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

101 N = len(indices) 

102 # select neighbours based on number 

103 else: 

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

105 N = nNeighbours 

106 

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

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

109 for i, neighbour in enumerate(indices): 

110 pos = points[neighbour] 

111 X[i, 0] = 1 

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

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

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

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

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

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

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

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

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

121 

122 # keep point's index 

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

124 

125 # fill in disp 

126 u = displacements[indices, 0] 

127 v = displacements[indices, 1] 

128 w = displacements[indices, 2] 

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

130 

131 # deviation of each disp vector from local median 

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

133 

134 # coefficient for gaussian weighting 

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

136 K += epsilon 

137 

138 # fill in gaussian weighting diag components 

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

140 # create the diag matrix 

141 Wdiag = numpy.diag(Wg) 

142 

143 # create matrices to solve with least-squares 

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

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

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

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

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

149 

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

151 au = numpy.dot(XtWXInv, XtWUu) 

152 av = numpy.dot(XtWXInv, XtWUv) 

153 aw = numpy.dot(XtWXInv, XtWUw) 

154 

155 # evaluate parabolic surface at point's position 

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

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

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

159 

160 # compute normalised residuals 

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

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

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

164 

165 # return coherency as the average normalised residual 

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

167 

168 if verbose: 

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

170 finishedPoints = 0 

171 

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

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

174 if verbose: 

175 finishedPoints += 1 

176 pbar.update(finishedPoints) 

177 LQC[returns[0]] = returns[1] 

178 pool.close() 

179 pool.join() 

180 

181 if verbose: pbar.finish() 

182 

183 return LQC 

184 

185 

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

187 ''' 

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

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

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

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

192 

193 Parameters 

194 ---------- 

195 fieldCoords : n x 3 numpy array of floats 

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

197 

198 displacements : n x 3 numpy array of floats 

199 Displacements of the points 

200 

201 pointsToEstimate : m x 3 numpy array of floats 

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

203 

204 neighbourRadius: float, optional 

205 Distance in pixels around the point to extract neighbours. 

206 This OR nNeighbours must be set. 

207 Default = None 

208 

209 nNeighbours : int, optional 

210 Number of the nearest neighbours to consider 

211 This OR neighbourRadius must be set. 

212 Default = None 

213 

214 epsilon: float, optional 

215 Background error as per (Westerweel and Scarano 2005) 

216 Default = 0.1 

217 

218 nProcesses : integer, optional 

219 Number of processes for multiprocessing 

220 Default = number of CPUs in the system 

221 

222 verbose : boolean, optional 

223 Print progress? 

224 Default = True 

225 

226 Returns 

227 ------- 

228 displacements: m x 3 array of floats 

229 The estimated displacements at the requested positions. 

230 

231 Note 

232 ----- 

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

234 ''' 

235 estimatedDisplacements = numpy.zeros_like(pointsToEstimate) 

236 

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

238 treeCoord = scipy.spatial.KDTree(fieldCoords) 

239 

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

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

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

243 

244 ball = None 

245 if nNeighbours is not None: 

246 ball = False 

247 elif neighbourRadius is not None: 

248 ball = True 

249 

250 # estimate disp for each incoherent point 

251 global dispOnePoint 

252 def dispOnePoint(pointToEstimate): 

253 # select neighbours based on radius 

254 if ball: 

255 radius = neighbourRadius 

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

257 # make sure that at least 27 neighbours are selected 

258 while len(indices) <= 27: 

259 radius *= 2 

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

261 N = len(indices) 

262 # select neighbours based on number 

263 else: 

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

265 N = nNeighbours 

266 

267 # fill in neighbours positions for the parabolic surface coefficients 

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

269 for i, neighbour in enumerate(indices): 

270 pos = fieldCoords[neighbour] 

271 X[i, 0] = 1 

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

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

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

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

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

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

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

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

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

281 

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

283 pos0 = pointsToEstimate[pointToEstimate] 

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

285 X0[0] = 1 

286 X0[1] = pos0[0] 

287 X0[2] = pos0[1] 

288 X0[3] = pos0[2] 

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

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

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

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

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

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

295 

296 # fill in disp of neighbours 

297 u = displacements[indices, 0] 

298 v = displacements[indices, 1] 

299 w = displacements[indices, 2] 

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

301 

302 # deviation of each disp vector from local median 

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

304 

305 # coefficient for gaussian weighting 

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

307 K += epsilon 

308 

309 # fill in gaussian weighting diag components 

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

311 # create the diag matrix 

312 Wdiag = numpy.diag(Wg) 

313 

314 # create matrices to solve with least-squares 

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

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

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

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

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

320 

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

322 au = numpy.dot(XtWXInv, XtWUu) 

323 av = numpy.dot(XtWXInv, XtWUv) 

324 aw = numpy.dot(XtWXInv, XtWUw) 

325 

326 # evaluate parabolic surface at incoherent point's position 

327 phiu = numpy.dot(au, X0) 

328 phiv = numpy.dot(av, X0) 

329 phiw = numpy.dot(aw, X0) 

330 

331 return pointToEstimate, [phiu, phiv, phiw] 

332 

333 # Iterate through flat field of Fs 

334 if verbose: 

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

336 finishedPoints = 0 

337 

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

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

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

341 if verbose: 

342 finishedPoints += 1 

343 pbar.update(finishedPoints) 

344 estimatedDisplacements[returns[0]] = returns[1] 

345 pool.close() 

346 pool.join() 

347 

348 if verbose: pbar.finish() 

349 

350 # overwrite bad points displacements 

351 return estimatedDisplacements 

352 

353 

354def interpolatePhiField(fieldCoords, PhiField, pointsToInterpolate, nNeighbours=None, neighbourRadius=None, interpolateF="all", neighbourDistanceWeight='inverse', checkPointSurrounded=False, nProcesses=nProcessesDefault, verbose=False): 

355 """ 

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

357 

358 Parameters 

359 ---------- 

360 fieldCoords : 2D array 

361 nx3 array of n points coordinates (ZYX) 

362 centre where each deformation function Phi has been measured 

363 

364 PhiField : 3D array 

365 nx4x4 array of n points deformation functions 

366 

367 pointsToInterpolate : 2D array 

368 mx3 array of m points coordinates (ZYX) 

369 Points where the deformation function Phi should be interpolated 

370 

371 nNeighbours : int, optional 

372 Number of the nearest neighbours to consider 

373 This OR neighbourRadius must be set. 

374 Default = None 

375 

376 neighbourRadius: float, optional 

377 Distance in pixels around the point to extract neighbours. 

378 This OR nNeighbours must be set. 

379 Default = None 

380 

381 interpolateF : string, optional 

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

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

384 Default = "all" 

385 

386 neighbourDistanceWeight : string, optional 

387 How to weight neigbouring points? 

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

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

390 

391 checkPointSurrounded : bool, optional 

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

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

394 Default = False 

395 

396 nProcesses : integer, optional 

397 Number of processes for multiprocessing 

398 Default = number of CPUs in the system 

399 

400 verbose : bool, optional 

401 follow the progress of the function. 

402 Default = False 

403 

404 Returns 

405 ------- 

406 PhiField : mx4x4 array 

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

408 """ 

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

410 

411 numberOfPointsToInterpolate = pointsToInterpolate.shape[0] 

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

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

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

415 treeCoord = scipy.spatial.KDTree(fieldCoords) 

416 

417 # extract the Phi matrices of the bad points 

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

419 

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

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

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

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

424 

425 if nNeighbours is not None: 

426 ball = False 

427 elif neighbourRadius is not None: 

428 ball = True 

429 

430 global interpolateOnePoint 

431 def interpolateOnePoint(pointNumber): 

432 pointToInterpolate = pointsToInterpolate[pointNumber] 

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

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

435 if interpolateF == 'no': 

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

437 

438 ####################################################### 

439 ### Find neighbours 

440 ####################################################### 

441 if ball: 

442 # Ball lookup 

443 indices = treeCoord.query_ball_point(pointToInterpolate, neighbourRadius) 

444 if len(indices) == 0: 

445 # No point! 

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

447 else: 

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

449 else: 

450 # Number of Neighbour lookup 

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

452 indices = numpy.array(indices) 

453 distances = numpy.array(distances) 

454 

455 ####################################################### 

456 ### Check if there is only one neighbour 

457 ####################################################### 

458 if indices.size == 1: 

459 if checkPointSurrounded: 

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

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

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

463 

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

465 outputPhi = PhiField[indices].copy() 

466 if interpolateF == "rigid": 

467 outputPhi = spam.deformation.computeRigidPhi(outputPhi) 

468 else: # interpolate only displacements 

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

470 

471 return pointNumber, outputPhi 

472 

473 ####################################################### 

474 ### If > 1 neighbour, interpolate Phi or displacements 

475 ####################################################### 

476 else: 

477 if neighbourDistanceWeight == 'inverse': 

478 weights = (1/(distances+tol)) 

479 elif neighbourDistanceWeight == 'gaussian': 

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

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

482 elif neighbourDistanceWeight == 'mean': 

483 weights = numpy.ones_like(distances) 

484 elif neighbourDistanceWeight == 'median': 

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

486 weights = numpy.zeros_like(distances) 

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

488 

489 if checkPointSurrounded: 

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

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

492 if not numpy.logical_and(numpy.all(pointToInterpolate >= posMin), 

493 numpy.all(pointToInterpolate <= posMax)): 

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

495 

496 if interpolateF == 'no': 

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

498 else: 

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

500 

501 if interpolateF == "rigid": 

502 outputPhi = spam.deformation.computeRigidPhi(outputPhi) 

503 

504 return pointNumber, outputPhi 

505 

506 

507 if verbose: 

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

509 widgets = [progressbar.Bar(), ' ', progressbar.AdaptiveETA()] 

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

511 pbar.start() 

512 finishedNodes = 0 

513 

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

515 for returns in pool.imap_unordered(interpolateOnePoint, range(numberOfPointsToInterpolate)): 

516 # Update progres bar if point is not skipped 

517 if verbose: 

518 pbar.update(finishedNodes) 

519 finishedNodes += 1 

520 

521 outputPhiField[returns[0]] = returns[1] 

522 pool.close() 

523 pool.join() 

524 

525 if verbose: pbar.finish() 

526 

527 return outputPhiField 

528 

529 

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

531 """ 

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

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

534 by applying the discrete measurements to the grid points. 

535 

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

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

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

539 

540 Parameters 

541 ----------- 

542 regularGrid : dictionary 

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

544 Default = None 

545 

546 discrete : dictionary or list of dictionaries 

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

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

549 Default = None 

550 

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

552 Labelled volume used for discrete computation 

553 Default = None 

554 

555 binningLabelled : int 

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

557 the regular field? 

558 Default = 1 

559 

560 alwaysLabel : bool 

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

562 label displacement automatically? 

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

564 Default = True (always use Label displacement) 

565 

566 Returns 

567 -------- 

568 Either dictionary or TSV file 

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

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

571 """ 

572 import spam.helpers 

573 

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

575 if type(discrete) == list: 

576 if type(labelledImage) != list: 

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

578 return 

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

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

581 return 

582 nDiscrete = len(discrete) 

583 

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

585 else: 

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

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

588 return 

589 discrete = [discrete] 

590 labelledImage = [labelledImage] 

591 nDiscrete = 1 

592 

593 output = {} 

594 

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

596 output['fieldDims'] = regularGrid['fieldDims'] 

597 output['fieldCoords'] = regularGrid['fieldCoords'] 

598 

599 output['PhiField'] = numpy.zeros_like(regularGrid['PhiField']) 

600 output['iterations'] = numpy.zeros_like(regularGrid['iterations']) 

601 output['deltaPhiNorm'] = numpy.zeros_like(regularGrid['deltaPhiNorm']) 

602 output['returnStatus'] = numpy.zeros_like(regularGrid['returnStatus']) 

603 output['pixelSearchCC'] = numpy.zeros_like(regularGrid['returnStatus']) 

604 output['error'] = numpy.zeros_like(regularGrid['returnStatus']) 

605 output['mergeSource'] = numpy.zeros_like(regularGrid['iterations']) 

606 

607 # from progressbar import ProgressBar 

608 # pbar = ProgressBar() 

609 

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

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

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

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

614 labels = [] 

615 for m in range(nDiscrete): 

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

617 int(gridPoint[1] / float(binningLabelled)), 

618 int(gridPoint[2] / float(binningLabelled))])) 

619 labels = numpy.array(labels) 

620 

621 # Is the point inside a discrete label? 

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

623 ### Use the REGULAR GRID MEASUREMENT 

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

625 output['PhiField'][n] = regularGrid['PhiField'][n] 

626 output['deltaPhiNorm'][n] = regularGrid['deltaPhiNorm'][n] 

627 output['returnStatus'][n] = regularGrid['returnStatus'][n] 

628 output['iterations'][n] = regularGrid['iterations'][n] 

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

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

631 else: 

632 ### Use the DISCRETE MEASUREMENT 

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

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

635 label = labels[m] 

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

637 tmp = discrete[m]['PhiField'][label].copy() 

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

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

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

641 phi = discrete[m]['PhiField'][label].copy() 

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

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

644 

645 output['PhiField'][n] = phi 

646 output['deltaPhiNorm'][n] = discrete[m]['deltaPhiNorm'][label] 

647 output['returnStatus'][n] = discrete[m]['returnStatus'][label] 

648 output['iterations'][n] = discrete[m]['iterations'][label] 

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

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

651 output['mergeSource'][n] = m+1 

652 

653 #if fileName is not None: 

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

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

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

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

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

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

660 

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

662 #output['returnStatus'], 

663 #output['deltaPhiNorm'], 

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

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

666 

667 #numpy.savetxt(fileName, 

668 #outMatrix, 

669 #fmt='%.7f', 

670 #delimiter='\t', 

671 #newline='\n', 

672 #comments='', 

673 #header=TSVheader) 

674 #else: 

675 return output 

676 

677 

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

679 """ 

680 This function computes the displacement as the mean displacement from the neighbours, for non-converged grains using a TSV file obtained from `spam-ddic` script. Returns a new TSV file with the new Phi (composed only by the displacement part).  

681 

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

683 

684 Parameters 

685 ----------- 

686 lab : 3D array of integers 

687 Labelled volume, with lab.max() labels 

688 

689 DVC : dictionary 

690 Dictionary with deformation field, obtained from `spam-ddic` script, and read using `spam.helpers.tsvio.readCorrelationTSV()` with `readConvergence=True, readPSCC=True, readLabelDilate=True` 

691 

692 fileName : string 

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

694 

695 method : string, optional 

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

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

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

699 Default = 'getLabel' 

700  

701 neighboursRange : int 

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

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

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

705 

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

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

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

709 

710 previousDVC : dictionary, optional 

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

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

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

714 Default = None 

715 

716 Returns 

717 -------- 

718 Dictionary 

719 TSV file with the same columns as the input 

720 """ 

721 import spam.label 

722 

723 # Compute centreOfMass if needed 

724 if centresOfMass == None: 

725 centresOfMass = spam.label.centresOfMass(labIm) 

726 # Get number of labels 

727 numberOfLabels = (labIm.max() + 1).astype('u4') 

728 # Create Phi field 

729 PhiField = numpy.zeros((numberOfLabels, 4, 4), dtype='<f4') 

730 # Rest of arrays 

731 try: 

732 iterations = DVC['iterations'] 

733 returnStatus = DVC['returnStatus'] 

734 deltaPhiNorm = DVC['deltaPhiNorm'] 

735 PSCC = DVC['pixelSearchCC'] 

736 labelDilateList = DVC['LabelDilate'] 

737 error = DVC['error'] 

738 

739 # Get the problematic labels 

740 probLab = numpy.where(DVC['returnStatus']!= 2)[0] 

741 # Remove the 0 from the wrongLab list 

742 probLab = probLab[~numpy.in1d(probLab, 0)] 

743 # Get neighbours 

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

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

746 for i in range(numberOfLabels): 

747 PhiField[i] = DVC['PhiField'][i] 

748 # Work on the problematic labels 

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

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

751 pbar.start() 

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

753 wrongLab = probLab[i] 

754 neighboursLabel = neighbours[i] 

755 t = [] 

756 for n in neighboursLabel: 

757 # Check the return status of each neighbour  

758 if DVC['returnStatus'][n] == 2: 

759 # We dont have a previous DVC file loaded 

760 if previousDVC is None: 

761 # Get translation, rotation and zoom from Phi 

762 nPhi = spam.deformation.deformationFunction.decomposePhi(DVC['PhiField'][n]) 

763 # Append the results 

764 t.append(nPhi['t']) 

765 # We have a previous DVC file loaded 

766 else: 

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

768 nPhi = spam.deformation.deformationFunction.decomposePhi(DVC['PhiField'][n]) 

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

770 nPhiP = spam.deformation.deformationFunction.decomposePhi(previousDVC['PhiField'][n]) 

771 # Append the incremental results 

772 t.append(nPhi['t']-nPhiP['t']) 

773 # Transform list to array 

774 if not t: 

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

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

777 PhiField[wrongLab] = Phi 

778 else: 

779 t = numpy.asarray(t) 

780 # Compute mean 

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

782 # Reconstruct  

783 transformation = {'t': meanT} 

784 Phi = spam.deformation.computePhi(transformation) 

785 # Save 

786 if previousDVC is None: 

787 PhiField[wrongLab] = Phi 

788 else: 

789 PhiField[wrongLab] = previousDVC['PhiField'][wrongLab] 

790 # Add the incremental displacement 

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

792 # Update the progressbar 

793 widgets[0] = progressbar.FormatLabel("{}/{} ".format(i, numberOfLabels)) 

794 pbar.update(i) 

795 # Save 

796 outMatrix = numpy.array([numpy.array(range(numberOfLabels)), 

797 centresOfMass[:, 0], centresOfMass[:, 1], centresOfMass[:, 2], 

798 PhiField[:, 0, 3], PhiField[:, 1, 3], PhiField[:, 2, 3],