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

1# Library of SPAM functions for deforming images.

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

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

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

7# any later version.

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

10# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or

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

12# more details.

13#

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

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

17import multiprocessing

19try:

20 multiprocessing.set_start_method("fork")

21except RuntimeError:

22 pass

23import numpy

24import scipy.ndimage

25from spambind.DIC.DICToolkit import applyMeshTransformation as applyMeshTransformationCPP

26from spambind.DIC.DICToolkit import applyPhi as applyPhiCPP

27from spambind.DIC.DICToolkit import binningChar, binningFloat, binningUInt

29nProcessesDefault = multiprocessing.cpu_count()

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

33###########################################################

34# Take an Phi and apply it (C++) to an image

35###########################################################

36def applyPhi(im, Phi=None, PhiCentre=None, interpolationOrder=1):

37 """

38 Deform a 3D image using a deformation function "Phi", applied using spam's C++ interpolator.

39 Only interpolation order = 1 is implemented.

41 Parameters

42 ----------

43 im : 3D numpy array

44 3D numpy array of grey levels to be deformed

46 Phi : 4x4 array, optional

47 "Phi" deformation function.

48 Highly recommended additional argument (why are you calling this function otherwise?)

50 PhiCentre : 3x1 array of floats, optional

51 Centre of application of Phi.

52 Default = (numpy.array(im1.shape)-1)/2.0

53 i.e., the centre of the image

55 interpolationOrder : int, optional

56 Order of image interpolation to use, options are either 0 (strict nearest neighbour) or 1 (trilinear interpolation)

57 Default = 1

59 Returns

60 -------

61 imDef : 3D array

62 Deformed greyscales by Phi

63 """

65 # Detect 2D images, and bail, doesn't work with our interpolator

66 if len(im.shape) == 2 or (numpy.array(im.shape) == 1).any():

67 print("DIC.deformationFunction.applyPhi(): looks like a 2D image which cannot be handled. Please use DIC.deformationFunction.applyPhiPython")

68 return

70 # Sort out Phi and calculate inverse

71 if Phi is None:

72 PhiInv = numpy.eye(4, dtype="<f4")

73 else:

74 try:

75 PhiInv = numpy.linalg.inv(Phi).astype("<f4")

76 except numpy.linalg.LinAlgError:

77 # print( "\tapplyPhi(): Can't inverse Phi, setting it to identity matrix. Phi is:\n{}".format( Phi ) )

78 PhiInv = numpy.eye(4, dtype="<f4")

80 if PhiCentre is None:

81 PhiCentre = (numpy.array(im.shape) - 1) / 2.0

83 if interpolationOrder > 1:

84 print("DIC.deformationFunction.applyPhi(): interpolation Order > 1 not implemented")

85 return

87 im = im.astype("<f4")

88 PhiCentre = numpy.array(PhiCentre).astype("<f4")

89 # We need to inverse Phi for question of direction

90 imDef = numpy.zeros_like(im, dtype="<f4")

91 applyPhiCPP(

92 im.astype("<f4"),

93 imDef,

94 PhiInv.astype("<f4"),

95 PhiCentre.astype("<f4"),

96 int(interpolationOrder),

97 )

99 return imDef

100

101

102###########################################################

103# Take an Phi and apply it to an image

104###########################################################

105def applyPhiPython(im, Phi=None, PhiCentre=None, interpolationOrder=3):

106 """

107 Deform a 3D image using a deformation function "Phi", applied using scipy.ndimage.map_coordinates

108 Can have orders > 1 but is hungry in memory.

109

110 Parameters

111 ----------

112 im : 2D/3D numpy array

113 2D/3D numpy array of grey levels to be deformed

114

115 Phi : 4x4 array, optional

116 "Phi" linear deformation function.

117 Highly recommended additional argument (why are you calling this function otherwise?)

118

119 PhiCentre : 3x1 array of floats, optional

120 Centre of application of Phi.

121 Default = (numpy.array(im1.shape)-1)/2.0

122 i.e., the centre of the image

123

124 interpolationOrder : int, optional

125 Order of image interpolation to use. This value is passed directly to ``scipy.ndimage.map_coordinates`` as "order".

126 Default = 3

127

128 Returns

129 -------

130 imSub : 3D array

131 Deformed greyscales by Phi

132 """

