Filters API

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"``, ``"gaussian"``, ``"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

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,
    padding_value: float = 0.0,
) -> 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.
        padding_value: The constant of constant value padding (Z3VE), with the ZERO
            (constant) boundary. Default 0.

    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
        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_mask = _bool_mask(source_mask)  # a mask of another type: != 0
    # Validate support
    if support < 1 or support % 2 == 0:
        raise ValueError(f"Support must be an odd positive integer, got {support}")

    # The passes run in float32. The normalized convolution of a source mask keeps the
    # image type: it divides two filtered sums.
    image = ensure_float32(image) if source_mask is not None else as_float32(image)

    # Handle string boundary condition
    boundary = resolve_boundary(boundary)

    problem = _padding_value_problem(boundary, padding_value)
    if problem:
        raise ValueError(problem)
    if padding_value and source_mask is not None:
        return _constant_padded(
            mean_filter, image, source_mask, padding_value, support // 2, support=support, boundary=boundary
        )  # fmt: skip
    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:
        # A float32 response: each pass sums its lines in double and writes float32. The
        # passes read the padding value beyond the image (no padded copy).
        cval = float(image.dtype.type(padding_value))
        return _uniform_filter(image, support, mode, cval=cval)

def gaussian_filter(
    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,
    padding_value: float = 0.0,
) -> Union[
    npt.NDArray[np.floating[Any]],
    Tuple[npt.NDArray[np.floating[Any]], npt.NDArray[np.bool_]],
]:
    """
    Apply a 3D Gaussian filter (IBSI code: 8BC3): a smoothing low-pass filter.

    Args:
        image: 3D input image array
        sigma_mm: Scale σ* in mm (41LN)
        spacing_mm: Voxel spacing in mm (scalar or a value for each axis)
        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,
            normalized convolution leaves out the invalid voxels.
        padding_value: The constant of constant value padding (Z3VE), with the ZERO
            (constant) boundary. Default 0.

    Returns:
        If source_mask is None: Response map (float32) with the dimensions of the input.
        If source_mask is given: (response map, output valid mask).

    Raises:
        ValueError: If `padding_value` is not 0 with another boundary than ZERO.

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

        image = np.random.rand(50, 50, 50)
        smoothed = gaussian_filter(image, sigma_mm=2.0, spacing_mm=(0.8, 0.8, 2.0))
        ```

    Note:
        - σ is converted from mm to voxels: σ_voxels = σ_mm / spacing_mm
        - The kernel of each axis has the radius ⌊truncate × σ_voxels + 0.5⌋ and sums to 1
    """
    source_mask = _bool_mask(source_mask)  # a mask of another type: != 0
    # The passes run in float32. The normalized convolution of a source mask keeps the
    # image type: it divides two filtered sums.
    image = ensure_float32(image) if source_mask is not None else as_float32(image)
    if isinstance(spacing_mm, (int, float)):
        spacing_mm = (float(spacing_mm),) * 3
    sigma_voxels = tuple(sigma_mm / s for s in spacing_mm)
    boundary = resolve_boundary(boundary)
    problem = _padding_value_problem(boundary, padding_value)
    if problem:
        raise ValueError(problem)
    if padding_value and source_mask is not None:
        reach = tuple(int(truncate * s + 0.5) for s in sigma_voxels)
        return _constant_padded(
            gaussian_filter, image, source_mask, padding_value, reach,
            sigma_mm=sigma_mm, spacing_mm=spacing_mm, truncate=truncate, boundary=boundary,
        )  # fmt: skip
    mode = get_scipy_mode(boundary)
    if source_mask is not None:
        return _normalized_gaussian(image, source_mask, sigma_voxels, mode, truncate)
    # A float32 response: each pass sums its lines in double and writes float32. The passes
    # read the padding value beyond the image (no padded copy).
    cval = float(image.dtype.type(padding_value))
    return _gaussian_filter(image, sigma_voxels, mode, truncate, cval=cval)

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,
    padding_value: float = 0.0,
) -> 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.
        padding_value: The constant of constant value padding (Z3VE), with the ZERO
            (constant) boundary. Default 0.

    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 `sigma_mm`, an entry of `spacing_mm` or `truncate` is not a finite
            positive number, or if `padding_value` is not 0 with another boundary than
            ZERO.

    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
        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)
        - A float64 image of 131,072 voxels (about 51³) or more runs as one FFT
          convolution with the same truncated kernels and boundary: about 2 times faster
          than the separable passes, with values within about 1e-15 of the largest
          response value.
        - A smaller float64 image runs the passes in float32, as a float32 image does:
          faster, with values within about 1e-7 of the largest response value.
        - With a `padding_value`, the padded image counts.
    """
    source_mask = _bool_mask(source_mask)  # a mask of another type: != 0
    # Handle scalar spacing
    if isinstance(spacing_mm, (int, float)):
        spacing_mm = (float(spacing_mm),) * 3
    problem = _log_parameter_problem(sigma_mm, spacing_mm, truncate)
    if problem:
        raise ValueError(problem)

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

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

    # Handle string boundary condition
    boundary = resolve_boundary(boundary)

    problem = _padding_value_problem(boundary, padding_value)
    if problem:
        raise ValueError(problem)
    if padding_value:
        reach = tuple(int(truncate * s + 0.5) for s in sigma_voxels)
        return _constant_padded(
            laplacian_of_gaussian, image, source_mask, padding_value, reach,
            sigma_mm=sigma_mm, spacing_mm=spacing_mm, truncate=truncate, boundary=boundary,
        )  # fmt: skip
    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:
        # A float64 image below _LOG_FFT_MIN voxels runs its 8 passes in float32 (faster;
        # the FFT path of a larger float64 image keeps float64). The response is float32.
        if image.dtype == np.float64 and image.size < _LOG_FFT_MIN:
            image = as_float32(image)
        result = _gaussian_laplace(image, sigma_voxels, mode, truncate, fft_dtype=np.float32)
        if result.dtype == np.float32:
            return result
        return _float32_cut(result, None)

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,
    padding_value: float = 0.0,
) -> 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 > 20,000 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).
        padding_value: The constant of constant value padding (Z3VE), with the ZERO
            (constant) boundary. Default 0.

    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
    """

    source_mask = _bool_mask(source_mask)  # a mask of another type: != 0
    # The passes run in float32. The normalized convolution of a source mask (without
    # rotation invariance) keeps the image type: it divides two filtered sums.
    masked_sums = source_mask is not None and not rotation_invariant
    image = ensure_float32(image) if masked_sums else as_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
    boundary = resolve_boundary(boundary)
    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
    problem = _padding_value_problem(boundary, padding_value)
    if problem:
        raise ValueError(problem)
    # The convolution reads the constant beyond the image; the energy step then pads the
    # response with 0 by the boundary, as with a padding value of 0. A source mask takes a
    # padded copy.
    cval = 0.0
    if padding_value and source_mask is None:
        cval = float(image.dtype.type(padding_value))
    elif padding_value:
        padded = _constant_padded(
            laws_filter, image, source_mask, padding_value, max(g1.size, g2.size, g3.size) // 2,
            kernels=kernels, boundary=boundary, rotation_invariant=rotation_invariant,
            pooling=pooling, use_parallel=use_parallel,
        )  # fmt: skip
        if not compute_energy:
            return padded
        response, valid_part = padded if isinstance(padded, tuple) else (padded, None)
        energy = _laws_energy(response, energy_distance, mode)
        return energy if valid_part is None else (energy, valid_part)

    # 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]

        # One step per rotation: the permuted kernels (the key) and the sign.
        steps: List[Tuple[Tuple[str, str, str], int]] = []
        base_perms: Dict[Tuple[str, str, str], Tuple[int, int, int]] = {}
        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
            steps.append((key, sign))
            base_perms.setdefault(key, perm)
        if pooling != "average":
            # A (key, sign) step that comes again leaves a max or a min unchanged.
            steps = list(dict.fromkeys(steps))

        # The bases run in a pool (at most _BASES_IN_FLIGHT and about 2 GB of them at a
        # time), and each one has its share of the threads for its passes.
        threads = get_num_threads() if use_parallel else 1
        # A small image has no slab threads in its passes, so all its bases run at once
        in_flight = _BASES_IN_FLIGHT if image.size >= _SLAB_MIN_SIZE else threads
        workers = max(1, min(len(base_perms), in_flight, threads, (2 << 30) // (4 * image.size)))

        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,
                threads=threads // workers,
                cval=cval,
            )

        bases = _ordered_map(_base, base_perms.values(), workers)
        result = _pool_rotations(bases, steps, pooling)

        # Finalize average pooling (all bases are done: the slab threads are free)
        if pooling == "average" and result is not None:
            _slab_ufunc(np.true_divide, (result, len(rotations)), result)
    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:
            # The passes of the float32 image give the float32 response
            result = _separable_convolve_3d(image, g1, g2, g3, mode, cval=cval)
            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")

        result = _laws_energy(result, energy_distance, mode)

    if result is None:  # pragma: no cover
        raise RuntimeError("Result should not be None")
    # float32, as the other filters: the rotation-invariant and masked paths pool or divide
    # in float64, and a large result copies in threads
    if result.dtype != np.float32:
        result = _float32_cut(result, None)

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