ごあいさつ

こんにちは！建設業でSEをしています、Seikyoと申します。
今回で2回目の記事投稿となります。
プログラムというより数学色強めで記事を書いていきます

前回の振り返り

前回の記事では、3次元空間上の平面の上下左右を区別させる方法をご紹介しました。

▼ 前回の記事
[Unity] 3次元空間で定義した任意の平面を自由に扱う方法

前回の記事に対して、こんな声が聞こえてきそうです。

「いやいや、そんな理想的な平面なんて、現実世界には存在しないでしょうよ」

まさにその通りです。現実で相手にする面なんか全部デコボコしてますよね。
では、そんな現実世界相手にどのように我々は戦っていきましょうか。

今回も、数学的な問いを立てて、解決策を考えていきましょう。

Q1. 3次元空間上に散らばった N 個の点群から、それらの「平均的な面」を求めるには？

現実の壁を3Dスキャンしたとき、返ってくるのは滑らかな平面ではなく散らばった点群。
この点群を最もよく代表する1枚の平面が取り出せれば、あとは前回の方法にそのまま接続できる。

方針はシンプルで、各点と面との距離（の二乗和）が最小となるような面を求めればよい。
（座標系での距離計算なので、+と-で打消しあってしまうのを防ぎたい）

記法について:
ベクトルはすべて3次元の縦ベクトルとし、矢印表記は省略して $n, c, p_i, q_i$ と書く。「$\cdot$」は内積、右肩の $T$ は転置を表している

平面の定義

この平面を、各点の重心 $c$ を通り、法線 $n$ を持つ面と定義する。各点を $p_i$ とする。

p_i = (x_i,\, y_i,\, z_i)

法線は単位ベクトルとする。

|n| = 1

重心 $c$ は次式。

c = \frac{1}{N}\sum_{i=1}^{N} p_i

点と面の距離

点 $p_i$ から面までの符号付き距離を $d_i$ とすると、

d_i = n \cdot (p_i - c)

補足:
内積の定義より、 $n \cdot (p_i - c) = |p_i - c||n|\cos\theta = |p_i - c|\cos\theta$（$\because |n| = 1$）
重心から各点へ向かうベクトル $p_i - c$ を単位法線 $n$ の方向へ投影した長さであり、これが点と平面の符号付き距離そのものになっている。

誤差関数の定義

$ d_i$の二乗和が最小であれば良かったので、
$\sum_{i=1}^{N} d_i^2$ が最小となればよい。これを誤差関数 $E(c, n)$ とすると、

E(c, n) = \sum_{i=1}^{N} d_i^2 = \sum_{i=1}^{N} \bigl(n \cdot (p_i - c)\bigr)^2

未知数は通過点 $c$ と単位法線 $n$ の2つ。ここから1つずつ潰していく。

2. 最適な平面は必ず重心を通るはずという直感の確認

まず $c$ から。誤差関数 $E$ を $c$ で偏微分し、$0$ となる条件（停留条件）を求める。
合成関数の微分（${f(x)^2}' = 2f(x) \cdot f'(x)$、後で検算する）より、

\frac{\partial E}{\partial c} = \sum_{i=1}^{N} 2\bigl(n \cdot (p_i - c)\bigr) \cdot (-n) = -2\sum_{i=1}^{N} \bigl(n \cdot (p_i - c)\bigr)\, n

これを $0$ とおく。

0 = -2\sum_{i=1}^{N} \bigl(n \cdot (p_i - c)\bigr)\, n

和の中で $n$ は共通なので括り出せて、

-2n\left(n \cdot \sum_{i=1}^{N} (p_i - c)\right) = 0

$n \neq 0$ だから、これが任意の $n$ について成り立つためには、

\sum_{i=1}^{N} (p_i - c) = 0

和を展開すると、

\sum_{i=1}^{N} p_i - Nc = 0

Nc = \sum_{i=1}^{N} p_i

c = \frac{1}{N}\sum_{i=1}^{N} p_i

結論: $n$ がどんな向きであれ、その向きのまま誤差を最小にしたければ、とにかく点群の重心を通るように面を置くのがベストである。

これで未知数が1つ減った。以降は $c$ を重心に固定して、法線 $n$ だけを追いかける。