133

134 if Phi is None:

135 PhiInv = numpy.eye(4, dtype="<f4")

136 else:

137 try:

138 PhiInv = numpy.linalg.inv(Phi).astype("<f4")

139 except numpy.linalg.LinAlgError:

140 # print( "\tapplyPhiPython(): Can't inverse Phi, setting it to identity matrix. Phi is:\n{}".format( Phi ) )

141 PhiInv = numpy.eye(4)

142

143 if im.ndim == 2:

144 twoD = True

145 im = im[numpy.newaxis, ...]

146 else:

147 twoD = False

148

149 if PhiCentre is None:

150 PhiCentre = (numpy.array(im.shape) - 1) / 2.0

151

152 imDef = numpy.zeros_like(im, dtype="<f4")

153

154 coordinatesInitial = numpy.ones((4, im.shape[0] * im.shape[1] * im.shape[2]), dtype="<f4")

155

156 coordinates_mgrid = numpy.mgrid[0 : im.shape[0], 0 : im.shape[1], 0 : im.shape[2]]

157

158 # Copy into coordinatesInitial

159 coordinatesInitial[0, :] = coordinates_mgrid[0].ravel() - PhiCentre[0]

160 coordinatesInitial[1, :] = coordinates_mgrid[1].ravel() - PhiCentre[1]

161 coordinatesInitial[2, :] = coordinates_mgrid[2].ravel() - PhiCentre[2]

162

163 # Apply Phi to coordinates

164 coordinatesDef = numpy.dot(PhiInv, coordinatesInitial)

165

166 coordinatesDef[0, :] += PhiCentre[0]

167 coordinatesDef[1, :] += PhiCentre[1]

168 coordinatesDef[2, :] += PhiCentre[2]

169

170 imDef += scipy.ndimage.map_coordinates(im, coordinatesDef[0:3], order=interpolationOrder).reshape(imDef.shape).astype("<f4")

171

172 if twoD:

173 imDef = imDef[0]

174

175 return imDef

176

177

178###############################################################

179# Take a field of Phi and apply it (quite slowly) to an image

180###############################################################

181def applyPhiField(

182 im,

183 fieldCoordsRef,

184 PhiField,

185 imMaskDef=None,

186 displacementMode="applyPhi",

187 nNeighbours=8,

188 interpolationOrder=1,

189 nProcesses=nProcessesDefault,

190 verbose=False,

191):

192 """

193 Deform a 3D image using a field of deformation functions "Phi" coming from a regularGrid,

194 applied using scipy.ndimage.map_coordinates.

195

196 Parameters

197 ----------

198 im : 3D array

199 3D array of grey levels to be deformed

200

201 fieldCoordsRef: 2D array, optional

202 nx3 array of n points coordinates (ZYX)

203 centre where each deformation function "Phi" has been calculated

204

205 PhiField: 3D array, optional

206 nx4x4 array of n points deformation functions

207

208 imMaskDef: 3D array of bools, optional

209 3D array same size as im but DEFINED IN THE DEFORMED CONFIGURATION

210 which should be True in the pixels to fill in in the deformed configuration.

211 Default = None

212

213 displacementMode : string, optional

214 How do you want to calculate displacements?

215 With "interpolate" they are just interpolated from the PhiField

216 With "applyPhi" each neighbour's Phi function is applied to the pixel position

217 and the resulting translations weighted and summed.

218 Default = "applyPhi"

219

220 nNeighbours : int, optional

221 Number of the nearest neighbours to consider

222 #This OR neighbourRadius must be set.

223 Default = 8

224

225 interpolationOrder : int, optional

226 Order of image interpolation to use. This value is passed directly to ``scipy.ndimage.map_coordinates`` as "order".

227 Default = 1

228

229 nProcesses : integer, optional

230 Number of processes for multiprocessing

231 Default = number of CPUs in the system

232

233 verbose : boolean, optional

234 Print progress?

235 Default = True

236

237 Returns

238 -------

239 imDef : 3D array

240 deformed greylevels by a field of deformation functions "Phi"

241 """

