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

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

19import numpy

20import scipy.ndimage

21from spam.DIC.DICToolkit import applyMeshTransformation as applyMeshTransformationCPP

22from spam.DIC.DICToolkit import applyPhi as applyPhiCPP

23from spam.DIC.DICToolkit import binningChar, binningFloat, binningUInt

25nProcessesDefault = multiprocessing.cpu_count()

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

29###########################################################

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

31###########################################################

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

33 """

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

35 Only interpolation order = 1 is implemented.

37 Parameters

38 ----------

39 im : 3D numpy array

40 3D numpy array of grey levels to be deformed

42 Phi : 4x4 array, optional

43 "Phi" deformation function.

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

46 PhiCentre : 3x1 array of floats, optional

47 Centre of application of Phi.

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

49 i.e., the centre of the image

51 interpolationOrder : int, optional

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

53 Default = 1

55 Returns

56 -------

57 imDef : 3D array

58 Deformed greyscales by Phi

59 """

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

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

63 print(

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

65 )

66 return

68 # Sort out Phi and calculate inverse

69 if Phi is None:

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

71 else:

72 try:

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

74 except numpy.linalg.linalg.LinAlgError:

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

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

78 if PhiCentre is None:

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

81 if interpolationOrder > 1:

82 print(

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

84 )

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 : 3D numpy array

113 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.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 PhiCentre is None:

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

145

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

147

148 coordinatesInitial = numpy.ones(

149 (4, im.shape[0] * im.shape[1] * im.shape[2]), dtype="<f4"

150 )

151

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

153

154 # Copy into coordinatesInitial

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

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

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

158

159 # Apply Phi to coordinates

160 coordinatesDef = numpy.dot(PhiInv, coordinatesInitial)

161

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

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

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

165

166 imDef += (

167 scipy.ndimage.map_coordinates(im, coordinatesDef[0:3], order=interpolationOrder)

168 .reshape(imDef.shape)

169 .astype("<f4")

170 )

171 return imDef

172

173

174###############################################################

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

176###############################################################

177def applyPhiField(

178 im,

179 fieldCoordsRef,

180 PhiField,

181 imMaskDef=None,

182 displacementMode="applyPhi",

183 nNeighbours=8,

184 interpolationOrder=1,

185 nProcesses=nProcessesDefault,

186 verbose=False,

187):

188 """

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

190 applied using scipy.ndimage.map_coordinates.

191

192 Parameters

193 ----------

194 im : 3D array

195 3D array of grey levels to be deformed

196

197 fieldCoordsRef: 2D array, optional

198 nx3 array of n points coordinates (ZYX)

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

200

201 PhiField: 3D array, optional

202 nx4x4 array of n points deformation functions

203

204 imMaskDef: 3D array of bools, optional

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

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

207 Default = None

208

209 displacementMode : string, optional

210 How do you want to calculate displacements?

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

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

213 and the resulting translations weighted and summed.

214 Default = "applyPhi"

215

216 nNeighbours : int, optional

217 Number of the nearest neighbours to consider

218 #This OR neighbourRadius must be set.

219 Default = 8

220

221 interpolationOrder : int, optional

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

223 Default = 1

224

225 nProcesses : integer, optional

226 Number of processes for multiprocessing

227 Default = number of CPUs in the system

228

229 verbose : boolean, optional

230 Print progress?

231 Default = True

232

233 Returns

234 -------

235 imDef : 3D array

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

237 """

238 import multiprocessing

239

240 import progressbar

241 import spam.deformation

242 import spam.DIC

243

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

245

246 # print("making pixel grid")

247 if imMaskDef is not None:

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

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

250 else:

251 # Create the grid of the input image

252 imSize = im.shape

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

254

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

256

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

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

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

260 # print(pixCoordsDef.shape)

261

262 # Initialise deformed coordinates

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

264 # print("done")

265

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

267

268 if displacementMode == "interpolate":

269 """

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

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

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

273 """

274 backwardsDisplacementsPhi = numpy.zeros_like(PhiField)

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

276

277 pixDispsDef = spam.DIC.interpolatePhiField(

278 fieldCoordsDef[maskPhiField],

279 backwardsDisplacementsPhi[maskPhiField],

280 pixCoordsDef,

281 nNeighbours=nNeighbours,

282 interpolateF="no",

283 nProcesses=nProcesses,

284 verbose=verbose,

285 )

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

287

288 elif displacementMode == "applyPhi":

289 """

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

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

292 at the resulting displacement of the pixel.

293 Those different displacements are weighted as a function of distance

294 and averaged into the point's final displacement.

295

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

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

298 this can become interesting.

299 """

300 # print("inversing PhiField")

301 PhiFieldInv = numpy.zeros_like(PhiField)

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

303 try:

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

305 except numpy.linalg.LinAlgError:

306 maskPhiField[n] = False

307 # print("done")

308

309 # mask everything

310 PhiFieldInvMasked = PhiFieldInv[maskPhiField]

311 fieldCoordsRefMasked = fieldCoordsRef[maskPhiField]

312 fieldCoordsDefMasked = fieldCoordsDef[maskPhiField]

313

314 # build KD-tree

315 treeCoordDef = scipy.spatial.KDTree(fieldCoordsDefMasked)

316

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

318

319 """

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

321 """

322 global computeDisplacementForSeriesOfPixels

323

