From ec71d60d76d64ab93e2a5e1031517ca81bf0e6c7 Mon Sep 17 00:00:00 2001 From: Saskia Gloy Date: Sun, 21 Dec 2025 11:19:06 +0100 Subject: [PATCH 01/15] Add surfaces to Discretization and update related documentation --- src/i2pp/core/discretization_helpers.py | 8 +++---- .../discretization_reader.py | 21 +++++++++++++++++++ .../fourc_yaml_reader.py | 13 ++++++++++++ .../discretization_readers/mesh_reader.py | 5 ++++- src/i2pp/core/export_data.py | 4 ++-- src/i2pp/core/interpolate_element_data.py | 2 +- src/i2pp/core/interpolators/interpolator.py | 2 +- .../interpolators/interpolator_all_voxel.py | 2 +- .../core/interpolators/interpolator_center.py | 8 +++---- .../core/interpolators/interpolator_nodes.py | 4 ++-- src/i2pp/core/visualize_results.py | 4 ++-- .../visualizers/discretization_visualizer.py | 4 ++-- .../test_fourc_yaml_reader.py | 12 +++++++++++ .../interpolators/test_interpolator_center.py | 4 ++-- .../interpolators/test_interpolator_nodes.py | 2 +- .../i2pp/core/test_discretization_helpers.py | 5 ++++- 16 files changed, 76 insertions(+), 24 deletions(-) diff --git a/src/i2pp/core/discretization_helpers.py b/src/i2pp/core/discretization_helpers.py index ecbda81..44e3157 100644 --- a/src/i2pp/core/discretization_helpers.py +++ b/src/i2pp/core/discretization_helpers.py @@ -94,8 +94,8 @@ def initialize_unstructured_grid( values. pixel_type (PixelValueType): The type of pixel values (e.g., RGB, MRT, CT). - dis (Discretization): The discretization object containing nodes and - elements. + dis (Discretization): The discretization object containing nodes, + elements and surfaces. Returns: tuple[pv.UnstructuredGrid, np.ndarray]: A tuple containing the PyVista `UnstructuredGrid` and a boolean array indicating which @@ -150,8 +150,8 @@ def get_elementwise_image_values( Arguments: elements_with_values (list[Element]): List of elements with assigned values. - dis (Discretization): The discretization object containing nodes and - elements. + dis (Discretization): The discretization object containing nodes, + elements and surfaces. pixel_type (PixelValueType): The type of pixel values (e.g., RGB, MRT, CT). Returns: diff --git a/src/i2pp/core/discretization_readers/discretization_reader.py b/src/i2pp/core/discretization_readers/discretization_reader.py index 2b94c49..67243f2 100644 --- a/src/i2pp/core/discretization_readers/discretization_reader.py +++ b/src/i2pp/core/discretization_readers/discretization_reader.py @@ -62,6 +62,23 @@ class Element: data: Optional[np.ndarray] = None +@dataclass +class Surface: + """Class for storing information about a surface in the Discretization. + + This class represents a surface in the mesh Discretization, defined + by its node IDs, an surface ID, and associated data. + + Attributes: + node_ids (np.ndarray): An array of node IDs that define the nodes of + the surface. + id (int): A unique identifier for the surface. + """ + + node_ids: np.ndarray + id: int + + @dataclass class Discretization: """Class for storing Discretization data. @@ -77,12 +94,16 @@ class Discretization: elements (list[Element]): A list of elements, each representing a part of the Discretization, containing information such as node IDs, element ID, center coordinates, and element data. + surfaces (list[Surface]): A list of surfaces, each representing a part + of the Discretization, containing information such as node IDs + and surface ID. bounding_box (Optional[np.ndarray]): Boundary limits of the Discretization. """ nodes: Nodes elements: list[Element] + surfaces: list[Surface] bounding_box: Optional[BoundingBox] = None diff --git a/src/i2pp/core/discretization_readers/fourc_yaml_reader.py b/src/i2pp/core/discretization_readers/fourc_yaml_reader.py index 24242ea..da0179c 100644 --- a/src/i2pp/core/discretization_readers/fourc_yaml_reader.py +++ b/src/i2pp/core/discretization_readers/fourc_yaml_reader.py @@ -10,6 +10,7 @@ DiscretizationReader, Element, Nodes, + Surface, ) from lnmmeshio import Discretization as FourCDiscretization from tqdm import tqdm @@ -119,9 +120,21 @@ def load_discretization( Element(node_ids=np.array(ele_node_ids), id=ele.id) ) + surfaces = [] + + for surf in raw_dis.surfacenodesets: + surf_node_ids = [] + for node in surf.nodes: + surf_node_ids.append(node.id) + + surfaces.append( + Surface(node_ids=np.array(surf_node_ids), id=surf.id) + ) + dis = Discretization( nodes=Nodes(coords=np.array(nodes_coords), ids=np.array(node_ids)), elements=elements, + surfaces=surfaces, ) return dis diff --git a/src/i2pp/core/discretization_readers/mesh_reader.py b/src/i2pp/core/discretization_readers/mesh_reader.py index 603dcd2..09a9848 100644 --- a/src/i2pp/core/discretization_readers/mesh_reader.py +++ b/src/i2pp/core/discretization_readers/mesh_reader.py @@ -65,6 +65,9 @@ def load_discretization( for i, face in enumerate(raw_dis.faces): elements.append(Element(node_ids=face, id=i)) - dis = Discretization(nodes=nodes, elements=elements) + # TODO: Add surfaces if needed in the future + # surfaces = [] + + dis = Discretization(nodes=nodes, elements=elements, surfaces=[]) return dis diff --git a/src/i2pp/core/export_data.py b/src/i2pp/core/export_data.py index 5a7c333..44d08f5 100644 --- a/src/i2pp/core/export_data.py +++ b/src/i2pp/core/export_data.py @@ -82,8 +82,8 @@ def export_data( Arguments: transformed_data (Any): The result from the transform_data function. elements (List[Element]): List of elements with IDs and data. - dis (Discretization): The discretization object containing nodes - and elements. + dis (Discretization): The discretization object containing nodes, + elements and surfaces. export_format (str): The format in which the data will be exported (e.g., "json", "txt"). property_output_file (Path): Path to the output file where the diff --git a/src/i2pp/core/interpolate_element_data.py b/src/i2pp/core/interpolate_element_data.py index 580ed41..b00f830 100644 --- a/src/i2pp/core/interpolate_element_data.py +++ b/src/i2pp/core/interpolate_element_data.py @@ -27,7 +27,7 @@ def interpolate_image_to_discretization( Arguments: dis (Discretization): The Discretization object containing FEM - elements and node coordinates. + surfaces, elements and node coordinates. image_data (ImageData): A structured representation containing 3D pixel data, grid coordinates, orientation, and metadata. interpolation_method (str): The type of interpolation to perform for diff --git a/src/i2pp/core/interpolators/interpolator.py b/src/i2pp/core/interpolators/interpolator.py index bc64af8..e732ae2 100644 --- a/src/i2pp/core/interpolators/interpolator.py +++ b/src/i2pp/core/interpolators/interpolator.py @@ -142,7 +142,7 @@ def compute_element_data( Arguments: dis (Discretization): The Discretization object containing FEM - elements and node coordinates. + surfaces, elements and node coordinates. image_data (ImageData): The 3D image data containing voxel coordinates and intensity values. diff --git a/src/i2pp/core/interpolators/interpolator_all_voxel.py b/src/i2pp/core/interpolators/interpolator_all_voxel.py index e251c1e..3b6447c 100644 --- a/src/i2pp/core/interpolators/interpolator_all_voxel.py +++ b/src/i2pp/core/interpolators/interpolator_all_voxel.py @@ -174,7 +174,7 @@ def compute_element_data( Arguments: dis (Discretization): The Discretization object containing FEM - elements and node coordinates. + surfaces, elements and node coordinates. image_data (ImageData): A structured representation containing 3D pixel data, grid coordinates, orientation, and metadata. diff --git a/src/i2pp/core/interpolators/interpolator_center.py b/src/i2pp/core/interpolators/interpolator_center.py index 5408ebc..53b688d 100644 --- a/src/i2pp/core/interpolators/interpolator_center.py +++ b/src/i2pp/core/interpolators/interpolator_center.py @@ -30,8 +30,8 @@ def compute_element_centers(self, dis: Discretization) -> Discretization: element. Arguments: - dis (Discretization): A finite element mesh containing elements - and their associated nodes. + dis (Discretization): A finite element mesh containing surfaces, + elements and their associated nodes. Returns: Discretization: The input Discretization object with updated @@ -66,8 +66,8 @@ def compute_element_data( retrieves the corresponding pixel values through interpolation. Arguments: - dis (Discretization): A finite element mesh containing elements - and their associated nodes. + dis (Discretization): A finite element mesh containing surfaces + elements and their associated nodes. image_data (ImageData): A structured representation of 3D image data, including pixel intensities, grid coordinates, and orientation. diff --git a/src/i2pp/core/interpolators/interpolator_nodes.py b/src/i2pp/core/interpolators/interpolator_nodes.py index 9320dda..584ad39 100644 --- a/src/i2pp/core/interpolators/interpolator_nodes.py +++ b/src/i2pp/core/interpolators/interpolator_nodes.py @@ -34,8 +34,8 @@ def compute_element_data( representative value to the element. Arguments: - dis (Discretization): The FEM discretization containing elements - and node coordinate data. + dis (Discretization): The FEM discretization containing surfaces, + elements and node coordinate data. image_data (ImageData): The 3D image dataset, including voxel values, spatial positioning, and metadata. diff --git a/src/i2pp/core/visualize_results.py b/src/i2pp/core/visualize_results.py index cd620e1..18c7a50 100644 --- a/src/i2pp/core/visualize_results.py +++ b/src/i2pp/core/visualize_results.py @@ -54,8 +54,8 @@ def visualize_results( assigned values. image_data (ImageData): Image data containing pixel values and grid coordinates. - dis (Discretization): The discretization object containing nodes - and elements. + dis (Discretization): The discretization object containing nodes, + elements and surfaces. Returns: None """ diff --git a/src/i2pp/core/visualizers/discretization_visualizer.py b/src/i2pp/core/visualizers/discretization_visualizer.py index 2cf5f5f..83a9658 100644 --- a/src/i2pp/core/visualizers/discretization_visualizer.py +++ b/src/i2pp/core/visualizers/discretization_visualizer.py @@ -43,8 +43,8 @@ def compute_grid( Arguments: elements_with_values (list[Element]): List of elements with assigned values. - dis (Discretization): The discretization object containing nodes - and elements. + dis (Discretization): The discretization object containing nodes, + elements and surfaces. Returns: None """ diff --git a/tests/unittests/i2pp/core/discretization_readers/test_fourc_yaml_reader.py b/tests/unittests/i2pp/core/discretization_readers/test_fourc_yaml_reader.py index afafb97..7a8b7a4 100644 --- a/tests/unittests/i2pp/core/discretization_readers/test_fourc_yaml_reader.py +++ b/tests/unittests/i2pp/core/discretization_readers/test_fourc_yaml_reader.py @@ -77,6 +77,10 @@ def test_load_discretization_fourc_yaml_without_filter(tmp_path: Path) -> None: ele1 = MagicMock(nodes=[node3, node1]) ele2 = MagicMock(nodes=[node2, node4]) + surface1 = MagicMock(nodes=[node1, node2]) + surface2 = MagicMock(nodes=[node3, node4]) + + mock_dis.surfacenodesets = [surface1, surface2] mock_dis.elements.structure = [ele1, ele2] mock_dis.nodes = [node1, node2, node3, node4] @@ -96,3 +100,11 @@ def test_load_discretization_fourc_yaml_without_filter(tmp_path: Path) -> None: assert np.array_equal( dis_loaded.elements[1].node_ids, np.array([2, 4]) ) + + assert np.array_equal( + dis_loaded.surfaces[0].node_ids, np.array([1, 2]) + ) + + assert np.array_equal( + dis_loaded.surfaces[1].node_ids, np.array([3, 4]) + ) diff --git a/tests/unittests/i2pp/core/interpolators/test_interpolator_center.py b/tests/unittests/i2pp/core/interpolators/test_interpolator_center.py index 57aebba..8d58008 100644 --- a/tests/unittests/i2pp/core/interpolators/test_interpolator_center.py +++ b/tests/unittests/i2pp/core/interpolators/test_interpolator_center.py @@ -31,7 +31,7 @@ def test_get_center(): nodes = Nodes(nodes_coords, node_ids) - test_dis = Discretization(nodes=nodes, elements=elements) + test_dis = Discretization(nodes=nodes, elements=elements, surfaces=[]) interpol = InterpolatorCenter() @@ -53,7 +53,7 @@ def test_compute_element_data(): elements = [element1, element2] nodes = Nodes(coords=node_coords, ids=[0, 1, 2, 3]) - dis = Discretization(nodes=nodes, elements=elements) + dis = Discretization(nodes=nodes, elements=elements, surfaces=[]) pixel_data = np.random.randint(0, 256, size=(4, 4, 4, 3)) diff --git a/tests/unittests/i2pp/core/interpolators/test_interpolator_nodes.py b/tests/unittests/i2pp/core/interpolators/test_interpolator_nodes.py index 7dd6aaa..6fb6d45 100644 --- a/tests/unittests/i2pp/core/interpolators/test_interpolator_nodes.py +++ b/tests/unittests/i2pp/core/interpolators/test_interpolator_nodes.py @@ -25,7 +25,7 @@ def test_compute_element_data(): elements = [element1, element2] nodes = Nodes(coords=node_coords, ids=[0, 1, 2, 3]) - dis = Discretization(nodes=nodes, elements=elements) + dis = Discretization(nodes=nodes, elements=elements, surfaces=[]) pixel_data = np.random.randint(0, 256, size=(4, 4, 4, 3)) diff --git a/tests/unittests/i2pp/core/test_discretization_helpers.py b/tests/unittests/i2pp/core/test_discretization_helpers.py index bcfe5f8..632d19d 100644 --- a/tests/unittests/i2pp/core/test_discretization_helpers.py +++ b/tests/unittests/i2pp/core/test_discretization_helpers.py @@ -75,7 +75,7 @@ def test_verify_and_load_discretization(): options_dict = {} mock_bounding = tuple([[0, 0, 0], [1, 1, 1]]) nodes = Nodes([0, 0, 0], 0) - mock_dis = Discretization(nodes, []) + mock_dis = Discretization(nodes, [], []) with patch( "i2pp.core.discretization_helpers.determine_discretization_format", @@ -122,6 +122,7 @@ def test_initialize_unstructured_grid(node_ids, expected_cell_type): elements=[ Element(id=1, node_ids=node_ids), ], + surfaces=[], ) mock_elements_with_values = [ Element(node_ids=node_ids, id=1, data=np.array([255, 0, 0])), @@ -195,6 +196,7 @@ def test_get_elementwise_image_values_rgb(): Element(id=2, node_ids=[1, 2]), Element(id=3, node_ids=[1, 2]), ], + surfaces=[], ) pixel_type = PixelValueType.RGB values, ele_has_value = get_elementwise_image_values( @@ -222,6 +224,7 @@ def test_get_elementwise_image_values_mrt(): Element(id=2, node_ids=[1, 2]), Element(id=3, node_ids=[1, 2]), ], + surfaces=[], ) pixel_type = PixelValueType.MRT values, ele_has_value = get_elementwise_image_values( From b7e93ae65e5e17c412cf04ec77bd942b3bdb4205 Mon Sep 17 00:00:00 2001 From: Saskia Gloy Date: Mon, 29 Dec 2025 10:26:17 +0100 Subject: [PATCH 02/15] Enhance smoothing function with masking support --- src/i2pp/core/run.py | 9 +- src/i2pp/core/utilities.py | 207 ++++++++++++++++++++++++++++++++++++- 2 files changed, 211 insertions(+), 5 deletions(-) diff --git a/src/i2pp/core/run.py b/src/i2pp/core/run.py index 3eda9d5..d126f20 100644 --- a/src/i2pp/core/run.py +++ b/src/i2pp/core/run.py @@ -13,7 +13,7 @@ interpolate_image_to_discretization, ) from i2pp.core.transform_data import transform_data -from i2pp.core.utilities import smooth_data +from i2pp.core.utilities import create_mesh_mask, smooth_data from i2pp.core.visualize_results import visualize_results, visualize_smoothing @@ -54,6 +54,9 @@ def run_i2pp(config_i2pp): discretization.bounding_box, ) + # Build exact mesh mask + image.mask = create_mesh_mask(discretization, image) + # If smoothing is enabled, smooth the image data if config.processing.smoothing: # Create a copy of the image data for visualization purposes @@ -61,7 +64,9 @@ def run_i2pp(config_i2pp): image_raw = copy.deepcopy(image) image.pixel_data = smooth_data( - image.pixel_data, config.processing.smoothing.smoothing_area + image.pixel_data, + config.processing.smoothing.smoothing_area, + mask=image.mask, ) if config.processing.smoothing.visualize: diff --git a/src/i2pp/core/utilities.py b/src/i2pp/core/utilities.py index f4690f3..a484835 100644 --- a/src/i2pp/core/utilities.py +++ b/src/i2pp/core/utilities.py @@ -5,6 +5,7 @@ import numpy as np from scipy.ndimage import uniform_filter +from scipy.spatial import ConvexHull def find_mins_maxs( @@ -90,6 +91,7 @@ def get_node_position_of_element( def smooth_data( data: np.ndarray, smoothing_window: int, + mask: np.ndarray, ) -> np.ndarray: """Applies a smoothing filter to 3D image data by averaging pixel values. @@ -98,23 +100,73 @@ def smooth_data( within a neighborhood defined by the `smoothing_window` parameter, which helps to smooth out irregularities in the data. + If a mask is provided, performs masked smoothing by: + - Zero data outside mask → data_masked + - Smooth the zeroed data → smoothed_data (artificially darkened at edges) + - Smooth the mask → smoothed_mask (tells you the "weight" of valid data) + - Divide smoothed data by smoothed mask → This "undoes" the darkening by renormalizing + - Only update pixels that were originally inside the mask → Preserves original data outside + Args: data (np.ndarray): A 3D array containing the pixel data to be smoothed. smoothing_window (int): The size of the neighborhood (in points) used to compute the average. Larger values result in smoother data, but may reduce fine details. + mask (np.ndarray): A boolean mask array of the same shape as `data`. + If provided, smoothing is only applied within the + masked region. Returns: np.ndarray: The smoothed image data as a 3D array, where each pixel's value has been replaced by the average of its neighbors within the defined smoothing window. """ - logging.info("Smooth data!") - return uniform_filter( - data, size=smoothing_window, mode="nearest", axes=(0, 1, 2) + if mask is None or not np.any(mask): + logging.warning( + "No mask provided for smoothing. " + "Smoothing will be applied to entire data." + ) + return uniform_filter( + data, size=smoothing_window, mode="nearest", axes=(0, 1, 2) + ) + + orig_dtype = data.dtype + data_f = data.astype(np.float32, copy=False) + + # Create a copy of the data where all values outside the mask are zero. + data_masked = data_f.copy() + data_masked[~mask, ...] = 0 + + # Apply a uniform filter to data_masked, which darkens values near edges. + smoothed_data = uniform_filter( + data_masked, size=smoothing_window, mode="nearest", axes=(0, 1, 2) ) + # Apply the same filter to the mask to find the proportion of valid data. + smoothed_mask = uniform_filter( + mask.astype(np.float32), + size=smoothing_window, + mode="nearest", + axes=(0, 1, 2), + ) + + smoothed_mask_exp = smoothed_mask[..., None] + valid_mask = (mask) & (smoothed_mask > 1e-8) + result_f = data_f.copy() + + normalized = np.zeros_like(smoothed_data, dtype=np.float32) + # Normalize smoothed data to correct darkening at the edges. + np.divide( + smoothed_data, + smoothed_mask_exp, + out=normalized, + where=smoothed_mask_exp > 1e-8, + ) + # Update pixels inside the original mask with the smoothed values. + result_f[valid_mask, ...] = normalized[valid_mask, ...] + + return result_f.astype(orig_dtype, copy=False) def make_json_serializable(obj: Any) -> Any: @@ -148,3 +200,152 @@ def make_json_serializable(obj: Any) -> Any: elif isinstance(obj, (np.float64, np.float32, np.float16)): return float(obj) return obj + + +def world_to_grid_coords( + points_world: np.ndarray, orientation: np.ndarray, position: np.ndarray +) -> np.ndarray: + """Convert world coordinates to image grid coordinates using image + orientation and origin. + + Args: + points_world (np.ndarray): A NumPy array of shape (N, 3) representing + N points in world coordinates. + orientation (np.ndarray): A 3x3 matrix representing the orientation + of the image grid in world space. + position (np.ndarray): A 1D array of shape (3,) representing the origin + of the image grid in world space. + + Returns: + np.ndarray: A NumPy array of shape (N, 3) representing the points in + image grid coordinates. + """ + return np.linalg.solve(orientation, (points_world - position).T).T + + +def _grid_bbox_indices_for_points( + grid_points: np.ndarray, grid_coords +) -> Tuple[int, int, int, int, int, int]: + """Compute axis-aligned grid bounding box indices for given grid points. + + Args: + grid_points (np.ndarray): A NumPy array of shape (N, 3) + representing N points in grid coordinates. + grid_coords: An object containing the grid's coordinate + arrays (slice, row, col). + + Returns: + Tuple[int, int, int, int, int, int]: A tuple of six integers representing + the bounding box indices (z0, z1, y0, y1, x0, x1), where + - z0, z1: Start and end indices along the Z-axis. + - y0, y1: Start and end indices along the Y-axis. + - x0, x1: Start and end indices along the X-axis. + """ + mins = np.min(grid_points, axis=0) + maxs = np.max(grid_points, axis=0) + + # Use np.clip to avoid redundant max/min calls + z0 = np.clip( + np.searchsorted(grid_coords.slice, mins[0], side="left"), + 0, + len(grid_coords.slice), + ) + z1 = np.clip( + np.searchsorted(grid_coords.slice, maxs[0], side="right"), + 0, + len(grid_coords.slice), + ) + y0 = np.clip( + np.searchsorted(grid_coords.row, mins[1], side="left"), + 0, + len(grid_coords.row), + ) + y1 = np.clip( + np.searchsorted(grid_coords.row, maxs[1], side="right"), + 0, + len(grid_coords.row), + ) + x0 = np.clip( + np.searchsorted(grid_coords.col, mins[2], side="left"), + 0, + len(grid_coords.col), + ) + x1 = np.clip( + np.searchsorted(grid_coords.col, maxs[2], side="right"), + 0, + len(grid_coords.col), + ) + + return z0, z1, y0, y1, x0, x1 + + +def create_mesh_mask(discretization, image) -> np.ndarray: + """Create a voxel-based mask of the discretization in the image's grid + space. + + This function generates a 3D boolean mask that represents a discretized + mesh within the grid space of a given image. + The mask is computed by determining which voxels in the image grid are + enclosed by the convex hulls of the mesh elements. + + Args: + discretization: A discretized mesh object containing nodes, elements + and surfaces. + image: An image object that provides the grid space and pixel data. + + Returns: + np.ndarray: A 3D boolean array where `True` indicates that the + corresponding voxel is part of the mesh mask, and `False` otherwise. + + Notes: + - The function uses convex hulls to approximate the spatial extent of + each mesh element in the image grid space. + - The convex hull computation requires at least three points. Elements + with fewer than three nodes are skipped. + - The mask is constructed iteratively by updating the relevant voxels + for each mesh element. + """ + logging.info("Create mesh mask in image grid space") + + mask = np.zeros(image.pixel_data.shape[:3], dtype=bool) + + node_id_to_idx = {nid: i for i, nid in enumerate(discretization.nodes.ids)} + nodes_grid = world_to_grid_coords( + discretization.nodes.coords, image.orientation, image.position + ) + + for ele in discretization.elements: + node_idxs = np.fromiter( + (node_id_to_idx[nid] for nid in ele.node_ids), + dtype=int, + count=len(ele.node_ids), + ) + ele_node_grid = nodes_grid[node_idxs] + + if ele_node_grid.shape[0] < 3: + continue + + try: + hull = ConvexHull(ele_node_grid) + except Exception: + continue + + z0, z1, y0, y1, x0, x1 = _grid_bbox_indices_for_points( + ele_node_grid, image.grid_coords + ) + + zz, yy, xx = np.meshgrid( + image.grid_coords.slice[z0:z1], + image.grid_coords.row[y0:y1], + image.grid_coords.col[x0:x1], + indexing="ij", + ) + grid_points = np.stack((zz.ravel(), yy.ravel(), xx.ravel()), axis=-1) + + A, b = hull.equations[:, :-1], hull.equations[:, -1] + inside_mask = np.all(A @ grid_points.T + b[:, None] <= 0, axis=0) + + inside_mask = inside_mask.reshape(zz.shape) + mask[z0:z1, y0:y1, x0:x1] |= inside_mask + + return mask From 838a7c45f78dc42bb6b94a194abf26ed4e2837bd Mon Sep 17 00:00:00 2001 From: Saskia Gloy Date: Mon, 29 Dec 2025 11:57:35 +0100 Subject: [PATCH 03/15] Add tests for create_mesh_mask and update smoothing function calls --- .../test_data/expected_physical_property.json | 1834 ++++++++--------- tests/unittests/i2pp/core/test_run.py | 11 +- tests/unittests/i2pp/core/test_utilities.py | 95 +- 3 files changed, 1019 insertions(+), 921 deletions(-) diff --git a/tests/integration/test_data/expected_physical_property.json b/tests/integration/test_data/expected_physical_property.json index 5d4b811..9ff6c86 100644 --- a/tests/integration/test_data/expected_physical_property.json +++ b/tests/integration/test_data/expected_physical_property.json @@ -1,1004 +1,1004 @@ { "STIFFNESS": { - "1": 0.4532679617404938, - "2": 0.47418299317359924, - "3": 0.4960784316062927, - "4": 0.503104567527771, - "5": 0.5076797604560852, - "6": 0.5302287340164185, + "1": 0.45457515120506287, + "2": 0.4749999940395355, + "3": 0.49656862020492554, + "4": 0.5034313797950745, + "5": 0.5078431367874146, + "6": 0.5303921699523926, "7": 0.0, - "8": 0.4794117510318756, - "9": 0.4426470696926117, - "10": 0.3562091588973999, - "11": 0.45130717754364014, - "12": 0.5006536245346069, - "13": 0.5472221970558167, - "14": 0.5300653576850891, - "15": 0.5346405506134033, - "16": 0.5800653696060181, + "8": 0.4803921580314636, + "9": 0.4436274468898773, + "10": 0.3570261299610138, + "11": 0.45228758454322815, + "12": 0.5016340017318726, + "13": 0.5480391979217529, + "14": 0.5303921699523926, + "15": 0.5352941155433655, + "16": 0.5807189345359802, "17": 0.0, - "18": 0.506209135055542, - "19": 0.46209150552749634, - "20": 0.45261436700820923, - "21": 0.4624182879924774, - "22": 0.5013071894645691, + "18": 0.5071895718574524, + "19": 0.46339869499206543, + "20": 0.4539215564727783, + "21": 0.46339869499206543, + "22": 0.5026143789291382, "23": 0.0, - "24": 0.5148692727088928, - "25": 0.5009803771972656, - "26": 0.5215686559677124, - "27": 0.4960784316062927, - "28": 0.46388888359069824, - "29": 0.48643791675567627, - "30": 0.5375816822052002, - "31": 0.5075163245201111, - "32": 0.5408496856689453, - "33": 0.5495098233222961, - "34": 0.49836599826812744, + "24": 0.5155228972434998, + "25": 0.5021241903305054, + "26": 0.522549033164978, + "27": 0.4964052140712738, + "28": 0.46437907218933105, + "29": 0.4870915114879608, + "30": 0.5383986830711365, + "31": 0.5088235139846802, + "32": 0.5423202514648438, + "33": 0.5511437654495239, + "34": 0.4995098114013672, "35": 0.0, - "36": 0.49084967374801636, - "37": 0.44281044602394104, - "38": 0.48186275362968445, - "39": 0.5519607663154602, - "40": 0.5565359592437744, - "41": 0.5029411911964417, - "42": 0.535620927810669, - "43": 0.5277777910232544, - "44": 0.4732026159763336, + "36": 0.4919934570789337, + "37": 0.4434640407562256, + "38": 0.48235294222831726, + "39": 0.5527777671813965, + "40": 0.5571895241737366, + "41": 0.5045751929283142, + "42": 0.5367646813392639, + "43": 0.52875816822052, + "44": 0.47418299317359924, "45": 0.0, "46": 0.0, - "47": 0.529411792755127, - "48": 0.5439542531967163, - "49": 0.5870915055274963, - "50": 0.5800653696060181, - "51": 0.4816993474960327, - "52": 0.48986926674842834, - "53": 0.48905229568481445, - "54": 0.4740196168422699, - "55": 0.5080065131187439, + "47": 0.5299019813537598, + "48": 0.5444444417953491, + "49": 0.5882353186607361, + "50": 0.5807189345359802, + "51": 0.4830065369606018, + "52": 0.4906862676143646, + "53": 0.48954248428344727, + "54": 0.47467321157455444, + "55": 0.5086601376533508, "56": 0.0, - "57": 0.5986928343772888, - "58": 0.5367646813392639, - "59": 0.5300653576850891, - "60": 0.578758180141449, - "61": 0.4771241843700409, - "62": 0.46977123618125916, - "63": 0.4781045615673065, - "64": 0.5006536245346069, - "65": 0.49885621666908264, - "66": 0.5274509787559509, - "67": 0.5467320084571838, - "68": 0.4797385632991791, - "69": 0.4575163424015045, - "70": 0.4977124035358429, - "71": 0.42549020051956177, - "72": 0.4120914936065674, - "73": 0.41192811727523804, - "74": 0.4377450942993164, - "75": 0.45996731519699097, - "76": 0.48758170008659363, - "77": 0.4874182939529419, - "78": 0.4472222328186035, - "79": 0.4320261478424072, - "80": 0.44068628549575806, - "81": 0.35947713255882263, - "82": 0.349346399307251, - "83": 0.36666667461395264, - "84": 0.3928104639053345, - "85": 0.4220588207244873, - "86": 0.48267972469329834, - "87": 0.5037581920623779, - "88": 0.4905228614807129, - "89": 0.47058823704719543, - "90": 0.44558823108673096, - "91": 0.32679739594459534, - "92": 0.32647058367729187, - "93": 0.37385621666908264, - "94": 0.4117647111415863, - "95": 0.43039214611053467, + "57": 0.5988562107086182, + "58": 0.5377451181411743, + "59": 0.5312091708183289, + "60": 0.5792483687400818, + "61": 0.47777777910232544, + "62": 0.4704248309135437, + "63": 0.4785947799682617, + "64": 0.5013071894645691, + "65": 0.4995098114013672, + "66": 0.5279411673545837, + "67": 0.5475490093231201, + "68": 0.48104575276374817, + "69": 0.4586601257324219, + "70": 0.4982026219367981, + "71": 0.42679738998413086, + "72": 0.41290849447250366, + "73": 0.4125817120075226, + "74": 0.43921568989753723, + "75": 0.4614379107952118, + "76": 0.488725483417511, + "77": 0.488725483417511, + "78": 0.4485294222831726, + "79": 0.4333333373069763, + "80": 0.4413398802280426, + "81": 0.3611111044883728, + "82": 0.35049018263816833, + "83": 0.36781045794487, + "84": 0.3942810595035553, + "85": 0.4233660101890564, + "86": 0.48398691415786743, + "87": 0.5047385692596436, + "88": 0.4915032684803009, + "89": 0.47156861424446106, + "90": 0.4459150433540344, + "91": 0.32810458540916443, + "92": 0.32777777314186096, + "93": 0.37516340613365173, + "94": 0.41241830587387085, + "95": 0.4310457408428192, "96": 0.0, "97": 0.0, "98": 0.0, - "99": 0.5098039507865906, + "99": 0.5104575157165527, "100": 0.4519607722759247, - "101": 0.4454248249530792, - "102": 0.4629085063934326, - "103": 0.4870915114879608, - "104": 0.4740196168422699, - "105": 0.48627451062202454, - "106": 0.5398693084716797, + "101": 0.4480392038822174, + "102": 0.46486929059028625, + "103": 0.4882352948188782, + "104": 0.47467321157455444, + "105": 0.4870915114879608, + "106": 0.5408496856689453, "107": 0.0, - "108": 0.4812091588973999, - "109": 0.4573529362678528, - "110": 0.4117647111415863, - "111": 0.43562090396881104, - "112": 0.46666666865348816, - "113": 0.49934640526771545, - "114": 0.4905228614807129, - "115": 0.5297385454177856, + "108": 0.48235294222831726, + "109": 0.45800653100013733, + "110": 0.41241830587387085, + "111": 0.43790850043296814, + "112": 0.46846404671669006, + "113": 0.5003268122673035, + "114": 0.49117645621299744, + "115": 0.5310457348823547, "116": 0.0, "117": 0.0, - "118": 0.54313725233078, - "119": 0.5058823823928833, - "120": 0.47075164318084717, - "121": 0.44738560914993286, - "122": 0.46568626165390015, - "123": 0.489705890417099, - "124": 0.48725488781929016, - "125": 0.5302287340164185, - "126": 0.584967315196991, - "127": 0.5485293865203857, - "128": 0.5150327086448669, - "129": 0.5299019813537598, - "130": 0.5292483568191528, - "131": 0.4709150195121765, - "132": 0.49215686321258545, - "133": 0.5187908411026001, - "134": 0.49787580966949463, - "135": 0.5063725709915161, - "136": 0.5357843041419983, - "137": 0.5016340017318726, - "138": 0.5145424604415894, - "139": 0.568790853023529, - "140": 0.5598039031028748, - "141": 0.4588235318660736, - "142": 0.4866012930870056, - "143": 0.5212418437004089, - "144": 0.4995098114013672, - "145": 0.499019593000412, - "146": 0.5470588207244873, - "147": 0.5633987188339233, - "148": 0.5797385573387146, + "118": 0.5444444417953491, + "119": 0.5066993236541748, + "120": 0.4718954265117645, + "121": 0.4490196108818054, + "122": 0.4673202633857727, + "123": 0.4910130798816681, + "124": 0.48807188868522644, + "125": 0.5312091708183289, + "126": 0.5859476923942566, + "127": 0.5495098233222961, + "128": 0.5160130858421326, + "129": 0.5310457348823547, + "130": 0.5305555462837219, + "131": 0.4725490212440491, + "132": 0.49395424127578735, + "133": 0.520588219165802, + "134": 0.499019593000412, + "135": 0.5070261359214783, + "136": 0.5364379286766052, + "137": 0.5022875666618347, + "138": 0.5155228972434998, + "139": 0.5704248547554016, + "140": 0.5609477162361145, + "141": 0.4611110985279083, + "142": 0.4882352948188782, + "143": 0.522549033164978, + "144": 0.5006536245346069, + "145": 0.49967318773269653, + "146": 0.5478758215904236, + "147": 0.5640522837638855, + "148": 0.5807189345359802, "149": 0.0, - "150": 0.5604575276374817, - "151": 0.45359477400779724, - "152": 0.4779411852359772, - "153": 0.5165032744407654, - "154": 0.5153594613075256, - "155": 0.5140522718429565, - "156": 0.556209146976471, - "157": 0.5821895599365234, - "158": 0.5537581443786621, - "159": 0.5326797366142273, - "160": 0.5454248189926147, - "161": 0.4627451002597809, - "162": 0.4883987009525299, + "150": 0.561274528503418, + "151": 0.4552287459373474, + "152": 0.47924837470054626, + "153": 0.5178104639053345, + "154": 0.5166666507720947, + "155": 0.5147058963775635, + "156": 0.5566993355751038, + "157": 0.5828431248664856, + "158": 0.5550653338432312, + "159": 0.5341503024101257, + "160": 0.546241819858551, + "161": 0.4637254774570465, + "162": 0.489705890417099, "163": 0.0, - "164": 0.5616012811660767, - "165": 0.5424836874008179, - "166": 0.5397058725357056, - "167": 0.5362744927406311, - "168": 0.48186275362968445, - "169": 0.4547385573387146, - "170": 0.49754902720451355, - "171": 0.40294116735458374, - "172": 0.43839868903160095, - "173": 0.5094771385192871, + "164": 0.5625817179679871, + "165": 0.5433006286621094, + "166": 0.5406862497329712, + "167": 0.5377451181411743, + "168": 0.48366013169288635, + "169": 0.4562091529369354, + "170": 0.49836599826812744, + "171": 0.4045751690864563, + "172": 0.43970587849617004, + "173": 0.5104575157165527, "174": 0.0, - "175": 0.5011438131332397, - "176": 0.496895432472229, - "177": 0.5099673271179199, - "178": 0.46405228972435, - "179": 0.4086601436138153, - "180": 0.4187908470630646, - "181": 0.3336601257324219, - "182": 0.38153594732284546, - "183": 0.460947722196579, - "184": 0.47516340017318726, - "185": 0.44918301701545715, - "186": 0.4619280993938446, - "187": 0.4915032684803009, - "188": 0.4699346423149109, - "189": 0.4138889014720917, - "190": 0.3870915174484253, - "191": 0.31699347496032715, - "192": 0.3686274588108063, - "193": 0.4503268003463745, - "194": 0.46666666865348816, - "195": 0.4516339898109436, - "196": 0.4660130739212036, - "197": 0.4843137264251709, - "198": 0.47287580370903015, - "199": 0.43431371450424194, - "200": 0.39183005690574646, - "201": 0.4718954265117645, - "202": 0.5065359473228455, + "175": 0.5029411911964417, + "176": 0.49836599826812744, + "177": 0.5116013288497925, + "178": 0.46568626165390015, + "179": 0.4099673330783844, + "180": 0.4192810356616974, + "181": 0.3356209099292755, + "182": 0.3826797306537628, + "183": 0.46176469326019287, + "184": 0.4771241843700409, + "185": 0.45130717754364014, + "186": 0.4632352888584137, + "187": 0.4923202693462372, + "188": 0.47075164318084717, + "189": 0.41470587253570557, + "190": 0.38741829991340637, + "191": 0.3189542591571808, + "192": 0.36993464827537537, + "193": 0.45098039507865906, + "194": 0.4676470458507538, + "195": 0.4532679617404938, + "196": 0.4673202633857727, + "197": 0.484640508890152, + "198": 0.4732026159763336, + "199": 0.4349673092365265, + "200": 0.3921568691730499, + "201": 0.47385621070861816, + "202": 0.5080065131187439, "203": 0.0, "204": 0.0, - "205": 0.5223855972290039, - "206": 0.5196078419685364, - "207": 0.4732026159763336, - "208": 0.48316994309425354, + "205": 0.5233660340309143, + "206": 0.520588219165802, + "207": 0.47418299317359924, + "208": 0.48415032029151917, "209": 0.0, - "210": 0.5460784435272217, - "211": 0.48627451062202454, - "212": 0.5034313797950745, - "213": 0.5251634120941162, - "214": 0.5240195989608765, - "215": 0.552450954914093, - "216": 0.5740196108818054, - "217": 0.533823549747467, - "218": 0.5411764979362488, + "210": 0.5464052557945251, + "211": 0.48856207728385925, + "212": 0.5049019455909729, + "213": 0.5259804129600525, + "214": 0.5249999761581421, + "215": 0.5537581443786621, + "216": 0.5753268003463745, + "217": 0.5352941155433655, + "218": 0.5421568751335144, "219": 0.0, - "220": 0.49983659386634827, - "221": 0.4761437773704529, - "222": 0.4838235378265381, - "223": 0.5117647051811218, - "224": 0.5305555462837219, - "225": 0.5611110925674438, - "226": 0.5774509906768799, - "227": 0.5392156839370728, - "228": 0.5398693084716797, - "229": 0.5413398742675781, - "230": 0.4735293984413147, - "231": 0.4650326669216156, - "232": 0.46486929059028625, - "233": 0.4982026219367981, - "234": 0.5375816822052002, - "235": 0.5426470637321472, - "236": 0.5222222208976746, - "237": 0.503104567527771, - "238": 0.5334967374801636, + "220": 0.5006536245346069, + "221": 0.4771241843700409, + "222": 0.484640508890152, + "223": 0.5125817060470581, + "224": 0.5313725471496582, + "225": 0.5619280934333801, + "226": 0.5782679915428162, + "227": 0.540032684803009, + "228": 0.5408496856689453, + "229": 0.5428104400634766, + "230": 0.47516340017318726, + "231": 0.4660130739212036, + "232": 0.4660130739212036, + "233": 0.4991829991340637, + "234": 0.538071870803833, + "235": 0.5429738759994507, + "236": 0.5227124094963074, + "237": 0.5037581920623779, + "238": 0.534967303276062, "239": 0.0, - "240": 0.5016340017318726, - "241": 0.44411763548851013, - "242": 0.4429738521575928, - "243": 0.47058823704719543, - "244": 0.5279411673545837, - "245": 0.5354574918746948, - "246": 0.5052287578582764, - "247": 0.529411792755127, + "240": 0.5034313797950745, + "241": 0.44640523195266724, + "242": 0.44460785388946533, + "243": 0.4717320203781128, + "244": 0.5290849804878235, + "245": 0.5361111164093018, + "246": 0.5058823823928833, + "247": 0.5302287340164185, "248": 0.0, "249": 0.0, - "250": 0.5140522718429565, - "251": 0.4156862795352936, - "252": 0.44199347496032715, - "253": 0.48104575276374817, - "254": 0.5199346542358398, - "255": 0.5343137383460999, - "256": 0.5119280815124512, - "257": 0.525653600692749, - "258": 0.5372549295425415, - "259": 0.5114378929138184, - "260": 0.5042483806610107, - "261": 0.43790850043296814, - "262": 0.4678104519844055, + "250": 0.5151960849761963, + "251": 0.41699346899986267, + "252": 0.4431372582912445, + "253": 0.482516348361969, + "254": 0.5215686559677124, + "255": 0.5351307392120361, + "256": 0.5127450823783875, + "257": 0.5264706015586853, + "258": 0.5385621190071106, + "259": 0.5132352709770203, + "260": 0.505065381526947, + "261": 0.4382352828979492, + "262": 0.46879085898399353, "263": 0.0, - "264": 0.5388888716697693, - "265": 0.5428104400634766, - "266": 0.5199346542358398, - "267": 0.49624183773994446, - "268": 0.4718954265117645, - "269": 0.45179739594459534, - "270": 0.4696078300476074, - "271": 0.4542483687400818, - "272": 0.46535947918891907, - "273": 0.5026143789291382, - "274": 0.5326797366142273, - "275": 0.516339898109436, - "276": 0.48774510622024536, - "277": 0.47957515716552734, - "278": 0.47107842564582825, - "279": 0.44411763548851013, - "280": 0.4328431487083435, - "281": 0.44150325655937195, - "282": 0.4583333432674408, - "283": 0.4905228614807129, - "284": 0.5068627595901489, - "285": 0.4851307272911072, - "286": 0.46437907218933105, - "287": 0.4709150195121765, - "288": 0.4709150195121765, - "289": 0.4503268003463745, - "290": 0.42990195751190186, - "291": 0.44117647409439087, - "292": 0.4663398563861847, - "293": 0.4973856210708618, - "294": 0.5042483806610107, - "295": 0.48986926674842834, - "296": 0.4797385632991791, - "297": 0.47908496856689453, - "298": 0.46405228972435, - "299": 0.4493463933467865, - "300": 0.4336601197719574, - "301": 0.484640508890152, - "302": 0.4910130798816681, + "264": 0.540032684803009, + "265": 0.5439542531967163, + "266": 0.5214052200317383, + "267": 0.49787580966949463, + "268": 0.47369280457496643, + "269": 0.4534313678741455, + "270": 0.47026142477989197, + "271": 0.45490196347236633, + "272": 0.4665032625198364, + "273": 0.5037581920623779, + "274": 0.5341503024101257, + "275": 0.5184640288352966, + "276": 0.489705890417099, + "277": 0.481535941362381, + "278": 0.47287580370903015, + "279": 0.44509804248809814, + "280": 0.43300652503967285, + "281": 0.4431372582912445, + "282": 0.45980390906333923, + "283": 0.4919934570789337, + "284": 0.5089869499206543, + "285": 0.48725488781929016, + "286": 0.4658496677875519, + "287": 0.47238561511039734, + "288": 0.47238561511039734, + "289": 0.4511438012123108, + "290": 0.4300653636455536, + "291": 0.44379085302352905, + "292": 0.46830064058303833, + "293": 0.499019593000412, + "294": 0.506209135055542, + "295": 0.4915032684803009, + "296": 0.4807189404964447, + "297": 0.48006534576416016, + "298": 0.46535947918891907, + "299": 0.4503268003463745, + "300": 0.43398693203926086, + "301": 0.4852941036224365, + "302": 0.491830050945282, "303": 0.0, "304": 0.0, - "305": 0.5442810654640198, - "306": 0.49346405267715454, - "307": 0.4418300688266754, - "308": 0.4931372404098511, + "305": 0.545098066329956, + "306": 0.49444442987442017, + "307": 0.44248366355895996, + "308": 0.4937908351421356, "309": 0.0, - "310": 0.59624183177948, - "311": 0.5127450823783875, - "312": 0.5176470875740051, + "310": 0.5965686440467834, + "311": 0.5140522718429565, + "312": 0.5186274647712708, "313": 0.0, "314": 0.0, - "315": 0.5661764740943909, - "316": 0.5388888716697693, - "317": 0.49346405267715454, - "318": 0.5235294103622437, - "319": 0.5457516312599182, - "320": 0.4647058844566345, - "321": 0.4833333194255829, - "322": 0.4937908351421356, - "323": 0.5385621190071106, - "324": 0.5774509906768799, - "325": 0.5735294222831726, - "326": 0.5370914936065674, - "327": 0.5071895718574524, - "328": 0.5156862735748291, - "329": 0.4802287518978119, - "330": 0.3952614367008209, - "331": 0.4362744987010956, - "332": 0.4287581741809845, - "333": 0.46928104758262634, - "334": 0.5511437654495239, - "335": 0.560620903968811, - "336": 0.5047385692596436, - "337": 0.4802287518978119, + "315": 0.5669934749603271, + "316": 0.540032684803009, + "317": 0.49444442987442017, + "318": 0.5243464112281799, + "319": 0.5464052557945251, + "320": 0.46519607305526733, + "321": 0.48398691415786743, + "322": 0.49444442987442017, + "323": 0.5392156839370728, + "324": 0.5779411792755127, + "325": 0.5741829872131348, + "326": 0.5377451181411743, + "327": 0.5078431367874146, + "328": 0.5169934630393982, + "329": 0.4816993474960327, + "330": 0.39656862616539, + "331": 0.4375816881656647, + "332": 0.4302287697792053, + "333": 0.47026142477989197, + "334": 0.5517973899841309, + "335": 0.561274528503418, + "336": 0.5053921341896057, + "337": 0.48137253522872925, "338": 0.0, - "339": 0.48807188868522644, - "340": 0.4498366117477417, - "341": 0.4326797425746918, - "342": 0.42189541459083557, - "343": 0.4318627417087555, - "344": 0.5153594613075256, - "345": 0.5563725233078003, - "346": 0.4986928105354309, - "347": 0.48758170008659363, - "348": 0.5192810297012329, - "349": 0.5114378929138184, - "350": 0.4905228614807129, - "351": 0.42941176891326904, - "352": 0.44869279861450195, - "353": 0.46928104758262634, - "354": 0.5037581920623779, - "355": 0.5254902243614197, - "356": 0.4977124035358429, - "357": 0.49526143074035645, - "358": 0.49787580966949463, - "359": 0.4696078300476074, - "360": 0.4552287459373474, - "361": 0.4483660161495209, - "362": 0.4645424783229828, - "363": 0.48594769835472107, - "364": 0.49166667461395264, - "365": 0.506209135055542, - "366": 0.5218954086303711, - "367": 0.5008170008659363, - "368": 0.4665032625198364, - "369": 0.4439542591571808, - "370": 0.4264705777168274, - "371": 0.5192810297012329, - "372": 0.49395424127578735, - "373": 0.4645424783229828, - "374": 0.48055556416511536, - "375": 0.5217320322990417, - "376": 0.516339898109436, - "377": 0.475326806306839, - "378": 0.48006534576416016, - "379": 0.4964052140712738, - "380": 0.49656862020492554, - "381": 0.564379096031189, - "382": 0.530718982219696, - "383": 0.48104575276374817, - "384": 0.49166667461395264, - "385": 0.5263071656227112, - "386": 0.496895432472229, - "387": 0.45359477400779724, - "388": 0.47565358877182007, - "389": 0.519444465637207, - "390": 0.5570261478424072, - "391": 0.5699346661567688, - "392": 0.5490196347236633, - "393": 0.5120915174484253, - "394": 0.5101307034492493, - "395": 0.5267974138259888, - "396": 0.5055555701255798, - "397": 0.46666666865348816, - "398": 0.46405228972435, - "399": 0.5022875666618347, - "400": 0.5532679557800293, - "401": 0.5235294103622437, - "402": 0.49035948514938354, - "403": 0.4870915114879608, - "404": 0.5289215445518494, - "405": 0.5218954086303711, - "406": 0.4604575037956238, - "407": 0.4031045734882355, - "408": 0.46176469326019287, + "339": 0.4900326728820801, + "340": 0.45147058367729187, + "341": 0.4349673092365265, + "342": 0.42352941632270813, + "343": 0.43300652503967285, + "344": 0.5169934630393982, + "345": 0.5576797127723694, + "346": 0.49967318773269653, + "347": 0.488725483417511, + "348": 0.520588219165802, + "349": 0.5129085183143616, + "350": 0.4915032684803009, + "351": 0.43039214611053467, + "352": 0.44967320561408997, + "353": 0.47058823704719543, + "354": 0.5055555701255798, + "355": 0.5269607901573181, + "356": 0.499019593000412, + "357": 0.4964052140712738, + "358": 0.49885621666908264, + "359": 0.4709150195121765, + "360": 0.45588234066963196, + "361": 0.44869279861450195, + "362": 0.4655228853225708, + "363": 0.4870915114879608, + "364": 0.4931372404098511, + "365": 0.5080065131187439, + "366": 0.523692786693573, + "367": 0.5027777552604675, + "368": 0.46830064058303833, + "369": 0.4452614486217499, + "370": 0.42679738998413086, + "371": 0.5202614665031433, + "372": 0.49526143074035645, + "373": 0.4660130739212036, + "374": 0.48235294222831726, + "375": 0.5238562226295471, + "376": 0.5184640288352966, + "377": 0.4776143729686737, + "378": 0.48235294222831726, + "379": 0.4977124035358429, + "380": 0.49673202633857727, + "381": 0.5660130977630615, + "382": 0.532516360282898, + "383": 0.48316994309425354, + "384": 0.4936274588108063, + "385": 0.5279411673545837, + "386": 0.4985294044017792, + "387": 0.4557189643383026, + "388": 0.4776143729686737, + "389": 0.5202614665031433, + "390": 0.5573529601097107, + "391": 0.5718954205513, + "392": 0.5509803891181946, + "393": 0.5140522718429565, + "394": 0.5117647051811218, + "395": 0.5281046032905579, + "396": 0.5068627595901489, + "397": 0.4686274528503418, + "398": 0.46568626165390015, + "399": 0.5026143789291382, + "400": 0.5535947680473328, + "401": 0.5238562226295471, + "402": 0.49117645621299744, + "403": 0.48807188868522644, + "404": 0.529411792755127, + "405": 0.522549033164978, + "406": 0.4619280993938446, + "407": 0.4045751690864563, + "408": 0.4627451002597809, "409": 0.0, - "410": 0.6008169651031494, - "411": 0.49673202633857727, - "412": 0.5060457587242126, - "413": 0.5406862497329712, - "414": 0.5663398504257202, - "415": 0.5517973899841309, - "416": 0.5200980305671692, - "417": 0.5055555701255798, - "418": 0.5271241664886475, - "419": 0.5290849804878235, - "420": 0.48986926674842834, - "421": 0.4467320144176483, - "422": 0.4606209099292755, - "423": 0.5130718946456909, - "424": 0.561274528503418, - "425": 0.546895444393158, - "426": 0.5142157077789307, - "427": 0.525653600692749, + "410": 0.6011437773704529, + "411": 0.4973856210708618, + "412": 0.5068627595901489, + "413": 0.5418300628662109, + "414": 0.5669934749603271, + "415": 0.552450954914093, + "416": 0.5217320322990417, + "417": 0.5070261359214783, + "418": 0.5284313559532166, + "419": 0.5300653576850891, + "420": 0.4900326728820801, + "421": 0.44738560914993286, + "422": 0.4616013169288635, + "423": 0.5140522718429565, + "424": 0.5619280934333801, + "425": 0.5477124452590942, + "426": 0.5150327086448669, + "427": 0.0, "428": 0.0, - "429": 0.48725488781929016, - "430": 0.4217320382595062, - "431": 0.40980392694473267, - "432": 0.3973856270313263, - "433": 0.4256536066532135, - "434": 0.5058823823928833, - "435": 0.5200980305671692, - "436": 0.4861111044883728, + "429": 0.48856207728385925, + "430": 0.42271241545677185, + "431": 0.41111111640930176, + "432": 0.39901959896087646, + "433": 0.42728757858276367, + "434": 0.5075163245201111, + "435": 0.5214052200317383, + "436": 0.4866012930870056, "437": 0.0, "438": 0.0, - "439": 0.49983659386634827, - "440": 0.47156861424446106, - "441": 0.4454248249530792, - "442": 0.43529412150382996, - "443": 0.4385620951652527, - "444": 0.48235294222831726, - "445": 0.5022875666618347, - "446": 0.49803921580314636, - "447": 0.5122548937797546, - "448": 0.5098039507865906, - "449": 0.5053921341896057, - "450": 0.5130718946456909, + "439": 0.5016340017318726, + "440": 0.4730392098426819, + "441": 0.44738560914993286, + "442": 0.4369280934333801, + "443": 0.4400326907634735, + "444": 0.48415032029151917, + "445": 0.5037581920623779, + "446": 0.499019593000412, + "447": 0.5137255191802979, + "448": 0.5114378929138184, + "449": 0.5066993236541748, + "450": 0.5142157077789307, "451": 0.0, "452": 0.0, - "453": 0.5124183297157288, - "454": 0.5014705657958984, - "455": 0.4869281053543091, - "456": 0.5214052200317383, - "457": 0.5382353067398071, - "458": 0.5040849447250366, - "459": 0.48186275362968445, - "460": 0.4803921580314636, + "453": 0.5133987069129944, + "454": 0.5024510025978088, + "455": 0.4883987009525299, + "456": 0.5228758454322815, + "457": 0.5395424962043762, + "458": 0.5052287578582764, + "459": 0.4830065369606018, + "460": 0.48137253522872925, "461": 0.0, - "462": 0.49803921580314636, - "463": 0.5155228972434998, - "464": 0.48758170008659363, - "465": 0.48055556416511536, - "466": 0.5426470637321472, - "467": 0.5357843041419983, - "468": 0.489705890417099, - "469": 0.47696077823638916, - "470": 0.4519607722759247, - "471": 0.5218954086303711, - "472": 0.5099673271179199, - "473": 0.48398691415786743, - "474": 0.47336602210998535, - "475": 0.5125817060470581, - "476": 0.5375816822052002, - "477": 0.5045751929283142, - "478": 0.4955882430076599, - "479": 0.5014705657958984, - "480": 0.4915032684803009, - "481": 0.5653594732284546, - "482": 0.5454248189926147, - "483": 0.5075163245201111, - "484": 0.511274516582489, + "462": 0.4991829991340637, + "463": 0.5168300867080688, + "464": 0.48905229568481445, + "465": 0.4821895360946655, + "466": 0.543790876865387, + "467": 0.5370914936065674, + "468": 0.4915032684803009, + "469": 0.47826796770095825, + "470": 0.45228758454322815, + "471": 0.5232025980949402, + "472": 0.511274516582489, + "473": 0.4852941036224365, + "474": 0.47516340017318726, + "475": 0.5142157077789307, + "476": 0.5385621190071106, + "477": 0.5060457587242126, + "478": 0.49754902720451355, + "479": 0.5026143789291382, + "480": 0.491830050945282, + "481": 0.5666666626930237, + "482": 0.5467320084571838, + "483": 0.5084967613220215, + "484": 0.5122548937797546, "485": 0.0, - "486": 0.5204248428344727, - "487": 0.46519607305526733, - "488": 0.4727124273777008, - "489": 0.5120915174484253, - "490": 0.5568627715110779, - "491": 0.5732026100158691, - "492": 0.5620915293693542, - "493": 0.5424836874008179, - "494": 0.5447712540626526, + "486": 0.5215686559677124, + "487": 0.46666666865348816, + "488": 0.474346399307251, + "489": 0.5129085183143616, + "490": 0.55751633644104, + "491": 0.5745097994804382, + "492": 0.5633987188339233, + "493": 0.54313725233078, + "494": 0.545098066329956, "495": 0.0, - "496": 0.5179738402366638, - "497": 0.45359477400779724, - "498": 0.44640523195266724, - "499": 0.5039215683937073, - "500": 0.5797385573387146, - "501": 0.5866013169288635, - "502": 0.5570261478424072, - "503": 0.505065381526947, - "504": 0.4802287518978119, - "505": 0.481535941362381, - "506": 0.4413398802280426, - "507": 0.38415032625198364, - "508": 0.45898693799972534, - "509": 0.5844771265983582, - "510": 0.6161764860153198, - "511": 0.5006536245346069, - "512": 0.5238562226295471, - "513": 0.5540849566459656, - "514": 0.5416666865348816, - "515": 0.5271241664886475, - "516": 0.531862735748291, - "517": 0.5230392217636108, - "518": 0.5385621190071106, - "519": 0.5650326609611511, - "520": 0.5517973899841309, - "521": 0.44869279861450195, - "522": 0.4771241843700409, - "523": 0.5259804129600525, + "496": 0.5192810297012329, + "497": 0.4552287459373474, + "498": 0.4480392038822174, + "499": 0.5045751929283142, + "500": 0.5803921818733215, + "501": 0.5875816941261292, + "502": 0.5578431487083435, + "503": 0.5060457587242126, + "504": 0.48088234663009644, + "505": 0.4816993474960327, + "506": 0.4421568512916565, + "507": 0.385294109582901, + "508": 0.4601307213306427, + "509": 0.5857843160629272, + "510": 0.616830050945282, + "511": 0.5013071894645691, + "512": 0.5243464112281799, + "513": 0.5550653338432312, + "514": 0.5426470637321472, + "515": 0.5276143550872803, + "516": 0.5331699252128601, + "517": 0.5245097875595093, + "518": 0.5398693084716797, + "519": 0.5666666626930237, + "520": 0.5526143908500671, + "521": 0.4493463933467865, + "522": 0.4779411852359772, + "523": 0.5272876024246216, "524": 0.0, - "525": 0.5080065131187439, - "526": 0.5186274647712708, - "527": 0.5529412031173706, + "525": 0.5089869499206543, + "526": 0.5196078419685364, + "527": 0.554411768913269, "528": 0.0, "529": 0.0, - "530": 0.48104575276374817, - "531": 0.44607841968536377, - "532": 0.45457515120506287, - "533": 0.4749999940395355, - "534": 0.4933006465435028, - "535": 0.4807189404964447, - "536": 0.47107842564582825, + "530": 0.4821895360946655, + "531": 0.44738560914993286, + "532": 0.4557189643383026, + "533": 0.47647058963775635, + "534": 0.4950980246067047, + "535": 0.48235294222831726, + "536": 0.4718954265117645, "537": 0.0, "538": 0.0, "539": 0.0, - "540": 0.48627451062202454, - "541": 0.49346405267715454, - "542": 0.482516348361969, - "543": 0.4838235378265381, - "544": 0.49444442987442017, - "545": 0.4833333194255829, - "546": 0.5076797604560852, - "547": 0.5447712540626526, + "540": 0.48774510622024536, + "541": 0.4954248368740082, + "542": 0.48398691415786743, + "543": 0.48545750975608826, + "544": 0.4964052140712738, + "545": 0.48496732115745544, + "546": 0.5086601376533508, + "547": 0.5457516312599182, "548": 0.0, - "549": 0.5253267884254456, - "550": 0.5106208920478821, + "549": 0.5267974138259888, + "550": 0.5120915174484253, "551": 0.0, - "552": 0.5116013288497925, - "553": 0.5248365998268127, - "554": 0.5197712182998657, - "555": 0.491830050945282, - "556": 0.5374183058738708, - "557": 0.5565359592437744, - "558": 0.5165032744407654, - "559": 0.5057189464569092, - "560": 0.5166666507720947, - "561": 0.49673202633857727, - "562": 0.4964052140712738, - "563": 0.5161764621734619, - "564": 0.5032680034637451, - "565": 0.48758170008659363, - "566": 0.5207516551017761, - "567": 0.5088235139846802, - "568": 0.49656862020492554, - "569": 0.5078431367874146, - "570": 0.495915025472641, - "571": 0.5160130858421326, - "572": 0.5098039507865906, - "573": 0.48562091588974, - "574": 0.46519607305526733, - "575": 0.4919934570789337, - "576": 0.5202614665031433, - "577": 0.5114378929138184, - "578": 0.5259804129600525, - "579": 0.5116013288497925, - "580": 0.4470588266849518, - "581": 0.5624182820320129, - "582": 0.5444444417953491, - "583": 0.5057189464569092, - "584": 0.5013071894645691, - "585": 0.5326797366142273, - "586": 0.5292483568191528, - "587": 0.5009803771972656, - "588": 0.5114378929138184, - "589": 0.5078431367874146, - "590": 0.4694444537162781, - "591": 0.5718954205513, - "592": 0.5529412031173706, - "593": 0.5320261716842651, - "594": 0.5464052557945251, - "595": 0.5617647171020508, - "596": 0.5300653576850891, - "597": 0.4748365879058838, + "552": 0.5132352709770203, + "553": 0.5261437892913818, + "554": 0.5207516551017761, + "555": 0.4931372404098511, + "556": 0.5385621190071106, + "557": 0.5573529601097107, + "558": 0.517483651638031, + "559": 0.5070261359214783, + "560": 0.5179738402366638, + "561": 0.4973856210708618, + "562": 0.4973856210708618, + "563": 0.5171568393707275, + "564": 0.5042483806610107, + "565": 0.48905229568481445, + "566": 0.5217320322990417, + "567": 0.5096405148506165, + "568": 0.49803921580314636, + "569": 0.5091503262519836, + "570": 0.4964052140712738, + "571": 0.5169934630393982, + "572": 0.5106208920478821, + "573": 0.48627451062202454, + "574": 0.4665032625198364, + "575": 0.4936274588108063, + "576": 0.5210784077644348, + "577": 0.5124183297157288, + "578": 0.5274509787559509, + "579": 0.5125817060470581, + "580": 0.44771242141723633, + "581": 0.5633987188339233, + "582": 0.5455882549285889, + "583": 0.5065359473228455, + "584": 0.5021241903305054, + "585": 0.5339869260787964, + "586": 0.5302287340164185, + "587": 0.501960813999176, + "588": 0.5125817060470581, + "589": 0.5088235139846802, + "590": 0.47026142477989197, + "591": 0.572549045085907, + "592": 0.5542483925819397, + "593": 0.5333333611488342, + "594": 0.5470588207244873, + "595": 0.5627450942993164, + "596": 0.5310457348823547, + "597": 0.4761437773704529, "598": 0.0, - "599": 0.5042483806610107, - "600": 0.5120915174484253, - "601": 0.5483660101890564, - "602": 0.5279411673545837, - "603": 0.4506535828113556, - "604": 0.3978758156299591, - "605": 0.44379085302352905, - "606": 0.4660130739212036, - "607": 0.4220588207244873, - "608": 0.4431372582912445, - "609": 0.51437908411026, - "610": 0.545098066329956, - "611": 0.4732026159763336, - "612": 0.49248364567756653, - "613": 0.49705880880355835, - "614": 0.4629085063934326, - "615": 0.4821895360946655, - "616": 0.5297385454177856, - "617": 0.5302287340164185, - "618": 0.538071870803833, - "619": 0.5387254953384399, - "620": 0.5111111402511597, - "621": 0.43529412150382996, - "622": 0.4732026159763336, - "623": 0.5202614665031433, + "599": 0.5055555701255798, + "600": 0.5127450823783875, + "601": 0.5490196347236633, + "602": 0.5284313559532166, + "603": 0.4516339898109436, + "604": 0.3986928164958954, + "605": 0.44428104162216187, + "606": 0.46666666865348816, + "607": 0.4225490093231201, + "608": 0.44379085302352905, + "609": 0.5153594613075256, + "610": 0.5457516312599182, + "611": 0.47385621070861816, + "612": 0.49281045794487, + "613": 0.49803921580314636, + "614": 0.46437907218933105, + "615": 0.4830065369606018, + "616": 0.530718982219696, + "617": 0.5312091708183289, + "618": 0.5387254953384399, + "619": 0.5403594970703125, + "620": 0.5124183297157288, + "621": 0.43660131096839905, + "622": 0.47385621070861816, + "623": 0.5214052200317383, "624": 0.0, - "625": 0.4986928105354309, - "626": 0.5276143550872803, - "627": 0.5547385811805725, - "628": 0.5614379048347473, - "629": 0.5186274647712708, - "630": 0.47549018263816833, - "631": 0.4650326669216156, - "632": 0.4906862676143646, - "633": 0.5217320322990417, - "634": 0.5187908411026001, - "635": 0.5145424604415894, - "636": 0.5080065131187439, - "637": 0.5027777552604675, - "638": 0.51437908411026, - "639": 0.49983659386634827, - "640": 0.4627451002597809, - "641": 0.5147058963775635, - "642": 0.4964052140712738, - "643": 0.4923202693462372, - "644": 0.5124183297157288, - "645": 0.5267974138259888, - "646": 0.5434640645980835, - "647": 0.535620927810669, - "648": 0.5179738402366638, - "649": 0.4954248368740082, - "650": 0.4434640407562256, - "651": 0.5150327086448669, - "652": 0.48267972469329834, - "653": 0.484640508890152, - "654": 0.5248365998268127, - "655": 0.5220588445663452, - "656": 0.5316993594169617, - "657": 0.5405228734016418, - "658": 0.5204248428344727, - "659": 0.48725488781929016, - "660": 0.4434640407562256, - "661": 0.4931372404098511, - "662": 0.4763071835041046, - "663": 0.48398691415786743, - "664": 0.5047385692596436, - "665": 0.4941176474094391, - "666": 0.4828431308269501, - "667": 0.47957515716552734, - "668": 0.49967318773269653, - "669": 0.4986928105354309, - "670": 0.46029412746429443, - "671": 0.5421568751335144, - "672": 0.5197712182998657, - "673": 0.46928104758262634, - "674": 0.4382352828979492, - "675": 0.46225491166114807, - "676": 0.481535941362381, - "677": 0.4901960790157318, - "678": 0.5292483568191528, - "679": 0.5120915174484253, - "680": 0.4506535828113556, - "681": 0.6032679677009583, - "682": 0.5681372284889221, - "683": 0.49395424127578735, - "684": 0.4627451002597809, - "685": 0.49395424127578735, - "686": 0.4923202693462372, - "687": 0.48643791675567627, - "688": 0.5263071656227112, + "625": 0.49983659386634827, + "626": 0.5289215445518494, + "627": 0.5560457706451416, + "628": 0.5624182820320129, + "629": 0.5202614665031433, + "630": 0.47696077823638916, + "631": 0.46699345111846924, + "632": 0.49215686321258545, + "633": 0.5233660340309143, + "634": 0.520588219165802, + "635": 0.5160130858421326, + "636": 0.5091503262519836, + "637": 0.5034313797950745, + "638": 0.5151960849761963, + "639": 0.5009803771972656, + "640": 0.4637254774570465, + "641": 0.5173202753067017, + "642": 0.4985294044017792, + "643": 0.49444442987442017, + "644": 0.5145424604415894, + "645": 0.5281046032905579, + "646": 0.5442810654640198, + "647": 0.5362744927406311, + "648": 0.5186274647712708, + "649": 0.4964052140712738, + "650": 0.44428104162216187, + "651": 0.5173202753067017, + "652": 0.4848039150238037, + "653": 0.48643791675567627, + "654": 0.5261437892913818, + "655": 0.5232025980949402, + "656": 0.5331699252128601, + "657": 0.5418300628662109, + "658": 0.0, + "659": 0.4883987009525299, + "660": 0.4444444477558136, + "661": 0.4937908351421356, + "662": 0.477450966835022, + "663": 0.4848039150238037, + "664": 0.5052287578582764, + "665": 0.4955882430076599, + "666": 0.4843137264251709, + "667": 0.48055556416511536, + "668": 0.5004901885986328, + "669": 0.4995098114013672, + "670": 0.460947722196579, + "671": 0.5434640645980835, + "672": 0.520588219165802, + "673": 0.4696078300476074, + "674": 0.43888887763023376, + "675": 0.46339869499206543, + "676": 0.482516348361969, + "677": 0.4910130798816681, + "678": 0.5302287340164185, + "679": 0.5130718946456909, + "680": 0.45147058367729187, + "681": 0.6049019694328308, + "682": 0.5692810416221619, + "683": 0.4946078360080719, + "684": 0.46339869499206543, + "685": 0.49477124214172363, + "686": 0.49346405267715454, + "687": 0.4879084825515747, + "688": 0.5276143550872803, "689": 0.0, - "690": 0.47107842564582825, - "691": 0.615686297416687, - "692": 0.5764706134796143, - "693": 0.5196078419685364, - "694": 0.5160130858421326, - "695": 0.5359477400779724, - "696": 0.4954248368740082, - "697": 0.4565359354019165, + "690": 0.4718954265117645, + "691": 0.6163398623466492, + "692": 0.5774509906768799, + "693": 0.520588219165802, + "694": 0.5166666507720947, + "695": 0.536928117275238, + "696": 0.49673202633857727, + "697": 0.45849671959877014, "698": 0.0, "699": 0.0, "700": 0.0, - "701": 0.4261437952518463, - "702": 0.4297385513782501, - "703": 0.40375816822052, - "704": 0.37875816226005554, - "705": 0.42532679438591003, - "706": 0.4838235378265381, - "707": 0.4833333194255829, - "708": 0.42679738998413086, - "709": 0.4122548997402191, - "710": 0.4870915114879608, - "711": 0.41143789887428284, - "712": 0.42091503739356995, - "713": 0.42581698298454285, - "714": 0.41862744092941284, - "715": 0.4555555582046509, - "716": 0.5109477043151855, - "717": 0.5410130620002747, - "718": 0.5320261716842651, - "719": 0.47647058963775635, - "720": 0.46307188272476196, - "721": 0.3839869201183319, - "722": 0.40163397789001465, - "723": 0.44607841968536377, - "724": 0.474346399307251, - "725": 0.5047385692596436, - "726": 0.5429738759994507, - "727": 0.5611110925674438, - "728": 0.5542483925819397, - "729": 0.5073529481887817, - "730": 0.48398691415786743, - "731": 0.43921568989753723, - "732": 0.4315359592437744, - "733": 0.4593137204647064, - "734": 0.5137255191802979, - "735": 0.548202633857727, - "736": 0.5539215803146362, - "737": 0.525653600692749, - "738": 0.505065381526947, - "739": 0.49656862020492554, - "740": 0.4892156720161438, - "741": 0.5009803771972656, - "742": 0.45849671959877014, - "743": 0.4506535828113556, - "744": 0.5245097875595093, - "745": 0.558169960975647, - "746": 0.5392156839370728, - "747": 0.5168300867080688, - "748": 0.5013071894645691, - "749": 0.46977123618125916, - "750": 0.4228758215904236, - "751": 0.4888888895511627, - "752": 0.4444444477558136, - "753": 0.4387255012989044, - "754": 0.5176470875740051, - "755": 0.5284313559532166, - "756": 0.49983659386634827, - "757": 0.5178104639053345, - "758": 0.525653600692749, - "759": 0.4732026159763336, - "760": 0.3915032744407654, - "761": 0.48594769835472107, - "762": 0.45310458540916443, - "763": 0.44558823108673096, - "764": 0.48366013169288635, - "765": 0.4766339957714081, - "766": 0.4570261538028717, - "767": 0.46846404671669006, - "768": 0.4964052140712738, - "769": 0.48627451062202454, - "770": 0.4192810356616974, - "771": 0.5326797366142273, - "772": 0.505065381526947, - "773": 0.47287580370903015, - "774": 0.4547385573387146, - "775": 0.44999998807907104, - "776": 0.45996731519699097, - "777": 0.4660130739212036, - "778": 0.4812091588973999, - "779": 0.48137253522872925, - "780": 0.46176469326019287, - "781": 0.5882353186607361, - "782": 0.5591503381729126, - "783": 0.5181372761726379, - "784": 0.484640508890152, - "785": 0.4658496677875519, - "786": 0.46437907218933105, - "787": 0.4722222089767456, - "788": 0.4851307272911072, - "789": 0.488725483417511, - "790": 0.516339898109436, - "791": 0.6111111044883728, - "792": 0.578104555606842, - "793": 0.5408496856689453, - "794": 0.5160130858421326, - "795": 0.48267972469329834, - "796": 0.45359477400779724, - "797": 0.4503268003463745, - "798": 0.47418299317359924, + "701": 0.42712417244911194, + "702": 0.43039214611053467, + "703": 0.40441176295280457, + "704": 0.3795751631259918, + "705": 0.4264705777168274, + "706": 0.484640508890152, + "707": 0.4843137264251709, + "708": 0.4277777671813965, + "709": 0.4125817120075226, + "710": 0.4874182939529419, + "711": 0.41274508833885193, + "712": 0.4217320382595062, + "713": 0.42679738998413086, + "714": 0.42009803652763367, + "715": 0.45669934153556824, + "716": 0.5117647051811218, + "717": 0.5416666865348816, + "718": 0.532516360282898, + "719": 0.47777777910232544, + "720": 0.4642156958580017, + "721": 0.38562092185020447, + "722": 0.4024509787559509, + "723": 0.4470588266849518, + "724": 0.47598040103912354, + "725": 0.5060457587242126, + "726": 0.5446078181266785, + "727": 0.5620915293693542, + "728": 0.5547385811805725, + "729": 0.5089869499206543, + "730": 0.48562091588974, + "731": 0.44117647409439087, + "732": 0.4331699311733246, + "733": 0.4611110985279083, + "734": 0.5153594613075256, + "735": 0.549346387386322, + "736": 0.5553921461105347, + "737": 0.5266339778900146, + "738": 0.5055555701255798, + "739": 0.4972222149372101, + "740": 0.4900326728820801, + "741": 0.5039215683937073, + "742": 0.46078431606292725, + "743": 0.4524509906768799, + "744": 0.5261437892913818, + "745": 0.5588235259056091, + "746": 0.5401960611343384, + "747": 0.5181372761726379, + "748": 0.5017973780632019, + "749": 0.4700980484485626, + "750": 0.4230392277240753, + "751": 0.4915032684803009, + "752": 0.4462418258190155, + "753": 0.4400326907634735, + "754": 0.5189542770385742, + "755": 0.529411792755127, + "756": 0.5014705657958984, + "757": 0.519444465637207, + "758": 0.5264706015586853, + "759": 0.4740196168422699, + "760": 0.3919934630393982, + "761": 0.4869281053543091, + "762": 0.45408496260643005, + "763": 0.4467320144176483, + "764": 0.4848039150238037, + "765": 0.47777777910232544, + "766": 0.45816993713378906, + "767": 0.4696078300476074, + "768": 0.4973856210708618, + "769": 0.4869281053543091, + "770": 0.41993463039398193, + "771": 0.5336601138114929, + "772": 0.5057189464569092, + "773": 0.47369280457496643, + "774": 0.4560457468032837, + "775": 0.45098039507865906, + "776": 0.46078431606292725, + "777": 0.4673202633857727, + "778": 0.482516348361969, + "779": 0.4820261299610138, + "780": 0.4624182879924774, + "781": 0.5898692607879639, + "782": 0.5602940917015076, + "783": 0.5187908411026001, + "784": 0.4852941036224365, + "785": 0.46666666865348816, + "786": 0.4655228853225708, + "787": 0.47418299317359924, + "788": 0.4869281053543091, + "789": 0.489705890417099, + "790": 0.5171568393707275, + "791": 0.6130719184875488, + "792": 0.5797385573387146, + "793": 0.5415032505989075, + "794": 0.516339898109436, + "795": 0.48366013169288635, + "796": 0.45490196347236633, + "797": 0.45261436700820923, + "798": 0.4761437773704529, "799": 0.0, "800": 0.0, - "801": 0.4019607901573181, - "802": 0.4205882251262665, - "803": 0.4488562047481537, - "804": 0.4647058844566345, - "805": 0.4758169949054718, - "806": 0.4936274588108063, - "807": 0.5212418437004089, - "808": 0.5006536245346069, - "809": 0.4936274588108063, + "801": 0.40294116735458374, + "802": 0.42140522599220276, + "803": 0.44950979948043823, + "804": 0.46535947918891907, + "805": 0.4767973721027374, + "806": 0.4946078360080719, + "807": 0.5223855972290039, + "808": 0.501960813999176, + "809": 0.49444442987442017, "810": 0.0, - "811": 0.41470587253570557, - "812": 0.41503268480300903, - "813": 0.43251633644104004, - "814": 0.47156861424446106, - "815": 0.5086601376533508, - "816": 0.5186274647712708, - "817": 0.5308823585510254, - "818": 0.530718982219696, - "819": 0.49297386407852173, - "820": 0.49166667461395264, - "821": 0.4058823585510254, - "822": 0.3982026278972626, - "823": 0.4341503381729126, - "824": 0.4991829991340637, - "825": 0.5434640645980835, - "826": 0.563725471496582, - "827": 0.5375816822052002, - "828": 0.5109477043151855, - "829": 0.49787580966949463, - "830": 0.4838235378265381, - "831": 0.4418300688266754, - "832": 0.41601306200027466, - "833": 0.45310458540916443, - "834": 0.5473856329917908, - "835": 0.5805555582046509, + "811": 0.4156862795352936, + "812": 0.4161764681339264, + "813": 0.43349674344062805, + "814": 0.4725490212440491, + "815": 0.5099673271179199, + "816": 0.5200980305671692, + "817": 0.5321895480155945, + "818": 0.5316993594169617, + "819": 0.4942810535430908, + "820": 0.49281045794487, + "821": 0.406862735748291, + "822": 0.39901959896087646, + "823": 0.4349673092365265, + "824": 0.5004901885986328, + "825": 0.5449346303939819, + "826": 0.5653594732284546, + "827": 0.5388888716697693, + "828": 0.5117647051811218, + "829": 0.499019593000412, + "830": 0.48496732115745544, + "831": 0.4431372582912445, + "832": 0.41683006286621094, + "833": 0.4539215564727783, + "834": 0.5486928224563599, + "835": 0.5816993713378906, "836": 0.0, - "837": 0.5519607663154602, - "838": 0.5135620832443237, - "839": 0.5099673271179199, - "840": 0.5047385692596436, - "841": 0.49575161933898926, - "842": 0.4601307213306427, - "843": 0.48006534576416016, - "844": 0.5872548818588257, - "845": 0.5975490212440491, - "846": 0.5377451181411743, - "847": 0.540032684803009, - "848": 0.5405228734016418, - "849": 0.4937908351421356, - "850": 0.4573529362678528, - "851": 0.4937908351421356, - "852": 0.4665032625198364, - "853": 0.4766339957714081, - "854": 0.5532679557800293, - "855": 0.5442810654640198, - "856": 0.4820261299610138, - "857": 0.5049019455909729, - "858": 0.5305555462837219, - "859": 0.48905229568481445, - "860": 0.43790850043296814, - "861": 0.4673202633857727, - "862": 0.4552287459373474, - "863": 0.47287580370903015, - "864": 0.5006536245346069, - "865": 0.4748365879058838, - "866": 0.45800653100013733, - "867": 0.46666666865348816, - "868": 0.4776143729686737, - "869": 0.4766339957714081, - "870": 0.4369280934333801, - "871": 0.47875815629959106, - "872": 0.4745098054409027, - "873": 0.5034313797950745, - "874": 0.5080065131187439, - "875": 0.4776143729686737, - "876": 0.49787580966949463, - "877": 0.4982026219367981, - "878": 0.4735293984413147, - "879": 0.4650326669216156, - "880": 0.4459150433540344, + "837": 0.0, + "838": 0.5142157077789307, + "839": 0.5106208920478821, + "840": 0.5053921341896057, + "841": 0.4977124035358429, + "842": 0.4614379107952118, + "843": 0.48104575276374817, + "844": 0.5885620713233948, + "845": 0.5983660221099854, + "846": 0.5388888716697693, + "847": 0.5411764979362488, + "848": 0.5410130620002747, + "849": 0.4946078360080719, + "850": 0.45800653100013733, + "851": 0.4954248368740082, + "852": 0.4676470458507538, + "853": 0.4779411852359772, + "854": 0.5547385811805725, + "855": 0.5452614426612854, + "856": 0.4834967255592346, + "857": 0.5063725709915161, + "858": 0.5315359234809875, + "859": 0.4900326728820801, + "860": 0.4385620951652527, + "861": 0.46830064058303833, + "862": 0.4560457468032837, + "863": 0.474346399307251, + "864": 0.5021241903305054, + "865": 0.47565358877182007, + "866": 0.45898693799972534, + "867": 0.468137264251709, + "868": 0.47908496856689453, + "869": 0.477450966835022, + "870": 0.4375816881656647, + "871": 0.4794117510318756, + "872": 0.47516340017318726, + "873": 0.5045751929283142, + "874": 0.5094771385192871, + "875": 0.4785947799682617, + "876": 0.4986928105354309, + "877": 0.0, + "878": 0.4748365879058838, + "879": 0.46568626165390015, + "880": 0.4467320144176483, "881": 0.0, - "882": 0.5060457587242126, - "883": 0.5120915174484253, - "884": 0.508169949054718, - "885": 0.48725488781929016, - "886": 0.5160130858421326, - "887": 0.5210784077644348, - "888": 0.49133986234664917, - "889": 0.5001633763313293, - "890": 0.536928117275238, + "882": 0.5071895718574524, + "883": 0.5129085183143616, + "884": 0.5089869499206543, + "885": 0.48807188868522644, + "886": 0.5169934630393982, + "887": 0.522549033164978, + "888": 0.49248364567756653, + "889": 0.5011438131332397, + "890": 0.538071870803833, "891": 0.0, - "892": 0.5169934630393982, - "893": 0.49803921580314636, - "894": 0.48954248428344727, - "895": 0.4735293984413147, - "896": 0.491830050945282, - "897": 0.4977124035358429, - "898": 0.48496732115745544, - "899": 0.5271241664886475, + "892": 0.5189542770385742, + "893": 0.499019593000412, + "894": 0.4901960790157318, + "895": 0.47418299317359924, + "896": 0.49281045794487, + "897": 0.49934640526771545, + "898": 0.48594769835472107, + "899": 0.5284313559532166, "900": 0.0, "901": 0.44509804248809814, - "902": 0.4575163424015045, - "903": 0.4941176474094391, - "904": 0.5343137383460999, - "905": 0.5362744927406311, - "906": 0.5111111402511597, - "907": 0.5248365998268127, - "908": 0.5539215803146362, - "909": 0.5882353186607361, + "902": 0.4578431248664856, + "903": 0.49477124214172363, + "904": 0.5346405506134033, + "905": 0.536928117275238, + "906": 0.5124183297157288, + "907": 0.5258169770240784, + "908": 0.5552287697792053, + "909": 0.5895425081253052, "910": 0.0, "911": 0.4490196108818054, - "912": 0.44248366355895996, - "913": 0.46209150552749634, - "914": 0.5245097875595093, - "915": 0.5604575276374817, - "916": 0.5310457348823547, - "917": 0.5140522718429565, - "918": 0.5300653576850891, - "919": 0.5232025980949402, - "920": 0.5189542770385742, - "921": 0.44967320561408997, - "922": 0.43300652503967285, - "923": 0.46339869499206543, - "924": 0.5392156839370728, + "912": 0.4434640407562256, + "913": 0.46339869499206543, + "914": 0.5254902243614197, + "915": 0.5620915293693542, + "916": 0.5333333611488342, + "917": 0.5160130858421326, + "918": 0.5313725471496582, + "919": 0.5245097875595093, + "920": 0.5202614665031433, + "921": 0.4503268003463745, + "922": 0.43398693203926086, + "923": 0.46437907218933105, + "924": 0.5408496856689453, "925": 0.0, - "926": 0.5669934749603271, - "927": 0.506209135055542, - "928": 0.47385621070861816, - "929": 0.477450966835022, - "930": 0.4614379107952118, - "931": 0.45816993713378906, - "932": 0.43921568989753723, - "933": 0.48954248428344727, - "934": 0.5898692607879639, + "926": 0.5686274766921997, + "927": 0.5078431367874146, + "928": 0.4748365879058838, + "929": 0.47843137383461, + "930": 0.4624182879924774, + "931": 0.45947712659835815, + "932": 0.4398692846298218, + "933": 0.48986926674842834, + "934": 0.591176450252533, "935": 0.0, "936": 0.0, "937": 0.0, - "938": 0.5094771385192871, - "939": 0.501960813999176, - "940": 0.484640508890152, + "938": 0.5101307034492493, + "939": 0.5029411911964417, + "940": 0.4852941036224365, "941": 0.0, "942": 0.0, "943": 0.0, "944": 0.0, - "945": 0.628104567527771, - "946": 0.556209146976471, + "945": 0.6290849447250366, + "946": 0.5571895241737366, "947": 0.0, "948": 0.0, - "949": 0.5166666507720947, - "950": 0.4718954265117645, + "949": 0.5179738402366638, + "950": 0.4732026159763336, "951": 0.0, - "952": 0.5039215683937073, - "953": 0.5274509787559509, - "954": 0.5846405029296875, - "955": 0.5607843399047852, - "956": 0.4892156720161438, - "957": 0.5107843279838562, - "958": 0.5385621190071106, - "959": 0.5052287578582764, - "960": 0.4748365879058838, - "961": 0.45816993713378906, - "962": 0.4732026159763336, - "963": 0.516339898109436, - "964": 0.5264706015586853, - "965": 0.48366013169288635, - "966": 0.47058823704719543, - "967": 0.47843137383461, - "968": 0.4732026159763336, - "969": 0.4722222089767456, - "970": 0.45686274766921997, - "971": 0.4503268003463745, - "972": 0.46568626165390015, - "973": 0.5248365998268127, - "974": 0.530718982219696, - "975": 0.4950980246067047, - "976": 0.5303921699523926, + "952": 0.5049019455909729, + "953": 0.5290849804878235, + "954": 0.5859476923942566, + "955": 0.5617647171020508, + "956": 0.49084967374801636, + "957": 0.5124183297157288, + "958": 0.5398693084716797, + "959": 0.5065359473228455, + "960": 0.4758169949054718, + "961": 0.4588235318660736, + "962": 0.47418299317359924, + "963": 0.5179738402366638, + "964": 0.5274509787559509, + "965": 0.4843137264251709, + "966": 0.4718954265117645, + "967": 0.48006534576416016, + "968": 0.4748365879058838, + "969": 0.4732026159763336, + "970": 0.4575163424015045, + "971": 0.45098039507865906, + "972": 0.46666666865348816, + "973": 0.5261437892913818, + "974": 0.5316993594169617, + "975": 0.49575161933898926, + "976": 0.5310457348823547, "977": 0.0, - "978": 0.49575161933898926, - "979": 0.4650326669216156, - "980": 0.4320261478424072, + "978": 0.4964052140712738, + "979": 0.46568626165390015, + "980": 0.43300652503967285, "981": 0.0, - "982": 0.4771241843700409, - "983": 0.4915032684803009, - "984": 0.5022875666618347, - "985": 0.5107843279838562, + "982": 0.4781045615673065, + "983": 0.49215686321258545, + "984": 0.5029411911964417, + "985": 0.5114378929138184, "986": 0.0, "987": 0.0, - "988": 0.5199346542358398, - "989": 0.5218954086303711, - "990": 0.543790876865387, + "988": 0.5202614665031433, + "989": 0.5228758454322815, + "990": 0.545098066329956, "991": 0.0, "992": 0.0, - "993": 0.45098039507865906, - "994": 0.4601307213306427, - "995": 0.49477124214172363, - "996": 0.5320261716842651, - "997": 0.5228758454322815, + "993": 0.4516339898109436, + "994": 0.46078431606292725, + "995": 0.4954248368740082, + "996": 0.5333333611488342, + "997": 0.5241830348968506, "998": 0.5078431367874146, - "999": 0.5620915293693542, + "999": 0.5633987188339233, "1000": 0.0 } } diff --git a/tests/unittests/i2pp/core/test_run.py b/tests/unittests/i2pp/core/test_run.py index cd2c2e4..ee49b11 100644 --- a/tests/unittests/i2pp/core/test_run.py +++ b/tests/unittests/i2pp/core/test_run.py @@ -2,6 +2,7 @@ from unittest.mock import MagicMock, patch +import numpy as np import pytest from i2pp.core.run import run_i2pp @@ -62,7 +63,9 @@ def minimal_valid_config(tmp_path): @patch("i2pp.core.run.export_data") @patch("i2pp.core.run.visualize_smoothing") @patch("i2pp.core.run.visualize_results") +@patch("i2pp.core.run.create_mesh_mask") def test_run_i2pp_with_minimal_valid_config( + mock_create_mesh_mask, mock_visualize_results, mock_visualize_smoothing, mock_export_data, @@ -95,6 +98,7 @@ def test_run_i2pp_with_minimal_valid_config( # Assertions to ensure key steps were called mock_verify_discretization.assert_called_once() mock_verify_image.assert_called_once() + mock_create_mesh_mask.assert_called_once() mock_smooth_data.assert_not_called() mock_visualize_smoothing.assert_not_called() mock_interpolate.assert_called_once() @@ -172,7 +176,9 @@ def maximal_valid_config(tmp_path): @patch("i2pp.core.run.export_data") @patch("i2pp.core.run.visualize_smoothing") @patch("i2pp.core.run.visualize_results") +@patch("i2pp.core.run.create_mesh_mask") def test_run_i2pp_with_maximal_valid_config( + mock_create_mesh_mask, mock_visualize_results, mock_visualize_smoothing, mock_export_data, @@ -189,12 +195,12 @@ def test_run_i2pp_with_maximal_valid_config( mock_verify_discretization.return_value = mock_discretization mock_image_data = MagicMock() - mock_image_data.pixel_data = [[0]] + mock_image_data.pixel_data = np.array([[0, 1], [2, 3]]) mock_image_data.pixel_range = (0, 255) mock_image_data.pixel_type = "uint8" mock_verify_image.return_value = mock_image_data - mock_smooth_data.side_effect = lambda pixel_data, area: pixel_data + mock_smooth_data.side_effect = lambda *args, **kwargs: args[0] mock_interpolate.return_value = "interpolated_elements" mock_transform_data.return_value = "transformed_data" @@ -205,6 +211,7 @@ def test_run_i2pp_with_maximal_valid_config( # Assertions to ensure key steps were called mock_verify_discretization.assert_called_once() mock_verify_image.assert_called_once() + mock_create_mesh_mask.assert_called_once() mock_smooth_data.assert_called_once() mock_visualize_smoothing.assert_called_once() mock_interpolate.assert_called_once() diff --git a/tests/unittests/i2pp/core/test_utilities.py b/tests/unittests/i2pp/core/test_utilities.py index 1fbf221..9908880 100644 --- a/tests/unittests/i2pp/core/test_utilities.py +++ b/tests/unittests/i2pp/core/test_utilities.py @@ -3,6 +3,7 @@ import numpy as np from i2pp.core.discretization_readers.discretization_reader import BoundingBox from i2pp.core.utilities import ( + create_mesh_mask, find_mins_maxs, get_node_position_of_element, make_json_serializable, @@ -88,7 +89,7 @@ def test_smoothing_dicom(): ) pixel_data = np.array([array_slice1, array_slice2, array_slice3]) - pixel_data_smoothed = smooth_data(pixel_data, 3) + pixel_data_smoothed = smooth_data(pixel_data, 3, mask=None) assert pixel_data_smoothed[1][1][1] == 5 @@ -125,7 +126,7 @@ def test_smoothing_rgb(): pixel_data = np.array([array_slice1, array_slice2, array_slice3]) - pixel_data_smoothed = smooth_data(pixel_data, 3) + pixel_data_smoothed = smooth_data(pixel_data, 3, mask=None) assert np.array_equal(pixel_data_smoothed[1][1][1], np.array([3, 3, 2])) @@ -173,3 +174,93 @@ def test_make_json_serializable_non_numpy(): value = {"x": 1, "y": [2, 3]} result = make_json_serializable(value) assert result == value + + +def test_create_mesh_mask_cube_covers_entire_grid(): + """Create a cube mesh that spans the whole 3x3x3 grid and expect a full + True mask.""" + + # Dummy discretization: nodes are the 8 cube corners at coordinates 0 or 2 + class DummyNodes: + """Dummy class representing nodes in a discretization.""" + + def __init__(self, ids, coords): + """Initialize DummyNodes. + + Args: + ids (list): List of node IDs. + coords (list): List of node coordinates. + """ + self.ids = np.array(ids) + self.coords = np.array(coords, dtype=float) + + class DummyElement: + """Dummy class representing an element in a discretization.""" + + def __init__(self, node_ids): + """Initialize DummyElement. + + Args: + node_ids (list): List of node IDs forming the element. + """ + self.node_ids = np.array(node_ids) + + class DummyDiscretization: + """Dummy class representing a discretization.""" + + def __init__(self, nodes, elements): + """Initialize DummyDiscretization. + + Args: + nodes (DummyNodes): Nodes in the discretization. + elements (list): Elements in the discretization. + """ + self.nodes = nodes + self.elements = elements + + class DummyGridCoords: + """Dummy class representing grid coordinates.""" + + def __init__(self, coords): + """Initialize DummyGridCoords. + + Args: + coords (list): List of grid coordinates. + """ + self.slice = coords # z + self.row = coords # y + self.col = coords # x + + class DummyImage: + """Dummy class representing an image.""" + + def __init__(self): + """Initialize DummyImage.""" + self.pixel_data = np.zeros((3, 3, 3), dtype=np.float32) + self.orientation = np.eye(3, dtype=float) + self.position = np.zeros(3, dtype=float) + coords = np.array([0.0, 1.0, 2.0], dtype=float) + self.grid_coords = DummyGridCoords(coords) + + # Define cube corners + corners = [ + [0, 0, 0], + [2, 0, 0], + [0, 2, 0], + [0, 0, 2], + [2, 2, 0], + [2, 0, 2], + [0, 2, 2], + [2, 2, 2], + ] + node_ids = list(range(1, 9)) + nodes = DummyNodes(node_ids, corners) + # Single element with all 8 nodes (convex hull == cube) + elements = [DummyElement(node_ids)] + discretization = DummyDiscretization(nodes, elements) + image = DummyImage() + + mask = create_mesh_mask(discretization, image) + + assert mask.shape == (3, 3, 3) + assert np.all(mask) From 2b40e6a5b8ce6bcd21ae87f28d3e107e8fa083ef Mon Sep 17 00:00:00 2001 From: Saskia Gloy Date: Tue, 30 Dec 2025 11:14:41 +0100 Subject: [PATCH 04/15] Add weighted interpolation and tests --- .../discretization_reader.py | 3 + .../fourc_yaml_reader.py | 11 +- src/i2pp/core/interpolate_element_data.py | 14 +- .../interpolators/interpolator_all_voxel.py | 219 ++++++++++-- .../core/interpolators/interpolator_types.py | 54 +-- .../test_interpolator_all_voxel.py | 317 +++++++++++++++++- .../interpolators/test_interpolator_types.py | 47 ++- 7 files changed, 604 insertions(+), 61 deletions(-) diff --git a/src/i2pp/core/discretization_readers/discretization_reader.py b/src/i2pp/core/discretization_readers/discretization_reader.py index 67243f2..6a5a045 100644 --- a/src/i2pp/core/discretization_readers/discretization_reader.py +++ b/src/i2pp/core/discretization_readers/discretization_reader.py @@ -31,10 +31,13 @@ class Nodes: coords (np.ndarray): An (N, 3) array containing the (x, y, z) world coordinates of each node. ids (np.ndarray): An array of IDs for each node. + weights (Optional[np.ndarray]): Optional per-node weights, + e.g., low weights for surface nodes. """ coords: np.ndarray ids: np.ndarray + weights: Optional[np.ndarray] = None @dataclass diff --git a/src/i2pp/core/discretization_readers/fourc_yaml_reader.py b/src/i2pp/core/discretization_readers/fourc_yaml_reader.py index da0179c..a2a2160 100644 --- a/src/i2pp/core/discretization_readers/fourc_yaml_reader.py +++ b/src/i2pp/core/discretization_readers/fourc_yaml_reader.py @@ -104,10 +104,12 @@ def load_discretization( nodes_coords = [] node_ids = [] + nodes_weights = [] for node in raw_dis.nodes: nodes_coords.append(node.coords) node_ids.append(node.id) + nodes_weights.append(5.0) elements = [] @@ -126,13 +128,20 @@ def load_discretization( surf_node_ids = [] for node in surf.nodes: surf_node_ids.append(node.id) + # get position of node in node_ids list and set weight + position = node_ids.index(node.id) + nodes_weights[position] = 0.0 surfaces.append( Surface(node_ids=np.array(surf_node_ids), id=surf.id) ) dis = Discretization( - nodes=Nodes(coords=np.array(nodes_coords), ids=np.array(node_ids)), + nodes=Nodes( + coords=np.array(nodes_coords), + ids=np.array(node_ids), + weights=np.array(nodes_weights), + ), elements=elements, surfaces=surfaces, ) diff --git a/src/i2pp/core/interpolate_element_data.py b/src/i2pp/core/interpolate_element_data.py index b00f830..f27c5ce 100644 --- a/src/i2pp/core/interpolate_element_data.py +++ b/src/i2pp/core/interpolate_element_data.py @@ -9,7 +9,11 @@ def interpolate_image_to_discretization( - dis: Discretization, image_data: ImageData, interpolation_method: str + dis: Discretization, + image_data: ImageData, + interpolation_method: str, + *, + filter_outliers: bool = False, ) -> list[Element]: """Performs interpolation of image data onto the FEM Discretization based on the specified interpolation method. @@ -34,13 +38,17 @@ def interpolate_image_to_discretization( assigning pixel values to the elements. This should match one of the `InterpolationType` enum values (e.g., "nodes", "elementcenter", "allvoxels"). + filter_outliers (bool, optional): Whether to filter outliers during + interpolation. Defaults to False. Returns: list[Element]: A list of FEM elements with interpolated pixel data. """ - enum_interpolation_method = InterpolationType(interpolation_method) - interpolator = enum_interpolation_method.get_interpolator()() + # Use factory to configure weighted vs unweighted and outlier filtering + interpolator = enum_interpolation_method.create_interpolator( + filter_outliers=filter_outliers + ) return interpolator.compute_element_data(dis, image_data) diff --git a/src/i2pp/core/interpolators/interpolator_all_voxel.py b/src/i2pp/core/interpolators/interpolator_all_voxel.py index 3b6447c..2c7f74a 100644 --- a/src/i2pp/core/interpolators/interpolator_all_voxel.py +++ b/src/i2pp/core/interpolators/interpolator_all_voxel.py @@ -15,15 +15,27 @@ class InterpolatorAllVoxel(Interpolator): - """Subclass of Interpolator for mapping 3D image data to finite element - mesh elements. - - This class extends Interpolator and specializes in assigning pixel values - from 3D image data to finite element mesh elements by computing the mean - of all voxels contained within each element. This approach is used when - `interpolation_method` is set to "allvoxels". + """Interpolator for mapping 3D image data to finite element mesh elements. + + This class supports both unweighted and node-weighted voxel mean + calculations, controlled via the 'mode' parameter. It assigns pixel + values from 3D image data to finite element mesh elements by computing + the mean of all voxels within each element. + This functionality is used when the `interpolation_method` is set to + "allvoxels" or "allvoxels_weighted". """ + def __init__( + self, + *args, + filter_outliers: bool = False, + mode: str = "allvoxels", + **kwargs, + ): + super().__init__(*args, **kwargs) + self._filter_outliers_enabled = filter_outliers + self._mode = mode # "allvoxels" or "allvoxels_weighted" + def _search_bounding_box( self, grid_coords: GridCoords, element_grid_coords: np.ndarray ) -> Tuple[np.ndarray, np.ndarray, np.ndarray]: @@ -98,8 +110,95 @@ def _is_inside_element(self, point: np.ndarray, hull: ConvexHull): return np.all(A @ point + b <= 0) + def _filter_outliers_modified_zscore( + self, values: np.ndarray, threshold: float = 3.5 + ) -> np.ndarray: + """Identify outliers using the Modified Z-Score method. + + Args: + values (np.ndarray): + Array of voxel values (1D or 2D, axis=0 = voxels) + threshold (float): + Cutoff for Modified Z-Score, default 3.5 + + Returns: + mask (np.ndarray): Boolean array, True = value is not an outlier + """ + values = np.asarray(values) + median = np.median(values, axis=0) + + # Median Absolute Deviation (MAD) + mad = np.median(np.abs(values - median), axis=0) + mad = np.maximum(mad, 1e-10) # avoid division by zero + + # Modified Z-Score + modified_z = 0.6745 * (values - median) / mad + + # Mask: True = not an outlier + mask = np.abs(modified_z) <= threshold + + return mask + + def _weighted_voxel_mean( + self, + voxels_phys: np.ndarray, + values: np.ndarray, + element_node_phys: np.ndarray, + node_weights: np.ndarray, + ) -> float: + """Calculate a node-weighted mean of voxel values, with optional + outlier filtering. + + Computes the weighted mean of voxel values within an element, + where the weight is determined based on the distances between the voxel + and the element's nodes, scaled by the node weights. Outliers can be + excluded using the Modified Z-Score method. + + Args: + voxels_phys (np.ndarray): + Physical coordinate of voxels within the element (N_voxel x 3) + values (np.ndarray): + Corresponding voxel values (N_voxels x ...). + element_node_phys (np.ndarray): + Physical coordinates of the element's nodes (N_nodes x 3). + node_weights (np.ndarray): + Weights assigned to the nodes (N_nodes,). + + Returns: + float: The weighted mean of the voxel values. + """ + + distances = np.linalg.norm( + element_node_phys[:, np.newaxis, :] + - voxels_phys[np.newaxis, :, :], + axis=2, + ) + distances = np.maximum(distances, 1e-10) + + voxel_weights = np.sum(node_weights[:, np.newaxis] / distances, axis=0) + + if len(values) > 5: + mask = self._filter_outliers_modified_zscore(values) + filtered_values = values[mask] + filtered_weights = voxel_weights[mask] + + # Fallback if all voxels are filtered out + if len(filtered_values) == 0: + filtered_values = values + filtered_weights = voxel_weights + + else: + # No outlier filtering for small voxel counts + return np.average(values, weights=voxel_weights, axis=0) + + # Compute weighted mean + return np.average(filtered_values, weights=filtered_weights, axis=0) + def _get_data_of_element( - self, element_node_grid_coords: np.ndarray, image_data: ImageData + self, + element_node_grid_coords: np.ndarray, + image_data: ImageData, + node_weights_current: np.ndarray | None = None, ) -> np.ndarray: """Computes the representative pixel value for a given element based on its nodes in grid coordinates. @@ -118,6 +217,8 @@ def _get_data_of_element( the element's nodes. image_data (ImageData): Image data containing voxel coordinates and pixel values. + node_weights_current (np.ndarray | None): + Weights assigned to the nodes. Returns: np.ndarray: The mean pixel value of all voxels inside the element. @@ -126,12 +227,13 @@ def _get_data_of_element( If all nodes are outside the grid, returns `np.nan`. """ - data = [] + # Lists for collection (typed for mypy) + voxels_phys: list[np.ndarray] = [] + data_list: list[np.ndarray] = [] slice_indices, row_indices, col_indices = self._search_bounding_box( image_data.grid_coords, element_node_grid_coords ) - hull = ConvexHull(element_node_grid_coords) for i in slice_indices: @@ -144,18 +246,89 @@ def _get_data_of_element( image_data.grid_coords.col[k], ] ) - if self._is_inside_element(grid_coord, hull): + voxels_phys.append( + np.array( + [ + image_data.grid_coords.slice[i], + image_data.grid_coords.row[j], + image_data.grid_coords.col[k], + ] + ) + ) + data_list.append(image_data.pixel_data[i, j, k]) + + if len(voxels_phys) > 0: + voxels_phys_np = np.asarray(voxels_phys) + data = np.asarray( + data_list + ) # ensure ndarray for boolean masks and vector ops + + if self._mode == "allvoxels": + if self._filter_outliers_enabled and len(data) > 5: + mask = self._filter_outliers_modified_zscore(data) + filtered = data[mask] + mean_val = ( + np.mean(filtered, axis=0) + if len(filtered) > 0 + else np.mean(data, axis=0) + ) + else: + mean_val = np.mean(data, axis=0) + return ( + np.array([mean_val]) + if image_data.pixel_type.num_values == 1 + else np.full(image_data.pixel_type.num_values, mean_val) + ) + + if self._mode == "allvoxels_weighted": + element_node_phys = element_node_grid_coords + if node_weights_current is None: + weighted = np.mean(data, axis=0) + else: + weighted = ( + self._weighted_voxel_mean( + voxels_phys_np, + data, + element_node_phys, + node_weights_current, + ) + if self._filter_outliers_enabled + else np.average( + data, + weights=np.sum( + node_weights_current[:, np.newaxis] + / np.maximum( + np.linalg.norm( + element_node_phys[:, np.newaxis, :] + - voxels_phys_np[np.newaxis, :, :], + axis=2, + ), + 1e-10, + ), + axis=0, + ), + axis=0, + ) + ) - data.append(image_data.pixel_data[i, j, k]) + return ( + np.array([weighted]) + if image_data.pixel_type.num_values == 1 + else np.full(image_data.pixel_type.num_values, weighted) + ) + + # Unknown mode fallback + mean_val = np.mean(data, axis=0) + return ( + np.array([mean_val]) + if image_data.pixel_type.num_values == 1 + else np.full(image_data.pixel_type.num_values, mean_val) + ) - if data: - return np.mean(data, axis=0) else: self.backup_interpolation += 1 - element_center = np.mean(element_node_grid_coords, axis=0) - return self.interpolate_image_values_to_points( element_center, image_data )[0] @@ -199,16 +372,24 @@ def compute_element_data( total=len(dis.elements), desc="Element values", ): - element_node_grid_coords = node_grid_coords[node_positions[i]] + # Cast weights to ndarray when present to satisfy type checker + node_weights_current = None + if getattr(dis.nodes, "weights", None) is not None: + # mypy: explicitly cast Optional[Any] to np.ndarray + from typing import cast + + weights_nd = cast(np.ndarray, dis.nodes.weights) + node_weights_current = weights_nd[node_positions[i]] ele.data = self._get_data_of_element( - element_node_grid_coords, image_data + element_node_grid_coords, + image_data, + node_weights_current=node_weights_current, ) if np.all(np.isnan(ele.data)): self.nan_elements += 1 self._log_interpolation_warnings() - return dis.elements diff --git a/src/i2pp/core/interpolators/interpolator_types.py b/src/i2pp/core/interpolators/interpolator_types.py index d4b0651..e7400de 100644 --- a/src/i2pp/core/interpolators/interpolator_types.py +++ b/src/i2pp/core/interpolators/interpolator_types.py @@ -1,7 +1,6 @@ """Interpolation type definitions and handling.""" from enum import Enum -from typing import Type from i2pp.core.interpolators.interpolator import Interpolator from i2pp.core.interpolators.interpolator_all_voxel import InterpolatorAllVoxel @@ -21,40 +20,43 @@ class InterpolationType(Enum): Attributes: NODES (str): Represents the interpolation method where the pixel value is averaged over the nodes of the element. + NODES_WEIGHTED (str): Represents the interpolation method where the + pixel value is weighted over the nodes of the element. CENTER (str): Represents the interpolation method where the pixel value is based on the center of the element. ALLVOXELS (str): Represents the interpolation method where the pixel value is averaged over all voxels inside the element. + ALLVOXELS_WEIGHTED (str): Represents the interpolation method where the + pixel value is weighted over all voxels inside the element. + + Use create_interpolator to obtain a configured interpolator instance. """ NODES = "nodes" + NODES_WEIGHTED = "nodes_weighted" CENTER = "elementcenter" ALLVOXELS = "allvoxels" + ALLVOXELS_WEIGHTED = "allvoxels_weighted" - def get_interpolator(self) -> Type[Interpolator]: - """Retrieves the appropriate interpolation class based on the selected - interpolation method. - - This method returns the corresponding interpolator class for assigning - pixel values to FEM elements, depending on the current interpolation - method. Supported methods include interpolation at element nodes, - element centers, or averaging all voxels within an element. - - Returns: - Type[Interpolator]: The interpolator class that matches the - specified interpolation method. + def create_interpolator( + self, *, filter_outliers: bool = False + ) -> Interpolator: + """Factory that returns a configured interpolator instance. - Raises: - ValueError: If the interpolation method is not supported. + For ALLVOXELS/ALLVOXELS_WEIGHTED, this method configures the + unified InterpolatorAllVoxel with the appropriate mode and + optional outlier filtering. """ - - interpolator_map = { - InterpolationType.NODES: InterpolatorNodes, - InterpolationType.ALLVOXELS: InterpolatorAllVoxel, - InterpolationType.CENTER: InterpolatorCenter, - } - - if self not in interpolator_map: - raise ValueError(f"Unsupported interpolation method: {self}") - - return interpolator_map[self] + if self == InterpolationType.ALLVOXELS: + return InterpolatorAllVoxel( + mode="allvoxels", filter_outliers=filter_outliers + ) + if self == InterpolationType.ALLVOXELS_WEIGHTED: + return InterpolatorAllVoxel( + mode="allvoxels_weighted", filter_outliers=filter_outliers + ) + if self in (InterpolationType.NODES, InterpolationType.NODES_WEIGHTED): + return InterpolatorNodes() + if self == InterpolationType.CENTER: + return InterpolatorCenter() + raise ValueError(f"Unsupported interpolation method: {self}") diff --git a/tests/unittests/i2pp/core/interpolators/test_interpolator_all_voxel.py b/tests/unittests/i2pp/core/interpolators/test_interpolator_all_voxel.py index 67906c9..a670d3a 100644 --- a/tests/unittests/i2pp/core/interpolators/test_interpolator_all_voxel.py +++ b/tests/unittests/i2pp/core/interpolators/test_interpolator_all_voxel.py @@ -1,6 +1,16 @@ """Test Interpolator Routine.""" import numpy as np +from i2pp.core.discretization_readers.discretization_reader import ( + Discretization, +) +from i2pp.core.discretization_readers.discretization_reader import ( + Element as DisElement, +) +from i2pp.core.discretization_readers.discretization_reader import ( + Nodes, + Surface, +) from i2pp.core.image_readers.image_reader import ( GridCoords, ImageData, @@ -144,7 +154,9 @@ def test_get_data_of_element_element_in_grid_scalar(): col_coords = np.arange(5) grid_coords = GridCoords(slice_coords, row_coords, col_coords) - image_data = ImageData(pixel_data, grid_coords, [], [], PixelValueType.CT) + image_data = ImageData( + pixel_data, grid_coords, np.eye(3), np.zeros(3), PixelValueType.CT + ) interpolator = InterpolatorAllVoxel() @@ -178,7 +190,9 @@ def test_get_data_of_element_element_in_grid_RGB(): ) interpolator = InterpolatorAllVoxel() - image_data = ImageData(pixel_data, grid_coords, [], [], PixelValueType.RGB) + image_data = ImageData( + pixel_data, grid_coords, np.eye(3), np.zeros(3), PixelValueType.RGB + ) data = interpolator._get_data_of_element( element_node_grid_coords, image_data @@ -217,7 +231,9 @@ def test_get_data_of_element_element_not_in_grid(): col_coords = np.arange(5) grid_coords = GridCoords(slice_coords, row_coords, col_coords) - image_data = ImageData(pixel_data, grid_coords, [], [], PixelValueType.CT) + image_data = ImageData( + pixel_data, grid_coords, np.eye(3), np.zeros(3), PixelValueType.CT + ) interpolator = InterpolatorAllVoxel() @@ -254,7 +270,9 @@ def test_get_data_of_element_element_low_resolution_image(): col_coords = np.arange(5) * 3 grid_coords = GridCoords(slice_coords, row_coords, col_coords) - image_data = ImageData(pixel_data, grid_coords, [], [], PixelValueType.CT) + image_data = ImageData( + pixel_data, grid_coords, np.eye(3), np.zeros(3), PixelValueType.CT + ) interpolator = InterpolatorAllVoxel() @@ -264,3 +282,294 @@ def test_get_data_of_element_element_low_resolution_image(): interpolator._get_data_of_element(element, image_data), np.mean(interpol_point), ) + + +def _make_simple_discretization( + node_coords, elements_node_ids, node_weights=None +): + """Helper to create a minimal Discretization with nodes/elements and + weights.""" + node_ids = np.arange(len(node_coords)) + nodes = Nodes(ids=node_ids, coords=np.asarray(node_coords)) + if node_weights is None: + node_weights = np.ones(len(node_coords), dtype=float) + setattr(nodes, "weights", np.asarray(node_weights)) + + # Create Elements with required 'id' field; + # center_coords/data left as defaults + elements = [ + DisElement(node_ids=np.asarray(nids), id=idx) + for idx, nids in enumerate(elements_node_ids) + ] + + # Single Surface referencing all element ids + surfaces = [Surface(node_ids=np.array([], dtype=int), id=0)] + + return Discretization(nodes=nodes, elements=elements, surfaces=surfaces) + + +def test_unweighted_vs_weighted_mode_mean_differs_when_weights_bias(): + """Weighted mode should differ from unweighted mean when node weights bias + proximity.""" + # Image: small grid 3x3x3 with deterministic values + pixel_data = np.arange(27, dtype=float).reshape((3, 3, 3)) + grid_coords = GridCoords( + slice=np.arange(3), row=np.arange(3), col=np.arange(3) + ) + image_data = ImageData( + pixel_data, grid_coords, np.eye(3), np.zeros(3), PixelValueType.CT + ) + + # Element: cube corners spanning [0,1] in each axis + ele_nodes = np.array( + [ + [0, 0, 0], + [1, 0, 0], + [0, 1, 0], + [1, 1, 0], + [0, 0, 1], + [1, 0, 1], + [0, 1, 1], + [1, 1, 1], + ] + ) + + # Discretization with 8 nodes, one element using all nodes + dis = _make_simple_discretization( + node_coords=ele_nodes, + elements_node_ids=[np.arange(8)], + node_weights=np.array( + [10, 10, 10, 10, 1, 1, 1, 1], dtype=float + ), # bias towards lower z nodes + ) + + # Unweighted + interpolator_unweighted = InterpolatorAllVoxel( + mode="allvoxels", filter_outliers=False + ) + elems_unweighted = interpolator_unweighted.compute_element_data( + dis, image_data + ) + val_unweighted = elems_unweighted[0].data + + # Weighted + interpolator_weighted = InterpolatorAllVoxel( + mode="allvoxels_weighted", filter_outliers=False + ) + elems_weighted = interpolator_weighted.compute_element_data( + dis, image_data + ) + val_weighted = elems_weighted[0].data + + assert np.isfinite(val_unweighted).all() + assert np.isfinite(val_weighted).all() + # Expect a difference due to weight bias toward voxels near lower-z nodes + assert not np.allclose(val_unweighted, val_weighted) + + +def test_filter_outliers_reduces_extreme_values_in_unweighted_mode(): + """Outlier filtering should reduce influence of extreme voxel values.""" + # Construct image with an element that will include voxels; + # insert an extreme outlier + N = 5 + pixel_data = np.ones((N, N, N), dtype=float) + pixel_data[2, 2, 2] = 1e6 # extreme outlier inside the element bbox + grid_coords = GridCoords( + slice=np.arange(N), row=np.arange(N), col=np.arange(N) + ) + image_data = ImageData( + pixel_data, grid_coords, np.eye(3), np.zeros(3), PixelValueType.CT + ) + + element = np.array( + [ + [1, 1, 1], + [3, 1, 1], + [1, 3, 1], + [3, 3, 1], + [1, 1, 3], + [3, 1, 3], + [1, 3, 3], + [3, 3, 3], + ] + ) + + # Unfiltered + interp_no_filter = InterpolatorAllVoxel( + mode="allvoxels", filter_outliers=False + ) + val_no_filter = interp_no_filter._get_data_of_element(element, image_data) + + # Filtered + interp_filter = InterpolatorAllVoxel( + mode="allvoxels", filter_outliers=True + ) + val_filter = interp_filter._get_data_of_element(element, image_data) + + # Without filtering, mean should be much larger due to the outlier + assert val_no_filter > 10 # arbitrary threshold above normal mean of ones + # With filtering, mean should be close to 1 + assert np.isclose(val_filter, 1.0, rtol=1e-3) + + +def test_filter_outliers_applies_in_weighted_mode(): + """Outlier filtering should apply in weighted mode as well.""" + N = 5 + pixel_data = np.ones((N, N, N), dtype=float) + pixel_data[2, 2, 2] = 1e6 + grid_coords = GridCoords( + slice=np.arange(N), row=np.arange(N), col=np.arange(N) + ) + image_data = ImageData( + pixel_data, grid_coords, np.eye(3), np.zeros(3), PixelValueType.CT + ) + + element = np.array( + [ + [1, 1, 1], + [3, 1, 1], + [1, 3, 1], + [3, 3, 1], + [1, 1, 3], + [3, 1, 3], + [1, 3, 3], + [3, 3, 3], + ] + ) + + # Node weights biased toward corners near + # the outlier but filtering should mitigate + weights = np.array([5, 5, 5, 5, 5, 5, 5, 5], dtype=float) + + interp_w_no_filter = InterpolatorAllVoxel( + mode="allvoxels_weighted", filter_outliers=False + ) + val_w_no_filter = interp_w_no_filter._get_data_of_element( + element, image_data, node_weights_current=weights + ) + + interp_w_filter = InterpolatorAllVoxel( + mode="allvoxels_weighted", filter_outliers=True + ) + val_w_filter = interp_w_filter._get_data_of_element( + element, image_data, node_weights_current=weights + ) + + # Unfiltered should be larger than ~1 due to outlier influence + assert np.all(val_w_no_filter > 1.0) + # Filtered result should be close to baseline of ones + assert np.allclose(val_w_filter, 1.0, rtol=1e-3) + # And filtering should reduce the value relative to unfiltered + assert np.all(val_w_filter < val_w_no_filter) + + +def test_unknown_mode_defaults_to_unweighted_mean(): + """Unknown mode should fall back to unweighted mean.""" + element = np.array( + [ + [0, 2, 1], + [0, 3, 1], + [1, 3, 1], + [1, 2, 1], + [0, 2, 2], + [0, 3, 2], + [1, 3, 2], + [1, 2, 2], + ] + ) + N_slice, N_row, N_col = 5, 5, 5 + pixel_data = ( + np.arange(N_slice * N_row * N_col) + .reshape((N_slice, N_row, N_col)) + .astype(float) + ) + grid_coords = GridCoords(np.arange(5), np.arange(5), np.arange(5)) + image_data = ImageData( + pixel_data, grid_coords, np.eye(3), np.zeros(3), PixelValueType.CT + ) + + # Expected unweighted mean of voxels inside the element + data_inside = [11, 12, 16, 17, 36, 37, 41, 42] + expected_mean = np.mean(data_inside) + + # Use unknown mode to trigger fallback + interpolator = InterpolatorAllVoxel(mode="unknown", filter_outliers=False) + result = interpolator._get_data_of_element(element, image_data) + + # Compare scalar value robustly regardless of return shape + result_val = float(np.asarray(result).mean()) + assert np.isclose(result_val, expected_mean) + + +def test_weighted_small_voxel_count_uses_simple_weighted_average(): + """When voxel count is small (<=5), _weighted_voxel_mean should return + simple weighted average.""" + # Grid 2x2x2 + grid_coords = GridCoords( + slice=np.array([0.0, 1.0]), + row=np.array([0.0, 1.0]), + col=np.array([0.0, 1.0]), + ) + # Deterministic pixel values: s + 10*r + 100*c + s = grid_coords.slice[:, None, None] + r = grid_coords.row[None, :, None] + c = grid_coords.col[None, None, :] + pixel_data = (s + 10 * r + 100 * c).astype(float) + image_data = ImageData( + pixel_data, grid_coords, np.eye(3), np.zeros(3), PixelValueType.CT + ) + + # Element hull tightly around x,y ~ 0 and spanning z from 0 to 1 + element_nodes = np.array( + [ + [0.0, 0.0, 0.0], + [0.0, 0.0, 1.0], + [0.1, 0.0, 0.5], + [0.0, 0.1, 0.5], + ] + ) + # One element uses the 4 nodes, with explicit node weights + node_ids = np.arange(len(element_nodes)) + nodes = Nodes( + ids=node_ids, + coords=element_nodes, + weights=np.array([2.0, 2.0, 1.0, 1.0]), + ) + elements = [DisElement(node_ids=node_ids, id=0)] + surfaces = [Surface(node_ids=np.array([], dtype=int), id=0)] + dis = Discretization(nodes=nodes, elements=elements, surfaces=surfaces) + + interpolator = InterpolatorAllVoxel( + mode="allvoxels_weighted", filter_outliers=True + ) + elems = interpolator.compute_element_data(dis, image_data) + result = float(np.asarray(elems[0].data).mean()) + + # Manually compute included voxels and expected weighted average + hull = ConvexHull(element_nodes) + included_points = [] + included_vals = [] + for zi, zv in enumerate(grid_coords.slice): + for yi, yv in enumerate(grid_coords.row): + for xi, xv in enumerate(grid_coords.col): + p = np.array([zv, yv, xv]) + A, b = hull.equations[:, :-1], hull.equations[:, -1] + if np.all(A @ p + b <= 0): + included_points.append(p) + included_vals.append(pixel_data[zi, yi, xi]) + + voxels_phys = np.asarray(included_points) + values = np.asarray(included_vals) + + # Expect small voxel count (<=5) to use simple weighted average + assert len(values) <= 5 + + # Recompute voxel weights exactly as in implementation + distances = np.linalg.norm( + element_nodes[:, np.newaxis, :] - voxels_phys[np.newaxis, :, :], axis=2 + ) + distances = np.maximum(distances, 1e-10) + voxel_weights = np.sum(nodes.weights[:, np.newaxis] / distances, axis=0) + + expected = np.average(values, weights=voxel_weights, axis=0) + assert np.isclose(result, expected) diff --git a/tests/unittests/i2pp/core/interpolators/test_interpolator_types.py b/tests/unittests/i2pp/core/interpolators/test_interpolator_types.py index 11ac902..d83befa 100644 --- a/tests/unittests/i2pp/core/interpolators/test_interpolator_types.py +++ b/tests/unittests/i2pp/core/interpolators/test_interpolator_types.py @@ -11,25 +11,55 @@ "enum_value, expected_class", [ (InterpolationType.NODES, InterpolatorNodes), + (InterpolationType.NODES_WEIGHTED, InterpolatorNodes), (InterpolationType.CENTER, InterpolatorCenter), (InterpolationType.ALLVOXELS, InterpolatorAllVoxel), + (InterpolationType.ALLVOXELS_WEIGHTED, InterpolatorAllVoxel), ], ) -def test_get_interpolator(enum_value, expected_class): - """Test that get_interpolator returns the correct interpolator class.""" - assert enum_value.get_interpolator() == expected_class +def test_create_interpolator_returns_expected_class( + enum_value, expected_class +): + """Test create_interpolator returns instance of the expected class.""" + inst = enum_value.create_interpolator() + assert isinstance(inst, expected_class) def test_enum_values(): - """Test that enum values are correctly defined.""" - + """Enum values are correctly defined.""" assert InterpolationType.NODES.value == "nodes" + assert InterpolationType.NODES_WEIGHTED.value == "nodes_weighted" assert InterpolationType.CENTER.value == "elementcenter" assert InterpolationType.ALLVOXELS.value == "allvoxels" + assert InterpolationType.ALLVOXELS_WEIGHTED.value == "allvoxels_weighted" + + +@pytest.mark.parametrize( + "enum_value, expected_mode", + [ + (InterpolationType.ALLVOXELS, "allvoxels"), + (InterpolationType.ALLVOXELS_WEIGHTED, "allvoxels_weighted"), + ], +) +def test_create_interpolator_configures_mode(enum_value, expected_mode): + """InterpolatorAllVoxel is configured with correct mode.""" + inst = enum_value.create_interpolator() + assert isinstance(inst, InterpolatorAllVoxel) + assert getattr(inst, "_mode") == expected_mode + + +def test_create_interpolator_filter_outliers_flag(): + """filter_outliers flag is propagated to InterpolatorAllVoxel.""" + inst = InterpolationType.ALLVOXELS.create_interpolator( + filter_outliers=True + ) + assert isinstance(inst, InterpolatorAllVoxel) + assert getattr(inst, "_filter_outliers_enabled") is True -def test_get_interpolator_unsupported_member(): - """Test get_interpolator raises ValueError for unsupported enum members.""" +def test_create_interpolator_unsupported_member(): + """Test create_interpolator raises ValueError for unsupported enum + members.""" class FakeEnum: """Fake enum class to simulate an unsupported interpolation type.""" @@ -37,4 +67,5 @@ class FakeEnum: pass with pytest.raises(ValueError, match="Unsupported interpolation method"): - InterpolationType.get_interpolator(FakeEnum()) + # Directly call method on enum class to mimic bad usage + InterpolationType.create_interpolator(FakeEnum()) From d70600c6347053a9e4cda9151f3759c1abdcbb1d Mon Sep 17 00:00:00 2001 From: Saskia Gloy Date: Tue, 30 Dec 2025 12:22:57 +0100 Subject: [PATCH 05/15] Add filter_outliers in Processing configuration --- .../core/configuration_validator/validator.py | 21 +++++++++++++++++-- src/i2pp/core/run.py | 1 + 2 files changed, 20 insertions(+), 2 deletions(-) diff --git a/src/i2pp/core/configuration_validator/validator.py b/src/i2pp/core/configuration_validator/validator.py index 690808e..a413e54 100644 --- a/src/i2pp/core/configuration_validator/validator.py +++ b/src/i2pp/core/configuration_validator/validator.py @@ -108,18 +108,35 @@ class Processing: interpolation_method: str smoothing: Optional[Smoothing] transformation: Transformation + filter_outliers: Optional[bool] = field(default=False) @staticmethod def from_dict(d: Dict[str, Any]) -> "Processing": """Creates a Processing instance from a dictionary.""" + # support both top-level 'interpolation_method' + # and nested 'interpolation.method' + interp_method = d.get("interpolation_method") + if interp_method is None and "interpolation" in d: + interp_method = d["interpolation"].get("method") + if not isinstance(interp_method, str) or not interp_method: + raise KeyError("interpolation_method") + + # support filter_outliers at top-level or nested + filter_out = d.get("filter_outliers") + if filter_out is None: + filter_out = d.get("interpolation", {}).get( + "filter_outliers", False + ) + return Processing( - interpolation_method=d["interpolation_method"], + interpolation_method=interp_method, smoothing=( - Smoothing.from_dict(d["smoothing"]) + Smoothing.from_dict(d.get("smoothing")) if d.get("smoothing") is not None else None ), transformation=Transformation.from_dict(d["transformation"]), + filter_outliers=filter_out, ) diff --git a/src/i2pp/core/run.py b/src/i2pp/core/run.py index d126f20..703424b 100644 --- a/src/i2pp/core/run.py +++ b/src/i2pp/core/run.py @@ -77,6 +77,7 @@ def run_i2pp(config_i2pp): discretization, image, interpolation_method=config.processing.interpolation_method, + filter_outliers=config.processing.filter_outliers, ) # Transform the data using the user-defined python function From 9ddd9e380c02e7fde7400c5c3519a3c05192a524 Mon Sep 17 00:00:00 2001 From: Saskia Gloy Date: Tue, 30 Dec 2025 13:00:46 +0100 Subject: [PATCH 06/15] Add weighted and unweighted mean computations to InterpolatorNodes --- .../core/interpolators/interpolator_nodes.py | 133 +++++++++++++++++- .../core/interpolators/interpolator_types.py | 27 +++- .../interpolators/test_interpolator_nodes.py | 104 ++++++++++++++ 3 files changed, 252 insertions(+), 12 deletions(-) diff --git a/src/i2pp/core/interpolators/interpolator_nodes.py b/src/i2pp/core/interpolators/interpolator_nodes.py index 584ad39..fdced2e 100644 --- a/src/i2pp/core/interpolators/interpolator_nodes.py +++ b/src/i2pp/core/interpolators/interpolator_nodes.py @@ -21,6 +21,114 @@ class InterpolatorNodes(Interpolator): "nodes". """ + # Add mode to control whether to use node weights + def __init__(self, *args, mode: str = "nodes", **kwargs): + """Initializes the InterpolatorNodes.""" + super().__init__() + self._mode = mode # "nodes" or "nodes_weighted" + + # Helpers for readability and error handling + def _compute_unweighted_mean( + self, ele_nodes: np.ndarray, num_values: int + ) -> np.ndarray: + """Computes the unweighted mean of pixel values for an element's nodes. + + This method calculates the arithmetic mean of pixel values associated + with an element's nodes. It handles both single-channel and multi- + channel image data. A key feature is its handling of NaN values: if a + node has a NaN value in any channel, that node is excluded from the + mean calculation. If all nodes for an element have NaN values, the + element's resulting value will be NaN. + + Args: + ele_nodes (np.ndarray): An array of interpolated pixel values for + the nodes of a single element. Shape can be (num_nodes,) for + single-channel data or (num_nodes, num_values) for multi- + channel data. + num_values (int): The number of values per pixel (e.g., 1 for + grayscale, 3 for RGB). + + Returns: + np.ndarray: An array containing the mean pixel value(s) for the + element. The shape is (num_values,). + """ + # Mask out nodes that contain any NaN across channels + if ele_nodes.ndim == 1: + nan_mask = ~np.isnan(ele_nodes) + else: + nan_mask = ~np.isnan(ele_nodes).any(axis=1) + + if not np.any(nan_mask): + return np.full(num_values, np.nan) + + valid = ele_nodes[nan_mask] + mean_val = np.mean(valid, axis=0) + return mean_val if num_values > 1 else np.array([mean_val]) + + def _compute_weighted_mean( + self, + ele_nodes: np.ndarray, + weights_current: np.ndarray, + num_values: int, + ) -> np.ndarray: + """Computes the weighted mean of pixel values for an element's nodes. + + This method calculates the weighted average of pixel values from an + element's nodes. It is used when a more nuanced contribution of each + node is desired, based on pre-assigned weights. Nodes with NaN values + are excluded from the calculation. The weights of the valid nodes are + normalized to sum to 1 before computing the average to prevent scaling + issues. + + Args: + ele_nodes (np.ndarray): An array of interpolated pixel values for + the nodes of a single element. Shape can be (num_nodes,) for + single-channel data or (num_nodes, num_values) for multi- + channel data. + weights_current (np.ndarray): An array of weights corresponding to + each node in `ele_nodes`. Shape is (num_nodes,). + num_values (int): The number of values per pixel (e.g., 1 for + grayscale, 3 for RGB). + + Returns: + np.ndarray: An array containing the weighted mean pixel value(s) + for the element. The shape is (num_values,). + + Raises: + ValueError: If `weights_current` is None or if the shape of + `ele_nodes` and `weights_current` are incompatible. + """ + if weights_current is None: + raise ValueError("Node weights are required for weighted mode.") + if ele_nodes.shape[0] != weights_current.shape[0]: + raise ValueError( + "Incompatible shapes:" + f" ele_nodes has {ele_nodes.shape[0]} nodes, " + f"weights_current has {weights_current.shape[0]}." + ) + + # Mask out nodes that contain any NaN across channels + if ele_nodes.ndim == 1: + nan_mask = ~np.isnan(ele_nodes) + else: + nan_mask = ~np.isnan(ele_nodes).any(axis=1) + + if not np.any(nan_mask): + return np.full(num_values, np.nan) + + vals = ele_nodes[nan_mask] + w = weights_current[nan_mask] + + wsum = float(np.sum(w)) + if wsum > 0: + w = w / wsum # normalize to avoid scaling issues + + if num_values == 1: + val = np.average(vals.reshape(-1), weights=w) + return np.array([val]) + else: + return np.average(vals, weights=w, axis=0) + def compute_element_data( self, dis: Discretization, image_data: ImageData ) -> list[Element]: @@ -59,18 +167,31 @@ def compute_element_data( ] ) + # Only prepare node weights if we're in weighted mode + node_weights = None + if ( + self._mode == "nodes_weighted" + and getattr(dis.nodes, "weights", None) is not None + ): + from typing import cast + + node_weights = cast(np.ndarray, dis.nodes.weights) + for i, ele in tqdm( enumerate(dis.elements), total=len(dis.elements), desc="Processing Elements", ): ele_nodes = node_values[node_positions[i]] - - ele.data = ( - np.nanmean(ele_nodes, axis=0) - if not np.all(np.isnan(ele_nodes)) - else np.full(image_data.pixel_type.num_values, np.nan) - ) + num_values = image_data.pixel_type.num_values + + if node_weights is not None: + weights_current = node_weights[node_positions[i]] + ele.data = self._compute_weighted_mean( + ele_nodes, weights_current, num_values + ) + else: + ele.data = self._compute_unweighted_mean(ele_nodes, num_values) if np.all(np.isnan(ele.data)): self.nan_elements += 1 diff --git a/src/i2pp/core/interpolators/interpolator_types.py b/src/i2pp/core/interpolators/interpolator_types.py index e7400de..42dd024 100644 --- a/src/i2pp/core/interpolators/interpolator_types.py +++ b/src/i2pp/core/interpolators/interpolator_types.py @@ -41,11 +41,24 @@ class InterpolationType(Enum): def create_interpolator( self, *, filter_outliers: bool = False ) -> Interpolator: - """Factory that returns a configured interpolator instance. + """Creates and returns a configured interpolator instance based on the + selected interpolation method. - For ALLVOXELS/ALLVOXELS_WEIGHTED, this method configures the - unified InterpolatorAllVoxel with the appropriate mode and - optional outlier filtering. + This method returns a configured interpolator instance for assigning + pixel values to FEM elements, depending on the current interpolation + method. Supported methods include interpolation at element nodes, + element centers, or averaging all voxels within an element. + + Args: + filter_outliers (bool): If True, outliers will be filtered during + interpolation. Defaults to False. + + Returns: + Interpolator: An instance of the interpolator that matches the + specified interpolation method. + + Raises: + ValueError: If the interpolation method is not supported. """ if self == InterpolationType.ALLVOXELS: return InterpolatorAllVoxel( @@ -55,8 +68,10 @@ def create_interpolator( return InterpolatorAllVoxel( mode="allvoxels_weighted", filter_outliers=filter_outliers ) - if self in (InterpolationType.NODES, InterpolationType.NODES_WEIGHTED): - return InterpolatorNodes() + if self == InterpolationType.NODES: + return InterpolatorNodes(mode="nodes") + if self == InterpolationType.NODES_WEIGHTED: + return InterpolatorNodes(mode="nodes_weighted") if self == InterpolationType.CENTER: return InterpolatorCenter() raise ValueError(f"Unsupported interpolation method: {self}") diff --git a/tests/unittests/i2pp/core/interpolators/test_interpolator_nodes.py b/tests/unittests/i2pp/core/interpolators/test_interpolator_nodes.py index 6fb6d45..87cbace 100644 --- a/tests/unittests/i2pp/core/interpolators/test_interpolator_nodes.py +++ b/tests/unittests/i2pp/core/interpolators/test_interpolator_nodes.py @@ -80,3 +80,107 @@ def test_compute_element_data(): return_grid_coords, image_data ) mock_get_node_position_of_element.call_count == 2 + + +# Helper +def _build_simple_setup_rgb(): + """Build simple discretization and image data for RGB tests.""" + node_coords = np.array([[0, 0, 0], [1, 0, 0]]) + nodes = Nodes(coords=node_coords, ids=[0, 1]) + # Provide weights for nodes + nodes.weights = np.array([0.2, 0.8]) + + elements = [Element(node_ids=[0, 1], id=0)] + dis = Discretization(nodes=nodes, elements=elements, surfaces=[]) + + pixel_data = np.zeros((2, 2, 2, 3)) + grid_coords = GridCoords(np.arange(2), np.arange(2), np.arange(2)) + image_data = ImageData( + pixel_data=pixel_data, + grid_coords=grid_coords, + orientation=np.eye(3), + position=np.array([0, 0, 0]), + pixel_type=PixelValueType.RGB, + ) + return dis, image_data + + +def test_nodes_unweighted_rgb_ignores_weights(): + """Unweighted mode computes nanmean and ignores weights for RGB.""" + dis, image_data = _build_simple_setup_rgb() + interpolator = InterpolatorNodes(mode="nodes") + + node_values = np.array([[100, 150, 200], [200, 250, 300]]) + with ( + patch.object( + InterpolatorNodes, + "world_to_grid_coords", + return_value=np.array([[0.5, 0.5, 0.5], [1.5, 1.5, 1.5]]), + ), + patch.object( + InterpolatorNodes, + "interpolate_image_values_to_points", + return_value=node_values, + ), + patch( + "i2pp.core.utilities.get_node_position_of_element", + return_value=np.array([0, 1]), + ), + ): + result = interpolator.compute_element_data(dis, image_data) + expected = np.mean(node_values, axis=0) + assert np.allclose(result[0].data, expected) + + +def test_nodes_weighted_rgb(): + """Weighted mode uses node weights for vector-valued (RGB) data.""" + dis, image_data = _build_simple_setup_rgb() + interpolator = InterpolatorNodes(mode="nodes_weighted") + + node_values = np.array([[100, 150, 200], [200, 250, 300]]) + with ( + patch.object( + InterpolatorNodes, + "world_to_grid_coords", + return_value=np.array([[0.5, 0.5, 0.5], [1.5, 1.5, 1.5]]), + ), + patch.object( + InterpolatorNodes, + "interpolate_image_values_to_points", + return_value=node_values, + ), + patch( + "i2pp.core.utilities.get_node_position_of_element", + return_value=np.array([0, 1]), + ), + ): + result = interpolator.compute_element_data(dis, image_data) + expected = 0.2 * node_values[0] + 0.8 * node_values[1] + assert np.allclose(result[0].data, expected) + + +def test_nodes_weighted_rgb_handles_nans(): + """Weighted mode handles NaNs by masking them out for RGB.""" + dis, image_data = _build_simple_setup_rgb() + interpolator = InterpolatorNodes(mode="nodes_weighted") + + node_values = np.array([[np.nan, np.nan, np.nan], [200, 250, 300]]) + with ( + patch.object( + InterpolatorNodes, + "world_to_grid_coords", + return_value=np.array([[0.5, 0.5, 0.5], [1.5, 1.5, 1.5]]), + ), + patch.object( + InterpolatorNodes, + "interpolate_image_values_to_points", + return_value=node_values, + ), + patch( + "i2pp.core.utilities.get_node_position_of_element", + return_value=np.array([0, 1]), + ), + ): + result = interpolator.compute_element_data(dis, image_data) + # Only the non-NaN node should contribute + assert np.allclose(result[0].data, np.array([200, 250, 300])) From 87b0da96412f9c3a9a780e5cfcb522ff743d403b Mon Sep 17 00:00:00 2001 From: Saskia Gloy Date: Tue, 30 Dec 2025 14:47:36 +0100 Subject: [PATCH 07/15] enable setting weights in config file --- .../core/configuration_validator/validator.py | 59 +++++++++++-------- src/i2pp/core/discretization_helpers.py | 9 ++- .../discretization_reader.py | 2 +- .../fourc_yaml_reader.py | 14 +++-- .../core/interpolators/interpolator_nodes.py | 13 ++-- src/i2pp/core/run.py | 7 ++- .../test_data/random_colors.yaml.jinja | 5 +- tests/integration/test_random_colors.py | 1 + .../configuration_validator/test_validator.py | 13 +++- .../test_fourc_yaml_reader.py | 12 +++- .../i2pp/core/test_discretization_helpers.py | 5 +- tests/unittests/i2pp/core/test_run.py | 4 +- 12 files changed, 92 insertions(+), 52 deletions(-) diff --git a/src/i2pp/core/configuration_validator/validator.py b/src/i2pp/core/configuration_validator/validator.py index a413e54..2d7f4ca 100644 --- a/src/i2pp/core/configuration_validator/validator.py +++ b/src/i2pp/core/configuration_validator/validator.py @@ -67,7 +67,7 @@ def from_dict(d: Dict[str, Any]) -> "Import": class Smoothing: """Class representing the smoothing configuration.""" - smoothing_area: int = 3 + area: int = 3 visualize: bool = False @staticmethod @@ -76,7 +76,7 @@ def from_dict(d: Optional[Dict[str, Any]]) -> Optional["Smoothing"]: if d is None: return None return Smoothing( - smoothing_area=d.get("smoothing_area", 3), + area=d.get("area", 3), visualize=d.get("visualize", False), ) @@ -101,42 +101,49 @@ def from_dict(d: Dict[str, Any]) -> "Transformation": ) +@dataclass +class Interpolation: + """Class representing the interpolation configuration.""" + + method: str + interior_node_weight: float = 1.0 + surface_node_weight: float = 0.0 + filter_outliers: bool = False + + @staticmethod + def from_dict(d: Dict[str, Any]) -> "Interpolation": + """Creates an Interpolation instance from a dictionary.""" + method = d.get("method") + if not isinstance(method, str) or not method: + raise KeyError("interpolation.method") + + node_weight_cfg = d.get("node_weight", {}) + interior_weight = float(node_weight_cfg.get("interior", 1.0)) + surface_weight = float(node_weight_cfg.get("surface", 0.0)) + + return Interpolation( + method=method, + interior_node_weight=interior_weight, + surface_node_weight=surface_weight, + filter_outliers=d.get("filter_outliers", False), + ) + + @dataclass class Processing: """Class representing the processing configuration.""" - interpolation_method: str + interpolation: Interpolation smoothing: Optional[Smoothing] transformation: Transformation - filter_outliers: Optional[bool] = field(default=False) @staticmethod def from_dict(d: Dict[str, Any]) -> "Processing": """Creates a Processing instance from a dictionary.""" - # support both top-level 'interpolation_method' - # and nested 'interpolation.method' - interp_method = d.get("interpolation_method") - if interp_method is None and "interpolation" in d: - interp_method = d["interpolation"].get("method") - if not isinstance(interp_method, str) or not interp_method: - raise KeyError("interpolation_method") - - # support filter_outliers at top-level or nested - filter_out = d.get("filter_outliers") - if filter_out is None: - filter_out = d.get("interpolation", {}).get( - "filter_outliers", False - ) - return Processing( - interpolation_method=interp_method, - smoothing=( - Smoothing.from_dict(d.get("smoothing")) - if d.get("smoothing") is not None - else None - ), + interpolation=Interpolation.from_dict(d["interpolation"]), + smoothing=Smoothing.from_dict(d.get("smoothing")), transformation=Transformation.from_dict(d["transformation"]), - filter_outliers=filter_out, ) diff --git a/src/i2pp/core/discretization_helpers.py b/src/i2pp/core/discretization_helpers.py index 44e3157..a4f0649 100644 --- a/src/i2pp/core/discretization_helpers.py +++ b/src/i2pp/core/discretization_helpers.py @@ -50,7 +50,9 @@ def determine_discretization_format(file_path: Path) -> DiscretizationFormat: def verify_and_load_discretization( - discretization_path: Path, options: dict + discretization_path: Path, + options: dict, + processing: dict, ) -> Discretization: """Loads and processes mesh data. @@ -62,6 +64,7 @@ def verify_and_load_discretization( discretization_path (Path): Path to the discretization file. options (dict): Options for loading the discretization that are passed to the reader classes. + processing (dict): Processing options are passed to the reader classes. Returns: DiscretizationData: The loaded and processed mesh data. @@ -70,7 +73,9 @@ def verify_and_load_discretization( dis_reader = dis_format.get_reader()() - dis = dis_reader.load_discretization(discretization_path, options) + dis = dis_reader.load_discretization( + discretization_path, options, processing + ) bounding_box = find_mins_maxs(points=dis.nodes.coords, enlargement=2) diff --git a/src/i2pp/core/discretization_readers/discretization_reader.py b/src/i2pp/core/discretization_readers/discretization_reader.py index 6a5a045..3cc8bff 100644 --- a/src/i2pp/core/discretization_readers/discretization_reader.py +++ b/src/i2pp/core/discretization_readers/discretization_reader.py @@ -126,7 +126,7 @@ def __init__(self): @abstractmethod def load_discretization( - self, file_path: Path, options: dict + self, file_path: Path, options: dict, processing: dict ) -> Discretization: """Abstract method to load discretization data from a file path. diff --git a/src/i2pp/core/discretization_readers/fourc_yaml_reader.py b/src/i2pp/core/discretization_readers/fourc_yaml_reader.py index a2a2160..5c25ad9 100644 --- a/src/i2pp/core/discretization_readers/fourc_yaml_reader.py +++ b/src/i2pp/core/discretization_readers/fourc_yaml_reader.py @@ -71,7 +71,10 @@ def _filter_discretization( return dis def load_discretization( - self, file_path: Path, options: dict + self, + file_path: Path, + options: dict, + processing: "I2PPConfig.processing", ) -> Discretization: """Loads and processes a finite element discretization from a .4C.yaml file. @@ -85,6 +88,7 @@ def load_discretization( options (dict): Options for loading the discretization. Filtering for material ids can be enabled by specifying `material_ids` in the options dictionary. + processing (I2PPConfig.processing): Processing configuration object. Returns: Discretization: The finite element discretization including nodes @@ -102,6 +106,9 @@ def load_discretization( raw_dis, np.array(options["material_ids"]) ) + interior_weight = processing.interpolation.interior_node_weight + surface_weight = processing.interpolation.surface_node_weight + nodes_coords = [] node_ids = [] nodes_weights = [] @@ -109,7 +116,7 @@ def load_discretization( for node in raw_dis.nodes: nodes_coords.append(node.coords) node_ids.append(node.id) - nodes_weights.append(5.0) + nodes_weights.append(interior_weight) elements = [] @@ -128,9 +135,8 @@ def load_discretization( surf_node_ids = [] for node in surf.nodes: surf_node_ids.append(node.id) - # get position of node in node_ids list and set weight position = node_ids.index(node.id) - nodes_weights[position] = 0.0 + nodes_weights[position] = surface_weight surfaces.append( Surface(node_ids=np.array(surf_node_ids), id=surf.id) diff --git a/src/i2pp/core/interpolators/interpolator_nodes.py b/src/i2pp/core/interpolators/interpolator_nodes.py index fdced2e..7b49833 100644 --- a/src/i2pp/core/interpolators/interpolator_nodes.py +++ b/src/i2pp/core/interpolators/interpolator_nodes.py @@ -1,5 +1,7 @@ """Interpolates pixel values from image-data to mesh-data.""" +from typing import cast + import numpy as np from i2pp.core.discretization_readers.discretization_reader import ( Discretization, @@ -120,8 +122,13 @@ def _compute_weighted_mean( w = weights_current[nan_mask] wsum = float(np.sum(w)) - if wsum > 0: - w = w / wsum # normalize to avoid scaling issues + if wsum == 0: + if num_values > 1: + return np.mean(vals, axis=0) + else: + return np.array([np.mean(vals)]) + + w = w / wsum # Normalize weights if num_values == 1: val = np.average(vals.reshape(-1), weights=w) @@ -173,8 +180,6 @@ def compute_element_data( self._mode == "nodes_weighted" and getattr(dis.nodes, "weights", None) is not None ): - from typing import cast - node_weights = cast(np.ndarray, dis.nodes.weights) for i, ele in tqdm( diff --git a/src/i2pp/core/run.py b/src/i2pp/core/run.py index 703424b..5ebc5fe 100644 --- a/src/i2pp/core/run.py +++ b/src/i2pp/core/run.py @@ -45,6 +45,7 @@ def run_i2pp(config_i2pp): discretization = verify_and_load_discretization( config.import_.discretization.path, config.import_.discretization.options, + config.processing, ) # Load the image data @@ -65,7 +66,7 @@ def run_i2pp(config_i2pp): image.pixel_data = smooth_data( image.pixel_data, - config.processing.smoothing.smoothing_area, + config.processing.smoothing.area, mask=image.mask, ) @@ -76,8 +77,8 @@ def run_i2pp(config_i2pp): elements = interpolate_image_to_discretization( discretization, image, - interpolation_method=config.processing.interpolation_method, - filter_outliers=config.processing.filter_outliers, + interpolation_method=config.processing.interpolation.method, + filter_outliers=config.processing.interpolation.filter_outliers, ) # Transform the data using the user-defined python function diff --git a/tests/integration/test_data/random_colors.yaml.jinja b/tests/integration/test_data/random_colors.yaml.jinja index b01d542..7f9bba3 100644 --- a/tests/integration/test_data/random_colors.yaml.jinja +++ b/tests/integration/test_data/random_colors.yaml.jinja @@ -14,9 +14,10 @@ import: image_position: [-0.25, -0.1, -0.5] processing: smoothing: - smoothing_area: 2 + area: 2 visualize: false - interpolation_method: nodes + interpolation: + method: nodes transformation: user_script: {{ user_function }} user_function: user_function diff --git a/tests/integration/test_random_colors.py b/tests/integration/test_random_colors.py index 42c51a1..111bdfd 100644 --- a/tests/integration/test_random_colors.py +++ b/tests/integration/test_random_colors.py @@ -42,6 +42,7 @@ def test_i2pp_integration_random_colors(): image_folder=image_folder, user_function=user_function, ) + input_file_filled = os.path.join(tmpdir, "random_colors.yaml") with open(input_file_filled, "w") as file: file.write(rendered_content) diff --git a/tests/unittests/i2pp/core/configuration_validator/test_validator.py b/tests/unittests/i2pp/core/configuration_validator/test_validator.py index 07136c8..a5f2526 100644 --- a/tests/unittests/i2pp/core/configuration_validator/test_validator.py +++ b/tests/unittests/i2pp/core/configuration_validator/test_validator.py @@ -28,7 +28,9 @@ def minimal_valid_config(tmp_path): "image": {"path": str(image_path), "type": "dicom"}, }, "processing": { - "interpolation_method": "nodes", + "interpolation":{ + "method": "nodes", + }, "transformation": { "user_script": str(script_path), "user_function": "process_image_data", @@ -86,8 +88,13 @@ def large_valid_config(tmp_path): }, }, "processing": { - "smoothing": {"smoothing_area": 5, "visualize": True}, - "interpolation_method": "elements", + "smoothing": { + "area": 5, + "visualize": True + }, + "interpolation": { + "method": "elements", + }, "transformation": { "user_script": str(script_path), "user_function": "process_image_data", diff --git a/tests/unittests/i2pp/core/discretization_readers/test_fourc_yaml_reader.py b/tests/unittests/i2pp/core/discretization_readers/test_fourc_yaml_reader.py index 7a8b7a4..f3852c8 100644 --- a/tests/unittests/i2pp/core/discretization_readers/test_fourc_yaml_reader.py +++ b/tests/unittests/i2pp/core/discretization_readers/test_fourc_yaml_reader.py @@ -84,9 +84,15 @@ def test_load_discretization_fourc_yaml_without_filter(tmp_path: Path) -> None: mock_dis.elements.structure = [ele1, ele2] mock_dis.nodes = [node1, node2, node3, node4] - test_config = {"material_ids": None} - - dis_loaded = test_dis.load_discretization(Path(test_path), test_config) + # test_config + test_options = {"material_ids": None} + mock_processing = MagicMock() + mock_processing.interpolation.interior_node_weight = 1.0 + mock_processing.interpolation.surface_node_weight = 0.0 + + dis_loaded = test_dis.load_discretization( + Path(test_path), test_options, mock_processing + ) assert np.array_equal( dis_loaded.nodes.coords, diff --git a/tests/unittests/i2pp/core/test_discretization_helpers.py b/tests/unittests/i2pp/core/test_discretization_helpers.py index 632d19d..da15423 100644 --- a/tests/unittests/i2pp/core/test_discretization_helpers.py +++ b/tests/unittests/i2pp/core/test_discretization_helpers.py @@ -73,6 +73,7 @@ def test_verify_and_load_discretization(): absolut_path = Path.cwd() / "test_path.4C.yaml" options_dict = {} + processing_dict = {} mock_bounding = tuple([[0, 0, 0], [1, 1, 1]]) nodes = Nodes([0, 0, 0], 0) mock_dis = Discretization(nodes, [], []) @@ -89,14 +90,14 @@ def test_verify_and_load_discretization(): return_value=mock_bounding, ) as mock_find_mins_maxs: dis = verify_and_load_discretization( - absolut_path, options_dict + absolut_path, options_dict, processing_dict ) mock_determine_discretization_format.assert_called_once_with( absolut_path ) mock_load_discretization.assert_called_once_with( - absolut_path, options_dict + absolut_path, options_dict, processing_dict ) mock_find_mins_maxs.assert_called_once_with( points=[0, 0, 0], enlargement=2 diff --git a/tests/unittests/i2pp/core/test_run.py b/tests/unittests/i2pp/core/test_run.py index ee49b11..519f8a0 100644 --- a/tests/unittests/i2pp/core/test_run.py +++ b/tests/unittests/i2pp/core/test_run.py @@ -38,7 +38,7 @@ def minimal_valid_config(tmp_path): }, }, "processing": { - "interpolation_method": "nodes", + "interpolation": {"method": "nodes"}, "transformation": { "user_script": str(script_path), "user_function": "process_image_data", @@ -147,7 +147,7 @@ def maximal_valid_config(tmp_path): }, }, "processing": { - "interpolation_method": "nodes", + "interpolation": {"method": "nodes"}, "transformation": { "user_script": str(script_path), "user_function": "process_image_data", From bc630d1f2d9e639a9c0fa4b85c184eba9967c8bc Mon Sep 17 00:00:00 2001 From: Saskia Gloy Date: Sun, 4 Jan 2026 14:59:43 +0100 Subject: [PATCH 08/15] Refactor interpolation configuration, enable setting surface node value prior node intepolation --- .../core/configuration_validator/validator.py | 29 ++++---- src/i2pp/core/interpolate_element_data.py | 18 ++--- .../interpolators/interpolator_all_voxel.py | 5 +- .../core/interpolators/interpolator_nodes.py | 66 ++++++++++++++++++- .../core/interpolators/interpolator_types.py | 15 ++++- src/i2pp/core/run.py | 3 +- .../configuration_validator/test_validator.py | 9 +-- tests/unittests/i2pp/core/test_run.py | 11 +++- 8 files changed, 108 insertions(+), 48 deletions(-) diff --git a/src/i2pp/core/configuration_validator/validator.py b/src/i2pp/core/configuration_validator/validator.py index 2d7f4ca..6c44eff 100644 --- a/src/i2pp/core/configuration_validator/validator.py +++ b/src/i2pp/core/configuration_validator/validator.py @@ -106,26 +106,16 @@ class Interpolation: """Class representing the interpolation configuration.""" method: str - interior_node_weight: float = 1.0 - surface_node_weight: float = 0.0 - filter_outliers: bool = False + filter_outliers: Optional[bool] = field(default=False) + set_node_value: Optional[float] = field(default=None) @staticmethod def from_dict(d: Dict[str, Any]) -> "Interpolation": """Creates an Interpolation instance from a dictionary.""" - method = d.get("method") - if not isinstance(method, str) or not method: - raise KeyError("interpolation.method") - - node_weight_cfg = d.get("node_weight", {}) - interior_weight = float(node_weight_cfg.get("interior", 1.0)) - surface_weight = float(node_weight_cfg.get("surface", 0.0)) - return Interpolation( - method=method, - interior_node_weight=interior_weight, - surface_node_weight=surface_weight, + method=d["method"], filter_outliers=d.get("filter_outliers", False), + set_node_value=d.get("set_surface_node_value", None), ) @@ -133,17 +123,22 @@ def from_dict(d: Dict[str, Any]) -> "Interpolation": class Processing: """Class representing the processing configuration.""" - interpolation: Interpolation smoothing: Optional[Smoothing] transformation: Transformation + interpolation: Interpolation @staticmethod def from_dict(d: Dict[str, Any]) -> "Processing": """Creates a Processing instance from a dictionary.""" + return Processing( - interpolation=Interpolation.from_dict(d["interpolation"]), - smoothing=Smoothing.from_dict(d.get("smoothing")), + smoothing=( + Smoothing.from_dict(d.get("smoothing")) + if d.get("smoothing") is not None + else None + ), transformation=Transformation.from_dict(d["transformation"]), + interpolation=Interpolation.from_dict(d["interpolation"]), ) diff --git a/src/i2pp/core/interpolate_element_data.py b/src/i2pp/core/interpolate_element_data.py index f27c5ce..a2cd00a 100644 --- a/src/i2pp/core/interpolate_element_data.py +++ b/src/i2pp/core/interpolate_element_data.py @@ -1,5 +1,6 @@ """Interpolate image data to FEM-Elements.""" +from i2pp.core.configuration_validator.validator import Interpolation from i2pp.core.discretization_readers.discretization_reader import ( Discretization, Element, @@ -11,9 +12,7 @@ def interpolate_image_to_discretization( dis: Discretization, image_data: ImageData, - interpolation_method: str, - *, - filter_outliers: bool = False, + interpolation: Interpolation, ) -> list[Element]: """Performs interpolation of image data onto the FEM Discretization based on the specified interpolation method. @@ -34,21 +33,16 @@ def interpolate_image_to_discretization( surfaces, elements and node coordinates. image_data (ImageData): A structured representation containing 3D pixel data, grid coordinates, orientation, and metadata. - interpolation_method (str): The type of interpolation to perform for - assigning pixel values to the elements. This should match one of - the `InterpolationType` enum values (e.g., "nodes", - "elementcenter", "allvoxels"). - filter_outliers (bool, optional): Whether to filter outliers during - interpolation. Defaults to False. + interpolation (Interpolation): The interpolation configuration object. Returns: list[Element]: A list of FEM elements with interpolated pixel data. """ - enum_interpolation_method = InterpolationType(interpolation_method) + enum_interpolation_method = InterpolationType(interpolation.method) - # Use factory to configure weighted vs unweighted and outlier filtering interpolator = enum_interpolation_method.create_interpolator( - filter_outliers=filter_outliers + filter_outliers=interpolation.filter_outliers, + set_surf_node_value=interpolation.set_surface_node_value, ) return interpolator.compute_element_data(dis, image_data) diff --git a/src/i2pp/core/interpolators/interpolator_all_voxel.py b/src/i2pp/core/interpolators/interpolator_all_voxel.py index 2c7f74a..ea7c050 100644 --- a/src/i2pp/core/interpolators/interpolator_all_voxel.py +++ b/src/i2pp/core/interpolators/interpolator_all_voxel.py @@ -1,6 +1,6 @@ """Interpolates pixel values from image-data to mesh-data.""" -from typing import Tuple +from typing import Tuple, cast import numpy as np from i2pp.core.discretization_readers.discretization_reader import ( @@ -376,9 +376,6 @@ def compute_element_data( # Cast weights to ndarray when present to satisfy type checker node_weights_current = None if getattr(dis.nodes, "weights", None) is not None: - # mypy: explicitly cast Optional[Any] to np.ndarray - from typing import cast - weights_nd = cast(np.ndarray, dis.nodes.weights) node_weights_current = weights_nd[node_positions[i]] diff --git a/src/i2pp/core/interpolators/interpolator_nodes.py b/src/i2pp/core/interpolators/interpolator_nodes.py index 7b49833..a891386 100644 --- a/src/i2pp/core/interpolators/interpolator_nodes.py +++ b/src/i2pp/core/interpolators/interpolator_nodes.py @@ -1,6 +1,6 @@ """Interpolates pixel values from image-data to mesh-data.""" -from typing import cast +from typing import Optional, Union, cast import numpy as np from i2pp.core.discretization_readers.discretization_reader import ( @@ -24,10 +24,17 @@ class InterpolatorNodes(Interpolator): """ # Add mode to control whether to use node weights - def __init__(self, *args, mode: str = "nodes", **kwargs): + def __init__( + self, + *args, + mode: str = "nodes", + surf_val: Optional[Union[float, np.ndarray]] = None, + **kwargs, + ): """Initializes the InterpolatorNodes.""" super().__init__() self._mode = mode # "nodes" or "nodes_weighted" + self.set_surf_node_value = surf_val # Helpers for readability and error handling def _compute_unweighted_mean( @@ -136,6 +143,53 @@ def _compute_weighted_mean( else: return np.average(vals, weights=w, axis=0) + def _override_surface_nodes( + self, + *, + node_values: np.ndarray, + dis: Discretization, + surf_node_value: Union[np.ndarray, float], + num_values: int, + ) -> None: + """Overrides the pixel values at surface nodes with a specified value. + + This method modifies the `node_values` array in place, setting the + values of nodes that belong to any surface in the discretization to + the provided `surf_node_value`. + + Arguments: + node_values (np.ndarray): Array of pixel values at each node. + dis (Discretization): The FEM discretization containing surfaces + and node data. + surf_node_value (np.ndarray | float): The value to assign to + surface nodes. Can be a single float or an array matching + the number of pixel values. + num_values (int): The number of pixel values per node. + """ + + surf_val = np.asarray(surf_node_value) + + if surf_val.size != num_values: + raise ValueError( + f"set_surf_node_value must have {num_values} value(s), " + f"got {surf_val.size}" + ) + + surface_node_ids = { + nid for surf in dis.surfaces for nid in surf.node_ids + } + + node_id_to_index = {nid: i for i, nid in enumerate(dis.nodes.ids)} + + surface_indices = [ + node_id_to_index[nid] + for nid in surface_node_ids + if nid in node_id_to_index + ] + + if surface_indices: + node_values[surface_indices] = surf_node_value + def compute_element_data( self, dis: Discretization, image_data: ImageData ) -> list[Element]: @@ -167,6 +221,14 @@ def compute_element_data( node_grid_coords, image_data ) + if self.set_surf_node_value is not None and dis.surfaces: + self._override_surface_nodes( + node_values=node_values, + dis=dis, + surf_node_value=self.set_surf_node_value, + num_values=image_data.pixel_type.num_values, + ) + node_positions = np.array( [ get_node_position_of_element(ele.node_ids, dis.nodes.ids) diff --git a/src/i2pp/core/interpolators/interpolator_types.py b/src/i2pp/core/interpolators/interpolator_types.py index 42dd024..0045ea6 100644 --- a/src/i2pp/core/interpolators/interpolator_types.py +++ b/src/i2pp/core/interpolators/interpolator_types.py @@ -39,7 +39,10 @@ class InterpolationType(Enum): ALLVOXELS_WEIGHTED = "allvoxels_weighted" def create_interpolator( - self, *, filter_outliers: bool = False + self, + *, + filter_outliers: bool = False, + set_surf_node_value: float | None = None, ) -> Interpolator: """Creates and returns a configured interpolator instance based on the selected interpolation method. @@ -52,6 +55,8 @@ def create_interpolator( Args: filter_outliers (bool): If True, outliers will be filtered during interpolation. Defaults to False. + set_surf_node_value (float | None): Value to set for surface nodes. + Defaults to None. Returns: Interpolator: An instance of the interpolator that matches the @@ -69,9 +74,13 @@ def create_interpolator( mode="allvoxels_weighted", filter_outliers=filter_outliers ) if self == InterpolationType.NODES: - return InterpolatorNodes(mode="nodes") + return InterpolatorNodes( + mode="nodes", surf_val=set_surf_node_value + ) if self == InterpolationType.NODES_WEIGHTED: - return InterpolatorNodes(mode="nodes_weighted") + return InterpolatorNodes( + mode="nodes_weighted", surf_val=set_surf_node_value + ) if self == InterpolationType.CENTER: return InterpolatorCenter() raise ValueError(f"Unsupported interpolation method: {self}") diff --git a/src/i2pp/core/run.py b/src/i2pp/core/run.py index 5ebc5fe..ef77def 100644 --- a/src/i2pp/core/run.py +++ b/src/i2pp/core/run.py @@ -77,8 +77,7 @@ def run_i2pp(config_i2pp): elements = interpolate_image_to_discretization( discretization, image, - interpolation_method=config.processing.interpolation.method, - filter_outliers=config.processing.interpolation.filter_outliers, + interpolation=config.processing.interpolation, ) # Transform the data using the user-defined python function diff --git a/tests/unittests/i2pp/core/configuration_validator/test_validator.py b/tests/unittests/i2pp/core/configuration_validator/test_validator.py index a5f2526..f0a11db 100644 --- a/tests/unittests/i2pp/core/configuration_validator/test_validator.py +++ b/tests/unittests/i2pp/core/configuration_validator/test_validator.py @@ -28,7 +28,7 @@ def minimal_valid_config(tmp_path): "image": {"path": str(image_path), "type": "dicom"}, }, "processing": { - "interpolation":{ + "interpolation": { "method": "nodes", }, "transformation": { @@ -88,12 +88,9 @@ def large_valid_config(tmp_path): }, }, "processing": { - "smoothing": { - "area": 5, - "visualize": True - }, + "smoothing": {"smoothing_area": 5, "visualize": True}, "interpolation": { - "method": "elements", + "method": "element", }, "transformation": { "user_script": str(script_path), diff --git a/tests/unittests/i2pp/core/test_run.py b/tests/unittests/i2pp/core/test_run.py index 519f8a0..50b2a0b 100644 --- a/tests/unittests/i2pp/core/test_run.py +++ b/tests/unittests/i2pp/core/test_run.py @@ -38,7 +38,9 @@ def minimal_valid_config(tmp_path): }, }, "processing": { - "interpolation": {"method": "nodes"}, + "interpolation": { + "method": "nodes", + }, "transformation": { "user_script": str(script_path), "user_function": "process_image_data", @@ -147,7 +149,12 @@ def maximal_valid_config(tmp_path): }, }, "processing": { - "interpolation": {"method": "nodes"}, + "interpolation": { + "method": "nodes_weighted", + "node_weight": {"interior": 0.8, "surface": 0.2}, + "filter_outliers": True, + "set_surface_node_value": 100.0, + }, "transformation": { "user_script": str(script_path), "user_function": "process_image_data", From 209f25ca0b117ba24eb89b6aa08d9565fb761cc8 Mon Sep 17 00:00:00 2001 From: Saskia Gloy Date: Sun, 4 Jan 2026 16:13:30 +0100 Subject: [PATCH 09/15] Add setting fixed values to boundary elements after interpolation --- .../core/configuration_validator/validator.py | 7 + src/i2pp/core/interpolate_element_data.py | 44 +++++-- .../core/interpolators/interpolator_nodes.py | 10 +- .../core/interpolators/interpolator_types.py | 4 +- .../core/test_interpolate_element_data.py | 121 +++++++++++++++++- 5 files changed, 167 insertions(+), 19 deletions(-) diff --git a/src/i2pp/core/configuration_validator/validator.py b/src/i2pp/core/configuration_validator/validator.py index 6c44eff..f21bd55 100644 --- a/src/i2pp/core/configuration_validator/validator.py +++ b/src/i2pp/core/configuration_validator/validator.py @@ -108,14 +108,21 @@ class Interpolation: method: str filter_outliers: Optional[bool] = field(default=False) set_node_value: Optional[float] = field(default=None) + set_ele_value: Optional[float] = field(default=None) @staticmethod def from_dict(d: Dict[str, Any]) -> "Interpolation": """Creates an Interpolation instance from a dictionary.""" + # check if both set node and element values are provided, which is not allowed + if d.get("set_surface_node_value") is not None and d.get("set_surface_element_value") is not None: + raise ValueError( + "Both 'set_surface_node_value' and 'set_surface_element_value' cannot be set at the same time." + ) return Interpolation( method=d["method"], filter_outliers=d.get("filter_outliers", False), set_node_value=d.get("set_surface_node_value", None), + set_ele_value=d.get("set_surface_element_value", None), ) diff --git a/src/i2pp/core/interpolate_element_data.py b/src/i2pp/core/interpolate_element_data.py index a2cd00a..b022e2b 100644 --- a/src/i2pp/core/interpolate_element_data.py +++ b/src/i2pp/core/interpolate_element_data.py @@ -9,17 +9,33 @@ from i2pp.core.interpolators.interpolator_types import InterpolationType +def _set_surf_element_value( + dis: Discretization, elements: list[Element], value: float | list[float] +) -> None: + """Sets a specified value to all elements containing at least one surface + node.""" + surface_node_ids: set[int] = set() + for surface in dis.surfaces: + surface_node_ids.update(surface.node_ids) + + if not surface_node_ids: + return + + for ele in elements: + if not surface_node_ids.isdisjoint(ele.node_ids): + ele.data = value + + def interpolate_image_to_discretization( dis: Discretization, image_data: ImageData, interpolation: Interpolation, ) -> list[Element]: - """Performs interpolation of image data onto the FEM Discretization based - on the specified interpolation method. + """Performs interpolation of image data onto the FEM Discretization. - This function applies different interpolation methods depending on the - user configuration. The pixel values are assigned to the FEM elements - using one of the following approaches: + This function uses the settings provided in the `interpolation` object + to control the interpolation process. The pixel values are assigned to + the FEM elements using one of the following approaches: - "nodes": Computes the mean pixel value for each element based on its node values. @@ -29,14 +45,16 @@ def interpolate_image_to_discretization( element. Arguments: - dis (Discretization): The Discretization object containing FEM + dis: The Discretization object containing FEM surfaces, elements and node coordinates. - image_data (ImageData): A structured representation containing 3D + image_data: A structured representation containing 3D pixel data, grid coordinates, orientation, and metadata. - interpolation (Interpolation): The interpolation configuration object. + interpolation: The interpolation configuration object, containing + settings like the interpolation method, outlier filtering, and + surface value assignments. Returns: - list[Element]: A list of FEM elements with interpolated pixel data. + A list of FEM elements with interpolated pixel data. """ enum_interpolation_method = InterpolationType(interpolation.method) @@ -45,4 +63,10 @@ def interpolate_image_to_discretization( set_surf_node_value=interpolation.set_surface_node_value, ) - return interpolator.compute_element_data(dis, image_data) + elements = interpolator.compute_element_data(dis, image_data) + + # apply fixed value to elements at boundary if configured + if interpolation.set_ele_value is not None: + _set_surf_element_value(dis, elements, interpolation.set_ele_value) + + return elements diff --git a/src/i2pp/core/interpolators/interpolator_nodes.py b/src/i2pp/core/interpolators/interpolator_nodes.py index a891386..dfe35df 100644 --- a/src/i2pp/core/interpolators/interpolator_nodes.py +++ b/src/i2pp/core/interpolators/interpolator_nodes.py @@ -28,13 +28,13 @@ def __init__( self, *args, mode: str = "nodes", - surf_val: Optional[Union[float, np.ndarray]] = None, + surf_node_val: Optional[Union[float, np.ndarray]] = None, **kwargs, ): """Initializes the InterpolatorNodes.""" super().__init__() self._mode = mode # "nodes" or "nodes_weighted" - self.set_surf_node_value = surf_val + self.set_surf_node_value = surf_node_val # Helpers for readability and error handling def _compute_unweighted_mean( @@ -167,12 +167,12 @@ def _override_surface_nodes( num_values (int): The number of pixel values per node. """ - surf_val = np.asarray(surf_node_value) + surf_node_val = np.asarray(surf_node_value) - if surf_val.size != num_values: + if surf_node_val.size != num_values: raise ValueError( f"set_surf_node_value must have {num_values} value(s), " - f"got {surf_val.size}" + f"got {surf_node_val.size}" ) surface_node_ids = { diff --git a/src/i2pp/core/interpolators/interpolator_types.py b/src/i2pp/core/interpolators/interpolator_types.py index 0045ea6..873b92b 100644 --- a/src/i2pp/core/interpolators/interpolator_types.py +++ b/src/i2pp/core/interpolators/interpolator_types.py @@ -75,11 +75,11 @@ def create_interpolator( ) if self == InterpolationType.NODES: return InterpolatorNodes( - mode="nodes", surf_val=set_surf_node_value + mode="nodes", surf_node_val=set_surf_node_value ) if self == InterpolationType.NODES_WEIGHTED: return InterpolatorNodes( - mode="nodes_weighted", surf_val=set_surf_node_value + mode="nodes_weighted", surf_node_val=set_surf_node_value ) if self == InterpolationType.CENTER: return InterpolatorCenter() diff --git a/tests/unittests/i2pp/core/test_interpolate_element_data.py b/tests/unittests/i2pp/core/test_interpolate_element_data.py index cb5f236..4d22d54 100644 --- a/tests/unittests/i2pp/core/test_interpolate_element_data.py +++ b/tests/unittests/i2pp/core/test_interpolate_element_data.py @@ -1,3 +1,120 @@ -"""Test Interpolator Routine.""" +"""Test setting fixed values for elements at the boundary.""" -pass +import numpy as np +from i2pp.core.configuration_validator.validator import Interpolation +from i2pp.core.discretization_readers.discretization_reader import ( + Discretization, + Element, + Nodes, + Surface, +) +from i2pp.core.image_readers.image_reader import ( + GridCoords, + ImageData, + PixelValueType, +) +from i2pp.core.interpolate_element_data import ( + interpolate_image_to_discretization, +) + + +def _make_simple_discretization(): + """Create a minimal discretization with two elements touching a boundary + surface.""" + # Nodes: 0..3, with coordinates forming a unit square in grid space + node_coords = np.array( + [ + [0, 0, 0], # node 0 + [1, 0, 0], # node 1 + [0, 1, 0], # node 2 + [1, 1, 0], # node 3 + ] + ) + node_ids = np.arange(4) + nodes = Nodes(coords=node_coords, ids=node_ids) + + # Elements: E0 uses nodes [0,1]; E1 uses nodes [2,3] + elements = [ + Element(node_ids=[0, 1], id=0), + Element(node_ids=[2, 3], id=1), + ] + + # Surface contains node 1 and 2 -> E0 and E1 touch boundary + surfaces = [Surface(node_ids=np.array([1, 2], dtype=int), id=0)] + + return Discretization(nodes=nodes, elements=elements, surfaces=surfaces) + + +def _make_ct_image(): + """Create a small 3D scalar CT image with zero values for interpolation + tests.""" + # Simple 3x3x3 scalar image with zeros + pixel_data = np.zeros((3, 3, 3), dtype=float) + grid_coords = GridCoords( + slice=np.arange(3), row=np.arange(3), col=np.arange(3) + ) + return ImageData( + pixel_data=pixel_data, + grid_coords=grid_coords, + orientation=np.eye(3), + position=np.zeros(3), + pixel_type=PixelValueType.CT, + ) + + +def _make_rgb_image(): + """Create a small 3D RGB image with zero values for interpolation tests.""" + # Simple 3x3x3x3 RGB image with zeros + pixel_data = np.zeros((3, 3, 3, 3), dtype=np.uint8) + grid_coords = GridCoords( + slice=np.arange(3), row=np.arange(3), col=np.arange(3) + ) + return ImageData( + pixel_data=pixel_data, + grid_coords=grid_coords, + orientation=np.eye(3), + position=np.zeros(3), + pixel_type=PixelValueType.RGB, + ) + + +def test_boundary_elements_receive_fixed_scalar_value(): + """Verify boundary-touching elements receive the configured fixed scalar + value.""" + dis = _make_simple_discretization() + image = _make_ct_image() + + # Interpolation: nodes, set surface element value to 999.0 + interp_cfg = Interpolation( + method="nodes", + filter_outliers=False, + set_node_value=None, + set_ele_value=999.0, + ) + + elements = interpolate_image_to_discretization(dis, image, interp_cfg) + + # Both elements touch the boundary -> both should receive fixed value + assert float(elements[0].data) == 999.0 + assert float(elements[1].data) == 999.0 + + +def test_boundary_elements_receive_fixed_rgb_value_vector(): + """Verify boundary-touching elements receive the configured fixed RGB + vector.""" + dis = _make_simple_discretization() + image = _make_rgb_image() + + fixed_rgb = [10, 20, 30] + interp_cfg = Interpolation( + method="nodes", + filter_outliers=False, + set_node_value=None, + set_ele_value=fixed_rgb, + ) + + elements = interpolate_image_to_discretization(dis, image, interp_cfg) + + # Both elements touch the boundary -> both should receive fixed RGB vector + assert np.array_equal(np.asarray(elements[0].data), np.asarray(fixed_rgb)) + assert np.array_equal(np.asarray(elements[1].data), np.asarray(fixed_rgb)) From 089b3a3a835b0bb235f58f87c8c954d9e9b2ef04 Mon Sep 17 00:00:00 2001 From: Saskia Gloy Date: Sun, 4 Jan 2026 18:28:07 +0100 Subject: [PATCH 10/15] Enhance element value assignment and handling of zero-sum weights in interpolation --- .../core/configuration_validator/validator.py | 26 +++- .../fourc_yaml_reader.py | 4 +- src/i2pp/core/interpolate_element_data.py | 64 +++++++++- .../interpolators/interpolator_all_voxel.py | 50 ++++++-- .../core/interpolators/interpolator_nodes.py | 6 +- .../core/interpolators/interpolator_types.py | 8 +- src/i2pp/core/run.py | 2 +- templates/config/example_config.yaml | 5 +- tests/integration/test_random_colors.py | 2 - .../test_interpolator_all_voxel.py | 112 ++++++++++++++++++ .../core/test_interpolate_element_data.py | 45 +++++++ 11 files changed, 290 insertions(+), 34 deletions(-) diff --git a/src/i2pp/core/configuration_validator/validator.py b/src/i2pp/core/configuration_validator/validator.py index f21bd55..d4174a5 100644 --- a/src/i2pp/core/configuration_validator/validator.py +++ b/src/i2pp/core/configuration_validator/validator.py @@ -101,6 +101,24 @@ def from_dict(d: Dict[str, Any]) -> "Transformation": ) +@dataclass +class NodeWeight: + """Class representing node weight configuration for interpolation.""" + + surface: float + interior: float + + @staticmethod + def from_dict(d: Optional[Dict[str, Any]]) -> Optional["node_weight"]: + """Creates a NodeWeight instance from a dictionary.""" + if d is None: + return NodeWeight(surface = 1.0, interior = 1.0) + return NodeWeight( + surface=d.get("surface", 1.0), + interior=d.get("interior", 1.0), + ) + + @dataclass class Interpolation: """Class representing the interpolation configuration.""" @@ -108,7 +126,8 @@ class Interpolation: method: str filter_outliers: Optional[bool] = field(default=False) set_node_value: Optional[float] = field(default=None) - set_ele_value: Optional[float] = field(default=None) + set_ele_value: Optional[float | list[float]] = field(default=None) + node_weight: Optional[NodeWeight] = field(default=None) @staticmethod def from_dict(d: Dict[str, Any]) -> "Interpolation": @@ -121,8 +140,9 @@ def from_dict(d: Dict[str, Any]) -> "Interpolation": return Interpolation( method=d["method"], filter_outliers=d.get("filter_outliers", False), - set_node_value=d.get("set_surface_node_value", None), - set_ele_value=d.get("set_surface_element_value", None), + set_node_value=d.get("set_surface_node_value"), + set_ele_value=d.get("set_surface_element_value"), + node_weight=NodeWeight.from_dict(d.get("node_weight")), ) diff --git a/src/i2pp/core/discretization_readers/fourc_yaml_reader.py b/src/i2pp/core/discretization_readers/fourc_yaml_reader.py index 5c25ad9..223a1df 100644 --- a/src/i2pp/core/discretization_readers/fourc_yaml_reader.py +++ b/src/i2pp/core/discretization_readers/fourc_yaml_reader.py @@ -106,8 +106,8 @@ def load_discretization( raw_dis, np.array(options["material_ids"]) ) - interior_weight = processing.interpolation.interior_node_weight - surface_weight = processing.interpolation.surface_node_weight + interior_weight = processing.interpolation.node_weight.interior + surface_weight = processing.interpolation.node_weight.surface nodes_coords = [] node_ids = [] diff --git a/src/i2pp/core/interpolate_element_data.py b/src/i2pp/core/interpolate_element_data.py index b022e2b..c695e94 100644 --- a/src/i2pp/core/interpolate_element_data.py +++ b/src/i2pp/core/interpolate_element_data.py @@ -1,5 +1,8 @@ """Interpolate image data to FEM-Elements.""" +import logging + +import numpy as np from i2pp.core.configuration_validator.validator import Interpolation from i2pp.core.discretization_readers.discretization_reader import ( Discretization, @@ -9,11 +12,41 @@ from i2pp.core.interpolators.interpolator_types import InterpolationType +def _value_compatible(expected: object, value: float | list[float]) -> bool: + """Check if value has the same 'format' as expected element data. + + - Scalars must remain scalars. + - Arrays/lists must match shape and dtype. + """ + # Determine expected format + exp_arr = np.asarray(expected) + # expected scalar + if exp_arr.shape == (): + # value must be scalar-like + val_arr = np.asarray(value) + return val_arr.shape == () + # expected vector/array + val_arr = np.asarray(value) + if val_arr.shape != exp_arr.shape: + return False + # dtype compatibility: allow safe casting + try: + np.asarray(value, dtype=exp_arr.dtype) + return True + except Exception: + return False + + def _set_surf_element_value( dis: Discretization, elements: list[Element], value: float | list[float] ) -> None: """Sets a specified value to all elements containing at least one surface - node.""" + node. + + Applies the provided scalar or vector value to all elements that + touch any surface, regardless of whether element.data was already + computed. + """ surface_node_ids: set[int] = set() for surface in dis.surfaces: surface_node_ids.update(surface.node_ids) @@ -21,9 +54,18 @@ def _set_surf_element_value( if not surface_node_ids: return - for ele in elements: - if not surface_node_ids.isdisjoint(ele.node_ids): - ele.data = value + # Coerce provided value: scalar -> float; array-like -> np.ndarray + val_arr = np.asarray(value) + if val_arr.shape == (): + coerced_scalar = float(val_arr) + for ele in elements: + if not surface_node_ids.isdisjoint(ele.node_ids): + ele.data = coerced_scalar + else: + coerced_vec = val_arr.astype(val_arr.dtype, copy=False) + for ele in elements: + if not surface_node_ids.isdisjoint(ele.node_ids): + ele.data = coerced_vec def interpolate_image_to_discretization( @@ -56,11 +98,23 @@ def interpolate_image_to_discretization( Returns: A list of FEM elements with interpolated pixel data. """ + + if ( + interpolation.set_ele_value is not None + and interpolation.set_node_value is not None + ): + logging.warning( + "Both 'set_surface_element_value' and " + "'set_surface_node_value' are configured; " + "'set_surface_node_value' will be ignored." + ) + interpolation.set_node_value = None + enum_interpolation_method = InterpolationType(interpolation.method) interpolator = enum_interpolation_method.create_interpolator( filter_outliers=interpolation.filter_outliers, - set_surf_node_value=interpolation.set_surface_node_value, + set_node_value=interpolation.set_node_value, ) elements = interpolator.compute_element_data(dis, image_data) diff --git a/src/i2pp/core/interpolators/interpolator_all_voxel.py b/src/i2pp/core/interpolators/interpolator_all_voxel.py index ea7c050..c4e1fc1 100644 --- a/src/i2pp/core/interpolators/interpolator_all_voxel.py +++ b/src/i2pp/core/interpolators/interpolator_all_voxel.py @@ -187,8 +187,15 @@ def _weighted_voxel_mean( filtered_values = values filtered_weights = voxel_weights + # Guard against zero-sum weights + if np.sum(filtered_weights) <= 0: + return np.mean(filtered_values, axis=0) + else: # No outlier filtering for small voxel counts + # Guard against zero-sum weights + if np.sum(voxel_weights) <= 0: + return np.mean(values, axis=0) return np.average(values, weights=voxel_weights, axis=0) # Compute weighted mean @@ -294,21 +301,40 @@ def _get_data_of_element( node_weights_current, ) if self._filter_outliers_enabled - else np.average( - data, - weights=np.sum( - node_weights_current[:, np.newaxis] - / np.maximum( - np.linalg.norm( - element_node_phys[:, np.newaxis, :] - - voxels_phys_np[np.newaxis, :, :], - axis=2, + else ( + # Non-filtered weighted path with zero-sum guard + (np.mean(data, axis=0)) + if np.sum( + np.sum( + node_weights_current[:, np.newaxis] + / np.maximum( + np.linalg.norm( + element_node_phys[:, np.newaxis, :] + - voxels_phys_np[np.newaxis, :, :], + axis=2, + ), + 1e-10, ), - 1e-10, + axis=0, + ) + ) + <= 0 + else np.average( + data, + weights=np.sum( + node_weights_current[:, np.newaxis] + / np.maximum( + np.linalg.norm( + element_node_phys[:, np.newaxis, :] + - voxels_phys_np[np.newaxis, :, :], + axis=2, + ), + 1e-10, + ), + axis=0, ), axis=0, - ), - axis=0, + ) ) ) diff --git a/src/i2pp/core/interpolators/interpolator_nodes.py b/src/i2pp/core/interpolators/interpolator_nodes.py index dfe35df..093fe61 100644 --- a/src/i2pp/core/interpolators/interpolator_nodes.py +++ b/src/i2pp/core/interpolators/interpolator_nodes.py @@ -34,7 +34,7 @@ def __init__( """Initializes the InterpolatorNodes.""" super().__init__() self._mode = mode # "nodes" or "nodes_weighted" - self.set_surf_node_value = surf_node_val + self.set_node_value = surf_node_val # Helpers for readability and error handling def _compute_unweighted_mean( @@ -221,11 +221,11 @@ def compute_element_data( node_grid_coords, image_data ) - if self.set_surf_node_value is not None and dis.surfaces: + if self.set_node_value is not None and dis.surfaces: self._override_surface_nodes( node_values=node_values, dis=dis, - surf_node_value=self.set_surf_node_value, + surf_node_value=self.set_node_value, num_values=image_data.pixel_type.num_values, ) diff --git a/src/i2pp/core/interpolators/interpolator_types.py b/src/i2pp/core/interpolators/interpolator_types.py index 873b92b..459ccc0 100644 --- a/src/i2pp/core/interpolators/interpolator_types.py +++ b/src/i2pp/core/interpolators/interpolator_types.py @@ -42,7 +42,7 @@ def create_interpolator( self, *, filter_outliers: bool = False, - set_surf_node_value: float | None = None, + set_node_value: float | None = None, ) -> Interpolator: """Creates and returns a configured interpolator instance based on the selected interpolation method. @@ -55,7 +55,7 @@ def create_interpolator( Args: filter_outliers (bool): If True, outliers will be filtered during interpolation. Defaults to False. - set_surf_node_value (float | None): Value to set for surface nodes. + set_node_value (float | None): Value to set for surface nodes. Defaults to None. Returns: @@ -75,11 +75,11 @@ def create_interpolator( ) if self == InterpolationType.NODES: return InterpolatorNodes( - mode="nodes", surf_node_val=set_surf_node_value + mode="nodes", surf_node_val=set_node_value ) if self == InterpolationType.NODES_WEIGHTED: return InterpolatorNodes( - mode="nodes_weighted", surf_node_val=set_surf_node_value + mode="nodes_weighted", surf_node_val=set_node_value ) if self == InterpolationType.CENTER: return InterpolatorCenter() diff --git a/src/i2pp/core/run.py b/src/i2pp/core/run.py index ef77def..47abee1 100644 --- a/src/i2pp/core/run.py +++ b/src/i2pp/core/run.py @@ -77,7 +77,7 @@ def run_i2pp(config_i2pp): elements = interpolate_image_to_discretization( discretization, image, - interpolation=config.processing.interpolation, + config.processing.interpolation, ) # Transform the data using the user-defined python function diff --git a/templates/config/example_config.yaml b/templates/config/example_config.yaml index e01aded..95ae208 100644 --- a/templates/config/example_config.yaml +++ b/templates/config/example_config.yaml @@ -16,9 +16,10 @@ import: image_position: [0, 0, 0] processing: smoothing: - smoothing_area: 5 + area: 5 visualize: true - interpolation_method: elements + interpolation: + method: elements transformation: user_script: tests/testdata/user_script.py user_function: process_image_data diff --git a/tests/integration/test_random_colors.py b/tests/integration/test_random_colors.py index 111bdfd..f246a0c 100644 --- a/tests/integration/test_random_colors.py +++ b/tests/integration/test_random_colors.py @@ -6,10 +6,8 @@ import tempfile import jinja2 -import pytest -@pytest.mark.integration def test_i2pp_integration_random_colors(): """Test that i2pp generates the expected output for a sample image data set with random colors. diff --git a/tests/unittests/i2pp/core/interpolators/test_interpolator_all_voxel.py b/tests/unittests/i2pp/core/interpolators/test_interpolator_all_voxel.py index a670d3a..2603164 100644 --- a/tests/unittests/i2pp/core/interpolators/test_interpolator_all_voxel.py +++ b/tests/unittests/i2pp/core/interpolators/test_interpolator_all_voxel.py @@ -573,3 +573,115 @@ def test_weighted_small_voxel_count_uses_simple_weighted_average(): expected = np.average(values, weights=voxel_weights, axis=0) assert np.isclose(result, expected) + + +def test_weighted_mode_no_filter_falls_back_when_weights_sum_zero(): + """Weighted/no-filter path should fall back to unweighted mean when weights + sum to zero.""" + # Simple 3x3x3 grid with deterministic values + pixel_data = np.arange(27, dtype=float).reshape((3, 3, 3)) + grid_coords = GridCoords( + slice=np.arange(3), row=np.arange(3), col=np.arange(3) + ) + image_data = ImageData( + pixel_data, grid_coords, np.eye(3), np.zeros(3), PixelValueType.CT + ) + + # Element spanning [0,1] cube corners + element_nodes = np.array( + [ + [0, 0, 0], + [1, 0, 0], + [0, 1, 0], + [1, 1, 0], + [0, 0, 1], + [1, 0, 1], + [0, 1, 1], + [1, 1, 1], + ] + ) + + # Discretization: one element using these nodes; + # set all node weights to zero + dis = _make_simple_discretization( + node_coords=element_nodes, + elements_node_ids=[np.arange(8)], + node_weights=np.zeros(8, dtype=float), + ) + + # Compute with weighted mode, no outlier filtering + interp_weighted = InterpolatorAllVoxel( + mode="allvoxels_weighted", filter_outliers=False + ) + elems = interp_weighted.compute_element_data(dis, image_data) + val_weighted = float(np.asarray(elems[0].data).mean()) + + # Compute unweighted reference + interp_unweighted = InterpolatorAllVoxel( + mode="allvoxels", filter_outliers=False + ) + val_unweighted = float( + np.asarray( + interp_unweighted.compute_element_data(dis, image_data)[0].data + ).mean() + ) + + # Should match unweighted mean (fallback) and not raise ZeroDivisionError + assert np.isclose(val_weighted, val_unweighted) + + +def test_weighted_mode_with_filter_falls_back_when_weights_sum_zero(): + """Weighted/filtered path should fall back to unweighted mean when weights + sum to zero.""" + # Grid with ones and an outlier to exercise filtering path + N = 5 + pixel_data = np.ones((N, N, N), dtype=float) + pixel_data[2, 2, 2] = 1000.0 + grid_coords = GridCoords( + slice=np.arange(N), row=np.arange(N), col=np.arange(N) + ) + image_data = ImageData( + pixel_data, grid_coords, np.eye(3), np.zeros(3), PixelValueType.CT + ) + + # Element covering a cube [1,3] in each axis to include the outlier voxel + element_nodes = np.array( + [ + [1, 1, 1], + [3, 1, 1], + [1, 3, 1], + [3, 3, 1], + [1, 1, 3], + [3, 1, 3], + [1, 3, 3], + [3, 3, 3], + ] + ) + + # Zero node weights to trigger zero-sum voxel weights + dis = _make_simple_discretization( + node_coords=element_nodes, + elements_node_ids=[np.arange(8)], + node_weights=np.zeros(8, dtype=float), + ) + + # Weighted with filtering; should not raise + # and should approximate unweighted filtered mean + interp_weighted_filter = InterpolatorAllVoxel( + mode="allvoxels_weighted", filter_outliers=True + ) + elems_weighted = interp_weighted_filter.compute_element_data( + dis, image_data + ) + val_weighted = float(np.asarray(elems_weighted[0].data).mean()) + + # Unweighted with filtering as reference + interp_unweighted_filter = InterpolatorAllVoxel( + mode="allvoxels", filter_outliers=True + ) + elems_unweighted = interp_unweighted_filter.compute_element_data( + dis, image_data + ) + val_unweighted = float(np.asarray(elems_unweighted[0].data).mean()) + + assert np.isclose(val_weighted, val_unweighted, rtol=1e-6) diff --git a/tests/unittests/i2pp/core/test_interpolate_element_data.py b/tests/unittests/i2pp/core/test_interpolate_element_data.py index 4d22d54..c9ec09e 100644 --- a/tests/unittests/i2pp/core/test_interpolate_element_data.py +++ b/tests/unittests/i2pp/core/test_interpolate_element_data.py @@ -90,6 +90,7 @@ def test_boundary_elements_receive_fixed_scalar_value(): filter_outliers=False, set_node_value=None, set_ele_value=999.0, + node_weight=None, ) elements = interpolate_image_to_discretization(dis, image, interp_cfg) @@ -111,6 +112,7 @@ def test_boundary_elements_receive_fixed_rgb_value_vector(): filter_outliers=False, set_node_value=None, set_ele_value=fixed_rgb, + node_weight=None, ) elements = interpolate_image_to_discretization(dis, image, interp_cfg) @@ -118,3 +120,46 @@ def test_boundary_elements_receive_fixed_rgb_value_vector(): # Both elements touch the boundary -> both should receive fixed RGB vector assert np.array_equal(np.asarray(elements[0].data), np.asarray(fixed_rgb)) assert np.array_equal(np.asarray(elements[1].data), np.asarray(fixed_rgb)) + + +def test_boundary_elements_receive_fixed_value_without_prior_data(): + """Boundary-touching elements get fixed values even if element.data was not + set before.""" + dis = _make_simple_discretization() + image = _make_ct_image() + + # Configure interpolation with method that will compute, + # but we simulate unset data by overwriting after + interp_cfg = Interpolation( + method="nodes", + filter_outliers=False, + set_node_value=None, + set_ele_value=123.0, + node_weight=None, + ) + + # Perform interpolation + elements = interpolate_image_to_discretization(dis, image, interp_cfg) + + # Verify fixed scalar applied + assert float(elements[0].data) == 123.0 + assert float(elements[1].data) == 123.0 + + # Now test vector assignment on RGB image + image_rgb = _make_rgb_image() + interp_cfg_vec = Interpolation( + method="nodes", + filter_outliers=False, + set_node_value=None, + set_ele_value=[7, 8, 9], + node_weight=None, + ) + elements_vec = interpolate_image_to_discretization( + dis, image_rgb, interp_cfg_vec + ) + assert np.array_equal( + np.asarray(elements_vec[0].data), np.array([7, 8, 9]) + ) + assert np.array_equal( + np.asarray(elements_vec[1].data), np.array([7, 8, 9]) + ) From 9d1e5aae1a54ff7ea4f9c03880901430008059ee Mon Sep 17 00:00:00 2001 From: Saskia Gloy Date: Sun, 4 Jan 2026 19:10:59 +0100 Subject: [PATCH 11/15] Update README and example config for interpolation methods and value settings --- README.md | 25 ++++++++++++++++++++++++- 1 file changed, 24 insertions(+), 1 deletion(-) diff --git a/README.md b/README.md index d8fbff5..9268862 100644 --- a/README.md +++ b/README.md @@ -24,6 +24,7 @@ The remaining parts of the README are structured as follows: - [Execute img2physiprop](#execute-img2physiprop) - [Run testing framework and create coverage report](#run-testing-framework-and-create-coverage-report) - [Create documentation](#create-documentation) + - [Interpolation and value settings](#interpolation-and-value-settings) - [Dependency Management](#dependency-management) - [Contributing](#contributing) - [License](#license) @@ -75,7 +76,7 @@ To execute img2physiprop run ``` i2pp --config path/to/config.yaml -```` +``` with your custom configuration file. A template configuration file containing all possible input configurations can be found in the folder `templates/config`. @@ -95,6 +96,28 @@ To locally create the documentation from the provided docstrings simply run pdoc --html --output-dir docs src/i2pp ``` +### Interpolation and value settings + +- Interpolation methods (processing.interpolation.method): + - nodes: Interpolates values at the element’s nodes and assigns the element mean (ignoring NaN nodes). Fast and robust; respects node sampling. + - nodes_weighted: Like `nodes`, but computes a weighted mean using node-specific weights (`dis.nodes.weights`), which are set via `processing.interpolation.node_weight.surface` and `processing.interpolation.node_weight.interior`. This reduces the influence of low-weight nodes. + - elementcenter: Interpolates at each element centroid and assigns that value. Reliable fallback if no voxels lie inside the element. + - allvoxels: Collects all voxels whose grid coordinates lie inside the convex hull of the element nodes; assigns the mean value; optionally filters outliers. + - allvoxels_weighted: Computes a voxel-weighted mean where voxel weights derive from node weights and inverse node-to-voxel distances; optionally filters outliers. + +- Element and node value overrides: + - set_surface_node_value: If provided, all nodes that belong to any surface receive the fixed value (vector size must match the number of pixel channels); only relevant for `nodes`and `nodes_weighted` interpolation method. + - set_surface_element_value: If provided, all elements touching any surface node receive the fixed value (scalar or vector); applicable to all interpolation methods. + - When both are set, surface node override is ignored (element override wins). + +- Outlier filtering (processing.interpolation.filter_outliers): + - In `allvoxels` and `allvoxels_weighted`, if enabled and enough voxels are present (>5), outliers are removed using a modified Z-score (median/MAD-based, threshold=3.5) before averaging. + +- Fallbacks and warnings: + - If outlier filtering removes all voxels, the method falls back to the unfiltered mean. + - If an element contains no voxels (allvoxels modes), interpolation falls back to the element center. + - If interpolated points fall outside the image grid, element data is NaN and a warning summary is logged after processing. + ## Dependency Management To ease the dependency update process [`pip-tools`](https://github.com/jazzband/pip-tools) is utilized. To create the necessary [`requirements.txt`](./requirements.txt) file simply execute From 81e812d407af28b9ccdf38ace2a434966aa3526f Mon Sep 17 00:00:00 2001 From: Saskia Gloy Date: Thu, 29 Jan 2026 12:22:58 +0100 Subject: [PATCH 12/15] Refine interpolation methods description in README Removed redundant text from interpolation methods section. --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index 9268862..fd10f01 100644 --- a/README.md +++ b/README.md @@ -101,7 +101,7 @@ pdoc --html --output-dir docs src/i2pp - Interpolation methods (processing.interpolation.method): - nodes: Interpolates values at the element’s nodes and assigns the element mean (ignoring NaN nodes). Fast and robust; respects node sampling. - nodes_weighted: Like `nodes`, but computes a weighted mean using node-specific weights (`dis.nodes.weights`), which are set via `processing.interpolation.node_weight.surface` and `processing.interpolation.node_weight.interior`. This reduces the influence of low-weight nodes. - - elementcenter: Interpolates at each element centroid and assigns that value. Reliable fallback if no voxels lie inside the element. + - elementcenter: Interpolates at each element centroid and assigns that value. - allvoxels: Collects all voxels whose grid coordinates lie inside the convex hull of the element nodes; assigns the mean value; optionally filters outliers. - allvoxels_weighted: Computes a voxel-weighted mean where voxel weights derive from node weights and inverse node-to-voxel distances; optionally filters outliers. From ebf8a79ac249dfe4ea9a77da6d75b041db90c1f0 Mon Sep 17 00:00:00 2001 From: Saskia Gloy Date: Thu, 12 Feb 2026 10:51:19 +0100 Subject: [PATCH 13/15] implement copilot PR draft suggestions --- .../core/configuration_validator/validator.py | 38 ++++++-- src/i2pp/core/discretization_helpers.py | 4 +- .../discretization_reader.py | 6 +- .../fourc_yaml_reader.py | 9 +- .../discretization_readers/mesh_reader.py | 9 +- src/i2pp/core/interpolate_element_data.py | 38 -------- .../interpolators/interpolator_all_voxel.py | 96 +++++++++---------- .../core/interpolators/interpolator_types.py | 2 +- src/i2pp/core/run.py | 6 +- src/i2pp/core/utilities.py | 15 ++- templates/config/example_config.yaml | 2 +- tests/integration/test_random_colors.py | 2 + .../configuration_validator/test_validator.py | 4 +- .../test_mesh_reader.py | 25 ++++- tests/unittests/i2pp/core/test_run.py | 2 +- 15 files changed, 136 insertions(+), 122 deletions(-) diff --git a/src/i2pp/core/configuration_validator/validator.py b/src/i2pp/core/configuration_validator/validator.py index d4174a5..e68d27d 100644 --- a/src/i2pp/core/configuration_validator/validator.py +++ b/src/i2pp/core/configuration_validator/validator.py @@ -8,6 +8,7 @@ from i2pp.core.configuration_validator.validation_helpers import ( resolve_and_validate_path, ) +from i2pp.core.interpolators.interpolator_types import InterpolationType @dataclass @@ -75,6 +76,13 @@ def from_dict(d: Optional[Dict[str, Any]]) -> Optional["Smoothing"]: """Creates a Smoothing instance from a dictionary.""" if d is None: return None + if d.get("area", 3) <= 0: + raise ValueError("Smoothing area must be a positive integer.") + if "smoothing_area" in d: + raise ValueError( + "The key 'smoothing_area' is deprecated. " + "Please use 'area' instead." + ) return Smoothing( area=d.get("area", 3), visualize=d.get("visualize", False), @@ -109,10 +117,10 @@ class NodeWeight: interior: float @staticmethod - def from_dict(d: Optional[Dict[str, Any]]) -> Optional["node_weight"]: + def from_dict(d: Optional[Dict[str, Any]]) -> "NodeWeight": """Creates a NodeWeight instance from a dictionary.""" if d is None: - return NodeWeight(surface = 1.0, interior = 1.0) + return NodeWeight(surface=1.0, interior=1.0) return NodeWeight( surface=d.get("surface", 1.0), interior=d.get("interior", 1.0), @@ -124,18 +132,32 @@ class Interpolation: """Class representing the interpolation configuration.""" method: str - filter_outliers: Optional[bool] = field(default=False) - set_node_value: Optional[float] = field(default=None) + node_weight: NodeWeight + filter_outliers: bool = False + set_node_value: Optional[float | list[float]] = field(default=None) set_ele_value: Optional[float | list[float]] = field(default=None) - node_weight: Optional[NodeWeight] = field(default=None) @staticmethod def from_dict(d: Dict[str, Any]) -> "Interpolation": """Creates an Interpolation instance from a dictionary.""" - # check if both set node and element values are provided, which is not allowed - if d.get("set_surface_node_value") is not None and d.get("set_surface_element_value") is not None: + + method = d["method"] + try: + InterpolationType(method) + except ValueError as error: + allowed = ", ".join(t.value for t in InterpolationType) + raise ValueError( + f"Unsupported interpolation method '{method}'. " + f"Supported methods are: {allowed}." + ) from error + if ( + d.get("set_surface_node_value") is not None + and d.get("set_surface_element_value") is not None + ): raise ValueError( - "Both 'set_surface_node_value' and 'set_surface_element_value' cannot be set at the same time." + "Both 'set_surface_node_value' " + "and 'set_surface_element_value' " + "cannot be set at the same time." ) return Interpolation( method=d["method"], diff --git a/src/i2pp/core/discretization_helpers.py b/src/i2pp/core/discretization_helpers.py index a4f0649..51cddf2 100644 --- a/src/i2pp/core/discretization_helpers.py +++ b/src/i2pp/core/discretization_helpers.py @@ -4,6 +4,7 @@ import numpy as np import pyvista as pv +from i2pp.core.configuration_validator.validator import Processing from i2pp.core.discretization_readers.discretization_format import ( DiscretizationFormat, ) @@ -52,7 +53,7 @@ def determine_discretization_format(file_path: Path) -> DiscretizationFormat: def verify_and_load_discretization( discretization_path: Path, options: dict, - processing: dict, + processing: Processing, ) -> Discretization: """Loads and processes mesh data. @@ -64,7 +65,6 @@ def verify_and_load_discretization( discretization_path (Path): Path to the discretization file. options (dict): Options for loading the discretization that are passed to the reader classes. - processing (dict): Processing options are passed to the reader classes. Returns: DiscretizationData: The loaded and processed mesh data. diff --git a/src/i2pp/core/discretization_readers/discretization_reader.py b/src/i2pp/core/discretization_readers/discretization_reader.py index 3cc8bff..6e37637 100644 --- a/src/i2pp/core/discretization_readers/discretization_reader.py +++ b/src/i2pp/core/discretization_readers/discretization_reader.py @@ -6,6 +6,7 @@ from typing import Optional import numpy as np +from i2pp.core.configuration_validator.validator import Processing @dataclass @@ -70,7 +71,7 @@ class Surface: """Class for storing information about a surface in the Discretization. This class represents a surface in the mesh Discretization, defined - by its node IDs, an surface ID, and associated data. + by its node IDs, a surface ID, and associated data. Attributes: node_ids (np.ndarray): An array of node IDs that define the nodes of @@ -126,7 +127,7 @@ def __init__(self): @abstractmethod def load_discretization( - self, file_path: Path, options: dict, processing: dict + self, file_path: Path, options: dict, processing: Processing ) -> Discretization: """Abstract method to load discretization data from a file path. @@ -140,6 +141,7 @@ def load_discretization( file_path (Path): Path to the discretization file. options (dict): A dictionary containing options for loading the discretization. + processing (Processing): Processing configuration object. Returns: Discretization: An instance of Discretization containing diff --git a/src/i2pp/core/discretization_readers/fourc_yaml_reader.py b/src/i2pp/core/discretization_readers/fourc_yaml_reader.py index 223a1df..fc85907 100644 --- a/src/i2pp/core/discretization_readers/fourc_yaml_reader.py +++ b/src/i2pp/core/discretization_readers/fourc_yaml_reader.py @@ -5,6 +5,7 @@ import lnmmeshio import numpy as np +from i2pp.core.configuration_validator.validator import Processing from i2pp.core.discretization_readers.discretization_reader import ( Discretization, DiscretizationReader, @@ -74,7 +75,7 @@ def load_discretization( self, file_path: Path, options: dict, - processing: "I2PPConfig.processing", + processing: Processing, ) -> Discretization: """Loads and processes a finite element discretization from a .4C.yaml file. @@ -88,7 +89,8 @@ def load_discretization( options (dict): Options for loading the discretization. Filtering for material ids can be enabled by specifying `material_ids` in the options dictionary. - processing (I2PPConfig.processing): Processing configuration object. + processing (I2PPConfig.processing): + Processing configuration object. Returns: Discretization: The finite element discretization including nodes @@ -131,11 +133,12 @@ def load_discretization( surfaces = [] + node_id_to_idx = {nid: i for i, nid in enumerate(node_ids)} for surf in raw_dis.surfacenodesets: surf_node_ids = [] for node in surf.nodes: surf_node_ids.append(node.id) - position = node_ids.index(node.id) + position = node_id_to_idx[node.id] nodes_weights[position] = surface_weight surfaces.append( diff --git a/src/i2pp/core/discretization_readers/mesh_reader.py b/src/i2pp/core/discretization_readers/mesh_reader.py index 09a9848..53f448f 100644 --- a/src/i2pp/core/discretization_readers/mesh_reader.py +++ b/src/i2pp/core/discretization_readers/mesh_reader.py @@ -4,6 +4,7 @@ from pathlib import Path import trimesh +from i2pp.core.configuration_validator.validator import Processing from i2pp.core.discretization_readers.discretization_reader import ( Discretization, DiscretizationReader, @@ -33,7 +34,10 @@ def _filter_discretization(self) -> None: raise RuntimeError("This function is not implemented yet.") def load_discretization( - self, file_path: Path, options: dict + self, + file_path: Path, + options: dict, + processing: Processing, ) -> Discretization: """Loads and processes a finite element model from a .mesh file. @@ -65,8 +69,7 @@ def load_discretization( for i, face in enumerate(raw_dis.faces): elements.append(Element(node_ids=face, id=i)) - # TODO: Add surfaces if needed in the future - # surfaces = [] + # Add surfaces here if needed in the future dis = Discretization(nodes=nodes, elements=elements, surfaces=[]) diff --git a/src/i2pp/core/interpolate_element_data.py b/src/i2pp/core/interpolate_element_data.py index c695e94..9b145c0 100644 --- a/src/i2pp/core/interpolate_element_data.py +++ b/src/i2pp/core/interpolate_element_data.py @@ -1,7 +1,5 @@ """Interpolate image data to FEM-Elements.""" -import logging - import numpy as np from i2pp.core.configuration_validator.validator import Interpolation from i2pp.core.discretization_readers.discretization_reader import ( @@ -12,31 +10,6 @@ from i2pp.core.interpolators.interpolator_types import InterpolationType -def _value_compatible(expected: object, value: float | list[float]) -> bool: - """Check if value has the same 'format' as expected element data. - - - Scalars must remain scalars. - - Arrays/lists must match shape and dtype. - """ - # Determine expected format - exp_arr = np.asarray(expected) - # expected scalar - if exp_arr.shape == (): - # value must be scalar-like - val_arr = np.asarray(value) - return val_arr.shape == () - # expected vector/array - val_arr = np.asarray(value) - if val_arr.shape != exp_arr.shape: - return False - # dtype compatibility: allow safe casting - try: - np.asarray(value, dtype=exp_arr.dtype) - return True - except Exception: - return False - - def _set_surf_element_value( dis: Discretization, elements: list[Element], value: float | list[float] ) -> None: @@ -99,17 +72,6 @@ def interpolate_image_to_discretization( A list of FEM elements with interpolated pixel data. """ - if ( - interpolation.set_ele_value is not None - and interpolation.set_node_value is not None - ): - logging.warning( - "Both 'set_surface_element_value' and " - "'set_surface_node_value' are configured; " - "'set_surface_node_value' will be ignored." - ) - interpolation.set_node_value = None - enum_interpolation_method = InterpolationType(interpolation.method) interpolator = enum_interpolation_method.create_interpolator( diff --git a/src/i2pp/core/interpolators/interpolator_all_voxel.py b/src/i2pp/core/interpolators/interpolator_all_voxel.py index c4e1fc1..47ce82d 100644 --- a/src/i2pp/core/interpolators/interpolator_all_voxel.py +++ b/src/i2pp/core/interpolators/interpolator_all_voxel.py @@ -201,6 +201,26 @@ def _weighted_voxel_mean( # Compute weighted mean return np.average(filtered_values, weights=filtered_weights, axis=0) + def _format_output_value( + self, value: float | np.ndarray, image_data: ImageData + ) -> np.ndarray: + """Formats the output value based on the pixel type. + + If the pixel type has one value, it wraps the value in a NumPy array. + If it has multiple values, it creates a NumPy array filled with the + value. + + Args: + value (float | np.ndarray): The value to format. + image_data (ImageData): Image data containing pixel type info. + + Returns: + np.ndarray: The formatted output value. + """ + if image_data.pixel_type.num_values == 1: + return np.array([value]) + return np.full(image_data.pixel_type.num_values, value) + def _get_data_of_element( self, element_node_grid_coords: np.ndarray, @@ -282,75 +302,45 @@ def _get_data_of_element( ) else: mean_val = np.mean(data, axis=0) - return ( - np.array([mean_val]) - if image_data.pixel_type.num_values == 1 - else np.full(image_data.pixel_type.num_values, mean_val) - ) + return self._format_output_value(mean_val, image_data) if self._mode == "allvoxels_weighted": element_node_phys = element_node_grid_coords if node_weights_current is None: weighted = np.mean(data, axis=0) else: - weighted = ( - self._weighted_voxel_mean( + if self._filter_outliers_enabled: + weighted = self._weighted_voxel_mean( voxels_phys_np, data, element_node_phys, node_weights_current, ) - if self._filter_outliers_enabled - else ( - # Non-filtered weighted path with zero-sum guard - (np.mean(data, axis=0)) - if np.sum( - np.sum( - node_weights_current[:, np.newaxis] - / np.maximum( - np.linalg.norm( - element_node_phys[:, np.newaxis, :] - - voxels_phys_np[np.newaxis, :, :], - axis=2, - ), - 1e-10, - ), - axis=0, - ) - ) - <= 0 - else np.average( - data, - weights=np.sum( - node_weights_current[:, np.newaxis] - / np.maximum( - np.linalg.norm( - element_node_phys[:, np.newaxis, :] - - voxels_phys_np[np.newaxis, :, :], - axis=2, - ), - 1e-10, - ), - axis=0, - ), - axis=0, - ) + else: + # Non-filtered weighted path with zero-sum guard + distances = np.linalg.norm( + element_node_phys[:, np.newaxis, :] + - voxels_phys_np[np.newaxis, :, :], + axis=2, + ) + distances = np.maximum(distances, 1e-10) + voxel_weights = np.sum( + node_weights_current[:, np.newaxis] / distances, + axis=0, ) - ) - return ( - np.array([weighted]) - if image_data.pixel_type.num_values == 1 - else np.full(image_data.pixel_type.num_values, weighted) - ) + if np.sum(voxel_weights) <= 0: + weighted = np.mean(data, axis=0) + else: + weighted = np.average( + data, weights=voxel_weights, axis=0 + ) + + return self._format_output_value(weighted, image_data) # Unknown mode fallback mean_val = np.mean(data, axis=0) - return ( - np.array([mean_val]) - if image_data.pixel_type.num_values == 1 - else np.full(image_data.pixel_type.num_values, mean_val) - ) + return self._format_output_value(mean_val, image_data) else: self.backup_interpolation += 1 diff --git a/src/i2pp/core/interpolators/interpolator_types.py b/src/i2pp/core/interpolators/interpolator_types.py index 459ccc0..78a7bc2 100644 --- a/src/i2pp/core/interpolators/interpolator_types.py +++ b/src/i2pp/core/interpolators/interpolator_types.py @@ -42,7 +42,7 @@ def create_interpolator( self, *, filter_outliers: bool = False, - set_node_value: float | None = None, + set_node_value: float | list[float] | None = None, ) -> Interpolator: """Creates and returns a configured interpolator instance based on the selected interpolation method. diff --git a/src/i2pp/core/run.py b/src/i2pp/core/run.py index 47abee1..092d327 100644 --- a/src/i2pp/core/run.py +++ b/src/i2pp/core/run.py @@ -55,15 +55,15 @@ def run_i2pp(config_i2pp): discretization.bounding_box, ) - # Build exact mesh mask - image.mask = create_mesh_mask(discretization, image) - # If smoothing is enabled, smooth the image data if config.processing.smoothing: # Create a copy of the image data for visualization purposes if config.processing.smoothing.visualize: image_raw = copy.deepcopy(image) + # Build exact mesh mask + image.mask = create_mesh_mask(discretization, image) + image.pixel_data = smooth_data( image.pixel_data, config.processing.smoothing.area, diff --git a/src/i2pp/core/utilities.py b/src/i2pp/core/utilities.py index a484835..afe436f 100644 --- a/src/i2pp/core/utilities.py +++ b/src/i2pp/core/utilities.py @@ -104,8 +104,10 @@ def smooth_data( - Zero data outside mask → data_masked - Smooth the zeroed data → smoothed_data (artificially darkened at edges) - Smooth the mask → smoothed_mask (tells you the "weight" of valid data) - - Divide smoothed data by smoothed mask → This "undoes" the darkening by renormalizing - - Only update pixels that were originally inside the mask → Preserves original data outside + - Divide smoothed data by smoothed mask + → This "undoes" the darkening by renormalizing + - Only update pixels that were originally inside the mask + → Preserves original data outside Args: data (np.ndarray): A 3D array containing the pixel data to be smoothed. @@ -134,7 +136,7 @@ def smooth_data( orig_dtype = data.dtype data_f = data.astype(np.float32, copy=False) - + # Create a copy of the data where all values outside the mask are zero. data_masked = data_f.copy() data_masked[~mask, ...] = 0 @@ -151,7 +153,12 @@ def smooth_data( axes=(0, 1, 2), ) - smoothed_mask_exp = smoothed_mask[..., None] + # Reshape smoothed_mask to match smoothed_data for broadcasting + target_shape = list(smoothed_mask.shape) + [1] * ( + smoothed_data.ndim - smoothed_mask.ndim + ) + smoothed_mask_exp = smoothed_mask.reshape(target_shape) + valid_mask = (mask) & (smoothed_mask > 1e-8) result_f = data_f.copy() diff --git a/templates/config/example_config.yaml b/templates/config/example_config.yaml index 95ae208..2c362b3 100644 --- a/templates/config/example_config.yaml +++ b/templates/config/example_config.yaml @@ -19,7 +19,7 @@ processing: area: 5 visualize: true interpolation: - method: elements + method: elementcenter transformation: user_script: tests/testdata/user_script.py user_function: process_image_data diff --git a/tests/integration/test_random_colors.py b/tests/integration/test_random_colors.py index f246a0c..111bdfd 100644 --- a/tests/integration/test_random_colors.py +++ b/tests/integration/test_random_colors.py @@ -6,8 +6,10 @@ import tempfile import jinja2 +import pytest +@pytest.mark.integration def test_i2pp_integration_random_colors(): """Test that i2pp generates the expected output for a sample image data set with random colors. diff --git a/tests/unittests/i2pp/core/configuration_validator/test_validator.py b/tests/unittests/i2pp/core/configuration_validator/test_validator.py index f0a11db..bf91512 100644 --- a/tests/unittests/i2pp/core/configuration_validator/test_validator.py +++ b/tests/unittests/i2pp/core/configuration_validator/test_validator.py @@ -88,9 +88,9 @@ def large_valid_config(tmp_path): }, }, "processing": { - "smoothing": {"smoothing_area": 5, "visualize": True}, + "smoothing": {"area": 5, "visualize": True}, "interpolation": { - "method": "element", + "method": "elementcenter", }, "transformation": { "user_script": str(script_path), diff --git a/tests/unittests/i2pp/core/discretization_readers/test_mesh_reader.py b/tests/unittests/i2pp/core/discretization_readers/test_mesh_reader.py index e39f5b1..cb61afc 100644 --- a/tests/unittests/i2pp/core/discretization_readers/test_mesh_reader.py +++ b/tests/unittests/i2pp/core/discretization_readers/test_mesh_reader.py @@ -3,6 +3,12 @@ from pathlib import Path from unittest.mock import patch +from i2pp.core.configuration_validator.validator import ( + Interpolation, + NodeWeight, + Processing, + Transformation, +) from i2pp.core.discretization_readers.mesh_reader import MeshReader @@ -13,6 +19,21 @@ def test_load_discretization_mesh(tmp_path: Path) -> None: test_dis = MeshReader() test_config = {"material_ids": None} + # create dummy processing config + test_processing = Processing( + smoothing=None, + interpolation=Interpolation( + method="nodes", + filter_outliers=False, + node_weight=NodeWeight(interior=0.5, surface=0.5), + ), + transformation=Transformation( + user_script=Path(""), + user_function="", + normalize_values=False, + visualize=False, + ), + ) with patch( "i2pp.core.discretization_readers.mesh_reader.Discretization", @@ -20,6 +41,8 @@ def test_load_discretization_mesh(tmp_path: Path) -> None: ) as MockClass: with patch("trimesh.load", returnValue=None) as mock_trimesh: - test_dis.load_discretization(Path(test_path), test_config) + test_dis.load_discretization( + Path(test_path), test_config, test_processing + ) assert mock_trimesh.call_count == 1 assert MockClass.call_count == 1 diff --git a/tests/unittests/i2pp/core/test_run.py b/tests/unittests/i2pp/core/test_run.py index 50b2a0b..735fc3c 100644 --- a/tests/unittests/i2pp/core/test_run.py +++ b/tests/unittests/i2pp/core/test_run.py @@ -100,7 +100,7 @@ def test_run_i2pp_with_minimal_valid_config( # Assertions to ensure key steps were called mock_verify_discretization.assert_called_once() mock_verify_image.assert_called_once() - mock_create_mesh_mask.assert_called_once() + mock_create_mesh_mask.assert_not_called() mock_smooth_data.assert_not_called() mock_visualize_smoothing.assert_not_called() mock_interpolate.assert_called_once() From e36699c14fd8644082f65e8a24c3844634d6b025 Mon Sep 17 00:00:00 2001 From: Saskia Gloy Date: Tue, 24 Feb 2026 13:31:11 +0100 Subject: [PATCH 14/15] refactor/rename in interpolation methods to use scaling factors instead of weights --- README.md | 8 +- .../core/configuration_validator/validator.py | 30 ++-- .../discretization_reader.py | 21 ++- .../fourc_yaml_reader.py | 27 ++- .../discretization_readers/mesh_reader.py | 7 +- .../interpolators/interpolator_all_voxel.py | 52 +++--- .../core/interpolators/interpolator_nodes.py | 60 ++++--- .../core/interpolators/interpolator_types.py | 20 +-- .../test_fourc_yaml_reader.py | 7 +- .../test_mesh_reader.py | 6 +- .../test_interpolator_all_voxel.py | 159 +++++++++--------- .../interpolators/test_interpolator_nodes.py | 20 +-- .../interpolators/test_interpolator_types.py | 10 +- .../core/test_interpolate_element_data.py | 8 +- tests/unittests/i2pp/core/test_run.py | 4 +- 15 files changed, 245 insertions(+), 194 deletions(-) diff --git a/README.md b/README.md index fd10f01..fac6eec 100644 --- a/README.md +++ b/README.md @@ -100,18 +100,18 @@ pdoc --html --output-dir docs src/i2pp - Interpolation methods (processing.interpolation.method): - nodes: Interpolates values at the element’s nodes and assigns the element mean (ignoring NaN nodes). Fast and robust; respects node sampling. - - nodes_weighted: Like `nodes`, but computes a weighted mean using node-specific weights (`dis.nodes.weights`), which are set via `processing.interpolation.node_weight.surface` and `processing.interpolation.node_weight.interior`. This reduces the influence of low-weight nodes. + - nodes_scaled: Like `nodes`, but computes a scaled mean using node-specific scaling factors (`dis.nodes.scaling_factors`), which are set via `processing.interpolation.node_scaling_factors.surface` and `processing.interpolation.node_scaling_factors.interior`. This adjusts the influence of specific nodes. - elementcenter: Interpolates at each element centroid and assigns that value. - allvoxels: Collects all voxels whose grid coordinates lie inside the convex hull of the element nodes; assigns the mean value; optionally filters outliers. - - allvoxels_weighted: Computes a voxel-weighted mean where voxel weights derive from node weights and inverse node-to-voxel distances; optionally filters outliers. + - allvoxels_scaled: Computes a voxel-weighted mean where voxel weights derive from node scaling factors and inverse node-to-voxel distances; optionally filters outliers. - Element and node value overrides: - - set_surface_node_value: If provided, all nodes that belong to any surface receive the fixed value (vector size must match the number of pixel channels); only relevant for `nodes`and `nodes_weighted` interpolation method. + - set_surface_node_value: If provided, all nodes that belong to any surface receive the fixed value (vector size must match the number of pixel channels); only relevant for `nodes` and `nodes_scaled` interpolation methods. - set_surface_element_value: If provided, all elements touching any surface node receive the fixed value (scalar or vector); applicable to all interpolation methods. - When both are set, surface node override is ignored (element override wins). - Outlier filtering (processing.interpolation.filter_outliers): - - In `allvoxels` and `allvoxels_weighted`, if enabled and enough voxels are present (>5), outliers are removed using a modified Z-score (median/MAD-based, threshold=3.5) before averaging. + - In `allvoxels` and `allvoxels_scaled`, if enabled and enough voxels are present (>5), outliers are removed using a modified Z-score (median/MAD-based, threshold=3.5) before averaging. - Fallbacks and warnings: - If outlier filtering removes all voxels, the method falls back to the unfiltered mean. diff --git a/src/i2pp/core/configuration_validator/validator.py b/src/i2pp/core/configuration_validator/validator.py index e68d27d..6f247ff 100644 --- a/src/i2pp/core/configuration_validator/validator.py +++ b/src/i2pp/core/configuration_validator/validator.py @@ -110,20 +110,24 @@ def from_dict(d: Dict[str, Any]) -> "Transformation": @dataclass -class NodeWeight: - """Class representing node weight configuration for interpolation.""" +class NodeScaling: + """Class representing node scaling configuration for interpolation.""" - surface: float - interior: float + surface_node_scaling: float = 1.0 + interior_node_scaling: float = 1.0 @staticmethod - def from_dict(d: Optional[Dict[str, Any]]) -> "NodeWeight": - """Creates a NodeWeight instance from a dictionary.""" + def from_dict( + d: Optional[Dict[str, Any]], + ) -> "NodeScaling": + """Creates a NodeScaling instance from a dictionary.""" if d is None: - return NodeWeight(surface=1.0, interior=1.0) - return NodeWeight( - surface=d.get("surface", 1.0), - interior=d.get("interior", 1.0), + return NodeScaling( + surface_node_scaling=1.0, interior_node_scaling=1.0 + ) + return NodeScaling( + surface_node_scaling=d.get("surface", 1.0), + interior_node_scaling=d.get("interior", 1.0), ) @@ -132,7 +136,7 @@ class Interpolation: """Class representing the interpolation configuration.""" method: str - node_weight: NodeWeight + node_scaling_factors: NodeScaling filter_outliers: bool = False set_node_value: Optional[float | list[float]] = field(default=None) set_ele_value: Optional[float | list[float]] = field(default=None) @@ -164,7 +168,9 @@ def from_dict(d: Dict[str, Any]) -> "Interpolation": filter_outliers=d.get("filter_outliers", False), set_node_value=d.get("set_surface_node_value"), set_ele_value=d.get("set_surface_element_value"), - node_weight=NodeWeight.from_dict(d.get("node_weight")), + node_scaling_factors=NodeScaling.from_dict( + d.get("node_scaling_factors") + ), ) diff --git a/src/i2pp/core/discretization_readers/discretization_reader.py b/src/i2pp/core/discretization_readers/discretization_reader.py index 6e37637..0832dfa 100644 --- a/src/i2pp/core/discretization_readers/discretization_reader.py +++ b/src/i2pp/core/discretization_readers/discretization_reader.py @@ -1,12 +1,16 @@ """Import Discretization.""" +from __future__ import annotations + from abc import ABC, abstractmethod from dataclasses import dataclass from pathlib import Path -from typing import Optional +from typing import TYPE_CHECKING, Optional import numpy as np -from i2pp.core.configuration_validator.validator import Processing + +if TYPE_CHECKING: + from i2pp.core.configuration_validator.validator import Processing @dataclass @@ -32,13 +36,15 @@ class Nodes: coords (np.ndarray): An (N, 3) array containing the (x, y, z) world coordinates of each node. ids (np.ndarray): An array of IDs for each node. - weights (Optional[np.ndarray]): Optional per-node weights, - e.g., low weights for surface nodes. + scaling_factors (Optional[np.ndarray]): + Optional per-node scaling_factors that can be used to + adjust the contribution of each node during interpolation, + e.g., low scaling factor for surface nodes. """ coords: np.ndarray ids: np.ndarray - weights: Optional[np.ndarray] = None + scaling_factors: Optional[np.ndarray] = None @dataclass @@ -127,7 +133,10 @@ def __init__(self): @abstractmethod def load_discretization( - self, file_path: Path, options: dict, processing: Processing + self, + file_path: Path, + options: dict, + processing: Processing, ) -> Discretization: """Abstract method to load discretization data from a file path. diff --git a/src/i2pp/core/discretization_readers/fourc_yaml_reader.py b/src/i2pp/core/discretization_readers/fourc_yaml_reader.py index fc85907..9f1421b 100644 --- a/src/i2pp/core/discretization_readers/fourc_yaml_reader.py +++ b/src/i2pp/core/discretization_readers/fourc_yaml_reader.py @@ -1,11 +1,13 @@ """Import 4C.yaml data.""" +from __future__ import annotations + import logging from pathlib import Path +from typing import TYPE_CHECKING import lnmmeshio import numpy as np -from i2pp.core.configuration_validator.validator import Processing from i2pp.core.discretization_readers.discretization_reader import ( Discretization, DiscretizationReader, @@ -16,6 +18,9 @@ from lnmmeshio import Discretization as FourCDiscretization from tqdm import tqdm +if TYPE_CHECKING: + from i2pp.core.configuration_validator.validator import Processing + class FourCYamlReader(DiscretizationReader): """Class for reading and processing finite element models from .4C.yaml @@ -99,6 +104,12 @@ def load_discretization( logging.info("Importing discretization data") + if processing is None: + raise ValueError( + "Processing configuration is required" + "for loading the discretization." + ) + raw_dis = lnmmeshio.read(str(file_path)) raw_dis.compute_ids(zero_based=True) @@ -108,17 +119,19 @@ def load_discretization( raw_dis, np.array(options["material_ids"]) ) - interior_weight = processing.interpolation.node_weight.interior - surface_weight = processing.interpolation.node_weight.surface + scaling_factors = processing.interpolation.node_scaling_factors + + interior_node_scaling = scaling_factors.interior_node_scaling + surface_node_scaling = scaling_factors.surface_node_scaling nodes_coords = [] node_ids = [] - nodes_weights = [] + nodes_scaling = [] for node in raw_dis.nodes: nodes_coords.append(node.coords) node_ids.append(node.id) - nodes_weights.append(interior_weight) + nodes_scaling.append(interior_node_scaling) elements = [] @@ -139,7 +152,7 @@ def load_discretization( for node in surf.nodes: surf_node_ids.append(node.id) position = node_id_to_idx[node.id] - nodes_weights[position] = surface_weight + nodes_scaling[position] = surface_node_scaling surfaces.append( Surface(node_ids=np.array(surf_node_ids), id=surf.id) @@ -149,7 +162,7 @@ def load_discretization( nodes=Nodes( coords=np.array(nodes_coords), ids=np.array(node_ids), - weights=np.array(nodes_weights), + scaling_factors=np.array(nodes_scaling), ), elements=elements, surfaces=surfaces, diff --git a/src/i2pp/core/discretization_readers/mesh_reader.py b/src/i2pp/core/discretization_readers/mesh_reader.py index 53f448f..04542da 100644 --- a/src/i2pp/core/discretization_readers/mesh_reader.py +++ b/src/i2pp/core/discretization_readers/mesh_reader.py @@ -1,10 +1,12 @@ """Import Mesh data.""" +from __future__ import annotations + import logging from pathlib import Path +from typing import TYPE_CHECKING import trimesh -from i2pp.core.configuration_validator.validator import Processing from i2pp.core.discretization_readers.discretization_reader import ( Discretization, DiscretizationReader, @@ -12,6 +14,9 @@ Nodes, ) +if TYPE_CHECKING: + from i2pp.core.configuration_validator.validator import Processing + class MeshReader(DiscretizationReader): """Class for reading and processing finite element models from .mesh files. diff --git a/src/i2pp/core/interpolators/interpolator_all_voxel.py b/src/i2pp/core/interpolators/interpolator_all_voxel.py index 47ce82d..2c0c7b0 100644 --- a/src/i2pp/core/interpolators/interpolator_all_voxel.py +++ b/src/i2pp/core/interpolators/interpolator_all_voxel.py @@ -17,12 +17,12 @@ class InterpolatorAllVoxel(Interpolator): """Interpolator for mapping 3D image data to finite element mesh elements. - This class supports both unweighted and node-weighted voxel mean + This class supports both unscaled and node-scaled voxel mean calculations, controlled via the 'mode' parameter. It assigns pixel values from 3D image data to finite element mesh elements by computing the mean of all voxels within each element. This functionality is used when the `interpolation_method` is set to - "allvoxels" or "allvoxels_weighted". + "allvoxels" or "allvoxels_scaled". """ def __init__( @@ -34,7 +34,7 @@ def __init__( ): super().__init__(*args, **kwargs) self._filter_outliers_enabled = filter_outliers - self._mode = mode # "allvoxels" or "allvoxels_weighted" + self._mode = mode # "allvoxels" or "allvoxels_scaled" def _search_bounding_box( self, grid_coords: GridCoords, element_grid_coords: np.ndarray @@ -144,15 +144,15 @@ def _weighted_voxel_mean( voxels_phys: np.ndarray, values: np.ndarray, element_node_phys: np.ndarray, - node_weights: np.ndarray, + node_scaling_factors: np.ndarray, ) -> float: """Calculate a node-weighted mean of voxel values, with optional outlier filtering. Computes the weighted mean of voxel values within an element, where the weight is determined based on the distances between the voxel - and the element's nodes, scaled by the node weights. Outliers can be - excluded using the Modified Z-Score method. + and the element's nodes, scaled by the node scaling factors. + Outliers can be excluded using the Modified Z-Score method. Args: voxels_phys (np.ndarray): @@ -161,8 +161,8 @@ def _weighted_voxel_mean( Corresponding voxel values (N_voxels x ...). element_node_phys (np.ndarray): Physical coordinates of the element's nodes (N_nodes x 3). - node_weights (np.ndarray): - Weights assigned to the nodes (N_nodes,). + node_scaling_factors (np.ndarray): + Scaling factors assigned to the nodes (N_nodes,). Returns: float: The weighted mean of the voxel values. @@ -175,7 +175,9 @@ def _weighted_voxel_mean( ) distances = np.maximum(distances, 1e-10) - voxel_weights = np.sum(node_weights[:, np.newaxis] / distances, axis=0) + voxel_weights = np.sum( + node_scaling_factors[:, np.newaxis] / distances, axis=0 + ) if len(values) > 5: mask = self._filter_outliers_modified_zscore(values) @@ -225,7 +227,7 @@ def _get_data_of_element( self, element_node_grid_coords: np.ndarray, image_data: ImageData, - node_weights_current: np.ndarray | None = None, + node_scaling_factors_current: np.ndarray | None = None, ) -> np.ndarray: """Computes the representative pixel value for a given element based on its nodes in grid coordinates. @@ -244,8 +246,8 @@ def _get_data_of_element( the element's nodes. image_data (ImageData): Image data containing voxel coordinates and pixel values. - node_weights_current (np.ndarray | None): - Weights assigned to the nodes. + node_scaling_factors_current (np.ndarray | None): + Scaling factors assigned to the nodes. Returns: np.ndarray: The mean pixel value of all voxels inside the element. @@ -304,9 +306,9 @@ def _get_data_of_element( mean_val = np.mean(data, axis=0) return self._format_output_value(mean_val, image_data) - if self._mode == "allvoxels_weighted": + if self._mode == "allvoxels_scaled": element_node_phys = element_node_grid_coords - if node_weights_current is None: + if node_scaling_factors_current is None: weighted = np.mean(data, axis=0) else: if self._filter_outliers_enabled: @@ -314,7 +316,7 @@ def _get_data_of_element( voxels_phys_np, data, element_node_phys, - node_weights_current, + node_scaling_factors_current, ) else: # Non-filtered weighted path with zero-sum guard @@ -325,7 +327,8 @@ def _get_data_of_element( ) distances = np.maximum(distances, 1e-10) voxel_weights = np.sum( - node_weights_current[:, np.newaxis] / distances, + node_scaling_factors_current[:, np.newaxis] + / distances, axis=0, ) @@ -389,16 +392,21 @@ def compute_element_data( desc="Element values", ): element_node_grid_coords = node_grid_coords[node_positions[i]] - # Cast weights to ndarray when present to satisfy type checker - node_weights_current = None - if getattr(dis.nodes, "weights", None) is not None: - weights_nd = cast(np.ndarray, dis.nodes.weights) - node_weights_current = weights_nd[node_positions[i]] + # Cast scaling factors to ndarray when present + # to satisfy type checker + node_scaling_factors_current = None + if getattr(dis.nodes, "scaling_factors", None) is not None: + scaling_factors_nd = cast( + np.ndarray, dis.nodes.scaling_factors + ) + node_scaling_factors_current = scaling_factors_nd[ + node_positions[i] + ] ele.data = self._get_data_of_element( element_node_grid_coords, image_data, - node_weights_current=node_weights_current, + node_scaling_factors_current=node_scaling_factors_current, ) if np.all(np.isnan(ele.data)): diff --git a/src/i2pp/core/interpolators/interpolator_nodes.py b/src/i2pp/core/interpolators/interpolator_nodes.py index 093fe61..611395a 100644 --- a/src/i2pp/core/interpolators/interpolator_nodes.py +++ b/src/i2pp/core/interpolators/interpolator_nodes.py @@ -23,7 +23,7 @@ class InterpolatorNodes(Interpolator): "nodes". """ - # Add mode to control whether to use node weights + # Add mode to control whether to use node scaling factors def __init__( self, *args, @@ -33,7 +33,7 @@ def __init__( ): """Initializes the InterpolatorNodes.""" super().__init__() - self._mode = mode # "nodes" or "nodes_weighted" + self._mode = mode # "nodes" or "nodes_scaled" self.set_node_value = surf_node_val # Helpers for readability and error handling @@ -77,43 +77,47 @@ def _compute_unweighted_mean( def _compute_weighted_mean( self, ele_nodes: np.ndarray, - weights_current: np.ndarray, + scaling_factors_current: np.ndarray, num_values: int, ) -> np.ndarray: """Computes the weighted mean of pixel values for an element's nodes. This method calculates the weighted average of pixel values from an element's nodes. It is used when a more nuanced contribution of each - node is desired, based on pre-assigned weights. Nodes with NaN values - are excluded from the calculation. The weights of the valid nodes are - normalized to sum to 1 before computing the average to prevent scaling - issues. + node is desired, based on pre-assigned scaling factor. + Nodes with NaN values are excluded from the calculation. + The factors of the valid nodes are normalized to sum to 1 (weighted) + before computing the average. Args: ele_nodes (np.ndarray): An array of interpolated pixel values for the nodes of a single element. Shape can be (num_nodes,) for single-channel data or (num_nodes, num_values) for multi- channel data. - weights_current (np.ndarray): An array of weights corresponding to + scaling_factors_current (np.ndarray): + An array of scaling factors corresponding to each node in `ele_nodes`. Shape is (num_nodes,). num_values (int): The number of values per pixel (e.g., 1 for grayscale, 3 for RGB). Returns: - np.ndarray: An array containing the weighted mean pixel value(s) + np.ndarray: An array containing the scaled mean pixel value(s) for the element. The shape is (num_values,). Raises: - ValueError: If `weights_current` is None or if the shape of - `ele_nodes` and `weights_current` are incompatible. + ValueError: If `scaling_factors_current` is None or if the shape of + `ele_nodes` and `scaling_factors_current` are incompatible. """ - if weights_current is None: - raise ValueError("Node weights are required for weighted mode.") - if ele_nodes.shape[0] != weights_current.shape[0]: + if scaling_factors_current is None: + raise ValueError( + "Node scaling factors are required for scaled mode." + ) + if ele_nodes.shape[0] != scaling_factors_current.shape[0]: raise ValueError( "Incompatible shapes:" f" ele_nodes has {ele_nodes.shape[0]} nodes, " - f"weights_current has {weights_current.shape[0]}." + "scaling_factors_current has" + f"{scaling_factors_current.shape[0]}." ) # Mask out nodes that contain any NaN across channels @@ -126,16 +130,16 @@ def _compute_weighted_mean( return np.full(num_values, np.nan) vals = ele_nodes[nan_mask] - w = weights_current[nan_mask] + sf = scaling_factors_current[nan_mask] - wsum = float(np.sum(w)) - if wsum == 0: + sf_sum = float(np.sum(sf)) + if sf_sum == 0: if num_values > 1: return np.mean(vals, axis=0) else: return np.array([np.mean(vals)]) - w = w / wsum # Normalize weights + w = sf / sf_sum # Normalize scaling factors if num_values == 1: val = np.average(vals.reshape(-1), weights=w) @@ -236,13 +240,13 @@ def compute_element_data( ] ) - # Only prepare node weights if we're in weighted mode - node_weights = None + # Only prepare node scaling factors if we're in scaled mode + node_scaling_factors = None if ( - self._mode == "nodes_weighted" - and getattr(dis.nodes, "weights", None) is not None + self._mode == "nodes_scaled" + and getattr(dis.nodes, "scaling_factors", None) is not None ): - node_weights = cast(np.ndarray, dis.nodes.weights) + node_scaling_factors = cast(np.ndarray, dis.nodes.scaling_factors) for i, ele in tqdm( enumerate(dis.elements), @@ -252,10 +256,12 @@ def compute_element_data( ele_nodes = node_values[node_positions[i]] num_values = image_data.pixel_type.num_values - if node_weights is not None: - weights_current = node_weights[node_positions[i]] + if node_scaling_factors is not None: + scaling_factors_current = node_scaling_factors[ + node_positions[i] + ] ele.data = self._compute_weighted_mean( - ele_nodes, weights_current, num_values + ele_nodes, scaling_factors_current, num_values ) else: ele.data = self._compute_unweighted_mean(ele_nodes, num_values) diff --git a/src/i2pp/core/interpolators/interpolator_types.py b/src/i2pp/core/interpolators/interpolator_types.py index 78a7bc2..67fd504 100644 --- a/src/i2pp/core/interpolators/interpolator_types.py +++ b/src/i2pp/core/interpolators/interpolator_types.py @@ -20,23 +20,23 @@ class InterpolationType(Enum): Attributes: NODES (str): Represents the interpolation method where the pixel value is averaged over the nodes of the element. - NODES_WEIGHTED (str): Represents the interpolation method where the - pixel value is weighted over the nodes of the element. + NODES_SCALED (str): Represents the interpolation method where the + pixel value is scaled over the nodes of the element. CENTER (str): Represents the interpolation method where the pixel value is based on the center of the element. ALLVOXELS (str): Represents the interpolation method where the pixel value is averaged over all voxels inside the element. - ALLVOXELS_WEIGHTED (str): Represents the interpolation method where the - pixel value is weighted over all voxels inside the element. + ALLVOXELS_SCALED (str): Represents the interpolation method where the + pixel value is scaled over all voxels inside the element. Use create_interpolator to obtain a configured interpolator instance. """ NODES = "nodes" - NODES_WEIGHTED = "nodes_weighted" + NODES_SCALED = "nodes_scaled" CENTER = "elementcenter" ALLVOXELS = "allvoxels" - ALLVOXELS_WEIGHTED = "allvoxels_weighted" + ALLVOXELS_SCALED = "allvoxels_scaled" def create_interpolator( self, @@ -69,17 +69,17 @@ def create_interpolator( return InterpolatorAllVoxel( mode="allvoxels", filter_outliers=filter_outliers ) - if self == InterpolationType.ALLVOXELS_WEIGHTED: + if self == InterpolationType.ALLVOXELS_SCALED: return InterpolatorAllVoxel( - mode="allvoxels_weighted", filter_outliers=filter_outliers + mode="allvoxels_scaled", filter_outliers=filter_outliers ) if self == InterpolationType.NODES: return InterpolatorNodes( mode="nodes", surf_node_val=set_node_value ) - if self == InterpolationType.NODES_WEIGHTED: + if self == InterpolationType.NODES_SCALED: return InterpolatorNodes( - mode="nodes_weighted", surf_node_val=set_node_value + mode="nodes_scaled", surf_node_val=set_node_value ) if self == InterpolationType.CENTER: return InterpolatorCenter() diff --git a/tests/unittests/i2pp/core/discretization_readers/test_fourc_yaml_reader.py b/tests/unittests/i2pp/core/discretization_readers/test_fourc_yaml_reader.py index f3852c8..7ee338a 100644 --- a/tests/unittests/i2pp/core/discretization_readers/test_fourc_yaml_reader.py +++ b/tests/unittests/i2pp/core/discretization_readers/test_fourc_yaml_reader.py @@ -87,8 +87,11 @@ def test_load_discretization_fourc_yaml_without_filter(tmp_path: Path) -> None: # test_config test_options = {"material_ids": None} mock_processing = MagicMock() - mock_processing.interpolation.interior_node_weight = 1.0 - mock_processing.interpolation.surface_node_weight = 0.0 + mock_scaling_factors = ( + mock_processing.interpolation.node_scaling_factors + ) + mock_scaling_factors.interior_node_scaling = 1.0 + mock_scaling_factors.surface_node_scaling = 0.0 dis_loaded = test_dis.load_discretization( Path(test_path), test_options, mock_processing diff --git a/tests/unittests/i2pp/core/discretization_readers/test_mesh_reader.py b/tests/unittests/i2pp/core/discretization_readers/test_mesh_reader.py index cb61afc..a9d3028 100644 --- a/tests/unittests/i2pp/core/discretization_readers/test_mesh_reader.py +++ b/tests/unittests/i2pp/core/discretization_readers/test_mesh_reader.py @@ -5,7 +5,7 @@ from i2pp.core.configuration_validator.validator import ( Interpolation, - NodeWeight, + NodeScaling, Processing, Transformation, ) @@ -25,7 +25,9 @@ def test_load_discretization_mesh(tmp_path: Path) -> None: interpolation=Interpolation( method="nodes", filter_outliers=False, - node_weight=NodeWeight(interior=0.5, surface=0.5), + node_scaling_factors=NodeScaling( + interior_node_scaling=0.5, surface_node_scaling=0.5 + ), ), transformation=Transformation( user_script=Path(""), diff --git a/tests/unittests/i2pp/core/interpolators/test_interpolator_all_voxel.py b/tests/unittests/i2pp/core/interpolators/test_interpolator_all_voxel.py index 2603164..fe8fdc8 100644 --- a/tests/unittests/i2pp/core/interpolators/test_interpolator_all_voxel.py +++ b/tests/unittests/i2pp/core/interpolators/test_interpolator_all_voxel.py @@ -285,15 +285,15 @@ def test_get_data_of_element_element_low_resolution_image(): def _make_simple_discretization( - node_coords, elements_node_ids, node_weights=None + node_coords, elements_node_ids, node_scaling_factors=None ): """Helper to create a minimal Discretization with nodes/elements and - weights.""" + scaling factors.""" node_ids = np.arange(len(node_coords)) nodes = Nodes(ids=node_ids, coords=np.asarray(node_coords)) - if node_weights is None: - node_weights = np.ones(len(node_coords), dtype=float) - setattr(nodes, "weights", np.asarray(node_weights)) + if node_scaling_factors is None: + node_scaling_factors = np.ones(len(node_coords), dtype=float) + setattr(nodes, "scaling_factors", np.asarray(node_scaling_factors)) # Create Elements with required 'id' field; # center_coords/data left as defaults @@ -308,9 +308,9 @@ def _make_simple_discretization( return Discretization(nodes=nodes, elements=elements, surfaces=surfaces) -def test_unweighted_vs_weighted_mode_mean_differs_when_weights_bias(): - """Weighted mode should differ from unweighted mean when node weights bias - proximity.""" +def test_unscaled_vs_scaled_mode_mean_differs_when_scaling_factors_bias(): + """Scaled mode should differ from unscaled mean when node scaling factors + bias proximity.""" # Image: small grid 3x3x3 with deterministic values pixel_data = np.arange(27, dtype=float).reshape((3, 3, 3)) grid_coords = GridCoords( @@ -338,36 +338,35 @@ def test_unweighted_vs_weighted_mode_mean_differs_when_weights_bias(): dis = _make_simple_discretization( node_coords=ele_nodes, elements_node_ids=[np.arange(8)], - node_weights=np.array( + node_scaling_factors=np.array( [10, 10, 10, 10, 1, 1, 1, 1], dtype=float ), # bias towards lower z nodes ) - # Unweighted - interpolator_unweighted = InterpolatorAllVoxel( + # Unscaled + interpolator_unscaled = InterpolatorAllVoxel( mode="allvoxels", filter_outliers=False ) - elems_unweighted = interpolator_unweighted.compute_element_data( + elems_unscaled = interpolator_unscaled.compute_element_data( dis, image_data ) - val_unweighted = elems_unweighted[0].data + val_unscaled = elems_unscaled[0].data - # Weighted - interpolator_weighted = InterpolatorAllVoxel( - mode="allvoxels_weighted", filter_outliers=False + # Scaled + interpolator_scaled = InterpolatorAllVoxel( + mode="allvoxels_scaled", filter_outliers=False ) - elems_weighted = interpolator_weighted.compute_element_data( - dis, image_data - ) - val_weighted = elems_weighted[0].data + elems_scaled = interpolator_scaled.compute_element_data(dis, image_data) + val_scaled = elems_scaled[0].data - assert np.isfinite(val_unweighted).all() - assert np.isfinite(val_weighted).all() - # Expect a difference due to weight bias toward voxels near lower-z nodes - assert not np.allclose(val_unweighted, val_weighted) + assert np.isfinite(val_unscaled).all() + assert np.isfinite(val_scaled).all() + # Expect a difference due to scaling factor + # bias toward voxels near lower-z nodes + assert not np.allclose(val_unscaled, val_scaled) -def test_filter_outliers_reduces_extreme_values_in_unweighted_mode(): +def test_filter_outliers_reduces_extreme_values_in_unscaled_mode(): """Outlier filtering should reduce influence of extreme voxel values.""" # Construct image with an element that will include voxels; # insert an extreme outlier @@ -412,8 +411,8 @@ def test_filter_outliers_reduces_extreme_values_in_unweighted_mode(): assert np.isclose(val_filter, 1.0, rtol=1e-3) -def test_filter_outliers_applies_in_weighted_mode(): - """Outlier filtering should apply in weighted mode as well.""" +def test_filter_outliers_applies_in_scaled_mode(): + """Outlier filtering should apply in scaled mode as well.""" N = 5 pixel_data = np.ones((N, N, N), dtype=float) pixel_data[2, 2, 2] = 1e6 @@ -437,22 +436,22 @@ def test_filter_outliers_applies_in_weighted_mode(): ] ) - # Node weights biased toward corners near + # Node scaling factors biased toward corners near # the outlier but filtering should mitigate - weights = np.array([5, 5, 5, 5, 5, 5, 5, 5], dtype=float) + scaling_factors = np.array([5, 5, 5, 5, 5, 5, 5, 5], dtype=float) interp_w_no_filter = InterpolatorAllVoxel( - mode="allvoxels_weighted", filter_outliers=False + mode="allvoxels_scaled", filter_outliers=False ) val_w_no_filter = interp_w_no_filter._get_data_of_element( - element, image_data, node_weights_current=weights + element, image_data, node_scaling_factors_current=scaling_factors ) interp_w_filter = InterpolatorAllVoxel( - mode="allvoxels_weighted", filter_outliers=True + mode="allvoxels_scaled", filter_outliers=True ) val_w_filter = interp_w_filter._get_data_of_element( - element, image_data, node_weights_current=weights + element, image_data, node_scaling_factors_current=scaling_factors ) # Unfiltered should be larger than ~1 due to outlier influence @@ -463,8 +462,8 @@ def test_filter_outliers_applies_in_weighted_mode(): assert np.all(val_w_filter < val_w_no_filter) -def test_unknown_mode_defaults_to_unweighted_mean(): - """Unknown mode should fall back to unweighted mean.""" +def test_unknown_mode_defaults_to_unscaled_mean(): + """Unknown mode should fall back to unscaled mean.""" element = np.array( [ [0, 2, 1], @@ -488,7 +487,7 @@ def test_unknown_mode_defaults_to_unweighted_mean(): pixel_data, grid_coords, np.eye(3), np.zeros(3), PixelValueType.CT ) - # Expected unweighted mean of voxels inside the element + # Expected unscaled mean of voxels inside the element data_inside = [11, 12, 16, 17, 36, 37, 41, 42] expected_mean = np.mean(data_inside) @@ -501,9 +500,9 @@ def test_unknown_mode_defaults_to_unweighted_mean(): assert np.isclose(result_val, expected_mean) -def test_weighted_small_voxel_count_uses_simple_weighted_average(): +def test_scaled_small_voxel_count_uses_simple_scaled_average(): """When voxel count is small (<=5), _weighted_voxel_mean should return - simple weighted average.""" + simple scaled average.""" # Grid 2x2x2 grid_coords = GridCoords( slice=np.array([0.0, 1.0]), @@ -528,24 +527,24 @@ def test_weighted_small_voxel_count_uses_simple_weighted_average(): [0.0, 0.1, 0.5], ] ) - # One element uses the 4 nodes, with explicit node weights + # One element uses the 4 nodes, with explicit node scaling factors node_ids = np.arange(len(element_nodes)) nodes = Nodes( ids=node_ids, coords=element_nodes, - weights=np.array([2.0, 2.0, 1.0, 1.0]), + scaling_factors=np.array([2.0, 2.0, 1.0, 1.0]), ) elements = [DisElement(node_ids=node_ids, id=0)] surfaces = [Surface(node_ids=np.array([], dtype=int), id=0)] dis = Discretization(nodes=nodes, elements=elements, surfaces=surfaces) interpolator = InterpolatorAllVoxel( - mode="allvoxels_weighted", filter_outliers=True + mode="allvoxels_scaled", filter_outliers=True ) elems = interpolator.compute_element_data(dis, image_data) result = float(np.asarray(elems[0].data).mean()) - # Manually compute included voxels and expected weighted average + # Manually compute included voxels and expected scaled average hull = ConvexHull(element_nodes) included_points = [] included_vals = [] @@ -561,23 +560,25 @@ def test_weighted_small_voxel_count_uses_simple_weighted_average(): voxels_phys = np.asarray(included_points) values = np.asarray(included_vals) - # Expect small voxel count (<=5) to use simple weighted average + # Expect small voxel count (<=5) to use simple scaled average assert len(values) <= 5 - # Recompute voxel weights exactly as in implementation + # Recompute voxel scaling factors exactly as in implementation distances = np.linalg.norm( element_nodes[:, np.newaxis, :] - voxels_phys[np.newaxis, :, :], axis=2 ) distances = np.maximum(distances, 1e-10) - voxel_weights = np.sum(nodes.weights[:, np.newaxis] / distances, axis=0) + voxel_scaling_factors = np.sum( + nodes.scaling_factors[:, np.newaxis] / distances, axis=0 + ) - expected = np.average(values, weights=voxel_weights, axis=0) + expected = np.average(values, weights=voxel_scaling_factors, axis=0) assert np.isclose(result, expected) -def test_weighted_mode_no_filter_falls_back_when_weights_sum_zero(): - """Weighted/no-filter path should fall back to unweighted mean when weights - sum to zero.""" +def test_scaled_mode_no_filter_falls_back_when_scaling_factors_sum_zero(): + """Scaled/no-filter path should fall back to unscaled mean when scaling + factors sum to zero.""" # Simple 3x3x3 grid with deterministic values pixel_data = np.arange(27, dtype=float).reshape((3, 3, 3)) grid_coords = GridCoords( @@ -602,37 +603,37 @@ def test_weighted_mode_no_filter_falls_back_when_weights_sum_zero(): ) # Discretization: one element using these nodes; - # set all node weights to zero + # set all node scaling factors to zero dis = _make_simple_discretization( node_coords=element_nodes, elements_node_ids=[np.arange(8)], - node_weights=np.zeros(8, dtype=float), + node_scaling_factors=np.zeros(8, dtype=float), ) - # Compute with weighted mode, no outlier filtering - interp_weighted = InterpolatorAllVoxel( - mode="allvoxels_weighted", filter_outliers=False + # Compute with scaled mode, no outlier filtering + interp_scaled = InterpolatorAllVoxel( + mode="allvoxels_scaled", filter_outliers=False ) - elems = interp_weighted.compute_element_data(dis, image_data) - val_weighted = float(np.asarray(elems[0].data).mean()) + elems = interp_scaled.compute_element_data(dis, image_data) + val_scaled = float(np.asarray(elems[0].data).mean()) - # Compute unweighted reference - interp_unweighted = InterpolatorAllVoxel( + # Compute unscaled reference + interp_unscaled = InterpolatorAllVoxel( mode="allvoxels", filter_outliers=False ) - val_unweighted = float( + val_unscaled = float( np.asarray( - interp_unweighted.compute_element_data(dis, image_data)[0].data + interp_unscaled.compute_element_data(dis, image_data)[0].data ).mean() ) - # Should match unweighted mean (fallback) and not raise ZeroDivisionError - assert np.isclose(val_weighted, val_unweighted) + # Should match unscaled mean (fallback) and not raise ZeroDivisionError + assert np.isclose(val_scaled, val_unscaled) -def test_weighted_mode_with_filter_falls_back_when_weights_sum_zero(): - """Weighted/filtered path should fall back to unweighted mean when weights - sum to zero.""" +def test_scaled_mode_with_filter_falls_back_when_scaling_factors_sum_zero(): + """Scaled/filtered path should fall back to unscaled mean when scaling + factors sum to zero.""" # Grid with ones and an outlier to exercise filtering path N = 5 pixel_data = np.ones((N, N, N), dtype=float) @@ -658,30 +659,28 @@ def test_weighted_mode_with_filter_falls_back_when_weights_sum_zero(): ] ) - # Zero node weights to trigger zero-sum voxel weights + # Zero node scaling factors to trigger zero-sum voxel scaling factors dis = _make_simple_discretization( node_coords=element_nodes, elements_node_ids=[np.arange(8)], - node_weights=np.zeros(8, dtype=float), + node_scaling_factors=np.zeros(8, dtype=float), ) - # Weighted with filtering; should not raise - # and should approximate unweighted filtered mean - interp_weighted_filter = InterpolatorAllVoxel( - mode="allvoxels_weighted", filter_outliers=True - ) - elems_weighted = interp_weighted_filter.compute_element_data( - dis, image_data + # Scaled with filtering; should not raise + # and should approximate unscaled filtered mean + interp_scaled_filter = InterpolatorAllVoxel( + mode="allvoxels_scaled", filter_outliers=True ) - val_weighted = float(np.asarray(elems_weighted[0].data).mean()) + elems_scaled = interp_scaled_filter.compute_element_data(dis, image_data) + val_scaled = float(np.asarray(elems_scaled[0].data).mean()) - # Unweighted with filtering as reference - interp_unweighted_filter = InterpolatorAllVoxel( + # Unscaled with filtering as reference + interp_unscaled_filter = InterpolatorAllVoxel( mode="allvoxels", filter_outliers=True ) - elems_unweighted = interp_unweighted_filter.compute_element_data( + elems_unscaled = interp_unscaled_filter.compute_element_data( dis, image_data ) - val_unweighted = float(np.asarray(elems_unweighted[0].data).mean()) + val_unscaled = float(np.asarray(elems_unscaled[0].data).mean()) - assert np.isclose(val_weighted, val_unweighted, rtol=1e-6) + assert np.isclose(val_scaled, val_unscaled, rtol=1e-6) diff --git a/tests/unittests/i2pp/core/interpolators/test_interpolator_nodes.py b/tests/unittests/i2pp/core/interpolators/test_interpolator_nodes.py index 87cbace..cc621df 100644 --- a/tests/unittests/i2pp/core/interpolators/test_interpolator_nodes.py +++ b/tests/unittests/i2pp/core/interpolators/test_interpolator_nodes.py @@ -87,8 +87,8 @@ def _build_simple_setup_rgb(): """Build simple discretization and image data for RGB tests.""" node_coords = np.array([[0, 0, 0], [1, 0, 0]]) nodes = Nodes(coords=node_coords, ids=[0, 1]) - # Provide weights for nodes - nodes.weights = np.array([0.2, 0.8]) + # Provide scaling factors for nodes + nodes.scaling_factors = np.array([0.2, 0.8]) elements = [Element(node_ids=[0, 1], id=0)] dis = Discretization(nodes=nodes, elements=elements, surfaces=[]) @@ -105,8 +105,8 @@ def _build_simple_setup_rgb(): return dis, image_data -def test_nodes_unweighted_rgb_ignores_weights(): - """Unweighted mode computes nanmean and ignores weights for RGB.""" +def test_nodes_unscaled_rgb_ignores_scaling_factors(): + """Unscaled mode computes nanmean and ignores scaling factors for RGB.""" dis, image_data = _build_simple_setup_rgb() interpolator = InterpolatorNodes(mode="nodes") @@ -132,10 +132,10 @@ def test_nodes_unweighted_rgb_ignores_weights(): assert np.allclose(result[0].data, expected) -def test_nodes_weighted_rgb(): - """Weighted mode uses node weights for vector-valued (RGB) data.""" +def test_nodes_scaled_rgb(): + """Scaled mode uses node scaling factors for vector-valued (RGB) data.""" dis, image_data = _build_simple_setup_rgb() - interpolator = InterpolatorNodes(mode="nodes_weighted") + interpolator = InterpolatorNodes(mode="nodes_scaled") node_values = np.array([[100, 150, 200], [200, 250, 300]]) with ( @@ -159,10 +159,10 @@ def test_nodes_weighted_rgb(): assert np.allclose(result[0].data, expected) -def test_nodes_weighted_rgb_handles_nans(): - """Weighted mode handles NaNs by masking them out for RGB.""" +def test_nodes_scaled_rgb_handles_nans(): + """Scaled mode handles NaNs by masking them out for RGB.""" dis, image_data = _build_simple_setup_rgb() - interpolator = InterpolatorNodes(mode="nodes_weighted") + interpolator = InterpolatorNodes(mode="nodes_scaled") node_values = np.array([[np.nan, np.nan, np.nan], [200, 250, 300]]) with ( diff --git a/tests/unittests/i2pp/core/interpolators/test_interpolator_types.py b/tests/unittests/i2pp/core/interpolators/test_interpolator_types.py index d83befa..e910b94 100644 --- a/tests/unittests/i2pp/core/interpolators/test_interpolator_types.py +++ b/tests/unittests/i2pp/core/interpolators/test_interpolator_types.py @@ -11,10 +11,10 @@ "enum_value, expected_class", [ (InterpolationType.NODES, InterpolatorNodes), - (InterpolationType.NODES_WEIGHTED, InterpolatorNodes), + (InterpolationType.NODES_SCALED, InterpolatorNodes), (InterpolationType.CENTER, InterpolatorCenter), (InterpolationType.ALLVOXELS, InterpolatorAllVoxel), - (InterpolationType.ALLVOXELS_WEIGHTED, InterpolatorAllVoxel), + (InterpolationType.ALLVOXELS_SCALED, InterpolatorAllVoxel), ], ) def test_create_interpolator_returns_expected_class( @@ -28,17 +28,17 @@ def test_create_interpolator_returns_expected_class( def test_enum_values(): """Enum values are correctly defined.""" assert InterpolationType.NODES.value == "nodes" - assert InterpolationType.NODES_WEIGHTED.value == "nodes_weighted" + assert InterpolationType.NODES_SCALED.value == "nodes_scaled" assert InterpolationType.CENTER.value == "elementcenter" assert InterpolationType.ALLVOXELS.value == "allvoxels" - assert InterpolationType.ALLVOXELS_WEIGHTED.value == "allvoxels_weighted" + assert InterpolationType.ALLVOXELS_SCALED.value == "allvoxels_scaled" @pytest.mark.parametrize( "enum_value, expected_mode", [ (InterpolationType.ALLVOXELS, "allvoxels"), - (InterpolationType.ALLVOXELS_WEIGHTED, "allvoxels_weighted"), + (InterpolationType.ALLVOXELS_SCALED, "allvoxels_scaled"), ], ) def test_create_interpolator_configures_mode(enum_value, expected_mode): diff --git a/tests/unittests/i2pp/core/test_interpolate_element_data.py b/tests/unittests/i2pp/core/test_interpolate_element_data.py index c9ec09e..1f6390b 100644 --- a/tests/unittests/i2pp/core/test_interpolate_element_data.py +++ b/tests/unittests/i2pp/core/test_interpolate_element_data.py @@ -90,7 +90,7 @@ def test_boundary_elements_receive_fixed_scalar_value(): filter_outliers=False, set_node_value=None, set_ele_value=999.0, - node_weight=None, + node_scaling_factors=None, ) elements = interpolate_image_to_discretization(dis, image, interp_cfg) @@ -112,7 +112,7 @@ def test_boundary_elements_receive_fixed_rgb_value_vector(): filter_outliers=False, set_node_value=None, set_ele_value=fixed_rgb, - node_weight=None, + node_scaling_factors=None, ) elements = interpolate_image_to_discretization(dis, image, interp_cfg) @@ -135,7 +135,7 @@ def test_boundary_elements_receive_fixed_value_without_prior_data(): filter_outliers=False, set_node_value=None, set_ele_value=123.0, - node_weight=None, + node_scaling_factors=None, ) # Perform interpolation @@ -152,7 +152,7 @@ def test_boundary_elements_receive_fixed_value_without_prior_data(): filter_outliers=False, set_node_value=None, set_ele_value=[7, 8, 9], - node_weight=None, + node_scaling_factors=None, ) elements_vec = interpolate_image_to_discretization( dis, image_rgb, interp_cfg_vec diff --git a/tests/unittests/i2pp/core/test_run.py b/tests/unittests/i2pp/core/test_run.py index 735fc3c..6594497 100644 --- a/tests/unittests/i2pp/core/test_run.py +++ b/tests/unittests/i2pp/core/test_run.py @@ -150,8 +150,8 @@ def maximal_valid_config(tmp_path): }, "processing": { "interpolation": { - "method": "nodes_weighted", - "node_weight": {"interior": 0.8, "surface": 0.2}, + "method": "nodes_scaled", + "node_scaling_factors": {"interior": 0.8, "surface": 0.2}, "filter_outliers": True, "set_surface_node_value": 100.0, }, From 417ea97f5391c37103287aed680f20b9ca129579 Mon Sep 17 00:00:00 2001 From: Saskia Gloy Date: Tue, 24 Feb 2026 19:39:07 +0100 Subject: [PATCH 15/15] Add inverse distance weighting power parameter to to allvoxel_scaled interpolation method --- README.md | 19 +- .../core/configuration_validator/validator.py | 18 +- src/i2pp/core/interpolate_element_data.py | 1 + .../interpolators/interpolator_all_voxel.py | 129 +++++--- .../core/interpolators/interpolator_types.py | 7 +- .../configuration_validator/test_validator.py | 18 +- .../test_interpolator_all_voxel.py | 284 ++++++++++++++++++ tests/unittests/i2pp/core/test_run.py | 3 +- 8 files changed, 417 insertions(+), 62 deletions(-) diff --git a/README.md b/README.md index fac6eec..96ca01c 100644 --- a/README.md +++ b/README.md @@ -98,19 +98,18 @@ pdoc --html --output-dir docs src/i2pp ### Interpolation and value settings -- Interpolation methods (processing.interpolation.method): - - nodes: Interpolates values at the element’s nodes and assigns the element mean (ignoring NaN nodes). Fast and robust; respects node sampling. - - nodes_scaled: Like `nodes`, but computes a scaled mean using node-specific scaling factors (`dis.nodes.scaling_factors`), which are set via `processing.interpolation.node_scaling_factors.surface` and `processing.interpolation.node_scaling_factors.interior`. This adjusts the influence of specific nodes. - - elementcenter: Interpolates at each element centroid and assigns that value. - - allvoxels: Collects all voxels whose grid coordinates lie inside the convex hull of the element nodes; assigns the mean value; optionally filters outliers. - - allvoxels_scaled: Computes a voxel-weighted mean where voxel weights derive from node scaling factors and inverse node-to-voxel distances; optionally filters outliers. +- Interpolation methods (`processing.interpolation.method`): + - `nodes`: Interpolates values at the element’s nodes and assigns the element mean (ignoring NaN nodes). Fast and robust; respects node sampling. + - `nodes_scaled`: Like `nodes`, but computes a scaled mean using node-specific scaling factors (`dis.nodes.scaling_factors`), which are set via `processing.interpolation.node_scaling_factors.surface` and `processing.interpolation.node_scaling_factors.interior`. This adjusts the influence of specific nodes. + - `elementcenter`: Interpolates at each element centroid and assigns that value. + - `allvoxels`: Collects all voxels whose grid coordinates lie inside the convex hull of the element nodes; assigns the mean value; optionally filters outliers. + - `allvoxels_scaled`: Computes a voxel-weighted mean where voxel weights derive from node scaling factors and inverse node-to-voxel distances. The influence of the distance (decay) is controlled by `processing.interpolation.inverse_distance_power` p (p=1 linear, p=2 quadratic (default), p>=3 increasingly like step function); Optionally filters outliers. - Element and node value overrides: - - set_surface_node_value: If provided, all nodes that belong to any surface receive the fixed value (vector size must match the number of pixel channels); only relevant for `nodes` and `nodes_scaled` interpolation methods. - - set_surface_element_value: If provided, all elements touching any surface node receive the fixed value (scalar or vector); applicable to all interpolation methods. - - When both are set, surface node override is ignored (element override wins). + - `processing.interpolation.set_surface_node_value`: If provided, all nodes that belong to any surface receive the fixed value (vector size must match the number of pixel channels); only relevant for `nodes` and `nodes_scaled` interpolation methods. + - `processing.interpolation.set_surface_element_value`: If provided, all elements touching any surface node receive the fixed value (scalar or vector); applicable to all interpolation methods. -- Outlier filtering (processing.interpolation.filter_outliers): +- Outlier filtering (`processing.interpolation.filter_outliers`): - In `allvoxels` and `allvoxels_scaled`, if enabled and enough voxels are present (>5), outliers are removed using a modified Z-score (median/MAD-based, threshold=3.5) before averaging. - Fallbacks and warnings: diff --git a/src/i2pp/core/configuration_validator/validator.py b/src/i2pp/core/configuration_validator/validator.py index 6f247ff..475b5f1 100644 --- a/src/i2pp/core/configuration_validator/validator.py +++ b/src/i2pp/core/configuration_validator/validator.py @@ -125,9 +125,17 @@ def from_dict( return NodeScaling( surface_node_scaling=1.0, interior_node_scaling=1.0 ) + surface = d.get("surface", 1.0) + interior = d.get("interior", 1.0) + + if surface < 0 or interior < 0: + raise ValueError( + "Node scaling factors (surface and interior) must be " + "non-negative." + ) return NodeScaling( - surface_node_scaling=d.get("surface", 1.0), - interior_node_scaling=d.get("interior", 1.0), + surface_node_scaling=surface, + interior_node_scaling=interior, ) @@ -138,6 +146,7 @@ class Interpolation: method: str node_scaling_factors: NodeScaling filter_outliers: bool = False + idw_power: int = 2 set_node_value: Optional[float | list[float]] = field(default=None) set_ele_value: Optional[float | list[float]] = field(default=None) @@ -163,6 +172,10 @@ def from_dict(d: Dict[str, Any]) -> "Interpolation": "and 'set_surface_element_value' " "cannot be set at the same time." ) + if d.get("inverse_distance_power", 2) <= 0: + raise ValueError( + "Inverse distance power must be a positive integer." + ) return Interpolation( method=d["method"], filter_outliers=d.get("filter_outliers", False), @@ -171,6 +184,7 @@ def from_dict(d: Dict[str, Any]) -> "Interpolation": node_scaling_factors=NodeScaling.from_dict( d.get("node_scaling_factors") ), + idw_power=d.get("inverse_distance_power", 2), ) diff --git a/src/i2pp/core/interpolate_element_data.py b/src/i2pp/core/interpolate_element_data.py index 9b145c0..a792c0e 100644 --- a/src/i2pp/core/interpolate_element_data.py +++ b/src/i2pp/core/interpolate_element_data.py @@ -77,6 +77,7 @@ def interpolate_image_to_discretization( interpolator = enum_interpolation_method.create_interpolator( filter_outliers=interpolation.filter_outliers, set_node_value=interpolation.set_node_value, + idw_power=interpolation.idw_power, ) elements = interpolator.compute_element_data(dis, image_data) diff --git a/src/i2pp/core/interpolators/interpolator_all_voxel.py b/src/i2pp/core/interpolators/interpolator_all_voxel.py index 2c0c7b0..5e39289 100644 --- a/src/i2pp/core/interpolators/interpolator_all_voxel.py +++ b/src/i2pp/core/interpolators/interpolator_all_voxel.py @@ -30,11 +30,15 @@ def __init__( *args, filter_outliers: bool = False, mode: str = "allvoxels", + idw_power: int = 2, **kwargs, ): super().__init__(*args, **kwargs) self._filter_outliers_enabled = filter_outliers self._mode = mode # "allvoxels" or "allvoxels_scaled" + self._idw_power = ( + idw_power # Power parameter for inverse distance weighting + ) def _search_bounding_box( self, grid_coords: GridCoords, element_grid_coords: np.ndarray @@ -115,37 +119,83 @@ def _filter_outliers_modified_zscore( ) -> np.ndarray: """Identify outliers using the Modified Z-Score method. + Supports multi-channel data (e.g., RGB) by evaluating outliers across + all components. + Args: values (np.ndarray): - Array of voxel values (1D or 2D, axis=0 = voxels) + Array of voxel values (N_voxels x ...). threshold (float): - Cutoff for Modified Z-Score, default 3.5 + Cutoff for Modified Z-Score, default 3.5. Returns: - mask (np.ndarray): Boolean array, True = value is not an outlier + mask (np.ndarray): 1D boolean array of shape (N_voxels,), where + True indicates the pixel is not an outlier and should be used. """ values = np.asarray(values) + if values.ndim == 1: + values = values[:, np.newaxis] + + # Calculate median and MAD per channel (axis=0) median = np.median(values, axis=0) # Median Absolute Deviation (MAD) mad = np.median(np.abs(values - median), axis=0) mad = np.maximum(mad, 1e-10) # avoid division by zero - # Modified Z-Score + # Modified Z-Score for each component modified_z = 0.6745 * (values - median) / mad - # Mask: True = not an outlier - mask = np.abs(modified_z) <= threshold + # Mask: True if NOT an outlier across all components/channels + # A voxel is filtered out if any channel exceeds the threshold + mask = np.all(np.abs(modified_z) <= threshold, axis=1) return mask + def _compute_idw_voxel_weights( + self, + element_node_phys: np.ndarray, + voxels_phys: np.ndarray, + node_scaling_factors: np.ndarray, + ) -> np.ndarray: + """Computes inverse distance weighting (IDW) voxel weights. + + Args: + element_node_phys (np.ndarray): + Physical coordinates of the element's nodes. + voxels_phys (np.ndarray): + Physical coordinates of voxels within the element. + node_scaling_factors (np.ndarray): + Scaling factors assigned to the nodes. + + Returns: + np.ndarray: The computed weights for each voxel. + """ + + distances = np.linalg.norm( + element_node_phys[:, np.newaxis, :] + - voxels_phys[np.newaxis, :, :], + axis=2, + ) + + # Avoid division by zero by setting a minimum distance threshold + eps = 1e-9 + distances = np.maximum(distances, eps) + + inv_distances = 1.0 / distances**self._idw_power + voxel_weights = np.sum( + node_scaling_factors[:, np.newaxis] * inv_distances, axis=0 + ) / np.sum(inv_distances, axis=0) + + return voxel_weights + def _weighted_voxel_mean( self, voxels_phys: np.ndarray, values: np.ndarray, element_node_phys: np.ndarray, node_scaling_factors: np.ndarray, - ) -> float: + ) -> float | np.ndarray: """Calculate a node-weighted mean of voxel values, with optional outlier filtering. @@ -168,18 +218,19 @@ def _weighted_voxel_mean( float: The weighted mean of the voxel values. """ - distances = np.linalg.norm( - element_node_phys[:, np.newaxis, :] - - voxels_phys[np.newaxis, :, :], - axis=2, - ) - distances = np.maximum(distances, 1e-10) + # if node_scaling factors are the same for all nodes, + # this reduces to a standard IDW weighted mean + if np.all(node_scaling_factors == node_scaling_factors[0]): + voxel_weights = ( + np.ones(voxels_phys.shape[0]) * node_scaling_factors[0] + ) - voxel_weights = np.sum( - node_scaling_factors[:, np.newaxis] / distances, axis=0 - ) + else: + voxel_weights = self._compute_idw_voxel_weights( + element_node_phys, voxels_phys, node_scaling_factors + ) - if len(values) > 5: + if self._filter_outliers_enabled and len(values) > 5: mask = self._filter_outliers_modified_zscore(values) filtered_values = values[mask] filtered_weights = voxel_weights[mask] @@ -208,9 +259,10 @@ def _format_output_value( ) -> np.ndarray: """Formats the output value based on the pixel type. - If the pixel type has one value, it wraps the value in a NumPy array. - If it has multiple values, it creates a NumPy array filled with the - value. + If the value is already a vector matching the pixel type components, it + is returned. Otherwise, if the pixel type has one value, it wraps the + value in a NumPy array. If it has multiple values and input is scalar, + it creates a NumPy array filled with the value. Args: value (float | np.ndarray): The value to format. @@ -219,6 +271,10 @@ def _format_output_value( Returns: np.ndarray: The formatted output value. """ + val_arr = np.atleast_1d(value) + if val_arr.size == image_data.pixel_type.num_values: + return val_arr + if image_data.pixel_type.num_values == 1: return np.array([value]) return np.full(image_data.pixel_type.num_values, value) @@ -311,33 +367,12 @@ def _get_data_of_element( if node_scaling_factors_current is None: weighted = np.mean(data, axis=0) else: - if self._filter_outliers_enabled: - weighted = self._weighted_voxel_mean( - voxels_phys_np, - data, - element_node_phys, - node_scaling_factors_current, - ) - else: - # Non-filtered weighted path with zero-sum guard - distances = np.linalg.norm( - element_node_phys[:, np.newaxis, :] - - voxels_phys_np[np.newaxis, :, :], - axis=2, - ) - distances = np.maximum(distances, 1e-10) - voxel_weights = np.sum( - node_scaling_factors_current[:, np.newaxis] - / distances, - axis=0, - ) - - if np.sum(voxel_weights) <= 0: - weighted = np.mean(data, axis=0) - else: - weighted = np.average( - data, weights=voxel_weights, axis=0 - ) + weighted = self._weighted_voxel_mean( + voxels_phys_np, + data, + element_node_phys, + node_scaling_factors_current, + ) return self._format_output_value(weighted, image_data) diff --git a/src/i2pp/core/interpolators/interpolator_types.py b/src/i2pp/core/interpolators/interpolator_types.py index 67fd504..5dbe12d 100644 --- a/src/i2pp/core/interpolators/interpolator_types.py +++ b/src/i2pp/core/interpolators/interpolator_types.py @@ -43,6 +43,7 @@ def create_interpolator( *, filter_outliers: bool = False, set_node_value: float | list[float] | None = None, + idw_power: int = 2, ) -> Interpolator: """Creates and returns a configured interpolator instance based on the selected interpolation method. @@ -57,6 +58,8 @@ def create_interpolator( interpolation. Defaults to False. set_node_value (float | None): Value to set for surface nodes. Defaults to None. + idw_power (int): Power parameter for inverse distance weighting. + Only applicable for ALLVOXELS_SCALED method. Defaults to 2. Returns: Interpolator: An instance of the interpolator that matches the @@ -71,7 +74,9 @@ def create_interpolator( ) if self == InterpolationType.ALLVOXELS_SCALED: return InterpolatorAllVoxel( - mode="allvoxels_scaled", filter_outliers=filter_outliers + mode="allvoxels_scaled", + filter_outliers=filter_outliers, + idw_power=idw_power, ) if self == InterpolationType.NODES: return InterpolatorNodes( diff --git a/tests/unittests/i2pp/core/configuration_validator/test_validator.py b/tests/unittests/i2pp/core/configuration_validator/test_validator.py index bf91512..6ed9c8b 100644 --- a/tests/unittests/i2pp/core/configuration_validator/test_validator.py +++ b/tests/unittests/i2pp/core/configuration_validator/test_validator.py @@ -156,6 +156,22 @@ def test_config_validation_non_existing_path(tmp_path): with pytest.raises( FileNotFoundError, - match=f"Path does not exist: {str(tmp_path / "non_existing.mesh")}", + match=f"Path does not exist: {str(tmp_path / 'non_existing.mesh')}", + ): + _ = I2PPConfig.from_dict(invalid_config) + + +def test_config_validation_both_surface_values_set(minimal_valid_config): + """Test that an error is raised if both surface node and element values are + set.""" + invalid_config = minimal_valid_config.copy() + invalid_int_config = invalid_config["processing"]["interpolation"] + invalid_int_config["set_surface_node_value"] = 1.0 + invalid_int_config["set_surface_element_value"] = 2.0 + + with pytest.raises( + ValueError, + match="Both 'set_surface_node_value' and 'set_surface_element_value' " + "cannot be set at the same time.", ): _ = I2PPConfig.from_dict(invalid_config) diff --git a/tests/unittests/i2pp/core/interpolators/test_interpolator_all_voxel.py b/tests/unittests/i2pp/core/interpolators/test_interpolator_all_voxel.py index fe8fdc8..cfa61c6 100644 --- a/tests/unittests/i2pp/core/interpolators/test_interpolator_all_voxel.py +++ b/tests/unittests/i2pp/core/interpolators/test_interpolator_all_voxel.py @@ -684,3 +684,287 @@ def test_scaled_mode_with_filter_falls_back_when_scaling_factors_sum_zero(): val_unscaled = float(np.asarray(elems_unscaled[0].data).mean()) assert np.isclose(val_scaled, val_unscaled, rtol=1e-6) + + +def test_different_idw_power_setups(): + """Test that changing idw_power affects the interpolated values.""" + # Image: 12x12x12 grid, values depend on x-coordinate + pixel_data = np.zeros((12, 12, 12), dtype=float) + # make pixel data random but reproducible + # and deterministic based on x-coordinate + for z in range(12): + for y in range(12): + for x in range(12): + pixel_data[z, y, x] = x + 0.1 * (z + y) + + grid_coords = GridCoords( + slice=np.arange(12), row=np.arange(12), col=np.arange(12) + ) + image_data = ImageData( + pixel_data, grid_coords, np.eye(3), np.zeros(3), PixelValueType.CT + ) + + # Element: cube corners spanning [0, 10] in each axis + ele_nodes = np.array( + [ + [0, 0, 0], + [10, 0, 0], + [0, 10, 0], + [10, 10, 0], + [0, 0, 10], + [10, 0, 10], + [0, 10, 10], + [10, 10, 10], + ], + dtype=float, + ) + + # Node scaling factors: heavily weight the nodes at x=0 + scaling_factors = np.array( + [10.0, 1.0, 10.0, 1.0, 10.0, 1.0, 10.0, 1.0], dtype=float + ) + + dis = _make_simple_discretization( + node_coords=ele_nodes, + elements_node_ids=[np.arange(8)], + node_scaling_factors=scaling_factors, + ) + + # Test with idw_power = 1 + interp_p1 = InterpolatorAllVoxel( + mode="allvoxels_scaled", idw_power=1, filter_outliers=False + ) + val_p1 = float( + np.asarray( + interp_p1.compute_element_data(dis, image_data)[0].data + ).mean() + ) + + # Test with idw_power = 2 + interp_p2 = InterpolatorAllVoxel( + mode="allvoxels_scaled", idw_power=2, filter_outliers=False + ) + val_p2 = float( + np.asarray( + interp_p2.compute_element_data(dis, image_data)[0].data + ).mean() + ) + + # Test with idw_power = 3 + interp_p3 = InterpolatorAllVoxel( + mode="allvoxels_scaled", idw_power=3, filter_outliers=False + ) + val_p3 = float( + np.asarray( + interp_p3.compute_element_data(dis, image_data)[0].data + ).mean() + ) + + # Ensure they are all finite + assert np.isfinite(val_p1) + assert np.isfinite(val_p2) + assert np.isfinite(val_p3) + + # Ensure they are different due to different distance weighting + assert not np.isclose(val_p1, val_p2) + assert not np.isclose(val_p2, val_p3) + assert not np.isclose(val_p1, val_p3) + + +def test_large_element_many_voxels(): + """Test interpolation for a large element containing many voxels.""" + # Grid 20x20x20 + grid_coords = GridCoords( + slice=np.arange(20), row=np.arange(20), col=np.arange(20) + ) + + # Pixel data varies linearly along the x-axis (col) + # pixel_data[z, y, x] = x + pixel_data = np.broadcast_to(np.arange(20), (20, 20, 20)).astype(float) + image_data = ImageData( + pixel_data, grid_coords, np.eye(3), np.zeros(3), PixelValueType.CT + ) + + # Element: large cube from coordinates 2 to 17 + # This will enclose 16x16x16 = 4096 voxels + ele_nodes = np.array( + [ + [2, 2, 2], + [17, 2, 2], + [2, 17, 2], + [17, 17, 2], + [2, 2, 17], + [17, 2, 17], + [2, 17, 17], + [17, 17, 17], + ], + dtype=float, + ) + + dis = _make_simple_discretization( + node_coords=ele_nodes, + elements_node_ids=[np.arange(8)], + node_scaling_factors=np.ones(8, dtype=float), + ) + + # Unscaled mode + interp_unscaled = InterpolatorAllVoxel( + mode="allvoxels", filter_outliers=False + ) + val_unscaled = float( + np.asarray( + interp_unscaled.compute_element_data(dis, image_data)[0].data + ).mean() + ) + + # Scaled mode, uniformly scaled + interp_uni_scaled = InterpolatorAllVoxel( + mode="allvoxels_scaled", filter_outliers=False + ) + val_uni_scaled = float( + np.asarray( + interp_uni_scaled.compute_element_data(dis, image_data)[0].data + ).mean() + ) + + # Scaled mode, different node scaling factors + dis.nodes.scaling_factors = np.array([0, 0, 1, 0, 1, 0, 1, 0], dtype=float) + interp_diff_scaled = InterpolatorAllVoxel( + mode="allvoxels_scaled", filter_outliers=False + ) + val_diff_scaled = float( + np.asarray( + interp_diff_scaled.compute_element_data(dis, image_data)[0].data + ).mean() + ) + + # The x-coordinates inside the inclusive range [2, 17] are 2, 3, ..., 17. + # The mean of this sequence is (2 + 17) / 2 = 9.5. + expected_mean = 9.5 + + assert np.isclose(val_unscaled, expected_mean) + + # With uniform scaling factors and a symmetric element/data distribution, + # the scaled mean should also match the expected mean. + assert np.isclose(val_uni_scaled, expected_mean) + + # With different node scaling factors, the mean should differ. + assert not np.isclose(val_diff_scaled, expected_mean) + + +def test__compute_idw_voxel_weights_different_scaling_factors(): + """Test that different node scaling factors produce different voxel + weights.""" + interpolator = InterpolatorAllVoxel(mode="allvoxels_scaled", idw_power=2) + + # 4 nodes (e.g., a tetrahedron) + element_node_phys = np.array( + [ + [0.0, 0.0, 0.0], + [1.0, 0.0, 0.0], + [0.0, 1.0, 0.0], + [0.0, 0.0, 1.0], + ] + ) + + # 2 voxels + voxels_phys = np.array( + [ + [0.1, 0.1, 0.1], # Closer to node 0 + [0.8, 0.1, 0.1], # Closer to node 1 + ] + ) + + # Uniform scaling factors + scaling_uniform = np.array([1.0, 1.0, 1.0, 1.0]) + weights_uniform = interpolator._compute_idw_voxel_weights( + element_node_phys, voxels_phys, scaling_uniform + ) + + # Biased towards node 0 + scaling_biased_0 = np.array([10.0, 1.0, 1.0, 1.0]) + weights_biased_0 = interpolator._compute_idw_voxel_weights( + element_node_phys, voxels_phys, scaling_biased_0 + ) + + # Biased towards node 1 + scaling_biased_1 = np.array([1.0, 10.0, 1.0, 1.0]) + weights_biased_1 = interpolator._compute_idw_voxel_weights( + element_node_phys, voxels_phys, scaling_biased_1 + ) + + assert weights_uniform.shape == (2,) + assert weights_biased_0.shape == (2,) + assert weights_biased_1.shape == (2,) + + # Uniform scaling factors should result in uniform weights (all 1.0) + assert np.allclose(weights_uniform, 1.0) + + # Biased weights should differ from uniform + assert not np.allclose(weights_uniform, weights_biased_0) + assert not np.allclose(weights_uniform, weights_biased_1) + + # Biased weights should differ from each other + assert not np.allclose(weights_biased_0, weights_biased_1) + + # Voxel 0 is closer to Node 0, so its weight should be + # higher when Node 0 is biased + assert weights_biased_0[0] > weights_biased_0[1] + + # Voxel 1 is closer to Node 1, so its weight should be + # higher when Node 1 is biased + assert weights_biased_1[1] > weights_biased_1[0] + + +def test__compute_idw_voxel_weights_different_idw_powers(): + """Test that different idw_power values produce different voxel weights.""" + # 4 nodes (e.g., a tetrahedron) + element_node_phys = np.array( + [ + [0.0, 0.0, 0.0], + [1.0, 0.0, 0.0], + [0.0, 1.0, 0.0], + [0.0, 0.0, 1.0], + ] + ) + + # 2 voxels + voxels_phys = np.array( + [ + [0.1, 0.1, 0.1], # Closer to node 0 + [0.8, 0.1, 0.1], # Closer to node 1 + ] + ) + + # Must use non-uniform scaling factors, otherwise weights are always 1.0 + scaling_factors = np.array([10.0, 1.0, 1.0, 1.0]) + + interpolator_p1 = InterpolatorAllVoxel( + mode="allvoxels_scaled", idw_power=1 + ) + weights_p1 = interpolator_p1._compute_idw_voxel_weights( + element_node_phys, voxels_phys, scaling_factors + ) + + interpolator_p2 = InterpolatorAllVoxel( + mode="allvoxels_scaled", idw_power=2 + ) + weights_p2 = interpolator_p2._compute_idw_voxel_weights( + element_node_phys, voxels_phys, scaling_factors + ) + + interpolator_p3 = InterpolatorAllVoxel( + mode="allvoxels_scaled", idw_power=3 + ) + weights_p3 = interpolator_p3._compute_idw_voxel_weights( + element_node_phys, voxels_phys, scaling_factors + ) + + assert weights_p1.shape == (2,) + assert weights_p2.shape == (2,) + assert weights_p3.shape == (2,) + + # Ensure the weights differ due to different distance weighting powers + assert not np.allclose(weights_p1, weights_p2) + assert not np.allclose(weights_p2, weights_p3) + assert not np.allclose(weights_p1, weights_p3) diff --git a/tests/unittests/i2pp/core/test_run.py b/tests/unittests/i2pp/core/test_run.py index 6594497..cd7e27a 100644 --- a/tests/unittests/i2pp/core/test_run.py +++ b/tests/unittests/i2pp/core/test_run.py @@ -150,8 +150,9 @@ def maximal_valid_config(tmp_path): }, "processing": { "interpolation": { - "method": "nodes_scaled", + "method": "allvoxels_scaled", "node_scaling_factors": {"interior": 0.8, "surface": 0.2}, + "inverse_distance_power": 3, "filter_outliers": True, "set_surface_node_value": 100.0, },