diff --git a/README.md b/README.md index 68d27b0f..ff545638 100644 --- a/README.md +++ b/README.md @@ -171,6 +171,17 @@ In the GIS world, rasters are used for representing continuous phenomena (e.g. e ------- +### **Morphological** + +| Name | Description | NumPy xr.DataArray | Dask xr.DataArray | CuPy GPU xr.DataArray | Dask GPU xr.DataArray | +|:----------:|:------------|:----------------------:|:--------------------:|:-------------------:|:------:| +| [Erode](xrspatial/morphology.py) | Morphological erosion (local minimum over structuring element) | ✅️ | ✅️ | ✅️ | ✅️ | +| [Dilate](xrspatial/morphology.py) | Morphological dilation (local maximum over structuring element) | ✅️ | ✅️ | ✅️ | ✅️ | +| [Opening](xrspatial/morphology.py) | Erosion then dilation (removes small bright features) | ✅️ | ✅️ | ✅️ | ✅️ | +| [Closing](xrspatial/morphology.py) | Dilation then erosion (fills small dark gaps) | ✅️ | ✅️ | ✅️ | ✅️ | + +------- + ### **Fire** | Name | Description | NumPy xr.DataArray | Dask xr.DataArray | CuPy GPU xr.DataArray | Dask GPU xr.DataArray | diff --git a/docs/source/reference/index.rst b/docs/source/reference/index.rst index 5f8c36e3..9de20872 100644 --- a/docs/source/reference/index.rst +++ b/docs/source/reference/index.rst @@ -15,6 +15,7 @@ Reference focal hydrology interpolation + morphology multispectral pathfinding proximity diff --git a/docs/source/reference/morphology.rst b/docs/source/reference/morphology.rst new file mode 100644 index 00000000..c6625a15 --- /dev/null +++ b/docs/source/reference/morphology.rst @@ -0,0 +1,40 @@ +.. _reference.morphology: + +********** +Morphology +********** + +Erode +===== +.. autosummary:: + :toctree: _autosummary + + xrspatial.morphology.morph_erode + +Dilate +====== +.. autosummary:: + :toctree: _autosummary + + xrspatial.morphology.morph_dilate + +Opening +======= +.. autosummary:: + :toctree: _autosummary + + xrspatial.morphology.morph_opening + +Closing +======= +.. autosummary:: + :toctree: _autosummary + + xrspatial.morphology.morph_closing + +Kernel Construction +=================== +.. autosummary:: + :toctree: _autosummary + + xrspatial.morphology._circle_kernel diff --git a/examples/user_guide/17_Morphological_Operators.ipynb b/examples/user_guide/17_Morphological_Operators.ipynb new file mode 100644 index 00000000..d1793ce8 --- /dev/null +++ b/examples/user_guide/17_Morphological_Operators.ipynb @@ -0,0 +1,217 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Morphological Operators\n", + "\n", + "This notebook demonstrates the four morphological raster operators\n", + "in `xrspatial.morphology`:\n", + "\n", + "- **`morph_erode`** -- local minimum (shrinks bright regions)\n", + "- **`morph_dilate`** -- local maximum (expands bright regions)\n", + "- **`morph_opening`** -- erosion then dilation (removes small bright features)\n", + "- **`morph_closing`** -- dilation then erosion (fills small dark gaps)\n", + "\n", + "These operations work on any numeric raster, including grayscale\n", + "elevation data and binary classification masks. All four backends\n", + "(numpy, cupy, dask+numpy, dask+cupy) are supported." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import xarray as xr\n", + "import matplotlib.pyplot as plt\n", + "\n", + "from xrspatial.morphology import (\n", + " _circle_kernel,\n", + " morph_closing,\n", + " morph_dilate,\n", + " morph_erode,\n", + " morph_opening,\n", + ")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 1. Synthetic elevation data\n", + "\n", + "Create a small raster with a peak and a pit to show how erosion and\n", + "dilation behave." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "rng = np.random.default_rng(42)\n", + "y, x = np.mgrid[-3:3:0.05, -3:3:0.05]\n", + "z = np.exp(-(x**2 + y**2)) - 0.5 * np.exp(-((x - 1.5)**2 + (y - 1)**2) / 0.3)\n", + "z += 0.05 * rng.standard_normal(z.shape) # add noise\n", + "\n", + "raster = xr.DataArray(z, dims=['y', 'x'], name='elevation')\n", + "\n", + "fig, ax = plt.subplots(figsize=(6, 5))\n", + "raster.plot(ax=ax, cmap='terrain')\n", + "ax.set_title('Original raster')\n", + "plt.tight_layout()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 2. Erosion and dilation" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "kernel = np.ones((5, 5), dtype=np.uint8)\n", + "\n", + "eroded = morph_erode(raster, kernel=kernel, boundary='nearest')\n", + "dilated = morph_dilate(raster, kernel=kernel, boundary='nearest')\n", + "\n", + "fig, axes = plt.subplots(1, 3, figsize=(15, 4))\n", + "for ax, data, title in zip(\n", + " axes,\n", + " [raster, eroded, dilated],\n", + " ['Original', 'Eroded (local min)', 'Dilated (local max)'],\n", + "):\n", + " data.plot(ax=ax, cmap='terrain', vmin=z.min(), vmax=z.max())\n", + " ax.set_title(title)\n", + "plt.tight_layout()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 3. Opening and closing\n", + "\n", + "Opening removes small bright spikes; closing fills small dark pits." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "opened = morph_opening(raster, kernel=kernel, boundary='nearest')\n", + "closed = morph_closing(raster, kernel=kernel, boundary='nearest')\n", + "\n", + "fig, axes = plt.subplots(1, 3, figsize=(15, 4))\n", + "for ax, data, title in zip(\n", + " axes,\n", + " [raster, opened, closed],\n", + " ['Original', 'Opening', 'Closing'],\n", + "):\n", + " data.plot(ax=ax, cmap='terrain', vmin=z.min(), vmax=z.max())\n", + " ax.set_title(title)\n", + "plt.tight_layout()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 4. Binary mask cleanup\n", + "\n", + "A common use case: clean up noisy classification results. Opening\n", + "removes salt noise; closing removes pepper noise." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Generate a binary mask with noise\n", + "mask = np.zeros((100, 100), dtype=np.float64)\n", + "mask[20:80, 20:80] = 1.0 # large square\n", + "\n", + "# Sprinkle salt-and-pepper noise\n", + "noise = rng.random(mask.shape)\n", + "mask[noise < 0.02] = 1.0 # salt\n", + "mask[noise > 0.98] = 0.0 # pepper\n", + "\n", + "noisy = xr.DataArray(mask, dims=['y', 'x'], name='mask')\n", + "cleaned = morph_opening(\n", + " morph_closing(noisy, kernel=kernel, boundary='nearest'),\n", + " kernel=kernel,\n", + " boundary='nearest',\n", + ")\n", + "\n", + "fig, axes = plt.subplots(1, 2, figsize=(10, 4))\n", + "noisy.plot(ax=axes[0], cmap='gray')\n", + "axes[0].set_title('Noisy mask')\n", + "cleaned.plot(ax=axes[1], cmap='gray')\n", + "axes[1].set_title('After closing + opening')\n", + "plt.tight_layout()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 5. Circular structuring element\n", + "\n", + "Use `_circle_kernel(radius)` to create a disk-shaped kernel\n", + "instead of a square." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "disk = _circle_kernel(3)\n", + "print('Circular kernel (radius=3):')\n", + "print(disk)\n", + "\n", + "eroded_disk = morph_erode(raster, kernel=disk, boundary='nearest')\n", + "\n", + "fig, axes = plt.subplots(1, 2, figsize=(10, 4))\n", + "raster.plot(ax=axes[0], cmap='terrain')\n", + "axes[0].set_title('Original')\n", + "eroded_disk.plot(ax=axes[1], cmap='terrain')\n", + "axes[1].set_title('Eroded with circular kernel (r=3)')\n", + "plt.tight_layout()\n", + "plt.show()" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "name": "python", + "version": "3.11.0" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/xrspatial/__init__.py b/xrspatial/__init__.py index 7f2ba46e..38b3b823 100644 --- a/xrspatial/__init__.py +++ b/xrspatial/__init__.py @@ -41,6 +41,10 @@ from xrspatial.flow_length import flow_length # noqa from xrspatial.flow_path import flow_path # noqa from xrspatial.focal import mean # noqa +from xrspatial.morphology import morph_closing # noqa +from xrspatial.morphology import morph_dilate # noqa +from xrspatial.morphology import morph_erode # noqa +from xrspatial.morphology import morph_opening # noqa from xrspatial.hand import hand # noqa from xrspatial.hillshade import hillshade # noqa from xrspatial.mahalanobis import mahalanobis # noqa diff --git a/xrspatial/accessor.py b/xrspatial/accessor.py index 837b3864..a97f6886 100644 --- a/xrspatial/accessor.py +++ b/xrspatial/accessor.py @@ -173,6 +173,24 @@ def focal_mean(self, **kwargs): from .focal import mean return mean(self._obj, **kwargs) + # ---- Morphological ---- + + def morph_erode(self, **kwargs): + from .morphology import morph_erode + return morph_erode(self._obj, **kwargs) + + def morph_dilate(self, **kwargs): + from .morphology import morph_dilate + return morph_dilate(self._obj, **kwargs) + + def morph_opening(self, **kwargs): + from .morphology import morph_opening + return morph_opening(self._obj, **kwargs) + + def morph_closing(self, **kwargs): + from .morphology import morph_closing + return morph_closing(self._obj, **kwargs) + # ---- Proximity / Distance ---- def proximity(self, **kwargs): @@ -501,6 +519,24 @@ def focal_mean(self, **kwargs): from .focal import mean return mean(self._obj, **kwargs) + # ---- Morphological ---- + + def morph_erode(self, **kwargs): + from .morphology import morph_erode + return morph_erode(self._obj, **kwargs) + + def morph_dilate(self, **kwargs): + from .morphology import morph_dilate + return morph_dilate(self._obj, **kwargs) + + def morph_opening(self, **kwargs): + from .morphology import morph_opening + return morph_opening(self._obj, **kwargs) + + def morph_closing(self, **kwargs): + from .morphology import morph_closing + return morph_closing(self._obj, **kwargs) + # ---- Diffusion ---- def diffuse(self, **kwargs): diff --git a/xrspatial/morphology.py b/xrspatial/morphology.py new file mode 100644 index 00000000..4aba31b5 --- /dev/null +++ b/xrspatial/morphology.py @@ -0,0 +1,562 @@ +"""Morphological raster operators: erode, dilate, opening, closing. + +Applies grayscale morphological operations using a structuring element +(kernel) over a 2D raster. Erosion computes the local minimum and +dilation the local maximum within the kernel footprint. Opening +(erosion then dilation) removes small bright features; closing +(dilation then erosion) fills small dark gaps. + +Supports all four backends: numpy, cupy, dask+numpy, dask+cupy. +""" + +from __future__ import annotations + +from functools import partial + +import numpy as np +import xarray as xr +from xarray import DataArray + +from numba import cuda, prange + +try: + import cupy +except ImportError: + class cupy: + ndarray = False + +try: + import dask.array as da +except ImportError: + da = None + +from xrspatial.utils import ( + ArrayTypeFunctionMapping, + _boundary_to_dask, + _pad_array, + _validate_boundary, + _validate_raster, + calc_cuda_dims, + has_cuda_and_cupy, + ngjit, + not_implemented_func, +) +from xrspatial.dataset_support import supports_dataset + + +# --------------------------------------------------------------------------- +# Kernel construction +# --------------------------------------------------------------------------- + +def _circle_kernel(radius): + """Return a boolean 2D array for a circular structuring element. + + The kernel has shape ``(2*radius+1, 2*radius+1)`` with ``True`` + for cells whose centre is within *radius* cells of the centre. + """ + size = 2 * radius + 1 + y, x = np.ogrid[-radius:radius + 1, -radius:radius + 1] + return (x * x + y * y <= radius * radius).astype(np.uint8) + + +def _validate_kernel(kernel, func_name): + """Ensure *kernel* is a 2D uint8 numpy array with odd dimensions.""" + if not isinstance(kernel, np.ndarray): + raise TypeError( + f"{func_name}(): `kernel` must be a numpy ndarray, " + f"got {type(kernel).__name__}" + ) + if kernel.ndim != 2: + raise ValueError( + f"{func_name}(): `kernel` must be 2D, got {kernel.ndim}D" + ) + if kernel.shape[0] % 2 == 0 or kernel.shape[1] % 2 == 0: + raise ValueError( + f"{func_name}(): `kernel` dimensions must be odd, " + f"got {kernel.shape}" + ) + return kernel.astype(np.uint8) + + +# --------------------------------------------------------------------------- +# numpy backend +# --------------------------------------------------------------------------- + +@ngjit +def _erode_kernel_numpy(data, kernel, rows, cols, ky, kx): + """Erosion (local minimum) on a padded array.""" + out = np.empty((rows, cols), dtype=data.dtype) + hy = ky // 2 + hx = kx // 2 + for i in prange(rows): + for j in range(cols): + val = data[i + hy, j + hx] + if val != val: # NaN + out[i, j] = val + continue + mn = val + for dy in range(ky): + for dx in range(kx): + if kernel[dy, dx] == 0: + continue + v = data[i + dy, j + dx] + if v != v: # NaN neighbour + mn = v + break + if v < mn: + mn = v + if mn != mn: # propagate NaN + break + out[i, j] = mn + return out + + +@ngjit +def _dilate_kernel_numpy(data, kernel, rows, cols, ky, kx): + """Dilation (local maximum) on a padded array.""" + out = np.empty((rows, cols), dtype=data.dtype) + hy = ky // 2 + hx = kx // 2 + for i in prange(rows): + for j in range(cols): + val = data[i + hy, j + hx] + if val != val: # NaN + out[i, j] = val + continue + mx = val + for dy in range(ky): + for dx in range(kx): + if kernel[dy, dx] == 0: + continue + v = data[i + dy, j + dx] + if v != v: # NaN neighbour + mx = v + break + if v > mx: + mx = v + if mx != mx: + break + out[i, j] = mx + return out + + +def _morph_numpy(data, kernel, op, boundary): + """Run erosion or dilation on a numpy array.""" + ky, kx = kernel.shape + hy, hx = ky // 2, kx // 2 + rows, cols = data.shape + arr = data.astype(np.float64) + if boundary == 'nan': + padded = np.pad(arr, ((hy, hy), (hx, hx)), + mode='constant', constant_values=np.nan) + else: + padded = _pad_array(arr, (hy, hx), boundary) + if op == 'erode': + return _erode_kernel_numpy(padded, kernel, rows, cols, ky, kx) + else: + return _dilate_kernel_numpy(padded, kernel, rows, cols, ky, kx) + + +def _erode_numpy(data, kernel, boundary): + return _morph_numpy(data, kernel, 'erode', boundary) + + +def _dilate_numpy(data, kernel, boundary): + return _morph_numpy(data, kernel, 'dilate', boundary) + + +def _opening_numpy(data, kernel, boundary): + tmp = _morph_numpy(data, kernel, 'erode', boundary) + return _morph_numpy(tmp, kernel, 'dilate', boundary) + + +def _closing_numpy(data, kernel, boundary): + tmp = _morph_numpy(data, kernel, 'dilate', boundary) + return _morph_numpy(tmp, kernel, 'erode', boundary) + + +# --------------------------------------------------------------------------- +# cupy backend +# --------------------------------------------------------------------------- + +@cuda.jit +def _erode_gpu(data, kernel, out, hy, hx, ky, kx): + i, j = cuda.grid(2) + rows = out.shape[0] + cols = out.shape[1] + if i < rows and j < cols: + val = data[i + hy, j + hx] + if val != val: + out[i, j] = val + return + mn = val + for dy in range(ky): + for dx in range(kx): + if kernel[dy, dx] == 0: + continue + v = data[i + dy, j + dx] + if v != v: + mn = v + break + if v < mn: + mn = v + if mn != mn: + break + out[i, j] = mn + + +@cuda.jit +def _dilate_gpu(data, kernel, out, hy, hx, ky, kx): + i, j = cuda.grid(2) + rows = out.shape[0] + cols = out.shape[1] + if i < rows and j < cols: + val = data[i + hy, j + hx] + if val != val: + out[i, j] = val + return + mx = val + for dy in range(ky): + for dx in range(kx): + if kernel[dy, dx] == 0: + continue + v = data[i + dy, j + dx] + if v != v: + mx = v + break + if v > mx: + mx = v + if mx != mx: + break + out[i, j] = mx + + +def _morph_cupy(data, kernel, op, boundary): + import cupy as cp + + ky, kx = kernel.shape + hy, hx = ky // 2, kx // 2 + rows, cols = data.shape + arr = cp.asarray(data, dtype=cp.float64) + kern_dev = cp.asarray(kernel) + if boundary == 'nan': + padded = cp.pad(arr, ((hy, hy), (hx, hx)), + mode='constant', constant_values=cp.nan) + else: + padded = _pad_array(arr, (hy, hx), boundary) + out = cp.empty((rows, cols), dtype=cp.float64) + bpg, tpb = calc_cuda_dims((rows, cols)) + if op == 'erode': + _erode_gpu[bpg, tpb](padded, kern_dev, out, hy, hx, ky, kx) + else: + _dilate_gpu[bpg, tpb](padded, kern_dev, out, hy, hx, ky, kx) + return out + + +def _erode_cupy(data, kernel, boundary): + return _morph_cupy(data, kernel, 'erode', boundary) + + +def _dilate_cupy(data, kernel, boundary): + return _morph_cupy(data, kernel, 'dilate', boundary) + + +def _opening_cupy(data, kernel, boundary): + tmp = _morph_cupy(data, kernel, 'erode', boundary) + return _morph_cupy(tmp, kernel, 'dilate', boundary) + + +def _closing_cupy(data, kernel, boundary): + tmp = _morph_cupy(data, kernel, 'dilate', boundary) + return _morph_cupy(tmp, kernel, 'erode', boundary) + + +# --------------------------------------------------------------------------- +# dask + numpy backend +# --------------------------------------------------------------------------- + +def _morph_chunk_numpy(chunk, kernel, op, block_info=None): + """Process a single overlapped dask chunk.""" + ky, kx = kernel.shape + hy, hx = ky // 2, kx // 2 + rows = chunk.shape[0] - 2 * hy + cols = chunk.shape[1] - 2 * hx + if op == 'erode': + interior = _erode_kernel_numpy(chunk, kernel, rows, cols, ky, kx) + else: + interior = _dilate_kernel_numpy(chunk, kernel, rows, cols, ky, kx) + result = chunk.copy() + result[hy:-hy, hx:-hx] = interior + return result + + +def _morph_dask_numpy(data, kernel, op, boundary): + ky, kx = kernel.shape + hy, hx = ky // 2, kx // 2 + _func = partial(_morph_chunk_numpy, kernel=kernel, op=op) + return data.astype(np.float64).map_overlap( + _func, + depth=(hy, hx), + boundary=_boundary_to_dask(boundary), + meta=np.array(()), + ) + + +def _erode_dask_numpy(data, kernel, boundary): + return _morph_dask_numpy(data, kernel, 'erode', boundary) + + +def _dilate_dask_numpy(data, kernel, boundary): + return _morph_dask_numpy(data, kernel, 'dilate', boundary) + + +def _opening_dask_numpy(data, kernel, boundary): + tmp = _morph_dask_numpy(data, kernel, 'erode', boundary) + return _morph_dask_numpy(tmp, kernel, 'dilate', boundary) + + +def _closing_dask_numpy(data, kernel, boundary): + tmp = _morph_dask_numpy(data, kernel, 'dilate', boundary) + return _morph_dask_numpy(tmp, kernel, 'erode', boundary) + + +# --------------------------------------------------------------------------- +# dask + cupy backend +# --------------------------------------------------------------------------- + +def _morph_chunk_cupy(chunk, kernel, op, block_info=None): + """Process a single overlapped dask+cupy chunk.""" + import cupy as cp + + ky, kx = kernel.shape + hy, hx = ky // 2, kx // 2 + rows = chunk.shape[0] - 2 * hy + cols = chunk.shape[1] - 2 * hx + kern_dev = cp.asarray(kernel) + out = cp.empty((rows, cols), dtype=cp.float64) + bpg, tpb = calc_cuda_dims((rows, cols)) + if op == 'erode': + _erode_gpu[bpg, tpb](chunk, kern_dev, out, hy, hx, ky, kx) + else: + _dilate_gpu[bpg, tpb](chunk, kern_dev, out, hy, hx, ky, kx) + result = chunk.copy() + result[hy:-hy, hx:-hx] = out + return result + + +def _morph_dask_cupy(data, kernel, op, boundary): + import cupy as cp + + ky, kx = kernel.shape + hy, hx = ky // 2, kx // 2 + _func = partial(_morph_chunk_cupy, kernel=kernel, op=op) + return data.astype(cp.float64).map_overlap( + _func, + depth=(hy, hx), + boundary=_boundary_to_dask(boundary, is_cupy=True), + meta=cp.array(()), + ) + + +def _erode_dask_cupy(data, kernel, boundary): + return _morph_dask_cupy(data, kernel, 'erode', boundary) + + +def _dilate_dask_cupy(data, kernel, boundary): + return _morph_dask_cupy(data, kernel, 'dilate', boundary) + + +def _opening_dask_cupy(data, kernel, boundary): + tmp = _morph_dask_cupy(data, kernel, 'erode', boundary) + return _morph_dask_cupy(tmp, kernel, 'dilate', boundary) + + +def _closing_dask_cupy(data, kernel, boundary): + tmp = _morph_dask_cupy(data, kernel, 'dilate', boundary) + return _morph_dask_cupy(tmp, kernel, 'erode', boundary) + + +# --------------------------------------------------------------------------- +# Dispatcher helpers +# --------------------------------------------------------------------------- + +def _dispatch(agg, kernel, boundary, name, numpy_fn, cupy_fn, dask_fn, dask_cupy_fn): + """Common dispatch logic for all four morphological functions.""" + _validate_raster(agg, func_name=name, name='agg', ndim=2) + kernel = _validate_kernel(kernel, name) + _validate_boundary(boundary) + + mapper = ArrayTypeFunctionMapping( + numpy_func=partial(numpy_fn, kernel=kernel, boundary=boundary), + cupy_func=partial(cupy_fn, kernel=kernel, boundary=boundary), + dask_func=partial(dask_fn, kernel=kernel, boundary=boundary), + dask_cupy_func=partial(dask_cupy_fn, kernel=kernel, boundary=boundary), + ) + + out = mapper(agg)(agg.data) + + return DataArray( + out, + name=name, + dims=agg.dims, + coords=agg.coords, + attrs=agg.attrs, + ) + + +# --------------------------------------------------------------------------- +# Public API +# --------------------------------------------------------------------------- + +@supports_dataset +def morph_erode(agg, kernel=None, boundary='nan', name='erode'): + """Apply morphological erosion (local minimum) to a 2D raster. + + Each output cell receives the minimum value found within the + structuring element footprint centred on that cell. + + Parameters + ---------- + agg : xarray.DataArray + 2D raster of numeric values. + kernel : numpy.ndarray or None + 2D structuring element with odd dimensions. Non-zero entries + mark the footprint. Defaults to a 3x3 square + (``np.ones((3, 3), dtype=np.uint8)``). + boundary : str, default ``'nan'`` + Edge handling: ``'nan'``, ``'nearest'``, ``'reflect'``, or + ``'wrap'``. + name : str, default ``'erode'`` + Name for the output DataArray. + + Returns + ------- + xarray.DataArray + Eroded raster with the same shape, coords, dims, and attrs + as the input. + + Examples + -------- + >>> from xrspatial.morphology import morph_erode + >>> eroded = morph_erode(raster, kernel=np.ones((5, 5), dtype=np.uint8)) + """ + if kernel is None: + kernel = np.ones((3, 3), dtype=np.uint8) + return _dispatch( + agg, kernel, boundary, name, + _erode_numpy, _erode_cupy, _erode_dask_numpy, _erode_dask_cupy, + ) + + +@supports_dataset +def morph_dilate(agg, kernel=None, boundary='nan', name='dilate'): + """Apply morphological dilation (local maximum) to a 2D raster. + + Each output cell receives the maximum value found within the + structuring element footprint centred on that cell. + + Parameters + ---------- + agg : xarray.DataArray + 2D raster of numeric values. + kernel : numpy.ndarray or None + 2D structuring element with odd dimensions. Non-zero entries + mark the footprint. Defaults to a 3x3 square. + boundary : str, default ``'nan'`` + Edge handling: ``'nan'``, ``'nearest'``, ``'reflect'``, or + ``'wrap'``. + name : str, default ``'dilate'`` + Name for the output DataArray. + + Returns + ------- + xarray.DataArray + Dilated raster. + + Examples + -------- + >>> from xrspatial.morphology import morph_dilate + >>> dilated = morph_dilate(raster, kernel=np.ones((5, 5), dtype=np.uint8)) + """ + if kernel is None: + kernel = np.ones((3, 3), dtype=np.uint8) + return _dispatch( + agg, kernel, boundary, name, + _dilate_numpy, _dilate_cupy, _dilate_dask_numpy, _dilate_dask_cupy, + ) + + +@supports_dataset +def morph_opening(agg, kernel=None, boundary='nan', name='opening'): + """Apply morphological opening (erosion then dilation) to a 2D raster. + + Opening removes small bright features and salt noise while + preserving the overall shape of larger structures. + + Parameters + ---------- + agg : xarray.DataArray + 2D raster of numeric values. + kernel : numpy.ndarray or None + 2D structuring element with odd dimensions. Defaults to a + 3x3 square. + boundary : str, default ``'nan'`` + Edge handling: ``'nan'``, ``'nearest'``, ``'reflect'``, or + ``'wrap'``. + name : str, default ``'opening'`` + Name for the output DataArray. + + Returns + ------- + xarray.DataArray + Opened raster. + + Examples + -------- + >>> from xrspatial.morphology import morph_opening + >>> cleaned = morph_opening(binary_mask) + """ + if kernel is None: + kernel = np.ones((3, 3), dtype=np.uint8) + return _dispatch( + agg, kernel, boundary, name, + _opening_numpy, _opening_cupy, _opening_dask_numpy, _opening_dask_cupy, + ) + + +@supports_dataset +def morph_closing(agg, kernel=None, boundary='nan', name='closing'): + """Apply morphological closing (dilation then erosion) to a 2D raster. + + Closing fills small dark gaps and pepper noise while preserving + the overall shape of larger structures. + + Parameters + ---------- + agg : xarray.DataArray + 2D raster of numeric values. + kernel : numpy.ndarray or None + 2D structuring element with odd dimensions. Defaults to a + 3x3 square. + boundary : str, default ``'nan'`` + Edge handling: ``'nan'``, ``'nearest'``, ``'reflect'``, or + ``'wrap'``. + name : str, default ``'closing'`` + Name for the output DataArray. + + Returns + ------- + xarray.DataArray + Closed raster. + + Examples + -------- + >>> from xrspatial.morphology import morph_closing + >>> filled = morph_closing(binary_mask) + """ + if kernel is None: + kernel = np.ones((3, 3), dtype=np.uint8) + return _dispatch( + agg, kernel, boundary, name, + _closing_numpy, _closing_cupy, _closing_dask_numpy, _closing_dask_cupy, + ) diff --git a/xrspatial/tests/test_morphology.py b/xrspatial/tests/test_morphology.py new file mode 100644 index 00000000..503a9851 --- /dev/null +++ b/xrspatial/tests/test_morphology.py @@ -0,0 +1,375 @@ +"""Tests for xrspatial.morphology (erode, dilate, opening, closing).""" + +from __future__ import annotations + +from functools import partial + +import numpy as np +import pytest +import xarray as xr + +from xrspatial.morphology import ( + _circle_kernel, + morph_closing, + morph_dilate, + morph_erode, + morph_opening, +) +from xrspatial.tests.general_checks import ( + assert_boundary_mode_correctness, + assert_numpy_equals_cupy, + assert_numpy_equals_dask_cupy, + assert_numpy_equals_dask_numpy, + create_test_raster, + cuda_and_cupy_available, + dask_array_available, + general_output_checks, +) + + +# --------------------------------------------------------------------------- +# Fixtures / test data +# --------------------------------------------------------------------------- + +_DATA = np.array([ + [1., 5., 3., 2., 7.], + [4., 8., 6., 1., 3.], + [2., 9., 7., 5., 4.], + [3., 1., 4., 8., 6.], + [6., 2., 3., 7., 9.], +], dtype=np.float64) + +_DATA_NAN = _DATA.copy() +_DATA_NAN[0, 2] = np.nan +_DATA_NAN[3, 1] = np.nan + +_KERNEL_3x3 = np.ones((3, 3), dtype=np.uint8) + +_KERNEL_CROSS = np.array([ + [0, 1, 0], + [1, 1, 1], + [0, 1, 0], +], dtype=np.uint8) + + +# --------------------------------------------------------------------------- +# Kernel construction +# --------------------------------------------------------------------------- + +def test_circle_kernel_radius1(): + k = _circle_kernel(1) + expected = np.array([ + [0, 1, 0], + [1, 1, 1], + [0, 1, 0], + ], dtype=np.uint8) + np.testing.assert_array_equal(k, expected) + + +def test_circle_kernel_radius2(): + k = _circle_kernel(2) + assert k.shape == (5, 5) + # corners should be zero + assert k[0, 0] == 0 + assert k[0, 4] == 0 + # centre and cardinal neighbours should be one + assert k[2, 2] == 1 + assert k[0, 2] == 1 + + +# --------------------------------------------------------------------------- +# Validation +# --------------------------------------------------------------------------- + +def test_erode_rejects_non_dataarray(): + with pytest.raises(TypeError, match="must be an xarray.DataArray"): + morph_erode(np.ones((5, 5))) + + +def test_erode_rejects_3d_input(): + agg = xr.DataArray(np.ones((3, 5, 5)), dims=['z', 'y', 'x']) + with pytest.raises(ValueError, match="must be 2D"): + morph_erode(agg) + + +def test_erode_rejects_even_kernel(): + agg = create_test_raster(_DATA) + with pytest.raises(ValueError, match="must be odd"): + morph_erode(agg, kernel=np.ones((4, 4), dtype=np.uint8)) + + +def test_erode_rejects_bad_boundary(): + agg = create_test_raster(_DATA) + with pytest.raises(ValueError, match="boundary must be one of"): + morph_erode(agg, boundary='invalid') + + +# --------------------------------------------------------------------------- +# Correctness -- numpy backend +# --------------------------------------------------------------------------- + +def test_erode_numpy_correctness(): + """Verify erosion matches hand-computed local minima.""" + agg = create_test_raster(_DATA) + result = morph_erode(agg, kernel=_KERNEL_3x3) + general_output_checks(agg, result, verify_attrs=True) + + # Interior cell [2, 2]: min of 3x3 neighbourhood centred at (2,2) + # Neighbourhood: [[8,6,1],[9,7,5],[1,4,8]] + assert result.data[2, 2] == 1.0 + + # Interior cell [2, 1]: min of [[1,5,3],[4,8,6],[2,9,7]] + assert result.data[2, 1] == 1.0 + + +def test_dilate_numpy_correctness(): + """Verify dilation matches hand-computed local maxima.""" + agg = create_test_raster(_DATA) + result = morph_dilate(agg, kernel=_KERNEL_3x3) + general_output_checks(agg, result, verify_attrs=True) + + # Interior cell [2, 2]: max of 3x3 neighbourhood centred at (2,2) + # Neighbourhood: [[8,6,1],[9,7,5],[1,4,8]] + assert result.data[2, 2] == 9.0 + + +def test_erode_cross_kernel(): + """Verify cross kernel only considers cardinal neighbours + centre.""" + agg = create_test_raster(_DATA) + result = morph_erode(agg, kernel=_KERNEL_CROSS) + + # Interior cell [2, 2]: min of {data[1,2], data[2,1], data[2,2], data[2,3], data[3,2]} + # = min(6, 9, 7, 5, 4) = 4 + assert result.data[2, 2] == 4.0 + + +def test_opening_removes_bright_spike(): + """Opening should suppress isolated bright peaks.""" + data = np.zeros((7, 7), dtype=np.float64) + data[3, 3] = 100.0 # single bright spike + agg = create_test_raster(data) + result = morph_opening(agg, kernel=_KERNEL_3x3) + # After erosion the spike disappears; after dilation it stays gone + assert result.data[3, 3] == 0.0 + + +def test_closing_fills_dark_hole(): + """Closing should fill isolated dark pits.""" + data = np.full((7, 7), 100.0, dtype=np.float64) + data[3, 3] = 0.0 # single dark pit + agg = create_test_raster(data) + result = morph_closing(agg, kernel=_KERNEL_3x3) + assert result.data[3, 3] == 100.0 + + +def test_opening_idempotent(): + """Applying opening twice should give the same result as once.""" + agg = create_test_raster(_DATA) + once = morph_opening(agg, kernel=_KERNEL_3x3, boundary='nearest') + twice = morph_opening(once, kernel=_KERNEL_3x3, boundary='nearest') + np.testing.assert_allclose(once.data, twice.data, equal_nan=True) + + +def test_closing_idempotent(): + """Applying closing twice should give the same result as once.""" + agg = create_test_raster(_DATA) + once = morph_closing(agg, kernel=_KERNEL_3x3, boundary='nearest') + twice = morph_closing(once, kernel=_KERNEL_3x3, boundary='nearest') + np.testing.assert_allclose(once.data, twice.data, equal_nan=True) + + +# --------------------------------------------------------------------------- +# NaN handling +# --------------------------------------------------------------------------- + +def test_erode_nan_propagation(): + """NaN cells in the source should remain NaN; NaN neighbours propagate.""" + agg = create_test_raster(_DATA_NAN) + result = morph_erode(agg, kernel=_KERNEL_3x3) + # Cell [0, 2] is NaN in source -> stays NaN + assert np.isnan(result.data[0, 2]) + # Cell [0, 1] touches NaN at [0, 2] -> becomes NaN + assert np.isnan(result.data[0, 1]) + + +def test_dilate_nan_propagation(): + agg = create_test_raster(_DATA_NAN) + result = morph_dilate(agg, kernel=_KERNEL_3x3) + assert np.isnan(result.data[0, 2]) + assert np.isnan(result.data[0, 1]) + + +# --------------------------------------------------------------------------- +# Edge cases +# --------------------------------------------------------------------------- + +def test_single_cell_raster(): + data = np.array([[42.0]], dtype=np.float64) + agg = create_test_raster(data) + # With boundary='nan', the NaN padding propagates into the result. + assert np.isnan(morph_erode(agg).data[0, 0]) + # With boundary='nearest', padding replicates the single value. + assert morph_erode(agg, boundary='nearest').data[0, 0] == 42.0 + assert morph_dilate(agg, boundary='nearest').data[0, 0] == 42.0 + + +def test_uniform_raster(): + """All values identical -> all operations should be identity.""" + data = np.full((5, 5), 7.0, dtype=np.float64) + agg = create_test_raster(data) + for func in [morph_erode, morph_dilate, morph_opening, morph_closing]: + result = func(agg, kernel=_KERNEL_3x3, boundary='nearest') + np.testing.assert_allclose(result.data, 7.0) + + +def test_default_kernel(): + """Functions work with the default (None) kernel argument.""" + agg = create_test_raster(_DATA) + result = morph_erode(agg) + general_output_checks(agg, result, verify_attrs=True) + + +# --------------------------------------------------------------------------- +# Boundary modes -- numpy vs dask +# --------------------------------------------------------------------------- + +@dask_array_available +def test_erode_boundary_modes(): + numpy_agg = create_test_raster(_DATA, backend='numpy') + dask_agg = create_test_raster(_DATA, backend='dask') + assert_boundary_mode_correctness( + numpy_agg, dask_agg, morph_erode, depth=1, nan_edges=True, + ) + + +@dask_array_available +def test_dilate_boundary_modes(): + numpy_agg = create_test_raster(_DATA, backend='numpy') + dask_agg = create_test_raster(_DATA, backend='dask') + assert_boundary_mode_correctness( + numpy_agg, dask_agg, morph_dilate, depth=1, nan_edges=True, + ) + + +# --------------------------------------------------------------------------- +# Dask backend +# --------------------------------------------------------------------------- + +@dask_array_available +def test_erode_dask_equals_numpy(): + numpy_agg = create_test_raster(_DATA, backend='numpy') + dask_agg = create_test_raster(_DATA, backend='dask') + assert_numpy_equals_dask_numpy(numpy_agg, dask_agg, morph_erode, nan_edges=True) + + +@dask_array_available +def test_dilate_dask_equals_numpy(): + numpy_agg = create_test_raster(_DATA, backend='numpy') + dask_agg = create_test_raster(_DATA, backend='dask') + assert_numpy_equals_dask_numpy(numpy_agg, dask_agg, morph_dilate, nan_edges=True) + + +@dask_array_available +def test_opening_dask_equals_numpy(): + numpy_agg = create_test_raster(_DATA, backend='numpy') + dask_agg = create_test_raster(_DATA, backend='dask') + assert_numpy_equals_dask_numpy(numpy_agg, dask_agg, morph_opening, nan_edges=True) + + +@dask_array_available +def test_closing_dask_equals_numpy(): + numpy_agg = create_test_raster(_DATA, backend='numpy') + dask_agg = create_test_raster(_DATA, backend='dask') + assert_numpy_equals_dask_numpy(numpy_agg, dask_agg, morph_closing, nan_edges=True) + + +# --------------------------------------------------------------------------- +# CuPy backend +# --------------------------------------------------------------------------- + +@cuda_and_cupy_available +def test_erode_cupy_equals_numpy(): + numpy_agg = create_test_raster(_DATA, backend='numpy') + cupy_agg = create_test_raster(_DATA, backend='cupy') + assert_numpy_equals_cupy(numpy_agg, cupy_agg, morph_erode, nan_edges=True) + + +@cuda_and_cupy_available +def test_dilate_cupy_equals_numpy(): + numpy_agg = create_test_raster(_DATA, backend='numpy') + cupy_agg = create_test_raster(_DATA, backend='cupy') + assert_numpy_equals_cupy(numpy_agg, cupy_agg, morph_dilate, nan_edges=True) + + +@cuda_and_cupy_available +def test_opening_cupy_equals_numpy(): + numpy_agg = create_test_raster(_DATA, backend='numpy') + cupy_agg = create_test_raster(_DATA, backend='cupy') + assert_numpy_equals_cupy(numpy_agg, cupy_agg, morph_opening, nan_edges=True) + + +@cuda_and_cupy_available +def test_closing_cupy_equals_numpy(): + numpy_agg = create_test_raster(_DATA, backend='numpy') + cupy_agg = create_test_raster(_DATA, backend='cupy') + assert_numpy_equals_cupy(numpy_agg, cupy_agg, morph_closing, nan_edges=True) + + +# --------------------------------------------------------------------------- +# Dask + CuPy backend +# --------------------------------------------------------------------------- + +@cuda_and_cupy_available +@dask_array_available +def test_erode_dask_cupy_equals_numpy(): + numpy_agg = create_test_raster(_DATA, backend='numpy') + dask_cupy_agg = create_test_raster(_DATA, backend='dask+cupy') + assert_numpy_equals_dask_cupy(numpy_agg, dask_cupy_agg, morph_erode, nan_edges=True) + + +@cuda_and_cupy_available +@dask_array_available +def test_dilate_dask_cupy_equals_numpy(): + numpy_agg = create_test_raster(_DATA, backend='numpy') + dask_cupy_agg = create_test_raster(_DATA, backend='dask+cupy') + assert_numpy_equals_dask_cupy(numpy_agg, dask_cupy_agg, morph_dilate, nan_edges=True) + + +# --------------------------------------------------------------------------- +# 5x5 kernel +# --------------------------------------------------------------------------- + +def test_erode_5x5_kernel(): + """Larger kernel should produce a wider erosion effect.""" + data = np.full((9, 9), 10.0, dtype=np.float64) + data[4, 4] = 0.0 # single dark pit + agg = create_test_raster(data) + k3 = np.ones((3, 3), dtype=np.uint8) + k5 = np.ones((5, 5), dtype=np.uint8) + + r3 = morph_erode(agg, kernel=k3, boundary='nearest') + r5 = morph_erode(agg, kernel=k5, boundary='nearest') + + # 3x3 kernel: cell [3,3] should be 0 (touches the pit) + assert r3.data[3, 3] == 0.0 + # but cell [2,2] should remain 10 + assert r3.data[2, 2] == 10.0 + + # 5x5 kernel: cell [2,2] should now see the pit too + assert r5.data[2, 2] == 0.0 + + +# --------------------------------------------------------------------------- +# Dataset support +# --------------------------------------------------------------------------- + +def test_erode_dataset(): + """morph_erode should work on xr.Dataset via @supports_dataset.""" + data = np.random.default_rng(42).random((5, 5)).astype(np.float64) + ds = xr.Dataset({ + 'a': xr.DataArray(data, dims=['y', 'x']), + 'b': xr.DataArray(data * 2, dims=['y', 'x']), + }) + result = morph_erode(ds, boundary='nearest') + assert isinstance(result, xr.Dataset) + assert set(result.data_vars) == {'a', 'b'}