Coverage for /usr/local/lib/python3.8/site-packages/spam/deformation/deformationFunction.py: 90%

1# Library of SPAM functions for manipulating a deformation function Phi, which is a 4x4 matrix.

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/>.

17# 2020-02-24 Olga Stamati and Edward Ando

19import numpy

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

23# Value at which to consider no rotation to avoid numerical issues

24rotationAngleDegThreshold = 0.0001

26# Value at which to consider no stretch to avoid numerical issues

27distanceFromIdentity = 0.00001

30###########################################################

31# From components (translation, rotation, zoom, shear) compute Phi

32###########################################################

33def computePhi(transformation, PhiCentre=[0.0, 0.0, 0.0], PhiPoint=[0.0, 0.0, 0.0]):

34 """

35 Builds "Phi", a 4x4 deformation function from a dictionary of transformation parameters (translation, rotation, zoom, shear).

36 Phi can be used to deform coordinates as follows:

37 $$ Phi.x = x'$$

39 Parameters

40 ----------

41 transformation : dictionary of 3x1 arrays

42 Input to computeTransformationOperator is a "transformation" dictionary where all items are optional.

44 Keys

45 t = translation (z,y,x). Note: (0, 0, 0) does nothing

46 r = rotation in "rotation vector" format. Note: (0, 0, 0) does nothing

47 z = zoom. Note: (1, 1, 1) does nothing

48 s = "shear". Note: (0, 0, 0) does nothing

49 U = Right stretch tensor

51 PhiCentre : 3x1 array, optional

52 Point where Phi is centered (centre of rotation)

54 PhiPoint : 3x1 array, optional

55 Point where Phi is going to be applied

57 Returns

58 -------

59 Phi : 4x4 array of floats

60 Phi, deformation function

62 Note

63 ----

64 Useful reference: Chapter 2 -- Rigid Body Registration -- John Ashburner & Karl J. Friston, although we use a symmetric shear.

65 Here we follow the common choice of applying components to Phi in the following order: U or (z or s), r, t

66 """

67 Phi = numpy.eye(4)

69 # Stretch tensor overrides 'z' and 's', hence elif next

70 if "U" in transformation:

71 # symmetric check from https://stackoverflow.com/questions/42908334/checking-if-a-matrix-is-symmetric-in-numpy

72 assert numpy.allclose(transformation["U"], transformation["U"].T), "U not symmetric"

73 tmp = numpy.eye(4)

74 tmp[:3, :3] = transformation["U"]

75 Phi = numpy.dot(tmp, Phi)

76 if "z" in transformation.keys():

77 print("spam.deform.computePhi(): You passed U and z, z is being ignored")

78 if "s" in transformation.keys():

79 print("spam.deform.computePhi(): You passed U and s, s is being ignored")

81 # Zoom + Shear

82 elif "z" in transformation or "s" in transformation:

83 tmp = numpy.eye(4)

85 if "z" in transformation:

86 # Zoom components

87 tmp[0, 0] = transformation["z"][0]

88 tmp[1, 1] = transformation["z"][1]

89 tmp[2, 2] = transformation["z"][2]

91 if "s" in transformation:

92 # Shear components

93 tmp[0, 1] = transformation["s"][0]

94 tmp[0, 2] = transformation["s"][1]

95 tmp[1, 2] = transformation["s"][2]

96 # Shear components

97 tmp[1, 0] = transformation["s"][0]

98 tmp[2, 0] = transformation["s"][1]

99 tmp[2, 1] = transformation["s"][2]

100 Phi = numpy.dot(tmp, Phi)

101

102 # Rotation

103 if "r" in transformation:

104 # https://en.wikipedia.org/wiki/Rodrigues'_rotation_formula

105 # its length is the rotation angle

106 rotationAngleDeg = numpy.linalg.norm(transformation["r"])

107

108 if rotationAngleDeg > rotationAngleDegThreshold:

109 # its direction is the rotation axis.

110 rotationAxis = transformation["r"] / rotationAngleDeg

111

112 # positive angle is clockwise

113 K = numpy.array(

114 [

115 [0, -rotationAxis[2], rotationAxis[1]],

116 [rotationAxis[2], 0, -rotationAxis[0]],

117 [-rotationAxis[1], rotationAxis[0], 0],

118 ]

119 )

120

121 # Note the numpy.dot is very important.

122 R = numpy.eye(3) + (numpy.sin(numpy.deg2rad(rotationAngleDeg)) * K) + ((1.0 - numpy.cos(numpy.deg2rad(rotationAngleDeg))) * numpy.dot(K, K))

123

124 tmp = numpy.eye(4)

