これは何？

ガウシアンフィルター（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