commit bcaae5ec228b844902b88c520f6ce5a68bda7b12 Author: mohamad Date: Thu Sep 10 01:58:07 2026 +0330 feat: initialize standalone image-filter Rename the distribution to image-filter and imports to image_filter. Keep medical fixtures local and document their setup for fresh clones. Include composite pipelines, consistent examples, and class diagrams. diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..080ead1 --- /dev/null +++ b/.gitignore @@ -0,0 +1,20 @@ +/.pytest_cache/ +/build/ +/dist/ +*.egg-info/ +__pycache__/ +*.py[cod] + +/example-output/ +/data/ +.DS_Store +.venv/ +venv/ +.env +.env.* +.coverage +htmlcov/ +.mypy_cache/ +.ruff_cache/ +*.log +/filter-qa.png diff --git a/README.md b/README.md new file mode 100644 index 0000000..4e24aee --- /dev/null +++ b/README.md @@ -0,0 +1,280 @@ +# Image Filter + +`image-filter` is a standalone, typed Python library for reproducible +filtering of in-memory NumPy images. It has no dependency on the Radiuma or +ViSERA application and performs no file or folder I/O. + +The repository and distribution are named `image-filter`. Python imports use +`image_filter` because hyphens are not valid Python identifiers. Source lives +in `src/image_filter/`. Commands below run from this repository's root. + +See the [class diagrams](docs/architecture.md) for the Composite architecture +and domain model, with editable Mermaid source files. + +## Local medical data + +Real scans are kept locally in ignored `data/samples/` and `data/tests/`; +they are not included in Git or the wheel. This Mac has the existing fixtures +copied into those directories. A fresh clone needs MRI.nii.gz, PET.nii.gz, +CT.nii.gz, SPECT.nii.gz and Ultrasound.nii.gz in `data/samples/`, plus +test4.nii.gz and 0001_MRI_Image_a.nii.gz in `data/tests/` to reproduce the full +real-data suite. Alternatively, examples accept `--input` for your own files. +DICOM demo fixtures come from the installed pydicom package. + +## Install + +```bash +python -m pip install . +``` + +For development: + +```bash +python -m pip install -e '.[test,examples]' +python -m pytest tests +``` + +## Quick start + +```python +import numpy as np +from image_filter import MeanFilter, Pipeline, WaveletFilter + +image = np.asarray(...) # a 2D or 3D numeric image +pipeline = Pipeline([ + MeanFilter(size=3), + WaveletFilter(configuration="LH", wavelet="db2"), + MeanFilter(size=4, origin=-1, mode="nearest"), +]) +result = pipeline.apply(image, {"case_id": "anonymous-17"}) + +filtered = result.image +history = result.metadata["processing_history"] +``` + +Configuration objects and pipelines are immutable. Repeating a filter class is +supported because a pipeline stores ordered filter instances, not a registry +keyed by type. + +## Composite API and runnable code examples + +`ImageFilter` is the abstract Component with `apply(image, metadata)` and +`apply_batch(items)`. The ten filters are Leaves. `Pipeline(ImageFilter)` is +the Composite: its children can be filters, adapters, or nested pipelines. +Each child receives the preceding child's image and metadata. A nested +pipeline records its leaves in execution order, without duplicate pipeline +history entries. Empty pipelines are rejected. Shared instances are immutable. + +```python +from image_filter import GaussianFilter, ImageFilter, Pipeline, SobelFilter + +def process(image, operation: ImageFilter): + return operation.apply(image) + +preprocessing = Pipeline([GaussianFilter(sigma=1.0)]) +edges = Pipeline([preprocessing, SobelFilter()]) +# process(image, GaussianFilter()) and process(image, edges) use the same API. +``` + +All five files have the same structure: `build_operation() -> ImageFilter`, +then `run_example(build_operation(), name=...)`. Only the operation changes. +The shared runner loads the image, applies the component, and saves the result. +It defaults to local MRI data; every example accepts the same input options. + +| File | Demonstrates | +| --- | --- | +| `examples/01_single_filter.py` | One-filter pipeline | +| `examples/02_sequential_filters.py` | Mean → wavelet → mean with different settings | +| `examples/03_nested_pipeline.py` | Reusable preprocessing composite followed by Sobel | +| `examples/04_component_api.py` | A direct leaf passed to the same runner | +| `examples/05_slicewise_pipeline.py` | Gaussian → Canny; use `--mode slice-wise` for a volume | + +From the library root, with example dependencies installed: + +```bash +cd /Users/mohamad/vollabor/image-filter +for example in examples/0[1-5]_*.py; do + PYTHONDONTWRITEBYTECODE=1 PYTHONPATH=src python3 "$example" || break +done +``` + +Outputs are NumPy arrays and JSON processing metadata in +`example-output/pipelines/`. Matching filenames are overwritten; +use different `--output-dir` values to retain multiple runs. + +All five scripts accept these options: + +| Option | Behavior | +| --- | --- | +| `--case mri` | Local case: mri, pet, ct, spect, ultrasound, mr-dicom, ct-series | +| `--input PATH` | NIfTI, NRRD, single DICOM, or DICOM series directory | +| `--series-uid UID` | Select a series when the directory contains several | +| `--mode plane` | Default: representative 2D plane, 1–99% window normalized to [0, 1] | +| `--mode volume` | Original intensities and dimensions; every step must support the input | +| `--mode slice-wise` | Original intensities; automatically wrap a 3D input in `SliceWise` | +| `--axis 2` | Plane selection or slice-wise axis (0, 1, or 2) | +| `--modality CT` | Label for a user-supplied image | +| `--output-dir PATH` | Destination for .npy and .json files | + +```bash +PYTHONPATH=src python3 examples/03_nested_pipeline.py --case ct +PYTHONPATH=src python3 examples/01_single_filter.py --case mri --mode volume +PYTHONPATH=src python3 examples/05_slicewise_pipeline.py --case mri --mode slice-wise +PYTHONPATH=src python3 examples/03_nested_pipeline.py --input /path/to/scan.nii.gz +``` + +`slice-wise` is an explicit user choice; within that mode a 2D input runs +directly and a 3D input gets the adapter automatically. The runner does not +change wavelet configuration or thresholds to accommodate incompatible data. +For example, `LH` in example 02 requires 2D execution: use `plane` or +`slice-wise`, or explicitly configure a 3D wavelet for `volume` mode. + +The existing `visual_compare.py` still produces all ten independent medical +comparisons; it now executes each operation through a one-step composite. +These comparison panels are independent, whereas examples 02, 03 and 05 are +sequential pipelines. + +## Dimensionality and channels + +Dimension is inferred for every image independently. A singleton spatial axis +does not change a 2D image into a volume and its original shape is retained: +`(1, H, W)`, `(H, 1, W)`, and `(H, W, 1)` are all 2D. + +Always declare a channel axis: + +```python +result = MeanFilter(3).apply(rgb, {"channel_axis": -1}) +``` + +A three-axis `(H, W, 3)` or `(H, W, 4)` array without metadata is deliberately +rejected because it may be either color data or a short volume. Set +`channel_axis=2` for color, or `spatial_dims=3` for a volume. Four-axis data +must declare `channel_axis`; time-series inference is intentionally outside the +core contract. + +## Explicit slice-wise filtering + +Gabor and Canny are intrinsically 2D. It raises `UnsupportedDimensionError` for a volume. +Slice processing is available only through an explicit adapter with a chosen +axis: + +```python +from image_filter import GaborFilter, SliceWise + +filter_3d = SliceWise(GaborFilter(frequency=0.2), axis=0) +result = filter_3d.apply(volume) +``` + +See [docs/semantics.md](docs/semantics.md) for numerical contracts, boundary +behavior, capability declarations, and intentionally unsupported options. + +## Collection execution + +```python +from image_filter import ImageItem + +results = Pipeline([MeanFilter(3)]).apply_batch([ + ImageItem(two_dimensional_image, {"id": "2d"}), + ImageItem(volume, {"id": "3d", "spatial_dims": 3}), +]) +``` + +No caller metadata dictionary is mutated. Each output receives a deep-copied +mapping with deterministic, ordered processing provenance. + + +## Ten public filters + +| Filter | Spatial dimensions | Output | +| --- | --- | --- | +| MeanFilter | 2D, 3D | float32/float64 | +| GaussianFilter | 2D, 3D | float32/float64 | +| MedianFilter | 2D, 3D | float32/float64 | +| SobelFilter | 2D, 3D | float32/float64 magnitude by default | +| LaplacianFilter | 2D, 3D | float32/float64 signed response | +| LoGFilter | 2D, 3D | float32/float64 | +| LawsFilter | 2D, 3D | float32/float64 | +| GaborFilter | 2D | float32/float64 | +| CannyFilter | 2D | bool | +| WaveletFilter | 2D, 3D | float32/float64 | + +```python +from skimage.data import camera +from image_filter import GaussianFilter, MedianFilter, CannyFilter, Pipeline + +edges = Pipeline([ + GaussianFilter(sigma=0.8), + MedianFilter(size=3), + CannyFilter(low_threshold=0.1, high_threshold=0.3, use_quantiles=True), +]).apply(camera(), {"source": "public-domain camera"}) +assert edges.image.dtype == bool +``` + +Canny follows scikit-image dtype scaling; float images are not automatically +normalized. See the [threshold and boundary contracts](docs/semantics.md). + +## Reproduce validation on macOS + +From the library root, with Python 3.10+ and the dependencies installed: + +```bash +cd /Users/mohamad/vollabor/image-filter +python3 -m pip install -e '.[test,examples]' build +PYTHONDONTWRITEBYTECODE=1 PYTHONPATH=src python3 -m pytest tests -p no:cacheprovider +PYTHONDONTWRITEBYTECODE=1 PYTHONPATH=src MPLBACKEND=Agg python3 examples/visual_compare.py --all --output-dir example-output --contact-sheet +python3 -m build --wheel --no-isolation --outdir dist . +open example-output/source-overview.png +``` + +Real-data tests read local NIfTI assets and medical fixtures. They do not +import application runtime code. + +## Generate the presentation examples + +Generate seven consistently styled medical comparison panels and a source +overview with one command, from the library root (development dependencies +installed): + +```bash +cd /Users/mohamad/vollabor/image-filter +PYTHONDONTWRITEBYTECODE=1 PYTHONPATH=src python3 examples/visual_compare.py --all --output-dir example-output --contact-sheet +open example-output/source-overview.png +``` + +The bundled command uses local local NIfTI volumes and installed pydicom +fixtures: + +| Source | Input | Panel filename | +| --- | --- | --- | +| mri | local MRI NIfTI volume | mri-filters.png | +| pet | local PET NIfTI volume | pet-filters.png | +| ct | local CT NIfTI volume | ct-filters.png | +| spect | local SPECT NIfTI volume | spect-filters.png | +| ultrasound | local ultrasound NIfTI volume | ultrasound-filters.png | +| mr-dicom | installed pydicom single MR DICOM | mr-dicom-filters.png | +| ct-series | installed pydicom five-slice CT series | ct-series-filters.png | + +All sources are loaded locally; no network download is attempted. Volumes are +canonicalized to RAS axes where the reader supports it. Each panel selects a +deterministic high-contrast non-border slice and labels modality, format, +source, axis/index, orientation, and file count. All ten filters operate +independently on that selected 2D slice. Intensities are windowed to the 1st– +99th percentile and normalized to [0, 1] for the demo only; filter outputs are +displayed with independent contrast. + +All panels are 2400×1800 PNGs. The overview, `source-overview.png`, is +1800×1350 for the seven-case set. Outputs remain in the ignored +`example-output/` directory, also the default directory +for `--all` or `--images` when `--output-dir` is omitted. Rerunning overwrites +matching filenames. The existing root `filter-qa.png` is preserved by this +multi-image command. + +Select a subset with `--images mri pet ct`, or supply user files with +`--input /path/case.nii.gz /path/case.nrrd /path/image.dcm` and +`--dicom-series /path/series-directory`. NRRD is supported through SimpleITK, +but no bundled NRRD fixture exists in this checkout, so it is not claimed in +the default demonstration. A DICOM directory is read as one ordered volume +using SimpleITK's series reader. `--slice-axis` and `--slice-index` override +the deterministic selection. `--contact-sheet` is optional. `--output` and +`--output-dir` are mutually exclusive. `cases.json` records the exact source, +format, dimensions, slice choice, and window for each generated panel. diff --git a/docs/architecture.md b/docs/architecture.md new file mode 100644 index 0000000..179faee --- /dev/null +++ b/docs/architecture.md @@ -0,0 +1,123 @@ +# Image Filter class diagrams + +These UML-style Mermaid diagrams describe the implemented classes in +`src/image_filter`. Open this Markdown in a Mermaid-capable viewer, or use +the adjacent `.mmd` source files with Mermaid CLI. + +## Composite and numerical filters + +```mermaid +classDiagram +direction TB +class ImageFilter { + <> + +Capabilities capabilities + +apply(image, metadata) FilterResult + +apply_batch(items) list~FilterResult~ +} +class NumericFilter { + <> + +apply(image, metadata) FilterResult + #_apply_with_layout(image, metadata, layout) FilterResult + #_compute(image, layout) ndarray + #_history_entry(image, output, layout) dict +} +class Pipeline { + +tuple~ImageFilter~ steps + +Capabilities capabilities + +apply(image, metadata) FilterResult +} +class SliceWise { + +ImageFilter filter + +int axis + +apply(image, metadata) FilterResult + #_compute(image, layout) ndarray +} +ImageFilter <|-- NumericFilter +ImageFilter <|-- Pipeline +NumericFilter <|-- SliceWise +Pipeline "0..*" o-- "1..*" ImageFilter : ordered children +SliceWise "0..*" o-- "1" ImageFilter : adapts +NumericFilter <|-- MeanFilter +NumericFilter <|-- GaussianFilter +NumericFilter <|-- MedianFilter +NumericFilter <|-- SobelFilter +NumericFilter <|-- LaplacianFilter +NumericFilter <|-- LoGFilter +NumericFilter <|-- LawsFilter +NumericFilter <|-- GaborFilter +NumericFilter <|-- CannyFilter +NumericFilter <|-- WaveletFilter +``` + +ImageFilter is the public Component, Pipeline the Composite, and the ten +concrete numerical filters the Leaves. NumericFilter shares validation and +history behavior; it is an internal abstract implementation, not a separate +public filtering operation. SliceWise is an adapter that can wrap a Leaf or +a Pipeline. + +The hollow diamond denotes aggregation: immutable children can be shared +between pipelines and remain valid independently of any parent. A Pipeline +contains at least one child and executes children in order. Inheritance arrows +point toward the base class. Public methods use +; internal methods use #. + +## Input, output, and capability model + +```mermaid +classDiagram +direction LR +class ImageItem { + +ndarray image + +Mapping metadata +} +class FilterResult { + +ndarray image + +Mapping metadata +} +class SpatialLayout { + +int dimension + +tuple spatial_axes + +tuple singleton_axes + +int channel_axis +} +class Capabilities { + +frozenset spatial_dimensions + +str accepted_dtypes + +str input_range + +str output_dtype + +bool preserves_shape + +bool supports_channels +} +class ImageFilter { + <> + +apply(image, metadata) FilterResult + +apply_batch(items) list~FilterResult~ +} +class NumericFilter { + <> +} +ImageFilter <|-- NumericFilter +ImageFilter ..> ImageItem : consumes in batches +ImageFilter ..> FilterResult : returns +ImageFilter ..> Capabilities : declares +NumericFilter ..> SpatialLayout : validates and computes +``` + +ImageItem groups an array with metadata for batch execution. FilterResult +contains the resulting array and copied metadata. SpatialLayout describes the +effective spatial axes, ignoring singleton axes and excluding an explicitly +declared channel axis. Capabilities describes supported input dimensions and +other contracts. Pipeline derives its capabilities from its children. + +## Source and presentation order + +1. `src/image_filter/base.py`: Component and common numerical behavior. +2. `src/image_filter/filters.py`: the ten Leaves. +3. `src/image_filter/pipeline.py`: Composite and slice adapter. +4. `src/image_filter/model.py`: domain objects and dimension inference. +5. `examples/01_single_filter.py` through `05_slicewise_pipeline.py`: + a shared runner with different component configurations. + +Editable sources: [components.mmd](components.mmd) and +[domain-model.mmd](domain-model.mmd). Keep their diagrams identical to the +embedded blocks here when changing class relationships. diff --git a/docs/components.mmd b/docs/components.mmd new file mode 100644 index 0000000..d49058d --- /dev/null +++ b/docs/components.mmd @@ -0,0 +1,41 @@ +classDiagram +direction TB +class ImageFilter { + <> + +Capabilities capabilities + +apply(image, metadata) FilterResult + +apply_batch(items) list~FilterResult~ +} +class NumericFilter { + <> + +apply(image, metadata) FilterResult + #_apply_with_layout(image, metadata, layout) FilterResult + #_compute(image, layout) ndarray + #_history_entry(image, output, layout) dict +} +class Pipeline { + +tuple~ImageFilter~ steps + +Capabilities capabilities + +apply(image, metadata) FilterResult +} +class SliceWise { + +ImageFilter filter + +int axis + +apply(image, metadata) FilterResult + #_compute(image, layout) ndarray +} +ImageFilter <|-- NumericFilter +ImageFilter <|-- Pipeline +NumericFilter <|-- SliceWise +Pipeline "0..*" o-- "1..*" ImageFilter : ordered children +SliceWise "0..*" o-- "1" ImageFilter : adapts +NumericFilter <|-- MeanFilter +NumericFilter <|-- GaussianFilter +NumericFilter <|-- MedianFilter +NumericFilter <|-- SobelFilter +NumericFilter <|-- LaplacianFilter +NumericFilter <|-- LoGFilter +NumericFilter <|-- LawsFilter +NumericFilter <|-- GaborFilter +NumericFilter <|-- CannyFilter +NumericFilter <|-- WaveletFilter diff --git a/docs/domain-model.mmd b/docs/domain-model.mmd new file mode 100644 index 0000000..d2cc2a9 --- /dev/null +++ b/docs/domain-model.mmd @@ -0,0 +1,37 @@ +classDiagram +direction LR +class ImageItem { + +ndarray image + +Mapping metadata +} +class FilterResult { + +ndarray image + +Mapping metadata +} +class SpatialLayout { + +int dimension + +tuple spatial_axes + +tuple singleton_axes + +int channel_axis +} +class Capabilities { + +frozenset spatial_dimensions + +str accepted_dtypes + +str input_range + +str output_dtype + +bool preserves_shape + +bool supports_channels +} +class ImageFilter { + <> + +apply(image, metadata) FilterResult + +apply_batch(items) list~FilterResult~ +} +class NumericFilter { + <> +} +ImageFilter <|-- NumericFilter +ImageFilter ..> ImageItem : consumes in batches +ImageFilter ..> FilterResult : returns +ImageFilter ..> Capabilities : declares +NumericFilter ..> SpatialLayout : validates and computes diff --git a/docs/semantics.md b/docs/semantics.md new file mode 100644 index 0000000..5ca5800 --- /dev/null +++ b/docs/semantics.md @@ -0,0 +1,174 @@ +# Numerical semantics and scope + +This package is a greenfield implementation. Its numerical contracts are the +documented behavior of SciPy, scikit-image, PyWavelets, and NumPy—not any +application-specific filter implementation. + +## Shared contract + +- Inputs are finite, real numeric NumPy arrays with two or three inferred + non-singleton spatial axes. +- Filters preserve the complete input shape, including singleton and channel + axes. +- Numeric output is explicitly `float32` or `float64` (default `float64`). + Canny always returns a boolean mask. +- SciPy-backed filters accept `reflect`, `constant`, `nearest`, `mirror`, and + `wrap`. `cval` is used only by `constant` mode. +- Each filter exposes a `Capabilities` object. Unsupported dimensions and + recognized-but-unimplemented options use domain-specific exceptions. + +## MeanFilter + +`MeanFilter` calls `scipy.ndimage.uniform_filter` exactly, with a scalar kernel +size applied only to spatial axes. Channel and retained singleton axes use size +one. + +SciPy's `origin=0` alignment is the default. Odd kernels are centered on a +sample. Even kernels necessarily have one extra sample on one side of the +origin; the exact placement is SciPy's documented convention. Set `origin` +explicitly to shift the footprint. Valid origins are checked at construction. +This policy makes even-kernel behavior deliberate and testable. + +## LoGFilter + +`LoGFilter` sums second-order `scipy.ndimage.gaussian_filter` derivatives over +the inferred spatial axes. This is numerically equivalent to +`scipy.ndimage.gaussian_laplace` when every array axis is spatial, while also +preventing smoothing or differentiation across channels. Sigma is expressed in +pixels and may be scalar or one value per effective spatial axis. + +## LawsFilter + +Laws responses use normalized standard L3/E3/S3 and L5/E5/S5/W5/R5 vectors +with `scipy.ndimage.convolve1d`. A configuration contains one two-character +kernel name per spatial axis, such as `L5E5` or `L5E5S5`. Optional texture +energy is the local uniform mean of the absolute response with window width +`2 * energy_radius + 1`. + +Rotation-invariant pooling is intentionally unsupported. Construct explicit +orientation filters and combine their outputs according to the domain's chosen +statistic instead of relying on an undocumented convention. + +## GaborFilter + +Gabor responses call `skimage.filters.gabor` for each orientation. Frequency is +in cycles per pixel, angles are radians, and `sigma_y = sigma / gamma`. +Available complex response projections are magnitude, real, and imaginary; +multiple orientations use max, min, mean, or sum pooling. + +The filter supports 2D data, including 2D multichannel arrays with an explicit +channel axis. A 3D volume requires `SliceWise`, whose axis is part of the public +configuration and provenance. + +## WaveletFilter + +Wavelet coefficients use PyWavelets' shift-invariant `pywt.swtn`. Configuration +letters map directly to PyWavelets keys: `L` to approximation (`a`) and `H` to +detail (`d`) in the order of inferred spatial axes. The selected array is the +requested coefficient at the final decomposition level. + +SWT requires transformed lengths divisible by `2**level`. Each spatial axis is +padded only at its high edge to the next valid length using `numpy.pad`, then +the coefficient is cropped to the original shape. `pad_mode` is therefore a +NumPy padding mode rather than a PyWavelets signal-extension mode. The optional +`normalize=True` is passed directly to `pywt.swtn`. + +Rotation-invariant coefficient pooling and Riesz steering are not implemented; +requesting rotation invariance raises `UnsupportedOptionError`. Decimated DWT +reconstruction is a distinct operation and is outside this first API. + + +## GaussianFilter + +2D/3D Gaussian smoothing calls [scipy.ndimage.gaussian_filter](https://docs.scipy.org/doc/scipy/reference/generated/scipy.ndimage.gaussian_filter.html). +`sigma` is a finite positive scalar or immutable nonempty tuple, one value per +inferred spatial axis in array order. Tuple length is validated on application. +`truncate` is finite and positive (default 4 standard deviations). Input and +intermediate values use the selected floating output dtype; no intensity +normalization occurs. Channels and singleton axes are indexed out before the +backend call, preserving SciPy 1.10 compatibility. + +## MedianFilter + +2D/3D rank filtering calls [scipy.ndimage.median_filter](https://docs.scipy.org/doc/scipy/reference/generated/scipy.ndimage.median_filter.html). +`size` is a positive integer applied along each spatial axis. `origin` must be +an integer in `[-size//2, (size-1)//2]`, interpreting the lower bound as +`-(size//2)`. Even footprints select sorted element `n//2`, rather than averaging +the middle pair. Input is converted to the selected floating dtype before +filtering, including fractional `cval` padding for integer images and support +for float16 input. Large 3D windows can require substantial backend memory. + +## SobelFilter + +2D/3D gradients call [scipy.ndimage.sobel](https://docs.scipy.org/doc/scipy/reference/generated/scipy.ndimage.sobel.html) +on each spatial block. The default `axis=None` returns Euclidean magnitude +`sqrt(sum(sobel(image, axis=i)**2))` over spatial axes, accumulated with `hypot` +to avoid intermediate squared overflow. There is no division by dimensionality +or kernel gain. Optional integer `axis` selects a signed component, indexing the +inferred spatial axes (negative indices supported). Channels are neither +smoothed nor differentiated. Computation uses the selected floating dtype. + +## LaplacianFilter + +2D/3D [scipy.ndimage.laplace](https://docs.scipy.org/doc/scipy/reference/generated/scipy.ndimage.laplace.html) +sums discrete second differences over spatial axes. Input is cast to the +selected floating dtype, retaining signed responses. This operation adds no +Gaussian smoothing; use LoGFilter when that is desired. + +## CannyFilter + +[skimage.feature.canny](https://scikit-image.org/docs/stable/api/skimage.feature.html#skimage.feature.canny) +produces a boolean edge mask from each 2D plane/channel. Singleton axes are +restored afterward. A 3D image is rejected; use `SliceWise(CannyFilter(), axis=0)` +with an explicitly chosen array axis. There is no implicit projection or slicing. + +Canny rejects int64/uint64 with `InvalidImageError`, matching the backend +restriction declared in its capabilities. Callers may explicitly convert to float +and choose thresholds in the resulting intensity units. + +The original supported input dtype reaches Canny. Its image conversion scales unsigned +integers by their dtype maximum and signed integers to the signed unit range; +floating input keeps its values, with no min/max normalization or clipping. +Absolute thresholds are expressed in original dtype intensity units; the +backend divides integer thresholds by the dtype maximum. For floating input, +thresholds are already in the unchanged intensity units. Missing low/high +thresholds become 0.1/0.2 in the backend's normalized scale. `use_quantiles=True` +instead interprets thresholds as gradient-magnitude quantiles in [0, 1], with +defaults 0.1/0.2. Thresholds must be finite, nonnegative, and ordered; ordering +against missing defaults for absolute thresholds is checked by the backend at +application time and converted to `InvalidConfigurationError`. + +`sigma` must be finite and positive. `cval` is in the converted image's units. +`mode` defaults to `constant`; it governs Gaussian padding, and Canny excludes +the outermost edge pixels during mask erosion. There is no public mask option. + +## Boundary and precision details for the new filters + +All five accept `reflect` (half-sample reflection), `mirror` (whole-sample +reflection), `nearest` (edge replication), `wrap` (periodic), and `constant` +(finite `cval`, default 0). SciPy aliases are intentionally excluded. Numeric +filters default to `reflect`. Boolean configuration values are not numbers. +Like their backends, floating filters can lose precision or overflow for extreme +finite inputs; finite outputs are verified on normal image ranges, not promised +for arbitrary floating-point magnitudes. Provenance records parameters, effective +spatial axes, shapes and dtypes; SliceWise records the nested filter and axis. + +### Composite contract + +`ImageFilter` defines `apply` and shared `apply_batch` behavior. Numerical +leaves reuse the internal `NumericFilter` implementation for validation and +history. `Pipeline` is an `ImageFilter` containing a nonempty immutable tuple +of components. Nested pipelines are supported at arbitrary depth and execute +depth-first in child order. History contains each executed leaf once; composite +boundaries do not add history entries. Errors propagate from the failing child. + +Pipeline dimensions are the intersection of child dimensions. This declaration +does not guarantee that intermediate ranges or dtypes satisfy later children; +each child validates its actual input. The last child determines output dtype. +Batch processing uses the same component API for both leaves and pipelines. + +`SliceWise` accepts any component with 2D capability, including a pipeline, +and calls its public API on each plane with remapped channel metadata. The +actual first plane result determines output dtype (including boolean Canny +results). Its outer history records one adapter entry containing the nested +configuration, rather than repeating history for every slice. diff --git a/examples/01_single_filter.py b/examples/01_single_filter.py new file mode 100644 index 0000000..623871b --- /dev/null +++ b/examples/01_single_filter.py @@ -0,0 +1,11 @@ +"""One-filter pipeline.""" +from image_filter import GaussianFilter, ImageFilter, Pipeline +from example_support import run_example + + +def build_operation() -> ImageFilter: + return Pipeline([GaussianFilter(sigma=1.0)]) + + +if __name__ == "__main__": + run_example(build_operation(), name="01-single-filter") diff --git a/examples/02_sequential_filters.py b/examples/02_sequential_filters.py new file mode 100644 index 0000000..295d2f6 --- /dev/null +++ b/examples/02_sequential_filters.py @@ -0,0 +1,15 @@ +"""Mean → wavelet → mean: each step consumes the previous output.""" +from image_filter import ImageFilter, MeanFilter, Pipeline, WaveletFilter +from example_support import run_example + + +def build_operation() -> ImageFilter: + return Pipeline([ + MeanFilter(size=3), + WaveletFilter(configuration="LH", wavelet="db2"), + MeanFilter(size=4, origin=-1), + ]) + + +if __name__ == "__main__": + run_example(build_operation(), name="02-sequential-filters") diff --git a/examples/03_nested_pipeline.py b/examples/03_nested_pipeline.py new file mode 100644 index 0000000..53fd61d --- /dev/null +++ b/examples/03_nested_pipeline.py @@ -0,0 +1,12 @@ +"""Reuse preprocessing as a child of another pipeline.""" +from image_filter import GaussianFilter, ImageFilter, MedianFilter, Pipeline, SobelFilter +from example_support import run_example + + +def build_operation() -> ImageFilter: + preprocessing = Pipeline([MedianFilter(size=3), GaussianFilter(sigma=0.8)]) + return Pipeline([preprocessing, SobelFilter()]) + + +if __name__ == "__main__": + run_example(build_operation(), name="03-nested-pipeline") diff --git a/examples/04_component_api.py b/examples/04_component_api.py new file mode 100644 index 0000000..af75d7e --- /dev/null +++ b/examples/04_component_api.py @@ -0,0 +1,11 @@ +"""A direct leaf uses exactly the same runner as a composite.""" +from image_filter import GaussianFilter, ImageFilter +from example_support import run_example + + +def build_operation() -> ImageFilter: + return GaussianFilter(sigma=1.0) + + +if __name__ == "__main__": + run_example(build_operation(), name="04-component-api") diff --git a/examples/05_slicewise_pipeline.py b/examples/05_slicewise_pipeline.py new file mode 100644 index 0000000..e514f42 --- /dev/null +++ b/examples/05_slicewise_pipeline.py @@ -0,0 +1,14 @@ +"""Gaussian → Canny; use --mode slice-wise for every volume slice.""" +from image_filter import CannyFilter, GaussianFilter, ImageFilter, Pipeline +from example_support import run_example + + +def build_operation() -> ImageFilter: + return Pipeline([ + GaussianFilter(sigma=1.0), + CannyFilter(low_threshold=0.1, high_threshold=0.3, use_quantiles=True), + ]) + + +if __name__ == "__main__": + run_example(build_operation(), name="05-slicewise-pipeline") diff --git a/examples/example_support.py b/examples/example_support.py new file mode 100644 index 0000000..87049b6 --- /dev/null +++ b/examples/example_support.py @@ -0,0 +1,51 @@ +"""One execution path for every example; I/O stays outside the library.""" +from __future__ import annotations + +import argparse +import json +from pathlib import Path + +import numpy as np + +from medical_io import CASE_NAMES, bundled_cases, load_medical, representative_plane +from image_filter import ImageFilter, SliceWise + + +def run_example(operation: ImageFilter, name: str, argv: list[str] | None = None) -> None: + """Load → adapt explicitly → apply → save, for leaves and pipelines alike.""" + parser = argparse.ArgumentParser(description=__doc__) + source = parser.add_mutually_exclusive_group() + source.add_argument("--case", choices=CASE_NAMES, default="mri") + source.add_argument("--input", type=Path, help="NIfTI, NRRD, DICOM file or series directory") + parser.add_argument("--series-uid", help="Select among multiple DICOM series") + parser.add_argument("--modality", default="Unknown") + parser.add_argument("--mode", choices=("plane", "volume", "slice-wise"), default="plane") + parser.add_argument("--axis", type=int, choices=(0, 1, 2), default=2) + parser.add_argument("--output-dir", type=Path, + default=Path(__file__).resolve().parents[1] / "example-output" / "pipelines") + args = parser.parse_args(argv) + if args.input is not None: + path, modality = args.input, args.modality + else: + catalog = bundled_cases() + if args.case not in catalog: + parser.error(f"Local case unavailable: {args.case}") + path, modality = catalog[args.case] + case = load_medical(path, modality, args.series_uid) + if args.mode == "plane": + image, metadata = representative_plane(case, axis=args.axis) + metadata["processing"] = "sequential processing on selected 2D plane" + else: + image = case.pixels + metadata = {"modality": case.modality, "format": case.format, + "spatial_dims": image.ndim, "processing": args.mode} + if args.mode == "slice-wise" and image.ndim == 3: + operation = SliceWise(operation, axis=args.axis) + # Native 2D inputs need no adapter. Volume mode never silently slices. + result = operation.apply(image, metadata) + args.output_dir.mkdir(parents=True, exist_ok=True) + output = args.output_dir / f"{name}.npy" + np.save(output, result.image) + output.with_suffix(".json").write_text(json.dumps(dict(result.metadata), indent=2) + "\n") + history = [entry["filter"] for entry in result.metadata["processing_history"]] + print(f"{' -> '.join(history)} | shape={result.image.shape} | {output.resolve()}") diff --git a/examples/medical_io.py b/examples/medical_io.py new file mode 100644 index 0000000..dcddbfb --- /dev/null +++ b/examples/medical_io.py @@ -0,0 +1,130 @@ +"""Example-only local medical I/O; the filter library remains NumPy-only.""" +from __future__ import annotations + +from dataclasses import dataclass +from pathlib import Path + +import nibabel as nib +import numpy as np + +SAMPLES = Path(__file__).resolve().parents[1] / 'data/samples' +CASE_NAMES = ('mri', 'pet', 'ct', 'spect', 'ultrasound', 'mr-dicom', 'ct-series') + + +@dataclass +class MedicalImage: + pixels: np.ndarray + modality: str + format: str + source: Path + orientation: str + files: int = 1 + + +def bundled_cases() -> dict[str, tuple[Path, str]]: + cases = {name.lower(): (SAMPLES / f'{name}.nii.gz', name) for name in ['MRI', 'PET', 'CT', 'SPECT', 'Ultrasound']} + try: + import pydicom + root = Path(pydicom.__file__).parent / 'data/test_files' + cases['mr-dicom'] = (root / 'MR_small.dcm', 'MR') + cases['ct-series'] = (root / 'dicomdirtests/98892001/CT5N', 'CT') + except ImportError: + pass + return cases + + +def _validate(pixels: np.ndarray) -> np.ndarray: + # Only singleton extras are removed; time/channel dimensions are rejected. + pixels = np.squeeze(pixels) + if pixels.ndim not in (2, 3) or pixels.dtype.kind not in 'buif': + raise ValueError('Only scalar 2D/3D medical images are supported; select time/channel externally') + if not np.isfinite(pixels).all(): + raise ValueError('Medical image contains nonfinite values') + return pixels + + +def load_medical(path: Path, modality: str = 'Unknown', series_uid: str | None = None) -> MedicalImage: + path = path.resolve() + if not path.exists(): + raise FileNotFoundError(f'Local medical fixture unavailable: {path}; no download attempted') + if path.name.lower().endswith(('.nii', '.nii.gz')): + source = nib.as_closest_canonical(nib.load(path)) + pixels = _validate(source.get_fdata()) + return MedicalImage(pixels, modality, 'NIfTI', path, 'RAS voxel axes') + if path.is_dir() or path.suffix.lower() in ('.nrrd', '.nhdr'): + import SimpleITK as sitk + count = 1 + if path.is_dir(): + ids = sitk.ImageSeriesReader.GetGDCMSeriesIDs(str(path)) or () + if not ids: + raise ValueError(f'No DICOM series in {path.name}') + if series_uid is None and len(ids) != 1: + raise ValueError('Multiple DICOM series found; select one with --series-uid') + uid = series_uid or ids[0] + if uid not in ids: + raise ValueError('--series-uid not found in directory') + files = sitk.ImageSeriesReader.GetGDCMSeriesFileNames(str(path), uid) + if len(files) < 2: + raise ValueError('DICOM series directory must contain at least two slices') + import pydicom + modality = str(pydicom.dcmread(files[0], stop_before_pixels=True).get('Modality', modality)) + reader = sitk.ImageSeriesReader() + reader.SetFileNames(files) + source = reader.Execute() + count = len(files) + format_name = 'DICOM series' + else: + source = sitk.ReadImage(str(path)) + format_name = 'NRRD' + if source.GetNumberOfComponentsPerPixel() != 1: + raise ValueError('Vector/color images are not supported in this medical demo') + if source.GetDimension() == 3: + source = sitk.DICOMOrient(source, 'RAS') + pixels = sitk.GetArrayFromImage(source) + # SimpleITK arrays are z,y,x; expose x,y,z like canonical nibabel. + pixels = pixels.transpose(tuple(reversed(range(pixels.ndim)))) + return MedicalImage(_validate(pixels), modality, format_name, path, 'RAS voxel axes' if source.GetDimension() == 3 else 'native axes', count) + if path.suffix.lower() == '.dcm': + import pydicom + from pydicom.pixels import apply_modality_lut + source = pydicom.dcmread(path) + if int(source.get('SamplesPerPixel', 1)) != 1 or int(source.get('NumberOfFrames', 1)) != 1: + raise ValueError('Single DICOM accepts one monochrome frame; use a series directory for volumes') + pixels = apply_modality_lut(source.pixel_array, source).astype(np.float64) + if source.get('PhotometricInterpretation') == 'MONOCHROME1': + pixels = pixels.max() + pixels.min() - pixels + return MedicalImage(_validate(pixels), str(source.get('Modality', modality)), 'DICOM', path, 'native row/column') + raise ValueError('Expected .nii/.nii.gz, .nrrd/.nhdr, .dcm, or a DICOM series directory') + + +def representative_plane(case: MedicalImage, axis: int = 2, index: int | None = None) -> tuple[np.ndarray, dict]: + pixels = case.pixels + selected = None + if pixels.ndim == 3: + if index is not None and not 0 <= index < pixels.shape[axis]: + raise ValueError(f'Slice index outside axis length {pixels.shape[axis]}') + if index is None: + # Ignore border slices; maximize contrast among the central 80%. + lo = int(pixels.shape[axis] * .1) + hi = max(lo + 1, int(np.ceil(pixels.shape[axis] * .9))) + scores = np.std(pixels, axis=tuple(i for i in range(3) if i != axis)) + selected = lo + int(np.argmax(scores[lo:hi])) + else: + selected = index + plane = np.rot90(np.take(pixels, selected, axis=axis)) + else: + if index is not None: + raise ValueError('--slice-index requires a volume') + plane = pixels + low, high = np.percentile(plane, (1, 99)) + if high <= low: + low, high = float(plane.min()), float(plane.max()) + if high <= low: + raise ValueError('Selected slice has no intensity variation') + image = np.clip((plane - low) / (high - low), 0, 1).astype(np.float64) + info = {'modality': case.modality, 'format': case.format, 'source': case.source.name, + 'volume_shape': list(pixels.shape), 'orientation': case.orientation, + 'slice_axis': axis if selected is not None else None, 'slice_index': selected, + 'window': [float(low), float(high)], 'files': case.files, + 'input_shape': list(image.shape), 'processing': 'all filters independently on selected 2D plane'} + return image, info diff --git a/examples/visual_compare.py b/examples/visual_compare.py new file mode 100644 index 0000000..7d4b49c --- /dev/null +++ b/examples/visual_compare.py @@ -0,0 +1,155 @@ +"""Offline medical-image filter comparisons from local volumes and DICOM.""" +from __future__ import annotations + +import argparse +from pathlib import Path + +import matplotlib +matplotlib.use("Agg") +import matplotlib.pyplot as plt +import numpy as np +import json +import re +import textwrap +import sys + +sys.path.insert(0, str(Path(__file__).resolve().parent)) +from medical_io import CASE_NAMES, bundled_cases, load_medical, representative_plane + +from image_filter import ( + CannyFilter, GaussianFilter, MedianFilter, SobelFilter, LaplacianFilter, + GaborFilter, LawsFilter, LoGFilter, MeanFilter, WaveletFilter, Pipeline, +) + +DEFAULT_DIR = Path(__file__).resolve().parents[1] / "example-output" +STYLE = {"font.family": "DejaVu Sans", "font.size": 11, + "text.color": "#18283b", "axes.titlecolor": "#18283b", + "figure.facecolor": "#f4f6f9", "savefig.facecolor": "#f4f6f9"} + + +def render_panel(name: str, image: np.ndarray, info: dict, output: Path) -> None: + operations = [ + ("Mean 5×5", MeanFilter(5)), + ("Gaussian σ=1.5", GaussianFilter(1.5)), + ("Median 5×5", MedianFilter(5)), + ("Sobel magnitude", SobelFilter()), + ("Laplacian", LaplacianFilter()), + ("Canny σ=1.2 (quantiles .1/.3)", CannyFilter(1.2, .1, .3, True)), + ("Laws L5E5 energy r=3", LawsFilter("L5E5", energy_radius=3)), + ("LoG σ=1.5", LoGFilter(1.5)), + ("Gabor magnitude", GaborFilter(0.15)), + ("SWT LH (db2)", WaveletFilter("LH", "db2")), + ] + # One-step composites preserve the existing independent comparisons. + panels = [("Input", image)] + [ + (title, Pipeline([operation]).apply(image).image) + for title, operation in operations + ] + with plt.rc_context(STYLE): + figure, axes = plt.subplots(3, 4, figsize=(16, 12)) + figure.subplots_adjust(left=.025, right=.975, bottom=.035, top=.88, + hspace=.23, wspace=.10) + figure.text(.025, .965, "IMAGE FILTER", fontsize=12, weight="bold", color="#426080") + figure.text(.025, .925, f"{info['modality']} · {info['format']} · {name}", fontsize=23, weight="bold") + for index, (axis, (title, values)) in enumerate(zip(axes.flat, panels)): + limits = (0, 1) if values.dtype == np.bool_ else np.percentile(values, (1, 99)) + axis.imshow(values, cmap="gray", vmin=limits[0], vmax=limits[1]) + axis.set_title(f"{index:02d} {title}", fontsize=12, loc="left", pad=9) + axis.axis("off") + note = axes.flat[-1] + note.axis("off") + selection = f"Axis {info['slice_axis']} / slice {info['slice_index']}" if info['slice_index'] is not None else "Native 2D frame" + source_label = textwrap.fill(info['source'], width=30) + note.text(.02, .95, "CASE / PROCESSING", transform=note.transAxes, weight="bold", fontsize=12) + note.text(.02, .83, + f"Source: {source_label}\n{selection} · {info['files']} file(s)\nInput: {image.shape[1]} × {image.shape[0]}\n{info['orientation']}\n\n1–99% window → [0, 1] input.\nAll filters act on this 2D slice.\nIndependent output display contrast.\nCanny: boolean edges.", + transform=note.transAxes, fontsize=10, va="top", linespacing=1.55) + figure.savefig(output, dpi=150, metadata={"Title": f"Image Filter comparison: {name}"}) + plt.close(figure) + + +def render_overview(originals: dict[str, tuple[np.ndarray, dict]], output: Path) -> None: + rows = (len(originals) + 2) // 3 + with plt.rc_context(STYLE): + figure, axes = plt.subplots(rows, 3, figsize=(12, 4 * rows + 1), squeeze=False) + figure.subplots_adjust(left=.035, right=.965, bottom=.045, top=1 - 1.25 / (4 * rows + 1), + hspace=.25, wspace=.12) + figure.suptitle("Image Filter · medical imaging examples", fontsize=20, weight="bold", y=.975) + for axis in axes.flat: + axis.axis("off") + for axis, (name, (original, info)) in zip(axes.flat, originals.items()): + axis.imshow(original, cmap="gray", vmin=0, vmax=1) + axis.set_anchor("N") + axis.set_title(f"{info['modality']} · {info['format']}\n{name}-filters.png", fontsize=12, loc="left") + figure.savefig(output, dpi=150) + plt.close(figure) + + +def parse_args(argv: list[str] | None = None) -> argparse.Namespace: + parser = argparse.ArgumentParser(description=__doc__) + selection = parser.add_mutually_exclusive_group() + selection.add_argument("--all", action="store_true", help="All available local medical fixtures") + selection.add_argument("--images", nargs="+", choices=CASE_NAMES, help="Select local medical cases") + parser.add_argument("--input", type=Path, nargs="+", default=[], help="User NIfTI, NRRD or single DICOM paths") + parser.add_argument("--dicom-series", type=Path, nargs="+", default=[], help="DICOM series directories (one volume each)") + parser.add_argument("--series-uid", help="Select a series when a directory contains multiple series") + parser.add_argument("--modality", default="Unknown", help="Modality label for user NIfTI/NRRD files") + parser.add_argument("--slice-axis", type=int, choices=(0, 1, 2), default=2, help="Canonical voxel axis for volumes") + parser.add_argument("--slice-index", type=int, help="Override deterministic non-empty slice selection") + destination = parser.add_mutually_exclusive_group() + destination.add_argument("--output", type=Path, help="Single PNG output") + destination.add_argument("--output-dir", type=Path, default=None, help="Default: package example-output directory") + parser.add_argument("--contact-sheet", action="store_true", help="Write source-overview.png") + args = parser.parse_args(argv) + args.names = list(CASE_NAMES) if args.all else list(dict.fromkeys(args.images or ([] if args.input or args.dicom_series else ['mri']))) + count = len(args.names) + len(args.input) + len(args.dicom_series) + if args.output and (count != 1 or args.output.suffix.lower() != '.png'): + parser.error('--output requires exactly one source and a .png extension') + if args.output and args.output.name == 'source-overview.png': + parser.error('--output must not use the reserved overview filename') + return args + + +def main(argv: list[str] | None = None) -> None: + args = parse_args(argv) + catalog = bundled_cases() + selections = [] + for name in args.names: + if name not in catalog or not catalog[name][0].exists(): + if args.all: + print(f"Unavailable bundled case: {name}; no download attempted") + continue + raise SystemExit(f"Unavailable bundled case: {name}") + path, modality = catalog[name] + selections.append((name, path, modality)) + for path in args.input + args.dicom_series: + name = re.sub(r'[^a-z0-9_-]+', '-', path.name.lower().removesuffix('.gz').removesuffix('.nii').removesuffix('.nrrd').removesuffix('.dcm')).strip('-') + selections.append((name or 'case', path, args.modality)) + if len({name for name, _, _ in selections}) != len(selections): + raise SystemExit('Source basenames collide; run cases separately with distinct --output names') + if not selections: + raise SystemExit('No local medical fixtures available') + directory = args.output.parent if args.output else (args.output_dir or DEFAULT_DIR) + directory.mkdir(parents=True, exist_ok=True) + results = {} + for name, path, modality in selections: + try: + case = load_medical(path, modality, args.series_uid) + image, info = representative_plane(case, args.slice_axis, args.slice_index) + except (ValueError, FileNotFoundError, RuntimeError, ImportError) as exc: + raise SystemExit(f"{name}: {exc}") from exc + output = args.output or directory / f'{name}-filters.png' + render_panel(name, image, info, output) + results[name] = (image, info) + print(output.resolve()) + (directory / 'cases.json').write_text(json.dumps({name: info for name, (_, info) in results.items()}, indent=2) + '\n') + if args.contact_sheet: + output = directory / 'source-overview.png' + render_overview(results, output) + print(output.resolve()) + if args.all: + print('Bundled NRRD coverage unavailable; use --input /path/to/case.nrrd. No downloads performed.') + + +if __name__ == "__main__": + main() diff --git a/pyproject.toml b/pyproject.toml new file mode 100644 index 0000000..6009146 --- /dev/null +++ b/pyproject.toml @@ -0,0 +1,35 @@ +[build-system] +requires = ["setuptools>=68"] +build-backend = "setuptools.build_meta" + +[project] +name = "image-filter" +version = "0.1.0" +description = "Typed, application-independent scientific image filtering pipelines" +readme = "README.md" +requires-python = ">=3.10" +license = { text = "Proprietary" } +authors = [{ name = "Radiuma contributors" }] +dependencies = [ + "numpy>=1.24", + "scipy>=1.10", + "scikit-image>=0.21", + "PyWavelets>=1.4", +] + +[project.optional-dependencies] +test = ["pytest>=8", "nibabel>=5", "matplotlib>=3.7"] +examples = ["nibabel>=5", "matplotlib>=3.7", "pydicom>=3", "SimpleITK>=2.3", "pynrrd>=1.0"] + +[tool.setuptools.packages.find] +where = ["src"] + +[tool.setuptools.package-data] +image_filter = ["py.typed"] + +[tool.pytest.ini_options] +addopts = "-ra --strict-markers" +testpaths = ["tests"] +markers = [ + "real_data: tests which load non-generated image data", +] diff --git a/src/image_filter/__init__.py b/src/image_filter/__init__.py new file mode 100644 index 0000000..b975991 --- /dev/null +++ b/src/image_filter/__init__.py @@ -0,0 +1,36 @@ +"""Standalone scientific image filters for NumPy arrays.""" + +from .base import ImageFilter +from .errors import ( + DimensionInferenceError, + ImageFilterError, + InvalidConfigurationError, + InvalidImageError, + UnsupportedDimensionError, + UnsupportedOptionError, +) +from .filters import ( + CannyFilter, + GaborFilter, + GaussianFilter, + LaplacianFilter, + LawsFilter, + LoGFilter, + MeanFilter, + MedianFilter, + SobelFilter, + WaveletFilter, +) +from .model import Capabilities, FilterResult, ImageItem, SpatialLayout, infer_spatial_layout +from .pipeline import Pipeline, SliceWise + +__all__ = [ + "CannyFilter", "Capabilities", "DimensionInferenceError", "FilterResult", "GaborFilter", + "GaussianFilter", "LaplacianFilter", + "ImageFilter", "ImageFilterError", "ImageItem", "InvalidConfigurationError", + "InvalidImageError", "LawsFilter", "LoGFilter", "MeanFilter", "MedianFilter", "Pipeline", + "SliceWise", "SobelFilter", "SpatialLayout", "UnsupportedDimensionError", + "UnsupportedOptionError", "WaveletFilter", "infer_spatial_layout", +] + +__version__ = "0.1.0" diff --git a/src/image_filter/_ops.py b/src/image_filter/_ops.py new file mode 100644 index 0000000..efc949f --- /dev/null +++ b/src/image_filter/_ops.py @@ -0,0 +1,238 @@ +"""Stateless numerical operations used by public filter objects.""" + +from __future__ import annotations + +from typing import Any, Iterable, Iterator + +import numpy as np +import pywt +from numpy.typing import NDArray +from scipy import ndimage as ndi +from skimage.feature import canny +from skimage.filters import gabor + + +def axis_sizes(ndim: int, axes: tuple[int, ...], size: int) -> tuple[int, ...]: + return tuple(size if axis in axes else 1 for axis in range(ndim)) + + +def mean_filter( + image: NDArray[np.number], axes: tuple[int, ...], size: int, mode: str, + cval: float, origin: int, dtype: np.dtype, +) -> NDArray[np.floating]: + return ndi.uniform_filter( + image, size=axis_sizes(image.ndim, axes, size), mode=mode, + cval=cval, origin=origin, output=dtype, + ) + + +def _spatial_blocks( + image: NDArray[Any], axes: tuple[int, ...], +) -> Iterator[tuple[tuple[Any, ...], NDArray[Any]]]: + """Index away channels/singletons, preserving spatial order and shape.""" + excluded = tuple(axis for axis in range(image.ndim) if axis not in axes) + for indices in np.ndindex(*(image.shape[axis] for axis in excluded)): + key = [slice(None)] * image.ndim + for axis, index in zip(excluded, indices): + key[axis] = index + yield tuple(key), image[tuple(key)] + + +def gaussian_filter( + image: NDArray[Any], axes: tuple[int, ...], sigma: tuple[float, ...], + mode: str, cval: float, truncate: float, dtype: np.dtype, +) -> NDArray[np.floating]: + output = np.empty(image.shape, dtype=dtype) + for key, block in _spatial_blocks(image, axes): + output[key] = ndi.gaussian_filter( + np.asarray(block, dtype=dtype), sigma, output=dtype, + mode=mode, cval=cval, truncate=truncate, + ) + return output + + +def median_filter( + image: NDArray[Any], axes: tuple[int, ...], size: int, mode: str, + cval: float, origin: int, dtype: np.dtype, +) -> NDArray[np.floating]: + output = np.empty(image.shape, dtype=dtype) + for key, block in _spatial_blocks(image, axes): + output[key] = ndi.median_filter( + np.asarray(block, dtype=dtype), size, output=dtype, + mode=mode, cval=cval, origin=origin, + ) + return output + + +def sobel_filter( + image: NDArray[Any], axes: tuple[int, ...], spatial_axis: int | None, + mode: str, cval: float, dtype: np.dtype, +) -> NDArray[np.floating]: + output = np.empty(image.shape, dtype=dtype) + for key, block in _spatial_blocks(image, axes): + data = np.asarray(block, dtype=dtype) + if spatial_axis is not None: + result = ndi.sobel(data, axis=spatial_axis, mode=mode, cval=cval) + else: + result = np.zeros(data.shape, dtype=dtype) + for axis in range(data.ndim): + component = ndi.sobel(data, axis=axis, mode=mode, cval=cval) + # hypot avoids unnecessary overflow when squaring derivatives. + np.hypot(result, component, out=result) + output[key] = result + return output + + +def laplacian_filter( + image: NDArray[Any], axes: tuple[int, ...], mode: str, + cval: float, dtype: np.dtype, +) -> NDArray[np.floating]: + output = np.empty(image.shape, dtype=dtype) + for key, block in _spatial_blocks(image, axes): + output[key] = ndi.laplace( + np.asarray(block, dtype=dtype), output=dtype, mode=mode, cval=cval, + ) + return output + + +def log_filter( + image: NDArray[np.number], axes: tuple[int, ...], sigma: tuple[float, ...], + mode: str, cval: float, truncate: float, dtype: np.dtype, +) -> NDArray[np.floating]: + data = np.asarray(image, dtype=dtype) + result = np.zeros_like(data) + sigma_all = [0.0] * image.ndim + for axis, value in zip(axes, sigma, strict=True): + sigma_all[axis] = value + for axis in axes: + order = [0] * image.ndim + order[axis] = 2 + result += ndi.gaussian_filter( + data, sigma=sigma_all, order=order, mode=mode, cval=cval, truncate=truncate + ) + return result + + +LAWS_KERNELS: dict[str, NDArray[np.float64]] = { + "L5": np.array([1, 4, 6, 4, 1], dtype=np.float64), + "E5": np.array([-1, -2, 0, 2, 1], dtype=np.float64), + "S5": np.array([-1, 0, 2, 0, -1], dtype=np.float64), + "W5": np.array([-1, 2, 0, -2, 1], dtype=np.float64), + "R5": np.array([1, -4, 6, -4, 1], dtype=np.float64), + "L3": np.array([1, 2, 1], dtype=np.float64), + "E3": np.array([-1, 0, 1], dtype=np.float64), + "S3": np.array([-1, 2, -1], dtype=np.float64), +} + + +def parse_laws(configuration: str, dimension: int) -> tuple[NDArray[np.float64], ...]: + names = tuple(configuration[index:index + 2].upper() for index in range(0, len(configuration), 2)) + if len(names) != dimension or any(name not in LAWS_KERNELS for name in names): + expected = "".join(["L5"] * dimension) + raise ValueError(f"configuration must contain {dimension} valid kernels, for example {expected}") + return tuple(LAWS_KERNELS[name] / np.linalg.norm(LAWS_KERNELS[name]) for name in names) + + +def separable_convolve( + image: NDArray[np.number], axes: tuple[int, ...], kernels: Iterable[NDArray[np.floating]], + mode: str, cval: float, dtype: np.dtype, +) -> NDArray[np.floating]: + result = np.asarray(image, dtype=dtype) + for axis, kernel in zip(axes, kernels, strict=True): + result = ndi.convolve1d(result, kernel, axis=axis, mode=mode, cval=cval) + return result + + +def laws_filter( + image: NDArray[np.number], axes: tuple[int, ...], configuration: str, + mode: str, cval: float, energy_radius: int | None, dtype: np.dtype, +) -> NDArray[np.floating]: + response = separable_convolve(image, axes, parse_laws(configuration, len(axes)), mode, cval, dtype) + if energy_radius is None: + return response + size = 2 * energy_radius + 1 + return ndi.uniform_filter( + np.abs(response), size=axis_sizes(image.ndim, axes, size), mode=mode, + cval=cval, output=dtype, + ) + + +def _map_slices(image: NDArray[np.number], channel_axis: int | None): + if channel_axis is None: + yield (), image + return + for index in range(image.shape[channel_axis]): + key = [slice(None)] * image.ndim + key[channel_axis] = index + yield tuple(key), image[tuple(key)] + + +def gabor_filter( + image: NDArray[np.number], channel_axis: int | None, frequency: float, + orientations: tuple[float, ...], sigma: float, gamma: float, mode: str, + cval: float, response: str, pooling: str, dtype: np.dtype, +) -> NDArray[np.floating]: + output = np.empty(image.shape, dtype=dtype) + for key, plane in _map_slices(image, channel_axis): + responses: list[NDArray[np.floating]] = [] + for theta in orientations: + real, imag = gabor( + np.asarray(plane, dtype=dtype), frequency=frequency, theta=theta, + sigma_x=sigma, sigma_y=sigma / gamma, mode=mode, cval=cval, + ) + if response == "magnitude": + current = np.hypot(real, imag) + elif response == "real": + current = real + else: + current = imag + responses.append(current) + stack = np.stack(responses) + if pooling == "mean": + pooled = stack.mean(axis=0) + elif pooling == "sum": + pooled = stack.sum(axis=0) + elif pooling == "min": + pooled = stack.min(axis=0) + else: + pooled = stack.max(axis=0) + output[key] = pooled.astype(dtype, copy=False) + return output + + +def canny_filter( + image: NDArray[np.number], channel_axis: int | None, sigma: float, + low_threshold: float | None, high_threshold: float | None, + use_quantiles: bool, mode: str, cval: float, +) -> NDArray[np.bool_]: + """Apply scikit-image Canny independently to each declared channel.""" + + output = np.empty(image.shape, dtype=np.bool_) + for key, plane in _map_slices(image, channel_axis): + output[key] = canny( + np.squeeze(plane), sigma=sigma, low_threshold=low_threshold, + high_threshold=high_threshold, use_quantiles=use_quantiles, + mode=mode, cval=cval, + ).reshape(plane.shape) + return output + + +def wavelet_filter( + image: NDArray[np.number], axes: tuple[int, ...], configuration: str, + wavelet_name: str, level: int, pad_mode: str, normalize: bool, dtype: np.dtype, +) -> NDArray[np.floating]: + if len(configuration) != len(axes) or any(value not in "LHlh" for value in configuration): + raise ValueError(f"configuration must have {len(axes)} characters, each L or H") + multiple = 2 ** level + pad_width = [(0, 0)] * image.ndim + for axis in axes: + pad_width[axis] = (0, (-image.shape[axis]) % multiple) + padded = np.pad(np.asarray(image, dtype=dtype), pad_width, mode=pad_mode) + coefficients = pywt.swtn( + padded, wavelet=wavelet_name, level=level, axes=axes, + trim_approx=False, norm=normalize, + ) + key = configuration.lower().replace("l", "a").replace("h", "d") + result = coefficients[0][key] + crop = tuple(slice(0, size) for size in image.shape) + return np.asarray(result[crop], dtype=dtype) diff --git a/src/image_filter/base.py b/src/image_filter/base.py new file mode 100644 index 0000000..ef909f9 --- /dev/null +++ b/src/image_filter/base.py @@ -0,0 +1,101 @@ +"""Base filter protocol, validation, and provenance behavior.""" + +from __future__ import annotations + +from abc import ABC, abstractmethod +from dataclasses import dataclass, fields +from typing import Any, ClassVar, Iterable + +import numpy as np +from numpy.typing import NDArray + +from .errors import InvalidConfigurationError, InvalidImageError, UnsupportedDimensionError +from .model import Capabilities, FilterResult, ImageItem, Metadata, SpatialLayout, infer_spatial_layout, mutable_metadata_copy, readonly_metadata + +BOUNDARY_MODES = frozenset({"reflect", "constant", "nearest", "mirror", "wrap"}) +OUTPUT_DTYPES = frozenset({"float32", "float64"}) + + +@dataclass(frozen=True, slots=True) +class ImageFilter(ABC): + """Component interface shared by leaves, adapters, and composite pipelines.""" + + capabilities: ClassVar[Capabilities] + + @abstractmethod + def apply(self, image: NDArray[Any], metadata: Metadata | None = None) -> FilterResult: + """Return a filtered image and copied processing metadata.""" + + def apply_batch(self, items: Iterable[ImageItem]) -> list[FilterResult]: + return [self.apply(item.image, item.metadata) for item in items] + + +@dataclass(frozen=True, slots=True) +class NumericFilter(ImageFilter): + """Shared validation and provenance implementation for numerical leaves.""" + + def apply(self, image: NDArray[Any], metadata: Metadata | None = None) -> FilterResult: + layout = infer_spatial_layout(image, metadata) + return self._apply_with_layout(image, metadata, layout) + + def _apply_with_layout( + self, image: NDArray[Any], metadata: Metadata | None, layout: SpatialLayout + ) -> FilterResult: + if layout.dimension not in self.capabilities.spatial_dimensions: + supported = ", ".join(f"{value}D" for value in sorted(self.capabilities.spatial_dimensions)) + raise UnsupportedDimensionError( + f"{type(self).__name__} supports {supported}, received {layout.dimension}D; " + "wrap it in SliceWise(filter, axis=...) for explicit 3D slice processing" + ) + if not np.isfinite(image).all(): + raise InvalidImageError("image contains NaN or infinite values") + output = self._compute(image, layout) + if output.shape != image.shape and self.capabilities.preserves_shape: + raise RuntimeError(f"{type(self).__name__} violated its shape-preservation contract") + copied = mutable_metadata_copy(metadata) + history = list(copied.get("processing_history", [])) + history.append(self._history_entry(image, output, layout)) + copied["processing_history"] = history + return FilterResult(output, readonly_metadata(copied)) + + def _history_entry( + self, image: NDArray[Any], output: NDArray[Any], layout: SpatialLayout + ) -> dict[str, Any]: + def serialize(filter_object: ImageFilter) -> dict[str, Any]: + config: dict[str, Any] = {} + for item in fields(filter_object): + value = getattr(filter_object, item.name) + if isinstance(value, np.dtype): + value = value.name + elif isinstance(value, tuple): + value = [ + {"filter": type(child).__name__, "parameters": serialize(child)} + if isinstance(child, ImageFilter) else child for child in value + ] + elif isinstance(value, ImageFilter): + value = {"filter": type(value).__name__, "parameters": serialize(value)} + config[item.name] = value + return config + + return { + "filter": type(self).__name__, + "parameters": serialize(self), + "spatial_dimension": layout.dimension, + "spatial_axes": list(layout.spatial_axes), + "input_shape": list(image.shape), + "output_shape": list(output.shape), + "input_dtype": image.dtype.name, + "output_dtype": output.dtype.name, + } + + @abstractmethod + def _compute(self, image: NDArray[Any], layout: SpatialLayout) -> NDArray[Any]: + """Compute a same-shape result after public validation.""" + + +def validate_common(mode: str, output_dtype: str) -> np.dtype: + if mode not in BOUNDARY_MODES: + raise InvalidConfigurationError(f"mode must be one of {sorted(BOUNDARY_MODES)}, got {mode!r}") + if output_dtype not in OUTPUT_DTYPES: + raise InvalidConfigurationError("output_dtype must be 'float32' or 'float64'") + return np.dtype(output_dtype) diff --git a/src/image_filter/errors.py b/src/image_filter/errors.py new file mode 100644 index 0000000..6e76731 --- /dev/null +++ b/src/image_filter/errors.py @@ -0,0 +1,25 @@ +"""Domain-specific exceptions for :mod:`image_filter`.""" + + +class ImageFilterError(Exception): + """Base class for errors intentionally exposed by the public API.""" + + +class InvalidImageError(ImageFilterError, ValueError): + """The image array or its metadata violates the input contract.""" + + +class DimensionInferenceError(InvalidImageError): + """Spatial dimensionality cannot be determined without guessing.""" + + +class UnsupportedDimensionError(ImageFilterError, ValueError): + """A filter does not support the inferred spatial dimensionality.""" + + +class InvalidConfigurationError(ImageFilterError, ValueError): + """Filter configuration is invalid.""" + + +class UnsupportedOptionError(ImageFilterError, NotImplementedError): + """A recognized advanced option has no verified implementation.""" diff --git a/src/image_filter/filters.py b/src/image_filter/filters.py new file mode 100644 index 0000000..13785fc --- /dev/null +++ b/src/image_filter/filters.py @@ -0,0 +1,421 @@ +"""Production filter configurations backed by scientific Python libraries.""" + +from __future__ import annotations + +from dataclasses import dataclass +from typing import Any, ClassVar +from numbers import Real + +import numpy as np +import pywt +from numpy.typing import NDArray + +from . import _ops +from .base import NumericFilter, validate_common +from .errors import InvalidConfigurationError, InvalidImageError, UnsupportedOptionError +from .model import Capabilities, SpatialLayout + +_GENERAL = Capabilities(frozenset({2, 3})) +_TWO_DIMENSIONAL = Capabilities(frozenset({2})) +_BOOLEAN_TWO_DIMENSIONAL = Capabilities( + frozenset({2}), + output_dtype="boolean edge mask", + accepted_dtypes="boolean, integers up to 32 bits, or real floating-point", +) + + +def _finite_real(name: str, value: Any) -> None: + if isinstance(value, (bool, np.bool_)) or not isinstance(value, Real) or not np.isfinite(value): + raise InvalidConfigurationError(f"{name} must be a finite real number") + + +def _validate_numeric(mode: str, output_dtype: str, cval: float) -> None: + if not isinstance(mode, str) or not isinstance(output_dtype, str): + raise InvalidConfigurationError("mode and output_dtype must be strings") + validate_common(mode, output_dtype) + _finite_real("cval", cval) + + +@dataclass(frozen=True, slots=True) +class MeanFilter(NumericFilter): + """N-dimensional box mean using :func:`scipy.ndimage.uniform_filter`. + + ``origin=0`` uses SciPy's centered convention. For an even kernel, the + footprint is asymmetric by one sample: changing origin explicitly shifts + that footprint. Only inferred spatial axes are filtered. + """ + + size: int = 3 + mode: str = "reflect" + cval: float = 0.0 + origin: int = 0 + output_dtype: str = "float64" + capabilities = _GENERAL + + def __post_init__(self) -> None: + validate_common(self.mode, self.output_dtype) + if not isinstance(self.size, int) or isinstance(self.size, bool) or self.size < 1: + raise InvalidConfigurationError("size must be a positive integer") + low, high = -(self.size // 2), (self.size - 1) // 2 + if not isinstance(self.origin, int) or not low <= self.origin <= high: + raise InvalidConfigurationError( + f"origin must be an integer in [{low}, {high}] for size={self.size}" + ) + + def _compute(self, image: NDArray[Any], layout: SpatialLayout) -> NDArray[Any]: + return _ops.mean_filter( + image, layout.spatial_axes, self.size, self.mode, self.cval, + self.origin, np.dtype(self.output_dtype), + ) + + +@dataclass(frozen=True, slots=True) +class GaussianFilter(NumericFilter): + """N-dimensional Gaussian smoothing using ``scipy.ndimage``. + + Sigma is measured in pixels and may be scalar or one value per inferred + spatial axis. Channel and retained singleton axes are not smoothed. + """ + + sigma: float | tuple[float, ...] = 1.0 + truncate: float = 4.0 + mode: str = "reflect" + cval: float = 0.0 + output_dtype: str = "float64" + capabilities = _GENERAL + + def __post_init__(self) -> None: + _validate_numeric(self.mode, self.output_dtype, self.cval) + if isinstance(self.sigma, Real): + values = (self.sigma,) + elif isinstance(self.sigma, tuple) and self.sigma: + values = self.sigma + else: + raise InvalidConfigurationError("sigma must be a scalar or nonempty tuple") + for value in values: + _finite_real("sigma", value) + if value <= 0: + raise InvalidConfigurationError("sigma must be greater than zero") + _finite_real("truncate", self.truncate) + if self.truncate <= 0: + raise InvalidConfigurationError("truncate must be greater than zero") + + def _compute(self, image: NDArray[Any], layout: SpatialLayout) -> NDArray[Any]: + if isinstance(self.sigma, Real): + sigma = (float(self.sigma),) * layout.dimension + else: + sigma = tuple(float(value) for value in self.sigma) + if len(sigma) != layout.dimension: + raise InvalidConfigurationError( + f"sigma contains {len(sigma)} values but this image is {layout.dimension}D" + ) + return _ops.gaussian_filter( + image, layout.spatial_axes, sigma, self.mode, self.cval, + self.truncate, np.dtype(self.output_dtype), + ) + + +@dataclass(frozen=True, slots=True) +class MedianFilter(NumericFilter): + """N-dimensional median denoising using ``scipy.ndimage``. + + The scalar window size applies to spatial axes only. ``origin`` follows + SciPy's alignment convention and makes even-sized windows explicit. + """ + + size: int = 3 + mode: str = "reflect" + cval: float = 0.0 + origin: int = 0 + output_dtype: str = "float64" + capabilities = _GENERAL + + def __post_init__(self) -> None: + _validate_numeric(self.mode, self.output_dtype, self.cval) + if not isinstance(self.size, int) or isinstance(self.size, bool) or self.size < 1: + raise InvalidConfigurationError("size must be a positive integer") + low, high = -(self.size // 2), (self.size - 1) // 2 + if not isinstance(self.origin, int) or isinstance(self.origin, bool) or not low <= self.origin <= high: + raise InvalidConfigurationError( + f"origin must be an integer in [{low}, {high}] for size={self.size}" + ) + + def _compute(self, image: NDArray[Any], layout: SpatialLayout) -> NDArray[Any]: + return _ops.median_filter( + image, layout.spatial_axes, self.size, self.mode, self.cval, + self.origin, np.dtype(self.output_dtype), + ) + + +@dataclass(frozen=True, slots=True) +class SobelFilter(NumericFilter): + """Spatial Sobel derivative or gradient magnitude. + + ``axis=None`` returns Euclidean magnitude across all spatial axes. An + integer selects an index within the inferred spatial-axis order, not an + absolute array axis, so channels and singleton axes remain excluded. + """ + + axis: int | None = None + mode: str = "reflect" + cval: float = 0.0 + output_dtype: str = "float64" + capabilities = _GENERAL + + def __post_init__(self) -> None: + _validate_numeric(self.mode, self.output_dtype, self.cval) + if self.axis is not None and ( + not isinstance(self.axis, int) or isinstance(self.axis, bool) + ): + raise InvalidConfigurationError("axis must be an integer or None") + + if self.axis is not None and not -3 <= self.axis <= 2: + raise InvalidConfigurationError("axis must index a 2D or 3D spatial axis") + + def _compute(self, image: NDArray[Any], layout: SpatialLayout) -> NDArray[Any]: + spatial_axis = self.axis + if spatial_axis is not None: + if spatial_axis < 0: + spatial_axis += layout.dimension + if not 0 <= spatial_axis < layout.dimension: + raise InvalidConfigurationError( + f"axis {self.axis} is invalid for a {layout.dimension}D image" + ) + return _ops.sobel_filter( + image, layout.spatial_axes, spatial_axis, self.mode, self.cval, + np.dtype(self.output_dtype), + ) + + +@dataclass(frozen=True, slots=True) +class LaplacianFilter(NumericFilter): + """N-dimensional discrete Laplacian using ``scipy.ndimage``.""" + + mode: str = "reflect" + cval: float = 0.0 + output_dtype: str = "float64" + capabilities = _GENERAL + + def __post_init__(self) -> None: + _validate_numeric(self.mode, self.output_dtype, self.cval) + + def _compute(self, image: NDArray[Any], layout: SpatialLayout) -> NDArray[Any]: + return _ops.laplacian_filter( + image, layout.spatial_axes, self.mode, self.cval, + np.dtype(self.output_dtype), + ) + + +@dataclass(frozen=True, slots=True) +class LoGFilter(NumericFilter): + """Laplacian of Gaussian, summed over spatial axes in pixel units.""" + + sigma: float | tuple[float, ...] = 1.0 + truncate: float = 4.0 + mode: str = "reflect" + cval: float = 0.0 + output_dtype: str = "float64" + capabilities = _GENERAL + + def __post_init__(self) -> None: + validate_common(self.mode, self.output_dtype) + values = (self.sigma,) if isinstance(self.sigma, (int, float)) else self.sigma + if not values or any(not np.isfinite(value) or value <= 0 for value in values): + raise InvalidConfigurationError("sigma values must be finite and greater than zero") + if not np.isfinite(self.truncate) or self.truncate <= 0: + raise InvalidConfigurationError("truncate must be finite and greater than zero") + + def _compute(self, image: NDArray[Any], layout: SpatialLayout) -> NDArray[Any]: + if isinstance(self.sigma, (int, float)): + sigma = (float(self.sigma),) * layout.dimension + else: + sigma = tuple(float(value) for value in self.sigma) + if len(sigma) != layout.dimension: + raise InvalidConfigurationError( + f"sigma contains {len(sigma)} values but this image is {layout.dimension}D" + ) + return _ops.log_filter( + image, layout.spatial_axes, sigma, self.mode, self.cval, + self.truncate, np.dtype(self.output_dtype), + ) + + +@dataclass(frozen=True, slots=True) +class LawsFilter(NumericFilter): + """Separable Laws texture response with optional local absolute energy.""" + + configuration: str = "L5E5" + energy_radius: int | None = None + rotation_invariant: bool = False + mode: str = "reflect" + cval: float = 0.0 + output_dtype: str = "float64" + capabilities = _GENERAL + + def __post_init__(self) -> None: + validate_common(self.mode, self.output_dtype) + if self.energy_radius is not None and ( + not isinstance(self.energy_radius, int) or isinstance(self.energy_radius, bool) + or self.energy_radius < 0 + ): + raise InvalidConfigurationError("energy_radius must be a non-negative integer or None") + if self.rotation_invariant: + raise UnsupportedOptionError( + "rotation_invariant Laws pooling is not implemented; use an explicit pipeline of orientations" + ) + + def _compute(self, image: NDArray[Any], layout: SpatialLayout) -> NDArray[Any]: + try: + return _ops.laws_filter( + image, layout.spatial_axes, self.configuration, self.mode, + self.cval, self.energy_radius, np.dtype(self.output_dtype), + ) + except ValueError as exc: + raise InvalidConfigurationError(str(exc)) from exc + + +@dataclass(frozen=True, slots=True) +class GaborFilter(NumericFilter): + """2D Gabor response using :func:`skimage.filters.gabor`. + + Angles are radians. Multiple orientations are pooled explicitly. For 3D + images use :class:`~image_filter.SliceWise` and choose a slicing axis. + """ + + frequency: float = 0.25 + orientations: tuple[float, ...] = (0.0,) + sigma: float = 2.0 + gamma: float = 1.0 + response: str = "magnitude" + pooling: str = "max" + mode: str = "reflect" + cval: float = 0.0 + output_dtype: str = "float64" + capabilities = _TWO_DIMENSIONAL + + def __post_init__(self) -> None: + validate_common(self.mode, self.output_dtype) + if not np.isfinite(self.frequency) or not 0 < self.frequency < 0.5: + raise InvalidConfigurationError("frequency must be finite and in the open interval (0, 0.5)") + if not self.orientations or any(not np.isfinite(value) for value in self.orientations): + raise InvalidConfigurationError("orientations must contain finite angles in radians") + if not np.isfinite(self.sigma) or self.sigma <= 0: + raise InvalidConfigurationError("sigma must be finite and greater than zero") + if not np.isfinite(self.gamma) or self.gamma <= 0: + raise InvalidConfigurationError("gamma must be finite and greater than zero") + if self.response not in {"magnitude", "real", "imaginary"}: + raise InvalidConfigurationError("response must be magnitude, real, or imaginary") + if self.pooling not in {"max", "min", "mean", "sum"}: + raise InvalidConfigurationError("pooling must be max, min, mean, or sum") + + def _compute(self, image: NDArray[Any], layout: SpatialLayout) -> NDArray[Any]: + return _ops.gabor_filter( + image, layout.channel_axis, self.frequency, tuple(self.orientations), + self.sigma, self.gamma, self.mode, self.cval, self.response, + self.pooling, np.dtype(self.output_dtype), + ) + + +@dataclass(frozen=True, slots=True) +class CannyFilter(NumericFilter): + """2D Canny edges using :func:`skimage.feature.canny`. + + Output is always boolean. With an explicit channel axis, channels are + processed independently. Volumes require an explicit ``SliceWise`` policy. + """ + + sigma: float = 1.0 + low_threshold: float | None = None + high_threshold: float | None = None + use_quantiles: bool = False + mode: str = "constant" + cval: float = 0.0 + output_dtype: ClassVar[str] = "bool" + capabilities = _BOOLEAN_TWO_DIMENSIONAL + + def __post_init__(self) -> None: + if not isinstance(self.mode, str) or self.mode not in {"reflect", "constant", "nearest", "mirror", "wrap"}: + raise InvalidConfigurationError( + "mode must be reflect, constant, nearest, mirror, or wrap" + ) + _finite_real("cval", self.cval) + _finite_real("sigma", self.sigma) + if not isinstance(self.use_quantiles, bool): + raise InvalidConfigurationError("use_quantiles must be boolean") + if not np.isfinite(self.sigma) or self.sigma <= 0: + raise InvalidConfigurationError("sigma must be finite and greater than zero") + for name, value in ( + ("low_threshold", self.low_threshold), + ("high_threshold", self.high_threshold), + ): + if value is not None: + _finite_real(name, value) + if value < 0: + raise InvalidConfigurationError(f"{name} must be nonnegative") + if self.use_quantiles and value is not None and not 0 <= value <= 1: + raise InvalidConfigurationError(f"{name} must be in [0, 1] when use_quantiles=True") + if ( + (0.1 if self.low_threshold is None else self.low_threshold) + > (0.2 if self.high_threshold is None else self.high_threshold) + and (self.use_quantiles or (self.low_threshold is not None and self.high_threshold is not None)) + ): + raise InvalidConfigurationError("low_threshold must not exceed high_threshold") + + def _compute(self, image: NDArray[Any], layout: SpatialLayout) -> NDArray[Any]: + if image.dtype.kind in "iu" and image.dtype.itemsize == 8: + raise InvalidImageError( + "Canny does not support 64-bit integers; explicitly convert the image " + "to floating point and choose thresholds in those units" + ) + try: + return _ops.canny_filter( + image, layout.channel_axis, self.sigma, self.low_threshold, + self.high_threshold, self.use_quantiles, self.mode, self.cval, + ) + except ValueError as exc: + raise InvalidConfigurationError(str(exc)) from exc + + + +@dataclass(frozen=True, slots=True) +class WaveletFilter(NumericFilter): + """Same-shape stationary wavelet coefficient from ``pywt.swtn``. + + ``configuration`` uses ``L`` (approximation) and ``H`` (detail), one per + inferred spatial axis. Arrays are padded at their high edge to a multiple + of ``2**level`` and the selected final-level coefficient is cropped back. + """ + + configuration: str = "LL" + wavelet: str = "db1" + level: int = 1 + pad_mode: str = "reflect" + normalize: bool = False + rotation_invariant: bool = False + output_dtype: str = "float64" + capabilities = _GENERAL + + def __post_init__(self) -> None: + if self.output_dtype not in {"float32", "float64"}: + raise InvalidConfigurationError("output_dtype must be 'float32' or 'float64'") + if not isinstance(self.level, int) or isinstance(self.level, bool) or self.level < 1: + raise InvalidConfigurationError("level must be a positive integer") + if self.pad_mode not in {"constant", "edge", "reflect", "symmetric", "wrap"}: + raise InvalidConfigurationError("pad_mode is not a supported numpy.pad mode") + try: + pywt.Wavelet(self.wavelet) + except ValueError as exc: + raise InvalidConfigurationError(f"unknown discrete wavelet {self.wavelet!r}") from exc + if self.rotation_invariant: + raise UnsupportedOptionError( + "rotation_invariant coefficient pooling is not implemented; select coefficients explicitly" + ) + + def _compute(self, image: NDArray[Any], layout: SpatialLayout) -> NDArray[Any]: + try: + return _ops.wavelet_filter( + image, layout.spatial_axes, self.configuration, self.wavelet, + self.level, self.pad_mode, self.normalize, np.dtype(self.output_dtype), + ) + except (ValueError, TypeError) as exc: + raise InvalidConfigurationError(str(exc)) from exc diff --git a/src/image_filter/model.py b/src/image_filter/model.py new file mode 100644 index 0000000..05c3871 --- /dev/null +++ b/src/image_filter/model.py @@ -0,0 +1,133 @@ +"""Immutable public domain objects and spatial-layout inference.""" + +from __future__ import annotations + +from copy import deepcopy +from dataclasses import dataclass, field +from types import MappingProxyType +from typing import Any, Mapping + +import numpy as np +from numpy.typing import NDArray + +from .errors import DimensionInferenceError, InvalidImageError + +Metadata = Mapping[str, Any] + + +@dataclass(frozen=True, slots=True) +class Capabilities: + """Machine-readable promises made by a filter implementation.""" + + spatial_dimensions: frozenset[int] + accepted_dtypes: str = "boolean, integer, or real floating-point" + input_range: str = "unrestricted finite numeric values" + output_dtype: str = "configured floating-point dtype" + preserves_shape: bool = True + supports_channels: bool = True + + +@dataclass(frozen=True, slots=True) +class SpatialLayout: + """The axes participating in filtering for one image.""" + + dimension: int + spatial_axes: tuple[int, ...] + singleton_axes: tuple[int, ...] + channel_axis: int | None + + +@dataclass(frozen=True, slots=True) +class ImageItem: + """One in-memory image and caller-owned descriptive metadata.""" + + image: NDArray[Any] + metadata: Metadata = field(default_factory=dict) + + +@dataclass(frozen=True, slots=True) +class FilterResult: + """A filtered image and a new metadata mapping containing provenance.""" + + image: NDArray[Any] + metadata: Mapping[str, Any] + + +def mutable_metadata_copy(metadata: Metadata | None) -> dict[str, Any]: + """Deep-copy metadata so nested caller values cannot be mutated.""" + + if metadata is None: + return {} + if not isinstance(metadata, Mapping): + raise InvalidImageError("metadata must be a mapping") + return deepcopy(dict(metadata)) + + +def readonly_metadata(metadata: Mapping[str, Any]) -> Mapping[str, Any]: + """Expose a read-only top-level metadata mapping.""" + + return MappingProxyType(dict(metadata)) + + +def infer_spatial_layout(image: NDArray[Any], metadata: Metadata | None = None) -> SpatialLayout: + """Infer 2D versus 3D independently for a single image. + + Singleton spatial axes are ignored for the effective dimension but retained + in the array. ``metadata['channel_axis']`` excludes a color/channel axis. + ``metadata['spatial_dims']`` may explicitly disambiguate a short 3D volume. + """ + + if not isinstance(image, np.ndarray): + raise InvalidImageError("image must be a numpy.ndarray") + if image.ndim < 2 or image.ndim > 4: + raise DimensionInferenceError( + f"expected 2 to 4 array axes, received shape {image.shape}" + ) + if image.dtype.kind not in "buif": + raise InvalidImageError(f"image dtype {image.dtype} is not a real numeric dtype") + + meta = metadata or {} + raw_channel = meta.get("channel_axis") + channel_axis: int | None = None + if raw_channel is not None: + if not isinstance(raw_channel, (int, np.integer)): + raise DimensionInferenceError("metadata['channel_axis'] must be an integer") + channel_axis = int(raw_channel) + if channel_axis < 0: + channel_axis += image.ndim + if channel_axis < 0 or channel_axis >= image.ndim: + raise DimensionInferenceError( + f"channel_axis {raw_channel} is invalid for shape {image.shape}" + ) + if image.shape[channel_axis] < 1: + raise DimensionInferenceError("channel axis must not be empty") + + candidate_axes = tuple(axis for axis in range(image.ndim) if axis != channel_axis) + active_axes = tuple(axis for axis in candidate_axes if image.shape[axis] > 1) + singleton_axes = tuple(axis for axis in candidate_axes if image.shape[axis] == 1) + if image.ndim == 4 and channel_axis is None: + raise DimensionInferenceError( + "4D input requires metadata['channel_axis']; time-series dimensions are not inferred" + ) + explicit = meta.get("spatial_dims") + if explicit is not None: + if explicit not in (2, 3): + raise DimensionInferenceError("metadata['spatial_dims'] must be 2 or 3") + if len(active_axes) != explicit: + raise DimensionInferenceError( + f"spatial_dims={explicit} conflicts with {len(active_axes)} non-singleton " + f"spatial axes in shape {image.shape}" + ) + + if len(active_axes) not in (2, 3): + raise DimensionInferenceError( + f"shape {image.shape} has {len(active_axes)} non-singleton spatial axes; " + "exactly 2 or 3 are required" + ) + if image.ndim == 3 and channel_axis is None and explicit is None: + if image.shape[-1] in (3, 4) and not singleton_axes: + raise DimensionInferenceError( + f"shape {image.shape} is ambiguous between a short volume and a color image; " + "set metadata['channel_axis']=2 for color or metadata['spatial_dims']=3 for a volume" + ) + return SpatialLayout(len(active_axes), active_axes, singleton_axes, channel_axis) diff --git a/src/image_filter/pipeline.py b/src/image_filter/pipeline.py new file mode 100644 index 0000000..1eff80a --- /dev/null +++ b/src/image_filter/pipeline.py @@ -0,0 +1,93 @@ +"""Ordered pipelines, explicit slice-wise adaptation, and batch execution.""" + +from __future__ import annotations + +from dataclasses import dataclass +from typing import Any, Iterable + +import numpy as np +from numpy.typing import NDArray + +from .base import ImageFilter, NumericFilter +from .errors import InvalidConfigurationError, UnsupportedDimensionError +from .model import Capabilities, FilterResult, Metadata, SpatialLayout, infer_spatial_layout + + +@dataclass(frozen=True, slots=True) +class SliceWise(NumericFilter): + """Apply an intrinsically 2D filter across a caller-selected volume axis.""" + + filter: ImageFilter + axis: int + capabilities = Capabilities(frozenset({3})) + + def __post_init__(self) -> None: + if not isinstance(self.filter, ImageFilter): + raise InvalidConfigurationError("SliceWise requires an ImageFilter") + if 2 not in self.filter.capabilities.spatial_dimensions: + raise InvalidConfigurationError("SliceWise requires a filter with 2D capability") + if not isinstance(self.axis, int) or isinstance(self.axis, bool): + raise InvalidConfigurationError("axis must be an integer") + + def apply(self, image: NDArray[Any], metadata: Metadata | None = None) -> FilterResult: + # The explicit slicing axis resolves a 3-axis array that might otherwise + # be ambiguous with channel-last color data. + inference_metadata = dict(metadata or {}) + if image.ndim == 3 and "channel_axis" not in inference_metadata: + inference_metadata["spatial_dims"] = 3 + layout = infer_spatial_layout(image, inference_metadata) + return self._apply_with_layout(image, metadata, layout) + + def _compute(self, image: NDArray[Any], layout: SpatialLayout) -> NDArray[Any]: + axis = self.axis + image.ndim if self.axis < 0 else self.axis + if axis not in layout.spatial_axes: + raise UnsupportedDimensionError( + f"slice axis {self.axis} must be one of spatial axes {layout.spatial_axes}" + ) + output = None + for index in range(image.shape[axis]): + key = [slice(None)] * image.ndim + key[axis] = index + plane = image[tuple(key)] + def remap(value: int) -> int: + return value - int(value > axis) + plane_metadata = {"spatial_dims": 2} + if layout.channel_axis is not None: + plane_metadata["channel_axis"] = remap(layout.channel_axis) + filtered = self.filter.apply(plane, plane_metadata).image + if output is None: + output = np.empty(image.shape, dtype=filtered.dtype) + output[tuple(key)] = filtered + return output + + +@dataclass(frozen=True, slots=True) +class Pipeline(ImageFilter): + """Composite: an immutable ordered sequence of leaves or other pipelines.""" + + steps: tuple[ImageFilter, ...] + + def __init__(self, steps: Iterable[ImageFilter]) -> None: + values = tuple(steps) + if not values: + raise InvalidConfigurationError("a pipeline must contain at least one filter") + if any(not isinstance(step, ImageFilter) for step in values): + raise InvalidConfigurationError("every pipeline step must be an ImageFilter") + object.__setattr__(self, "steps", values) + + @property + def capabilities(self) -> Capabilities: + return Capabilities( + spatial_dimensions=frozenset.intersection(*(step.capabilities.spatial_dimensions for step in self.steps)), + accepted_dtypes="subject to each ordered child", + input_range="subject to each ordered child", + output_dtype=self.steps[-1].capabilities.output_dtype, + preserves_shape=all(step.capabilities.preserves_shape for step in self.steps), + supports_channels=all(step.capabilities.supports_channels for step in self.steps), + ) + + def apply(self, image: NDArray[Any], metadata: Metadata | None = None) -> FilterResult: + result = FilterResult(image, metadata or {}) + for step in self.steps: + result = step.apply(result.image, result.metadata) + return result diff --git a/src/image_filter/py.typed b/src/image_filter/py.typed new file mode 100644 index 0000000..8b13789 --- /dev/null +++ b/src/image_filter/py.typed @@ -0,0 +1 @@ + diff --git a/tests/test_additional_filters.py b/tests/test_additional_filters.py new file mode 100644 index 0000000..ddfae41 --- /dev/null +++ b/tests/test_additional_filters.py @@ -0,0 +1,259 @@ +from __future__ import annotations + +import numpy as np +import pytest +from scipy import ndimage as ndi +from skimage.feature import canny + +from image_filter import ( + CannyFilter, + GaussianFilter, + InvalidConfigurationError, + LaplacianFilter, + MedianFilter, + Pipeline, + SliceWise, + SobelFilter, + UnsupportedDimensionError, +) + + +@pytest.mark.parametrize("dimension", [2, 3]) +def test_gaussian_matches_scipy_for_2d_and_3d(dimension): + shape = (7, 8) if dimension == 2 else (6, 7, 8) + image = np.arange(np.prod(shape), dtype=np.int16).reshape(shape) + sigma = tuple(0.6 + 0.2 * index for index in range(dimension)) + metadata = None if dimension == 2 else {"spatial_dims": 3} + actual = GaussianFilter( + sigma=sigma, truncate=3.0, mode="mirror", output_dtype="float32" + ).apply(image, metadata).image + expected = ndi.gaussian_filter( + image, sigma=sigma, truncate=3.0, mode="mirror", output=np.float32 + ) + np.testing.assert_allclose(actual, expected, rtol=0, atol=0) + assert actual.dtype == np.float32 + + +@pytest.mark.parametrize("size,mode", [(3, "reflect"), (4, "constant")]) +def test_median_matches_scipy_for_window_and_boundary_modes(size, mode): + image = np.array( + [[9, 1, 7, 2, 8], [4, 6, 0, 5, 3], [8, 2, 9, 1, 7], [0, 5, 3, 6, 4]], + dtype=np.int16, + ) + origin = -1 if size == 4 else 0 + actual = MedianFilter( + size=size, origin=origin, mode=mode, cval=-3, output_dtype="float64" + ).apply(image).image + expected = ndi.median_filter( + image, size=size, origin=origin, mode=mode, cval=-3, output=np.float64 + ) + np.testing.assert_array_equal(actual, expected) + + +@pytest.mark.parametrize("axis", [1, None]) +def test_sobel_matches_scipy_components_and_magnitude(axis): + image = np.arange(8 * 9, dtype=np.float64).reshape(8, 9) ** 1.2 + actual = SobelFilter(axis=axis, mode="nearest").apply(image).image + components = [ndi.sobel(image, spatial_axis, mode="nearest", output=np.float64) for spatial_axis in (0, 1)] + expected = components[axis] if axis is not None else np.hypot(*components) + np.testing.assert_allclose(actual, expected, rtol=1e-15, atol=1e-15) + + +def test_sobel_matches_scipy_in_3d(): + image = np.arange(5 * 6 * 7, dtype=np.float32).reshape(5, 6, 7) + actual = SobelFilter(axis=-1, mode="wrap", output_dtype="float32").apply( + image, {"spatial_dims": 3} + ).image + expected = ndi.sobel(image, axis=2, mode="wrap", output=np.float32) + np.testing.assert_array_equal(actual, expected) + + +@pytest.mark.parametrize("dimension", [2, 3]) +def test_laplacian_matches_scipy_for_2d_and_3d(dimension): + shape = (8, 9) if dimension == 2 else (5, 6, 7) + image = np.linspace(-2, 3, np.prod(shape), dtype=np.float64).reshape(shape) ** 3 + metadata = None if dimension == 2 else {"spatial_dims": 3} + actual = LaplacianFilter(mode="constant", cval=1.5).apply(image, metadata).image + expected = ndi.laplace(image, mode="constant", cval=1.5, output=np.float64) + np.testing.assert_array_equal(actual, expected) + + +def test_canny_matches_scikit_image_reference(): + image = np.zeros((40, 44), dtype=np.float32) + image[9:31, 11:35] = 1 + configured = CannyFilter( + sigma=1.4, low_threshold=0.1, high_threshold=0.3, + mode="reflect", + ) + actual = configured.apply(image).image + expected = canny( + image, sigma=1.4, low_threshold=0.1, high_threshold=0.3, + mode="reflect", + ) + np.testing.assert_array_equal(actual, expected) + assert actual.dtype == np.bool_ + + +NEW_FILTERS = [GaussianFilter(.8), MedianFilter(4, origin=-1), SobelFilter(), LaplacianFilter(), CannyFilter(.8, .1, .3, True)] + + +def reference(configured, image): + if isinstance(configured, CannyFilter): + return canny(image, sigma=.8, low_threshold=.1, high_threshold=.3, use_quantiles=True) + image = image.astype(np.float64) + if isinstance(configured, GaussianFilter): + return ndi.gaussian_filter(image, .8) + if isinstance(configured, MedianFilter): + return ndi.median_filter(image, 4, origin=-1) + if isinstance(configured, SobelFilter): + return np.sqrt(sum(ndi.sobel(image, axis=i) ** 2 for i in range(image.ndim))) + return ndi.laplace(image) + + +# Representative layouts exercise each spatial-block implementation without a +# filter × singleton × channel Cartesian product. Canny has its own mapper. +@pytest.mark.parametrize('configured,singleton_axis,channel_axis', [ + (NEW_FILTERS[0], 0, None), (NEW_FILTERS[1], 1, 0), + (NEW_FILTERS[2], 2, -1), (NEW_FILTERS[3], 0, -1), + (NEW_FILTERS[4], 1, 0), (NEW_FILTERS[4], None, -1), + (NEW_FILTERS[0], None, 0), +]) +def test_numerics_singletons_channels_and_history(configured, singleton_axis, channel_axis): + base = np.random.default_rng(123).normal(size=(12, 13)) + def retain_singleton(values): + return values if singleton_axis is None else np.expand_dims(values, singleton_axis) + plane = retain_singleton(base) + expected = retain_singleton(reference(configured, base)) + metadata = {'case': {'ids': [12]}} + if channel_axis is not None: + image = np.stack([plane, plane * 2], axis=channel_axis) + expected = np.stack([expected, retain_singleton(reference(configured, base * 2))], axis=channel_axis) + metadata['channel_axis'] = channel_axis + else: + image = plane + before = image.copy() + result = configured.apply(image, metadata) + np.testing.assert_allclose(result.image, expected, rtol=1e-14, atol=1e-14) + np.testing.assert_array_equal(image, before) + assert result.image.shape == image.shape + assert result.image.dtype == (np.bool_ if isinstance(configured, CannyFilter) else np.float64) + assert np.isfinite(result.image).all() + entry = result.metadata['processing_history'][0] + assert entry['filter'] == type(configured).__name__ + assert entry['spatial_dimension'] == 2 + assert entry['input_shape'] == entry['output_shape'] == list(image.shape) + assert entry['output_dtype'] == result.image.dtype.name + assert entry['parameters']['mode'] == configured.mode + assert result.metadata['case'] == metadata['case'] + result.metadata['case']['ids'].append(42) + assert metadata['case']['ids'] == [12] + assert 'processing_history' not in metadata + with pytest.raises(TypeError): + result.metadata['case'] = {} + + +@pytest.mark.parametrize('configured,axis', list(zip(NEW_FILTERS[:4], [0, 1, 2, 3]))) +def test_multichannel_volumes_numerically(configured, axis): + volume = np.random.default_rng(40).normal(size=(6, 7, 8)) + image = np.stack([volume, volume * 3], axis=axis) + actual = configured.apply(image, {'channel_axis': axis}).image + expected = np.stack([reference(configured, volume), reference(configured, volume * 3)], axis=axis) + np.testing.assert_allclose(actual, expected, rtol=1e-14, atol=1e-14) + assert actual.shape == image.shape + assert actual.dtype == np.float64 + assert np.isfinite(actual).all() + + +@pytest.mark.parametrize('factory', [ + pytest.param(lambda: GaussianFilter(sigma=[1, 2]), id='immutable-sigma'), + pytest.param(lambda: GaussianFilter(sigma=0), id='positive-sigma'), + pytest.param(lambda: GaussianFilter(truncate=0), id='positive-truncate'), + pytest.param(lambda: MedianFilter(size=0), id='positive-window'), + pytest.param(lambda: MedianFilter(size=4, origin=2), id='origin-boundary'), + pytest.param(lambda: SobelFilter(axis=True), id='boolean-is-not-axis'), + pytest.param(lambda: LaplacianFilter(cval=np.inf), id='finite-padding'), + pytest.param(lambda: LaplacianFilter(mode='invalid'), id='boundary-mode'), + pytest.param(lambda: GaussianFilter(output_dtype='int16'), id='floating-output'), + pytest.param(lambda: CannyFilter(low_threshold=-1), id='nonnegative-threshold'), + pytest.param(lambda: CannyFilter(low_threshold=1.1, use_quantiles=True), id='quantile-range'), + pytest.param(lambda: CannyFilter(low_threshold=.8, high_threshold=.2), id='threshold-order'), +]) + +def test_invalid_types_and_nonfinite_parameters_use_domain_errors(factory): + with pytest.raises(InvalidConfigurationError): + factory() + + +@pytest.mark.parametrize('configured', [GaussianFilter((1, 2, 3)), SobelFilter(axis=2), CannyFilter(low_threshold=.8)]) +def test_image_dependent_configuration_errors(configured): + with pytest.raises(InvalidConfigurationError): + configured.apply(np.ones((10, 11))) + + +@pytest.mark.parametrize('dtype,quantiles', [(np.uint8, False), (np.int16, True), (np.float32, False), (np.bool_, True)]) +def test_canny_dtype_normalization_matches_backend(dtype, quantiles): + image = np.random.default_rng(12).integers(0, 100, (24, 25)).astype(dtype) + config = CannyFilter(use_quantiles=quantiles) + actual = config.apply(image).image + np.testing.assert_array_equal(actual, canny(image, use_quantiles=quantiles)) + assert actual.dtype == np.bool_ + + +@pytest.mark.parametrize('axis', [1, -1]) +def test_canny_slicewise_reference_and_nested_provenance(axis): + volume = np.random.default_rng(11).normal(size=(10, 11, 12)) + with pytest.raises(UnsupportedDimensionError, match='SliceWise'): + NEW_FILTERS[-1].apply(volume) + result = SliceWise(NEW_FILTERS[-1], axis=axis).apply(volume, {'source': 'test'}) + expected = np.stack([reference(NEW_FILTERS[-1], np.take(volume, i, axis=axis)) for i in range(volume.shape[axis])], axis=axis) + np.testing.assert_array_equal(result.image, expected) + assert result.image.shape == volume.shape + assert result.image.dtype == np.bool_ + entry = result.metadata['processing_history'][0] + assert entry['parameters']['filter']['filter'] == 'CannyFilter' + assert entry['parameters']['axis'] == axis + assert entry['output_dtype'] == 'bool' + assert result.metadata['source'] == 'test' + + +def test_combined_repeated_pipeline_batch_reference(): + from image_filter import ImageItem + images = [np.random.default_rng(i).normal(size=shape) for i, shape in enumerate([(13, 14), (7, 8, 9)])] + steps = [*NEW_FILTERS[:4], GaussianFilter(.8)] + pipeline = Pipeline(steps) + results = pipeline.apply_batch([ImageItem(image, {'id': i}) for i, image in enumerate(images)]) + for i, (image, result) in enumerate(zip(images, results)): + expected = image + for step in steps: + expected = reference(step, expected) + np.testing.assert_allclose(result.image, expected, rtol=1e-12, atol=1e-12) + assert result.image.shape == image.shape + assert result.image.dtype == np.float64 + assert np.isfinite(result.image).all() + assert result.metadata['id'] == i + assert all(entry['spatial_dimension'] == image.ndim for entry in result.metadata['processing_history']) + assert [entry['filter'] for entry in result.metadata['processing_history']] == [type(step).__name__ for step in steps] + # Boolean output remains valid numeric input for a subsequent filter. + image = images[0] + result = Pipeline([NEW_FILTERS[-1], GaussianFilter(.8), NEW_FILTERS[-1]]).apply(image) + expected = reference(NEW_FILTERS[-1], reference(NEW_FILTERS[0], reference(NEW_FILTERS[-1], image))) + np.testing.assert_array_equal(result.image, expected) + assert len(result.metadata['processing_history']) == 3 + + +def test_exactly_ten_public_concrete_filters(): + import image_filter as api + names = [name for name in api.__all__ if name.endswith('Filter') and name != 'ImageFilter'] + assert set(names) == {'MeanFilter', 'LoGFilter', 'LawsFilter', 'GaborFilter', 'WaveletFilter', 'GaussianFilter', 'MedianFilter', 'SobelFilter', 'LaplacianFilter', 'CannyFilter'} + assert api.CannyFilter.capabilities.spatial_dimensions == frozenset({2}) + + +@pytest.mark.parametrize('dtype', [np.int64, np.uint64]) +def test_canny_declares_backend_integer_restriction(dtype): + from image_filter import InvalidImageError + image = np.ones((10, 11), dtype=dtype) + with pytest.raises(ValueError, match='64-bit'): + canny(image) + with pytest.raises(InvalidImageError, match='64-bit'): + CannyFilter().apply(image) + assert '32 bits' in CannyFilter.capabilities.accepted_dtypes diff --git a/tests/test_composite.py b/tests/test_composite.py new file mode 100644 index 0000000..853f78d --- /dev/null +++ b/tests/test_composite.py @@ -0,0 +1,51 @@ +"""Composite behavior: substitution, ordering, metadata, and adapters.""" +import json + +import numpy as np +import pytest + +from image_filter import ( + CannyFilter, GaussianFilter, ImageFilter, ImageItem, InvalidConfigurationError, + MeanFilter, Pipeline, SliceWise, SobelFilter, UnsupportedDimensionError, +) + + +def test_nested_pipeline_matches_explicit_execution_and_preserves_metadata(): + image = np.arange(120, dtype=float).reshape(10, 12) + metadata = {"patient": {"label": "test"}} + first, second, third = MeanFilter(3), GaussianFilter(0.8), SobelFilter() + operation = Pipeline([Pipeline([first, second]), third]) + assert isinstance(operation, ImageFilter) + result = operation.apply(image, metadata) + expected = third.apply(second.apply(first.apply(image).image).image) + np.testing.assert_array_equal(result.image, expected.image) + assert [x["filter"] for x in result.metadata["processing_history"]] == [ + "MeanFilter", "GaussianFilter", "SobelFilter", + ] + assert metadata == {"patient": {"label": "test"}} + assert result.metadata["patient"] is not metadata["patient"] + np.testing.assert_array_equal(image, np.arange(120).reshape(10, 12)) + for component in (first, operation): + batch = component.apply_batch([ImageItem(image), ImageItem(image[None])]) + np.testing.assert_array_equal(batch[0].image, batch[1].image[0]) + + +def test_pipeline_in_slicewise_preserves_boolean_output_and_serializable_history(): + volume = np.random.default_rng(12).random((10, 12, 5)) + pipeline = Pipeline([GaussianFilter(0.8), CannyFilter(use_quantiles=True)]) + result = SliceWise(pipeline, axis=2).apply(volume) + expected = np.stack([pipeline.apply(volume[:, :, i]).image for i in range(5)], axis=2) + np.testing.assert_array_equal(result.image, expected) + assert result.image.dtype == bool + json.dumps(dict(result.metadata)) + assert pipeline.capabilities.spatial_dimensions == frozenset({2}) + with pytest.raises(UnsupportedDimensionError): + pipeline.apply(volume) + + +def test_invalid_composites_fail_at_construction(): + for steps in ([], [object()], [Pipeline([MeanFilter()]), "invalid"]): + with pytest.raises(InvalidConfigurationError): + Pipeline(steps) + with pytest.raises(InvalidConfigurationError): + SliceWise(object(), axis=0) diff --git a/tests/test_dimensions.py b/tests/test_dimensions.py new file mode 100644 index 0000000..fecad9c --- /dev/null +++ b/tests/test_dimensions.py @@ -0,0 +1,53 @@ +from __future__ import annotations + +import numpy as np +import pytest + +from image_filter import ( + DimensionInferenceError, + GaborFilter, + SliceWise, + UnsupportedDimensionError, + infer_spatial_layout, +) + + +def test_short_volume_or_rgb_without_metadata_is_ambiguous(): + with pytest.raises(DimensionInferenceError, match="ambiguous"): + infer_spatial_layout(np.zeros((16, 17, 3))) + assert infer_spatial_layout(np.zeros((16, 17, 3)), {"spatial_dims": 3}).dimension == 3 + + +def test_4d_requires_channel_axis(): + image = np.zeros((8, 9, 10, 2)) + with pytest.raises(DimensionInferenceError, match="requires"): + infer_spatial_layout(image) + assert infer_spatial_layout(image, {"channel_axis": 3}).dimension == 3 + + +def test_intrinsically_2d_filter_rejects_volume_in_strict_mode(): + with pytest.raises(UnsupportedDimensionError, match="SliceWise"): + GaborFilter().apply(np.ones((6, 7, 8)), {"spatial_dims": 3}) + + +def test_slice_wise_is_explicit_and_axis_sensitive(): + image = np.zeros((4, 12, 13), dtype=np.float32) + image[:, 5, :] = 1 + adapter = SliceWise(GaborFilter(frequency=0.2), axis=0) + result = adapter.apply(image) + from skimage.filters import gabor + expected = np.stack([np.hypot(*gabor(plane.astype(np.float64), frequency=.2, sigma_x=2, sigma_y=2)) for plane in image]) + np.testing.assert_allclose(result.image, expected) + assert result.image.shape == image.shape + assert result.image.dtype == np.float64 + assert np.isfinite(result.image).all() + assert result.metadata["processing_history"][-1]["filter"] == "SliceWise" + nested = result.metadata["processing_history"][-1]["parameters"]["filter"] + assert nested["filter"] == "GaborFilter" + assert nested["parameters"]["frequency"] == 0.2 + + +def test_slice_wise_rejects_non_spatial_axis(): + image = np.ones((5, 6, 7, 3)) + with pytest.raises(UnsupportedDimensionError, match="spatial axes"): + SliceWise(GaborFilter(), axis=3).apply(image, {"channel_axis": 3}) diff --git a/tests/test_pipeline_and_contracts.py b/tests/test_pipeline_and_contracts.py new file mode 100644 index 0000000..459b4f2 --- /dev/null +++ b/tests/test_pipeline_and_contracts.py @@ -0,0 +1,11 @@ +"""Unsupported features must fail explicitly instead of silently approximating.""" +import pytest + +from image_filter import LawsFilter, UnsupportedOptionError, WaveletFilter + + +def test_unsupported_advanced_options_fail_at_configuration_time(): + with pytest.raises(UnsupportedOptionError): + LawsFilter(rotation_invariant=True) + with pytest.raises(UnsupportedOptionError): + WaveletFilter(rotation_invariant=True) diff --git a/tests/test_real_data.py b/tests/test_real_data.py new file mode 100644 index 0000000..d41fd47 --- /dev/null +++ b/tests/test_real_data.py @@ -0,0 +1,112 @@ +"""Integration tests using real repository scans and scikit-image's camera. + +``skimage.data.camera`` is a public-domain photograph bundled/cached by the +installed scientific dependency; it supplies a genuine 2D image without +network access. NIfTI inputs below are existing Radiuma repository assets and +are used only as image data, never as an implementation reference. +""" + +from __future__ import annotations + +from pathlib import Path + +import nibabel as nib +import numpy as np +import pytest +from scipy import ndimage as ndi +from skimage import data +from skimage.filters import gabor + +from image_filter import ( + CannyFilter, + GaborFilter, + LoGFilter, + MeanFilter, + Pipeline, + SliceWise, +) + +REPOSITORY = Path(__file__).resolve().parents[1] / 'data/tests' + + +@pytest.mark.real_data +def test_public_api_on_real_2d_photograph(): + image = data.camera() + result = Pipeline([ + MeanFilter(3, output_dtype="float32"), + GaborFilter(frequency=0.16, orientations=(0.0, np.pi / 2), output_dtype="float32"), + ]).apply(image, {"source": "skimage.data.camera (public domain)"}) + assert result.image.shape == (512, 512) + assert result.image.dtype == np.float32 + assert np.isfinite(result.image).all() + smoothed = ndi.uniform_filter(image, 3, output=np.float32) + expected = np.maximum(*[ + np.hypot(*gabor(smoothed, frequency=.16, theta=theta, sigma_x=2, sigma_y=2)) + for theta in (0, np.pi / 2) + ]) + np.testing.assert_allclose(result.image, expected) + assert result.metadata['processing_history'][1]['filter'] == 'GaborFilter' + + +@pytest.mark.real_data +@pytest.mark.parametrize( + "relative_path", + ["test4.nii.gz", "0001_MRI_Image_a.nii.gz"], +) +def test_real_nifti_volume_end_to_end(relative_path): + nii = nib.load(REPOSITORY / relative_path) + full = np.asanyarray(nii.dataobj) + widths = (12, 24, 28) + peak = np.unravel_index(np.argmax(np.abs(full)), full.shape) + slices = tuple( + slice(start := min(max(0, center - width // 2), size - width), start + width) + for size, width, center in zip(full.shape, widths, peak, strict=True) + ) + source = full[slices] + metadata = {"spacing": tuple(float(v) for v in nii.header.get_zooms()[:3]), "spatial_dims": 3} + result = Pipeline([ + MeanFilter(3, mode="nearest", output_dtype="float32"), + LoGFilter(0.9, output_dtype="float32"), + MeanFilter(4, origin=-1, output_dtype="float32"), + ]).apply(source, metadata) + assert result.image.shape == source.shape + assert result.image.dtype == np.float32 + assert np.isfinite(result.image).all() + assert not np.shares_memory(result.image, source) + assert len(result.metadata["processing_history"]) == 3 + expected = ndi.uniform_filter(source, 3, mode="nearest", output=np.float32) + expected = ndi.gaussian_laplace(expected, .9) + expected = ndi.uniform_filter(expected, 4, origin=-1, output=np.float32) + np.testing.assert_allclose(result.image, expected, rtol=1e-6, atol=1e-6) + assert result.metadata['spacing'] == metadata['spacing'] + + +@pytest.mark.real_data +@pytest.mark.parametrize('source', ['camera', 'test4.nii.gz', '0001_MRI_Image_a.nii.gz']) +def test_all_new_filters_real_backend_numerics(source): + from test_additional_filters import NEW_FILTERS, reference + if source == 'camera': + image = data.camera() + else: + nii = nib.load(REPOSITORY / source) + full = np.asanyarray(nii.dataobj) + peak = np.unravel_index(np.argmax(np.abs(full)), full.shape) + region = tuple(slice(max(0, min(c - 8, n - 16)), max(0, min(c - 8, n - 16)) + 16) for c, n in zip(peak, full.shape)) + image = full[region] + metadata = {'source': source, 'spatial_dims': image.ndim} + for configured in NEW_FILTERS: + if isinstance(configured, CannyFilter) and image.ndim == 3: + applied = SliceWise(configured, axis=0) + expected = np.stack([reference(configured, plane) for plane in image]) + else: + applied = configured + expected = reference(configured, image) + result = applied.apply(image, metadata) + np.testing.assert_allclose(result.image, expected, rtol=1e-13, atol=1e-10) + assert result.image.shape == image.shape + assert result.image.dtype == expected.dtype + assert np.isfinite(result.image).all() + assert result.metadata['source'] == source + assert result.metadata['processing_history'][0]['filter'] == type(applied).__name__ + assert result.metadata['processing_history'][0]['spatial_dimension'] == image.ndim + assert 'processing_history' not in metadata diff --git a/tests/test_reference_numerics.py b/tests/test_reference_numerics.py new file mode 100644 index 0000000..81e946e --- /dev/null +++ b/tests/test_reference_numerics.py @@ -0,0 +1,95 @@ +from __future__ import annotations + +import numpy as np +import pywt +import pytest +from scipy import ndimage as ndi +from skimage.filters import gabor + +from image_filter import GaborFilter, LawsFilter, LoGFilter, MeanFilter, WaveletFilter + + +@pytest.mark.parametrize("size,mode", [(3, "reflect"), (4, "nearest"), (3, "mirror"), (4, "wrap"), (4, "constant")]) +def test_mean_matches_scipy_for_even_odd_and_boundaries(size, mode): + image = np.arange(42, dtype=np.int16).reshape(6, 7) + actual = MeanFilter(size=size, mode=mode, cval=-2, output_dtype="float32").apply(image).image + expected = ndi.uniform_filter(image, size=size, mode=mode, cval=-2, output=np.float32) + np.testing.assert_allclose(actual, expected, rtol=0, atol=0) + assert actual.dtype == np.float32 + + +def test_mean_origin_defines_even_kernel_alignment(): + image = np.arange(30, dtype=np.float64).reshape(5, 6) + actual = MeanFilter(size=4, origin=-1).apply(image).image + expected = ndi.uniform_filter(image, size=4, origin=-1, output=np.float64) + np.testing.assert_array_equal(actual, expected) + assert not np.array_equal(actual, MeanFilter(size=4, origin=0).apply(image).image) + + +def test_log_matches_scipy_gaussian_laplace(): + image = np.arange(72, dtype=np.float64).reshape(8, 9) + actual = LoGFilter(sigma=(1.2, 0.8), truncate=3, mode="mirror").apply(image).image + expected = ndi.gaussian_laplace(image, sigma=(1.2, 0.8), truncate=3, mode="mirror") + np.testing.assert_allclose(actual, expected, rtol=2e-14, atol=2e-14) + + +def test_log_rgb_matches_per_channel_scipy_reference(): + image = np.arange(8 * 9 * 3, dtype=np.float64).reshape(8, 9, 3) + actual = LoGFilter(sigma=1.1).apply(image, {"channel_axis": 2}).image + expected = np.stack([ndi.gaussian_laplace(image[..., c], 1.1) for c in range(3)], axis=2) + np.testing.assert_allclose(actual, expected, rtol=2e-14, atol=2e-14) + + +def test_laws_matches_separable_scipy_reference(): + image = np.arange(90, dtype=np.float64).reshape(9, 10) + l5 = np.array([1, 4, 6, 4, 1], dtype=float); l5 /= np.linalg.norm(l5) + e5 = np.array([-1, -2, 0, 2, 1], dtype=float); e5 /= np.linalg.norm(e5) + expected = ndi.convolve1d(image, l5, axis=0, mode="reflect") + expected = ndi.convolve1d(expected, e5, axis=1, mode="reflect") + actual = LawsFilter("L5E5").apply(image).image + np.testing.assert_allclose(actual, expected, rtol=1e-14, atol=1e-14) + + +def test_laws_energy_is_local_mean_absolute_response(): + image = np.arange(90, dtype=np.float64).reshape(9, 10) + response = LawsFilter("L5E5").apply(image).image + expected = ndi.uniform_filter(np.abs(response), size=3, mode="reflect") + actual = LawsFilter("L5E5", energy_radius=1).apply(image).image + np.testing.assert_allclose(actual, expected) + assert np.all(actual >= 0) + + +def test_gabor_matches_skimage_reference(): + y, x = np.mgrid[:32, :35] + image = np.sin(2 * np.pi * 0.2 * x) + actual = GaborFilter( + frequency=0.2, orientations=(0.3,), sigma=2.1, gamma=0.7, + response="magnitude", mode="reflect", + ).apply(image).image + real, imaginary = gabor( + image, frequency=0.2, theta=0.3, sigma_x=2.1, + sigma_y=2.1 / 0.7, mode="reflect", cval=0, + ) + np.testing.assert_allclose(actual, np.hypot(real, imaginary)) + + +def test_wavelet_matches_pywavelets_stationary_transform(): + image = np.arange(8 * 10, dtype=np.float64).reshape(8, 10) + padded = np.pad(image, ((0, 0), (0, 2)), mode="reflect") + expected = pywt.swtn(padded, "db2", level=2, axes=(0, 1), trim_approx=False)[0]["ad"] + actual = WaveletFilter("LH", wavelet="db2", level=2).apply(image).image + np.testing.assert_allclose(actual, expected[:, :10]) + assert actual.shape == image.shape + + + + +def test_mean_excludes_channel_and_singleton_axes(): + # Mean has a separate full-array implementation from the newer filters. + base = np.random.default_rng(14).normal(size=(8, 9)) + image = np.stack([base, base * 2 + 20], axis=0)[:, :, None, :] + result = MeanFilter(4).apply(image, {"channel_axis": 0}) + expected = np.stack([ndi.uniform_filter(plane[:, 0, :], 4) for plane in image])[:, :, None, :] + np.testing.assert_allclose(result.image, expected) + assert result.image.shape == image.shape + assert result.metadata['processing_history'][0]['spatial_axes'] == [1, 3] diff --git a/tests/test_usage_examples.py b/tests/test_usage_examples.py new file mode 100644 index 0000000..ea5bb60 --- /dev/null +++ b/tests/test_usage_examples.py @@ -0,0 +1,47 @@ +"""Exercise the common runner with real MRI data and different components.""" +import importlib.util +import json +from pathlib import Path +import sys + +import numpy as np +import pytest + +EXAMPLES = Path(__file__).resolve().parents[1] / "examples" +sys.path.insert(0, str(EXAMPLES)) +from example_support import run_example +from image_filter import UnsupportedDimensionError + + +def operation(filename): + spec = importlib.util.spec_from_file_location("usage_example", EXAMPLES / filename) + module = importlib.util.module_from_spec(spec) + spec.loader.exec_module(module) + return module.build_operation() + + +def test_all_examples_share_runner_and_leaf_matches_single_step_pipeline(tmp_path): + outputs = [] + for script in sorted(EXAMPLES.glob("0[1-5]_*.py")): + run_example(operation(script.name), script.stem, ["--output-dir", str(tmp_path)]) + outputs.append(np.load(tmp_path / f"{script.stem}.npy")) + metadata = json.loads((tmp_path / f"{script.stem}.json").read_text()) + assert metadata["processing_history"] + assert outputs[-1].ndim == 2 + np.testing.assert_array_equal(outputs[0], outputs[3]) + + +def test_volume_policy_and_user_file_input(tmp_path): + from medical_io import bundled_cases + path, _ = bundled_cases()["mri"] + args = ["--input", str(path), "--output-dir", str(tmp_path)] + run_example(operation("01_single_filter.py"), "volume", args + ["--mode", "volume"]) + volume = np.load(tmp_path / "volume.npy") + assert volume.ndim == 3 + edges = operation("05_slicewise_pipeline.py") + run_example(edges, "slices", args + ["--mode", "slice-wise"]) + result = np.load(tmp_path / "slices.npy") + assert result.shape == volume.shape and result.dtype == bool + with pytest.raises(UnsupportedDimensionError): + run_example(edges, "invalid", args + ["--mode", "volume"]) + assert not (tmp_path / "invalid.npy").exists() diff --git a/tests/test_visual_compare.py b/tests/test_visual_compare.py new file mode 100644 index 0000000..5fb3272 --- /dev/null +++ b/tests/test_visual_compare.py @@ -0,0 +1,46 @@ +"""Small CLI contract and real rendering checks; numerical filters are tested elsewhere.""" +import importlib.util +from pathlib import Path + +import numpy as np +from PIL import Image +import pytest + +SCRIPT = Path(__file__).resolve().parents[1] / 'examples' / 'visual_compare.py' +spec = importlib.util.spec_from_file_location('visual_compare', SCRIPT) +visual = importlib.util.module_from_spec(spec) +spec.loader.exec_module(visual) + + +def test_cli_selection_and_destination_contract(): + assert visual.parse_args([]).names == ['mri'] + assert visual.parse_args(['--all']).names == list(visual.CASE_NAMES) + args = visual.parse_args(['--images', 'mri', 'pet', '--output-dir', 'panels']) + assert args.names == ['mri', 'pet'] and args.output_dir == Path('panels') + with pytest.raises(SystemExit): + visual.parse_args(['--all', '--output', 'one.png']) + with pytest.raises(SystemExit): + visual.parse_args(['--images', 'unknown']) + + +def test_real_multi_image_rendering_offline(tmp_path, monkeypatch): + # Fail if any sample loader attempts to use networking. + import socket + def forbid_network(*args, **kwargs): + raise AssertionError('Example generation must stay offline') + monkeypatch.setattr(socket.socket, 'connect', forbid_network) + visual.main(['--images', 'mri', 'mr-dicom', '--output-dir', str(tmp_path), '--contact-sheet']) + expected = {'mri-filters.png': (2400, 1800), 'mr-dicom-filters.png': (2400, 1800), + 'source-overview.png': (1800, 750), 'cases.json': None} + assert {p.name for p in tmp_path.iterdir()} == set(expected) + for name, size in expected.items(): + if name == 'cases.json': + assert 'mri' in (tmp_path / name).read_text() + continue + with Image.open(tmp_path / name) as image: + image.load() + assert image.format == 'PNG' and image.size == size and image.mode == 'RGBA' + pixels = np.asarray(image.convert('RGB')) + assert pixels.std() > 20 + # Image area must contain actual content, beyond just titles. + assert pixels[size[1] // 3:2 * size[1] // 3].std() > 20