242

243 import progressbar

244 import spam.deformation

245 import spam.DIC

246

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

248

249 # print("making pixel grid")

250 if imMaskDef is not None:

251 # print("...from a passed mask, cool, this should shave time")

252 pixCoordsDef = numpy.array(numpy.where(imMaskDef)).T

253 else:

254 # Create the grid of the input image

255 imSize = im.shape

256 coordinates_mgrid = numpy.mgrid[0 : imSize[0], 0 : imSize[1], 0 : imSize[2]]

257

258 pixCoordsDef = numpy.ones((imSize[0] * imSize[1] * imSize[2], 3))

259

260 pixCoordsDef[:, 0] = coordinates_mgrid[0].ravel()

261 pixCoordsDef[:, 1] = coordinates_mgrid[1].ravel()

262 pixCoordsDef[:, 2] = coordinates_mgrid[2].ravel()

263 # print(pixCoordsDef.shape)

264

265 # Initialise deformed coordinates

266 fieldCoordsDef = fieldCoordsRef + PhiField[:, 0:3, -1]

267 # print("done")

268

269 maskPhiField = numpy.isfinite(PhiField[:, 0, 0])

270

271 if displacementMode == "interpolate":

272 """

273 In this mode we're only taking into account displacements.

274 We use interpolatePhiField in the deformed configuration, in displacements only,

275 and we don't feed PhiInv, but only the negative of the displacements

276 """

277 backwardsDisplacementsPhi = numpy.zeros_like(PhiField)

278 backwardsDisplacementsPhi[:, 0:3, -1] = -1 * PhiField[:, 0:3, -1]

279

280 pixDispsDef = spam.DIC.interpolatePhiField(

281 fieldCoordsDef[maskPhiField],

282 backwardsDisplacementsPhi[maskPhiField],

283 pixCoordsDef,

284 nNeighbours=nNeighbours,

285 interpolateF="no",

286 nProcesses=nProcesses,

287 verbose=verbose,

288 )

289 pixCoordsRef = pixCoordsDef + pixDispsDef[:, 0:3, -1]

290

291 elif displacementMode == "applyPhi":

292 """

293 In this mode we're NOT interpolating the displacement field.

294 For each pixel, we're applying the neighbouring Phis and looking

295 at the resulting displacement of the pixel.

296 Those different displacements are weighted as a function of distance

297 and averaged into the point's final displacement.

298

299 Obviously if your PhiField is only a displacement field, this changes

300 nothing from above (except for computation time), but with some stretches

301 this can become interesting.

302 """

303 # print("inversing PhiField")

304 PhiFieldInv = numpy.zeros_like(PhiField)

305 for n in range(PhiField.shape[0]):

306 try:

307 PhiFieldInv[n] = numpy.linalg.inv(PhiField[n])

308 except numpy.linalg.LinAlgError:

309 maskPhiField[n] = False

310 # print("done")

311

312 # mask everything

313 PhiFieldInvMasked = PhiFieldInv[maskPhiField]

314 fieldCoordsRefMasked = fieldCoordsRef[maskPhiField]

315 fieldCoordsDefMasked = fieldCoordsDef[maskPhiField]

316

317 # build KD-tree

318 treeCoordDef = scipy.spatial.KDTree(fieldCoordsDefMasked)

319

320 pixCoordsRef = numpy.zeros_like(pixCoordsDef, dtype="f4")

321

322 """

323 Define multiproc function only for displacementMode == "applyPhi"

324 """

325 global _multiprocessingComputeDisplacementForSeriesOfPixels

326

327 def _multiprocessingComputeDisplacementForSeriesOfPixels(seriesNumber):

328 pixelNumbers = splitPixNumbers[seriesNumber]

329

330 pixCoordsRefSeries = numpy.zeros((len(pixelNumbers), 3), dtype="f4")

331

332 # all jobs should take the same time, so just show progress bar in 0th process

333 if seriesNumber == 0 and verbose:

334 pbar = progressbar.ProgressBar(maxval=len(pixelNumbers)).start()

335

336 for localPixelNumber, globalPixelNumber in enumerate(pixelNumbers):

337 if seriesNumber == 0 and verbose:

338 pbar.update(localPixelNumber)

339

340 pixCoordDef = pixCoordsDef[globalPixelNumber]

341 # get nNeighbours and compute distance weights

342 distances, indices = treeCoordDef.query(pixCoordDef, k=nNeighbours)

343 weights = 1 / (distances + tol)

344

345 displacement = numpy.zeros(3, dtype="float")

346

347 # for each neighbour

348 for neighbour, index in enumerate(indices):

349 # apply PhiInv to current point with PhiCentre = fieldCoordsREF <- this is important

350 # -> this gives a displacement for each neighbour

351 PhiInv = PhiFieldInvMasked[index]

352 # print("PhiInv", PhiInv)

353 translationTmp = PhiInv[0:3, -1].copy()

354 dist = pixCoordDef - fieldCoordsRefMasked[index]

355 # print(f"dist = {dist}")

356 translationTmp -= dist - numpy.dot(PhiInv[0:3, 0:3], dist)

357 # print(f"translationTmp ({neighbour}): {translationTmp} (weight = {weights[neighbour]})")

358 displacement += translationTmp * weights[neighbour]

359

360 # compute resulting displacement as weights * displacements

361 # compute pixel coordinates in reference config

362 # print("pixCoordDef", pixCoordDef)

363 # print("displacement", displacement)

364 pixCoordsRefSeries[localPixelNumber] = pixCoordDef + displacement / weights.sum()

365

366 if seriesNumber == 0 and verbose:

367 pbar.finish()

368 return pixelNumbers, pixCoordsRefSeries

369 # pixCoordsRef[pixNumber] = pixCoordDef + numpy.sum(displacements*weights[:, numpy.newaxis], axis=0)

370

371 splitPixNumbers = numpy.array_split(numpy.arange(pixCoordsDef.shape[0]), nProcesses)

372

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

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

375 for returns in pool.imap_unordered(_multiprocessingComputeDisplacementForSeriesOfPixels, range(nProcesses)):

376 pixCoordsRef[returns[0]] = returns[1]

377 pool.close()

378 pool.join()

379

380 if imMaskDef is not None:

381 imDef = numpy.zeros_like(im)

382 imDef[imMaskDef] = scipy.ndimage.map_coordinates(im, pixCoordsRef.T, mode="constant", order=interpolationOrder)

383

384 else: # no pixel mask, image comes out directly

385 imDef = scipy.ndimage.map_coordinates(im, pixCoordsRef.T, mode="constant", order=interpolationOrder).reshape(im.shape)

386

387 return imDef

388

389

390def binning(im, binning, returnCropAndCentre=False):

391 """

392 This function downscales images by averaging NxNxN voxels together in 3D and NxN pixels in 2D.

393 This is useful for reducing data volumes, and denoising data (due to averaging procedure).

394

395 Parameters

396 ----------

397 im : 2D/3D numpy array

398 Input measurement field

399

400 binning : int

401 The number of pixels/voxels to average together

402

403 returnCropAndCentre: bool (optional)

404 Return the position of the centre of the binned image

405 in the coordinates of the original image, and the crop

406 Default = False

407

408 Returns

409 -------

410 imBin : 2/3D numpy array

411 `binning`-binned array

412

413 (otherwise if returnCropAndCentre): list containing:

414 imBin,

415 topCrop, bottomCrop

416 centre of imBin in im coordinates (useful for re-stitching)

417 Notes

418 -----

419 Here we will only bin pixels/voxels if they is a sufficient number of

420 neighbours to perform the binning. This means that the number of pixels that

421 will be rejected is the dimensions of the image, modulo the binning amount.

422

423 The returned volume is computed with only fully binned voxels, meaning that some voxels on the edges

424 may be excluded.

425 This means that the output volume size is the input volume size / binning or less (in fact the crop

426 in the input volume is the input volume size % binning

427 """

