Coverage for /usr/local/lib/python3.10/site-packages/spam/deformation/deformationFunction.py: 90%
125 statements
« prev ^ index » next coverage.py v7.8.0, created at 2025-04-24 17:26 +0000
1# Library of SPAM functions for manipulating a deformation function Phi, which is a 4x4 matrix.
2# Copyright (C) 2020 SPAM Contributors
3#
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.
8#
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)
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"])
108 if rotationAngleDeg > rotationAngleDegThreshold:
109 # its direction is the rotation axis.
110 rotationAxis = transformation["r"] / rotationAngleDeg
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 )
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))
124 tmp = numpy.eye(4)
125 tmp[0:3, 0:3] = R
127 Phi = numpy.dot(tmp, Phi)
129 # Translation:
130 if "t" in transformation:
131 Phi[0:3, 3] += transformation["t"]
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)
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)
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.")
143 return Phi
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
153 Parameters
154 ----------
155 F : 3x3 array
156 The transformation gradient tensor F (I + du/dx)
158 twoD : bool, optional
159 is the F defined in 2D? This applies corrections to the strain outputs
160 Default = False
162 verbose : boolean, optional
163 Print errors?
164 Default = True
166 Returns
167 -------
168 transformation : dictionary of arrays
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 - Up = 3x1 numpy array: Principal strain components (min (z), med (y), max (x)) from diagonalisation of U tensor
174 - e = 3x3 numpy array: Strain tensor in small strains (symmetric)
175 - vol = float: First invariant of the strain tensor ("Volumetric Strain"), det(F)-1
176 - dev = float: Second invariant of the strain tensor ("Deviatoric Strain")
177 - volss = float: First invariant of the strain tensor ("Volumetric Strain") in small strains, trace(e)/3
178 - devss = float: Second invariant of the strain tensor ("Deviatoric Strain") in small strains
180 """
181 # - G = 3x3 numpy array: Eigen vectors * eigenvalues of strains, from which principal directions of strain can be obtained
182 # Default non-success values to be over written if all goes well
183 transformation = {
184 "r": numpy.array([numpy.nan] * 3),
185 "z": numpy.array([numpy.nan] * 3),
186 "Up": numpy.array([numpy.nan] * 3),
187 "U": numpy.eye(3) * numpy.nan,
188 "e": numpy.eye(3) * numpy.nan,
189 "vol": numpy.nan,
190 "dev": numpy.nan,
191 "volss": numpy.nan,
192 "devss": numpy.nan
193 # 'G': 3x3
194 }
196 # Check for NaNs if any quit
197 if numpy.isnan(F).sum() > 0:
198 if verbose:
199 print("deformationFunction.decomposeF(): Nan value in F. Exiting")
200 return transformation
202 # Check for inf if any quit
203 if numpy.isinf(F).sum() > 0:
204 if verbose:
205 print("deformationFunction.decomposeF(): Inf value in F. Exiting")
206 return transformation
208 # Check for singular F if yes quit
209 try:
210 numpy.linalg.inv(F)
211 except numpy.linalg.LinAlgError:
212 if verbose:
213 print("deformationFunction.decomposeF(): F is singular. Exiting")
214 return transformation
216 ###########################################################
217 # Polar decomposition of F = RU
218 # U is the right stretch tensor
219 # R is the rotation tensor
220 ###########################################################
222 # Compute the Right Cauchy tensor
223 C = numpy.dot(F.T, F)
225 # 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)
226 # TODO: At least also catch the case when two eigenvales are very small
227 if numpy.abs(numpy.subtract(C, numpy.eye(3))).sum() < distanceFromIdentity:
228 # This forces the rest of the function to give trivial results
229 C = numpy.eye(3)
231 # Solve eigen problem
232 CeigVal, CeigVec = numpy.linalg.eig(C)
234 # 2018-06-29 OS & ER check for negative eigenvalues
235 # test "really" negative eigenvalues
236 if CeigVal.any() / CeigVal.mean() < -1:
237 print("deformationFunction.decomposeF(): negative eigenvalue in transpose(F). Exiting")
238 print("Eigenvalues are: {}".format(CeigVal))
239 exit()
240 # for negative eigen values but close to 0 we set it to 0
241 CeigVal[CeigVal < 0] = 0
243 # Diagonalise C --> which is U**2
244 diagUsqr = numpy.array([[CeigVal[0], 0, 0], [0, CeigVal[1], 0], [0, 0, CeigVal[2]]])
246 diagU = numpy.sqrt(diagUsqr)
248 # 2018-02-16 check for both issues with negative (Udiag)**2 values and inverse errors
249 try:
250 U = numpy.dot(numpy.dot(CeigVec, diagU), CeigVec.T)
251 R = numpy.dot(F, numpy.linalg.inv(U))
252 except numpy.linalg.LinAlgError:
253 return transformation
255 # normalisation of rotation matrix in order to respect basic properties
256 # otherwise it gives errors like trace(R) > 3
257 # this issue might come from numerical noise.
258 # ReigVal, ReigVec = numpy.linalg.eig(R)
259 for i in range(3):
260 R[i, :] /= numpy.linalg.norm(R[i, :])
261 # print("traceR - sumEig = {}".format(R.trace() - ReigVal.sum()))
262 # print("u.v = {}".format(numpy.dot(R[:, 0], R[:, 1])))
263 # print("detR = {}".format(numpy.linalg.det(R)))
265 # Calculate rotation angle
266 # Detect an identity -- zero rotation
267 # if numpy.allclose(R, numpy.eye(3), atol=1e-03):
268 # rotationAngleRad = 0.0
269 # rotationAngleDeg = 0.0
271 # Detect trace(R) > 3 (should not happen but happens)
272 arccosArg = 0.5 * (R.trace() - 1.0)
273 if arccosArg > 1.0:
274 rotationAngleRad = 0.0
275 else:
276 # https://en.wikipedia.org/wiki/Rotation_formalisms_in_three_dimensions#Rotation_matrix_.E2.86.94_Euler_axis.2Fangle
277 rotationAngleRad = numpy.arccos(arccosArg)
278 rotationAngleDeg = numpy.rad2deg(float(rotationAngleRad))
280 if rotationAngleDeg > rotationAngleDegThreshold:
281 rotationAxis = numpy.array([R[2, 1] - R[1, 2], R[0, 2] - R[2, 0], R[1, 0] - R[0, 1]])
282 rotationAxis /= 2.0 * numpy.sin(rotationAngleRad)
283 rot = rotationAngleDeg * rotationAxis
284 else:
285 rot = [0.0, 0.0, 0.0]
286 ###########################################################
288 # print "R is \n", R, "\n"
289 # print "|R| is ", numpy.linalg.norm(R), "\n"
290 # print "det(R) is ", numpy.linalg.det(R), "\n"
291 # print "R.T - R-1 is \n", R.T - numpy.linalg.inv( R ), "\n\n"
293 # print "U is \n", U, "\n"
294 # print "U-1 is \n", numpy.linalg.inv( U ), "\n\n"
296 # Also output eigenvectors * their eigenvalues as output:
297 # G = []
298 # for eigenvalue, eigenvector in zip(CeigVal, CeigVec):
299 # G.append(numpy.multiply(eigenvalue, eigenvector))
301 # Compute the volumetric strain from the determinant of F
302 vol = numpy.linalg.det(F) - 1
304 # Decompose U into an isotropic and deviatoric part
305 # and compute the deviatoric strain as the norm of the deviatoric part
306 if twoD:
307 Udev = U[1:, 1:] * (numpy.linalg.det(F[1:, 1:]) ** (-1 / 2.0))
308 dev = numpy.linalg.norm(Udev - numpy.eye(2))
309 else:
310 Udev = U * (numpy.linalg.det(F) ** (-1 / 3.0))
311 dev = numpy.linalg.norm(Udev - numpy.eye(3))
313 ###########################################################
314 # Small strains bit
315 ###########################################################
316 # to get rid of numerical noise in 2D
317 if twoD:
318 F[0, :] = [1.0, 0.0, 0.0]
319 F[:, 0] = [1.0, 0.0, 0.0]
321 # In small strains: 0.5(F+F.T)
322 e = 0.5 * (F + F.T) - numpy.eye(3)
324 # The volumetric strain is the trace of the strain matrix
325 volss = numpy.trace(e)
327 # The deviatoric in the norm of the matrix
328 if twoD:
329 devss = numpy.linalg.norm(e[1:, 1:] - numpy.eye(2) * volss / 2.0)
330 else:
331 devss = numpy.linalg.norm(e - numpy.eye(3) * volss / 3.0)
333 transformation = {
334 "r": rot,
335 "z": [U[i, i] for i in range(3)],
336 "U": U,
337 "Up": numpy.sort(numpy.diag(diagU - numpy.eye(3))),
338 # 'G': G,
339 "e": e,
340 "vol": vol,
341 "dev": dev,
342 "volss": volss,
343 "devss": devss,
344 }
346 return transformation
349def decomposePhi(Phi, PhiCentre=[0.0, 0.0, 0.0], PhiPoint=[0.0, 0.0, 0.0], twoD=False, verbose=False):
350 """
351 Get components out of a linear deformation function "Phi"
353 Parameters
354 ----------
355 Phi : 4x4 array
356 The deformation function operator "Phi"
358 PhiCentre : 3x1 array, optional
359 Point where Phi was calculated
361 PhiPoint : 3x1 array, optional
362 Point where Phi is going to be applied
364 twoD : bool, optional
365 is the Phi defined in 2D? This applies corrections to the strain outputs
366 Default = False
368 verbose : boolean, optional
369 Print errors?
370 Default = True
372 Returns
373 -------
374 transformation : dictionary of arrays
376 - t = 3x1 numpy array: Translation vector (z, y, x)
377 - r = 3x1 numpy array: Rotation in "rotation vector" format
378 - z = 3x1 numpy array: Zoom in "zoom vector" format (z, y, x)
379 - U = 3x3 numpy array: Right stretch tensor, U - numpy.eye(3) is the strain tensor in large strains
380 - Up = 3x1 numpy array: Principal strain components (min, med, max) from diagonalisation of U tensor
381 - e = 3x3 numpy array: Strain tensor in small strains (symmetric)
382 - vol = float: First invariant of the strain tensor ("Volumetric Strain"), det(F)-1
383 - dev = float: Second invariant of the strain tensor ("Deviatoric Strain")
384 - volss = float: First invariant of the strain tensor ("Volumetric Strain") in small strains, trace(e)/3
385 - devss = float: Second invariant of the strain tensor ("Deviatoric Strain") in small strains
387 """
388 # - G = 3x3 numpy array: Eigen vectors * eigenvalues of strains, from which principal directions of strain can be obtained
389 # Default non-success values to be over written if all goes well
390 transformation = {
391 "t": numpy.array([numpy.nan] * 3),
392 "r": numpy.array([numpy.nan] * 3),
393 "z": numpy.array([numpy.nan] * 3),
394 "Up": numpy.array([numpy.nan] * 3),
395 "U": numpy.eye(3) * numpy.nan,
396 "e": numpy.eye(3) * numpy.nan,
397 "vol": numpy.nan,
398 "dev": numpy.nan,
399 "volss": numpy.nan,
400 "devss": numpy.nan
401 # 'G': 3x3
402 }
404 # Check for singular Phi if yes quit
405 try:
406 numpy.linalg.inv(Phi)
407 except numpy.linalg.LinAlgError:
408 return transformation
410 # Check for NaNs if any quit
411 if numpy.isnan(Phi).sum() > 0:
412 return transformation
414 # Check for NaNs if any quit
415 if numpy.isinf(Phi).sum() > 0:
416 return transformation
418 ###########################################################
419 # F, the inside 3x3 displacement gradient
420 ###########################################################
421 F = Phi[0:3, 0:3].copy()
423 ###########################################################
424 # Calculate transformation by undoing F on the PhiPoint
425 ###########################################################
426 tra = Phi[0:3, 3].copy()
428 # compute distance between given point and the point where Phi was calculated
429 dist = numpy.array(PhiPoint) - numpy.array(PhiCentre)
431 # apply Phi to the given point and calculate its displacement
432 tra -= dist - numpy.dot(F, dist)
434 decomposedF = decomposeF(F, verbose=verbose)
436 transformation = {
437 "t": tra,
438 "r": decomposedF["r"],
439 "z": decomposedF["z"],
440 "U": decomposedF["U"],
441 "Up": decomposedF["Up"],
442 # 'G': G,
443 "e": decomposedF["e"],
444 "vol": decomposedF["vol"],
445 "dev": decomposedF["dev"],
446 "volss": decomposedF["volss"],
447 "devss": decomposedF["devss"],
448 }
450 return transformation
453def computeRigidPhi(Phi):
454 """
455 Returns only the rigid part of the passed Phi
457 Parameters
458 ----------
459 Phi : 4x4 numpy array of floats
460 Phi, deformation function
462 Returns
463 -------
464 PhiRigid : 4x4 numpy array of floats
465 Phi with only the rotation and translation parts on inputted Phi
466 """
467 decomposedPhi = decomposePhi(Phi)
468 PhiRigid = computePhi({"t": decomposedPhi["t"], "r": decomposedPhi["r"]})
469 return PhiRigid