点群からの法線の導出

はじめに

最近、法線ベクトルを求める機会があった。Point Cloud Library（PCL）が提供する関数を使えば簡単に求めることができる。法線ベクトルを求めるアルゴリズムについてもここに説明があるが、結果の式しか記載されていない。今回はこの式の導出を行う。

問題設定

いま、3次元空間内に存在する点群を考える。この点群はなんらかの物体の表面から取得されたものである。点群だけから物体の表面を推測し、表面上の各点における法線ベクトルを求めたい。

方針

点群から法線を求める方法には2つある。

  1. 点群から曲面を求め、その曲面から各点の法線を計算する。
  2. 点群から直接法線を計算する。

前者では、点を結ぶ小さな三角形をたくさん作る（ポリゴン化）。ノイズの載った点が混じると凸凹した表面が出来上がるため、ノイズに弱い手法である。一方、後者では、対象点近傍の点を複数個考慮し法線を求めるため、ノイズの効果は相対的に小さくなる（後述）。PCLが採用している手法は後者であり、その手順は以下の通りである。

  1. ある点Aを考える。
  2. 点Aの近傍K個の点を考える。
  3. 点Aと近傍K個の点を使い、最小自乗法により平面を決定する。
  4. その平面の法線ベクトルを求める。

以下この方針に従って法線ベクトルを求める。

導出

点A近傍のK個の点の集合をP=(\vec{p}_1,\cdots,\vec{p}_K)と書くことにする。点群Pが近似する平面Bの法線を\vec{n}と書くと、平面B上の点\vec{x}は次式を満たす。

(1)   \begin{equation*} \vec{n}^T \vec{x}=d \end{equation*}

ここで、dは定数である。点群Pが平面Bを近似するには次を最小にすれば良い。

(2)   \begin{equation*} F(\vec{n},d)=\sum_{k=1}^{K}\left(\vec{n}^T \vec{p}_k-d\right)^2 \end{equation*}

上式を最小にするような\vec{n},dを求める。平面の式(1)には定数倍だけの不定性がある。これを取り除くため、\vec{n}の大きさを1に限定する条件を、Lagrangeの未定係数法により導入する。

(3)   \begin{equation*} F(\vec{n},d,\lambda)=\sum_{k=1}^{K}\left(\vec{n}^T \vec{p}_k-d\right)^2-\lambda(\|\vec{n}\|^2-1) \end{equation*}

ここで、\lambdaもまたFを最小化するように決定すべきパラメータである。最初に、Fをdで偏微分し0とおくことにより次式を得る。

(4)   \begin{equation*} \sum_{k=1}^K(\vec{n}^T\vec{p}_k-d)=0 \end{equation*}

これより次式を得る。

(5)   \begin{eqnarray*} \vec{n}^T\vec{p}_a&=&d \\ \vec{p}_a&=&\frac{1}{N}\sum_{k=1}^K \vec{p}_k \end{eqnarray*}

\vec{p}_aは点群Pの重心である。上式を式(3)に代入すると

(6)   \begin{equation*} F(\vec{n},\lambda)=\sum_{k=1}^{K}\left[\vec{n}^T(\vec{p}_k-\vec{p}_a)\right]^2-\lambda(\|\vec{n}\|^2-1) \end{equation*}

を得る。これを今度は\vec{n}で偏微分し0とおくことにより次式を得る。

(7)   \begin{equation*} \left[\vec{n}^T\sum_{k=1}^K\left(\vec{p}_k-\vec{p}_a\right)\right]\left(\vec{p}_k-\vec{p}_a\right)=\lambda\vec{n} \end{equation*}

上の式は少し分かり難いので説明しておく。[]で囲われた部分は内積であり、スカラー量となる。従って上式左辺は全体としてベクトル量になる。また、添字kについての和は内積部分だけで閉じておらず、後ろの\left(\vec{p}_k-\vec{p}_a\right)にもかかっていることに注意する。さて、式(7)は以下のように変形できる。

(8)   \begin{eqnarray*} Q\vec{n}&=&\lambda\vec{n}\\ Q&=&\sum_{k=1}^K\left(\vec{p}_k-\vec{p}_a\right)\left(\vec{p}_k-\vec{p}_a\right)^T \end{eqnarray*}

この式は、点群Pの共分散行列Qの固有値を求める式である。つまり、法線の導出問題が共分散行列Qの固有値問題に帰着したことになる。これが、ここに記載されている式である。いま考えている空間は3次元であるからQは3\times 3の行列である。さて、式(2)に式(5)を代入すると

(9)   \begin{equation*} F(\vec{n})=\sum_{k=1}^{K}\left[\vec{n}^T\left(\vec{p}_k-\vec{p}_a\right)\right]^2 \end{equation*}

を得る。Qを使って上式を書き直すと

(10)   \begin{eqnarray*} F(\vec{n})&=&\vec{n}^T Q\vec{n} \\ &=&\lambda\vec{n}^T\vec{n} \\ &=&\lambda \end{eqnarray*}

となる。ここで、\|\vec{n}\|=1を用いた。つまり、F(\vec{n})を最小にするには、複数の固有値のうち最小のものを選択すれば良いことになる（いまの場合、非ゼロの固有値は最大で3つ求まる）。その固有値に対応する固有ベクトルが求めたい法線ベクトルである。

主成分分析との関係

式(8)は、点群Pの位置座標について主成分分析を行う際に出てくる式である。得られる3つの主成分のうち一番小さい寄与を持つ方向が法線に対応するわけである。

まとめ

今回はPoint Cloud Library（PCL）が提供する法線ベクトル導出アルゴリズムを定式化した。PCLを用いた法線導出のサンプルコードはここを見てほしい。ここでは適用結果だけを示す。

法線\vec{n}の3成分をRGBに対応させた画像である。今回用いた3Dモデルの場合、近傍点の数Kを20個程度にすると上のような綺麗な法線,法線ベクトル,点群（法線が滑らかに変化する）画像となった。

Kumada Seiya

Kumada Seiya

仕事であろうとなかろうと勉強し続ける、その結果”中身”を知ったエンジニアになれる

最近の記事

  • 関連記事
  • おすすめ記事
  • 特集記事

アーカイブ

カテゴリー

PAGE TOP