Coverage for /usr/local/lib/python3.8/site-packages/spam/DIC/globalDVC.py: 85%
279 statements
« prev ^ index » next coverage.py v7.2.3, created at 2023-11-22 13:26 +0000
1# Library of SPAM image correlation functions.
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/>.
17import time # for debug
19import numpy
20import spam.DIC
21import spam.label # for label tet
22import tifffile
24# 2017-05-29 ER and EA
25# This is spam's C++ DIC toolkit, but since we're in the tools/ directory we can import it directly
26from spam.DIC.DICToolkit import (
27 computeDICglobalMatrix,
28 computeDICglobalVector,
29 computeGradientPerTet,
30)
33def _errorCalc(im1, im2, im2ref, meshPaddingSlice):
34 errorInitial = numpy.sqrt(numpy.square(im2ref[meshPaddingSlice] - im1[meshPaddingSlice]).sum())
35 errorCurrent = numpy.sqrt(numpy.square(im2[meshPaddingSlice] - im1[meshPaddingSlice]).sum())
36 return errorCurrent / errorInitial
39def _makeRegularisationMatrix(regularisation, points, connectivity, Mc):
41 # 1. check minimal number of parameters
42 # 1.1: check mesh type
43 try:
44 MESH = regularisation["MESH"]
45 MESH_TYPE = MESH.get("type", "cuboid")
46 except KeyError as e:
47 raise KeyError(f"[regularisation] Missing parameter: {e}")
49 # 1.2: check mandatory BULK
50 try:
51 BULK = regularisation["BULK"]
52 BULK_YOUNG = BULK["young"]
53 BULK_POISSON = BULK["poisson"]
54 BULK_KSI = BULK["ksi"]
55 except KeyError as e:
56 raise KeyError(f"[regularisation] Missing parameter in BULK: {e}")
57 # 1.2 check cast not mandatory BULK
58 BULK_FEW = BULK.get("few_times", 3.0) # we still need to figure out what it means
59 print(f"[regularisation] [neumann] young = {BULK_YOUNG}")
60 print(f"[regularisation] [neumann] poisson = {BULK_POISSON}")
61 print(f"[regularisation] [neumann] ksi = {BULK_KSI}")
62 print(f"[regularisation] [neumann] few = {BULK_FEW}")
64 # 1.3: check dirichlet
65 DIRICHLET = regularisation.get("DIRICHLET")
66 try:
67 if DIRICHLET:
68 DIRICHLET_SURFACES = [(k, v) for k, v in DIRICHLET["ksi"].items() if isinstance(v, int)]
69 else:
70 DIRICHLET_SURFACES = []
71 except KeyError as e:
72 raise KeyError(f"Missing parameter in DIRICHLET: {e}")
74 print(f'[regularisation] [dirichlet] Surface regularisation activated: {"yes" if DIRICHLET else "no"}')
75 if DIRICHLET:
76 for k, v in DIRICHLET_SURFACES:
77 # k is where, v is ksi
78 print(f"[regularisation] [dirichlet]\t{k}: {v}")
80 # 2. compute complete connectivity matrix
81 K = spam.mesh.globalStiffnessMatrix(points, connectivity, BULK_YOUNG, BULK_POISSON)
83 # 3. compute projection matrices
84 # note: for cylinders we can't deduce the surface from the mesh
85 # so it will need to be a user input in `regularisation`
86 # diagonal matrix with 1 only on surface dof
87 DS = numpy.zeros_like(K)
88 # diagonal matrix with 1 only on surface with dirichlet dof for each
89 # surface
90 n_dirichlet = len(DIRICHLET_SURFACES)
91 print(f"[regularisation] [dirichlet] number of surfaces -> {n_dirichlet}")
92 DSd = []
93 DSd_A = []
94 DSd_B = []
95 DSd_nn = []
96 DSd_en = []
97 for si in range(n_dirichlet):
98 DSd.append(numpy.zeros_like(K))
99 DSd_A.append(numpy.zeros_like(K)) # A matrix (to be inverted)
100 DSd_B.append(numpy.zeros_like(K)) # B matrix (A.B = I)
101 DSd_nn.append([]) # list of all node numbers on dirichlet surfaces
102 DSd_en.append([]) # list of all equations on dirichlet surfaces
104 # diagonal matrix with 1 only on surface with neuman dof
105 DSn = numpy.zeros_like(K)
106 if MESH_TYPE == "cuboid":
107 # we assume it's a cuboid to define the surfaces
108 # 1 x 6 array with [min_z, min_y, min_x, max_z, max_y, max_x]
109 maxCoord = numpy.amax(points, axis=0).astype("<u2")
110 minCoord = numpy.amin(points, axis=0).astype("<u2")
111 surfaces_positions = list(minCoord) + list(maxCoord)
112 print(f"[regularisation] [dirichlet] compute surfaces based on min/max mesh coordinates: {surfaces_positions}")
114 # a list of [direction, positions fo reach surface]
115 def _dirpos(key, surfaces_positions):
116 """Simple parsing function that takes surface keys
117 z_start, z_end, y_start, y_end, x_start, x_end
118 and convert it to (direction, position) based on surfaces_positions
120 surface_positions = [min_z, min_y, min_x, max_z, max_y, max_x]
121 """
122 # split key
123 dir_str, pos_str = key.split("_")
125 # convert direction into integer
126 _str2int = {"z": 0, "y": 1, "x": 2}
127 direction = _str2int[dir_str]
129 # convert position into coordinate
130 _str2int = {"s": 0, "e": 1}
131 position = surfaces_positions[3 * _str2int[pos_str[0]] + direction]
133 print(f"[regularisation] [dirichlet] BC {key} -> {direction}, {position}")
135 return direction, position # needs to be a tuple for later lookup
137 bcd_lookup = [_dirpos(key, surfaces_positions) for key, _ in DIRICHLET_SURFACES]
138 if DIRICHLET:
139 print(f"[regularisation] [dirichlet] BC bcd_lookup = {bcd_lookup}")
141 # loop over the points
142 # PUT SOME NUMBA HERE
143 for A, point in enumerate(points):
144 # print(f"[regularisation] point {A} {point} ({A/points.shape[0]*100:.2f}%)")
145 # loop over the surfaces to find if current point belongs to it
146 # A: global point number
147 # point: [z, y, x] coordinate
149 # get all the dofs for surfaces: DS
150 for direction, position in enumerate(surfaces_positions):
151 # since surfaces_positions =
152 # [min_z, min_y, min_x, max_z, max_y, max_x]
153 # position is alternatively min_z, min_y, ...
154 # with mod 3 direction is 0, 1, 2, 0, 1, 2
155 direction = direction % 3 # due to concatenation of min/max
156 # check point position
157 if abs(point[direction] - position) < 0.00001:
158 # print(f"[regularisation] Point number {A}: {point} in surface N{direction} = {position}")
159 # loop over the 3 dofs / deduce global equation number
160 for i in range(3):
161 # P: global equation number
162 P = 3 * A + i
163 DS[P, P] = 1
165 # check if also dirichlet
166 if DIRICHLET:
167 try:
168 # WARNING: position is float so we need to be
169 # more careful about the comparison we make here
170 # as one comes from user input and the other from
171 # node coordinates
172 # maybe use numpy.isclose
173 si = bcd_lookup.index((direction, position))
174 DSd[si][P, P] = 1
175 DSd_nn[si].append(A)
176 DSd_en[si].append(P)
177 # print(f'Dirichlet Surface {si} point {A} {P}')
178 except ValueError:
179 DSn[P, P] = 1
180 else:
181 DSn[P, P] = 1
183 # with the break the point will be assigned to the
184 # first surface of the loop (ie we could miss some
185 # dirichlet bc if neuman comes first).
186 # And since we don't have the control over the ordering
187 # of the surface it's not a really good thing
188 # break
190 # WARNING/TODO: without the break above, points on
191 # edges can be assigne to multiple surfaces
192 # for edges between two surfaces of the same kind
193 # it's okay but edges between neuman and dirchlet
194 # pose a problem.
195 # Might be good to clean the neumann matrix where there overlap to give
196 # priority to dirichlet
198 elif MESH_TYPE == "cylinder":
199 # later check the good parameters needed (center/radius/orientation)
200 raise NotImplementedError("Cylindrical meshes are not implement yet")
201 else:
202 raise KeyError(f"Unknown mesh type {MESH_TYPE}")
204 # loop over all tetraedra to build up triangle meshes of dirichlet
205 # convert node numbers into sets
206 dirichlet_connectivity_all = []
207 dirichlet_points_all = []
208 for si in range(n_dirichlet):
209 DSd_nn[si] = set(DSd_nn[si])
210 dirichlet_connectivity_all.append([])
211 dirichlet_points_all.append([])
213 # loop over dirichlet matrix
214 for si, surface_node_numbers in enumerate(DSd_nn):
215 for tetra_connectivity in connectivity:
216 tri_con = list(surface_node_numbers.intersection(tetra_connectivity))
217 tri_nodes = [list(points[i]) for i in tri_con]
218 if len(tri_con) != 3:
219 # the element doesn't have a triangle a this surface
220 continue
222 # get 4th point of the tetrahedron to compute 3D shape functions
223 point_4_n = list(set(tri_con).symmetric_difference(set(tetra_connectivity)))[0]
224 _, tri_coefficients = spam.mesh.shapeFunctions(tri_nodes + [list(points[point_4_n])])
225 # we've got a triangle on the surface!!!
226 # dirichlet_direction = regularisation["bc_dirichlet"][si][2]
227 dirichlet_connectivity_all[si].append(list(tri_con))
228 dirichlet_points_all[si].append(tri_nodes)
230 # assemble the dirichlet connectivity matrix
231 a = numpy.array(tri_coefficients[:3, 0])
232 bcd = numpy.array(tri_coefficients[:3, 1:])
233 B = numpy.array(tri_nodes).T
235 def phi(i, L):
236 return a[i] + numpy.matmul(bcd[i], numpy.matmul(B, L.T))
238 def dphi(i):
239 return bcd[i]
241 # gauss points
242 L_gp = numpy.array([[0.5, 0.5, 0.0], [0.5, 0.0, 0.5], [0.0, 0.5, 0.5]])
244 # STEP 3: compute the area
245 area = 0.5 * numpy.linalg.norm(
246 numpy.cross(
247 numpy.subtract(tri_nodes[1], tri_nodes[0]),
248 numpy.subtract(tri_nodes[2], tri_nodes[0]),
249 )
250 )
252 # STEP 3: compute inner products of the shape functions
253 for i in range(3):
254 for j in range(3):
255 inner = 0.0
256 for L in L_gp:
257 # print(inner, i, j, L, phi(i, L), phi(j, L))
258 inner += phi(i, L) * phi(j, L)
259 inner *= area / 3.0 # the 1/3. comes from the weight
260 dinner = area * numpy.inner(dphi(i), dphi(j))
261 for dirichlet_direction in range(3):
262 P = 3 * tri_con[i] + dirichlet_direction
263 Q = 3 * tri_con[j] + dirichlet_direction
264 DSd_A[si][P, Q] += inner
265 DSd_B[si][P, Q] += dinner
267 for si in range(n_dirichlet):
268 # invert matrix A
269 # 1. extract submatrix A from full size DSd_A and invert
270 sub_A = DSd_A[si][:, DSd_en[si]]
271 sub_A = numpy.linalg.inv(sub_A[DSd_en[si], :])
272 # 2. push back inverted submatrix into full size DSd_A
273 for i, P in enumerate(DSd_en[si]):
274 for j, Q in enumerate(DSd_en[si]):
275 DSd_A[si][P, Q] = sub_A[i, j]
277 # DEBUG
278 # for si, _ in enumerate(dirichlet_connectivity_all):
279 # a = numpy.array(dirichlet_connectivity_all[si])
280 # meshio.write_points_cells(
281 # f'dirichlet_{si}.vtk',
282 # points,
283 # {'triangle': a},
284 # file_format='vtk',
285 # binary=False
286 # )
288 # diagonal matrix with 1 only on bulks dof
289 DB = numpy.eye(K.shape[0]) - DS
291 # 4.1 build bulk connectivity matrix
292 Km = numpy.matmul(DB + DSn, K)
293 # 4.2 build dirichlet connectivity matrices
294 Ks = [numpy.matmul(DSd_i, K) for DSd_i in DSd]
296 # 5. normalisation of each functionals
297 # 5.1: build v(x) = sin(2pi k.x)
298 # lookup for largest dimension to build k
299 lengths = [0, 0, 0]
300 for i in range(3):
301 lengths[i] = surfaces_positions[i + 3] - surfaces_positions[i]
302 largest_direction = numpy.argmax(lengths)
303 largest_size = numpy.max(lengths)
304 k = [0, 0, 0]
305 k_mag = BULK_FEW / float(largest_size)
306 k[largest_direction] = k_mag
307 # print(f"[regularisation] k = {k} |k| = {numpy.linalg.norm(k)} 1/|k| = {1.0 / numpy.linalg.norm(k)}")
309 # build pure shear field for normalisation
310 v = []
311 for point in points:
312 vector = [0, 0, 0]
313 norm = numpy.sin(2 * numpy.pi * numpy.dot(k, point))
314 vector[(largest_direction + 1) % 3] = numpy.sqrt(0.5) * norm
315 vector[(largest_direction + 2) % 3] = numpy.sqrt(0.5) * norm
316 # vector[(largest_direction + 1) % 3] = norm
317 v.append(vector)
318 # print(point, k, numpy.dot(k, point), norm, vector)
320 # debug build vtk with v(x) on current mesh
321 v = numpy.array(v)
322 spam.helpers.writeUnstructuredVTK(
323 points,
324 connectivity,
325 elementType="tetra",
326 pointData={"v": v},
327 cellData={},
328 fileName="debug-regularisation-v.vtk",
329 )
331 # 5.2 compute normalized energies
332 v = numpy.ravel(v)
333 Ec = numpy.matmul(v.T, numpy.matmul(Mc, v))
334 Em = numpy.matmul(v.T, numpy.matmul(Km.T, numpy.matmul(Km, v)))
335 Es = [
336 numpy.matmul(
337 v.T,
338 numpy.matmul(
339 Ks[i].T,
340 numpy.matmul(DSd_A[i], numpy.matmul(DSd_B[i], numpy.matmul(Ks[i], v))),
341 ),
342 )
343 for i in range(n_dirichlet)
344 ]
346 Es_string = " ".join([f"{_:.0f}" for _ in Es])
347 print(f"[regularisation] Ec = {Ec:.0f}, Em = {Em:.0f}, Es = {Es_string}")
349 # 5.3 compute functional weights
350 # 5.3.1 DVC weight
351 omega_c = 1.0
353 # 5.3.2 Neumann weight
354 ksi = BULK_KSI if BULK_KSI else 1.0 / k_mag
355 print(f"[regularisation] [neumann] ksi = {ksi}")
356 omega_m = (ksi * numpy.linalg.norm(k)) ** 4
358 # 5.3.3 Dirichlet weight
359 omega_s = []
360 for _, ksi in DIRICHLET_SURFACES:
361 if ksi == 0: # compute value based on k_mag
362 ksi = 1.0 / k_mag
363 print(f"[regularisation] [dirichlet] use 1/k_mag for ksi for surface {_}: {ksi}")
364 omega_s.append((ksi * numpy.linalg.norm(k)) ** 4)
366 # 5.3.4 Total weigth for normalisation
367 omega_t = omega_c + omega_m + sum(omega_s)
369 print(f"[regularisation] omega_t = {omega_t:.2f}")
370 print(f"[regularisation] omega_c = {omega_c:.2f} ({100 * (omega_c / omega_t):.2f}%)")
371 print(f"[regularisation] omega_m = {omega_m:.2f} ({100 * (omega_m / omega_t):.2f}%)")
372 for i, omega in enumerate(omega_s):
373 print(f"[regularisation] omega_s_{i} = {omega} ({100 * (omega / omega_t):.2f}%)")
375 # weight_m = omega_m * Ec / Em
376 # weight_s = [omega_s[i] * Ec / Es[i] for i in range(n_dirichlet)]
377 # print(f"[regularisation] weight_m = {weight_m}, weight_s = {weight_s}")
379 # 5.4 compute Mreg
380 Mreg = numpy.zeros_like(K)
381 Mreg = Ec * omega_m * numpy.matmul(Km.T, Km) / Em
382 for si in range(n_dirichlet):
383 Mreg += omega_s[si] * Ec * numpy.matmul(Ks[si].T, numpy.matmul(DSd_A[si], numpy.matmul(DSd_B[si], Ks[si]))) / Es[si]
385 # from matplotlib import pyplot as plt
386 # # plt.imshow(numpy.log(Mreg))
387 # plt.imshow(numpy.log(Mreg))
388 # plt.show()
389 # plt.imshow(numpy.log(Mc))
390 # plt.show()
392 # def phi_t(u):
393 # omega_t = omega_c + omega_m
394 # phi_c = omega_c * numpy.matmul(u.T, numpy.matmul(Mc, u)) / Ec
395 # phi_c += omega_m * numpy.matmul(u.T, numpy.matmul(Km.T, numpy.matmul(Km, u))) / Em
396 # phi_c /= omega_t
397 # return phi_c
398 # return
400 return Mreg
403def globalCorrelation(
404 im1,
405 im2,
406 points,
407 connectivity,
408 regularisation={},
409 initialDisplacements=None,
410 convergenceCriterion=0.01,
411 maxIterations=20,
412 medianFilterEachIteration=False,
413 debugFiles=False,
414 prefix="globalCorrelation",
415 nThreads=None,
416):
417 """
418 Global DVC (works only with 3D images).
420 Parameters
421 ----------
422 im1 : 3D numpy array
423 Reference image in which the mesh is defined
425 im2 : 3D numpy array
426 Deformed image, should be same size as im1
428 points : M x 3 numpy array
429 M nodal coordinates in reference configuration
431 connectivity : N x 4 numpy array
432 connectivityhedral connectivity generated by spam.mesh.triangulate() for example
434 regularisation : dict (optional)
435 regularisation parameters
437 .. code-block:: python
439 # ksi is a normalisation parameter to control the weight
440 # of each functionals (DVC, bulk and surface).
441 # If set to 0 it will be automatically computed.
443 # Example for a cuboid mesh
444 {
445 # Information about the mesh
446 "MESH": {"type": "cuboid"}, # mandatory
447 # Information about the elastic regularisation
448 # (Bulk + Neumann boundary conditions)
449 "BULK": {
450 "young": 10, # mandatory (Young modulus)
451 "poisson": 0.25, # mandatory (Poisson ratio)
452 "ksi": 30, # mandatory
453 "few_times": 3, # optionnal (whatever...)
454 },
455 # Information about the surface regularisation
456 # (Dirichlet boundary conditions)
457 # Each surface of the cuboid is labelled by the keywords
458 # z_start: z == 0, z_end, y_start, y_end, x_start and x_end)
459 # If a keyword is ommited the surface is not regularised.
460 "DIRICHLET": {
461 "ksi": { # mandatory
462 "z_start": 0, # automatically computed
463 "z_end": 30, # ksi normalisation is 30
464 }
465 },
466 }
468 initialDisplacements : M x 3 numpy array of floats (optional)
469 Initial guess for nodal displacements, must be coherent with input mesh
470 Default = None
472 convergenceCriterion : float
473 Convergence criterion for change in displacements in px
474 Default = 0.01
476 maxIterations : int
477 Number of iterations to stop after if convergence has not been reached
478 Default = 20
480 debugFiles : bool
481 Write temporary results to file for debugging?
482 Default = 'globalCorrelation'
484 prefix : string
485 Output file prefix for debugFiles
486 Default = None
488 Returns
489 -------
490 displacements : N x 3 numpy array of floats
491 (converged?) Nodal displacements
493 Example
494 -------
495 >>> import spam.DIC
496 >>> spam.DIC.globalCorrelation(
497 imRef,
498 imDef
499 )
500 """
501 import multiprocessing
503 import spam.helpers
504 import spam.mesh
506 # Global number of processes
507 nThreads = multiprocessing.cpu_count() if nThreads is None else nThreads
508 print(f"[globalCorrelation] C++ parallelisation on {nThreads} threads")
510 print(f"[globalCorrelation] Convergence criterion = {convergenceCriterion}")
511 print(f"[globalCorrelation] Max iterations = {maxIterations}")
512 print("[globalCorrelation] Converting im1 to 32-bit float")
513 im1 = im1.astype("<f4")
515 points = points.astype("<f8")
516 connectivity = connectivity.astype("<u4")
518 maxCoord = numpy.amax(points, axis=0).astype("<u2")
519 minCoord = numpy.amin(points, axis=0).astype("<u2")
520 print(f"[globalCorrelation] Mesh box: min = {minCoord} max = {maxCoord}")
522 meshPaddingSlice = (
523 slice(minCoord[0], maxCoord[0]),
524 slice(minCoord[1], maxCoord[1]),
525 slice(minCoord[2], maxCoord[2]),
526 )
528 displacements = numpy.zeros((points.shape[0], 3), dtype="<f8")
530 print(f"[globalCorrelation] Points: {points.shape}")
531 print(f"[globalCorrelation] Displacements: {displacements.shape}")
532 print(f"[globalCorrelation] Cells: {connectivity.shape}")
533 print(f"[globalCorrelation] Padding: {meshPaddingSlice}")
535 ###############################################################
536 # Step 2-1 Apply deformation and interpolate pixels
537 ###############################################################
539 print("[globalCorrelation] Allocating 3D data (deformed image)")
540 if initialDisplacements is None:
541 im1Def = im1.copy()
542 imTetLabel = spam.label.labelTetrahedra(im1.shape, points, connectivity, nThreads=nThreads)
543 else:
544 print("[globalCorrelation] Applying initial deformation to image")
545 displacements = initialDisplacements.copy()
546 tic = time.perf_counter()
547 imTetLabel = spam.label.labelTetrahedra(im1.shape, points + displacements, connectivity, nThreads=nThreads)
548 print(f"[globalCorrelation] Running labelTetrahedra: {time.perf_counter()-tic:.3f} seconds.")
550 im1Def = spam.DIC.applyMeshTransformation(
551 im1,
552 points,
553 connectivity,
554 displacements,
555 imTetLabel=imTetLabel,
556 nThreads=nThreads,
557 )
558 if debugFiles:
559 print("[globalCorrelation] Saving initial images")
560 for name, image in [
561 [f"{prefix}-def-init.tif", im1Def],
562 [f"{prefix}-imTetLabel-init.tif", imTetLabel],
563 ]:
564 print(f"[globalCorrelation]\t{name}: {image.shape}")
565 tifffile.imwrite(name, image)
567 # print("[globalCorrelation] Correlating (MF)!")
568 print("[globalCorrelation] Calculating gradient of IM TWO...")
569 im2Grad = numpy.array(numpy.gradient(im2), dtype="<f4")
571 print("[globalCorrelation] Computing global matrix")
572 # This generates the globalMatrix (big Mc matrix) with imGrad as input
573 Mc = numpy.zeros((3 * points.shape[0], 3 * points.shape[0]), dtype="<f8")
575 if debugFiles:
576 print("[globalCorrelation] Computing debug files fields")
577 gradientPerTet = numpy.zeros((connectivity.shape[0], 3), dtype="<f8")
578 IDPerTet = numpy.array([_ for _ in range(connectivity.shape[0])])
580 computeGradientPerTet(
581 imTetLabel.astype("<u4"),
582 im2Grad.astype("<f4"),
583 connectivity.astype("<u4"),
584 (points + displacements).astype("<f8"),
585 gradientPerTet,
586 )
588 spam.helpers.writeUnstructuredVTK(
589 (points + displacements),
590 connectivity,
591 cellData={"meanGradient": gradientPerTet, "id": IDPerTet},
592 fileName=f"{prefix}-gradient.vtk",
593 )
594 del gradientPerTet
596 computeDICglobalMatrix(
597 imTetLabel.astype("<u4"),
598 im2Grad.astype("<f4"),
599 connectivity.astype("<u4"),
600 (points + displacements).astype("<f8"),
601 Mc,
602 )
604 ###############################################################
605 # Setup left hand vector
606 ###############################################################
607 if len(regularisation):
608 print("[globalCorrelation] Entering regularisation")
609 Mreg = _makeRegularisationMatrix(regularisation, points, connectivity, Mc)
610 left_hand_inverse = numpy.linalg.inv(Mc + Mreg)
611 else:
612 print("[globalCorrelation] Skip regularisation")
613 left_hand_inverse = numpy.linalg.inv(Mc)
614 del Mc
616 # error = _errorCalc(im2, im1Def, im1, meshPaddingSlice)
617 # print("\[globalCorrelation] Initial Error (abs) = ", error)
619 # We try to solve Md=F
620 # while error > 0.1 and error < errorIn:
621 # while error > 0.1 and i <= maxIterations and error < errorIn:
622 dxNorm = numpy.inf
623 i = 0
624 while dxNorm > convergenceCriterion and i < maxIterations:
625 i += 1
627 # This function returns globalVector (F) taking in im1Def and im2 and the gradients
628 tic = time.perf_counter()
629 # print("[globalCorrelation] [newton] run computeDICglobalVector: ", end="")
630 right_hand_vector = numpy.zeros((3 * points.shape[0]), dtype="<f8")
631 computeDICglobalVector(
632 imTetLabel.astype("<u4"),
633 im2Grad.astype("<f4"),
634 im1Def.astype("<f4"),
635 im2.astype("<f4"),
636 connectivity.astype("<u4"),
637 (points + displacements).astype("<f8"),
638 right_hand_vector,
639 )
640 # print(f"{time.perf_counter()-tic:.3f} seconds.")
642 tic = time.perf_counter()
643 # print("[globalCorrelation] [newton] run solve: ", end="")
645 # solve: we can use solve here for sake of precision (over computing
646 # M^-1). However solve takes quite a lot of time for "small" meshes).
648 if len(regularisation):
649 right_hand_vector -= numpy.matmul(Mreg, displacements.ravel())
650 dx = numpy.matmul(left_hand_inverse, right_hand_vector).astype("<f8")
651 # dx_solve = numpy.linalg.solve(
652 # Mc,
653 # right_hand_vector
654 # ).astype('<f8')
655 # print(numpy.linalg.norm(dx - dx_solve))
657 displacements += dx.reshape(points.shape[0], 3)
658 dxNorm = numpy.linalg.norm(dx)
659 # print(f"{time.perf_counter()-tic:.3f} seconds.")
661 if medianFilterEachIteration:
662 # use connectivity to filter
663 print("[globalCorrelation] [newton] Median filter of displacements...")
664 for nodeNumber in range(points.shape[0]):
665 # get rows of connectivity (i.e., tets) which include this point
666 connectedTets = numpy.where(connectivity == nodeNumber)[0]
667 neighbourPoints = numpy.unique(connectivity[connectedTets])
668 diff = numpy.median(displacements[neighbourPoints], axis=0) - displacements[nodeNumber]
669 displacements[nodeNumber] += 0.5 * diff
671 tic = time.perf_counter()
672 # print("[globalCorrelation] [newton] run labelTetrahedra: ", end="")
674 imTetLabel = spam.label.labelTetrahedra(im1.shape, points + displacements, connectivity, nThreads=nThreads)
675 # print(f"{time.perf_counter()-tic:.3f} seconds.")
677 tic = time.perf_counter()
678 # print("[globalCorrelation] [newton] run applyMeshTransformation: ", end="")
679 im1Def = spam.DIC.applyMeshTransformation(
680 im1,
681 points,
682 connectivity,
683 displacements,
684 imTetLabel=imTetLabel,
685 nThreads=nThreads,
686 )
687 # print(f"{time.perf_counter()-tic:.3f} seconds.")
689 if debugFiles:
690 tifffile.imwrite(f"{prefix}-def-i{i:03d}.tif", im1Def)
691 tifffile.imwrite(
692 f"{prefix}-residual-cropped-i{i:03d}.tif",
693 im1Def[meshPaddingSlice] - im2[meshPaddingSlice],
694 )
695 # tifffile.imwrite(f"{prefix}-imTetLabel-i{i:03d}.tif", imTetLabel)
697 pointData = {"displacements": displacements, "initialDisplacements": initialDisplacements, "fluctuations": numpy.subtract(displacements, initialDisplacements)}
699 # compute strain for each fields
700 cellData = {}
701 components = ["vol", "dev", "volss", "devss"]
702 for fieldName, field in pointData.items():
703 Ffield = spam.deformation.FfieldBagi(points, connectivity, field, verbose=False)
704 decomposedFfield = spam.deformation.decomposeFfield(Ffield, components)
705 for c in components:
706 cellData[f"{fieldName}-{c}"] = decomposedFfield[c]
708 spam.helpers.writeUnstructuredVTK(
709 points.copy(),
710 connectivity.copy(),
711 pointData=pointData,
712 cellData=cellData,
713 fileName=f"{prefix}-displacementFE-i{i:03d}.vtk",
714 )
716 # print("\t\[globalCorrelation] Error Out = %0.5f%%" % (error))
717 # reshapedDispl = displacements.reshape(points.shape[0], 3)
718 dMin = numpy.min(displacements, axis=0)
719 dMed = numpy.median(displacements, axis=0)
720 dMax = numpy.max(displacements, axis=0)
721 strMin = f"Min={dMin[0]: .3f} {dMin[1]: .3f} {dMin[2]: .3f}"
722 strMed = f"Med={dMed[0]: .3f} {dMed[1]: .3f} {dMed[2]: .3f}"
723 strMax = f"Max={dMax[0]: .3f} {dMax[1]: .3f} {dMax[2]: .3f}"
724 print(f"[globalCorrelation] [newton] i={i:03d}, displacements {strMin}, {strMed}, {strMax}, dx={dxNorm:.2f}")
726 return displacements