Skip to content

Commit ec31dc6

Browse files
brendancolbrendancol
authored
Add bilateral filter for feature-preserving smoothing (#969) (#970)
1 parent b66c00e commit ec31dc6

File tree

7 files changed

+814
-0
lines changed

7 files changed

+814
-0
lines changed

‎README.md‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -168,6 +168,7 @@ In the GIS world, rasters are used for representing continuous phenomena (e.g. e
168168
| [Emerging Hotspots](xrspatial/emerging_hotspots.py) | Classifies time-series hot/cold spot trends using Gi* and Mann-Kendall | ✅️ | ✅️ | ✅️ | ✅️ |
169169
| [Mean](xrspatial/focal.py) | Computes the mean value within a sliding neighborhood window | ✅️ | ✅️ | ✅️ | ✅️ |
170170
| [Focal Statistics](xrspatial/focal.py) | Computes summary statistics over a sliding neighborhood window | ✅️ | ✅️ | ✅️ | ✅️ |
171+
| [Bilateral](xrspatial/bilateral.py) | Feature-preserving smoothing via bilateral filtering | ✅️ | ✅️ | ✅️ | ✅️ |
171172

172173
-------
173174

‎docs/source/reference/focal.rst‎

Lines changed: 7 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -25,6 +25,13 @@ Mean
2525

2626
xrspatial.focal.mean
2727

28+
Bilateral
29+
=========
30+
.. autosummary::
31+
:toctree: _autosummary
32+
33+
xrspatial.bilateral.bilateral
34+
2835

2936
Focal Statistics
3037
================

‎examples/user_guide/17_Bilateral_Filter.ipynb‎

Lines changed: 185 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,185 @@
1+
{
2+
"cells": [
3+
{
4+
"cell_type": "markdown",
5+
"metadata": {},
6+
"source": [
7+
"# Bilateral filtering\n",
8+
"\n",
9+
"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",
10+
"\n",
11+
"Two parameters control the behavior:\n",
12+
"- **sigma_spatial**: how far the spatial Gaussian reaches (kernel radius = ceil(2 * sigma_spatial))\n",
13+
"- **sigma_range**: how much value difference is tolerated before a neighbor gets downweighted"
14+
]
15+
},
16+
{
17+
"cell_type": "code",
18+
"execution_count": null,
19+
"metadata": {},
20+
"outputs": [],
21+
"source": [
22+
"import numpy as np\n",
23+
"import xarray as xr\n",
24+
"import matplotlib.pyplot as plt\n",
25+
"\n",
26+
"from xrspatial import bilateral\n",
27+
"from xrspatial import mean\n",
28+
"from xrspatial.terrain import generate_terrain"
29+
]
30+
},
31+
{
32+
"cell_type": "markdown",
33+
"metadata": {},
34+
"source": [
35+
"## Generate a synthetic terrain with noise\n",
36+
"\n",
37+
"We'll create a DEM, add Gaussian noise, and then compare bilateral filtering against the standard mean filter."
38+
]
39+
},
40+
{
41+
"cell_type": "code",
42+
"execution_count": null,
43+
"metadata": {},
44+
"outputs": [],
45+
"source": [
46+
"W, H = 600, 400\n",
47+
"cvs_terrain = xr.DataArray(\n",
48+
" np.zeros((H, W)),\n",
49+
" dims=['y', 'x'],\n",
50+
" coords={'y': np.linspace(0, 100, H), 'x': np.linspace(0, 150, W)},\n",
51+
")\n",
52+
"terrain = generate_terrain(cvs_terrain, seed=42)\n",
53+
"\n",
54+
"# Add noise\n",
55+
"rng = np.random.default_rng(123)\n",
56+
"noise = rng.normal(0, 15, terrain.shape)\n",
57+
"noisy_terrain = terrain.copy(data=terrain.values + noise)\n",
58+
"\n",
59+
"fig, axes = plt.subplots(1, 2, figsize=(14, 5))\n",
60+
"terrain.plot(ax=axes[0], cmap='terrain')\n",
61+
"axes[0].set_title('Clean terrain')\n",
62+
"noisy_terrain.plot(ax=axes[1], cmap='terrain')\n",
63+
"axes[1].set_title('Noisy terrain')\n",
64+
"plt.tight_layout()"
65+
]
66+
},
67+
{
68+
"cell_type": "markdown",
69+
"metadata": {},
70+
"source": [
71+
"## Compare bilateral vs. mean filter\n",
72+
"\n",
73+
"The mean filter blurs edges. The bilateral filter preserves them."
74+
]
75+
},
76+
{
77+
"cell_type": "code",
78+
"execution_count": null,
79+
"metadata": {},
80+
"outputs": [],
81+
"source": [
82+
"smoothed_bilateral = bilateral(noisy_terrain, sigma_spatial=2.0, sigma_range=20.0)\n",
83+
"smoothed_mean = mean(noisy_terrain, passes=3)\n",
84+
"\n",
85+
"fig, axes = plt.subplots(1, 3, figsize=(18, 5))\n",
86+
"\n",
87+
"noisy_terrain.plot(ax=axes[0], cmap='terrain')\n",
88+
"axes[0].set_title('Noisy input')\n",
89+
"\n",
90+
"smoothed_mean.plot(ax=axes[1], cmap='terrain')\n",
91+
"axes[1].set_title('Mean filter (3 passes)')\n",
92+
"\n",
93+
"smoothed_bilateral.plot(ax=axes[2], cmap='terrain')\n",
94+
"axes[2].set_title('Bilateral filter')\n",
95+
"\n",
96+
"plt.tight_layout()"
97+
]
98+
},
99+
{
100+
"cell_type": "markdown",
101+
"metadata": {},
102+
"source": [
103+
"## Effect of sigma_range\n",
104+
"\n",
105+
"Smaller `sigma_range` preserves more edges; larger values allow smoothing across bigger value differences."
106+
]
107+
},
108+
{
109+
"cell_type": "code",
110+
"execution_count": null,
111+
"metadata": {},
112+
"outputs": [],
113+
"source": [
114+
"sigma_ranges = [5.0, 20.0, 100.0]\n",
115+
"\n",
116+
"fig, axes = plt.subplots(1, len(sigma_ranges), figsize=(18, 5))\n",
117+
"for ax, sr in zip(axes, sigma_ranges):\n",
118+
" result = bilateral(noisy_terrain, sigma_spatial=2.0, sigma_range=sr)\n",
119+
" result.plot(ax=ax, cmap='terrain')\n",
120+
" ax.set_title(f'sigma_range = {sr}')\n",
121+
"plt.tight_layout()"
122+
]
123+
},
124+
{
125+
"cell_type": "markdown",
126+
"metadata": {},
127+
"source": [
128+
"## Step-edge preservation\n",
129+
"\n",
130+
"A clear demonstration: a raster with a sharp vertical edge. The bilateral filter keeps the boundary; the mean filter blurs it."
131+
]
132+
},
133+
{
134+
"cell_type": "code",
135+
"execution_count": null,
136+
"metadata": {},
137+
"outputs": [],
138+
"source": [
139+
"step = np.zeros((50, 100))\n",
140+
"step[:, 50:] = 100.0\n",
141+
"\n",
142+
"# Add a bit of noise\n",
143+
"step_noisy = step + rng.normal(0, 5, step.shape)\n",
144+
"step_agg = xr.DataArray(step_noisy, dims=['y', 'x'])\n",
145+
"\n",
146+
"step_bilateral = bilateral(step_agg, sigma_spatial=2.0, sigma_range=10.0)\n",
147+
"step_mean = mean(step_agg, passes=3)\n",
148+
"\n",
149+
"fig, axes = plt.subplots(1, 3, figsize=(15, 4))\n",
150+
"step_agg.plot(ax=axes[0], cmap='gray')\n",
151+
"axes[0].set_title('Noisy step edge')\n",
152+
"step_mean.plot(ax=axes[1], cmap='gray')\n",
153+
"axes[1].set_title('Mean filter')\n",
154+
"step_bilateral.plot(ax=axes[2], cmap='gray')\n",
155+
"axes[2].set_title('Bilateral filter')\n",
156+
"plt.tight_layout()\n",
157+
"\n",
158+
"# Cross-section\n",
159+
"row = 25\n",
160+
"fig, ax = plt.subplots(figsize=(10, 4))\n",
161+
"ax.plot(step_agg.data[row], label='Noisy', alpha=0.5)\n",
162+
"ax.plot(step_mean.data[row], label='Mean', linewidth=2)\n",
163+
"ax.plot(step_bilateral.data[row], label='Bilateral', linewidth=2)\n",
164+
"ax.legend()\n",
165+
"ax.set_xlabel('Column')\n",
166+
"ax.set_ylabel('Value')\n",
167+
"ax.set_title('Cross-section at row 25')\n",
168+
"plt.tight_layout()"
169+
]
170+
}
171+
],
172+
"metadata": {
173+
"kernelspec": {
174+
"display_name": "Python 3",
175+
"language": "python",
176+
"name": "python3"
177+
},
178+
"language_info": {
179+
"name": "python",
180+
"version": "3.10.0"
181+
}
182+
},
183+
"nbformat": 4,
184+
"nbformat_minor": 4
185+
}

‎xrspatial/__init__.py‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,6 @@
11
from xrspatial.aspect import aspect # noqa
22
from xrspatial.balanced_allocation import balanced_allocation # noqa
3+
from xrspatial.bilateral import bilateral # noqa
34
from xrspatial.bump import bump # noqa
45
from xrspatial.corridor import least_cost_corridor # noqa
56
from xrspatial.cost_distance import cost_distance # noqa

‎xrspatial/accessor.py‎

Lines changed: 8 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -181,6 +181,10 @@ def focal_mean(self, **kwargs):
181181
from .focal import mean
182182
return mean(self._obj, **kwargs)
183183

184+
def bilateral(self, **kwargs):
185+
from .bilateral import bilateral
186+
return bilateral(self._obj, **kwargs)
187+
184188
# ---- Morphological ----
185189

186190
def morph_erode(self, **kwargs):
@@ -545,6 +549,10 @@ def focal_mean(self, **kwargs):
545549
from .focal import mean
546550
return mean(self._obj, **kwargs)
547551

552+
def bilateral(self, **kwargs):
553+
from .bilateral import bilateral
554+
return bilateral(self._obj, **kwargs)
555+
548556
# ---- Morphological ----
549557

550558
def morph_erode(self, **kwargs):

0 commit comments

Comments
 (0)