Coverage for src/pygnd/io.py: 97%
241 statements
« prev ^ index » next coverage.py v7.16.0, created at 2026-09-03 23:45 +0000
1import os
2from pathlib import Path
4from PIL import Image
5import h5py
6import numpy as np
9# Define Dream3d data types
10_DREAM3D_DTYPES = {
11 np.uint8: "DataArray<uint8_t> ",
12 np.int8: "DataArray<int8_t> ",
13 np.uint16: "DataArray<uint16_t> ",
14 np.int16: "DataArray<int16_t> ",
15 np.uint32: "DataArray<uint32_t> ",
16 np.int32: "DataArray<int32_t> ",
17 np.uint64: "DataArray<uint64_t> ",
18 np.int64: "DataArray<int64_t> ",
19 np.float32: "DataArray<float> ",
20 np.float64: "DataArray<double> ",
21 bool: "DataArray<bool> ",
22}
23_XDMF_DTYPE_FORMATS = { # (NumberType, Precision)
24 np.uint8: ("UChar", "1"),
25 np.int8: ("Char", "1"),
26 np.uint16: ("UInt", "2"),
27 np.int16: ("Int", "2"),
28 np.uint32: ("UInt", "4"),
29 np.int32: ("Int", "4"),
30 np.uint64: ("UInt", "8"),
31 np.int64: ("Int", "8"),
32 np.float32: ("Float", "4"),
33 np.float64: ("Float", "8"),
34 bool: ("uchar", "1"),
35}
38def read_ang(
39 path: str | Path, ids_path: str | Path | None = None
40) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
41 """
42 Reads an .ang file into a numpy array
43 If ids_path is provided, reads the Feature IDs from that file as well.
44 If ids_path is None, all Feature IDs are set to 1.
45 The ids_path should be a path to a grain data file generated by OIM Analysis.
46 If not generated by OIM, the Feature IDs should be in the 9th column (index 8) of the file and be a integer grain ID for each point of the .ang file.
48 Parameters:
49 path (str): Path to the .ang file.
50 ids_path (Optional[str]): Path to the file containing Feature IDs.
51 Returns:
52 Tuple[np.ndarray, np.ndarray, np.ndarray]: Euler angles, Feature IDs, spacing.
53 """
54 num_header_lines = 0
55 col_names = None
56 ncols = nrows = res = None
57 with open(path, "r", encoding="utf-8") as f:
58 for line in f: 58 ↛ 71line 58 didn't jump to line 71
59 if line[0] == "#":
60 num_header_lines += 1
61 if "NCOLS_ODD" in line:
62 ncols = int(line.split(": ")[1].strip())
63 elif "NROWS" in line:
64 nrows = int(line.split(": ")[1].strip())
65 elif "COLUMN_HEADERS" in line:
66 col_names = line.split(": ")[1].strip().split(", ")
67 elif "XSTEP" in line:
68 res = float(line.split(": ")[1].strip())
69 else:
70 break
71 if ncols is None or nrows is None or res is None:
72 raise ValueError(
73 "The .ang header is missing one of the required NCOLS_ODD, NROWS, or XSTEP fields."
74 )
75 raw_data = np.genfromtxt(path, skip_header=num_header_lines)
76 n_entries = raw_data.shape[-1]
77 if col_names is None:
78 default_names = ["phi1", "PHI", "phi2", "x", "y", "IQ", "CI", "Phase index"]
79 col_names = default_names[:n_entries] + [
80 f"col_{i}" for i in range(len(default_names), n_entries)
81 ]
82 if raw_data.shape[0] != ncols * nrows:
83 raise ValueError(
84 f"The number of data points ({raw_data.shape[0]}) does not match the expected grid ({nrows} rows, {ncols} cols, {ncols * nrows} total points). "
85 )
86 data = raw_data.reshape((nrows, ncols, n_entries))
88 out = {col_names[i]: data[:, :, i] for i in range(n_entries)}
89 eulerangles = np.array([out["phi1"], out["PHI"], out["phi2"]]).T.astype(float)
90 eulerangles = eulerangles.reshape(1, *eulerangles.shape).transpose(0, 2, 1, 3)
91 if ids_path is not None:
92 grain_data = np.genfromtxt(ids_path, dtype=float, comments="#")
93 ids = grain_data[:, 8].reshape(eulerangles.shape[:-1]).astype(int)
94 else:
95 ids = np.ones(eulerangles.shape[:-1], dtype=int)
96 spacing = np.array([res, res, res]) * 1e-6 # Convert from microns to meters
97 return eulerangles, ids, spacing
100def read_dream3d(
101 path: str | Path,
102 ids_name: str = "FeatureIds",
103 euler_name: str = "EulerAngles",
104 spacing_units: str = "microns",
105) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
106 """
107 Reads Euler angles, Feature IDs, and spacing from a DREAM3D file.
109 Parameters:
110 path (str): Path to the DREAM3D file.
111 ids_name (str): Name of the DataArray containing the Feature IDs.
112 euler_name (str): Name of the DataArray containing the Euler angles.
113 spacing_units (str): Units of the spacing in the DREAM3D file. Options are
114 'nm', 'nanometer', 'nanometers', 'um', 'µm', 'micron', 'microns', 'micrometer',
115 'micrometers', 'mm', 'millimeter', 'millimeters', 'm', 'meter', 'meters'.
116 Returns:
117 Tuple[np.ndarray, np.ndarray, np.ndarray]: Euler angles, Feature IDs, spacing.
118 """
119 if not os.path.exists(path):
120 raise FileNotFoundError(f"DREAM3D file not found at path: {Path(path).absolute()}")
122 ids = extract_data_from_h5(path, ids_name)
123 if ids is None:
124 raise KeyError(
125 f"Could not find a data array with the name '{ids_name}' in the dream3d file."
126 )
127 ids = ids[..., 0] if ids.ndim == 4 and ids.shape[-1] == 1 else ids
129 eulerangles = extract_data_from_h5(path, euler_name)
130 if eulerangles is None:
131 raise KeyError(
132 f"Could not find a data array with the name '{euler_name}' in the dream3d file."
133 )
135 spacing = read_dream3d_spacing(path, spacing_units=spacing_units)
137 return eulerangles, ids, spacing
140def read_dream3d_spacing(path: str, spacing_units: str = "microns") -> np.ndarray:
141 """
142 Reads the spacing information from a DREAM3D file and converts it to meters.
144 Parameters:
145 path (str): Path to the DREAM3D file.
146 spacing_units (str): Units of the spacing in the DREAM3D file. Options are
147 'nm', 'nanometer', 'nanometers', 'um', 'µm', 'micron', 'microns', 'micrometer',
148 'micrometers', 'mm', 'millimeter', 'millimeters', 'm', 'meter', 'meters'.
149 Returns:
150 np.ndarray: Spacing in meters.
151 """
152 if not os.path.exists(path):
153 raise FileNotFoundError(f"DREAM3D file not found at path: {path}")
155 spacing = extract_data_from_h5(path, "SPACING")
157 # If the spacing path isn't found, it is likely a dream3dnx file and spacing is an attribute
158 if spacing is None:
159 spacing = extract_attribute_from_h5(path, "_SPACING")
160 if spacing is None: 160 ↛ 166line 160 didn't jump to line 166 because the condition on line 160 was always true
161 raise ValueError(
162 "Could not find spacing information in the dream3d file. Attempted to find a 'SPACING' dataset (DREAM3D format) and a '_SPACING' attribute (DREAM3DNX format)."
163 )
165 # Convert spacing to meters
166 if spacing_units == "nm" or spacing_units == "nanometer" or spacing_units == "nanometers":
167 spacing *= 1e-9
168 elif (
169 spacing_units == "um"
170 or spacing_units == "micron"
171 or spacing_units == "microns"
172 or spacing_units == "micrometer"
173 or spacing_units == "micrometers"
174 or spacing_units == "µm"
175 ):
176 spacing *= 1e-6
177 elif spacing_units == "mm" or spacing_units == "millimeter" or spacing_units == "millimeters":
178 spacing *= 1e-3
179 elif spacing_units == "m" or spacing_units == "meter" or spacing_units == "meters":
180 pass
181 else:
182 raise ValueError(
183 "units must be one of 'nm', 'nanometer', 'um', 'µm', 'micron', 'microns', 'micrometer', 'micrometers', 'mm', 'millimeter', 'millimeters', 'm', or 'meters'."
184 )
186 return spacing
189def add_dataset_to_h5(h5group: h5py.Group, name: str, data: np.ndarray) -> h5py.Dataset:
190 """
191 Adds a new dataset to an existing HDF5 group. Designed to be used with DREAM3D files.
193 Parameters:
194 h5group (h5py.Group): The HDF5 group to add the dataset to.
195 name (str): The name of the new dataset.
196 data (np.ndarray): The data to be stored in the dataset.
197 Returns:
198 h5py.Dataset: The created dataset.
199 """
200 # Check that h5group is indeed an h5py.Group
201 if not isinstance(h5group, h5py.Group):
202 raise TypeError("h5group must be an instance of h5py.Group")
204 # Check that the data is the same shape as other datasets in the group
205 for key in h5group.keys():
206 if isinstance(h5group[key], h5py.Dataset): 206 ↛ 205line 206 didn't jump to line 205 because the condition on line 206 was always true
207 if h5group[key].shape[:-1] != data.shape[:-1]:
208 raise ValueError(
209 f"Data shape {data.shape} does not match existing dataset shape {h5group[key].shape} in the group."
210 )
211 break
213 # Check to see if the dataset already exists, if so just overwrite it
214 if name in h5group:
215 print(f"Dataset '{name}' already exists in HDF5 group. Overwriting.")
216 h5group[name][...] = data
217 return h5group[name]
219 dtype = data.dtype.type
220 if dtype not in _DREAM3D_DTYPES:
221 raise TypeError(f"Unsupported data type for DREAM3D: {dtype}")
222 dset = h5group.create_dataset(name, data=data, dtype=dtype)
223 dset.attrs["ComponentDimensions"] = np.uint64([data.shape[-1]])
224 dset.attrs["Tuple Axis Dimensions"] = np.bytes_(
225 f"x={str(data.shape[2])},y={str(data.shape[1])},z={str(data.shape[0])} "
226 )
227 dset.attrs["DataArrayVersion"] = np.int32([2])
228 dset.attrs["ObjectType"] = np.bytes_(_DREAM3D_DTYPES[dtype])
229 dset.attrs["TupleDimensions"] = np.uint64(np.squeeze(data.shape[:-1][::-1]))
230 print(f"Added dataset '{name}' to HDF5 group.")
232 return dset
235def add_dataset_to_xdmf(xdmf_path: str | Path, dataset_name: str, data_array: np.ndarray) -> None:
236 """
237 Adds a new dataset to an existing XDMF file. Designed to be used with DREAM3D files.
239 Parameters:
240 xdmf_path (str): Path to the XDMF file.
241 dataset_name (str): The name of the new dataset.
242 data_array (np.ndarray): The data array to be referenced in the XDMF file
243 Returns:
244 None
245 """
246 # Read the existing XDMF file
247 with open(xdmf_path, "r") as file:
248 xdmf_content = file.readlines()
250 # Break the xdmf content into lines for easier manipulation
251 xdmf_content = [line.replace("\n", "") for line in xdmf_content]
253 # Make sure the shape of the data_array is compatible
254 if data_array.ndim == 3:
255 data_array = data_array.reshape(data_array.shape + (1,))
256 elif data_array.ndim < 3:
257 raise ValueError("data_array must be at least 3-dimensional")
258 elif data_array.ndim > 4: 258 ↛ 260line 258 didn't jump to line 260 because the condition on line 258 was always true
259 raise ValueError("data_array must be at most 4-dimensional")
260 dimensions = (
261 xdmf_content[["<Topology" in line for line in xdmf_content].index(True)]
262 .split("Dimensions=")[1]
263 .split('"')[1]
264 .strip()
265 )
266 data_array_dims = " ".join(map(str, np.array(data_array.shape[0:3]) + 1))
267 if dimensions != data_array_dims:
268 raise ValueError("data_array dimensions are not compatible with XDMF Topology dimensions")
270 # Make sure an entry with the same name does not already exist, if it does then just return
271 for line in xdmf_content:
272 if f'Attribute Name="{dataset_name}"' in line:
273 print(f"Dataset '{dataset_name}' already exists in XDMF file. Skipping addition.")
274 return
276 # Determine the insertion point (put the new entry at the end of the Grid section)
277 insertion_index = ["</Grid>" in line for line in xdmf_content].index(True)
279 # Gather relevant data for the new dataset
280 data_type, precision = _XDMF_DTYPE_FORMATS[data_array.dtype.type]
281 dimensions = " ".join(map(str, data_array.shape))
282 attribute_type = (
283 "Scalar"
284 if (data_array.ndim == 3) or ((data_array.ndim == 4) and (data_array.shape[-1] == 1))
285 else "Vector"
286 )
287 file_path = (
288 xdmf_content[[".dream3d:/" in line for line in xdmf_content].index(True)].strip().split("/")
289 )
290 file_path[-1] = dataset_name
291 file_path = "/".join(file_path)
293 # Create the new DataItem entry
294 xdmf_content.insert(
295 insertion_index,
296 f' <Attribute Name="{dataset_name}" AttributeType="{attribute_type}" Center="Cell">',
297 )
298 xdmf_content.insert(
299 insertion_index + 1,
300 f' <DataItem Format="HDF" Dimensions="{dimensions}" NumberType="{data_type}" Precision="{precision}" >',
301 )
302 xdmf_content.insert(
303 insertion_index + 2,
304 f" {file_path}",
305 )
306 xdmf_content.insert(insertion_index + 3, " </DataItem>")
307 xdmf_content.insert(insertion_index + 4, " </Attribute>")
309 # Write the modified content back to the XDMF file
310 with open(xdmf_path, "w") as file:
311 for line in xdmf_content:
312 file.write(line + "\n")
314 print(f"Added dataset '{dataset_name}' to XDMF file.")
315 return
318def save_to_dream3d(path: str | Path, ids_name: str, gnd_data: dict, fdm_data: np.ndarray) -> bool:
319 """
320 Saves GND and FDM data to a DREAM3D file.
322 Parameters:
323 path (str): Path to the DREAM3D file.
324 ids_name (str): Name of the DataArray containing the Feature IDs.
325 gnd_data (dict): Dictionary containing GND data arrays keyed by minimization method.
326 fdm_data (np.ndarray): FDM data array.
327 Returns:
328 bool: True if successful.
329 """
330 # Get the path to the ids array
331 ids_path = extract_path_from_h5(path, ids_name)
332 if ids_path is None:
333 raise KeyError(
334 f"Could not find Feature IDs data array with name '{ids_name}' in DREAM3D file. This is required to determine the cell data group."
335 )
337 # Use the ids path to find the cell data group and create paths for new data
338 cell_data_path = "/".join(ids_path.split("/")[:-1])
339 xdmf_path = str(Path(path).with_suffix(".xdmf"))
340 modify_xdmf = os.path.exists(xdmf_path)
341 if not modify_xdmf: 341 ↛ 347line 341 didn't jump to line 347 because the condition on line 341 was always true
342 print(
343 "WARNING: Could not find associated XDMF file. New datasets will only be added to the DREAM3D file. Please run DREAM3D to generate an updated XDMF file."
344 )
346 # Prep the data
347 data_shape = fdm_data.shape[1:] + (1,)
348 for m in gnd_data:
349 gnd_data[m] = gnd_data[m].sum(axis=0).reshape(data_shape).copy()
350 fdm_avg = fdm_data.mean(axis=0).reshape(data_shape)
351 fdm_max = fdm_data.max(axis=0).reshape(data_shape)
353 # Now open the new file and write the new data
354 with h5py.File(path, "r+") as h5:
355 cell_data = h5[cell_data_path]
357 # Add GND and FDM datasets
358 for m in gnd_data:
359 add_dataset_to_h5(cell_data, f"GND_{m}", gnd_data[m])
360 add_dataset_to_h5(cell_data, "FDM_avg", fdm_avg)
361 add_dataset_to_h5(cell_data, "FDM_max", fdm_max)
363 # Update the XDMF file if it exists
364 if modify_xdmf: 364 ↛ 365line 364 didn't jump to line 365 because the condition on line 364 was never true
365 for m in gnd_data:
366 add_dataset_to_xdmf(xdmf_path, f"GND_{m}", gnd_data[m])
367 add_dataset_to_xdmf(xdmf_path, "FDM_avg", fdm_avg)
368 add_dataset_to_xdmf(xdmf_path, "FDM_max", fdm_max)
370 return True
373def extract_path_from_h5(h5_file_path: str, target_name: str) -> str:
374 """
375 Extracts the path to a dataset from an HDF5 file given its name.
377 Parameters:
378 h5_file_path (str): Path to the HDF5 file.
379 target_name (str): Name of the target dataset to extract.
380 Returns:
381 str: The path to the dataset, or None if not found.
382 """
383 with h5py.File(h5_file_path, "r") as h5:
385 def recursive_search(name, obj):
386 if isinstance(obj, h5py.Dataset):
387 if name == target_name or name.endswith("/" + target_name):
388 return name
389 else:
390 return None
391 for key, item in obj.items():
392 result = recursive_search(f"{name}/{key}", item)
393 if result:
394 return result
395 return None
397 path = recursive_search("", h5)
398 return path
401def extract_attribute_from_h5(h5_file_path: str, attribute_name: str) -> any:
402 """
403 Extracts an attribute from an HDF5 file given its name.
405 Parameters:
406 h5_file_path (str): Path to the HDF5 file.
407 attribute_name (str): Name of the target attribute to extract.
408 Returns:
409 The value of the attribute, or None if not found.
410 """
411 if not os.path.exists(h5_file_path):
412 raise FileNotFoundError(f"HDF5 file not found at path: {h5_file_path}")
414 with h5py.File(h5_file_path, "r") as h5:
416 def recursive_search(name, obj):
417 print("Visiting:", name)
418 if isinstance(obj, h5py.Group):
419 print(" - attrs:", obj.attrs.keys())
420 if attribute_name in obj.attrs:
421 return obj.attrs[attribute_name]
422 else:
423 return None
424 print(" - Did not find attribute in group")
425 for key, item in obj.items():
426 result = recursive_search(f"{name}/{key}", item)
427 print(" - Result from", key, ":", result)
428 if result is not None:
429 return result
430 return None
432 attribute_value = recursive_search("", h5)
433 return attribute_value
436def extract_data_from_h5(h5_file_path: str, target_name: str) -> np.ndarray:
437 """
438 Extracts a data array from an HDF5 file given its name.
440 Parameters:
441 h5_file_path (str): Path to the HDF5 file.
442 target_name (str): Name of the target data array to extract.
443 Returns:
444 np.ndarray: The extracted data array, or None if not found.
445 """
446 if not os.path.exists(h5_file_path):
447 raise FileNotFoundError(f"HDF5 file not found at path: {h5_file_path}")
449 with h5py.File(h5_file_path, "r") as h5:
451 def recursive_search(name, obj):
452 if isinstance(obj, h5py.Dataset):
453 if name == target_name or name.endswith("/" + target_name):
454 return obj[...]
455 else:
456 return None
457 for key, item in obj.items():
458 result = recursive_search(f"{name}/{key}", item)
459 if result is not None:
460 return result
461 return None
463 data_array = recursive_search("", h5)
464 return data_array
467def save_npz(gnd_data: dict[np.ndarray], fdm_data: np.ndarray, folder: str | Path) -> None:
468 """
469 Saves GND and FDM data arrays to a .npz file.
471 Parameters:
472 gnd_data (dict[np.ndarray]): Dictionary containing GND data arrays keyed by minimization method.
473 fdm_data (np.ndarray): FDM data array.
474 folder (str): Folder to save the .npz file.
475 Returns:
476 None
477 """
478 for m in gnd_data:
479 np.save(Path(folder) / f"gnd_{m}.npy", gnd_data[m])
480 np.save(Path(folder) / "fdm.npy", fdm_data)
483def remove_npz(folder: str | Path) -> None:
484 """
485 Removes GND and FDM .npy files from the specified folder.
487 Parameters:
488 gnd_data (dict[np.ndarray]): Dictionary containing GND data arrays keyed by minimization method.
489 fdm_data (np.ndarray): FDM data array.
490 folder (str): Folder to remove the .npy files from.
491 Returns:
492 None
493 """
494 if (Path(folder) / "fdm.npy").exists():
495 Path(folder).joinpath("fdm.npy").unlink()
496 if (Path(folder) / "gnd_l1.npy").exists():
497 Path(folder).joinpath("gnd_l1.npy").unlink()
498 if (Path(folder) / "gnd_l2.npy").exists():
499 Path(folder).joinpath("gnd_l2.npy").unlink()
502def generate_images(gnd_data: dict[np.ndarray], fdm_data: np.ndarray, folder: str | Path) -> None:
503 """
504 Generates and saves images of GND and FDM data arrays as .npy files.
506 Parameters:
507 gnd_data (dict[np.ndarray]): Dictionary containing GND data arrays keyed by minimization method.
508 fdm_data (np.ndarray): FDM data array.
509 folder (str): Folder to save the .npy files.
510 Returns:
511 None
512 """
513 folder = Path(folder) / "images"
514 folder.mkdir(parents=True, exist_ok=True)
516 for m in gnd_data:
517 d = np.log10(np.clip(gnd_data[m], a_min=1e-10, a_max=None))
518 mx = d.max()
519 for i in range(d.shape[0]):
520 img_array = (d[i] / mx * 255).astype(np.uint8)
521 img = Image.fromarray(img_array)
522 img.save(folder / f"gnd_{m}_{i}.png")
524 d_sum = d.sum(axis=0)
525 img_array = (d_sum / d_sum.max() * 255).astype(np.uint8)
526 img = Image.fromarray(img_array)
527 img.save(folder / f"gnd_{m}_sum.png")
529 fdm_avg = fdm_data.mean(axis=0)
530 img_array = (fdm_avg / fdm_avg.max() * 255).astype(np.uint8)
531 img = Image.fromarray(img_array)
532 img.save(folder / "fdm_avg.png")
534 fdm_max = fdm_data.max(axis=0)
535 img_array = (fdm_max / fdm_max.max() * 255).astype(np.uint8)
536 img = Image.fromarray(img_array)
537 img.save(folder / "fdm_max.png")
538 return