これは何？
点群データには必ずノイズが混ざっている。鳥、電線、測量誤差……。
そのノイズに引っ張られず、「本当のモデル（平面・直線など）」を見つけるアルゴリズムが RANSAC だ。
森林LiDARでの主な使いどころ：
- 地面（DEM）の検出
- 個木の幹軸の推定
- 建物屋根面の検出
普通のフィッティングとの違い
最小二乗法: RANSAC:
● 外れ値 ● 外れ値（無視）
↗（引っ張られる）
● ● ● ● ● ● ● ●
────────── ← ずれる ────────── ← 正確
最小二乗法は全点を均等に扱うため、外れ値が多いと結果がズレる。
RANSACは「支持されなかった点は無視する」ことで外れ値に強くなる。
アルゴリズムの流れ
n_iter 回繰り返す:
1. 最小限の点をランダムに選ぶ（平面なら3点）
2. 仮のモデル（平面）を作る
3. 全点との誤差（距離）を計算
4. 距離 < threshold → インライア（モデルに合う点）
5. インライア数が最多なら採用
最後に:
6. 採用したインライア全体で再フィット（精度アップ）
Pythonで一から実装する
import numpy as np
def ransac_plane(points, threshold=0.5, n_iter=500, rng=None):
"""RANSAC による平面検出。(法線, 切片, インライアマスク) を返す。"""
if rng is None:
rng = np.random.default_rng(42)
best_mask = np.zeros(len(points), dtype=bool)
for _ in range(n_iter):
# 1. 3点をランダムに選ぶ
idx = rng.choice(len(points), 3, replace=False)
p1, p2, p3 = points[idx]
# 2. 法線ベクトルを計算（外積 → 長さ1に正規化）
n = np.cross(p2 - p1, p3 - p1)
norm = np.linalg.norm(n)
if norm < 1e-10:
continue # 3点が同一直線上 → スキップ
n = n / norm
# 3. 全点から平面までの距離
dist = np.abs((points - p1) @ n)
# 4. インライア判定
mask = dist < threshold
# 5. 最多インライアを記録
if mask.sum() > best_mask.sum():
best_mask = mask
# 6. 再フィット：インライア全体でSVDによる最適平面を求める
inliers = points[best_mask]
centroid = inliers.mean(axis=0)
_, _, Vt = np.linalg.svd(inliers - centroid, full_matrices=False)
normal = Vt[-1] # 分散最小方向 = 平面の法線
d = -normal @ centroid
return normal, d, best_mask
コアの計算を1行ずつ読む
法線ベクトルを求める
n = np.cross(p2 - p1, p3 - p1)
3点が張る2本のベクトルの外積を取ると、その平面に垂直なベクトルが得られる。
p3
↑ ↑ n（法線）= 平面に垂直
| |
p1 → p2 ───┘
長さを1に揃える（単位ベクトル化）
norm = np.linalg.norm(n) # ベクトルの長さ = √(nx²+ny²+nz²)
n = n / norm
長さを1にしないと、次の距離計算の単位がメートルにならない。
点から平面までの距離
dist = np.abs((points - p1) @ n)
平面の方程式 n·(x - p1) = 0 を使う。各点 x に対して |n·(x - p1)| が平面からの距離。
再フィット（SVD）
_, _, Vt = np.linalg.svd(inliers - centroid, full_matrices=False)
normal = Vt[-1]
重心を原点にずらした点群にSVDをかける。Vt[-1] は分散が最も小さい方向、つまり平面の薄さ方向 = 法線。
★full_matrices=Falseにしないと、巨大な行列の計算になり、クラッシュしてしまいます。
デフォルト（True）では U 行列が（点数 × 点数）で確保されるため、数十万点のデータに適用すると TB 単位のメモリが必要になります。
False にすると U は（点数 × 3）になり、Vt の値は変わりません。
なぜ何百回も試すのか（確率の話）
外れ値が30%混じっているとき、3点を選んで全部インライアである確率：
P(成功) = 0.7³ ≈ 34%
500回試したとき1回も成功しない確率 ≈ 10⁻¹¹⁰（事実上ゼロ）
試行回数を増やすことで、「偶然ノイズ点だけを引く」リスクを排除できる。
森林データへの適用例
import laspy
las = laspy.read('mitake.las')
xyz = np.vstack([las.x, las.y, las.z]).T
# 地盤点（classification=2）だけ使う
ground_pts = xyz[las.classification == 2]
normal, d, mask = ransac_plane(ground_pts, threshold=2.0, n_iter=300)
print(f'地面の法線: {normal.round(3)}')
print(f'インライア率: {mask.mean()*100:.1f}%')
山地では注意：御岳山のような起伏のある地形は1枚の平面では表せないため、インライア率が低くなる。なだらかな地形の検出か、Progressive RANSAC（検出→除去→繰り返し）で対応する。
DBSCANとの組み合わせ
| アルゴリズム | 役割 |
|---|---|
| RANSAC | 地面を検出・除去する |
| DBSCAN | 残った点（樹木など）を個別にクラスタリング |
# 地面以外の点だけDBSCANにかける
from sklearn.cluster import DBSCAN
above_ground = xyz[~mask]
labels = DBSCAN(eps=3.0, min_samples=5).fit(above_ground[:, :2]).labels_
まとめ
| 項目 | 内容 |
|---|---|
| 何をするか | ノイズ混じりの点群から平面を検出 |
| 核心 | ランダムに3点→仮平面→インライア集計→最多を採用 |
| 重要な仕上げ | インライア全体でSVD再フィット |
| パラメータ |
threshold（距離の許容幅）、n_iter（試行回数） |
| 弱点 | 曲面には向かない、計算コストが高い |
このシリーズは、東京都多摩地域の森林LiDARデータをPythonで分析する中で出会ったアルゴリズムを1つずつ解説しています。
実装コードは tama-forest notebooks で公開中。
作者のその他の実験 → https://niikun.net