ちょっと復習
偏微分して0になる値を探す行為は、他の変数を固定してある変数だけ変化させたときに、値が変化しないことを示す。
今回だと「法線ベクトル ($n$) はその向きで固定して動かさない」と決めた上で、面の置き場所（位置 ($c$)）だけを動かしたときに、最も誤差 $E(c,n)$が変化しなくなる（これ以上減らなくなる底の）状態を探す行為。

3. 固有値問題への帰着

転置行列で内積を書き換える

各点の重心からの偏差は、 $p_i - c$ と書ける。これを$q_i$と置くと,

E(c, n) = \sum_{i=1}^{N} d_i^2 = \sum_{i=1}^{N} \bigl(n \cdot (p_i - c)\bigr)^2

の、$Σ$に内包されている内積について、転置行列を用いて

n \cdot (p_i - c) = n \cdot q_i = n^T q_i

と書ける

二次形式への変形

これを使って距離の二乗和を変形していく。
ポイントは、スカラーの転置はそれ自身に等しいこと。$((n^T q_i)^T = q_i^T n)$
(内積の計算結果がスカラーなこと忘れがちなの私だけですかね...)

\sum_{i=1}^{N} d_i^2 = \sum_{i=1}^{N} \bigl(n^T q_i\bigr)^2 = \sum_{i=1}^{N} \bigl(n^T q_i\bigr)\bigl(n^T q_i\bigr)

= \sum_{i=1}^{N} \bigl(n^T q_i\bigr)\bigl(q_i^T n\bigr) = \sum_{i=1}^{N} n^T \bigl(q_i\, q_i^T\bigr)\, n

= n^T \left(\sum_{i=1}^{N} q_i\, q_i^T\right) n

括弧内の $3 \times 3$ 対称行列を $M$ とおく（今度は外積の話ですな）。

M = \sum_{i=1}^{N} q_i\, q_i^T

すると誤差関数は、法線 $n$ についての二次形式として簡潔に表される。

\sum_{i=1}^{N} d_i^2 = n^T M n = E(c, n)

ちょっと補足:
この $M$ は、データ数 $N$ で割ると「分散共分散行列」になります。今回は最小化する「向き（固有ベクトル）」だけが知りたいので、定数倍の $1/N$ は省略して、計算がシンプルな総和の形で進めています。

制約条件とラグランジュの未定乗数法

ここで、$|n| = 1$ より、

n^T n = 1 \quad\Rightarrow\quad n^T n - 1 = 0

制約条件付きの最小化問題になったので、ラグランジュの未定乗数法を使う。ラグランジュ関数 $L$ を次のように定義する。

L = n^T M n - \lambda\,(n^T n - 1)

ラグランジュ乗数法の直観:
第1項 $n^T M n$ は目的関数 $E(c, n)$ で、これを$n$で微分すれば勾配がわかる。
第2項の括弧内は制約条件。（ $n^T n = 1$）
最適解の場所では、制約条件の曲線と目的関数の等高線が接し、両者の勾配ベクトルが平行になると考える。
①向きが等しく ②長さは定数倍（その倍率が $\lambda$）。これを1本の式にまとめたものが $\frac{\partial}{\partial n}{E - \lambda(n^T n - 1)} = 0$、つまり $\frac{\partial L}{\partial n} = 0$ 。

4. 固有値問題を解く　∂L/∂n = 0

ラグランジュ関数 $L$ を $n$ で偏微分する（ベクトル微分の公式は5章で検算する）。

\frac{\partial L}{\partial n} = \frac{\partial}{\partial n}\Bigl(E - \lambda\,(n^T n - 1)\Bigr) = 2Mn - 2\lambda n

これを $0$ とおくと、

2Mn - 2\lambda n = 0

Mn = \lambda n

出てきたこの式、まさに行列 $M$ の固有値問題である（$\lambda$: 固有値、$n$: 固有ベクトル）。

誤差の二乗和は固有値そのもの

このとき誤差関数の値を計算すると、驚くほどきれいな結果になる。

E(c, n) = n^T M n = n^T \lambda n = \lambda\,(n^T n) = \lambda

n^T n = 1 より、E(c,n)=\lambda

結論:
点群データと最も誤差が少なくなる $n$ の条件は、$E(c, n)$ が最小のとき、つまり $\lambda$ が最小のとき。
求める法線 $n$ は「共分散行列 $M$ の最小固有値に対応する単位固有ベクトル」

難しいこと書いたけど、$M$をがいくつか求まるからその中の最小値が解だよってこと

固有値の求め方

固有値の求め方も確認しておく。単位行列を $I$ とすると、

