Add bilateral filter for feature-preserving smoothing (#969)

README.md

@@ -168,6 +168,7 @@ In the GIS world, rasters are used for representing continuous phenomena (e.g. e
 | [Emerging Hotspots](xrspatial/emerging_hotspots.py) | Classifies time-series hot/cold spot trends using Gi* and Mann-Kendall | ✅️ | ✅️ | ✅️ | ✅️ |
 | [Mean](xrspatial/focal.py) | Computes the mean value within a sliding neighborhood window | ✅️ | ✅️ | ✅️ | ✅️ |
 | [Focal Statistics](xrspatial/focal.py) | Computes summary statistics over a sliding neighborhood window | ✅️ | ✅️ | ✅️ | ✅️ |
+| [Bilateral](xrspatial/bilateral.py) | Feature-preserving smoothing via bilateral filtering | ✅️ | ✅️ | ✅️ | ✅️ |
 
 -------
 

docs/source/reference/focal.rst

@@ -25,6 +25,13 @@ Mean
 
    xrspatial.focal.mean
 
+Bilateral
+=========
+.. autosummary::
+   :toctree: _autosummary
+
+   xrspatial.bilateral.bilateral
+
 
 Focal Statistics
 ================

examples/user_guide/17_Bilateral_Filter.ipynb

@@ -0,0 +1,185 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Bilateral filtering\n",
+    "\n",
+    "The bilateral filter smooths a raster while preserving edges. Unlike a simple mean filter, it weights each neighbor by both spatial distance and value similarity. Pixels across a sharp boundary contribute very little, so edges stay sharp while flat areas get smoothed.\n",
+    "\n",
+    "Two parameters control the behavior:\n",
+    "- **sigma_spatial**: how far the spatial Gaussian reaches (kernel radius = ceil(2 * sigma_spatial))\n",
+    "- **sigma_range**: how much value difference is tolerated before a neighbor gets downweighted"
+   ]
+  },
+  {
+   "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 import bilateral\n",
+    "from xrspatial import mean\n",
+    "from xrspatial.terrain import generate_terrain"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Generate a synthetic terrain with noise\n",
+    "\n",
+    "We'll create a DEM, add Gaussian noise, and then compare bilateral filtering against the standard mean filter."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "W, H = 600, 400\n",
+    "cvs_terrain = xr.DataArray(\n",
+    "    np.zeros((H, W)),\n",
+    "    dims=['y', 'x'],\n",
+    "    coords={'y': np.linspace(0, 100, H), 'x': np.linspace(0, 150, W)},\n",
+    ")\n",
+    "terrain = generate_terrain(cvs_terrain, seed=42)\n",
+    "\n",
+    "# Add noise\n",
+    "rng = np.random.default_rng(123)\n",
+    "noise = rng.normal(0, 15, terrain.shape)\n",
+    "noisy_terrain = terrain.copy(data=terrain.values + noise)\n",
+    "\n",
+    "fig, axes = plt.subplots(1, 2, figsize=(14, 5))\n",
+    "terrain.plot(ax=axes[0], cmap='terrain')\n",
+    "axes[0].set_title('Clean terrain')\n",
+    "noisy_terrain.plot(ax=axes[1], cmap='terrain')\n",
+    "axes[1].set_title('Noisy terrain')\n",
+    "plt.tight_layout()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Compare bilateral vs. mean filter\n",
+    "\n",
+    "The mean filter blurs edges. The bilateral filter preserves them."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "smoothed_bilateral = bilateral(noisy_terrain, sigma_spatial=2.0, sigma_range=20.0)\n",
+    "smoothed_mean = mean(noisy_terrain, passes=3)\n",
+    "\n",
+    "fig, axes = plt.subplots(1, 3, figsize=(18, 5))\n",
+    "\n",
+    "noisy_terrain.plot(ax=axes[0], cmap='terrain')\n",
+    "axes[0].set_title('Noisy input')\n",
+    "\n",
+    "smoothed_mean.plot(ax=axes[1], cmap='terrain')\n",
+    "axes[1].set_title('Mean filter (3 passes)')\n",
+    "\n",
+    "smoothed_bilateral.plot(ax=axes[2], cmap='terrain')\n",
+    "axes[2].set_title('Bilateral filter')\n",
+    "\n",
+    "plt.tight_layout()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Effect of sigma_range\n",
+    "\n",
+    "Smaller `sigma_range` preserves more edges; larger values allow smoothing across bigger value differences."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "sigma_ranges = [5.0, 20.0, 100.0]\n",
+    "\n",
+    "fig, axes = plt.subplots(1, len(sigma_ranges), figsize=(18, 5))\n",
+    "for ax, sr in zip(axes, sigma_ranges):\n",
+    "    result = bilateral(noisy_terrain, sigma_spatial=2.0, sigma_range=sr)\n",
+    "    result.plot(ax=ax, cmap='terrain')\n",
+    "    ax.set_title(f'sigma_range = {sr}')\n",
+    "plt.tight_layout()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Step-edge preservation\n",
+    "\n",
+    "A clear demonstration: a raster with a sharp vertical edge. The bilateral filter keeps the boundary; the mean filter blurs it."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "step = np.zeros((50, 100))\n",
+    "step[:, 50:] = 100.0\n",
+    "\n",
+    "# Add a bit of noise\n",
+    "step_noisy = step + rng.normal(0, 5, step.shape)\n",
+    "step_agg = xr.DataArray(step_noisy, dims=['y', 'x'])\n",
+    "\n",
+    "step_bilateral = bilateral(step_agg, sigma_spatial=2.0, sigma_range=10.0)\n",
+    "step_mean = mean(step_agg, passes=3)\n",
+    "\n",
+    "fig, axes = plt.subplots(1, 3, figsize=(15, 4))\n",
+    "step_agg.plot(ax=axes[0], cmap='gray')\n",
+    "axes[0].set_title('Noisy step edge')\n",
+    "step_mean.plot(ax=axes[1], cmap='gray')\n",
+    "axes[1].set_title('Mean filter')\n",
+    "step_bilateral.plot(ax=axes[2], cmap='gray')\n",
+    "axes[2].set_title('Bilateral filter')\n",
+    "plt.tight_layout()\n",
+    "\n",
+    "# Cross-section\n",
+    "row = 25\n",
+    "fig, ax = plt.subplots(figsize=(10, 4))\n",
+    "ax.plot(step_agg.data[row], label='Noisy', alpha=0.5)\n",
+    "ax.plot(step_mean.data[row], label='Mean', linewidth=2)\n",
+    "ax.plot(step_bilateral.data[row], label='Bilateral', linewidth=2)\n",
+    "ax.legend()\n",
+    "ax.set_xlabel('Column')\n",
+    "ax.set_ylabel('Value')\n",
+    "ax.set_title('Cross-section at row 25')\n",
+    "plt.tight_layout()"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "name": "python",
+   "version": "3.10.0"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}

xrspatial/__init__.py

@@ -1,5 +1,6 @@
 from xrspatial.aspect import aspect  # noqa
 from xrspatial.balanced_allocation import balanced_allocation  # noqa
+from xrspatial.bilateral import bilateral  # noqa
 from xrspatial.bump import bump  # noqa
 from xrspatial.corridor import least_cost_corridor  # noqa
 from xrspatial.cost_distance import cost_distance  # noqa

xrspatial/accessor.py

@@ -181,6 +181,10 @@ def focal_mean(self, **kwargs):
         from .focal import mean
         return mean(self._obj, **kwargs)
 
+    def bilateral(self, **kwargs):
+        from .bilateral import bilateral
+        return bilateral(self._obj, **kwargs)
+
     # ---- Morphological ----
 
     def morph_erode(self, **kwargs):
@@ -545,6 +549,10 @@ def focal_mean(self, **kwargs):
         from .focal import mean
         return mean(self._obj, **kwargs)
 
+    def bilateral(self, **kwargs):
+        from .bilateral import bilateral
+        return bilateral(self._obj, **kwargs)
+
     # ---- Morphological ----
 
     def morph_erode(self, **kwargs):