これは何？
ガウシアンフィルター（Gaussian Filter） は、画像や格子データを「正規分布の重み」でなめらかにする処理だ。
森林LiDARでの使いどころ：
- 樹高マップ（CHM）のノイズ除去
- 樹頂点検出の前処理（これが主役）
- 地形モデルのスムージング
「ぼかすだけ」に見えて、次の検出ステップの精度を左右する重要な処理だ。
なぜ樹頂点検出の前にぼかすのか
1本の木でも、LiDAR点群の測量誤差や葉の揺れで 複数のピーク が生まれる。
【フィルターなし】 【フィルターあり】
高さ 高さ
↑ ↑
30│ ▲ ▲ 28│ ▲
│/|\ |\ ←誤ピーク │ /│\
25│ | \ | \ │ / │ \
│ | \| \ │ / │ \
└──────────→ X └──────────→ X
→ 1本の木を3本と誤検出 → 1本として正しく検出
ガウシアンフィルターで表面をなめらかにしてから探すことで、余分なピークを消せる。
アルゴリズム：重みが正規分布
各ピクセルの値を「周囲のピクセルの重み付き平均」に置き換える。
重みの式：
w(d) = exp( -d² / (2σ²) )
d = 中心からの距離（ピクセル）
σ = sigma パラメータ
これは正規分布の確率密度関数と同じ形だ。中心が最大で、離れるほど指数関数的に小さくなる。
sigma=1.5 のときの重み分布:
距離 d 重み w
0 1.00 ████████████████
1 0.64 ██████████
2 0.17 ███
3 0.01 ▏
4 0.00 （無視）
2Dに広げると：
0.01 0.12 0.20 0.12 0.01
0.12 0.64 1.00 0.64 0.12
0.20 1.00 [中心] 1.00 0.20
0.12 0.64 1.00 0.64 0.12
0.01 0.12 0.20 0.12 0.01
中心に近いピクセルほど影響が大きく、遠いほど無視される。
Pythonでの実装
import numpy as np
import rasterio
from scipy.ndimage import gaussian_filter
# CHMを読み込む
with rasterio.open('../data/processed/height_map.tif') as src:
chm = src.read(1).astype('float32')
chm[chm == src.nodata] = np.nan
# NaNを0に置換（NaNのまま渡すと周囲に伝染する）
chm_filled = np.where(np.isnan(chm), 0, chm)
# ガウシアンフィルター適用
sigma = 1.5
chm_smooth = gaussian_filter(chm_filled, sigma=sigma)
# 元のNaN箇所を戻す
chm_smooth[np.isnan(chm)] = np.nan
NaN の扱いが重要：gaussian_filter は NaN を「値0」として計算する。フィルター後に元の NaN 箇所を戻さないと、データなし領域が「樹高0m」として残り、境界部分の検出が乱れる。
sigma パラメータの意味
御岳山のデータは解像度 1m/pixel なので：
sigma = 1.5 → 1.5ピクセル = 1.5m の範囲でぼかす
→ 半径3σ = 4.5m 以上離れた点はほぼ無視（重み < 0.01）
統計の「平均±1σに68%のデータが入る」と同じ σ だ。
sigma の選び方
for sigma in [0.5, 1.0, 1.5, 2.0, 3.0]:
chm_smooth = gaussian_filter(chm_filled, sigma=sigma)
local_max = maximum_filter(chm_smooth, size=5)
tree_mask = (chm_smooth == local_max) & (chm_smooth >= 5.0)
print(f'sigma={sigma}: {tree_mask.sum()}本')
sigma=0.5: 1,203本 ← ノイズが残る、過検出
sigma=1.0: 1,018本
sigma=1.5: 941本 ← 御岳山での採用値
sigma=2.0: 762本
sigma=3.0: 498本 ← 隣の木と合体、過少検出
検出本数は sigma を大きくするにつれ単調に減る。「肘」は出ないので、現地の実測データや目視と照合して決める。
御岳山での処理結果
データ: 09KC7495.las（御岳山 400m × 300m）
解像度: 1m/pixel → 401 × 301 ピクセル
sigma: 1.5
フィルター前後のピーク数の違い：
フィルターなし: 局所最大値が 2,400+ 個（ノイズ込み）
sigma=1.5 後: 941 個（1本1ピーク に整理）
次のステップ：maximum_filter で樹頂点を探す
ガウシアンフィルターでなめらかにした CHM に対して、maximum_filter で「周囲より高いピクセル = 樹頂点」を探す。
from scipy.ndimage import maximum_filter
WINDOW_SIZE = 5 # 5×5ピクセル（= 5m×5m）の窓で探索
MIN_HEIGHT = 5.0
local_max = maximum_filter(chm_smooth, size=WINDOW_SIZE)
tree_mask = (
(chm_smooth == local_max) & # 周囲で最大
(chm_smooth >= MIN_HEIGHT) & # 最低樹高以上
~np.isnan(chm)
)
print(f'検出本数: {tree_mask.sum()}')
ガウシアンフィルターが前処理として機能して初めて、maximum_filter が1本1ピークを正しく拾えるようになる。
まとめ
| 項目 | 内容 |
|---|---|
| 何をするか | 正規分布の重みで周囲のピクセルを平均し、なめらかにする |
| 核心 |
w = exp(-d² / 2σ²) — 距離が遠いほど影響小 |
| sigma | ぼかしの半径（ピクセル単位）。解像度1m/pxなら sigma=1.5 → 1.5m |
| NaNの注意 | フィルター前に0置換、フィルター後に元のNaN箇所を戻す |
| sigma の決め方 | 単調減少するので「肘」なし。実測・目視で検証する |
chm_filled = np.where(np.isnan(chm), 0, chm)
chm_smooth = gaussian_filter(chm_filled, sigma=1.5)
chm_smooth[np.isnan(chm)] = np.nan
このシリーズは、東京都多摩地域の森林LiDARデータをPythonで分析する中で出会ったアルゴリズムを1つずつ解説しています。
実装コードは tama-forest notebooks で公開中。
作者のその他の実験 → https://niikun.net