125 tmp[0:3, 0:3] = R

126

127 Phi = numpy.dot(tmp, Phi)

128

129 # Translation:

130 if "t" in transformation:

131 Phi[0:3, 3] += transformation["t"]

132

133 # compute distance between point to apply Phi and the point where Phi is centered (centre of rotation)

134 dist = numpy.array(PhiPoint) - numpy.array(PhiCentre)

135

136 # apply Phi to the given point and calculate its displacement

137 Phi[0:3, 3] -= dist - numpy.dot(Phi[0:3, 0:3], dist)

138

139 # check that determinant of Phi is sound

140 if numpy.linalg.det(Phi) < 0.00001:

141 print("spam.deform.computePhi(): Determinant of Phi is very small, this is probably bad, transforming volume into a point.")

142

143 return Phi

144

145

146###########################################################

147# Polar Decomposition of a given F into human readable components

148###########################################################

149def decomposeF(F, twoD=False, verbose=False):

150 """

151 Get components out of a transformation gradient tensor F

152

153 Parameters

154 ----------

155 F : 3x3 array

156 The transformation gradient tensor F (I + du/dx)

157

158 twoD : bool, optional

159 is the F defined in 2D? This applies corrections to the strain outputs

160 Default = False

161

162 verbose : boolean, optional

163 Print errors?

164 Default = True

165

166 Returns

167 -------

168 transformation : dictionary of arrays

169

170 - r = 3x1 numpy array: Rotation in "rotation vector" format

171 - z = 3x1 numpy array: Zoom in "zoom vector" format (z, y, x)

172 - U = 3x3 numpy array: Right stretch tensor, U - numpy.eye(3) is the strain tensor in large strains

173 - e = 3x3 numpy array: Strain tensor in small strains (symmetric)

174 - vol = float: First invariant of the strain tensor ("Volumetric Strain"), det(F)-1

175 - dev = float: Second invariant of the strain tensor ("Deviatoric Strain")

176 - volss = float: First invariant of the strain tensor ("Volumetric Strain") in small strains, trace(e)/3

177 - devss = float: Second invariant of the strain tensor ("Deviatoric Strain") in small strains

178

179 """

180 # - G = 3x3 numpy array: Eigen vectors * eigenvalues of strains, from which principal directions of strain can be obtained

181 # Default non-success values to be over written if all goes well

182 transformation = {

183 "r": numpy.array([numpy.nan] * 3),

184 "z": numpy.array([numpy.nan] * 3),

185 "U": numpy.eye(3) * numpy.nan,

186 "e": numpy.eye(3) * numpy.nan,

187 "vol": numpy.nan,

188 "dev": numpy.nan,

189 "volss": numpy.nan,

190 "devss": numpy.nan

191 # 'G': 3x3

192 }

193

194 # Check for NaNs if any quit

195 if numpy.isnan(F).sum() > 0:

196 if verbose:

197 print("deformationFunction.decomposeF(): Nan value in F. Exiting")

198 return transformation

199

200 # Check for inf if any quit

201 if numpy.isinf(F).sum() > 0:

202 if verbose:

203 print("deformationFunction.decomposeF(): Inf value in F. Exiting")

204 return transformation

205

206 # Check for singular F if yes quit

207 try:

208 numpy.linalg.inv(F)

209 except numpy.linalg.linalg.LinAlgError:

210 if verbose:

211 print("deformationFunction.decomposeF(): F is singular. Exiting")

212 return transformation

213

214 ###########################################################

215 # Polar decomposition of F = RU

216 # U is the right stretch tensor

217 # R is the rotation tensor

218 ###########################################################

219

220 # Compute the Right Cauchy tensor

221 C = numpy.dot(F.T, F)

222

223 # 2020-02-24 EA and OS (day of the fire in 3SR): catch the case when C is practically the identity matrix (i.e., a rigid body motion)

224 # TODO: At least also catch the case when two eigenvales are very small

225 if numpy.abs(numpy.subtract(C, numpy.eye(3))).sum() < distanceFromIdentity:

226 # This forces the rest of the function to give trivial results

227 C = numpy.eye(3)

228

229 # Solve eigen problem

230 CeigVal, CeigVec = numpy.linalg.eig(C)

231

232 # 2018-06-29 OS & ER check for negative eigenvalues

233 # test "really" negative eigenvalues

234 if CeigVal.any() / CeigVal.mean() < -1:

235 print("deformationFunction.decomposeF(): negative eigenvalue in transpose(F). Exiting")

236 print("Eigenvalues are: {}".format(CeigVal))

237 exit()

238 # for negative eigen values but close to 0 we set it to 0

