Coverage for src/pygnd/core.py: 92%
385 statements
« prev ^ index » next coverage.py v7.16.0, created at 2026-09-03 23:45 +0000
1import os
2import warnings
3from pathlib import Path
4from collections.abc import Iterable
6import numpy as np
7from tqdm import tqdm
8from scipy import optimize
9from joblib import Parallel, delayed
11from pygnd import io, rotations, quaternions
12from pygnd.utils import tqdm_joblib
15warnings.filterwarnings("ignore", category=UserWarning)
17_PRECISION = np.float32
20def _resolve_n_cpus(n_cpus: int) -> int:
21 """Resolve a joblib-style n_cpus value (e.g. -1 for all cores, -2 for all but
22 one) to a positive core count, for use in chunk-size arithmetic."""
23 cpu_count = os.cpu_count() or 1
24 if n_cpus < 0:
25 return max(1, cpu_count + 1 + n_cpus)
26 return max(1, n_cpus)
29def get_linear_operator(cs: int, slip_systems: str = "all") -> tuple[np.ndarray, np.ndarray]:
30 """Pre-calculate the A matrix for the given crystal structure and desired slip systems.
32 Args:
33 cs (int): The crystal structure of the material. 1 for FCC, 2 for BCC, 3 for HCP.
34 slip_systems (str, optional): The slip systems to be used. Defaults to 'all'.
35 (FCC) - unused, always 'all'
36 (BCC) - 'screw+110', 'screw+112', 'screw+123', 'screw+110+112', 'screw+110+123', 'screw+112+123', 'all'
37 (HCP) - 'basal', 'prismatic', 'pyramidal', 'basal+prismatic', 'basal+pyramidal', 'prismatic+pyramidal', 'all'
39 Returns:
40 A (np.ndarray): The A matrix for the given crystal structure and slip systems. Shape (9, n_slip_systems)
41 B (np.ndarray): The B matrix (psuedo-inverse of A) for the given crystal structure. Shape (n_slip_systems, 9)
42 """
43 # Check the input values
44 if type(cs) != int:
45 raise ValueError("Crystal structure must be an integer value.")
46 if cs not in [1, 2, 3]:
47 raise ValueError("Crystal structure must be 1, 2, or 3.")
48 if type(slip_systems) != str:
49 raise ValueError("Slip systems must be a string.")
50 slip_systems = slip_systems.lower().strip()
51 if slip_systems not in [
52 "all",
53 "screw+110",
54 "screw+112",
55 "screw+123",
56 "screw+110+112",
57 "screw+110+123",
58 "screw+112+123",
59 "basal",
60 "prismatic",
61 "pyramidal",
62 "basal+prismatic",
63 "basal+pyramidal",
64 "prismatic+pyramidal",
65 ]:
66 raise ValueError(
67 "Slip systems must be 'all', 'screw+110', 'screw+112', 'screw+123', 'screw+110+112', 'basal', 'basal+prismatic', depending on the crystam structure."
68 )
70 # Create the A matrix for the given crystal structure
71 if cs == 1:
72 a = np.sqrt(3).astype(_PRECISION) / 9
73 c = np.sqrt(3).astype(_PRECISION) / 84
74 d = 1 / 18
75 f = 3 / 14
77 # See Arsenlis & Parks 1999
78 B = np.array(
79 [
80 [a, 7 * c, -13 * c, -7 * c, -a, 13 * c, c, -c, 0],
81 [-a, 13 * c, -7 * c, -c, 0, c, 7 * c, -13 * c, a],
82 [0, c, -c, -13 * c, a, 7 * c, 13 * c, -7 * c, -a],
83 [a, -7 * c, 13 * c, 7 * c, -a, 13 * c, -c, -c, 0],
84 [-a, -13 * c, 7 * c, c, 0, c, -7 * c, -13 * c, a],
85 [0, -c, c, 13 * c, a, 7 * c, -13 * c, -7 * c, -a],
86 [a, -7 * c, -13 * c, 7 * c, -a, -13 * c, c, c, 0],
87 [-a, -13 * c, -7 * c, c, 0, -c, 7 * c, 13 * c, a],
88 [0, -c, -c, 13 * c, a, -7 * c, 13 * c, 7 * c, -a],
89 [a, 7 * c, 13 * c, -7 * c, -a, -13 * c, -c, c, 0],
90 [-a, 13 * c, 7 * c, -c, 0, -c, -7 * c, 13 * c, -a],
91 [0, c, c, -13 * c, a, -7 * c, -13 * c, 7 * c, -a],
92 [5 * d, f, 0, f, 5 * d, 0, 0, 0, -d],
93 [5 * d, 0, f, 0, -d, 0, f, 0, 5 * d],
94 [-d, 0, 0, 0, 5 * d, f, 0, f, 5 * d],
95 [5 * d, -f, 0, -f, 5 * d, 0, 0, 0, -d],
96 [5 * d, 0, -f, 0, -d, 0, -f, 0, 5 * d],
97 [-d, 0, 0, 0, 5 * d, -f, 0, -f, 5 * d],
98 ]
99 ).astype(_PRECISION)
101 # FCC
102 A = pseudo_inverse(B)
104 elif cs == 2:
105 # BCC
106 A = generate_BCC_A_matrix()
107 if slip_systems == "screw+110":
108 A = A[:, :16]
109 elif slip_systems == "screw+112":
110 A = np.hstack((A[:, :4], A[:, 16:28]))
111 elif slip_systems == "screw+123":
112 A = np.hstack((A[:, :4], A[:, 28:]))
113 elif slip_systems == "screw+110+112":
114 A = A[:, :28]
115 elif slip_systems == "screw+110+123":
116 A = np.hstack((A[:, :16], A[:, 28:]))
117 elif slip_systems == "screw+112+123":
118 A = np.hstack((A[:, :4], A[:, 16:]))
119 elif slip_systems == "all": 119 ↛ 121line 119 didn't jump to line 121 because the condition on line 119 was always true
120 pass
121 B = pseudo_inverse(A)
123 elif cs == 3: 123 ↛ 147line 123 didn't jump to line 147 because the condition on line 123 was always true
124 # HCP
125 A = generate_HCP_A_matrix()
126 if slip_systems == "basal":
127 A = A[:, :6] # 3 edge basal and 3 screw basal slip systems
128 elif slip_systems == "prismatic":
129 A = A[:, 6:9] # 3 edge prismatic slip systems
130 elif slip_systems == "pyramidal":
131 A = A[:, 9:] # 12 edge pyramidal and 12 screw pyramidal slip systems
132 elif slip_systems == "basal+prismatic":
133 A = A[:, :9]
134 elif slip_systems == "basal+pyramidal":
135 A = np.hstack((A[:, :6], A[:, 9:]))
136 elif slip_systems == "prismatic+pyramidal":
137 A = A[:, 6:]
138 elif slip_systems == "all": 138 ↛ 140line 138 didn't jump to line 140 because the condition on line 138 was always true
139 pass
140 B = pseudo_inverse(A)
142 # Add this diagnostic
143 # print(f"Condition number of A: {np.linalg.cond(A)}")
144 # print(f"Rank of A: {np.linalg.matrix_rank(A)}")
145 # print(f"||B @ A - I||_F: {np.linalg.norm(B @ A - np.eye(A.shape[1]), 'fro')}")
146 # print(f"||A @ B - I||_F: {np.linalg.norm(A @ B - np.eye(A.shape[0]), 'fro')}")
147 return (A, B)
150def generate_BCC_A_matrix() -> np.ndarray:
151 """Generate the A matrix for BCC crystal structure."""
152 # Burgers vectors and slip plane normals for BCC
153 b_n = np.array(
154 [
155 [[1, 1, -1], [1, 1, -1]], # <111> screw
156 [[1, -1, -1], [1, -1, -1]],
157 [[1, -1, 1], [1, -1, 1]],
158 [[1, 1, 1], [1, 1, 1]],
159 [[1, 1, -1], [0, 1, 1]], # {110}<111> edge
160 [[1, 1, -1], [1, 0, 1]],
161 [[1, 1, -1], [1, -1, 0]],
162 [[1, -1, -1], [0, 1, -1]],
163 [[1, -1, -1], [1, 0, 1]],
164 [[1, -1, -1], [1, 1, 0]],
165 [[1, -1, 1], [0, 1, 1]],
166 [[1, -1, 1], [1, 0, -1]],
167 [[1, -1, 1], [1, 1, 0]],
168 [[1, 1, 1], [0, 1, -1]],
169 [[1, 1, 1], [1, 0, -1]],
170 [[1, 1, 1], [1, -1, 0]],
171 [[-1, -1, 1], [-2, 1, -1]], # {112}<111> edge
172 [[-1, -1, 1], [1, -2, -1]],
173 [[-1, -1, 1], [1, 1, 2]],
174 [[-1, 1, 1], [-2, -1, -1]],
175 [[-1, 1, 1], [1, 2, -1]],
176 [[-1, 1, 1], [1, -1, 2]],
177 [[1, -1, 1], [2, 1, -1]],
178 [[1, -1, 1], [-1, -2, -1]],
179 [[1, -1, 1], [-1, 1, 2]],
180 [[1, 1, 1], [2, -1, -1]],
181 [[1, 1, 1], [-1, 2, -1]],
182 [[1, 1, 1], [-1, -1, 2]],
183 [[1, 1, -1], [1, 2, 3]], # {123}<111> edge
184 [[1, 1, -1], [-1, 3, 2]],
185 [[1, 1, -1], [2, 1, 3]],
186 [[1, 1, -1], [-2, 3, 1]],
187 [[1, 1, -1], [3, -1, 2]],
188 [[1, 1, -1], [3, -2, 1]],
189 [[1, -1, -1], [-1, 2, -3]],
190 [[1, -1, -1], [1, 3, -2]],
191 [[1, -1, -1], [2, -1, 3]],
192 [[1, -1, -1], [2, 3, -1]],
193 [[1, -1, -1], [3, 1, 2]],
194 [[1, -1, -1], [3, 2, 1]],
195 [[1, -1, 1], [1, -2, -3]],
196 [[1, -1, 1], [1, 3, 2]],
197 [[1, -1, 1], [2, -1, -3]],
198 [[1, -1, 1], [2, 3, 1]],
199 [[1, -1, 1], [3, 1, -2]],
200 [[1, -1, 1], [3, 2, -1]],
201 [[1, 1, 1], [1, 2, -3]],
202 [[1, 1, 1], [1, -3, 2]],
203 [[1, 1, 1], [2, 1, -3]],
204 [[1, 1, 1], [2, -3, 1]],
205 [[1, 1, 1], [-3, 1, 2]],
206 [[1, 1, 1], [-3, 2, 1]],
207 ]
208 ).astype(float)
209 burgers = b_n[:, 0] / np.sqrt(3)
210 normals = b_n[:, 1] / np.linalg.norm(b_n[:, 1], axis=1)[:, None]
212 # Get the sense vectors
213 t = np.cross(normals, burgers)
215 # Fix the screw dislocations (sense vectors are the burgers vectors)
216 t[:4] = burgers[:4]
218 # Calculate the outer product of the two vectors
219 outer = np.einsum("...i,...j->...ij", burgers, t)
221 # Convert to the (n_slip_systems, 9) matrix
222 A_bcc = outer.reshape(-1, 9).T
224 return A_bcc
227def generate_HCP_A_matrix() -> np.ndarray:
228 """Generate the A matrix for HCP crystal structure."""
229 # Relevant Direcitons in [uvtw] notation
230 b_n_uvtw = np.array(
231 [
232 [[1, 1, -2, 0], [0, 0, 0, 1]], # Basal
233 [[1, -2, 1, 0], [0, 0, 0, 1]],
234 [[-2, 1, 1, 0], [0, 0, 0, 1]],
235 [[2, -1, -1, 0], [0, 1, -1, 0]], # Prismatic
236 [[-1, 2, -1, 0], [1, 0, -1, 0]],
237 [[1, 1, -2, 0], [1, -1, 0, 0]],
238 [[-1, -1, 2, 3], [1, 0, -1, 1]], # Pyramidal
239 [[-2, 1, 1, 3], [1, 0, -1, 1]],
240 [[1, 1, -2, 3], [0, -1, 1, 1]],
241 [[-1, 2, -1, 3], [0, -1, 1, 1]],
242 [[2, -1, -1, 3], [-1, 1, 0, 1]],
243 [[1, -2, 1, 3], [-1, 1, 0, 1]],
244 [[2, -1, -1, 3], [-1, 0, 1, 1]],
245 [[1, 1, -2, 3], [-1, 0, 1, 1]],
246 [[-1, -1, 2, 3], [0, 1, -1, 1]],
247 [[1, -2, 1, 3], [0, 1, -1, 1]],
248 [[-2, 1, 1, 3], [1, -1, 0, 1]],
249 [[-1, 2, -1, 3], [1, -1, 0, 1]],
250 ]
251 ).astype(float)
253 # Convert to uvw
254 u = b_n_uvtw[:, 0, 0]
255 v = b_n_uvtw[:, 0, 1]
256 t = b_n_uvtw[:, 0, 2]
257 w = b_n_uvtw[:, 0, 3]
258 burgers = np.array([u - t, v - t, w]).T
259 burgers /= np.linalg.norm(burgers, axis=1)[:, None]
261 u = b_n_uvtw[:, 1, 0]
262 v = b_n_uvtw[:, 1, 1]
263 t = b_n_uvtw[:, 1, 2]
264 w = b_n_uvtw[:, 1, 3]
265 normals = np.array([u - t, v - t, w]).T
266 normals /= np.linalg.norm(normals, axis=1)[:, None]
268 # Get the sense vectors
269 t = np.cross(normals, burgers)
271 # Put in the screw dislocations
272 burgers = np.vstack((burgers[:3], burgers)) # 3 screw basal dislocations
273 t = np.vstack((burgers[:3], t))
274 burgers = np.vstack((burgers, burgers[-12:])) # 12 screw pyramidal dislocations
275 t = np.vstack((t, burgers[-12:]))
277 # Calculate the outer product of the two vectors
278 outer = np.einsum("...i,...j->...ij", burgers, t)
279 outer = outer.reshape(-1, 9).T
280 outer = outer / np.linalg.norm(outer, axis=0)
282 return outer
285def pseudo_inverse(A: np.ndarray) -> np.ndarray:
286 """Calculate the B matrix (psuedo-inverse of A) for the given A matrix.
288 Args:
289 A (np.ndarray): The A matrix. Shape (9, n_slip_systems)
291 Returns:
292 np.ndarray: The B matrix. Shape (n_slip_systems, 9)
293 """
294 # return A.T.dot(np.linalg.inv(A.dot(A.T))).astype(PRECISION)
295 return np.linalg.pinv(A).astype(_PRECISION)
298def get_completeness(grain_ids: np.ndarray) -> np.ndarray:
299 """
300 Vectorized version of neighborhood analysis for 3D EBSD dataset.
302 Args:
303 grain_ids (np.ndarray): The grain ID map. Shape (n_x, n_y, n_z)
305 Returns:
306 np.ndarray: The completeness array. Shape (n_x, n_y, n_z, 3)
307 """
308 shape = grain_ids.shape
310 # Initialize output array
311 completeness = np.zeros((*shape, 3), dtype=np.int8)
313 # Create masks for valid grain IDs
314 valid_grains = grain_ids != 0
316 # Compute transitions (now using broadcasting)
317 x_trans = np.ones(shape, dtype=bool)
318 y_trans = np.ones(shape, dtype=bool)
319 z_trans = np.ones(shape, dtype=bool)
321 if shape[0] > 1:
322 x_trans = np.pad(grain_ids[:-1, ...] != grain_ids[1:, ...], ((0, 1), (0, 0), (0, 0)))
323 if shape[1] > 1: 323 ↛ 325line 323 didn't jump to line 325 because the condition on line 323 was always true
324 y_trans = np.pad(grain_ids[:, :-1, :] != grain_ids[:, 1:, :], ((0, 0), (0, 1), (0, 0)))
325 if shape[2] > 1: 325 ↛ 329line 325 didn't jump to line 329 because the condition on line 325 was always true
326 z_trans = np.pad(grain_ids[..., :-1] != grain_ids[..., 1:], ((0, 0), (0, 0), (0, 1)))
328 # Interior points
329 interior_mask = np.zeros_like(grain_ids, dtype=bool)
330 interior_mask[1:-1, :, :] = True
332 # X-direction vectorized analysis
333 if shape[0] > 1:
334 completeness[0, :, :, 0] = np.where(
335 x_trans[0, :, :], 0, 1
336 ) # Forward differences for first slice
337 completeness[-1, :, :, 0] = np.where(
338 x_trans[-2, :, :], 0, 2
339 ) # Backward differences for last slice
340 completeness[1:-1, :, :, 0] = np.select( # Central differences
341 [
342 (x_trans[:-2, :, :] & x_trans[1:-1, :, :]),
343 x_trans[:-2, :, :],
344 x_trans[1:-1, :, :],
345 ],
346 [0, 1, 2],
347 default=3,
348 )
350 # Y-direction vectorized analysis
351 if shape[1] > 1: 351 ↛ 365line 351 didn't jump to line 365 because the condition on line 351 was always true
352 completeness[:, 0, :, 1] = np.where(y_trans[:, 0, :], 0, 1)
353 completeness[:, -1, :, 1] = np.where(y_trans[:, -2, :], 0, 2)
354 completeness[:, 1:-1, :, 1] = np.select(
355 [
356 (y_trans[:, :-2, :] & y_trans[:, 1:-1, :]),
357 y_trans[:, :-2, :],
358 y_trans[:, 1:-1, :],
359 ],
360 [0, 1, 2],
361 default=3,
362 )
364 # Z-direction (similar logic)
365 if shape[2] > 1: 365 ↛ 379line 365 didn't jump to line 379 because the condition on line 365 was always true
366 completeness[:, :, 0, 2] = np.where(z_trans[:, :, 0], 0, 1)
367 completeness[:, :, -1, 2] = np.where(z_trans[:, :, -2], 0, 2)
368 completeness[:, :, 1:-1, 2] = np.select(
369 [
370 (z_trans[:, :, :-2] & z_trans[:, :, 1:-1]),
371 z_trans[:, :, :-2],
372 z_trans[:, :, 1:-1],
373 ],
374 [0, 1, 2],
375 default=3,
376 )
378 # Zero out values where grain_id is 0
379 completeness[~valid_grains] = 0
381 return completeness
384def get_neighbors(completeness: np.ndarray) -> np.ndarray:
385 """Convert the completeness array into a list of neighbor indices."""
386 # Get the shape of the completeness array
387 shape = completeness.shape[:-1]
389 # Create coordinate shifts for the pairs and the scale for the finite difference calculation
390 shifts0 = np.zeros((3,) + shape + (3,), dtype=np.int32)
391 shifts1 = np.zeros((3,) + shape + (3,), dtype=np.int32)
392 scale = np.zeros(shape + (3,), dtype=_PRECISION)
394 # Only central and backward differences will have a shift in the first point
395 shifts0[0][(completeness[..., 0] == 2) | (completeness[..., 0] == 3)] = [-1, 0, 0]
396 shifts0[1][(completeness[..., 1] == 2) | (completeness[..., 1] == 3)] = [0, -1, 0]
397 shifts0[2][(completeness[..., 2] == 2) | (completeness[..., 2] == 3)] = [0, 0, -1]
399 # Only central and forward differences will have a shift in the second point
400 shifts1[0][(completeness[..., 0] == 1) | (completeness[..., 0] == 3)] = [1, 0, 0]
401 shifts1[1][(completeness[..., 1] == 1) | (completeness[..., 1] == 3)] = [0, 1, 0]
402 shifts1[2][(completeness[..., 2] == 1) | (completeness[..., 2] == 3)] = [0, 0, 1]
404 # Create coordinate array
405 coords = np.indices(shape).transpose(1, 2, 3, 0) # (x, y, z, ndim)
406 coords0 = np.stack([coords + shift for shift in shifts0], axis=-2)
407 coords1 = np.stack([coords + shift for shift in shifts1], axis=-2)
409 # Handle forward and backward differences (scale is 1)
410 scale[..., 0][(completeness[..., 0] == 1) | (completeness[..., 0] == 2)] = 1
411 scale[..., 1][(completeness[..., 1] == 1) | (completeness[..., 1] == 2)] = 1
412 scale[..., 2][(completeness[..., 2] == 1) | (completeness[..., 2] == 2)] = 1
414 # Handle central differences (scale is 2)
415 scale[..., 0][(completeness[..., 0] == 3)] = 2
416 scale[..., 1][(completeness[..., 1] == 3)] = 2
417 scale[..., 2][(completeness[..., 2] == 3)] = 2
419 # Make sure everywhere that has a 0 completeness has a 0 scale
420 scale[completeness == 0] = 0
422 # coords dimensions
423 # 0: x index in volume
424 # 1: y index in volume
425 # 2: z index in volume
426 # 3: the axis of the difference (0: x, 1: y, 2: z)
427 # 4: the coordinate of the voxel (0: x, 1: y, 2: z)
428 # so coords0[4, 5, 6, 0, 1] is the y-coordinate of the first voxel in the finite difference pair along the x-direction for the voxel at (4, 5, 6)
429 # and coords1[4, 5, 6, 0, 1] is the y-coordinate of the second voxel in the pair
431 return (coords0, coords1, scale)
434def get_finite_difference_coordinates(
435 grain_ids: np.ndarray,
436) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
437 """Calculate the coordinates for finite difference pairs in a 3D EBSD dataset.
439 Args:
440 grain_ids: 3D numpy array containing grain IDs
442 Returns:
443 tuple of three 3D arrays containing the coordinates of the first voxel, the coordinates of the second voxel, and the scale factor
444 """
445 # Get the completeness array
446 completeness = get_completeness(grain_ids)
448 # Get the neighbors
449 coords0, coords1, scale = get_neighbors(completeness)
451 return (coords0, coords1, scale)
454def get_orientation_gradients(
455 quats: np.ndarray,
456 pts0: np.ndarray,
457 pts1: np.ndarray,
458 distances: np.ndarray,
459 cs: int,
460 n_cpus: int = 1,
461 chunk_size: int = None,
462 progress_bar: bool = False,
463) -> np.ndarray:
464 """Calculate the orientation gradients for a 3D EBSD dataset.
465 This is essentially the rotation vectors corresponding to the disorientation between neighboring voxels,
466 divided by the spacing along each dimension. The result is a 3x3 matrix for each voxel.
467 This function will call the private function _get_orientation_gradients to do the actual calculations.
468 Supports parallel processing.
470 Args:
471 quats: 3D numpy array containing quaternions, (X, Y, Z, 4)
472 pts0: 3D numpy array containing the coordinates of the first voxel in the finite difference pairs, (X, Y, Z, 3)
473 pts1: 3D numpy array containing the coordinates of the second voxel in the finite difference pairs, (X, Y, Z, 3)
474 distances: 3D numpy array containing the distances between the finite difference pairs, (X, Y, Z, 3)
475 cs: The crystal structure of the material. 1 for FCC, 2 for BCC, 3 for HCP.
476 n_cpus: The number of CPUs to use for parallel processing. If None, all available CPUs minus one will be used.
478 Returns:
479 3D numpy array containing the orientation gradients, (X, Y, Z, 3, 3)
480 This is essentially 3 rotation vectors corresponding to the disorientation between neighboring voxels,
481 each divided by the distance between the finite difference pair. The 3x3 matrix for each point is the rotation vector for each axis.
482 """
483 # Get the shape
484 out_shape = quats.shape[:-1]
486 # Reshape the data to be 1D
487 quats = quats.reshape(-1, 4).astype(_PRECISION)
488 N = quats.shape[0]
489 pts0 = pts0.reshape(-1, 3, 3)
490 pts1 = pts1.reshape(-1, 3, 3)
491 distances = distances.reshape(-1, 3)
493 # Convert points to raveled indices
494 pts0 = np.stack(
495 [
496 np.ravel_multi_index(pts0[:, 0].T, out_shape),
497 np.ravel_multi_index(pts0[:, 1].T, out_shape),
498 np.ravel_multi_index(pts0[:, 2].T, out_shape),
499 ],
500 axis=-1,
501 )
502 pts1 = np.stack(
503 [
504 np.ravel_multi_index(pts1[:, 0].T, out_shape),
505 np.ravel_multi_index(pts1[:, 1].T, out_shape),
506 np.ravel_multi_index(pts1[:, 2].T, out_shape),
507 ],
508 axis=-1,
509 )
511 # Get quaternion pairs
512 q0 = np.stack(
513 [quats[pts0[:, 0]], quats[pts0[:, 1]], quats[pts0[:, 2]]],
514 axis=1,
515 dtype=_PRECISION,
516 ) # (n_pairs, 3, 4)
517 q1 = np.stack(
518 [quats[pts1[:, 0]], quats[pts1[:, 1]], quats[pts1[:, 2]]],
519 axis=1,
520 dtype=_PRECISION,
521 ) # (n_pairs, 3, 4)
522 del quats, pts0, pts1 # Free memory
524 # Get laue_id
525 laue_id = 11 if cs == 1 or cs == 2 else 9
527 if n_cpus == 1:
528 quats_disorientation = quaternions.qu_disorientation(q0, q1, laue_id, laue_id).transpose(
529 1, 0, 2
530 )
531 else:
532 # Setup chunk size
533 if chunk_size is None:
534 resolved_cpus = _resolve_n_cpus(n_cpus)
535 chunk_size = max(1, min(N // resolved_cpus, N // 100))
537 # Split the data into chunks
538 q0 = np.array_split(q0, max(1, q0.shape[0] // chunk_size))
539 q1 = np.array_split(q1, max(1, q1.shape[0] // chunk_size))
540 n_chunks = len(q0)
541 chunks = zip(q0, q1)
543 # Run the calculations in parallel
544 quats_disorientation = np.empty((N, 3, 4), dtype=_PRECISION)
545 if progress_bar: 545 ↛ 546line 545 didn't jump to line 546 because the condition on line 545 was never true
546 with tqdm_joblib(
547 tqdm(total=n_chunks, desc="Calculating orientation gradients")
548 ) as progress_bar:
549 out = Parallel(n_jobs=n_cpus, timeout=9999999)(
550 delayed(quaternions.qu_disorientation)(q0, q1, laue_id, laue_id)
551 for q0, q1 in chunks
552 )
553 else:
554 out = Parallel(n_jobs=n_cpus, timeout=9999999)(
555 delayed(quaternions.qu_disorientation)(q0, q1, laue_id, laue_id)
556 for q0, q1 in chunks
557 )
558 del q0, q1 # Free memory
560 # Concatenate the results
561 start_idx = 0
562 for chunk in tqdm(out, desc="Unpacking orientation gradients"):
563 end_idx = start_idx + chunk.shape[0]
564 quats_disorientation[start_idx:end_idx] = chunk
565 start_idx = end_idx
566 del out, chunk
567 quats_disorientation = quats_disorientation.transpose(1, 0, 2)
569 # Convert quaternions to rotation vectors
570 rot_vectors = quaternions.qu_log(quats_disorientation) * 2
571 del quats_disorientation # Free memory
573 # Get the misorientations from the rotation vectors
574 misorientation = np.linalg.norm(rot_vectors, axis=-1)
576 # Get the orientation gradients
577 with np.errstate(divide="ignore", invalid="ignore"):
578 gradient_tensors = np.where(
579 (misorientation == 0).reshape(3, -1, 1),
580 0,
581 rot_vectors / distances.T[..., None],
582 )
584 # Reshape the output
585 gradient_tensors = gradient_tensors.transpose(1, 0, 2).reshape(out_shape + (3, 3))
586 misorientation = misorientation.T.reshape(out_shape + (3,))
587 return gradient_tensors, misorientation
590def _minimize_l2(Lambda: np.ndarray, B: np.ndarray, chunk_size: int = None) -> np.ndarray:
591 """Perform the minimization using the L2 norm.
592 Negative densities indicate left handed dislocations, positive densities indicate right handed dislocations.
594 Args:
595 Lambda: The Nye tensor components. Shape (n_voxels, 9)
596 B: The B matrix. Shape (n_slip_systems, 9)
597 chunk_size: The size of the chunks to process in parallel. If None, the entire array is processed at once.
599 Returns:
600 np.ndarray: The dislocation density. Shape (n_slip_systems, n_voxels)"""
601 if chunk_size is None: 601 ↛ 606line 601 didn't jump to line 606 because the condition on line 601 was always true
602 dd = B.dot(Lambda.T).reshape((-1,))
603 else:
604 # Split Lambda into chunks, rebuilding the full (n_slip_systems, n_voxels)
605 # matrix before the final flatten so the ordering matches the unchunked path
606 chunks = np.array_split(Lambda, max(1, Lambda.shape[0] // chunk_size))
607 dd = np.hstack([B.dot(chunk.T) for chunk in chunks]).reshape((-1,))
608 return dd
611def _minimize_l1(Lambda: np.ndarray, A: np.ndarray, tol: float = 1e-4) -> np.ndarray:
612 """Perform L1 minimization using expanded basis (separate left/right handed dislocations).
614 Following Arsenlis & Parks 1999, the basis is expanded so that left-handed and
615 right-handed dislocations each get their own index, all constrained to be non-negative.
617 Args:
618 Lambda: The Nye tensor components. Shape (n_voxels, 9)
619 A: The A matrix for compact basis. Shape (9, n_slip_systems)
620 tol: Tolerance for constraint satisfaction
622 Returns:
623 np.ndarray: The dislocation density in compact basis (allowing negative values).
624 Shape (n_slip_systems, n_voxels)
625 """
626 Lambda = Lambda.astype(np.float64)
627 A = A.astype(np.float64)
629 n_constraints = A.shape[0]
630 n_slip_systems = A.shape[1]
631 N = Lambda.shape[0]
633 # Expand the basis: [A, -A] so that each column represents one sign
634 # rho_expanded = [rho_plus, rho_minus] where both are ≥ 0
635 # Then A @ (rho_plus - rho_minus) = Lambda
636 A_expanded = np.hstack([A, -A]) # Shape (9, 2*n_slip_systems)
637 n_variables = 2 * n_slip_systems
639 dd_compact = np.zeros((n_slip_systems, N), dtype=np.float64)
640 bad_count = 0
642 for i in range(N):
643 # Minimize sum of all densities (L1 norm of expanded representation)
644 c = np.ones(n_variables, dtype=np.float64)
646 # Equality constraint: A_expanded @ rho_expanded = Lambda[i]
647 A_eq = A_expanded
648 b_eq = Lambda[i].reshape(-1).astype(np.float64)
650 # All variables non-negative (this is the key for L1)
651 bounds = [(0.0, None)] * n_variables
653 # Solve
654 result = optimize.linprog(
655 c,
656 A_eq=A_eq,
657 b_eq=b_eq,
658 bounds=bounds,
659 method="highs",
660 options={
661 "presolve": True,
662 "disp": False,
663 "dual_feasibility_tolerance": tol,
664 "primal_feasibility_tolerance": tol,
665 },
666 )
668 if not result.success: 668 ↛ 669line 668 didn't jump to line 669 because the condition on line 668 was never true
669 dd_compact[:, i] = 0
670 bad_count += 1
671 else:
672 # Convert back to compact basis: rho = rho_plus - rho_minus
673 rho_plus = result.x[:n_slip_systems]
674 rho_minus = result.x[n_slip_systems:]
675 dd_compact[:, i] = rho_plus - rho_minus
677 if bad_count > 0: 677 ↛ 678line 677 didn't jump to line 678 because the condition on line 677 was never true
678 print(f"L1 minimization: {bad_count} failed optimizations out of {N}")
679 return dd_compact
682def minimize(
683 alpha: np.ndarray,
684 cs: int,
685 A: np.ndarray,
686 B: np.ndarray,
687 burgers: np.ndarray,
688 minimization="l2",
689 n_cpus: int = 1,
690 chunk_size: int = None,
691 progress_bar: bool = False,
692) -> np.ndarray:
693 """Minimize the dislocation density using the given minimization scheme.
694 The equation to be solved is A*rho = Lambda, where A is the A matrix, rho is the dislocation density, and Lambda is the Nye tensor.
695 This is solved directly using L2 minimization with the pseudo-inverse of A.
696 This can also be solved using L1 minimization, which is done in parallel for each point in the Nye tensor.
698 Args:
699 alpha: The Nye tensor components. Shape (..., 3, 3)
700 cs: The crystal structure of the material. 1 for FCC, 2 for BCC, 3 for HCP.
701 A: The A matrix. Shape (9, n_slip_systems)
702 B: The B matrix, psuedo-inverse of A. Shape (n_slip_systems, 9)
703 burgers: The burgers vector for the material.
704 For HCP with pyramidal slip systems, this should be a tuple of the two burgers vectors.
705 minimization: The minimization scheme to use. Either 'l2' or 'l1'.
706 n_cpus: The number of CPUs to use for parallel processing during L1 minimization.
707 If None, all available CPUs minus one will be used. Not used for L2 minimization.
708 chunk_size: The size of the chunks to process in parallel during L1 minimization.
709 Default chunk size is max(1, n_voxels / (4 * n_cpus)).
710 progress_bar: Whether to display a progress bar during L1 minimization."""
711 # Equation to be solved -> A*rho[array form] = Lambda[Nye in array form]
712 # Solve: A*rho = Lambd
713 # Nye tensor must be converted into array form Lambda
714 # Get shape
715 shape = alpha.shape[:-2]
716 Lambda = alpha.reshape(-1, 9)
717 out_shape = (A.shape[1],) + shape
719 # DEBUG
720 # Lambda_i = Lambda[0]
721 # print(f"Debug: Lambda[0] = {Lambda_i}")
722 # print(f"Debug: A shape = {A.shape}, B shape = {B.shape}")
723 # print(f"A matrix:\n{A}")
724 # print(f"B matrix:\n{B}")
725 # print(f"Debug: A @ B (should be close to identity):\n{A @ B}")
726 # print(
727 # f"Debug: L1 minimization result for first point:\n{_minimize_l1(Lambda[0:1], A)}"
728 # )
729 # print(
730 # f"Debug: L2 minimization result for first point:\n{_minimize_l2(Lambda[0:1], B)}"
731 # )
732 # exit()
734 if minimization == "l2":
735 print("Performing L2 minimization...")
736 dd = _minimize_l2(Lambda, B, chunk_size).reshape(out_shape)
738 elif minimization == "l1":
739 # Setup chunk size
740 if chunk_size is None:
741 resolved_cpus = _resolve_n_cpus(n_cpus)
742 chunk_size = max(1, Lambda.shape[0] // (resolved_cpus * 4))
744 # Split into chunks
745 chunks = np.array_split(Lambda, max(1, Lambda.shape[0] // chunk_size))
747 # Add progress bar if desired
748 if progress_bar: 748 ↛ 749line 748 didn't jump to line 749 because the condition on line 748 was never true
749 chunks = tqdm(chunks, desc="Minimizing (L1) chunks")
751 # Process chunks in parallel
752 with Parallel(n_jobs=n_cpus, timeout=9999999) as parallel:
753 chunk_results = parallel(delayed(_minimize_l1)(chunk, A) for chunk in chunks)
755 # Combine the results
756 dd = np.hstack(chunk_results).reshape(out_shape)
758 else:
759 raise ValueError("Minimization scheme not recognized. Please choose either 'l1' or 'l2'")
761 # Divide by Burgers vector correctly based on crystal structure and slip systems
762 if cs == 1 or cs == 2:
763 dd = dd / burgers
764 else:
765 burgers_basal_prismatic = None
766 burgers_pyramidal = None
767 if not isinstance(burgers, (tuple, list, np.ndarray)):
768 burgers = (burgers,)
769 if len(burgers) == 2:
770 burgers_basal_prismatic = burgers[0]
771 burgers_pyramidal = burgers[1]
772 elif len(burgers) == 1 and dd.shape[0] <= 9:
773 burgers_basal_prismatic = burgers[0]
774 elif len(burgers) == 1 and dd.shape[0] == 24: 774 ↛ 777line 774 didn't jump to line 777 because the condition on line 774 was always true
775 burgers_pyramidal = burgers[0]
776 else:
777 raise ValueError(
778 "For HCP, when mixing basal/prismatic and pyramidal slip systems, the Burgers vector must be a tuple of (basa/prismatic, pyramidal) Burgers vectors."
779 )
781 # Basal slip
782 if dd.shape[0] == 6:
783 dd = dd / burgers_basal_prismatic
785 # Prismatic slip
786 elif dd.shape[0] == 3:
787 dd = dd / burgers_basal_prismatic
789 # Pyramidal slip
790 elif dd.shape[0] == 24:
791 dd = dd / burgers_pyramidal
793 # Basal + Prismatic slip
794 elif dd.shape[0] == 9:
795 dd = dd / burgers_basal_prismatic
797 # Basal + Pyramidal slip
798 elif dd.shape[0] == 30:
799 dd[:6] = dd[:6] / burgers_basal_prismatic
800 dd[6:] = dd[6:] / burgers_pyramidal
802 # Prismatic + Pyramidal slip
803 elif dd.shape[0] == 27:
804 dd[:3] = dd[:3] / burgers_basal_prismatic
805 dd[3:] = dd[3:] / burgers_pyramidal
807 # All slip systems
808 elif dd.shape[0] == 33: 808 ↛ 812line 808 didn't jump to line 812 because the condition on line 808 was always true
809 dd[:9] = dd[:9] / burgers_basal_prismatic
810 dd[9:] = dd[9:] / burgers_pyramidal
812 return dd
815def calculate(
816 euler: np.ndarray,
817 ids: np.ndarray,
818 cs: int,
819 slip_systems: str,
820 burgers: tuple[float, tuple],
821 spacing: tuple,
822 minimization: tuple[str, tuple] = "l2",
823 n_cpus: int = -1,
824 progress_bar: bool = True,
825 chunk_size: int = None,
826) -> tuple[np.ndarray, np.ndarray]:
827 """Calculate the GND density for a 2D or 3D EBSD dataset.
829 Args:
830 euler: The Euler angles for the dataset. Shape (..., 3)
831 ids: The grain IDs for the dataset. Shape (...)
832 cs: The crystal structure of the material. 1 for FCC, 2 for BCC, 3 for HCP.
833 slip_systems (str, optional): The slip systems to be used. Defaults to 'all'.
834 (FCC) - unused, always 'all'
835 (BCC) - 'screw+110', 'screw+112', 'screw+123', 'screw+110+112', 'screw+110+123', 'screw+112+123', 'all'
836 (HCP) - 'basal', 'prismatic', 'pyramidal', 'basal+prismatic', 'basal+pyramidal', 'prismatic+pyramidal', 'all'
837 burgers: The burgers vector for the material in meters.
838 For HCP with pyramidal slip systems, this should be a tuple of the two burgers vectors.
839 spacing: The spacing between voxels in meters. Needs to be a tuple of the same length as the ids.
840 minimization: The minimization scheme to use. Either 'l2' or 'l1' or a tuple containing both.
841 n_cpus: The number of CPUs to use for parallel processing during L1 minimization.
842 If None, all available CPUs minus one will be used. Not used for L2 minimization.
843 progress_bar: Whether to display a progress bar during L1 minimization.
845 Returns:
846 dd: The dislocation density. Shape (n_slip_systems, ...)
847 mis: The misorientation. Shape (3, ...)"""
849 # Check inputs
850 ndim = euler.ndim - 1
851 if euler.shape[-1] != 3:
852 raise ValueError("The Euler angles must have shape (..., 3)")
853 if ids.shape != euler.shape[:-1]:
854 raise ValueError("The grain IDs must have the same shape as the Euler angles")
855 if cs not in (1, 2, 3):
856 raise ValueError("The crystal structure must be 1 for FCC, 2 for BCC, or 3 for HCP")
857 if len(spacing) != ndim:
858 raise ValueError("The spacing must have the same number of dimensions as the Euler angles")
860 euler = euler.astype(_PRECISION)
862 # Handle minimization
863 if isinstance(minimization, tuple):
864 if len(minimization) > 2: 864 ↛ 865line 864 didn't jump to line 865 because the condition on line 864 was never true
865 raise ValueError(
866 "The minimization scheme must be either 'l1' or 'l2' or both, but cannot have more than two elements"
867 )
868 minimization = tuple(m.lower() for m in minimization)
869 elif isinstance(minimization, str):
870 minimization = (minimization.lower(),)
871 minimization = sorted(minimization)
872 for m in minimization:
873 if m not in ("l1", "l2"): 873 ↛ 874line 873 didn't jump to line 874 because the condition on line 873 was never true
874 raise ValueError("The minimization scheme must be either 'l1' or 'l2'")
876 # Get the linear operator
877 A, B = get_linear_operator(cs, slip_systems)
879 # Convert Euler angles to quaternions
880 quats = rotations.eu2qu(euler)
881 del euler # Free memory
883 # Get the finite difference coordinates
884 print("Getting finite difference coordinates...")
885 nbrs0, nbrs1, distances = get_finite_difference_coordinates(ids)
886 distances *= spacing
888 # Get the orientation gradients
889 dphi, mis = get_orientation_gradients(
890 quats,
891 nbrs0,
892 nbrs1,
893 distances,
894 cs,
895 n_cpus,
896 progress_bar=progress_bar,
897 chunk_size=chunk_size,
898 )
899 del quats, nbrs0, nbrs1, distances # Free memory
900 mis = np.rad2deg(mis)
901 mis = mis.transpose(3, 0, 1, 2) # (..., 3) -> (3, ...)
903 # Calculate the alpha tensor
904 trace = np.trace(dphi, axis1=3, axis2=4)
905 alpha = dphi.transpose(0, 1, 2, 4, 3) - trace[..., None, None] * np.eye(3).reshape(
906 1, 1, 1, 3, 3
907 )
908 del dphi, trace # Free memory
910 # Minimize the dislocation density
911 dd = {}
912 for m in minimization:
913 dd[m] = np.abs(minimize(alpha, cs, A, B, burgers, m, n_cpus, progress_bar=progress_bar))
915 return dd, mis
918def _run_calculation(
919 euler: np.ndarray,
920 ids: np.ndarray,
921 spacing: tuple,
922 cs: int,
923 burgers: np.ndarray,
924 minimization: str | Iterable[str],
925 n_cpus: int,
926 slip_systems: str,
927 progress_bar: bool,
928 chunk_size: int,
929) -> tuple[dict[str, np.ndarray], np.ndarray] | None:
930 """Shared calculation + results-summary step used by both
931 calculate_and_save_dream3d and calculate_and_save_ang.
933 Args:
934 euler: The Euler angles for the dataset. Shape (..., 3)
935 ids: The grain IDs for the dataset. Shape (...)
936 spacing: The spacing between voxels in meters.
937 cs: The crystal structure of the material. 1 for FCC, 2 for BCC, 3 for HCP.
938 burgers: The burgers vector for the material in meters.
939 minimization: The minimization scheme to use. Either 'l2' or 'l1' or both.
940 n_cpus: The number of CPUs to use for parallel processing during L1 minimization.
941 slip_systems: The slip systems to be used.
942 progress_bar: Whether to display a progress bar during L1 minimization.
943 chunk_size: The size of the chunks to process in parallel.
945 Returns:
946 `(dd, mis)` on success, or `None` if the calculation itself raised.
947 """
948 print("Euler angles shape:", euler.shape)
949 print("Grain IDs shape:", ids.shape)
950 print("Voxel spacing (m):", spacing)
952 try:
953 dd, mis = calculate(
954 euler,
955 ids,
956 cs,
957 slip_systems,
958 burgers,
959 spacing,
960 minimization,
961 n_cpus,
962 progress_bar,
963 chunk_size,
964 )
965 except Exception as e:
966 print(f"Error during GND calculation: {e}")
967 return None
969 print("\nResults summary:")
970 print("----------------")
971 for m in dd:
972 print(f"- {m} GND max: {dd[m].max():.3e} m\u207b\u00b2")
973 nonzero = dd[m][dd[m] > 0]
974 if nonzero.size > 0:
975 print(f"- {m} GND min (non-zero): {nonzero.min():.3e} m\u207b\u00b2")
976 else:
977 print(f"- {m} GND min (non-zero): none (all values are zero)")
978 print(f"- FDM_avg max: {mis.mean(axis=0).max():.3f}\u00b0")
979 print(f"- FDM_max max: {mis.max(axis=0).max():.3f}\u00b0")
980 print("----------------")
981 return dd, mis
984def calculate_and_save_dream3d(
985 dream3d_path: Path | str,
986 ids_name: str,
987 euler_name: str,
988 cs: int,
989 burgers: np.ndarray,
990 spacing_units: str = "um",
991 minimization: str | Iterable[str] = "l2",
992 n_cpus: int = -1,
993 slip_systems: str = "all",
994 progress_bar: bool = False,
995 chunk_size: int = 1000,
996) -> bool:
997 """Calculate the GND density for a DREAM3D EBSD dataset and save the
998 results back into the same DREAM3D file.
1000 Args:
1001 dream3d_path: Path to the DREAM3D file.
1002 ids_name: Name of the grain ID data array in the DREAM3D file.
1003 euler_name: Name of the Euler angles data array in the DREAM3D file.
1004 cs: The crystal structure of the material. 1 for FCC, 2 for BCC, 3 for HCP.
1005 burgers: The burgers vector for the material in meters. For HCP with
1006 pyramidal slip systems, this should be a tuple of the two burgers vectors.
1007 spacing_units: Units of the spacing in the DREAM3D file. Options are 'nm',
1008 'nanometer', 'nanometers', 'um', 'µm', 'micron', 'microns', 'micrometer',
1009 'micrometers', 'mm', 'millimeter', 'millimeters', 'm', 'meter', 'meters'.
1010 minimization: The minimization scheme to use. Either 'l2' or 'l1' or both.
1011 n_cpus: The number of CPUs to use for parallel processing during L1
1012 minimization. If -1, all available CPUs minus one will be used. Not
1013 used for L2 minimization.
1014 slip_systems: The slip systems to be used. Defaults to 'all'.
1015 (FCC) - unused, always 'all'.
1016 (BCC) - 'screw+110', 'screw+112', 'screw+123', 'screw+110+112',
1017 'screw+110+123', 'screw+112+123', 'all'.
1018 (HCP) - 'basal', 'prismatic', 'pyramidal', 'basal+prismatic',
1019 'basal+pyramidal', 'prismatic+pyramidal', 'all'.
1020 progress_bar: Whether to display a progress bar during L1 minimization.
1021 chunk_size: The size of the chunks to process in parallel.
1023 Returns:
1024 True if the results were saved back into the DREAM3D file. False if the
1025 calculation or save failed, in which case a `.npy` fallback is written
1026 next to the DREAM3D file instead.
1027 """
1028 folder = Path(dream3d_path).parent
1029 euler, ids, spacing = io.read_dream3d(dream3d_path, ids_name, euler_name, spacing_units)
1031 result = _run_calculation(
1032 euler, ids, spacing, cs, burgers, minimization, n_cpus, slip_systems, progress_bar, chunk_size
1033 )
1034 if result is None: 1034 ↛ 1035line 1034 didn't jump to line 1035 because the condition on line 1034 was never true
1035 return False
1036 dd, mis = result
1038 try:
1039 saved = io.save_to_dream3d(dream3d_path, ids_name, dd, mis)
1040 except Exception as e:
1041 print(f"Error saving data to DREAM3D file: {e}")
1042 saved = False
1044 if not saved: 1044 ↛ 1045line 1044 didn't jump to line 1045 because the condition on line 1044 was never true
1045 io.save_npz(dd, mis, folder)
1046 print("Failed to save results to DREAM3D file. See .npy files for raw data.")
1047 return False
1049 print(f"Results saved to DREAM3D file: {Path(dream3d_path).absolute()}")
1050 return True
1053def calculate_and_save_ang(
1054 ang_path: Path | str,
1055 cs: int,
1056 burgers: np.ndarray,
1057 grain_ids_path: Path | str | None = None,
1058 minimization: str | Iterable[str] = "l2",
1059 n_cpus: int = -1,
1060 slip_systems: str = "all",
1061 progress_bar: bool = False,
1062 chunk_size: int = 1000,
1063) -> bool:
1064 """Calculate the GND density for an .ang EBSD dataset and save the results
1065 as `.npy` files (and preview images) next to the .ang file.
1067 Args:
1068 ang_path: Path to the .ang file.
1069 cs: The crystal structure of the material. 1 for FCC, 2 for BCC, 3 for HCP.
1070 burgers: The burgers vector for the material in meters. For HCP with
1071 pyramidal slip systems, this should be a tuple of the two burgers vectors.
1072 grain_ids_path: Path to a grain-data file generated by OIM Analysis. If
1073 not provided, the entire dataset is assumed to be a single grain.
1074 minimization: The minimization scheme to use. Either 'l2' or 'l1' or both.
1075 n_cpus: The number of CPUs to use for parallel processing during L1
1076 minimization. If -1, all available CPUs minus one will be used. Not
1077 used for L2 minimization.
1078 slip_systems: The slip systems to be used. Defaults to 'all'.
1079 (FCC) - unused, always 'all'.
1080 (BCC) - 'screw+110', 'screw+112', 'screw+123', 'screw+110+112',
1081 'screw+110+123', 'screw+112+123', 'all'.
1082 (HCP) - 'basal', 'prismatic', 'pyramidal', 'basal+prismatic',
1083 'basal+pyramidal', 'prismatic+pyramidal', 'all'.
1084 progress_bar: Whether to display a progress bar during L1 minimization.
1085 chunk_size: The size of the chunks to process in parallel.
1087 Returns:
1088 True if the calculation succeeded and results were saved, False otherwise.
1089 """
1090 if grain_ids_path is None: 1090 ↛ 1096line 1090 didn't jump to line 1096 because the condition on line 1090 was always true
1091 warnings.warn(
1092 "grain_ids_path not provided. Assuming the entire dataset is a single grain.",
1093 Warning,
1094 )
1096 folder = Path(ang_path).parent
1097 euler, ids, spacing = io.read_ang(ang_path, grain_ids_path)
1099 result = _run_calculation(
1100 euler, ids, spacing, cs, burgers, minimization, n_cpus, slip_systems, progress_bar, chunk_size
1101 )
1102 if result is None: 1102 ↛ 1103line 1102 didn't jump to line 1103 because the condition on line 1102 was never true
1103 return False
1104 dd, mis = result
1106 # Remove z dimension if 2D, order is (C, Z, Y, X), with C either being 3 (mis) or n_slip_systems (dd)
1107 dd = {k: v[:, 0] for k, v in dd.items()}
1108 mis = mis[:, 0]
1109 io.save_npz(dd, mis, folder)
1110 print("Calculation complete. Check .npy files for raw data.")
1111 io.generate_images(gnd_data=dd, fdm_data=mis, folder=folder)
1112 return True
1115def calculate_and_save(
1116 cs: int,
1117 burgers: np.ndarray,
1118 dream3d_path: Path | str | None = None,
1119 ids_name: str | None = None,
1120 euler_name: str | None = None,
1121 ang_path: Path | str | None = None,
1122 grain_ids_path: Path | str | None = None,
1123 spacing_units: str = "um",
1124 minimization: str | Iterable[str] = "l2",
1125 n_cpus: int = -1,
1126 slip_systems: str = "all",
1127 progress_bar: bool = False,
1128 chunk_size: int = 1000,
1129) -> bool:
1130 """Deprecated: use `calculate_and_save_dream3d` or `calculate_and_save_ang` instead.
1132 This combined-argument entry point is kept only for backwards compatibility
1133 and dispatches to one of the two functions above based on whether
1134 `dream3d_path` or `ang_path` is provided.
1135 """
1136 warnings.warn(
1137 "calculate_and_save() is deprecated; use calculate_and_save_dream3d() "
1138 "or calculate_and_save_ang() instead.",
1139 DeprecationWarning,
1140 stacklevel=2,
1141 )
1142 if dream3d_path is not None and ang_path is not None:
1143 raise ValueError("Provide either dream3d_path or ang_path, not both.")
1144 if dream3d_path is not None:
1145 if ids_name is None or euler_name is None:
1146 raise ValueError(
1147 "If dream3d_path is provided, ids_name and euler_name must also be provided."
1148 )
1149 return calculate_and_save_dream3d(
1150 dream3d_path,
1151 ids_name,
1152 euler_name,
1153 cs,
1154 burgers,
1155 spacing_units,
1156 minimization,
1157 n_cpus,
1158 slip_systems,
1159 progress_bar,
1160 chunk_size,
1161 )
1162 if ang_path is not None:
1163 return calculate_and_save_ang(
1164 ang_path,
1165 cs,
1166 burgers,
1167 grain_ids_path,
1168 minimization,
1169 n_cpus,
1170 slip_systems,
1171 progress_bar,
1172 chunk_size,
1173 )
1174 raise ValueError("Either dream3d_path or ang_path must be provided.")