428 twoD = False

429

430 if im.dtype == "f8":

431 im = im.astype("<f4")

432

433 binning = int(binning)

434 # print("binning = ", binning)

435

436 dimsOrig = numpy.array(im.shape)

437 # print("dimsOrig = ", dimsOrig)

438

439 # Note: // is a floor-divide

440 imBin = numpy.zeros(dimsOrig // binning, dtype=im.dtype)

441 # print("imBin.shape = ", imBin.shape)

442

443 # Calculate number of pixels to throw away

444 offset = dimsOrig % binning

445 # print("offset = ", offset)

446

447 # Take less off the top corner than the bottom corner

448 topCrop = offset // 2

449 # print("topCrop = ", topCrop)

450 topCrop = topCrop.astype("<i2")

451

452 if len(im.shape) == 2:

453 # pad them

454 im = im[numpy.newaxis, ...]

455 imBin = imBin[numpy.newaxis, ...]

456 topCrop = numpy.array([0, topCrop[0], topCrop[1]]).astype("<i2")

457 offset = numpy.array([0, offset[0], offset[1]]).astype("<i2")

458 twoD = True

459

460 # Call C++

461 if im.dtype == "f4":

462 # print("Float binning")

463 binningFloat(im.astype("<f4"), imBin, topCrop.astype("<i4"), int(binning))

464 elif im.dtype == "u2":

465 # print("Uint 2 binning")

466 binningUInt(im.astype("<u2"), imBin, topCrop.astype("<i4"), int(binning))

467 elif im.dtype == "u1":

468 # print("Char binning")

469 binningChar(im.astype("<u1"), imBin, topCrop.astype("<i4"), int(binning))

470

471 if twoD:

472 imBin = imBin[0]

473

474 if returnCropAndCentre:

475 centreBinned = (numpy.array(imBin.shape) - 1) / 2.0

476 relCentOrig = offset + binning * centreBinned

477 return [imBin, [topCrop, offset - topCrop], relCentOrig]

478 else:

479 return imBin

480

481

482###############################################################

483# Take a tetrahedral mesh (defined by coords and conn) and use

484# it to deform an image

485###############################################################

486def applyMeshTransformation(im, points, connectivity, displacements, imTetLabel=None, nThreads=1):

487 """

488 This function deforms an image based on a tetrahedral mesh and

489 nodal displacements (normally from Global DVC),

490 using the mesh's shape functions to interpolate.

491

492 Parameters

493 ----------

494 im : 3D numpy array of greyvalues

495 Input image to be deformed

496

497 points : m x 3 numpy array

498 M nodal coordinates in reference configuration

499

500 connectivity : n x 4 numpy array

501 Tetrahedral connectivity generated by spam.mesh.triangulate() for example

502

503 displacements : m x 3 numpy array

504 M displacements defined at the nodes

505

506 imTetLabel : 3D numpy array of ints (optional)

507 Pixels labelled with the tetrahedron (i.e., line number in connectivity matrix) they belong to.

508 If this is not passed, it's calculated in this function (can be slow).

509 WARNING: This is in the deformed configuration.

510

511 nThreads: int (optional, default=1)

512 The number of threads used for the cpp parallelisation.

513

514 Returns

515 -------

516 imDef : 3D numpy array of greyvalues

517 Deformed image

518

519 """

520 # deformed tetrahedra

521 pointsDef = points + displacements

522

523 if imTetLabel is None:

524 import spam.label

525

526 print("spam.DIC.applyMeshTransformation(): imTetLabel not passed, recomputing it.")

527 imTetLabel = spam.label.labelTetrahedra(im.shape, pointsDef, connectivity, nThreads=nThreads)

528

529 # Allocate output array that will be painted in by C++

530 imDef = numpy.zeros_like(im, dtype="<f4")

531 applyMeshTransformationCPP(

532 im.astype("<f4"),

533 imTetLabel.astype("<u4"),

534 imDef,

535 connectivity.astype("<u4"),

536 pointsDef.astype("<f8"),

537 displacements.astype("<f8"),

538 nThreads,

539 )

540 return imDef