239 CeigVal[CeigVal < 0] = 0

240

241 # Diagonalise C --> which is U**2

242 diagUsqr = numpy.array([[CeigVal[0], 0, 0], [0, CeigVal[1], 0], [0, 0, CeigVal[2]]])

243

244 diagU = numpy.sqrt(diagUsqr)

245

246 # 2018-02-16 check for both issues with negative (Udiag)**2 values and inverse errors

247 try:

248 U = numpy.dot(numpy.dot(CeigVec, diagU), CeigVec.T)

249 R = numpy.dot(F, numpy.linalg.inv(U))

250 except numpy.linalg.LinAlgError:

251 return transformation

252

253 # normalisation of rotation matrix in order to respect basic properties

254 # otherwise it gives errors like trace(R) > 3

255 # this issue might come from numerical noise.

256 # ReigVal, ReigVec = numpy.linalg.eig(R)

257 for i in range(3):

258 R[i, :] /= numpy.linalg.norm(R[i, :])

259 # print("traceR - sumEig = {}".format(R.trace() - ReigVal.sum()))

260 # print("u.v = {}".format(numpy.dot(R[:, 0], R[:, 1])))

261 # print("detR = {}".format(numpy.linalg.det(R)))

262

263 # Calculate rotation angle

264 # Detect an identity -- zero rotation

265 # if numpy.allclose(R, numpy.eye(3), atol=1e-03):

266 # rotationAngleRad = 0.0

267 # rotationAngleDeg = 0.0

268

269 # Detect trace(R) > 3 (should not happen but happens)

270 arccosArg = 0.5 * (R.trace() - 1.0)

271 if arccosArg > 1.0:

272 rotationAngleRad = 0.0

273 else:

274 # https://en.wikipedia.org/wiki/Rotation_formalisms_in_three_dimensions#Rotation_matrix_.E2.86.94_Euler_axis.2Fangle

275 rotationAngleRad = numpy.arccos(arccosArg)

276 rotationAngleDeg = numpy.rad2deg(float(rotationAngleRad))

277

278 if rotationAngleDeg > rotationAngleDegThreshold:

279 rotationAxis = numpy.array([R[2, 1] - R[1, 2], R[0, 2] - R[2, 0], R[1, 0] - R[0, 1]])

280 rotationAxis /= 2.0 * numpy.sin(rotationAngleRad)

281 rot = rotationAngleDeg * rotationAxis

282 else:

283 rot = [0.0, 0.0, 0.0]

284 ###########################################################

285

286 # print "R is \n", R, "\n"

287 # print "|R| is ", numpy.linalg.norm(R), "\n"

288 # print "det(R) is ", numpy.linalg.det(R), "\n"

289 # print "R.T - R-1 is \n", R.T - numpy.linalg.inv( R ), "\n\n"

290

291 # print "U is \n", U, "\n"

292 # print "U-1 is \n", numpy.linalg.inv( U ), "\n\n"

293

294 # Also output eigenvectors * their eigenvalues as output:

295 # G = []

296 # for eigenvalue, eigenvector in zip(CeigVal, CeigVec):

297 # G.append(numpy.multiply(eigenvalue, eigenvector))

298

299 # Compute the volumetric strain from the determinant of F

300 vol = numpy.linalg.det(F) - 1

301

302 # Decompose U into an isotropic and deviatoric part

303 # and compute the deviatoric strain as the norm of the deviatoric part

304 if twoD:

305 Udev = U[1:, 1:] * (numpy.linalg.det(F[1:, 1:]) ** (-1 / 2.0))

306 dev = numpy.linalg.norm(Udev - numpy.eye(2))

307 else:

308 Udev = U * (numpy.linalg.det(F) ** (-1 / 3.0))

309 dev = numpy.linalg.norm(Udev - numpy.eye(3))

310

311 ###########################################################

312 # Small strains bit

313 ###########################################################

314 # to get rid of numerical noise in 2D

315 if twoD:

316 F[0, :] = [1.0, 0.0, 0.0]

317 F[:, 0] = [1.0, 0.0, 0.0]

318

319 # In small strains: 0.5(F+F.T)

320 e = 0.5 * (F + F.T) - numpy.eye(3)

321

322 # The volumetric strain is the trace of the strain matrix

323 volss = numpy.trace(e)

324

325 # The deviatoric in the norm of the matrix

326 if twoD:

327 devss = numpy.linalg.norm(e[1:, 1:] - numpy.eye(2) * volss / 2.0)

328 else:

329 devss = numpy.linalg.norm(e - numpy.eye(3) * volss / 3.0)

330