324 def computeDisplacementForSeriesOfPixels(seriesNumber):

325 pixelNumbers = splitPixNumbers[seriesNumber]

326

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

328

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

330 if seriesNumber == 0 and verbose:

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

332

333 for localPixelNumber, globalPixelNumber in enumerate(pixelNumbers):

334 if seriesNumber == 0 and verbose:

335 pbar.update(localPixelNumber)

336

337 pixCoordDef = pixCoordsDef[globalPixelNumber]

338 # get nNeighbours and compute distance weights

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

340 weights = 1 / (distances + tol)

341

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

343

344 # for each neighbour

345 for neighbour, index in enumerate(indices):

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

347 # -> this gives a displacement for each neighbour

348 PhiInv = PhiFieldInvMasked[index]

349 # print("PhiInv", PhiInv)

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

351 dist = pixCoordDef - fieldCoordsRefMasked[index]

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

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

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

355 displacement += translationTmp * weights[neighbour]

356

357 # compute resulting displacement as weights * displacements

358 # compute pixel coordinates in reference config

359 # print("pixCoordDef", pixCoordDef)

360 # print("displacement", displacement)

361 pixCoordsRefSeries[localPixelNumber] = (

362 pixCoordDef + displacement / weights.sum()

363 )

364

365 if seriesNumber == 0 and verbose:

366 pbar.finish()

367 return pixelNumbers, pixCoordsRefSeries

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

369

370 splitPixNumbers = numpy.array_split(

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

372 )

373

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

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

376 for returns in pool.imap_unordered(

377 computeDisplacementForSeriesOfPixels, range(nProcesses)

378 ):

379 pixCoordsRef[returns[0]] = returns[1]

380 pool.close()

381 pool.join()

382

383 if imMaskDef is not None:

384 imDef = numpy.zeros_like(im)

385 imDef[imMaskDef] = scipy.ndimage.map_coordinates(

386 im, pixCoordsRef.T, mode="constant", order=interpolationOrder

387 )

388

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

390 imDef = scipy.ndimage.map_coordinates(

391 im, pixCoordsRef.T, mode="constant", order=interpolationOrder

392 ).reshape(im.shape)

393

394 return imDef

395

396

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

398 """

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

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

401

402 Parameters

403 ----------

404 im : 2D/3D numpy array

405 Input measurement field

406

407 binning : int

408 The number of pixels/voxels to average together

409

410 returnCropAndCentre: bool (optional)

411 Return the position of the centre of the binned image

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

413 Default = False

414

415 Returns

416 -------

417 imBin : 2/3D numpy array

418 `binning`-binned array

419

420 (otherwise if returnCropAndCentre): list containing:

421 imBin,

422 topCrop, bottomCrop

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

424 Notes

425 -----

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

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

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

429

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

431 may be excluded.

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

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

434 """

435 twoD = False

436

437 if im.dtype == "f8":

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

439

440 binning = int(binning)

441 # print("binning = ", binning)

442

443 dimsOrig = numpy.array(im.shape)

444 # print("dimsOrig = ", dimsOrig)

445

446 # Note: // is a floor-divide

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

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

449

450 # Calculate number of pixels to throw away

451 offset = dimsOrig % binning

452 # print("offset = ", offset)

453

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

455 topCrop = offset // 2

456 # print("topCrop = ", topCrop)

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

458

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

460 # pad them

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

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

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

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

465 twoD = True

466

467 # Call C++

468 if im.dtype == "f4":

469 # print("Float binning")

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

471 elif im.dtype == "u2":

472 # print("Uint 2 binning")

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

474 elif im.dtype == "u1":

475 # print("Char binning")

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

477

478 if twoD:

479 imBin = imBin[0]

480

481 if returnCropAndCentre:

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

483 relCentOrig = offset + binning * centreBinned

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

485 else:

486 return imBin

487

488

489###############################################################

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

491# it to deform an image

492###############################################################

493def applyMeshTransformation(

494 im, points, connectivity, displacements, imTetLabel=None, nThreads=1

495):

496 """

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

498 nodal displacements (normally from Global DVC),

499 using the mesh's shape functions to interpolate.

500

501 Parameters

502 ----------

503 im : 3D numpy array of greyvalues

504 Input image to be deformed

505

506 points : m x 3 numpy array

507 M nodal coordinates in reference configuration

508

509 connectivity : n x 4 numpy array

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

511

512 displacements : m x 3 numpy array

513 M displacements defined at the nodes

514

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

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

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

518 WARNING: This is in the deformed configuration.

519

520 nThreads: int (optional, default=1)

521 The number of threads used for the cpp parallelisation.

522

523 Returns

524 -------

525 imDef : 3D numpy array of greyvalues

526 Deformed image

527

528 """

529 # deformed tetrahedra

530 pointsDef = points + displacements

531

532 if imTetLabel is None:

533 import spam.label

534

535 print(

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

537 )

538 imTetLabel = spam.label.labelTetrahedra(

539 im.shape, pointsDef, connectivity, nThreads=nThreads

540 )

541

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

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

544 applyMeshTransformationCPP(

545 im.astype("<f4"),

546 imTetLabel.astype("<u4"),

547 imDef,

548 connectivity.astype("<u4"),

549 pointsDef.astype("<f8"),

550 displacements.astype("<f8"),

551 nThreads,

552 )

553 return imDef