Skip to content

Filters API

Pictologics provides IBSI 2-compliant convolutional filters for image response map generation.

Overview

All filters can be used via the RadiomicsPipeline filter step or called directly:

# Pipeline usage
{"step": "filter", "params": {"type": "log", "sigma_mm": 1.5}}

# Direct usage
from pictologics.filters import laplacian_of_gaussian, BoundaryCondition
response = laplacian_of_gaussian(image.array, sigma_mm=1.5, spacing_mm=image.spacing)

Available Filters

Filter Function Use Case
Mean mean_filter Local averaging
LoG laplacian_of_gaussian Edge/blob detection
Laws laws_filter Texture energy
Gabor gabor_filter Directional patterns
Wavelet wavelet_transform Multi-resolution analysis
Simoncelli simoncelli_wavelet Non-separable wavelet
Riesz riesz_transform, riesz_log, riesz_simoncelli Rotation-equivariant transforms

Capability Metadata

Versioned, machine-readable descriptions of what each filter supports — input and kernel dimensionality, plane-wise execution, orthogonal-plane averaging, rotation pooling, supported and effective boundary handling, Riesz orders, structure-tensor steering, and anisotropic-spacing behaviour. Intended for compliance tooling, so support can be determined without inspecting signatures or private source.

Note the distinction between supported_boundaries (what the function accepts) and effective_boundary (what physically happens at the border): the FFT-based filters report as_specified_via_padding, because a requested non-periodic boundary is realised by a defined pad-filter-crop procedure and the transform itself remains periodic on the padded domain.

pictologics.filters.CAPABILITIES_SCHEMA_VERSION = '1.0.0' module-attribute

Semantic version of the :data:FILTER_CAPABILITIES schema.

Bump the major component on breaking changes (field removal/retyping), the minor component when adding fields or filter entries, and the patch component for corrections that don't change the schema shape.

pictologics.filters.FilterCapability dataclass

Declarative capability record for one IBSI 2 filter (or Riesz variant).

Every field is grounded in the corresponding filter implementation; see :data:FILTER_CAPABILITIES for the per-filter values and the source evidence behind each one.

Attributes:

Name Type Description
input_dimensionality tuple[int, ...]

Dimensionality/dimensionalities of the top-level image array the filter function accepts, e.g. (3,) for 3D-only.
kernel_dimensionality int

Dimensionality of the convolution/transform kernel itself (2 or 3). May differ from input_dimensionality, e.g. Gabor uses a 2D kernel on a 3D input.
slice_plane_execution bool

True if the filter applies its kernel plane-wise (slice-by-slice) rather than as a single volumetric operation.
orthogonal_plane_averaging bool

True if the filter supports averaging its plane-wise response over the 3 orthogonal anatomical planes.
rotation_pooling tuple[str, ...]

Pooling method names accepted for pseudo-rotational invariance (e.g. ("max", "average", "min")), or () if the filter has no rotation-invariance/pooling mechanism.
supported_boundaries tuple[str, ...]

BoundaryCondition member names the filter actually accepts via its boundary parameter, or () if it has no such parameter.
effective_boundary str

What actually happens at the image border. One of "as_specified" (a spatial-convolution filter that applies the requested boundary directly), "as_specified_via_padding" (FFT-domain filters: the requested boundary is realised by a defined pad-filter-crop procedure — the transform itself remains periodic on the padded domain, so the boundary is honoured approximately rather than exactly), or "periodic" (periodic regardless of what is requested).
supported_riesz_orders Optional[str]

Description of the Riesz derivative orders accepted (Riesz-family filters only), or None if the filter has no Riesz order concept.
structure_tensor_steering bool

True if the filter can steer its kernel using a structure tensor. Currently False for every filter (no tensor_sigma/steering exists anywhere in the package).
anisotropic_spacing str

How the filter handles anisotropic spacing_mm: "supported" (correctly converts per-axis), "warns_uses_first_axis" (warns and derives scale from the first axis only), or "not_applicable" (no spacing_mm parameter / concept).
Source code in pictologics/filters/capabilities.py 
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
@dataclass(frozen=True)
class FilterCapability:
    """
    Declarative capability record for one IBSI 2 filter (or Riesz variant).

    Every field is grounded in the corresponding filter implementation; see
    :data:`FILTER_CAPABILITIES` for the per-filter values and the source
    evidence behind each one.

    Attributes:
        input_dimensionality: Dimensionality/dimensionalities of the top-level
            image array the filter function accepts, e.g. ``(3,)`` for
            3D-only.
        kernel_dimensionality: Dimensionality of the convolution/transform
            kernel itself (2 or 3). May differ from `input_dimensionality`,
            e.g. Gabor uses a 2D kernel on a 3D input.
        slice_plane_execution: True if the filter applies its kernel
            plane-wise (slice-by-slice) rather than as a single volumetric
            operation.
        orthogonal_plane_averaging: True if the filter supports averaging its
            plane-wise response over the 3 orthogonal anatomical planes.
        rotation_pooling: Pooling method names accepted for pseudo-rotational
            invariance (e.g. ``("max", "average", "min")``), or ``()`` if the
            filter has no rotation-invariance/pooling mechanism.
        supported_boundaries: `BoundaryCondition` member names the filter
            actually accepts via its `boundary` parameter, or ``()`` if it has
            no such parameter.
        effective_boundary: What actually happens at the image border. One of
            ``"as_specified"`` (a spatial-convolution filter that applies the
            requested boundary directly), ``"as_specified_via_padding"``
            (FFT-domain filters: the requested boundary is realised by a
            defined pad-filter-crop procedure — the transform itself remains
            periodic on the padded domain, so the boundary is honoured
            approximately rather than exactly), or ``"periodic"`` (periodic
            regardless of what is requested).
        supported_riesz_orders: Description of the Riesz derivative orders
            accepted (Riesz-family filters only), or `None` if the filter has
            no Riesz order concept.
        structure_tensor_steering: True if the filter can steer its kernel
            using a structure tensor. Currently `False` for every filter (no
            `tensor_sigma`/steering exists anywhere in the package).
        anisotropic_spacing: How the filter handles anisotropic
            `spacing_mm`: ``"supported"`` (correctly converts per-axis),
            ``"warns_uses_first_axis"`` (warns and derives scale from the
            first axis only), or ``"not_applicable"`` (no `spacing_mm`
            parameter / concept).
    """

    input_dimensionality: tuple[int, ...]
    kernel_dimensionality: int
    slice_plane_execution: bool
    orthogonal_plane_averaging: bool
    rotation_pooling: tuple[str, ...]
    supported_boundaries: tuple[str, ...]
    effective_boundary: str
    supported_riesz_orders: Optional[str]
    structure_tensor_steering: bool
    anisotropic_spacing: str

pictologics.filters.FILTER_CAPABILITIES = {'mean': FilterCapability(input_dimensionality=(3,), kernel_dimensionality=3, slice_plane_execution=False, orthogonal_plane_averaging=False, rotation_pooling=(), supported_boundaries=_ALL_BOUNDARIES, effective_boundary='as_specified', supported_riesz_orders=None, structure_tensor_steering=False, anisotropic_spacing='not_applicable'), 'log': FilterCapability(input_dimensionality=(3,), kernel_dimensionality=3, slice_plane_execution=False, orthogonal_plane_averaging=False, rotation_pooling=(), supported_boundaries=_ALL_BOUNDARIES, effective_boundary='as_specified', supported_riesz_orders=None, structure_tensor_steering=False, anisotropic_spacing='supported'), 'laws': FilterCapability(input_dimensionality=(3,), kernel_dimensionality=3, slice_plane_execution=False, orthogonal_plane_averaging=False, rotation_pooling=('max', 'average', 'min'), supported_boundaries=_ALL_BOUNDARIES, effective_boundary='as_specified', supported_riesz_orders=None, structure_tensor_steering=False, anisotropic_spacing='not_applicable'), 'gabor': FilterCapability(input_dimensionality=(3,), kernel_dimensionality=2, slice_plane_execution=True, orthogonal_plane_averaging=True, rotation_pooling=('max', 'average', 'min'), supported_boundaries=_ALL_BOUNDARIES, effective_boundary='as_specified', supported_riesz_orders=None, structure_tensor_steering=False, anisotropic_spacing='supported'), 'wavelet': FilterCapability(input_dimensionality=(3,), kernel_dimensionality=3, slice_plane_execution=False, orthogonal_plane_averaging=False, rotation_pooling=('max', 'average', 'min'), supported_boundaries=_ALL_BOUNDARIES, effective_boundary='as_specified', supported_riesz_orders=None, structure_tensor_steering=False, anisotropic_spacing='not_applicable'), 'simoncelli': FilterCapability(input_dimensionality=(3,), kernel_dimensionality=3, slice_plane_execution=False, orthogonal_plane_averaging=False, rotation_pooling=(), supported_boundaries=_ALL_BOUNDARIES, effective_boundary='as_specified_via_padding', supported_riesz_orders=None, structure_tensor_steering=False, anisotropic_spacing='not_applicable'), 'riesz': FilterCapability(input_dimensionality=(3,), kernel_dimensionality=3, slice_plane_execution=False, orthogonal_plane_averaging=False, rotation_pooling=(), supported_boundaries=_ALL_BOUNDARIES, effective_boundary='as_specified_via_padding', supported_riesz_orders='Tuple[int, ...] (l1, ..., ld): any non-negative ints summing to L >= 1; see get_riesz_orders(max_order, ndim) to enumerate all combinations.', structure_tensor_steering=False, anisotropic_spacing='not_applicable'), 'riesz_log': FilterCapability(input_dimensionality=(3,), kernel_dimensionality=3, slice_plane_execution=False, orthogonal_plane_averaging=False, rotation_pooling=(), supported_boundaries=_ALL_BOUNDARIES, effective_boundary='as_specified_via_padding', supported_riesz_orders='Tuple[int, ...] (l1, ..., ld): any non-negative ints summing to L >= 1; see get_riesz_orders(max_order, ndim) to enumerate all combinations.', structure_tensor_steering=False, anisotropic_spacing='supported'), 'riesz_simoncelli': FilterCapability(input_dimensionality=(3,), kernel_dimensionality=3, slice_plane_execution=False, orthogonal_plane_averaging=False, rotation_pooling=(), supported_boundaries=_ALL_BOUNDARIES, effective_boundary='as_specified_via_padding', supported_riesz_orders='Tuple[int, ...] (l1, ..., ld): any non-negative ints summing to L >= 1; see get_riesz_orders(max_order, ndim) to enumerate all combinations.', structure_tensor_steering=False, anisotropic_spacing='not_applicable')} module-attribute