331 transformation = {

332 "r": rot,

333 "z": [U[i, i] for i in range(3)],

334 "U": U,

335 # 'G': G,

336 "e": e,

337 "vol": vol,

338 "dev": dev,

339 "volss": volss,

340 "devss": devss,

341 }

342

343 return transformation

344

345

346def decomposePhi(Phi, PhiCentre=[0.0, 0.0, 0.0], PhiPoint=[0.0, 0.0, 0.0], twoD=False, verbose=False):

347 """

348 Get components out of a linear deformation function "Phi"

349

350 Parameters

351 ----------

352 Phi : 4x4 array

353 The deformation function operator "Phi"

354

355 PhiCentre : 3x1 array, optional

356 Point where Phi was calculated

357

358 PhiPoint : 3x1 array, optional

359 Point where Phi is going to be applied

360

361 twoD : bool, optional

362 is the Phi defined in 2D? This applies corrections to the strain outputs

363 Default = False

364

365 verbose : boolean, optional

366 Print errors?

367 Default = True

368

369 Returns

370 -------

371 transformation : dictionary of arrays

372

373 - t = 3x1 numpy array: Translation vector (z, y, x)

374 - r = 3x1 numpy array: Rotation in "rotation vector" format

375 - z = 3x1 numpy array: Zoom in "zoom vector" format (z, y, x)

376 - U = 3x3 numpy array: Right stretch tensor, U - numpy.eye(3) is the strain tensor in large strains

377 - e = 3x3 numpy array: Strain tensor in small strains (symmetric)

378 - vol = float: First invariant of the strain tensor ("Volumetric Strain"), det(F)-1

379 - dev = float: Second invariant of the strain tensor ("Deviatoric Strain")

380 - volss = float: First invariant of the strain tensor ("Volumetric Strain") in small strains, trace(e)/3

381 - devss = float: Second invariant of the strain tensor ("Deviatoric Strain") in small strains

382

383 """

384 # - G = 3x3 numpy array: Eigen vectors * eigenvalues of strains, from which principal directions of strain can be obtained

385 # Default non-success values to be over written if all goes well

386 transformation = {

387 "t": numpy.array([numpy.nan] * 3),

388 "r": numpy.array([numpy.nan] * 3),

389 "z": numpy.array([numpy.nan] * 3),

390 "U": numpy.eye(3) * numpy.nan,

391 "e": numpy.eye(3) * numpy.nan,

392 "vol": numpy.nan,

393 "dev": numpy.nan,

394 "volss": numpy.nan,

395 "devss": numpy.nan

396 # 'G': 3x3

397 }

398

399 # Check for singular Phi if yes quit

400 try:

401 numpy.linalg.inv(Phi)

402 except numpy.linalg.linalg.LinAlgError:

403 return transformation

404

405 # Check for NaNs if any quit

406 if numpy.isnan(Phi).sum() > 0:

407 return transformation

408

409 # Check for NaNs if any quit

410 if numpy.isinf(Phi).sum() > 0:

411 return transformation

412

413 ###########################################################

414 # F, the inside 3x3 displacement gradient

415 ###########################################################

416 F = Phi[0:3, 0:3].copy()

417

418 ###########################################################

419 # Calculate transformation by undoing F on the PhiPoint

420 ###########################################################

421 tra = Phi[0:3, 3].copy()

422

423 # compute distance between given point and the point where Phi was calculated

424 dist = numpy.array(PhiPoint) - numpy.array(PhiCentre)

425

426 # apply Phi to the given point and calculate its displacement

427 tra -= dist - numpy.dot(F, dist)

428

429 decomposedF = decomposeF(F, verbose=verbose)

430

431 transformation = {

432 "t": tra,

433 "r": decomposedF["r"],

434 "z": decomposedF["z"],

435 "U": decomposedF["U"],

436 # 'G': G,

437 "e": decomposedF["e"],

438 "vol": decomposedF["vol"],

439 "dev": decomposedF["dev"],

440 "volss": decomposedF["volss"],

441 "devss": decomposedF["devss"],

442 }

443

444 return transformation

445

446

447def computeRigidPhi(Phi):

448 """

449 Returns only the rigid part of the passed Phi

450

451 Parameters

452 ----------

453 Phi : 4x4 numpy array of floats

454 Phi, deformation function

455

456 Returns

457 -------

458 PhiRigid : 4x4 numpy array of floats

459 Phi with only the rotation and translation parts on inputted Phi

460 """

461 decomposedPhi = decomposePhi(Phi)

462 PhiRigid = computePhi({"t": decomposedPhi["t"], "r": decomposedPhi["r"]})

463 return PhiRigid