Solving Issue: Handling interpolation at boundary

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,27 @@ 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_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:
+  - `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`):
+  - 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.
+  - 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

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
@@ -67,16 +68,23 @@ 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
     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(
-            smoothing_area=d.get("smoothing_area", 3),
+            area=d.get("area", 3),
             visualize=d.get("visualize", False),
         )
 
@@ -101,25 +109,105 @@ def from_dict(d: Dict[str, Any]) -> "Transformation":
         )
 
 
+@dataclass
+class NodeScaling:
+    """Class representing node scaling configuration for interpolation."""
+
+    surface_node_scaling: float = 1.0
+    interior_node_scaling: float = 1.0
+
+    @staticmethod
+    def from_dict(
+        d: Optional[Dict[str, Any]],
+    ) -> "NodeScaling":
+        """Creates a NodeScaling instance from a dictionary."""
+        if d is None:
+            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=surface,
+            interior_node_scaling=interior,
+        )
+
+
+@dataclass
+class Interpolation:
+    """Class representing the interpolation configuration."""
+
+    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)
+
+    @staticmethod
+    def from_dict(d: Dict[str, Any]) -> "Interpolation":
+        """Creates an Interpolation instance from a dictionary."""
+
+        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."
+            )
+        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),
+            set_node_value=d.get("set_surface_node_value"),
+            set_ele_value=d.get("set_surface_element_value"),
+            node_scaling_factors=NodeScaling.from_dict(
+                d.get("node_scaling_factors")
+            ),
+            idw_power=d.get("inverse_distance_power", 2),
+        )
+
+
 @dataclass
 class Processing:
     """Class representing the processing configuration."""
 
-    interpolation_method: str
     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_method=d["interpolation_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"]),
+            interpolation=Interpolation.from_dict(d["interpolation"]),
         )
 
 

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,
 )
@@ -50,7 +51,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: Processing,
 ) -> Discretization:
     """Loads and processes 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)
 
@@ -94,8 +99,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 +155,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:

src/i2pp/core/discretization_readers/discretization_reader.py

@@ -1,12 +1,17 @@
 """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
 
+if TYPE_CHECKING:
+    from i2pp.core.configuration_validator.validator import Processing
+
 
 @dataclass
 class BoundingBox:
@@ -31,10 +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.
+        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
+    scaling_factors: Optional[np.ndarray] = None
 
 
 @dataclass
@@ -62,6 +72,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, a 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 +104,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
 
 
@@ -102,7 +133,10 @@ def __init__(self):
 
     @abstractmethod
     def load_discretization(
-        self, file_path: Path, options: dict
+        self,
+        file_path: Path,
+        options: dict,
+        processing: Processing,
     ) -> Discretization:
         """Abstract method to load discretization data from a file path.
 
@@ -116,6 +150,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