diff --git a/docs/method.rst b/docs/method.rst index cb0c2932..a50f8988 100644 --- a/docs/method.rst +++ b/docs/method.rst @@ -127,75 +127,37 @@ The assignment algorithm begins by randomly assigning orientations obtained from EBSD reading ------------ -Kanapy can read two dimensional EBSD maps and convert the measured pixel -orientations into a graph representation of the microstructure. The graph is -intended to describe the grain structure on the EBSD map after region -segmentation, graph cleanup, and node merging. It can also be plotted directly -on the EBSD inverse pole figure map for visual inspection. - -The EBSD workflow is handled by :class:`kanapy.texture.EBSDmap`. During map -reading, Kanapy stores the EBSD pixel orientations, separates the retained -phases, and builds graph data for the selected phase. For each phase, the graph -data contains the final graph, diagnostic information from the initial graph, -and merge information from the cleanup procedure. +Kanapy can read two dimensional EBSD maps and convert the measured pixel orientations into a graph representation of the microstructure. The graph is intended to describe the grain structure on the EBSD map after region segmentation, graph cleanup, and node merging. It can also be plotted directly on the EBSD inverse pole figure map for visual inspection. + +The EBSD workflow is handled by :class:`kanapy.texture.EBSDmap`. During map reading, Kanapy stores the EBSD pixel orientations, separates the retained phases, and builds graph data for the selected phase. For each phase, the graph data contains the final graph, diagnostic information from the initial graph, and merge information from the cleanup procedure. Initial region construction ^^^^^^^^^^^^^^^^^^^^^^^^^^^ -Initial EBSD regions are identified by local pixel to pixel misorientation. The -procedure is implemented in :func:`kanapy.texture.find_similar_regions_by_misorientation`. -Starting from each unassigned pixel, a breadth first search grows a region by -adding neighboring pixels whose local misorientation is below the selected -threshold. The misorientation calculation uses the crystal symmetry of the -current phase. +Initial EBSD regions are identified by local pixel to pixel misorientation. The procedure is implemented in :func:`kanapy.texture.find_similar_regions_by_misorientation`. Starting from each unassigned pixel, a breadth first search grows a region by adding neighboring pixels whose local misorientation is below the selected threshold. The misorientation calculation uses the crystal symmetry of the current phase. -For the two dimensional EBSD graph path, neighboring pixels are evaluated with -4-neighbor connectivity. This gives each pixel direct horizontal and vertical -neighbors on the EBSD map. The resulting labeled regions form the initial graph -nodes. +For the two dimensional EBSD graph path, neighboring pixels are evaluated with 4-neighbor connectivity. This gives each pixel direct horizontal and vertical neighbors on the EBSD map. The resulting labeled regions form the initial graph nodes. Graph construction and node data ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ -The labeled EBSD regions are converted into graph nodes by -:func:`kanapy.texture.build_graph_from_labeled_pixels`. Each graph node -represents one connected labeled region. Node attributes include the pixel list, -the node center, the number of pixels, the original pixel rotations, and the -mean node orientation. +The labeled EBSD regions are converted into graph nodes by :func:`kanapy.texture.build_graph_from_labeled_pixels`. Each graph node represents one connected labeled region. Node attributes include the pixel list, the node center, the number of pixels, the original pixel rotations, and the mean node orientation. -The mean orientation of a node is computed by -:func:`kanapy.texture.mean_orientation_data`. Pixel quaternions are normalized, -symmetry equivalent orientations are aligned to a reference orientation, and the -``q`` versus ``-q`` quaternion ambiguity is handled before averaging. The final -mean quaternion is obtained from the dominant eigenvector of the quaternion -accumulator matrix. +The mean orientation of a node is computed by :func:`kanapy.texture.mean_orientation_data`. Pixel quaternions are normalized, symmetry equivalent orientations are aligned to a reference orientation, and the ``q`` versus ``-q`` quaternion ambiguity is handled before averaging. The final mean quaternion is obtained from the dominant eigenvector of the quaternion accumulator matrix. -Graph edges are added between neighboring labeled regions. Label ``0`` is -treated as background and is not used as a graph node. +Graph edges are added between neighboring labeled regions. Label ``0`` is treated as background and is not used as a graph node. Boundary artifact cleanup ^^^^^^^^^^^^^^^^^^^^^^^^^ -Some small regions can appear mainly along grain boundaries. Kanapy computes -node boundary information with :func:`kanapy.texture.get_node_boundary_stats`. -For each node, the function checks local neighboring pixels and reports -boundary pixel fraction, neighboring labels, map edge contact, bounding box -size, and bounding box fill fraction. +Some small regions can appear mainly along grain boundaries. Kanapy computes node boundary information with :func:`kanapy.texture.get_node_boundary_stats`. For each node, the function checks local neighboring pixels and reports boundary pixel fraction, neighboring labels, map edge contact, bounding box size, and bounding box fill fraction. -Small boundary artifact nodes can be merged with -:func:`kanapy.texture.merge_boundary_artifact_nodes`. The merge target is chosen -from neighboring nodes by comparing the node orientations. After node merging, -the mean orientation of the merged node is recomputed from the original pixel -orientations. +Small boundary artifact nodes can be merged with :func:`kanapy.texture.merge_boundary_artifact_nodes`. The merge target is chosen from neighboring nodes by comparing the node orientations. After node merging, the mean orientation of the merged node is recomputed from the original pixel orientations. Graph plotting ^^^^^^^^^^^^^^ -The final EBSD graph can be plotted on the EBSD inverse pole figure map by -setting ``show_graph=True`` when creating an :class:`kanapy.texture.EBSDmap` -object. The graph overlay is produced by -:meth:`kanapy.texture.EBSDmap.plot_graph_overlay`. Node centers are shown as -black points and graph edges are shown as black lines. +The final EBSD graph can be plotted on the EBSD inverse pole figure map by setting ``show_graph=True`` when creating an :class:`kanapy.texture.EBSDmap` object. The graph overlay is produced by :meth:`kanapy.texture.EBSDmap.plot_graph_overlay`. Node centers are shown as black points and graph edges are shown as black lines. Example usage: @@ -205,30 +167,28 @@ Example usage: ebsd = knpy.EBSDmap("p558_250x_1.ang", show_graph=True) -The same plotting method can also save the graph image to file when an output -path is provided. - -Example scripts -^^^^^^^^^^^^^^^ +The same plotting method can also save the graph image to file when an output path is provided. -Two example scripts are provided in ``examples/EBSD_graph_analysis``. +Example script +^^^^^^^^^^^^^^ -``plot_ebsd_graph_manually.py`` - Runs EBSD graph construction, saves the final graph overlay image, and - writes a short text summary of the graph. The output is written to - ``plot_ebsd_graph_manually_result``. +The graph workflow example is provided in ``examples/EBSD_graph_analysis``. ``run_ebsd_analysis_with_graph.py`` - Runs the standard EBSD analysis workflow and adds the final graph overlay - by using ``show_graph=True``. + Builds an in memory :class:`kanapy.graph_workflow.EBSDGraphResult` object and writes only the explicitly enabled graph handoff files. This example uses the graph only path and does not run grain statistics extraction or the legacy interactive plotting workflow. The output is written to ``2D_graph_result``. + +The graph workflow can be called directly from Python by using :func:`kanapy.graph_workflow.build_ebsd_graph`. The main data transfer object is :class:`kanapy.graph_workflow.EBSDGraphResult`; PKL, JSON, and PNG files are written only when enabled through :class:`kanapy.graph_workflow.EBSDGraphOutputOptions`. Users who need the classic EBSD statistics and interactive Kanapy plots can still create :class:`kanapy.texture.EBSDmap` directly with the usual plotting flags. + +The example keeps graph step figures and local zoom diagnostics disabled by default. These outputs are mainly useful for debugging or preparing publication figures, while routine workflows usually need only the final graph handoff files. -The manual plotting example produces: +The graph workflow example produces: .. code-block:: text - plot_ebsd_graph_manually_result/ebsd_graph_overlay.png - plot_ebsd_graph_manually_result/ebsd_graph_summary.txt + 2D_graph_result/ebsd_graph.pkl + 2D_graph_result/ebsd_graph_summary.json + 2D_graph_result/graph_overlay.png + 2D_graph_result/graph_ipf_map.png -The overlay image shows the final graph node centers and graph edges on the -EBSD inverse pole figure map. +The overlay image shows the final graph node centers and graph edges on the EBSD inverse pole figure map. diff --git a/examples/EBSD_graph_analysis/plot_ebsd_graph_manually.py b/examples/EBSD_graph_analysis/plot_ebsd_graph_manually.py deleted file mode 100644 index e065307e..00000000 --- a/examples/EBSD_graph_analysis/plot_ebsd_graph_manually.py +++ /dev/null @@ -1,143 +0,0 @@ -""" -Read an EBSD map and visualize the resulting 2D microstructure graph. - -The graph nodes represent detected grains or grain regions. Graph edges mark -neighboring nodes in the EBSD pixel map. -""" - -from pathlib import Path -import sys - -import numpy as np - -repo_src = Path(__file__).resolve().parents[2] / "src" -if repo_src.exists(): - sys.path.insert(0, str(repo_src)) - -import kanapy as knpy - - -def get_graph_summary_lines(ms_data): - """Create a concise summary of the EBSD graph. - - Parameters - ---------- - ms_data : dict - Phase data dictionary from ``ebsd.ms_data``. - - Returns - ------- - list[str] - Summary lines for terminal and text-file output. - """ - graph = ms_data["graph"] - degrees = np.array([degree for _, degree in graph.degree], dtype=float) - npix = np.array([node["npix"] for _, node in graph.nodes.items()], dtype=float) - - lines = [ - "*** EBSD graph information ***", - f"Phase: #{ms_data['index']} ({ms_data['name']})", - f"Nodes: {graph.number_of_nodes()}", - f"Edges: {graph.number_of_edges()}", - f"Mean node degree: {degrees.mean():.2f}", - f"Median node area: {np.median(npix):.1f} pixels", - ] - - if "graph_initial" in ms_data: - initial_graph = ms_data["graph_initial"] - lines.append( - f"Initial nodes before graph cleanup: {initial_graph.number_of_nodes()}" - ) - - if "merge_debug" in ms_data: - performed_merges = [] - for row in ms_data["merge_debug"]: - if row.get("will_merge", False): - performed_merges.append(row) - lines.append(f"Recorded graph merges: {len(performed_merges)}") - - return lines - - -def select_major_phase(ebsd): - """Return the retained phase with the largest volume fraction. - - Parameters - ---------- - ebsd : kanapy.texture.EBSDmap - EBSD map object. - - Returns - ------- - tuple[int, dict] - Phase list index and phase data dictionary. - """ - if len(ebsd.ms_data) == 0: - raise ValueError("No EBSD phase passed the volume fraction filter.") - - phase_rows = list(enumerate(ebsd.ms_data)) - return max(phase_rows, key=lambda row: row[1]["vf"]) - - -def write_graph_summary(ms_data, out_file): - """Print and save graph summary information. - - Parameters - ---------- - ms_data : dict - Phase data dictionary from ``ebsd.ms_data``. - out_file : pathlib.Path - Text file where the summary is saved. - """ - lines = get_graph_summary_lines(ms_data) - text = "\n".join(lines) - print(f"\n{text}") - out_file.write_text(text + "\n", encoding="utf-8") - - -def main(): - """Run the EBSD graph visualization example.""" - base_dir = Path(__file__).resolve().parent - input_file = base_dir / "p558_250x_1.ang" - result_dir = base_dir / "plot_ebsd_graph_manually_result" - result_dir.mkdir(exist_ok=True) - - overlay_file = result_dir / "ebsd_graph_overlay.png" - summary_file = result_dir / "ebsd_graph_summary.txt" - - vf_min = 0.03 # minimum phase volume fraction to keep - max_angle = 5.0 # maximum misorientation angle within one grain in degrees - min_size = 20.0 # minimum grain size in pixels - connectivity = 4 # use edge sharing pixel neighbors - - if not input_file.exists(): - raise FileNotFoundError( - f"EBSD input file not found: {input_file}. " - "Set input_file to the path of an ANG or CTF file before running this example." - ) - - ebsd = knpy.EBSDmap( - str(input_file), - vf_min=vf_min, - gs_min=min_size, - max_angle=max_angle, - connectivity=connectivity, - show_plot=False, - show_grains=False, - show_hist=False, - ) - - iphase, ms_data = select_major_phase(ebsd) - ebsd.plot_graph_overlay( - iphase=iphase, - save_path=overlay_file, - show=False, - ) - write_graph_summary(ms_data, summary_file) - - print(f"Graph overlay saved to: {overlay_file}") - print(f"Graph summary saved to: {summary_file}") - - -if __name__ == "__main__": - main() diff --git a/examples/EBSD_graph_analysis/plot_ebsd_graph_manually_result/ebsd_graph_overlay.png b/examples/EBSD_graph_analysis/plot_ebsd_graph_manually_result/ebsd_graph_overlay.png deleted file mode 100644 index a9528444..00000000 Binary files a/examples/EBSD_graph_analysis/plot_ebsd_graph_manually_result/ebsd_graph_overlay.png and /dev/null differ diff --git a/examples/EBSD_graph_analysis/plot_ebsd_graph_manually_result/ebsd_graph_summary.txt b/examples/EBSD_graph_analysis/plot_ebsd_graph_manually_result/ebsd_graph_summary.txt deleted file mode 100644 index d759f161..00000000 --- a/examples/EBSD_graph_analysis/plot_ebsd_graph_manually_result/ebsd_graph_summary.txt +++ /dev/null @@ -1,8 +0,0 @@ -*** EBSD graph information *** -Phase: #1 (Fe) -Nodes: 185 -Edges: 503 -Mean node degree: 5.44 -Median node area: 1645.0 pixels -Initial nodes before graph cleanup: 323 -Recorded graph merges: 137 diff --git a/examples/EBSD_graph_analysis/run_ebsd_analysis_with_graph.py b/examples/EBSD_graph_analysis/run_ebsd_analysis_with_graph.py index 4f6595aa..1d989a25 100644 --- a/examples/EBSD_graph_analysis/run_ebsd_analysis_with_graph.py +++ b/examples/EBSD_graph_analysis/run_ebsd_analysis_with_graph.py @@ -1,89 +1,95 @@ """ -Run the original EBSD analysis workflow with graph visualization enabled. +Build a reusable EBSD graph handoff from an EBSD map. -This script keeps the EBSD statistics and optional RVE-generation workflow from -the original example, but uses the local ANG file and the current 2D graph -construction settings. In addition to the standard EBSD maps, grain statistics, -and histogram plots, it displays the final microstructure graph as black node -centers and graph edges on the IPF map. +This example returns an ``EBSDGraphResult`` instance and writes only the graph +handoff files explicitly enabled below. It uses Kanapy's graph only path and +skips the grain statistics and legacy interactive plotting workflow. Users who +need the classic interactive Kanapy analysis plots should create an +``EBSDmap`` object directly with the usual plotting flags. Routine graph users +usually need only the four handoff files written by this example. """ from pathlib import Path import sys -import numpy as np -import matplotlib.pyplot as plt -from scipy.stats import lognorm - repo_src = Path(__file__).resolve().parents[2] / "src" if repo_src.exists(): sys.path.insert(0, str(repo_src)) import kanapy as knpy -base_dir = Path(__file__).resolve().parent -fname = base_dir / "p558_250x_1.ang" # cross_section - -nvox = 30 # number of voxels per side -box_length = 50 # side length of generated RVE in micron -periodic = False # create RVE with periodic structure - -vf_min = 0.03 # minimum volume fraction of phases to be considered -max_angle = 5.0 # maximum misorientation angle within one grain in degrees -min_size = 20.0 # minimum grain size in pixels -connectivity = 4 # use edge sharing neighbors - -# read EBSD map and evaluate statistics of microstructural features -ebsd = knpy.EBSDmap(str(fname), vf_min=vf_min, gs_min=min_size, max_angle=max_angle, show_plot=True, - connectivity=connectivity, show_grains=True, show_hist=True, show_graph=True) -# ebsd.showIPF() -ms_data = ebsd.ms_data[0] # analyse only data for majority phase with order parameter "0" -gs_param = ms_data['gs_param'] # lognorm distr of grain size: [std dev., location, scale] -ar_param = ms_data['ar_param'] # lognorm distr. of aspect ratios [std dev., loc., scale] -om_param = ms_data['om_param'] # normal distribution of tilt angles [std dev., mean] -matname = ms_data['name'] # material name -print('*** Statistical information on microstructure ***') -print(f'=== Phase: {matname} ===') -print('==== Grain size (equivalent grain diameter) ====') -print(f'scale: {gs_param[2].round(3)} micron, ' - f'std. deviation: {gs_param[0].round(3)}') -print('==== Aspect ratio ====') -print(f'scale: {ar_param[2].round(3)}, ' - f'std. deviation: {ar_param[0].round(3)}') -print('==== Tilt angle ====') -print(f'most frequent value: {0.5*(om_param[1] + np.pi):.3f} rad, ' + - f'std. deviation: {om_param[0]:.3f}') - -""" Analyze microstructure graph """ -G = ms_data['graph'] - -# plot components of major axis -va_ = [] -for nn in G.nodes: - va_.append(G.nodes[nn]['maj_ax']) -va_ = np.array(va_) -plt.hist(va_, bins=40, density=True, label=['x', 'y']) -plt.title('Components of major axis') -plt.legend() -plt.show() - -# plot angles from arctan2 -om1 = np.arctan2(np.abs(va_[:, 1]), va_[:, 0]) -om2 = np.arccos(va_[:, 0]) -plt.hist(np.c_[om1, om2], bins=40, density=True, label=['arctan2', 'arccos']) -plt.title('Tilt angle by arctan2 vs. acos') -plt.legend() -plt.show() - -# compare area of hull and pixels -dar = ms_data['delta_x'] * ms_data['delta_y'] -area = [] -for nn in G.nodes: - area.append(G.nodes[nn]['npix'] * dar) -eqd = 2.0 * (np.array(area) / np.pi) ** 0.5 -plt.hist(np.c_[ms_data['gs_data'], eqd], bins=15, density=True, label=['hull', 'pix']) -plt.title('Grain diameters obtained from pixels and from convex hull') -plt.xlabel('Grain diameter') -plt.ylabel('Frequency') -plt.legend() -plt.show() + +def print_graph_handoff_summary(result): + """ + Print a concise summary of the graph handoff. + + The terminal output confirms graph construction and lists only the files + explicitly requested by this example. + + Parameters + ---------- + result : kanapy.graph_workflow.EBSDGraphResult + In memory graph construction result. + """ + graph = result.graph + print("----------------------------------------") + print("EBSD Graph Result") + print(f"Phase: #{result.phase_index} ({result.phase_name})") + print(f"Graph nodes: {graph.number_of_nodes()}") + print(f"Graph edges: {graph.number_of_edges()}") + + if result.written_files: + print("Written files:") + for output_name, output_path in result.written_files.items(): + print(f" {output_name}: {output_path}") + else: + print("Written files: none") + + if result.elapsed_time_seconds is not None: + print(f"Elapsed time: {result.elapsed_time_seconds:.2f} s") + + print("----------------------------------------") + + +def main(): + """ + Build the EBSD graph handoff for downstream analysis. + + The API returns an in memory ``EBSDGraphResult`` and skips grain statistics + extraction. PKL, JSON, overlay PNG, and IPF map PNG files are written here + only because the example explicitly enables them. + """ + base_dir = Path(__file__).resolve().parent + input_file = base_dir / "p558_250x_1.ang" + output_dir = base_dir / "2D_graph_result" + + graph_config = knpy.EBSDGraphConfig( + min_phase_pixel_fraction=0.03, + max_misorientation_degrees=5.0, + min_grain_size_pixels=20.0, + connectivity=4, + boundary_artifact_size_factor=2.0, + boundary_artifact_fraction_min=0.85, + convex_hull_min_points=3, + similar_orientation_merge_degrees=5.0, + ) + output_options = knpy.EBSDGraphOutputOptions( + output_dir=output_dir, + clear_output_dir=True, + write_pickle=True, + write_summary_json=True, + write_graph_overlay=True, + write_ipf_map=True, + write_step_figures=False, + ) + + result = knpy.build_ebsd_graph( + input_file=input_file, + config=graph_config, + output_options=output_options, + ) + print_graph_handoff_summary(result) + + +if __name__ == "__main__": + main() diff --git a/src/kanapy/__init__.py b/src/kanapy/__init__.py index aa051bb5..44eba6bb 100644 --- a/src/kanapy/__init__.py +++ b/src/kanapy/__init__.py @@ -28,6 +28,14 @@ from .texture import EBSDmap, ODF, createOriset, createOrisetRandom, \ get_ipf_colors, plot_pole_figure, plot_pole_figure_proj +from .graph_workflow import ( + EBSDGraphConfig, + EBSDGraphOutputOptions, + EBSDGraphResult, + build_ebsd_graph, + load_ebsd_graph, + write_graph_result_outputs, +) log_level = 20 # Levels for logging: 10: DEBUG, 20: INFO, 30: WARNING, 40: ERROR logging.basicConfig(level=log_level) # set log level @@ -41,7 +49,9 @@ "import_stats", "write_stats", "start", "EBSDmap", "ODF", "createOriset", "createOrisetRandom", "get_ipf_colors", "plot_pole_figure", "plot_pole_figure_proj", "plot_voxels_3D", "plot_polygons_3D", "plot_mean_ellipsoids_from_stats", - "triple_surf"] + "triple_surf", "EBSDGraphConfig", "EBSDGraphOutputOptions", + "EBSDGraphResult", "build_ebsd_graph", "load_ebsd_graph", + "write_graph_result_outputs"] if triple_surf: from .core import create_ref_ell diff --git a/src/kanapy/graph_workflow.py b/src/kanapy/graph_workflow.py new file mode 100644 index 00000000..1b0406b2 --- /dev/null +++ b/src/kanapy/graph_workflow.py @@ -0,0 +1,688 @@ +""" +In memory EBSD graph workflow objects. + +This module keeps graph construction results as Python objects and makes file +outputs explicit choices by the caller. +""" + +from __future__ import annotations + +import json +import pickle +import shutil +import time +from dataclasses import dataclass, field, replace +from pathlib import Path +from typing import Any, Optional, Union + +import matplotlib.pyplot as plt +import numpy as np + +from .texture import ( + DEFAULT_BOUNDARY_ARTIFACT_FRACTION_MIN, + DEFAULT_BOUNDARY_ARTIFACT_SIZE_FACTOR, + DEFAULT_CONNECTIVITY, + DEFAULT_MAX_MISORIENTATION_DEGREES, + DEFAULT_MIN_VOLUME_FRACTION, + DEFAULT_SIMILAR_ORIENTATION_MERGE_DEGREES, + EBSDmap, + CONVEX_HULL_MIN_POINTS, + MIN_GRAIN_SIZE_PIXELS, +) + + +EBSD_GRAPH_FILENAME = "ebsd_graph.pkl" +GRAPH_SUMMARY_FILENAME = "ebsd_graph_summary.json" +GRAPH_OVERLAY_FILENAME = "graph_overlay.png" +GRAPH_IPF_MAP_FILENAME = "graph_ipf_map.png" + + +@dataclass +class EBSDGraphConfig: + """ + Store user configurable thresholds for EBSD graph construction. + + Defaults match the existing graph workflow while keeping each threshold + visible to callers. + """ + + min_phase_pixel_fraction: float = DEFAULT_MIN_VOLUME_FRACTION + max_misorientation_degrees: float = DEFAULT_MAX_MISORIENTATION_DEGREES + min_grain_size_pixels: float = MIN_GRAIN_SIZE_PIXELS + connectivity: int = DEFAULT_CONNECTIVITY + boundary_artifact_size_factor: float = DEFAULT_BOUNDARY_ARTIFACT_SIZE_FACTOR + boundary_artifact_fraction_min: float = DEFAULT_BOUNDARY_ARTIFACT_FRACTION_MIN + convex_hull_min_points: int = CONVEX_HULL_MIN_POINTS + similar_orientation_merge_degrees: float = ( + DEFAULT_SIMILAR_ORIENTATION_MERGE_DEGREES + ) + + def __post_init__(self) -> None: + """ + Validate graph threshold settings. + + Convex hull construction needs at least three points in two + dimensions. + + Raises + ------ + ValueError + If ``convex_hull_min_points`` is lower than three. + """ + if self.convex_hull_min_points < CONVEX_HULL_MIN_POINTS: + raise ValueError( + f"convex_hull_min_points must be at least {CONVEX_HULL_MIN_POINTS}" + ) + + def to_dict(self) -> dict[str, Union[float, int]]: + """ + Return graph construction thresholds as plain values. + + The dictionary is used in summaries and serialized handoff metadata. + + Returns + ------- + dict + Graph threshold and connectivity settings. + """ + return { + "min_phase_pixel_fraction": float(self.min_phase_pixel_fraction), + "max_misorientation_degrees": float(self.max_misorientation_degrees), + "min_grain_size_pixels": float(self.min_grain_size_pixels), + "connectivity": int(self.connectivity), + "boundary_artifact_size_factor": float( + self.boundary_artifact_size_factor + ), + "boundary_artifact_fraction_min": float( + self.boundary_artifact_fraction_min + ), + "convex_hull_min_points": int(self.convex_hull_min_points), + "similar_orientation_merge_degrees": float( + self.similar_orientation_merge_degrees + ), + } + + +@dataclass +class EBSDGraphOutputOptions: + """ + Store optional graph output settings. + + No files are written unless their corresponding option is enabled. The + example script opts into the four paper handoff files. + """ + + output_dir: Optional[Union[str, Path]] = None + write_pickle: bool = False + write_summary_json: bool = False + write_graph_overlay: bool = False + write_ipf_map: bool = False + write_step_figures: bool = False + clear_output_dir: bool = False + pickle_filename: str = EBSD_GRAPH_FILENAME + summary_filename: str = GRAPH_SUMMARY_FILENAME + graph_overlay_filename: str = GRAPH_OVERLAY_FILENAME + ipf_map_filename: str = GRAPH_IPF_MAP_FILENAME + + @property + def any_file_output_enabled(self) -> bool: + """ + Return whether any graph output file is enabled. + + Returns + ------- + bool + True when a configured output writes to disk. + """ + return ( + self.write_pickle + or self.write_summary_json + or self.write_graph_overlay + or self.write_ipf_map + or self.write_step_figures + ) + + @property + def requires_output_dir(self) -> bool: + """ + Return whether the configured options need an output directory. + + Returns + ------- + bool + True when an output directory is required. + """ + return self.any_file_output_enabled or self.clear_output_dir + + def resolved_output_dir(self) -> Path: + """ + Return the configured output directory as a path. + + Raises + ------ + ValueError + If an output directory is required but was not configured. + + Returns + ------- + pathlib.Path + Resolved output directory path. + """ + if self.output_dir is None: + raise ValueError("output_dir is required when graph file output is enabled") + + return Path(self.output_dir) + + +@dataclass +class EBSDGraphResult: + """ + Store an EBSD graph handoff in memory. + + Downstream Python functions can consume this object directly instead of + reading a PKL, JSON, or PNG file. + """ + + graph: Any + phase_data: dict[str, Any] + phase_index: int + phase_name: str + rgb_image: np.ndarray + npx: int + sh_x: int + sh_y: int + source_input_file: Optional[Path] = None + config: EBSDGraphConfig = field(default_factory=EBSDGraphConfig) + elapsed_time_seconds: Optional[float] = None + ebsd: Optional[EBSDmap] = None + written_files: dict[str, Path] = field(default_factory=dict) + + @classmethod + def from_phase_data( + cls, + phase_data: dict[str, Any], + npx: int, + sh_x: int, + sh_y: int, + source_input_file: Optional[Union[str, Path]], + config: EBSDGraphConfig, + elapsed_time_seconds: Optional[float] = None, + ebsd: Optional[EBSDmap] = None, + ) -> "EBSDGraphResult": + """ + Build a graph result from one ``EBSDmap.ms_data`` phase row. + + Parameters + ---------- + phase_data : dict + Phase row containing the generated graph and RGB IPF image. + npx : int + Total number of pixels in the EBSD map. + sh_x : int + Map height in pixels. + sh_y : int + Map width in pixels. + source_input_file : str or pathlib.Path or None + Input EBSD file used to build the graph. + config : EBSDGraphConfig + Graph construction thresholds. + elapsed_time_seconds : float, optional + Runtime for graph construction. + ebsd : EBSDmap, optional + Originating EBSD map object. + + Returns + ------- + EBSDGraphResult + In memory graph handoff object. + """ + rgb_image = _phase_rgb_image(phase_data, npx, sh_x, sh_y) + source_path = None + if source_input_file is not None: + source_path = Path(source_input_file) + + return cls( + graph=phase_data["graph"], + phase_data=phase_data, + phase_index=int(phase_data["index"]), + phase_name=str(phase_data["name"]), + rgb_image=rgb_image, + npx=int(npx), + sh_x=int(sh_x), + sh_y=int(sh_y), + source_input_file=source_path, + config=config, + elapsed_time_seconds=elapsed_time_seconds, + ebsd=ebsd, + ) + + def to_payload(self) -> dict[str, Any]: + """ + Return the full Python payload for optional pickle output. + + Returns + ------- + dict + Serializable handoff payload containing the graph object. + """ + return { + "graph": self.graph, + "phase": { + "index": self.phase_index, + "name": self.phase_name, + }, + "rgb_image": self.rgb_image, + "npx": self.npx, + "sh_x": self.sh_x, + "sh_y": self.sh_y, + "source_input_file": ( + str(self.source_input_file) if self.source_input_file else None + ), + "thresholds": self.config.to_dict(), + "elapsed_time_seconds": self.elapsed_time_seconds, + "graph_node_count": self.graph.number_of_nodes(), + "graph_edge_count": self.graph.number_of_edges(), + } + + def to_summary(self) -> dict[str, Any]: + """ + Return a compact JSON friendly graph summary. + + Returns + ------- + dict + JSON serializable graph metadata. + """ + summary = { + "phase_index": self.phase_index, + "phase_name": self.phase_name, + "graph_node_count": self.graph.number_of_nodes(), + "graph_edge_count": self.graph.number_of_edges(), + "connectivity": int(self.config.connectivity), + "source_input_file": ( + str(self.source_input_file) if self.source_input_file else None + ), + "thresholds": self.config.to_dict(), + "elapsed_time_seconds": self.elapsed_time_seconds, + } + + for key, path in self.written_files.items(): + summary[key] = Path(path).name + + return summary + + +def _phase_rgb_image( + phase_data: dict[str, Any], + npx: int, + sh_x: int, + sh_y: int, +) -> np.ndarray: + """ + Return the phase RGB image normalized to zero to one. + + Parameters + ---------- + phase_data : dict + Phase data containing an ``rgb_im`` field. + npx : int + Total number of map pixels. + sh_x : int + Map height in pixels. + sh_y : int + Map width in pixels. + + Returns + ------- + numpy.ndarray + RGB image with shape ``(sh_x, sh_y, 3)``. + """ + rgb_image = np.asarray(phase_data["rgb_im"], dtype=float) + if rgb_image.shape == (npx, 3): + rgb_image = rgb_image.reshape((sh_x, sh_y, 3)) + if rgb_image.max(initial=0.0) > 1.0: + rgb_image = rgb_image / 255.0 + + return rgb_image + + +def select_major_phase(ebsd: EBSDmap) -> tuple[int, dict[str, Any]]: + """ + Return the retained phase with the largest volume fraction. + + Parameters + ---------- + ebsd : EBSDmap + Analyzed EBSD map object. + + Raises + ------ + ValueError + If no phase passed the configured volume fraction filter. + + Returns + ------- + tuple + Phase list index and selected phase data dictionary. + """ + if len(ebsd.ms_data) == 0: + raise ValueError("No EBSD phase passed the volume fraction filter.") + + phase_rows = list(enumerate(ebsd.ms_data)) + return max(phase_rows, key=lambda row: row[1]["vf"]) + + +def write_graph_result_outputs( + result: EBSDGraphResult, + output_options: EBSDGraphOutputOptions, +) -> dict[str, Path]: + """ + Write selected graph handoff files. + + Parameters + ---------- + result : EBSDGraphResult + In memory graph result to serialize or plot. + output_options : EBSDGraphOutputOptions + File output settings. + + Returns + ------- + dict + Mapping from output kind to written path. + """ + if not output_options.any_file_output_enabled: + return {} + + output_dir = output_options.resolved_output_dir() + if output_options.clear_output_dir and output_dir.exists(): + shutil.rmtree(output_dir) + + output_dir.mkdir(parents=True, exist_ok=True) + written_files: dict[str, Path] = {} + + if output_options.write_pickle: + pickle_path = output_dir / output_options.pickle_filename + with pickle_path.open("wb") as file_handle: + pickle.dump( + result.to_payload(), + file_handle, + protocol=pickle.HIGHEST_PROTOCOL, + ) + written_files["pickle"] = pickle_path + result.written_files["pickle"] = pickle_path + + if output_options.write_graph_overlay: + overlay_path = output_dir / output_options.graph_overlay_filename + _write_graph_overlay(result, overlay_path) + written_files["graph_overlay"] = overlay_path + result.written_files["graph_overlay"] = overlay_path + + if output_options.write_ipf_map: + ipf_path = output_dir / output_options.ipf_map_filename + _write_ipf_map(result, ipf_path) + written_files["ipf_map"] = ipf_path + result.written_files["ipf_map"] = ipf_path + + if output_options.write_summary_json: + summary_path = output_dir / output_options.summary_filename + with summary_path.open("w", encoding="utf-8") as file_handle: + json.dump(result.to_summary(), file_handle, indent=2) + file_handle.write("\n") + written_files["summary_json"] = summary_path + result.written_files["summary_json"] = summary_path + + return written_files + + +def _node_center( + graph: Any, + node_id: int, + image_shape: tuple[int, int], +) -> Optional[tuple[float, float]]: + """ + Return a graph node center in image pixel coordinates. + + Parameters + ---------- + graph : networkx.Graph + Microstructure graph containing the node. + node_id : int + Graph node ID. + image_shape : tuple of int + EBSD map shape as ``(rows, columns)``. + + Returns + ------- + tuple or None + Center as ``(x, y)`` image coordinates, or None if no center exists. + """ + node = graph.nodes[node_id] + center = node.get("center") + if center is not None: + dx = graph.graph.get("dx", 1.0) + dy = graph.graph.get("dy", 1.0) + return float(center[1]) / dy, float(center[0]) / dx + + pixels = node.get("pixels") + if pixels is None or len(pixels) == 0: + return None + + rows, cols = np.unravel_index(np.asarray(pixels, dtype=int), image_shape) + return float(np.mean(cols)), float(np.mean(rows)) + + +def _figure_size_for_image(image_shape: tuple[int, ...]) -> tuple[float, float]: + """ + Return a figure size that follows the image aspect ratio. + + Parameters + ---------- + image_shape : tuple of int + Image shape as ``(rows, columns)`` or ``(rows, columns, channels)``. + + Returns + ------- + tuple + Figure size in inches. + """ + rows, cols = image_shape[:2] + return 8.0 * cols / rows, 8.0 + + +def _format_image_axes(ax: Any, image_shape: tuple[int, ...]) -> None: + """ + Format an EBSD image axis with compact pixel coordinates. + + Parameters + ---------- + ax : matplotlib.axes.Axes + Axis containing an EBSD image. + image_shape : tuple of int + Image shape as ``(rows, columns)`` or ``(rows, columns, channels)``. + """ + rows, cols = image_shape[:2] + ax.set_xlim(-0.5, cols - 0.5) + ax.set_ylim(rows - 0.5, -0.5) + + +def _save_tight_figure(fig: Any, output_path: Path) -> None: + """ + Save a figure with compact margins. + + Parameters + ---------- + fig : matplotlib.figure.Figure + Figure to save. + output_path : pathlib.Path + Output image path. + """ + output_path.parent.mkdir(parents=True, exist_ok=True) + fig.savefig(output_path, dpi=150, bbox_inches="tight", pad_inches=0.03) + + +def _write_ipf_map(result: EBSDGraphResult, output_path: Path) -> None: + """ + Write the EBSD IPF map PNG. + + Parameters + ---------- + result : EBSDGraphResult + Graph result containing the IPF RGB image. + output_path : pathlib.Path + PNG output path. + """ + fig, ax = plt.subplots(figsize=_figure_size_for_image(result.rgb_image.shape)) + ax.imshow(result.rgb_image, origin="upper") + ax.set_title("EBSD map", fontsize=16, fontfamily="DejaVu Sans", pad=14) + _format_image_axes(ax, result.rgb_image.shape) + fig.tight_layout(pad=0.3) + _save_tight_figure(fig, output_path) + plt.close(fig) + + +def _write_graph_overlay(result: EBSDGraphResult, output_path: Path) -> None: + """ + Write the final graph overlay PNG. + + Parameters + ---------- + result : EBSDGraphResult + Graph result containing the final graph and IPF RGB image. + output_path : pathlib.Path + PNG output path. + """ + graph = result.graph + image_shape = result.rgb_image.shape[:2] + fig, ax = plt.subplots(figsize=_figure_size_for_image(result.rgb_image.shape)) + ax.imshow(result.rgb_image, origin="upper") + + for node_a, node_b in graph.edges(): + center_a = _node_center(graph, node_a, image_shape) + center_b = _node_center(graph, node_b, image_shape) + if center_a is None or center_b is None: + continue + + xa, ya = center_a + xb, yb = center_b + ax.plot([xa, xb], [ya, yb], color="black", linewidth=0.7) + + for node_id in graph.nodes: + center = _node_center(graph, node_id, image_shape) + if center is None: + continue + + x, y = center + ax.plot(x, y, "ko", markersize=2.5) + + node_count = graph.number_of_nodes() + edge_count = graph.number_of_edges() + ax.set_title( + f"Final graph ({node_count} nodes, {edge_count} edges)", + fontsize=16, + fontfamily="DejaVu Sans", + pad=14, + ) + _format_image_axes(ax, result.rgb_image.shape) + fig.tight_layout(pad=0.3) + _save_tight_figure(fig, output_path) + plt.close(fig) + + +def load_ebsd_graph(path: Union[str, Path]) -> dict[str, Any]: + """ + Load an optional EBSD graph pickle handoff. + + Parameters + ---------- + path : str or pathlib.Path + Pickle file to read. + + Returns + ------- + dict + Graph handoff payload. + """ + with Path(path).open("rb") as file_handle: + return pickle.load(file_handle) + + +def build_ebsd_graph( + input_file: Union[str, Path], + config: Optional[EBSDGraphConfig] = None, + output_options: Optional[EBSDGraphOutputOptions] = None, +) -> EBSDGraphResult: + """ + Build an EBSD graph and return an in memory result object. + + Parameters + ---------- + input_file : str or pathlib.Path + EBSD input file used for graph construction. + config : EBSDGraphConfig, optional + Graph construction thresholds and connectivity. + output_options : EBSDGraphOutputOptions, optional + Optional file output settings. + + Raises + ------ + FileNotFoundError + If the EBSD input file does not exist. + + Returns + ------- + EBSDGraphResult + In memory graph result for downstream Python functions. + """ + graph_config = config or EBSDGraphConfig() + graph_output_options = output_options or EBSDGraphOutputOptions() + input_path = Path(input_file) + + if not input_path.exists(): + raise FileNotFoundError(f"EBSD input file not found: {input_path}") + + if graph_output_options.requires_output_dir: + output_dir = graph_output_options.resolved_output_dir() + if graph_output_options.clear_output_dir and output_dir.exists(): + shutil.rmtree(output_dir) + output_dir.mkdir(parents=True, exist_ok=True) + + start_time = time.perf_counter() + ebsd = EBSDmap( + str(input_path), + vf_min=graph_config.min_phase_pixel_fraction, + gs_min=graph_config.min_grain_size_pixels, + max_angle=graph_config.max_misorientation_degrees, + connectivity=graph_config.connectivity, + boundary_artifact_size_factor=graph_config.boundary_artifact_size_factor, + boundary_artifact_fraction_min=graph_config.boundary_artifact_fraction_min, + convex_hull_min_points=graph_config.convex_hull_min_points, + similar_orientation_merge_degrees=( + graph_config.similar_orientation_merge_degrees + ), + show_plot=False, + show_grains=False, + show_hist=False, + show_graph=False, + graph_only=True, + ) + + _, selected_phase_data = select_major_phase(ebsd) + elapsed_time_seconds = time.perf_counter() - start_time + result = EBSDGraphResult.from_phase_data( + phase_data=selected_phase_data, + npx=ebsd.npx, + sh_x=ebsd.sh_x, + sh_y=ebsd.sh_y, + source_input_file=input_path, + config=graph_config, + elapsed_time_seconds=elapsed_time_seconds, + ebsd=ebsd, + ) + + write_options = graph_output_options + if graph_output_options.clear_output_dir: + write_options = replace(graph_output_options, clear_output_dir=False) + + write_graph_result_outputs(result, write_options) + return result diff --git a/src/kanapy/texture.py b/src/kanapy/texture.py index 7b80a267..b5128de9 100644 --- a/src/kanapy/texture.py +++ b/src/kanapy/texture.py @@ -20,6 +20,7 @@ from scipy.integrate import quad from scipy.sparse import coo_matrix from scipy.sparse.csgraph import connected_components +from scipy.sparse.csgraph import connected_components as connected_components_csgraph from skimage.segmentation import mark_boundaries from orix import io from orix import plot as ox_plot @@ -33,6 +34,47 @@ from typing import Any, Sequence, Union, Dict, Tuple, List +MIN_GRAIN_SIZE_PIXELS = 10.0 +""" +Minimum number of pixels for a grain to be considered valid. +""" + +DEFAULT_MAX_MISORIENTATION_DEGREES = 5.0 +""" +Default maximum misorientation angle in degrees for grain segmentation. +""" + +DEFAULT_SIMILAR_ORIENTATION_MERGE_DEGREES = DEFAULT_MAX_MISORIENTATION_DEGREES +""" +Default maximum orientation difference for merging neighboring graph nodes. +""" + +DEFAULT_MIN_VOLUME_FRACTION = 0.03 +""" +Default minimum phase volume fraction for EBSD analysis. +""" + +DEFAULT_CONNECTIVITY = 4 +""" +Default pixel connectivity for graph construction. +""" + +DEFAULT_BOUNDARY_ARTIFACT_SIZE_FACTOR = 2.0 +""" +Default size factor for boundary artifact candidates. +""" + +DEFAULT_BOUNDARY_ARTIFACT_FRACTION_MIN = 0.85 +""" +Default minimum boundary pixel fraction for artifact candidates. +""" + +CONVEX_HULL_MIN_POINTS = 3 +""" +Minimum number of points required to compute a convex hull. +""" + + def get_distinct_colormap(N, cmap='prism'): """ Generate N visually distinct colors as an RGB colormap. @@ -162,7 +204,163 @@ def mean_orientation_data(pixel_orientations, sym): return mean_quat / np.linalg.norm(mean_quat) -def merge_nodes(G, node1, node2): +def get_proper_symmetry_quaternions(sym): + """ + Return proper symmetry operators as scalar first quaternions. + + Graph misorientation helpers use this compact array for vectorized + symmetry reduction. + + Parameters + ---------- + sym : object or numpy.ndarray + ORIX symmetry object or an existing quaternion array. + + Returns + ------- + numpy.ndarray + Proper symmetry operators with shape ``(n_ops, 4)``. + """ + if isinstance(sym, np.ndarray): + return np.asarray(sym, dtype=float) + + if getattr(sym, "contains_inversion", False): + sym = sym.proper_subgroup + + return np.asarray(sym.data, dtype=float) + + +def _quat_left_mul_batch(q, Q): + """ + Left multiply one quaternion with a batch of quaternions. + + Parameters + ---------- + q : numpy.ndarray + Single scalar first quaternion with shape ``(4,)``. + Q : numpy.ndarray + Batch of scalar first quaternions with shape ``(n, 4)``. + + Returns + ------- + numpy.ndarray + Quaternion products with shape ``(n, 4)``. + """ + w1, x1, y1, z1 = q[0], q[1], q[2], q[3] + w2, x2, y2, z2 = Q[:, 0], Q[:, 1], Q[:, 2], Q[:, 3] + return np.column_stack([ + w1 * w2 - x1 * x2 - y1 * y2 - z1 * z2, + w1 * x2 + x1 * w2 + y1 * z2 - z1 * y2, + w1 * y2 - x1 * z2 + y1 * w2 + z1 * x2, + w1 * z2 + x1 * y2 - y1 * x2 + z1 * w2, + ]) + + +def _batch_misorientation_angle(q_a, q_b, sym_ops): + """ + Return symmetry reduced misorientation angles for quaternion pairs. + + The implementation keeps memory bounded by looping over the first symmetry + side and vectorizing over the second. + + Parameters + ---------- + q_a : numpy.ndarray + First quaternion batch with shape ``(n, 4)``. + q_b : numpy.ndarray + Second quaternion batch with shape ``(n, 4)``. + sym_ops : numpy.ndarray + Proper symmetry operators with shape ``(n_ops, 4)``. + + Returns + ------- + numpy.ndarray + Misorientation angles in radians. + """ + q_a = np.asarray(q_a, dtype=float) + q_b = np.asarray(q_b, dtype=float) + sym_ops = np.asarray(sym_ops, dtype=float) + + norms_a = np.linalg.norm(q_a, axis=1, keepdims=True) + norms_b = np.linalg.norm(q_b, axis=1, keepdims=True) + q_a = q_a / np.where(norms_a > 0, norms_a, 1.0) + q_b = q_b / np.where(norms_b > 0, norms_b, 1.0) + + w1, x1, y1, z1 = q_b[:, 0], q_b[:, 1], q_b[:, 2], q_b[:, 3] + w2, x2, y2, z2 = q_a[:, 0], -q_a[:, 1], -q_a[:, 2], -q_a[:, 3] + raw = np.column_stack([ + w1 * w2 - x1 * x2 - y1 * y2 - z1 * z2, + w1 * x2 + x1 * w2 + y1 * z2 - z1 * y2, + w1 * y2 - x1 * z2 + y1 * w2 + z1 * x2, + w1 * z2 + x1 * y2 - y1 * x2 + z1 * w2, + ]) + + sign = np.array([1.0, -1.0, -1.0, -1.0]) + max_abs_w = np.zeros(len(q_a), dtype=np.float64) + for sym_left in sym_ops: + left = _quat_left_mul_batch(sym_left, raw) + candidate_w = (left * sign) @ sym_ops.T + np.maximum(max_abs_w, np.abs(candidate_w).max(axis=1), out=max_abs_w) + + np.clip(max_abs_w, 0.0, 1.0, out=max_abs_w) + return 2.0 * np.arccos(max_abs_w) + + +def _compute_all_edge_misorientations(G, sym_ops): + """ + Compute and cache misorientation values for every graph edge. + + Parameters + ---------- + G : networkx.Graph + Microstructure graph with average orientations on nodes. + sym_ops : numpy.ndarray + Proper symmetry operators used for misorientation reduction. + """ + edges = list(G.edges()) + if not edges: + return + + q_a = np.array([G.nodes[node_a]["ori_av"] for node_a, _node_b in edges]) + q_b = np.array([G.nodes[node_b]["ori_av"] for _node_a, node_b in edges]) + angles = _batch_misorientation_angle(q_a, q_b, sym_ops) + for (node_a, node_b), angle in zip(edges, angles): + G[node_a][node_b]["misorientation"] = float(angle) + + +def _refresh_node_misorientation_cache(G, node, sym_ops=None): + """ + Refresh cached misorientation values for one node's incident edges. + + Merging changes the surviving grain orientation, so cached values around + that node must be recomputed before later neighbor searches reuse them. + + Parameters + ---------- + G : networkx.Graph + Microstructure graph with average orientations on nodes. + node : int + Grain node whose incident edge cache values should be refreshed. + sym_ops : numpy.ndarray, optional + Proper symmetry operators. If omitted, they are built from graph + symmetry metadata. + """ + neighbors_sorted = sorted(G.adj[node]) + if not neighbors_sorted: + return + + if sym_ops is None: + sym_ops = get_proper_symmetry_quaternions(G.graph["symmetry"]) + + q_a = np.tile(G.nodes[node]["ori_av"], (len(neighbors_sorted), 1)) + q_b = np.array([G.nodes[neigh]["ori_av"] for neigh in neighbors_sorted]) + angles = _batch_misorientation_angle(q_a, q_b, sym_ops) + + for neigh, angle in zip(neighbors_sorted, angles): + G[node][neigh]["misorientation"] = float(angle) + + +def merge_nodes(G, node1, node2, preserve_target_orientation=False): """ Merge the pixel and attribute data of one node into another node. @@ -180,24 +378,32 @@ def merge_nodes(G, node1, node2): Node ID to be merged and removed. node2 : int Node ID to merge into and retain. + preserve_target_orientation : bool, optional + Keep the retained node orientation unchanged. Cleanup passes use this + when a small geometric artifact should not affect the larger node. Notes ----- This function modifies the graph in place. """ + has_misorientation_cache = any( + "misorientation" in data + for _node_a, _node_b, data in G.edges([node1, node2], data=True) + ) + # merge pixel lists before recomputing node attributes G.nodes[node2]['pixels'] = np.concatenate((G.nodes[node1]['pixels'], G.nodes[node2]['pixels'])) npix1 = G.nodes[node1]['npix'] npix2 = G.nodes[node2]['npix'] ntot = npix1 + npix2 - if 'rotations' in G.graph: + if not preserve_target_orientation and 'rotations' in G.graph: pix = G.nodes[node2]['pixels'] G.nodes[node2]['ori_av'] = mean_orientation_data( G.graph['rotations'][pix], G.graph['symmetry'], ) - else: + elif not preserve_target_orientation: # keep old graph objects usable when original pixel orientations are missing ori_av = (G.nodes[node1]['ori_av'] * G.nodes[node1]['npix'] + G.nodes[node2]['ori_av'] * G.nodes[node2]['npix']) / ntot @@ -211,10 +417,13 @@ def merge_nodes(G, node1, node2): # convert pixel distances to micron pts[:, 0] *= G.graph['dx'] pts[:, 1] *= G.graph['dy'] - G.nodes[node2]['hull'] = ConvexHull(pts) + if pts.shape[0] >= CONVEX_HULL_MIN_POINTS: + G.nodes[node2]['hull'] = ConvexHull(pts) + else: + G.nodes[node2].pop('hull', None) # add new edges and let the graph ignore existing ones - for neigh in G.adj[node1]: + for neigh in sorted(G.adj[node1]): if node2 != neigh: G.add_edge(node2, neigh) @@ -224,6 +433,9 @@ def merge_nodes(G, node1, node2): ix, iy = np.nonzero(G.graph['label_map'] == node1) G.graph['label_map'][ix, iy] = node2 + if has_misorientation_cache: + _refresh_node_misorientation_cache(G, node2) + def find_largest_neighbor(G, node): """ @@ -264,12 +476,12 @@ def find_largest_neighbor(G, node): return num_ln -def find_sim_neighbor(G, nn): +def find_sim_neighbor(G, nn, sym_ops=None): """ - Find the neighboring node most similar in orientation to the given node + Find the neighboring node most similar in orientation to the given node. This function computes the misorientation angles between the given node and - all its neighbors, and returns the neighbor with the smallest misorientation. + all its neighbors, or reads cached edge values when they are available. Parameters ---------- @@ -278,6 +490,9 @@ def find_sim_neighbor(G, nn): The graph also contains a 'symmetry' key in G.graph used for orientation calculations. nn : int Node ID for which the most similar neighbor is to be found. + sym_ops : numpy.ndarray, optional + Proper symmetry operators. Pass this when calling in a loop to avoid + rebuilding symmetry data. Returns ------- @@ -286,22 +501,31 @@ def find_sim_neighbor(G, nn): - neighbor_id : int, ID of the neighbor with smallest misorientation - angle : float, misorientation angle with the given node - Notes - ----- - - Uses the `Orientation` class for misorientation calculations - - Assumes all neighbors have a valid 'ori_av' attribute - - Does not modify the graph + Raises + ------ + ValueError + If the node has no neighbors. """ - sym = G.graph['symmetry'] - ori0 = Orientation(G.nodes[nn]['ori_av'], sym) - on = [] - ng = [] - for neigh in G.adj[nn]: - ang = ori0.angle_with(Orientation(G.nodes[neigh]['ori_av'], sym))[0] - on.append(ang) - ng.append(neigh) - nn = np.argmin(on) - return ng[nn], on[nn] + neighbors_sorted = sorted(G.adj[nn]) + if not neighbors_sorted: + raise ValueError(f"Grain {nn} has no neighbors.") + + if all("misorientation" in G[nn][neigh] for neigh in neighbors_sorted): + angles = [G[nn][neigh]["misorientation"] for neigh in neighbors_sorted] + else: + if sym_ops is None: + sym_ops = get_proper_symmetry_quaternions(G.graph["symmetry"]) + + q_a = np.tile(G.nodes[nn]["ori_av"], (len(neighbors_sorted), 1)) + q_b = np.array([G.nodes[neigh]["ori_av"] for neigh in neighbors_sorted]) + angles = _batch_misorientation_angle(q_a, q_b, sym_ops) + + rows = [] + for neigh, angle in zip(neighbors_sorted, angles): + rows.append((float(angle), int(neigh))) + + angle, neighbor = min(rows, key=lambda row: (row[0], row[1])) + return neighbor, angle def get_node_boundary_stats(G, node_id, connectivity=4): @@ -392,10 +616,12 @@ def merge_boundary_artifact_nodes( """ debug_records = [] npix_limit = max(int(npix_factor * gs_min), 40) + sym_ops = get_proper_symmetry_quaternions(G.graph['symmetry']) changed = True while changed: changed = False + _compute_all_edge_misorientations(G, sym_ops) candidate_rows = [] # collect small nodes that are mostly made of boundary pixels @@ -410,7 +636,7 @@ def merge_boundary_artifact_nodes( if stats['boundary_fraction'] < boundary_fraction_min: continue - sim_neigh, ori_neigh = find_sim_neighbor(G, node_id) + sim_neigh, ori_neigh = find_sim_neighbor(G, node_id, sym_ops=sym_ops) candidate_rows.append(( -stats['boundary_fraction'], npix, @@ -443,7 +669,7 @@ def merge_boundary_artifact_nodes( 'will_merge': True, } debug_records.append(record) - merge_nodes(G, node_id, sim_neigh) + merge_nodes(G, node_id, sim_neigh, preserve_target_orientation=True) changed = True break @@ -573,31 +799,55 @@ def build_graph_from_labeled_pixels(label_array, rot, sym, dx, dy, connectivity= - Loops over all pixels to add edges between neighboring grains - The function preserves the shape and indices of `label_array` """ - # t1 = time.time() nodes = summarize_labels(label_array, rot, sym) - # t2 = time.time() - # print(f'Time for extracting nodes: {t2 - t1}') - # print(f'Building microstructure graph with {len(nodes)} nodes (grains).') - G = nx.Graph(label_map=label_array, symmetry=sym, rotations=rot, + G = nx.Graph(label_map=label_array.copy(), symmetry=sym, rotations=rot, dx=dx, dy=dy) G.add_nodes_from(nodes) - # t3 = time.time() - # print(f'Time for building graph: {t3 - t2}') - - # print('Adding edges (grain boundaries) to microstructure graph.') - rows, cols = label_array.shape - for x in range(rows): - for y in range(cols): - label_here = label_array[x, y] - if label_here <= 0: - continue - for px, py in neighbors(x, y, connectivity=connectivity): - if 0 <= px < rows and 0 <= py < cols: - neighbor_label = label_array[px, py] - if neighbor_label > 0 and neighbor_label != label_here: - G.add_edge(label_here, neighbor_label) - # t4 = time.time() - # print(f'Time for adding edges: {t4 - t3}') + + edges = set() + + left = label_array[:, :-1].ravel() + right = label_array[:, 1:].ravel() + mask = left != right + if mask.any(): + pairs = np.stack((left[mask], right[mask]), axis=1) + for label_a, label_b in pairs: + if label_a > 0 and label_b > 0 and label_a != label_b: + edge = tuple(sorted((int(label_a), int(label_b)))) + edges.add(edge) + + top = label_array[:-1, :].ravel() + bottom = label_array[1:, :].ravel() + mask = top != bottom + if mask.any(): + pairs = np.stack((top[mask], bottom[mask]), axis=1) + for label_a, label_b in pairs: + if label_a > 0 and label_b > 0 and label_a != label_b: + edge = tuple(sorted((int(label_a), int(label_b)))) + edges.add(edge) + + if connectivity == 8: + top_left = label_array[:-1, :-1].ravel() + bottom_right = label_array[1:, 1:].ravel() + mask = top_left != bottom_right + if mask.any(): + pairs = np.stack((top_left[mask], bottom_right[mask]), axis=1) + for label_a, label_b in pairs: + if label_a > 0 and label_b > 0 and label_a != label_b: + edge = tuple(sorted((int(label_a), int(label_b)))) + edges.add(edge) + + top_right = label_array[:-1, 1:].ravel() + bottom_left = label_array[1:, :-1].ravel() + mask = top_right != bottom_left + if mask.any(): + pairs = np.stack((top_right[mask], bottom_left[mask]), axis=1) + for label_a, label_b in pairs: + if label_a > 0 and label_b > 0 and label_a != label_b: + edge = tuple(sorted((int(label_a), int(label_b)))) + edges.add(edge) + + G.add_edges_from(sorted(edges)) return G @@ -744,44 +994,61 @@ def find_similar_regions_by_misorientation(ori_map, phase_mask, sym, tolerance=0 raise ValueError("ori_map and phase_mask must have matching spatial dimensions.") rows, cols = phase_mask.shape - labeled_array = np.zeros((rows, cols), dtype=int) - visited = np.zeros((rows, cols), dtype=bool) - current_label = 1 - - for i in range(rows): - for j in range(cols): - if visited[i, j] or not phase_mask[i, j]: - continue - - # grow one region from this seed pixel - queue = deque([(i, j)]) - visited[i, j] = True - labeled_array[i, j] = current_label - - while queue: - x, y = queue.popleft() - ori_xy = Orientation(ori_map[x, y][None, :], sym) + sym_ops = get_proper_symmetry_quaternions(sym) + flat_index = np.arange(rows * cols, dtype=np.int32).reshape(rows, cols) + + directions = [ + (np.s_[:, :-1], np.s_[:, 1:]), + (np.s_[:-1, :], np.s_[1:, :]), + ] + if connectivity == 8: + directions += [ + (np.s_[:-1, :-1], np.s_[1:, 1:]), + (np.s_[:-1, 1:], np.s_[1:, :-1]), + ] + + source_rows = [] + target_rows = [] + for slice_a, slice_b in directions: + valid = phase_mask[slice_a] & phase_mask[slice_b] + if not np.any(valid): + continue - for xx, yy in neighbors(x, y, connectivity=connectivity): - if xx < 0 or xx >= rows or yy < 0 or yy >= cols: - continue - if visited[xx, yy] or not phase_mask[xx, yy]: - continue + q_a = ori_map[slice_a][valid] + q_b = ori_map[slice_b][valid] + same_region = _batch_misorientation_angle(q_a, q_b, sym_ops) <= tolerance + if not np.any(same_region): + continue - ori_nb = Orientation(ori_map[xx, yy][None, :], sym) - misorientation = ori_xy.angle_with(ori_nb)[0] + source_rows.append(flat_index[slice_a][valid][same_region]) + target_rows.append(flat_index[slice_b][valid][same_region]) - # only connect locally similar neighboring pixels - if misorientation > tolerance: - continue + pixel_count = rows * cols + if source_rows: + source = np.concatenate(source_rows) + target = np.concatenate(target_rows) + adjacency = coo_matrix( + (np.ones(len(source), dtype=np.float32), (source, target)), + shape=(pixel_count, pixel_count), + ) + adjacency = adjacency + adjacency.T + else: + adjacency = coo_matrix((pixel_count, pixel_count), dtype=np.float32) - queue.append((xx, yy)) - visited[xx, yy] = True - labeled_array[xx, yy] = current_label + _, labels_flat = connected_components_csgraph( + adjacency.tocsr(), + directed=False, + ) - current_label += 1 + labeled_array = np.zeros(pixel_count, dtype=int) + phase_flat = phase_mask.ravel() + phase_labels = labels_flat[phase_flat] + unique_components = np.unique(phase_labels) + remap = np.zeros(int(labels_flat.max()) + 1, dtype=int) + remap[unique_components] = np.arange(1, len(unique_components) + 1) + labeled_array[phase_flat] = remap[phase_labels] - return labeled_array, current_label - 1 + return labeled_array.reshape(rows, cols), len(unique_components) def calc_error(odf_ref, odf_test, res=10.): @@ -1908,6 +2175,15 @@ class EBSDmap: Default is 5.0. connectivity : int, optional Connectivity for grain identification. Default is 4. + boundary_artifact_size_factor : float, optional + Size factor used to identify small boundary artifact nodes. + boundary_artifact_fraction_min : float, optional + Minimum boundary pixel fraction for boundary artifact nodes. + convex_hull_min_points : int, optional + Minimum number of node points required for convex hull calculation. + similar_orientation_merge_degrees : float, optional + Maximum orientation difference in degrees for merging neighboring + graph nodes after invalid regions have been absorbed. show_plot : bool, optional If True, plots microstructure maps. Default is True. show_hist : bool, optional @@ -1918,6 +2194,9 @@ class EBSDmap: If True, plots grain labeling. Default is False. show_graph : bool, optional If True, plots graph nodes and edges on the EBSD IPF map. Default is False. + graph_only : bool, optional + If True, build the EBSD graph and skip grain statistics extraction. + Default is False. hist : bool, optional Deprecated. Use `show_hist` instead. plot : bool, optional @@ -1953,9 +2232,16 @@ class EBSDmap: Total number of grains after merging/pruning. """ - def __init__(self, fname, gs_min=10.0, vf_min=0.03, max_angle=5.0, connectivity=4, + def __init__(self, fname, gs_min=MIN_GRAIN_SIZE_PIXELS, + vf_min=DEFAULT_MIN_VOLUME_FRACTION, + max_angle=DEFAULT_MAX_MISORIENTATION_DEGREES, + connectivity=DEFAULT_CONNECTIVITY, + boundary_artifact_size_factor=DEFAULT_BOUNDARY_ARTIFACT_SIZE_FACTOR, + boundary_artifact_fraction_min=DEFAULT_BOUNDARY_ARTIFACT_FRACTION_MIN, + convex_hull_min_points=CONVEX_HULL_MIN_POINTS, + similar_orientation_merge_degrees=DEFAULT_SIMILAR_ORIENTATION_MERGE_DEGREES, show_plot=True, show_hist=None, felzenszwalb=False, show_grains=False, - show_graph=False, hist=None, plot=None): + show_graph=False, graph_only=False, hist=None, plot=None): def _reassign(pix, ori, phid, bads): @@ -2002,7 +2288,27 @@ def _reassign(pix, ori, phid, bads): logging.warning('Use of "hist" is depracted, use argument "show_plot" instead.') if show_hist is None: show_hist = show_plot - max_angle *= np.pi / 180 + if graph_only: + show_plot = False + show_hist = False + show_grains = False + show_graph = False + if convex_hull_min_points < CONVEX_HULL_MIN_POINTS: + raise ValueError( + f"convex_hull_min_points must be at least {CONVEX_HULL_MIN_POINTS}" + ) + segmentation_angle = max_angle * np.pi / 180 + similar_orientation_merge_angle = ( + similar_orientation_merge_degrees * np.pi / 180 + ) + self.gs_min = gs_min + self.vf_min = vf_min + self.max_angle = segmentation_angle + self.connectivity = connectivity + self.boundary_artifact_size_factor = boundary_artifact_size_factor + self.boundary_artifact_fraction_min = boundary_artifact_fraction_min + self.convex_hull_min_points = int(convex_hull_min_points) + self.similar_orientation_merge_angle = similar_orientation_merge_angle self.ms_data = [] # read EBSD map and return the orix object self.emap = io.load(fname) @@ -2145,7 +2451,7 @@ def _reassign(pix, ori, phid, bads): full_quat, phase_mask, data['cs'], - tolerance=max_angle, + tolerance=segmentation_angle, connectivity=connectivity, ) print(f"Phase #{data['index']} ({data['name']}): Identified Grains: {n_regions}") @@ -2172,6 +2478,8 @@ def _reassign(pix, ori, phid, bads): ms_graph, gs_min, connectivity=connectivity, + npix_factor=boundary_artifact_size_factor, + boundary_fraction_min=boundary_artifact_fraction_min, ) data["merge_debug"].extend(boundary_merge_records) print(f'After boundary-artifact pre-merge, {len(ms_graph.nodes)} grains left.') @@ -2179,6 +2487,8 @@ def _reassign(pix, ori, phid, bads): print('Starting to simplify microstructure graph.') + sym_ops = get_proper_symmetry_quaternions(ms_graph.graph['symmetry']) + _compute_all_edge_misorientations(ms_graph, sym_ops) grain_set = set(ms_graph.nodes) rem_grains = len(grain_set) while rem_grains > 0: @@ -2190,6 +2500,24 @@ def _reassign(pix, ori, phid, bads): pts[:, 0] *= self.dx pts[:, 1] *= self.dy + if pts.shape[0] < convex_hull_min_points: + num_ln = find_largest_neighbor(ms_graph, num) + data["merge_debug"].append({ + "node": int(num), + "reason": "no_hull", + "target": int(num_ln), + "node_npix": int(ms_graph.nodes[num]["npix"]), + "target_npix": int(ms_graph.nodes[num_ln]["npix"]), + "will_merge": True, + }) + merge_nodes( + ms_graph, + num, + num_ln, + preserve_target_orientation=True, + ) + continue + try: hull = ConvexHull(pts) ms_graph.nodes[num]['hull'] = hull @@ -2203,10 +2531,19 @@ def _reassign(pix, ori, phid, bads): "target_npix": int(ms_graph.nodes[num_ln]["npix"]), "will_merge": True, }) - merge_nodes(ms_graph, num, num_ln) + merge_nodes( + ms_graph, + num, + num_ln, + preserve_target_orientation=True, + ) continue - sim_neigh, ori_neigh = find_sim_neighbor(ms_graph, num) + sim_neigh, ori_neigh = find_sim_neighbor( + ms_graph, + num, + sym_ops=sym_ops, + ) data["merge_debug"].append({ "node": int(num), @@ -2214,13 +2551,13 @@ def _reassign(pix, ori, phid, bads): "target": int(sim_neigh), "angle_rad": float(ori_neigh), "angle_deg": float(np.degrees(ori_neigh)), - "threshold_deg": float(np.degrees(max_angle)), + "threshold_deg": float(np.degrees(similar_orientation_merge_angle)), "node_npix": int(ms_graph.nodes[num]["npix"]), "target_npix": int(ms_graph.nodes[sim_neigh]["npix"]), - "will_merge": bool(ori_neigh <= max_angle), + "will_merge": bool(ori_neigh <= similar_orientation_merge_angle), }) - if ori_neigh <= max_angle: + if ori_neigh <= similar_orientation_merge_angle: merge_nodes(ms_graph, num, sim_neigh) self.ngrains = len(ms_graph.nodes) @@ -2234,11 +2571,25 @@ def _reassign(pix, ori, phid, bads): rem_grains = len(grain_set) if ms_graph.nodes[num]['npix'] < gs_min: num_ln = find_largest_neighbor(ms_graph, num) - merge_nodes(ms_graph, num, num_ln) + merge_nodes( + ms_graph, + num, + num_ln, + preserve_target_orientation=True, + ) self.ngrains = len(ms_graph.nodes) data['ngrains'] = self.ngrains print(f'After elimination of small grains, {self.ngrains} grains left.') + if graph_only: + data['graph'] = ms_graph + self.ms_data.append(data) + print( + f"Built graph of phase #{data['index']} ({data['name']}) " + f"with {self.ngrains} grains." + ) + continue + # Extract grain statistics and axes arr_a = [] arr_b = [] diff --git a/tests/test_graph_operations.py b/tests/test_graph_operations.py new file mode 100644 index 00000000..6543cb7c --- /dev/null +++ b/tests/test_graph_operations.py @@ -0,0 +1,173 @@ +""" +Tests for EBSD graph helper functions. + +These tests cover graph misorientation caches and merge behavior used by the +EBSD cleanup workflow. +""" + +import networkx as nx +import numpy as np +from orix.crystal_map import Phase + +from kanapy.texture import ( + _batch_misorientation_angle, + _compute_all_edge_misorientations, + build_graph_from_labeled_pixels, + find_sim_neighbor, + get_proper_symmetry_quaternions, + merge_nodes, +) + + +def _make_orientation_graph(): + """ + Return a graph with four nodes at known misorientations. + + The graph is small enough to make direct and cached misorientation results + easy to compare. + + Returns + ------- + networkx.Graph + Test graph with average orientations and map metadata. + """ + sym = Phase(point_group="m-3m").point_group + graph = nx.Graph() + graph.graph["symmetry"] = sym + graph.graph["label_map"] = np.ones((8, 8), dtype=int) + graph.graph["dx"] = 1.0 + graph.graph["dy"] = 1.0 + + angles_deg = {1: 0.0, 2: 3.0, 3: 10.0, 4: 20.0} + for node_id, angle_deg in angles_deg.items(): + half_angle = np.deg2rad(angle_deg) / 2.0 + graph.add_node( + node_id, + npix=10, + pixels=np.arange((node_id - 1) * 10, node_id * 10), + ori_av=np.array( + [np.cos(half_angle), 0.0, 0.0, np.sin(half_angle)] + ), + ) + + graph.add_edge(1, 2) + graph.add_edge(1, 3) + graph.add_edge(2, 4) + graph.add_edge(3, 4) + return graph + + +def test_edge_misorientation_cache_matches_direct_batch_result(): + """ + Cached edge misorientation values match direct batch computations. + """ + graph = _make_orientation_graph() + sym_ops = get_proper_symmetry_quaternions(graph.graph["symmetry"]) + + expected = {} + for node_a, node_b in graph.edges(): + q_a = np.array([graph.nodes[node_a]["ori_av"]]) + q_b = np.array([graph.nodes[node_b]["ori_av"]]) + expected[(node_a, node_b)] = float( + _batch_misorientation_angle(q_a, q_b, sym_ops)[0] + ) + + _compute_all_edge_misorientations(graph, sym_ops) + + for node_a, node_b in graph.edges(): + assert "misorientation" in graph[node_a][node_b] + assert np.isclose( + graph[node_a][node_b]["misorientation"], + expected[(node_a, node_b)], + atol=1e-12, + ) + + +def test_find_sim_neighbor_uses_cached_values_without_changing_result(): + """ + Similar neighbor search returns the same result with or without cache. + """ + graph = _make_orientation_graph() + sym_ops = get_proper_symmetry_quaternions(graph.graph["symmetry"]) + + direct = { + node_id: find_sim_neighbor(graph, node_id, sym_ops=sym_ops) + for node_id in graph.nodes + } + _compute_all_edge_misorientations(graph, sym_ops) + cached = {node_id: find_sim_neighbor(graph, node_id) for node_id in graph.nodes} + + for node_id in graph.nodes: + assert direct[node_id][0] == cached[node_id][0] + assert np.isclose(direct[node_id][1], cached[node_id][1], atol=1e-12) + + +def test_merge_nodes_refreshes_surviving_node_misorientation_cache(): + """ + Merge updates cached edge misorientation values around the surviving node. + """ + graph = nx.Graph() + graph.graph["symmetry"] = Phase(point_group="m-3m").point_group + graph.graph["label_map"] = np.ones((6, 5), dtype=int) + graph.graph["dx"] = 1.0 + graph.graph["dy"] = 1.0 + + for node_id, angle_deg in [(1, 20.0), (2, 0.0), (3, 30.0)]: + half_angle = np.deg2rad(angle_deg) / 2.0 + graph.add_node( + node_id, + npix=10, + pixels=np.arange((node_id - 1) * 10, node_id * 10), + ori_av=np.array( + [np.cos(half_angle), 0.0, 0.0, np.sin(half_angle)] + ), + ) + + graph.add_edge(1, 2) + graph.add_edge(2, 3) + sym_ops = get_proper_symmetry_quaternions(graph.graph["symmetry"]) + _compute_all_edge_misorientations(graph, sym_ops) + + merge_nodes(graph, 1, 2) + + expected_angle = _batch_misorientation_angle( + np.array([graph.nodes[2]["ori_av"]]), + np.array([graph.nodes[3]["ori_av"]]), + sym_ops, + )[0] + + assert "misorientation" in graph[2][3] + assert np.isclose(graph[2][3]["misorientation"], expected_angle, atol=1e-12) + + +def test_build_graph_stores_independent_label_map_copy(): + """ + Graph metadata keeps its own label map copy. + + Downstream cleanup mutates graph labels during node merges, so callers must + not see their original label array modified through shared storage. + """ + label_map = np.array( + [ + [1, 1, 0], + [1, 2, 2], + [0, 2, 3], + ], + dtype=int, + ) + rotations = np.tile(np.array([[1.0, 0.0, 0.0, 0.0]]), (label_map.size, 1)) + symmetry = Phase(point_group="m-3m").point_group + + graph = build_graph_from_labeled_pixels( + label_map, + rotations, + symmetry, + dx=1.0, + dy=1.0, + connectivity=4, + ) + graph.graph["label_map"][0, 0] = 99 + + assert label_map[0, 0] == 1 + assert set(graph.nodes) == {1, 2, 3} + assert set(graph.edges) == {(1, 2), (2, 3)} diff --git a/tests/test_graph_workflow.py b/tests/test_graph_workflow.py new file mode 100644 index 00000000..e22e4f47 --- /dev/null +++ b/tests/test_graph_workflow.py @@ -0,0 +1,190 @@ +""" +Tests for the EBSD graph workflow objects. + +These tests keep the graph handoff in memory and make file outputs explicit. +""" + +import json + +import networkx as nx +import numpy as np + +from kanapy.graph_workflow import ( + EBSD_GRAPH_FILENAME, + GRAPH_IPF_MAP_FILENAME, + GRAPH_OVERLAY_FILENAME, + GRAPH_SUMMARY_FILENAME, + EBSDGraphConfig, + EBSDGraphOutputOptions, + EBSDGraphResult, + build_ebsd_graph, + write_graph_result_outputs, +) +import kanapy.graph_workflow as graph_workflow + + +def make_phase_data(): + """ + Build a minimal phase row for graph workflow tests. + + The row mirrors the fields produced by ``EBSDmap.ms_data`` that are needed + for plotting and serialization. + + Returns + ------- + dict + Minimal phase data containing a graph and an RGB image. + """ + graph = nx.Graph() + graph.graph["label_map"] = np.array([[1, 1], [2, 2]], dtype=int) + graph.graph["dx"] = 1.0 + graph.graph["dy"] = 1.0 + graph.add_node(1, pixels=np.array([0, 1]), npix=2) + graph.add_node(2, pixels=np.array([2, 3]), npix=2) + graph.add_edge(1, 2) + + return { + "index": 0, + "name": "Ni", + "vf": 1.0, + "graph": graph, + "rgb_im": np.array( + [ + [[1.0, 0.0, 0.0], [0.0, 1.0, 0.0]], + [[0.0, 0.0, 1.0], [1.0, 1.0, 0.0]], + ], + dtype=float, + ), + } + + +def test_graph_config_exposes_user_thresholds(): + """ + Graph configuration exposes defaults and user overrides. + """ + config = EBSDGraphConfig() + + assert config.min_phase_pixel_fraction == 0.03 + assert config.max_misorientation_degrees == 5.0 + assert config.min_grain_size_pixels == 10.0 + assert config.connectivity == 4 + assert config.boundary_artifact_size_factor == 2.0 + assert config.boundary_artifact_fraction_min == 0.85 + assert config.convex_hull_min_points == 3 + assert config.similar_orientation_merge_degrees == 5.0 + + custom = EBSDGraphConfig( + min_phase_pixel_fraction=0.05, + max_misorientation_degrees=4.0, + min_grain_size_pixels=25.0, + connectivity=8, + boundary_artifact_size_factor=3.0, + boundary_artifact_fraction_min=0.75, + convex_hull_min_points=4, + similar_orientation_merge_degrees=3.5, + ) + + assert custom.to_dict()["connectivity"] == 8 + assert custom.to_dict()["similar_orientation_merge_degrees"] == 3.5 + + +def test_graph_workflow_is_exported_from_kanapy(): + """ + Top level kanapy exports the graph workflow entry points. + """ + import kanapy as knpy + + assert knpy.EBSDGraphConfig is EBSDGraphConfig + assert knpy.EBSDGraphOutputOptions is EBSDGraphOutputOptions + assert knpy.EBSDGraphResult is EBSDGraphResult + assert callable(knpy.build_ebsd_graph) + + +def test_graph_outputs_are_explicit_and_exclude_step_figures(tmp_path): + """ + Selected graph outputs match the paper figure handoff files. + """ + result = EBSDGraphResult.from_phase_data( + phase_data=make_phase_data(), + npx=4, + sh_x=2, + sh_y=2, + source_input_file=tmp_path / "source.ang", + config=EBSDGraphConfig(), + ) + options = EBSDGraphOutputOptions( + output_dir=tmp_path, + write_pickle=True, + write_summary_json=True, + write_graph_overlay=True, + write_ipf_map=True, + write_step_figures=False, + ) + + written_files = write_graph_result_outputs(result, options) + + assert set(written_files) == { + "pickle", + "summary_json", + "graph_overlay", + "ipf_map", + } + assert written_files["pickle"].name == EBSD_GRAPH_FILENAME + assert written_files["summary_json"].name == GRAPH_SUMMARY_FILENAME + assert written_files["graph_overlay"].name == GRAPH_OVERLAY_FILENAME + assert written_files["ipf_map"].name == GRAPH_IPF_MAP_FILENAME + assert "step0" not in GRAPH_IPF_MAP_FILENAME + assert not list(tmp_path.glob("graph_step1*.png")) + + with written_files["summary_json"].open(encoding="utf-8") as file_handle: + summary = json.load(file_handle) + + assert summary["graph_node_count"] == 2 + assert summary["graph_edge_count"] == 1 + assert summary["graph_overlay"] == GRAPH_OVERLAY_FILENAME + assert summary["ipf_map"] == GRAPH_IPF_MAP_FILENAME + + +def test_build_ebsd_graph_uses_graph_only_map_analysis(tmp_path, monkeypatch): + """ + Graph workflow requests graph only EBSD analysis. + + This keeps the example focused on the graph handoff instead of running the + full statistical EBSD analysis path. + """ + calls = [] + + class FakeEBSDMap: + """ + Minimal EBSD map replacement for workflow tests. + """ + + def __init__(self, *_args, **kwargs): + """ + Record constructor options and expose one phase row. + + Parameters + ---------- + *_args : tuple + Positional arguments passed by the workflow. + **kwargs : dict + Keyword arguments passed by the workflow. + """ + calls.append(kwargs) + self.npx = 4 + self.sh_x = 2 + self.sh_y = 2 + self.ms_data = [make_phase_data()] + + input_file = tmp_path / "source.ang" + input_file.write_text("", encoding="utf-8") + monkeypatch.setattr(graph_workflow, "EBSDmap", FakeEBSDMap) + + result = build_ebsd_graph( + input_file=input_file, + config=EBSDGraphConfig(), + output_options=EBSDGraphOutputOptions(), + ) + + assert calls[0]["graph_only"] is True + assert result.graph.number_of_nodes() == 2