Capability record for every public filter, keyed by pipeline filter-type name.

See the module docstring for the "riesz" / "riesz_log" / "riesz_simoncelli" keying scheme.

pictologics.filters.get_filter_capabilities(name)

Look up the capability record for a filter by its pipeline filter-type name.

Parameters:

Name Type Description Default
name str

Pipeline filter-type name, e.g. "mean", "log", "laws", "gabor", "wavelet", "simoncelli", "riesz", "riesz_log", or "riesz_simoncelli".

 required

Returns:

Type Description
FilterCapability

The FilterCapability record registered under name.

Raises:

Type Description
ValueError

If name is not a key in FILTER_CAPABILITIES.
Example 
from pictologics.filters import get_filter_capabilities

capability = get_filter_capabilities("riesz_log")
print(capability.anisotropic_spacing)
# supported
Source code in pictologics/filters/capabilities.py 
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
def get_filter_capabilities(name: str) -> FilterCapability:
    """
    Look up the capability record for a filter by its pipeline filter-type name.

    Args:
        name: Pipeline filter-type name, e.g. ``"mean"``, ``"log"``, ``"laws"``,
            ``"gabor"``, ``"wavelet"``, ``"simoncelli"``, ``"riesz"``,
            ``"riesz_log"``, or ``"riesz_simoncelli"``.

    Returns:
        The `FilterCapability` record registered under `name`.

    Raises:
        ValueError: If `name` is not a key in `FILTER_CAPABILITIES`.

    Example:
        ```python
        from pictologics.filters import get_filter_capabilities

        capability = get_filter_capabilities("riesz_log")
        print(capability.anisotropic_spacing)
        # supported
        ```
    """
    try:
        return FILTER_CAPABILITIES[name]
    except KeyError as exc:
        valid = ", ".join(sorted(FILTER_CAPABILITIES))
        raise ValueError(f"Unknown filter name {name!r}; valid names are: {valid}") from exc

Boundary Conditions

pictologics.filters.BoundaryCondition

Bases: Enum

IBSI 2 boundary conditions for image padding (GBYQ).

Maps to scipy.ndimage mode parameter values.

Example 
from pictologics.filters import BoundaryCondition

boundary = BoundaryCondition.MIRROR
print(boundary.value)
# "reflect"

# Also constructible from the scipy mode string used by get_scipy_mode
boundary = BoundaryCondition["MIRROR"]
Source code in pictologics/filters/base.py 
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
class BoundaryCondition(Enum):
    """
    IBSI 2 boundary conditions for image padding (GBYQ).

    Maps to scipy.ndimage mode parameter values.

    Example:
        ```python
        from pictologics.filters import BoundaryCondition

        boundary = BoundaryCondition.MIRROR
        print(boundary.value)
        # "reflect"

        # Also constructible from the scipy mode string used by get_scipy_mode
        boundary = BoundaryCondition["MIRROR"]
        ```
    """

    ZERO = "constant"  # Zero padding (Z3VE)
    NEAREST = "nearest"  # Nearest value padding (SIJG)
    PERIODIC = "wrap"  # Periodic/wrap padding (Z7YO)
    MIRROR = "reflect"  # Mirror/symmetric padding (ZDTV)

pictologics.filters.FilterResult dataclass

Container for filter response maps and metadata.

Example 
import numpy as np
from pictologics.filters.base import FilterResult

result = FilterResult(
    response_map=np.zeros((4, 4, 4), dtype=np.float32),
    filter_name="mean",
    filter_params={"support": 3},
)
print(result.shape, result.dtype)
# (4, 4, 4) float32
Source code in pictologics/filters/base.py 
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
@dataclass
class FilterResult:
    """Container for filter response maps and metadata.

    Example:
        ```python
        import numpy as np
        from pictologics.filters.base import FilterResult

        result = FilterResult(
            response_map=np.zeros((4, 4, 4), dtype=np.float32),
            filter_name="mean",
            filter_params={"support": 3},
        )
        print(result.shape, result.dtype)
        # (4, 4, 4) float32
        ```
    """

    response_map: npt.NDArray[np.floating[Any]]
    filter_name: str
    filter_params: Dict[str, Any]

    @property
    def shape(self) -> tuple[int, ...]:
        """Shape of the response map."""
        return self.response_map.shape  # type: ignore[no-any-return]

    @property
    def dtype(self) -> np.dtype[Any]:
        """Data type of the response map."""
        return self.response_map.dtype  # type: ignore[no-any-return]

dtype property

Data type of the response map.

shape property

Shape of the response map.

pictologics.filters.LAWS_KERNELS = _LAWS_KERNELS module-attribute

Dictionary of normalized Laws kernels (IBSI 2 Table 6).

Filter Functions

pictologics.filters.mean_filter(image, support=15, boundary=BoundaryCondition.ZERO, source_mask=None) 

mean_filter(
    image: npt.NDArray[np.floating[Any]],
    support: int = ...,
    boundary: Union[BoundaryCondition, str] = ...,
    source_mask: None = ...,
) -> npt.NDArray[np.floating[Any]]

mean_filter(
    image: npt.NDArray[np.floating[Any]],
    support: int = ...,
    boundary: Union[BoundaryCondition, str] = ...,
    source_mask: npt.NDArray[np.bool_] = ...,
) -> tuple[
    npt.NDArray[np.floating[Any]], npt.NDArray[np.bool_]
]

Apply 3D mean filter (IBSI code: S60F).

The mean filter computes the average intensity over an M×M×M spatial support. Per IBSI 2 Eq. 2.

Parameters:

Name Type Description Default
image NDArray[floating[Any]]

3D input image array

 required
support int

Filter support M in voxels (must be odd, YNOF)

 15
boundary Union[BoundaryCondition, str]

Boundary condition for padding (GBYQ)

 ZERO
source_mask Optional[NDArray[bool_]]

Optional boolean mask where True = valid voxel. When provided, uses normalized convolution to exclude invalid (sentinel) voxels from mean computation.

 None

Returns:

Type Description
Union[NDArray[floating[Any]], tuple[NDArray[floating[Any]], NDArray[bool_]]]

If source_mask is None: Response map with same dimensions as input
Union[NDArray[floating[Any]], tuple[NDArray[floating[Any]], NDArray[bool_]]]

If source_mask provided: Tuple of (response_map, output_valid_mask)

Raises:

Type Description
ValueError

If support is not an odd positive integer
Example

Apply Mean filter with 15-voxel support:

import numpy as np
from pictologics.filters import mean_filter

# Create dummy 3D image
image = np.random.rand(50, 50, 50)

# Apply filter (original API)
response = mean_filter(image, support=15, boundary="zero")

# With source_mask for sentinel exclusion
mask = image > -1000  # Valid voxels
response, valid_mask = mean_filter(image, support=15, source_mask=mask)
Note

Support M is defined in voxel units as per IBSI specification.

Source code in pictologics/filters/mean.py 
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
def mean_filter(
    image: npt.NDArray[np.floating[Any]],
    support: int = 15,
    boundary: Union[BoundaryCondition, str] = BoundaryCondition.ZERO,
    source_mask: Optional[npt.NDArray[np.bool_]] = None,
) -> Union[
    npt.NDArray[np.floating[Any]],
    tuple[npt.NDArray[np.floating[Any]], npt.NDArray[np.bool_]],
]:
    """
    Apply 3D mean filter (IBSI code: S60F).

    The mean filter computes the average intensity over an M×M×M
    spatial support. Per IBSI 2 Eq. 2.

    Args:
        image: 3D input image array
        support: Filter support M in voxels (must be odd, YNOF)
        boundary: Boundary condition for padding (GBYQ)
        source_mask: Optional boolean mask where True = valid voxel.
            When provided, uses normalized convolution to exclude invalid
            (sentinel) voxels from mean computation.

    Returns:
        If source_mask is None: Response map with same dimensions as input
        If source_mask provided: Tuple of (response_map, output_valid_mask)

    Raises:
        ValueError: If support is not an odd positive integer

    Example:
        Apply Mean filter with 15-voxel support:

        ```python
        import numpy as np
        from pictologics.filters import mean_filter

        # Create dummy 3D image
        image = np.random.rand(50, 50, 50)

        # Apply filter (original API)
        response = mean_filter(image, support=15, boundary="zero")

        # With source_mask for sentinel exclusion
        mask = image > -1000  # Valid voxels
        response, valid_mask = mean_filter(image, support=15, source_mask=mask)
        ```

    Note:
        Support M is defined in voxel units as per IBSI specification.
    """
    # Validate support
    if support < 1 or support % 2 == 0:
        raise ValueError(f"Support must be an odd positive integer, got {support}")

    # Convert to float32 as required by IBSI
    image = ensure_float32(image)

    # Handle string boundary condition
    if isinstance(boundary, str):
        boundary = BoundaryCondition[boundary.upper()]

    mode = get_scipy_mode(boundary)

    if source_mask is not None:
        # Use normalized convolution for source masking
        return _normalized_uniform_filter(image, source_mask, size=support, mode=mode)
    else:
        # Cast to float32 for consistency with the masked path and the other filters
        # (uniform_filter accumulates in the input dtype, so a float64 image keeps
        # its precision through the running sum before the final downcast).
        return uniform_filter(image, size=support, mode=mode).astype(np.float32)  # type: ignore[no-any-return]