Mn = \lambda I n

Mn - \lambda I n = 0

(M - \lambda I)\, n = 0

$|n| = 1$ より $n \neq 0$ だから、この同次方程式が非自明な解を持つためには係数行列が正則であってはならない。すなわち、

\det(M - \lambda I) = 0

これが固有方程式（特性方程式）。$3 \times 3$ 行列なら $\lambda$ の解の候補は最大3個。その中から最小のものを選べばよい、というのが今回の問いの最終的な答え。

5. 検算：使った微分公式を成分計算で確かめる

途中で使ったベクトル微分の公式は教科書に載っている。
が、2次元・2×2行列で愚直に成分計算して公式を確認する。忘れやすいからね。

5.1 ∂(nᵀMn)/∂n = 2Mn の確認

$M$ を $2 \times 2$ の対称行列とする。

M = \begin{pmatrix} a & b \\ b & c \end{pmatrix}

二次形式を展開する。

n^T M n = \begin{pmatrix} x & y \end{pmatrix} \begin{pmatrix} a & b \\ b & c \end{pmatrix} \begin{pmatrix} x \\ y \end{pmatrix} = \begin{pmatrix} x & y \end{pmatrix} \begin{pmatrix} ax + by \\ bx + cy \end{pmatrix}

= ax^2 + bxy + bxy + cy^2 = ax^2 + 2bxy + cy^2

各成分で偏微分すると、

\frac{\partial}{\partial x}\,(ax^2 + 2bxy + cy^2) = 2ax + 2by

\frac{\partial}{\partial y}\,(ax^2 + 2bxy + cy^2) = 2bx + 2cy

ベクトルにまとめると、

\begin{pmatrix} 2ax + 2by \\ 2bx + 2cy \end{pmatrix} = 2\begin{pmatrix} a & b \\ b & c \end{pmatrix}\begin{pmatrix} x \\ y \end{pmatrix} = 2Mn

よって、$M$ が対称行列のとき $\frac{\partial (n^T M n)}{\partial n} = 2Mn$ が確認できた。

注意:
$M$ が非対称の場合、一般には $\frac{\partial (n^T M n)}{\partial n} = (M + M^T),n$ になる。
今回の共分散行列 $M = \sum q_i q_i^T$ は構成上必ず対称なので $2Mn$ で問題ない。

5.2 ∂(nᵀn)/∂n = 2n の確認

$n$ を2次元ベクトル $n = (x, y)^T$ として成分で確認する。

n^T n = \begin{pmatrix} x & y \end{pmatrix} \begin{pmatrix} x \\ y \end{pmatrix} = x^2 + y^2

各成分で偏微分すると（$\lambda$ 倍付きで確認）、

\frac{\partial}{\partial x}\,\lambda(x^2 + y^2) = 2\lambda x, \qquad \frac{\partial}{\partial y}\,\lambda(x^2 + y^2) = 2\lambda y

ベクトルにまとめると、

\begin{pmatrix} 2\lambda x \\ 2\lambda y \end{pmatrix} = 2\lambda \begin{pmatrix} x \\ y \end{pmatrix} = 2\lambda n

4章の $\frac{\partial L}{\partial n}$ はこの2つの公式の組み合わせでできるので是非お手元の紙で。

まとめ

ステップ	結果
問題設定	距離の二乗和 $E(c, n)$ を最小化する平面を探す
$c$ の最適化	最適な平面は必ず重心を通る
$n$ の最適化	共分散行列 $M$ の最小固有値の固有ベクトルが法線
誤差の意味	二乗和の最小値は最小固有値 $\lambda$ そのもの

これによって出来上がった平面を用いれば、
前回のプログラムをなんと現実世界から観測したデータに対して使用することができます。

次回は、せっかく自由に使える平面があるのでちょっと遊んでみましょうか。
UnityC#での実装方法も載せようかな。
今時AIにこの記事読ませたら書いてくれそうなもんだけど。

今回、数学的な記述のみでコードを書かなかった理由は、
全ての業界の人間の共通言語は数学である
と考えたためです。
ITに強い人間、現場仕事に強い人間、会社の偉い人たち、機械メーカー、学生、などなど
みんなとりあえずコードは書けなくても数学的に解釈できれば、
業界を絞らない発展に繋がるんじゃないかと考えたためです。
みんなで一緒に世界のサイエンスを盛り上げていきましょう！

ではでは、皆様！ご安全に！