3次元点群の地表面抽出の方法であるCSF (Cloth Simulation Filter) について

2025年6月12日に行われたオンラインの勉強会 #StudyCO にて3次元点群の地表面抽出の方法であるCSFについて紹介しました。

動画
https://youtu.be/pVe4oUVSBYQ?t=1237

Transcript

1. 3次元点群の地表面抽出の方法である CSFについて ImVisionLabs株式会社代表取締役 板倉健太 博士（農学）
2. CSFとは バックパック型スキャナー iPhone12 LiDAR 3次元点群から地表面を抽出する方法の1つ 以下の図では、地表面が黄色、その他が青色で示されている [動画] データ出典:上図はToronto3Dデータセットを利用しています （https://github.com/WeikaiTan/Toronto-3D)
3. CSFとは バックパック型スキャナー iPhone12 LiDAR  CSFとは、点群データを上下反転させ、その上から仮想の布を落とすことで、地面と非地面を 分類する手法である。  布の全体が固定された後、各点群と布との最短距離が計算され、その距離に基づいて地面か 非地面かの分類が行われる。 [動画] 図出典:Zhang, W., Qi, J., Wan, P., Wang, H., Xie, D., Wang, X., & Yan, G. (2016). An easy-to-use airborne LiDAR data filtering method based on cloth simulation. Remote Sensing, 8(6), 501. https://doi.org/10.3390/rs8060501
4. CSFの実行方法 iPhone12 LiDAR [動画]  CloudCompareの他にPDALのfilters.csfやMATLABのFile Exchangeのコードを利用可能  CloudCompareにて Plugins > CSF Filter の順に進むことでCSFを実行できる MATLAB版: https://jp.mathworks.com/matlabcentral/fileexchange/58139- csf-cloth-simulation-filter?s_tid=FX_rc3_behav
5. CSFの実行方法 iPhone12 LiDAR [動画]  本スライドでは主に以下のコードを利用してCSFを実行 MATLAB版: https://jp.mathworks.com/matlabcentral/fileexchange/58139- csf-cloth-simulation-filter?s_tid=FX_rc3_behav
6. パラメーターについて バックパック型スキャナー iPhone12 LiDAR  CSFアルゴリズムにおける主要なパラメーターは、以下の表に示す通りである [動画] パラメータ デフォルト値 説明 Rigidity Index(RI) 1.0 布の剛性：1 = 傾斜地, 2 = 緩 やかな斜面, 3 = 市街地の平 坦地。 Smooth Terrain(ST) False イテレーション終了後の後処 理が行われるかのフラグ (TrueまたはFalse) Grid resolution(GR) 0.5 地表を模擬するクロス（布）の グリッド間隔（単位: m）。大き くすると計算が速くなるが精度 が下がる。 デフォルト値はCloudCompareやMATLAB実装を参考にしています
7. パラメーターについて バックパック型スキャナー iPhone12 LiDAR  CSFアルゴリズムにおける主要なパラメーターは、以下の表に示す通りである [動画] パラメータ デフォルト値 説明 Height threshold for Cloud- to-Cloth distance(hcc ) 0.5 点が地面か非地面かを分け る高さ差のしきい値（単位: m）。 大きすぎると樹木や構造物も 地面と誤判定されやすい。 Maximum number of iterations 500 布が落下して安定するまでの 最大反復回数。 Duration of time(dT) 0.65 クロスシミュレーションでの時 間ステップ（動きの滑らかさに 関係）。単位はs。 Height threshold for Cloth Post-processing(hcp ) 0.3 布と点群の間の高さ差がこれ 以下なら、その点は布（地面） に吸着される。単位はm。 デフォルト値はCloudCompareやMATLAB実装を参考にしています
8. パラメーターについて バックパック型スキャナー iPhone12 LiDAR  CSFの主要な実装ファイルである csf_filtering.cpp 内の関数 csf_filtering にて以下のように パラメータが定義されている  関数の前半7つの引数のうち、1つ目を除くものをパラメータ調整可能 [動画] 参考ページ： https://github.com/jianboqi/CSF/blob/master/matlab/csf_filtering.cpp
9. イテレーション内における物理挙動のイメージ バックパック型スキャナー iPhone12 LiDAR [動画]  以下のように上下を逆さまにした状態から布を落としていくイメージ 山 山 人工構造物 天地を逆にしている
10. イテレーション内における物理挙動のイメージ バックパック型スキャナー iPhone12 LiDAR [動画]  前スライドでは平面的なアニメーションを示したが、実際は以下の図のように、2次元のグリッ ドにてシミュレーションを行う  布を落としていくイメージであるが、布は質量を持つ粒子とそれらをつなぐ接続からなるグリッ ドとしてモデル化され、これを質点・バネモデル（Mass-Spring Model）と呼ぶ 図出典:CSFの論文のFigure2を引用しています
11. イテレーション内における物理挙動のイメージ バックパック型スキャナー iPhone12 LiDAR [動画]  以下の図はCloudCompareにてCSFを実行し、布をメッシュモデルとしてエキスポートした様  この緑のメッシュを基準面として十分近い点を地表面であると分類する
12. イテレーション内における物理挙動のイメージ iPhone12 LiDAR [動画]  布の粒子の間隔をCloth Resolutionにて定義する  CloudCompareにおいてはAdvanced Settingにて変更可能 和訳：1 布の解像度（Cloth resolution）とは、地形を覆うために使 用される布グリッドのサイズを指す（単位は点群と同じ）。この 値が大きいほど、得られる数値地形モデル（DTM）は粗くなる。
13. イテレーション内における物理挙動 バックパック型スキャナー iPhone12 LiDAR [動画] ①重力落下 • 接触判定されていない点に対してのみ適用される • 以下の式に基づいて落下する • Xおよびmは布を構成する粒子の位置や質量を表し、Δt は時間ステップを表す • 時間ステップと初期位置が与えられれば、G が定数であるため現在の位置を 直接計算できる。
14. イテレーション内における物理挙動 バックパック型スキャナー [動画]  重力落下の式の導出 現在位置 X (t)の、時刻 t に対するテイラー展開は以下のように表すことができる。 n+1次以降の項をまとめた「誤差項」 この式の意味は次の通り： 未来の時刻における位置 速度（1階導関数） 加速度（2階導関数） n階導関数（n次の変化率） n+1次以降の項をまとめた「誤差項」
15. イテレーション内における物理挙動 iPhone12 LiDAR [動画]  重力落下の式の導出 この時刻 tに対するテイラー展開から、3次以上の高次項を無視して整理すると以下のようになる。
16. イテレーション内における物理挙動 iPhone12 LiDAR [動画]  重力落下の式の導出 両辺を加算して 整理して
17. イテレーション内における物理挙動 バックパック型スキャナー iPhone12 LiDAR [動画]  重力落下の式の導出 ሷ X (t) に代入して整理すると ニュートンの第二法則より、
18. イテレーション内における物理挙動 バックパック型スキャナー [動画] ②点の接触判定 • 布上の点が点群と接触判定された場合、その点はそれ以降動かない • 接触判定の条件は以下の式で表される 図出典:CSFの論文のFigure6を引用しています 水色: 布の点が点群と接触しているためこれ以 上下へは動かない
19. イテレーション内における物理挙動 バックパック型スキャナー [動画] ②点の接触判定 • 布上の点が点群と接触判定された場合、その点はそれ以降動かない • 接触判定の条件は以下の式で表される 図出典:CSFの論文のFigure6を引用しています 赤色: 布の点が点群と接触しておらず、次の試行 （イテレーション）でも移動する可能性がある
20. イテレーション内における物理挙動 iPhone12 LiDAR [動画] ②点の接触判定 • 布上の点が点群と接触判定された場合、その点はそれ以降動かない • 接触判定の条件は以下の式で表される 落下している点の高さを z c 、地表面の高さを z t として、以下の式を満たした場合、落下している点 は固定される。 h cp は高さ差のしきい値（Height threshold for Cloth Post-processing）を表している。
21. イテレーション内における物理挙動 バックパック型スキャナー iPhone12 LiDAR [動画] ③剛性による調整 重力落下、点の接触判定が終わると、布を構成する各点は剛性による調整が行われる。 調整をする点をp 0 、隣接する点を p i とすると、 p 0 に対しては以下の変異ベクトルdが適用される。 b は接触状態に基づくブール変数であり、点群と接触している場合はb=0となって d=0となる。接触してい ない場合はb=1である。 𝑛はz軸方向の垂直ベクトル (0,0,1) を表すことから、 ( Ԧ 𝑝 i - Ԧ 𝑝 0 ) ・𝑛は点 റ 𝑝 i が റ 𝑝 0 に対してどれだけ垂直方向 （高さ方向）に離れているかを表す。
22. イテレーション内における物理挙動 バックパック型スキャナー iPhone12 LiDAR [動画] ③剛性による調整 • この変位処理は、布の剛性を表すパラメータ RIの値に応じて複数回繰り返さ れる。
23. イテレーション内における物理挙動 バックパック型スキャナー iPhone12 LiDAR [動画] • CSFの論文では以下の図にて説明されている（Figure5） 図出典:CSFの論文のFigure5を引用しています
24. イテレーション内における物理挙動 iPhone12 LiDAR [動画] • 以下の式についての補足 • nは垂直成分の (0,0,1)であるため、2つの隣り合う粒子の位置ベクトルの差と の内積を利用することでz成分の値のみを取り出すことができる
25. イテレーション内における物理挙動 [動画] • これまでのスライドでのルールに従い、回数（イテレーション）を増やしていくこ とで、次第に布が地表面へフィットする • 以降のスライドにて、その様子をアニメーションにて示す • イテレーションについてはCloudCompareにおいてはAdvanced Settingにて変更可能 和訳：2. 最大反復回数（Max iterations）とは、地形シミュレー ションにおける反復計算の最大回数を指す。一般的なシーンで は、500回で十分である。
26. イテレーション100での布の落下過程 [a] ScanX2.0
27. イテレーション100での布の落下結果 [a] ScanX2.0 • 左右の2つのこぶが山であり、中央上部の点群が 建物である • 地面・非地面の区別が適切に行われていた
28. イテレーション80での布の落下過程 [a] ScanX2.0
29. イテレーション80での布の落下結果 [a] ScanX2.0 • イテレーションが80だと、布が十分に落ちるのに 時間が足りなかった • そのため、地面の多くが非地面に区別されていた
30. イテレーション120での布の落下過程 [a] ScanX2.0
31. イテレーション120での布の落下結果 [a] ScanX2.0 • イテレーションが120だと、布が地形に沈みすぎて いた • そのため、建物の多くが地面に区別されていた
32. 地面・非地面の判定について [a] ScanX2.0  後処理による点の位置調整 • この補正処理は、パラメータSmooth TerrainがTrueに設定されている場合に適用される。 • 特に、崖や急傾斜では、布が地表面に接触せず、非地面と誤判定されることがある。 この問題に対処するために、布を構成する一部の未接地点に次の調整を行う
33. 地面・非地面の判定について [a] ScanX2.0  後処理による点の位置調整 ①まだ接地していない点Aに対して、その4つの近傍点に接地済みの点Bが存在する かを確認する。点が見つからなければ調整は行われない。
34. 地面・非地面の判定について [a] ScanX2.0  後処理による点の位置調整 ②見つかった点Bと、点Aに対して、それぞれの直下の地形点を取得する。 ここで、点Aの直下の地形点を点Cとし、点Bの直下の地形点を点Dとする。 点Cと点Dの高さの差がℎ𝑐𝑝 の固定値である0.3mより小さい場合に処理は進 む。
35. 地面・非地面の判定について [a] ScanX2.0  後処理による点の位置調整 ③点Aの高さを調整し、点Cの高さと同じにする。この補正により、未接地点だった点A は接地済の点となる。その結果、地面として正しく分類されやすくなる。
36. 地面・非地面の判定について [a] ScanX2.0  後処理の伝搬について • 補正されて接地した点に隣接する未接地点へ補正が伝搬する。 • 接地済の点に隣接する未接地点群をSCC（強連結成分）とし、BFS（幅優先探索）に より外縁から内側へ順に補正を伝搬させる
37. 地面・非地面の判定について [a] ScanX2.0  点の地面・非地面の判定 • 点の接地・未接地の判定は以下の式によって行われる
38. 地面・非地面の判定について [a] ScanX2.0  hcc についてはCloudCompareにおいてはAdvanced Settingにて変更可能 ※hcc についての和訳： 分類しきい値（Classification threshold）とは、シミュレーション された地表面と各点群との距離に基づいて、点群を地表 （ground）と非地表（non-ground）に分類するためのしきい値で ある（単位は点群と同じ単位）。一般的なシーンにおいては、0.5 が適用されることが多い。
39. CSF適用後の地面分類例 [a] ScanX2.0  パラメーター設定について パラメータ 値 Rigidity Index(RI) 3.0(例2のみ1.0) Smooth Terrain(ST) True Grid resolution(GR) 1.0 Height threshold for Cloud-to-Cloth distance(ℎ𝑐𝑐) 以下の例によって異なる Maximum number of iterations 500 Duration of time(dT) 0.65 Height threshold for Cloth Post-processing(ℎ𝑐𝑝) 0.3 • 今回の分類は以下のようなパラメーター設定で行った
40. CSF適用後の地面分類例 [a] ScanX2.0  例1:hcc =0.5の場合 • 一般的な場合の結果で、 正しく分類されていた
41. CSF適用後の地面分類例 [a] ScanX2.0  例2:hcc=0.9,RI=1の場合 • 剛性が小さく、地面と判定される閾 値が大きいため、本来なら地面が多 くなるはずである。 • しかし今回の設定では、閾値が相対 的に小さかったため、地面として分 類される範囲はそれほど広がらな かった。
42. CSF適用後の地面分類例 [a] ScanX2.0  例3:hcc=0.1の場合 • 地面と判定する閾値を0.1と小さく設 定したことで、布との距離が0.1m以 内の点しか地面と判定されなくなる。 • これによって非地面として分類され る点の範囲が広がった。
43. CSF適用後の地面分類例 [a] ScanX2.0  極端な例1:hcc =0の場合 • 閾値を極端に小さく（たとえば0m）設 定すると、布との距離が完全に0で なければ地面と判定されない。 • その結果、すべての点が非地面と分 類されてしまう。
44. CSF適用後の地面分類例 [a] ScanX2.0  極端な例2:hcc =9の場合 • 閾値を9mと大きく設定することで、 布から大きく離れた点でも地面とみ なされる可能性が高くなる。 • これにより、実際には構造物や樹木 である点も、誤って地面として分類さ れるリスクが増加する。