pictologics.filters.laplacian_of_gaussian(image, sigma_mm, spacing_mm=1.0, truncate=4.0, boundary=BoundaryCondition.ZERO, source_mask=None) 

laplacian_of_gaussian(
    image: npt.NDArray[np.floating[Any]],
    sigma_mm: float,
    spacing_mm: Union[
        float, Tuple[float, float, float]
    ] = ...,
    truncate: float = ...,
    boundary: Union[BoundaryCondition, str] = ...,
    source_mask: None = ...,
) -> npt.NDArray[np.floating[Any]]

laplacian_of_gaussian(
    image: npt.NDArray[np.floating[Any]],
    sigma_mm: float,
    spacing_mm: Union[
        float, Tuple[float, float, float]
    ] = ...,
    truncate: float = ...,
    boundary: Union[BoundaryCondition, str] = ...,
    source_mask: npt.NDArray[np.bool_] = ...,
) -> tuple[
    npt.NDArray[np.floating[Any]], npt.NDArray[np.bool_]
]

Apply 3D Laplacian of Gaussian filter (IBSI code: L6PA).

The LoG is a band-pass, spherically symmetric operator. Per IBSI 2 Eq. 3.

Parameters:

Name Type Description Default
image NDArray[floating[Any]]

3D input image array

 required
sigma_mm float

Standard deviation in mm (σ*, 41LN)

 required
spacing_mm Union[float, Tuple[float, float, float]]

Voxel spacing in mm (scalar for isotropic, or tuple)

 1.0
truncate float

Filter size cutoff in σ units (default 4.0, WGPM)

 4.0
boundary Union[BoundaryCondition, str]

Boundary condition for padding (GBYQ)

 ZERO
source_mask Optional[NDArray[bool_]]

Optional boolean mask where True = valid voxel. When provided, uses normalized convolution to exclude invalid (sentinel) voxels from computation.

 None

Returns:

Type Description
Union[NDArray[floating[Any]], tuple[NDArray[floating[Any]], NDArray[bool_]]]

If source_mask is None: Response map with same dimensions as input
Union[NDArray[floating[Any]], tuple[NDArray[floating[Any]], NDArray[bool_]]]

If source_mask provided: Tuple of (response_map, output_valid_mask)
Example

Apply LoG filter with 5.0mm sigma on an image with 2.0mm spacing:

import numpy as np
from pictologics.filters import laplacian_of_gaussian

# Create dummy 3D image
image = np.random.rand(50, 50, 50)

# Apply filter (original API)
response = laplacian_of_gaussian(
    image,
    sigma_mm=5.0,
    spacing_mm=(2.0, 2.0, 2.0),
    truncate=4.0
)

# With source_mask for sentinel exclusion
mask = image > -1000
response, valid_mask = laplacian_of_gaussian(
    image, sigma_mm=5.0, spacing_mm=2.0, source_mask=mask
)
Note
  • σ is converted from mm to voxels: σ_voxels = σ_mm / spacing_mm
  • Filter size: M = 1 + 2⌊d×σ + 0.5⌋ where d=truncate
  • The kernel should sum to approximately 0 (zero-mean)
Source code in pictologics/filters/log.py 
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
def laplacian_of_gaussian(
    image: npt.NDArray[np.floating[Any]],
    sigma_mm: float,
    spacing_mm: Union[float, Tuple[float, float, float]] = 1.0,
    truncate: float = 4.0,
    boundary: Union[BoundaryCondition, str] = BoundaryCondition.ZERO,
    source_mask: Optional[npt.NDArray[np.bool_]] = None,
) -> Union[
    npt.NDArray[np.floating[Any]],
    tuple[npt.NDArray[np.floating[Any]], npt.NDArray[np.bool_]],
]:
    """
    Apply 3D Laplacian of Gaussian filter (IBSI code: L6PA).

    The LoG is a band-pass, spherically symmetric operator. Per IBSI 2 Eq. 3.

    Args:
        image: 3D input image array
        sigma_mm: Standard deviation in mm (σ*, 41LN)
        spacing_mm: Voxel spacing in mm (scalar for isotropic, or tuple)
        truncate: Filter size cutoff in σ units (default 4.0, WGPM)
        boundary: Boundary condition for padding (GBYQ)
        source_mask: Optional boolean mask where True = valid voxel.
            When provided, uses normalized convolution to exclude invalid
            (sentinel) voxels from computation.

    Returns:
        If source_mask is None: Response map with same dimensions as input
        If source_mask provided: Tuple of (response_map, output_valid_mask)

    Example:
        Apply LoG filter with 5.0mm sigma on an image with 2.0mm spacing:

        ```python
        import numpy as np
        from pictologics.filters import laplacian_of_gaussian

        # Create dummy 3D image
        image = np.random.rand(50, 50, 50)

        # Apply filter (original API)
        response = laplacian_of_gaussian(
            image,
            sigma_mm=5.0,
            spacing_mm=(2.0, 2.0, 2.0),
            truncate=4.0
        )

        # With source_mask for sentinel exclusion
        mask = image > -1000
        response, valid_mask = laplacian_of_gaussian(
            image, sigma_mm=5.0, spacing_mm=2.0, source_mask=mask
        )
        ```

    Note:
        - σ is converted from mm to voxels: σ_voxels = σ_mm / spacing_mm
        - Filter size: M = 1 + 2⌊d×σ + 0.5⌋ where d=truncate
        - The kernel should sum to approximately 0 (zero-mean)
    """
    # Convert to float32 as required by IBSI
    image = ensure_float32(image)

    # Handle scalar spacing
    if isinstance(spacing_mm, (int, float)):
        spacing_mm = (float(spacing_mm),) * 3

    # Convert sigma from mm to voxels for each axis
    sigma_voxels = tuple(sigma_mm / s for s in spacing_mm)

    # Handle string boundary condition
    if isinstance(boundary, str):
        boundary = BoundaryCondition[boundary.upper()]

    mode = get_scipy_mode(boundary)

    if source_mask is not None:
        # Use normalized convolution for source masking
        return _normalized_gaussian_laplace(
            image, source_mask, sigma=sigma_voxels, mode=mode, truncate=truncate
        )
    else:
        # Cast to float32 for consistency with the masked path and the other filters
        # (gaussian_laplace accumulates in the input dtype, so a float64 image keeps
        # its precision through the convolution before the final downcast).
        return gaussian_laplace(  # type: ignore[no-any-return]
            image, sigma=sigma_voxels, mode=mode, truncate=truncate
        ).astype(np.float32)

pictologics.filters.laws_filter(image, kernels, boundary=BoundaryCondition.ZERO, rotation_invariant=False, pooling='max', compute_energy=False, energy_distance=7, use_parallel=None, source_mask=None) 

laws_filter(
    image: npt.NDArray[np.floating[Any]],
    kernels: str,
    boundary: Union[BoundaryCondition, str] = ...,
    rotation_invariant: bool = ...,
    pooling: str = ...,
    compute_energy: bool = ...,
    energy_distance: int = ...,
    use_parallel: Union[bool, None] = ...,
    source_mask: None = ...,
) -> npt.NDArray[np.floating[Any]]

laws_filter(
    image: npt.NDArray[np.floating[Any]],
    kernels: str,
    boundary: Union[BoundaryCondition, str] = ...,
    rotation_invariant: bool = ...,
    pooling: str = ...,
    compute_energy: bool = ...,
    energy_distance: int = ...,
    use_parallel: Union[bool, None] = ...,
    source_mask: npt.NDArray[np.bool_] = ...,
) -> Tuple[
    npt.NDArray[np.floating[Any]], npt.NDArray[np.bool_]
]

Apply 3D Laws kernel filter (IBSI code: JTXT).

Laws kernels detect texture patterns via separable 1D filters combined into 2D/3D filters via outer products.

Parameters:

Name Type Description Default
image NDArray[floating[Any]]

3D input image array

 required
kernels str

Kernel specification as string, e.g., "E5L5S5" for 3D

 required
boundary Union[BoundaryCondition, str]

Boundary condition for padding (GBYQ)

 ZERO
rotation_invariant bool

If True, apply pseudo-rotational invariance (O1AQ) using max pooling over 24 right-angle rotations

 False
pooling str

Pooling method for rotation invariance ("max", "average", "min")

 'max'
compute_energy bool

If True, compute texture energy image (PQSD)

 False
energy_distance int

Chebyshev distance δ for energy computation (I176)

 7
use_parallel Union[bool, None]

If True, use parallel processing for rotation_invariant mode. If None (default), auto-enables for images > ~128³ voxels. Only affects rotation_invariant mode.

 None
source_mask Optional[NDArray[bool_]]

Optional boolean mask where True = valid voxel. In non-rotation-invariant mode, uses normalized separable convolution to exclude invalid (sentinel) voxels. In rotation-invariant mode, invalid voxels are zero-filled as a first-order approximation (the rotated kernels preclude normalized convolution).

 None

Returns:

Type Description
Union[NDArray[floating[Any]], Tuple[NDArray[floating[Any]], NDArray[bool_]]]

If source_mask is None: Response map (or energy image if compute_energy=True)
Union[NDArray[floating[Any]], Tuple[NDArray[floating[Any]], NDArray[bool_]]]

If source_mask provided: Tuple of (response_map, output_valid_mask)

Raises:

Type Description
ValueError

If kernels does not parse into exactly 3 Laws kernel codes (wrong count or malformed string), if any parsed kernel code is not a recognized name (see LAWS_KERNELS), or if rotation_invariant=True and pooling is not "max", "average", or "min".
RuntimeError

Defensive check raised if no response was computed; not expected to occur in normal use.
Example

Apply Laws E5L5S5 kernel with rotation invariance and texture energy:

import numpy as np
from pictologics.filters import laws_filter

# Create dummy 3D image
image = np.random.rand(50, 50, 50)

# Apply filter
response = laws_filter(
    image,
    "E5L5S5",
    rotation_invariant=True,
    pooling="max",
    compute_energy=True,
    energy_distance=7
)
Note
  • Kernels are normalized (deviate from Laws' original unnormalized)
  • Energy is computed as: mean(|h|) over δ neighborhood
  • For rotation invariance, energy is computed after pooling
  • Uses separable 1D convolutions for ~8x speedup over full 3D
Source code in pictologics/filters/laws.py 
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
def laws_filter(
    image: npt.NDArray[np.floating[Any]],
    kernels: str,
    boundary: Union[BoundaryCondition, str] = BoundaryCondition.ZERO,
    rotation_invariant: bool = False,
    pooling: str = "max",
    compute_energy: bool = False,
    energy_distance: int = 7,
    use_parallel: Union[bool, None] = None,
    source_mask: Optional[npt.NDArray[np.bool_]] = None,
) -> Union[
    npt.NDArray[np.floating[Any]],
    Tuple[npt.NDArray[np.floating[Any]], npt.NDArray[np.bool_]],
]:
    """
    Apply 3D Laws kernel filter (IBSI code: JTXT).

    Laws kernels detect texture patterns via separable 1D filters combined
    into 2D/3D filters via outer products.

    Args:
        image: 3D input image array
        kernels: Kernel specification as string, e.g., "E5L5S5" for 3D
        boundary: Boundary condition for padding (GBYQ)
        rotation_invariant: If True, apply pseudo-rotational invariance (O1AQ)
                            using max pooling over 24 right-angle rotations
        pooling: Pooling method for rotation invariance ("max", "average", "min")
        compute_energy: If True, compute texture energy image (PQSD)
        energy_distance: Chebyshev distance δ for energy computation (I176)
        use_parallel: If True, use parallel processing for rotation_invariant mode.
            If None (default), auto-enables for images > ~128³ voxels.
            Only affects rotation_invariant mode.
        source_mask: Optional boolean mask where True = valid voxel.
            In non-rotation-invariant mode, uses normalized separable convolution
            to exclude invalid (sentinel) voxels. In rotation-invariant mode,
            invalid voxels are zero-filled as a first-order approximation (the
            rotated kernels preclude normalized convolution).

    Returns:
        If source_mask is None: Response map (or energy image if compute_energy=True)
        If source_mask provided: Tuple of (response_map, output_valid_mask)

    Raises:
        ValueError: If `kernels` does not parse into exactly 3 Laws kernel
            codes (wrong count or malformed string), if any parsed kernel
            code is not a recognized name (see `LAWS_KERNELS`), or if
            `rotation_invariant=True` and `pooling` is not "max", "average",
            or "min".
        RuntimeError: Defensive check raised if no response was computed;
            not expected to occur in normal use.

    Example:
        Apply Laws E5L5S5 kernel with rotation invariance and texture energy:

        ```python
        import numpy as np
        from pictologics.filters import laws_filter

        # Create dummy 3D image
        image = np.random.rand(50, 50, 50)

        # Apply filter
        response = laws_filter(
            image,
            "E5L5S5",
            rotation_invariant=True,
            pooling="max",
            compute_energy=True,
            energy_distance=7
        )
        ```

    Note:
        - Kernels are normalized (deviate from Laws' original unnormalized)
        - Energy is computed as: mean(|h|) over δ neighborhood
        - For rotation invariance, energy is computed after pooling
        - Uses separable 1D convolutions for ~8x speedup over full 3D
    """

    # Convert to float32
    image = ensure_float32(image)

    # Parse kernel names (e.g., "E5L5S5" -> ["E5", "L5", "S5"])
    kernel_names = _parse_kernel_string(kernels)
    if len(kernel_names) != 3:
        raise ValueError(f"Expected 3 kernel names for 3D, got {len(kernel_names)}: {kernel_names}")

    # Handle boundary condition
    if isinstance(boundary, str):
        boundary = BoundaryCondition[boundary.upper()]
    mode = get_scipy_mode(boundary)

    # Validate pooling method if used
    if rotation_invariant and pooling not in ("max", "average", "min"):
        raise ValueError(f"Unknown pooling method: {pooling}")

    # Get 1D kernels for separable convolution
    try:
        g1 = LAWS_KERNELS[kernel_names[0]].astype(np.float32)
        g2 = LAWS_KERNELS[kernel_names[1]].astype(np.float32)
        g3 = LAWS_KERNELS[kernel_names[2]].astype(np.float32)
    except KeyError as exc:
        valid = ", ".join(sorted(LAWS_KERNELS))
        raise ValueError(
            f"Unknown Laws kernel {exc.args[0]!r}; valid kernels are: {valid}"
        ) from exc

    # Auto-detect parallel mode based on image size
    if use_parallel is None:
        use_parallel = image.size > _PARALLEL_THRESHOLD

    result: npt.NDArray[np.floating[Any]] | None = None

    if rotation_invariant:
        # Normalized convolution isn't available on the rotated kernels, so zero-fill
        # invalid voxels as a first-order approximation (same approach the FFT-based
        # filters use). The output validity mask is the input mask (zero-fill does
        # not shrink it, unlike normalized convolution).
        if source_mask is not None:
            image = _prepare_masked_image(image, source_mask)
            valid_mask = source_mask
        else:
            valid_mask = None

        rotations = _get_rotation_permutations_3d()

        # Every Laws 1D kernel is odd-length and either symmetric (L, S, R) or
        # antisymmetric (E, W). Reversing a symmetric kernel is a no-op; reversing an
        # antisymmetric one negates it. So each of the 24 rotated separable
        # convolutions equals ±(a convolution with the kernels permuted but not
        # flipped). We therefore compute only the (at most 6) distinct permutations
        # and recover every rotation's response with a sign flip.
        kernel_arrays = [g1, g2, g3]
        antisym = [bool(np.allclose(k, -k[::-1])) for k in kernel_arrays]

        base_keys: List[Tuple[Tuple[str, str, str], Tuple[int, int, int]]] = []
        seen: set[Tuple[str, str, str]] = set()
        for perm, _flips in rotations:
            key = (kernel_names[perm[0]], kernel_names[perm[1]], kernel_names[perm[2]])
            if key not in seen:
                seen.add(key)
                base_keys.append((key, perm))

        def _base(perm: Tuple[int, int, int]) -> npt.NDArray[np.floating[Any]]:
            return _separable_convolve_3d(
                image, kernel_arrays[perm[0]], kernel_arrays[perm[1]], kernel_arrays[perm[2]], mode
            )

        if use_parallel:
            with ThreadPoolExecutor() as executor:
                computed = list(executor.map(lambda kp: _base(kp[1]), base_keys))
        else:
            computed = [_base(perm) for _key, perm in base_keys]
        base_cache = {key: resp for (key, _perm), resp in zip(base_keys, computed, strict=True)}

        for perm, flips in rotations:
            key = (kernel_names[perm[0]], kernel_names[perm[1]], kernel_names[perm[2]])
            sign = 1
            for i, do_flip in enumerate(flips):
                if do_flip and antisym[perm[i]]:
                    sign = -sign
            base = base_cache[key]
            signed = base if sign > 0 else -base
            if result is None:
                result = signed.astype(np.float64) if pooling == "average" else signed.copy()
            elif pooling == "max":
                np.maximum(result, signed, out=result)
            elif pooling == "average":
                result += signed
            else:  # "min"
                np.minimum(result, signed, out=result)

        # Finalize average pooling
        if pooling == "average" and result is not None:
            result /= len(rotations)
    else:
        # Non-rotation-invariant: single separable convolution
        if source_mask is not None:
            result, valid_mask = _normalized_separable_convolve_3d(
                image, source_mask, g1, g2, g3, mode
            )
        else:
            result = _separable_convolve_3d(image, g1, g2, g3, mode)
            valid_mask = None

    # Compute energy image if requested
    if compute_energy:
        if result is None:  # pragma: no cover
            raise RuntimeError("Result should not be None")

        # Energy = mean of absolute values over δ neighborhood, i.e. uniform_filter
        # on |result|. Accumulate in float64: scipy's running moving-sum otherwise
        # drifts in float32 over long axes. Cast the result back to float32.
        abs_result = np.abs(result).astype(np.float64, copy=False)
        energy_support = 2 * energy_distance + 1
        result = uniform_filter(abs_result, size=energy_support, mode=mode).astype(np.float32)

    if result is None:  # pragma: no cover
        raise RuntimeError("Result should not be None")

    if source_mask is not None and valid_mask is not None:
        return result, valid_mask
    return result  # type: ignore[no-any-return]

pictologics.filters.gabor_filter(image, sigma_mm, lambda_mm, gamma=1.0, theta=0.0, spacing_mm=1.0, boundary=BoundaryCondition.ZERO, rotation_invariant=False, delta_theta=None, pooling='average', average_over_planes=False, use_parallel=None, source_mask=None)

Apply 2D Gabor filter to 3D image (IBSI code: Q88H).

The Gabor filter is applied in the axial plane (k1, k2) and optionally averaged over orthogonal planes. Per IBSI 2 Eq. 9.

Parameters:

Name Type Description Default
image NDArray[floating[Any]]

3D input image array

 required
sigma_mm float

Standard deviation of Gaussian envelope in mm (41LN)

 required
lambda_mm float

Wavelength in mm (S4N6)

 required
gamma float

Spatial aspect ratio (GDR5), typically 0.5 to 2.0

 1.0
theta float

Orientation angle in radians (FQER), clockwise in (k1,k2)

 0.0
spacing_mm Union[float, Tuple[float, float, float]]

Voxel spacing in mm (scalar or per-axis tuple). Each plane's kernel is built from that plane's own two in-plane axis spacings, so anisotropic spacing (including anisotropic z, relevant when average_over_planes=True) is handled correctly rather than approximated from a single axis.

 1.0
boundary Union[BoundaryCondition, str]

Boundary condition for padding (GBYQ)

 ZERO
rotation_invariant bool

If True, average over orientations

 False
delta_theta Optional[float]

Orientation step for rotation invariance (XTGK)

 None
pooling str

Pooling method ("average", "max", "min")

 'average'
average_over_planes bool

If True, average 2D responses over 3 orthogonal planes

 False
use_parallel Union[bool, None]

If True, process slices in parallel. If None (default), auto-enables for images > ~46³ voxels.

 None
source_mask Optional[NDArray[bool_]]

Optional boolean mask where True = valid voxel. When provided, zeros out invalid (sentinel) voxels before FFT-based convolution to prevent contamination.

 None

Returns:

Type Description
NDArray[floating[Any]]

Response map (modulus of complex response)

Raises:

Type Description
ValueError

If pooling is not "max", "average", or "min", or if rotation_invariant=True is set without providing delta_theta.
RuntimeError

Defensive check raised if plane averaging fails to produce a result; not expected to occur in normal use.
Example

Apply Gabor filter with rotation invariance over orthogonal planes:

import numpy as np
from pictologics.filters import gabor_filter

# Create dummy 3D image
image = np.random.rand(50, 50, 50)

# Apply filter
response = gabor_filter(
    image,
    sigma_mm=10.0,
    lambda_mm=4.0,
    gamma=0.5,
    rotation_invariant=True,
    delta_theta=0.7853981633974483,  # pi/4
    average_over_planes=True
)
Note
  • Returns modulus |h| = |g ⊗ f| for feature extraction
  • 2D filter applied slice-by-slice, then optionally over planes
  • Uses single complex FFT convolution for ~2x speedup
  • Each plane's kernel uses that plane's own two in-plane spacings. When they are equal (the isotropic-in-plane case, including the default axial-only plane under typical (x, y, z) spacing with x == y), the kernel is built on a voxel-unit grid. When they differ, the kernel is built on a physical-coordinate (mm) grid with a per-axis radius, giving a rectangular kernel that is physically correct rather than warning and guessing.
Source code in pictologics/filters/gabor.py 
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
def gabor_filter(
    image: npt.NDArray[np.floating[Any]],
    sigma_mm: float,
    lambda_mm: float,
    gamma: float = 1.0,
    theta: float = 0.0,
    spacing_mm: Union[float, Tuple[float, float, float]] = 1.0,
    boundary: Union[BoundaryCondition, str] = BoundaryCondition.ZERO,
    rotation_invariant: bool = False,
    delta_theta: Optional[float] = None,
    pooling: str = "average",
    average_over_planes: bool = False,
    use_parallel: Union[bool, None] = None,
    source_mask: Optional[npt.NDArray[np.bool_]] = None,
) -> npt.NDArray[np.floating[Any]]:
    """
    Apply 2D Gabor filter to 3D image (IBSI code: Q88H).

    The Gabor filter is applied in the axial plane (k1, k2) and optionally
    averaged over orthogonal planes. Per IBSI 2 Eq. 9.

    Args:
        image: 3D input image array
        sigma_mm: Standard deviation of Gaussian envelope in mm (41LN)
        lambda_mm: Wavelength in mm (S4N6)
        gamma: Spatial aspect ratio (GDR5), typically 0.5 to 2.0
        theta: Orientation angle in radians (FQER), clockwise in (k1,k2)
        spacing_mm: Voxel spacing in mm (scalar or per-axis tuple). Each
            plane's kernel is built from that plane's own two in-plane
            axis spacings, so anisotropic spacing (including anisotropic
            z, relevant when `average_over_planes=True`) is handled
            correctly rather than approximated from a single axis.
        boundary: Boundary condition for padding (GBYQ)
        rotation_invariant: If True, average over orientations
        delta_theta: Orientation step for rotation invariance (XTGK)
        pooling: Pooling method ("average", "max", "min")
        average_over_planes: If True, average 2D responses over 3 orthogonal planes
        use_parallel: If True, process slices in parallel. If None (default),
            auto-enables for images > ~46³ voxels.
        source_mask: Optional boolean mask where True = valid voxel.
            When provided, zeros out invalid (sentinel) voxels before
            FFT-based convolution to prevent contamination.

    Returns:
        Response map (modulus of complex response)

    Raises:
        ValueError: If `pooling` is not "max", "average", or "min", or if
            `rotation_invariant=True` is set without providing `delta_theta`.
        RuntimeError: Defensive check raised if plane averaging fails to
            produce a result; not expected to occur in normal use.

    Example:
        Apply Gabor filter with rotation invariance over orthogonal planes:

        ```python
        import numpy as np
        from pictologics.filters import gabor_filter

        # Create dummy 3D image
        image = np.random.rand(50, 50, 50)

        # Apply filter
        response = gabor_filter(
            image,
            sigma_mm=10.0,
            lambda_mm=4.0,
            gamma=0.5,
            rotation_invariant=True,
            delta_theta=0.7853981633974483,  # pi/4
            average_over_planes=True
        )
        ```

    Note:
        - Returns modulus |h| = |g ⊗ f| for feature extraction
        - 2D filter applied slice-by-slice, then optionally over planes
        - Uses single complex FFT convolution for ~2x speedup
        - Each plane's kernel uses that plane's own two in-plane spacings.
          When they are equal (the isotropic-in-plane case, including the
          default axial-only plane under typical (x, y, z) spacing with
          x == y), the kernel is built on a voxel-unit grid. When they
          differ, the kernel is built on a physical-coordinate (mm) grid
          with a per-axis radius, giving a rectangular kernel that is
          physically correct rather than warning and guessing.
    """
    # Convert to float32
    image = ensure_float32(image)

    # Apply source_mask preprocessing (zero out invalid voxels for FFT-based filter)
    if source_mask is not None:
        image = _prepare_masked_image(image, source_mask)

    # Handle spacing. The mm -> voxel/physical conversion is deferred to
    # _apply_gabor_to_plane, which is per-plane: each plane's in-plane axes
    # (and therefore in-plane spacings) depend on plane_axis.
    if isinstance(spacing_mm, (int, float)):
        spacing_mm = (float(spacing_mm),) * 3

    # Handle boundary
    if isinstance(boundary, str):
        boundary = BoundaryCondition[boundary.upper()]
    mode = get_scipy_mode(boundary)

    # Validate pooling parameter early
    valid_poolings = ("max", "average", "min")
    if pooling not in valid_poolings:
        raise ValueError(f"Unknown pooling: {pooling}. Must be one of {valid_poolings}")

    # Auto-detect parallel mode based on image size
    if use_parallel is None:
        use_parallel = image.size > _PARALLEL_THRESHOLD

    if rotation_invariant:
        if delta_theta is None:
            raise ValueError(
                "rotation_invariant=True requires delta_theta (the orientation step in radians)"
            )
        # Generate orientations from 0 to 2π
        n_orientations = int(np.ceil(2 * np.pi / delta_theta))
        thetas = [i * delta_theta for i in range(n_orientations)]
        # The Gabor response modulus is π-periodic in theta: kernel(θ+π) = conj(kernel(θ))
        # and the image is real, so |response| is identical for θ and θ+π. When the
        # orientation set is closed under +π (n even and spans exactly 2π), the second
        # half duplicates the first; drop it (max/min/average pooling are unchanged).
        if n_orientations % 2 == 0 and abs(n_orientations * delta_theta - 2 * np.pi) < 1e-9:
            thetas = thetas[: n_orientations // 2]
    else:
        thetas = [theta]

    if average_over_planes:
        # Apply to all 3 orthogonal planes and average with in-place aggregation
        result: npt.NDArray[np.floating[Any]] | None = None
        for plane_axis in range(3):
            plane_response = _apply_gabor_to_plane(
                image,
                sigma_mm,
                lambda_mm,
                gamma,
                thetas,
                plane_axis,
                spacing_mm,
                mode,
                pooling,
                use_parallel,
            )
            if result is None:
                result = plane_response.astype(np.float64)
            else:
                result += plane_response

        if result is None:  # pragma: no cover
            raise RuntimeError("Result should not be None after plane loop")

        return (result / 3.0).astype(np.float32)  # type: ignore[union-attr]
    else:
        # Apply only to axial plane (axis 2 = k3 slices)
        return _apply_gabor_to_plane(
            image,
            sigma_mm,
            lambda_mm,
            gamma,
            thetas,
            plane_axis=2,
            spacing_mm=spacing_mm,
            mode=mode,
            pooling=pooling,
            use_parallel=use_parallel,
        )

pictologics.filters.wavelet_transform(image, wavelet='db2', level=1, decomposition='LHL', boundary=BoundaryCondition.ZERO, rotation_invariant=False, pooling='average', use_parallel=None, source_mask=None)

Apply 3D separable wavelet transform (undecimated/stationary).

Uses the à trous algorithm for undecimated wavelet decomposition. The transform is translation-invariant (unlike decimated transform).

Supported wavelets
  • "haar" (UOUE): Haar wavelet
  • "db2", "db3": Daubechies wavelets
  • "coif1": Coiflet wavelet

Parameters:

Name Type Description Default
image NDArray[floating[Any]]

3D input image array

 required
wavelet str

Wavelet name (e.g., "db2", "coif1", "haar")

 'db2'
level int

Decomposition level (GCEK)

 1
decomposition str

Which response map to return, e.g., "LHL", "HHH"

 'LHL'
boundary Union[BoundaryCondition, str]

Boundary condition for padding

 ZERO
rotation_invariant bool

If True, average over 24 rotations

 False
pooling str

Pooling method for rotation invariance

 'average'
use_parallel Union[bool, None]

If True, use parallel processing for rotation_invariant mode. If None (default), auto-enables for images > ~128³ voxels.

 None
source_mask Optional[NDArray[bool_]]

Optional boolean mask where True = valid voxel. When provided, zeros out invalid (sentinel) voxels before wavelet decomposition to prevent contamination.

 None

Returns:

Type Description
NDArray[floating[Any]]

Response map for the specified decomposition

Raises:

Type Description
ValueError

If rotation_invariant=True and pooling is not "max", "average", or "min".
Example

Apply Daubechies 2 wavelet transform at level 1, returning LHL coefficients:

import numpy as np
from pictologics.filters import wavelet_transform

# Create dummy 3D image
image = np.random.rand(50, 50, 50)

# Apply transform
response = wavelet_transform(
    image,
    wavelet="db2",
    level=1,
    decomposition="LHL"
)
Source code in pictologics/filters/wavelets.py 
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
def wavelet_transform(
    image: npt.NDArray[np.floating[Any]],
    wavelet: str = "db2",
    level: int = 1,
    decomposition: str = "LHL",
    boundary: Union[BoundaryCondition, str] = BoundaryCondition.ZERO,
    rotation_invariant: bool = False,
    pooling: str = "average",
    use_parallel: Union[bool, None] = None,
    source_mask: Optional[npt.NDArray[np.bool_]] = None,
) -> npt.NDArray[np.floating[Any]]:
    """
    Apply 3D separable wavelet transform (undecimated/stationary).

    Uses the à trous algorithm for undecimated wavelet decomposition.
    The transform is translation-invariant (unlike decimated transform).

    Supported wavelets:
        - "haar" (UOUE): Haar wavelet
        - "db2", "db3": Daubechies wavelets
        - "coif1": Coiflet wavelet

    Args:
        image: 3D input image array
        wavelet: Wavelet name (e.g., "db2", "coif1", "haar")
        level: Decomposition level (GCEK)
        decomposition: Which response map to return, e.g., "LHL", "HHH"
        boundary: Boundary condition for padding
        rotation_invariant: If True, average over 24 rotations
        pooling: Pooling method for rotation invariance
        use_parallel: If True, use parallel processing for rotation_invariant mode.
            If None (default), auto-enables for images > ~128³ voxels.
        source_mask: Optional boolean mask where True = valid voxel.
            When provided, zeros out invalid (sentinel) voxels before
            wavelet decomposition to prevent contamination.

    Returns:
        Response map for the specified decomposition

    Raises:
        ValueError: If `rotation_invariant=True` and `pooling` is not "max",
            "average", or "min".

    Example:
        Apply Daubechies 2 wavelet transform at level 1, returning LHL coefficients:

        ```python
        import numpy as np
        from pictologics.filters import wavelet_transform

        # Create dummy 3D image
        image = np.random.rand(50, 50, 50)

        # Apply transform
        response = wavelet_transform(
            image,
            wavelet="db2",
            level=1,
            decomposition="LHL"
        )
        ```
    """
    # Convert to float32
    image = ensure_float32(image)

    # Apply source_mask preprocessing (zero out invalid voxels)
    if source_mask is not None:
        image = _prepare_masked_image(image, source_mask)

    # Handle boundary
    if isinstance(boundary, str):
        boundary = BoundaryCondition[boundary.upper()]
    mode = get_scipy_mode(boundary)

    # Get wavelet filters
    w = pywt.Wavelet(wavelet)
    lo = np.array(w.dec_lo, dtype=np.float32)  # Low-pass decomposition filter
    hi = np.array(w.dec_hi, dtype=np.float32)  # High-pass decomposition filter

    # Auto-detect parallel mode based on image size
    if use_parallel is None:
        use_parallel = image.size > _PARALLEL_THRESHOLD

    if rotation_invariant:
        if pooling not in ("max", "average", "min"):
            raise ValueError(f"Unknown pooling: {pooling}")

        rotations = _get_rotation_perms()

        def apply_rotated_wavelet(
            rotation: Tuple[Tuple[int, int, int], Tuple[bool, bool, bool]],
        ) -> npt.NDArray[np.floating[Any]]:
            """Apply wavelet transform with rotated image."""
            perm, flips = rotation
            # Permute and flip image
            rotated = np.transpose(image, perm)
            for axis, flip in enumerate(flips):
                if flip:
                    rotated = np.flip(rotated, axis=axis)

            # Apply wavelet
            response = _apply_undecimated_wavelet_3d(rotated, lo, hi, level, decomposition, mode)

            # Undo rotation for response
            for axis, flip in enumerate(flips):
                if flip:
                    response = np.flip(response, axis=axis)
            inv_perm = tuple(np.argsort(perm))
            return np.transpose(response, inv_perm)

        result: npt.NDArray[np.floating[Any]] | None = None

        def _pool(response: npt.NDArray[np.floating[Any]]) -> None:
            nonlocal result
            if result is None:
                result = response.astype(np.float64) if pooling == "average" else response
            elif pooling == "max":
                np.maximum(result, response, out=result)
            elif pooling == "average":
                result += response
            else:  # "min"
                np.minimum(result, response, out=result)

        if use_parallel:
            with ThreadPoolExecutor() as executor:
                future_to_rot = {
                    executor.submit(apply_rotated_wavelet, rot): rot for rot in rotations
                }
                # Pool responses as they complete to avoid holding all 24 at once.
                for future in as_completed(future_to_rot):
                    _pool(future.result())
        else:
            # Sequential processing for small images
            for rotation in rotations:
                _pool(apply_rotated_wavelet(rotation))

        # Finalize average pooling
        if pooling == "average" and result is not None:
            result /= len(rotations)
        return result.astype(np.float32)  # type: ignore[union-attr]
    else:
        return _apply_undecimated_wavelet_3d(image, lo, hi, level, decomposition, mode)

pictologics.filters.simoncelli_wavelet(image, level=1, boundary=BoundaryCondition.PERIODIC, source_mask=None)

Apply Simoncelli non-separable wavelet (IBSI code: PRT7).

The Simoncelli wavelet is isotropic (spherically symmetric) and implemented in the Fourier domain. Per IBSI 2 Eq. 27.

For decomposition level N, the frequency band is scaled by j = N-1: - Level 1 (j=0): band [π/4, π] (highest frequencies) - Level 2 (j=1): band [π/8, π/2] - Level 3 (j=2): band [π/16, π/4]

Parameters:

Name Type Description Default
image NDArray[floating[Any]]

3D input image array

 required
level int

Decomposition level (1 = highest frequency band)

 1
boundary Union[BoundaryCondition, str]

Boundary condition. The filter is inherently periodic (FFT-based), so BoundaryCondition.PERIODIC (the default) runs it directly on image. Any other condition is approximated via pad-filter-crop (see _apply_with_boundary_padding and _simoncelli_pad_width).

 PERIODIC
source_mask Optional[NDArray[bool_]]

Optional boolean mask where True = valid voxel

 None

Returns:

Type Description
NDArray[floating[Any]]

Band-pass response map (B map) for the specified level

Raises:

Type Description
ValueError

If boundary is a string that is not a valid BoundaryCondition member name.
Example

Apply first-level Simoncelli wavelet (highest frequency band):

import numpy as np
from pictologics.filters import simoncelli_wavelet

# Create dummy 3D image
image = np.random.rand(50, 50, 50)

# Apply wavelet
response = simoncelli_wavelet(image, level=1)
Source code in pictologics/filters/wavelets.py 
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
def simoncelli_wavelet(
    image: npt.NDArray[np.floating[Any]],
    level: int = 1,
    boundary: Union[BoundaryCondition, str] = BoundaryCondition.PERIODIC,
    source_mask: Optional[npt.NDArray[np.bool_]] = None,
) -> npt.NDArray[np.floating[Any]]:
    """
    Apply Simoncelli non-separable wavelet (IBSI code: PRT7).

    The Simoncelli wavelet is isotropic (spherically symmetric) and
    implemented in the Fourier domain. Per IBSI 2 Eq. 27.

    For decomposition level N, the frequency band is scaled by j = N-1:
        - Level 1 (j=0): band [π/4, π] (highest frequencies)
        - Level 2 (j=1): band [π/8, π/2]
        - Level 3 (j=2): band [π/16, π/4]

    Args:
        image: 3D input image array
        level: Decomposition level (1 = highest frequency band)
        boundary: Boundary condition. The filter is inherently periodic (FFT-based),
            so `BoundaryCondition.PERIODIC` (the default) runs it directly on
            `image`. Any other condition is approximated via pad-filter-crop (see
            `_apply_with_boundary_padding` and `_simoncelli_pad_width`).
        source_mask: Optional boolean mask where True = valid voxel

    Returns:
        Band-pass response map (B map) for the specified level

    Raises:
        ValueError: If `boundary` is a string that is not a valid
            `BoundaryCondition` member name.

    Example:
        Apply first-level Simoncelli wavelet (highest frequency band):

        ```python
        import numpy as np
        from pictologics.filters import simoncelli_wavelet

        # Create dummy 3D image
        image = np.random.rand(50, 50, 50)

        # Apply wavelet
        response = simoncelli_wavelet(image, level=1)
        ```
    """
    boundary = resolve_boundary(boundary)

    # Convert to float32
    image = ensure_float32(image)

    # Apply source_mask preprocessing (zero out invalid voxels for FFT-based filter)
    if source_mask is not None:
        image = _prepare_masked_image(image, source_mask)

    def _core(arr: npt.NDArray[np.floating[Any]]) -> npt.NDArray[np.floating[Any]]:
        shape = tuple(arr.shape)
        ndim = arr.ndim

        # Transfer function depends only on (shape, level) — never on image values or
        # the source mask — so it is built once and cached (see _simoncelli_transfer).
        g_sim = _simoncelli_transfer(shape, level)

        # Apply filter in frequency domain using full FFT (full FFT required because
        # the centered grid is non-symmetric for even N). scipy.fft with workers=-1
        # is multithreaded and matches np.fft to float32 precision.
        axes = tuple(range(ndim))
        F = scipy.fft.fftn(arr, workers=-1)
        response = scipy.fft.ifftn(F * g_sim, s=shape, axes=axes, workers=-1)

        return cast(npt.NDArray[np.floating[Any]], np.real(response).astype(np.float32))

    return _apply_with_boundary_padding(_core, image, boundary, _simoncelli_pad_width(level))

pictologics.filters.riesz_transform(image, order, boundary=BoundaryCondition.PERIODIC, source_mask=None)

Apply Riesz transform (IBSI code: AYRS).

The Riesz transform computes higher-order all-pass image derivatives in the Fourier domain. Per IBSI 2 Eq. 34.

Parameters:

Name Type Description Default
image NDArray[floating[Any]]

3D input image array

 required
order Tuple[int, ...]

Tuple (l1, l2, l3) specifying derivative order per axis e.g., (1,0,0) = first-order along k1 (gradient-like) (2,0,0), (1,1,0), (0,2,0) = second-order (Hessian-like)

 required
boundary Union[BoundaryCondition, str]

Boundary condition. The filter is inherently periodic (FFT-based), so BoundaryCondition.PERIODIC (the default) runs it directly on image. Any other condition is approximated via pad-filter-crop (see _apply_with_boundary_padding and _RIESZ_BASE_PAD).

 PERIODIC
source_mask Optional[NDArray[bool_]]

Optional boolean mask where True = valid voxel. When provided, zeros out invalid (sentinel) voxels before FFT-based transform to prevent contamination.

 None

Returns:

Type Description
NDArray[floating[Any]]

Riesz-transformed image (real part)

Raises:

Type Description
ValueError

If order sums to 0 (i.e. every component is 0), which would correspond to a zero-order (identity) transform, or if boundary is a string that is not a valid BoundaryCondition member name.
Example

Compute first-order Riesz transform along the k1 axis:

import numpy as np
from pictologics.filters import riesz_transform

# Create dummy 3D image
image = np.random.rand(50, 50, 50)

# Apply transform (gradient-like along axis 0)
response = riesz_transform(image, order=(1, 0, 0))
Note
  • First-order Riesz components form the image gradient
  • Second-order Riesz components form the image Hessian
  • All-pass: doesn't amplify high frequencies like regular derivatives
Source code in pictologics/filters/riesz.py 
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
def riesz_transform(
    image: npt.NDArray[np.floating[Any]],
    order: Tuple[int, ...],
    boundary: Union[BoundaryCondition, str] = BoundaryCondition.PERIODIC,
    source_mask: Optional[npt.NDArray[np.bool_]] = None,
) -> npt.NDArray[np.floating[Any]]:
    """
    Apply Riesz transform (IBSI code: AYRS).

    The Riesz transform computes higher-order all-pass image derivatives
    in the Fourier domain. Per IBSI 2 Eq. 34.

    Args:
        image: 3D input image array
        order: Tuple (l1, l2, l3) specifying derivative order per axis
               e.g., (1,0,0) = first-order along k1 (gradient-like)
                     (2,0,0), (1,1,0), (0,2,0) = second-order (Hessian-like)
        boundary: Boundary condition. The filter is inherently periodic (FFT-based),
            so `BoundaryCondition.PERIODIC` (the default) runs it directly on
            `image`. Any other condition is approximated via pad-filter-crop (see
            `_apply_with_boundary_padding` and `_RIESZ_BASE_PAD`).
        source_mask: Optional boolean mask where True = valid voxel.
            When provided, zeros out invalid (sentinel) voxels before
            FFT-based transform to prevent contamination.

    Returns:
        Riesz-transformed image (real part)

    Raises:
        ValueError: If `order` sums to 0 (i.e. every component is 0), which
            would correspond to a zero-order (identity) transform, or if
            `boundary` is a string that is not a valid `BoundaryCondition`
            member name.

    Example:
        Compute first-order Riesz transform along the k1 axis:

        ```python
        import numpy as np
        from pictologics.filters import riesz_transform

        # Create dummy 3D image
        image = np.random.rand(50, 50, 50)

        # Apply transform (gradient-like along axis 0)
        response = riesz_transform(image, order=(1, 0, 0))
        ```

    Note:
        - First-order Riesz components form the image gradient
        - Second-order Riesz components form the image Hessian
        - All-pass: doesn't amplify high frequencies like regular derivatives
    """
    boundary = resolve_boundary(boundary)

    # Convert to float32
    image = ensure_float32(image)

    # Apply source_mask preprocessing (zero out invalid voxels for FFT-based filter)
    if source_mask is not None:
        image = _prepare_masked_image(image, source_mask)

    L = sum(order)  # Total order

    if L == 0:
        raise ValueError("At least one order component must be > 0")

    # Coerce order to a tuple first so a list-typed order (e.g. from a YAML/JSON
    # pipeline config) stays hashable for the transfer-function cache key.
    order = tuple(order)

    def _core(arr: npt.NDArray[np.floating[Any]]) -> npt.NDArray[np.floating[Any]]:
        shape = tuple(arr.shape)
        ndim = arr.ndim

        # Transfer function depends only on (shape, order) — never on image values
        # or the source mask — so it is built once and cached (see _riesz_transfer).
        transfer = _riesz_transfer(shape, order)

        # Apply in frequency domain using Real FFT. scipy.fft (multithreaded via
        # workers=-1) is several times faster than the single-threaded np.fft and
        # matches it to float32 precision.
        axes = tuple(range(ndim))
        F = scipy.fft.rfftn(arr, workers=-1)

        # F has shape (N1, N2, N3//2 + 1); transfer is broadcastable to it.
        response = scipy.fft.irfftn(F * transfer, s=shape, axes=axes, workers=-1)

        return cast(npt.NDArray[np.floating[Any]], response.astype(np.float32))

    return _apply_with_boundary_padding(_core, image, boundary, _RIESZ_BASE_PAD)

pictologics.filters.riesz_log(image, sigma_mm, spacing_mm=1.0, order=(1, 0, 0), truncate=4.0, boundary=BoundaryCondition.PERIODIC, source_mask=None)

Apply Riesz transform to LoG-filtered image.

Combines multi-scale analysis (LoG) with directional analysis (Riesz). First applies LoG filtering, then applies Riesz transform.

Parameters:

Name Type Description Default
image NDArray[floating[Any]]

3D input image array

 required
sigma_mm float

LoG scale in mm

 required
spacing_mm Union[float, Tuple[float, float, float]]

Voxel spacing in mm

 1.0
order Tuple[int, ...]

Riesz order tuple (l1, l2, l3)

 (1, 0, 0)
truncate float

LoG truncation parameter

 4.0
boundary Union[BoundaryCondition, str]

Boundary condition for the whole LoG-then-Riesz chain. The default BoundaryCondition.PERIODIC reproduces today's exact behaviour: the internal LoG call keeps its own default (ZERO padding) and the Riesz stage stays periodic, with no outer padding at all. Any other condition pads image once (see _apply_with_boundary_padding and _riesz_log_pad_width), runs the LoG-then-Riesz chain on the padded array, and crops back — and is also forwarded to the internal LoG call so its own edge handling matches the requested condition instead of silently staying at ZERO.

 PERIODIC
source_mask Optional[NDArray[bool_]]

Optional boolean mask where True = valid voxel. Because source_mask shares image's (unpadded) shape, the mid-chain re-zeroing this function otherwise performs before the Riesz stage is only applied when no padding occurs (the default PERIODIC case); for any other boundary, the mask is instead applied once to the final, already-cropped response.

 None

Returns:

Type Description
NDArray[floating[Any]]

Riesz-transformed LoG response

Raises:

Type Description
ValueError

If boundary is a string that is not a valid BoundaryCondition member name.
Example

Compute first-order Riesz transform of LoG-filtered image at 5mm scale:

import numpy as np
from pictologics.filters import riesz_log

# Create dummy 3D image
image = np.random.rand(50, 50, 50)

# Apply filter
response = riesz_log(
    image,
    sigma_mm=5.0,
    spacing_mm=(2.0, 2.0, 2.0),
    order=(1, 0, 0)
)
Source code in pictologics/filters/riesz.py 
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
def riesz_log(
    image: npt.NDArray[np.floating[Any]],
    sigma_mm: float,
    spacing_mm: Union[float, Tuple[float, float, float]] = 1.0,
    order: Tuple[int, ...] = (1, 0, 0),
    truncate: float = 4.0,
    boundary: Union[BoundaryCondition, str] = BoundaryCondition.PERIODIC,
    source_mask: Optional[npt.NDArray[np.bool_]] = None,
) -> npt.NDArray[np.floating[Any]]:
    """
    Apply Riesz transform to LoG-filtered image.

    Combines multi-scale analysis (LoG) with directional analysis (Riesz).
    First applies LoG filtering, then applies Riesz transform.

    Args:
        image: 3D input image array
        sigma_mm: LoG scale in mm
        spacing_mm: Voxel spacing in mm
        order: Riesz order tuple (l1, l2, l3)
        truncate: LoG truncation parameter
        boundary: Boundary condition for the whole LoG-then-Riesz chain. The
            default `BoundaryCondition.PERIODIC` reproduces today's exact
            behaviour: the internal LoG call keeps its own default (ZERO padding)
            and the Riesz stage stays periodic, with no outer padding at all. Any
            other condition pads `image` once (see `_apply_with_boundary_padding`
            and `_riesz_log_pad_width`), runs the LoG-then-Riesz chain on the
            padded array, and crops back — and is *also* forwarded to the
            internal LoG call so its own edge handling matches the requested
            condition instead of silently staying at ZERO.
        source_mask: Optional boolean mask where True = valid voxel. Because
            `source_mask` shares `image`'s (unpadded) shape, the mid-chain
            re-zeroing this function otherwise performs before the Riesz stage is
            only applied when no padding occurs (the default `PERIODIC` case);
            for any other boundary, the mask is instead applied once to the final,
            already-cropped response.

    Returns:
        Riesz-transformed LoG response

    Raises:
        ValueError: If `boundary` is a string that is not a valid
            `BoundaryCondition` member name.

    Example:
        Compute first-order Riesz transform of LoG-filtered image at 5mm scale:

        ```python
        import numpy as np
        from pictologics.filters import riesz_log

        # Create dummy 3D image
        image = np.random.rand(50, 50, 50)

        # Apply filter
        response = riesz_log(
            image,
            sigma_mm=5.0,
            spacing_mm=(2.0, 2.0, 2.0),
            order=(1, 0, 0)
        )
        ```
    """
    from .log import laplacian_of_gaussian

    boundary = resolve_boundary(boundary)

    def _core(arr: npt.NDArray[np.floating[Any]]) -> npt.NDArray[np.floating[Any]]:
        # `_core` runs on `image` unchanged when boundary is PERIODIC (the default,
        # no padding), and on a *padded* array otherwise. `source_mask` always has
        # `image`'s original, unpadded shape, so it can only be forwarded to the
        # internal calls below in the PERIODIC case; the non-PERIODIC case masks
        # once, after cropping, below.
        if boundary is BoundaryCondition.PERIODIC:
            log_response = laplacian_of_gaussian(
                arr,
                sigma_mm=sigma_mm,
                spacing_mm=spacing_mm,
                truncate=truncate,
                source_mask=source_mask,
            )
            mask = source_mask
        else:
            # Also forward `boundary` here so LoG's own edge handling (ZERO by
            # default) matches the requested condition instead of silently
            # staying at ZERO.
            log_response = laplacian_of_gaussian(
                arr,
                sigma_mm=sigma_mm,
                spacing_mm=spacing_mm,
                truncate=truncate,
                boundary=boundary,
                source_mask=None,
            )
            mask = None

        # Handle tuple return from LoG if source_mask was used
        if isinstance(log_response, tuple):
            log_response = log_response[0]

        # Then apply Riesz transform. We pass the mask again (PERIODIC case only)
        # to enforce zeroing of invalid regions (though LoG normalized convolution
        # might have filled them, Riesz is global). The Riesz stage keeps its own
        # PERIODIC default: the outer pad-filter-crop below already accounts for
        # the boundary once for the whole chain.
        return riesz_transform(log_response, order=order, source_mask=mask)

    pad_width = _riesz_log_pad_width(sigma_mm, spacing_mm, truncate)
    result = _apply_with_boundary_padding(_core, image, boundary, pad_width)

    if source_mask is not None and boundary is not BoundaryCondition.PERIODIC:
        result = _prepare_masked_image(result, source_mask)

    return result

pictologics.filters.riesz_simoncelli(image, level=1, order=(1, 0, 0), boundary=BoundaryCondition.PERIODIC, source_mask=None)

Apply Riesz transform to Simoncelli wavelet-filtered image.

Combines isotropic multi-scale analysis (Simoncelli) with directional analysis (Riesz) for rotation-invariant directional features.

Parameters:

Name Type Description Default
image NDArray[floating[Any]]

3D input image array

 required
level int

Simoncelli decomposition level

 1
order Tuple[int, ...]

Riesz order tuple (l1, l2, l3)

 (1, 0, 0)
boundary Union[BoundaryCondition, str]

Boundary condition for the whole Simoncelli-then-Riesz chain. The default BoundaryCondition.PERIODIC reproduces today's exact behaviour (both stages run periodically, with no outer padding at all). Any other condition pads image once — covering both FFT stages with a single, uniform boundary treatment — runs the chain on the padded array (each stage keeping its own PERIODIC default), and crops back (see _apply_with_boundary_padding).

 PERIODIC
source_mask Optional[NDArray[bool_]]

Optional boolean mask where True = valid voxel. Because source_mask shares image's (unpadded) shape, the mid-chain re-zeroing this function otherwise performs before the Riesz stage is only applied when no padding occurs (the default PERIODIC case); for any other boundary, the mask is instead applied once to the final, already-cropped response.

 None

Returns:

Type Description
NDArray[floating[Any]]

Riesz-transformed Simoncelli response

Raises:

Type Description
ValueError

If boundary is a string that is not a valid BoundaryCondition member name.
Example

Compute second-order Riesz transform (Hessian-like) of Simoncelli level 2:

import numpy as np
from pictologics.filters import riesz_simoncelli

# Create dummy 3D image
image = np.random.rand(50, 50, 50)

# Apply filter
response = riesz_simoncelli(
    image,
    level=2,
    order=(2, 0, 0)
)
Source code in pictologics/filters/riesz.py 
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
def riesz_simoncelli(
    image: npt.NDArray[np.floating[Any]],
    level: int = 1,
    order: Tuple[int, ...] = (1, 0, 0),
    boundary: Union[BoundaryCondition, str] = BoundaryCondition.PERIODIC,
    source_mask: Optional[npt.NDArray[np.bool_]] = None,
) -> npt.NDArray[np.floating[Any]]:
    """
    Apply Riesz transform to Simoncelli wavelet-filtered image.

    Combines isotropic multi-scale analysis (Simoncelli) with
    directional analysis (Riesz) for rotation-invariant directional features.

    Args:
        image: 3D input image array
        level: Simoncelli decomposition level
        order: Riesz order tuple (l1, l2, l3)
        boundary: Boundary condition for the whole Simoncelli-then-Riesz chain.
            The default `BoundaryCondition.PERIODIC` reproduces today's exact
            behaviour (both stages run periodically, with no outer padding at
            all). Any other condition pads `image` once — covering both FFT
            stages with a single, uniform boundary treatment — runs the chain on
            the padded array (each stage keeping its own PERIODIC default), and
            crops back (see `_apply_with_boundary_padding`).
        source_mask: Optional boolean mask where True = valid voxel. Because
            `source_mask` shares `image`'s (unpadded) shape, the mid-chain
            re-zeroing this function otherwise performs before the Riesz stage is
            only applied when no padding occurs (the default `PERIODIC` case);
            for any other boundary, the mask is instead applied once to the final,
            already-cropped response.

    Returns:
        Riesz-transformed Simoncelli response

    Raises:
        ValueError: If `boundary` is a string that is not a valid
            `BoundaryCondition` member name.

    Example:
        Compute second-order Riesz transform (Hessian-like) of Simoncelli level 2:

        ```python
        import numpy as np
        from pictologics.filters import riesz_simoncelli

        # Create dummy 3D image
        image = np.random.rand(50, 50, 50)

        # Apply filter
        response = riesz_simoncelli(
            image,
            level=2,
            order=(2, 0, 0)
        )
        ```
    """
    from .wavelets import _simoncelli_pad_width, simoncelli_wavelet

    boundary = resolve_boundary(boundary)

    # Preprocess once: float32 conversion + source mask zeroing
    image = ensure_float32(image)
    if source_mask is not None:
        image = _prepare_masked_image(image, source_mask)

    def _core(arr: npt.NDArray[np.floating[Any]]) -> npt.NDArray[np.floating[Any]]:
        # Apply Simoncelli wavelet (already preprocessed, skip redundant work)
        sim_response = simoncelli_wavelet(arr, level=level)

        # Re-apply source_mask (PERIODIC case only, see docstring): Simoncelli's
        # global FFT spreads energy back into the invalid regions, and the Riesz
        # transform is likewise global, so re-zero before it (mirrors riesz_log).
        mask = source_mask if boundary is BoundaryCondition.PERIODIC else None
        return riesz_transform(sim_response, order=order, source_mask=mask)

    pad_width = _simoncelli_pad_width(level) + _RIESZ_BASE_PAD
    result = _apply_with_boundary_padding(_core, image, boundary, pad_width)

    if source_mask is not None and boundary is not BoundaryCondition.PERIODIC:
        result = _prepare_masked_image(result, source_mask)

    return result