ぱたへね

はてなダイアリーはrustの色分けができないのでこっちに来た

点群の法線ベクトルと主成分分析

詳解三次元点群処理に、法線ベクトルの求め方にPCA(主成分分析)を使うと書いてあり、少し驚いたので整理しました。

www.kspub.co.jp

点群における法線ベクトルとは

計測した点群のある点について、周辺の点群から面を推定しその面に対して垂直なベクトルのことです。点群の元が立方体であれば、点群の法線ベクトルを使ってどっちにどれくらい傾いているかが分かります。

PCAとは

今流行のデータサイエンスだと、多くの変数を持つデータに対して次元圧縮の手法として説明されることが多いです。次元圧縮した後、別の機械学習でクラスタリング等はよく見ます。

ある点周辺の点群に対してPCAを使うとこんなイメージになります。 点群周辺が平面になると仮定すると、その平面で楕円を作って楕円の半長軸と半短軸に相当する固有ベクトルが計算されます。次元圧縮であれば、大きな固有値に対応する2つの固有ベクトルを使えば終了です。

点群の面という観点から見ると、大きい固有値に対応する2つの固有ベクトル(緑)が平面を、最後の固有ベクトル(赤)が法線ベクトルを表します。

言われてみれば単純なのですが、仕組みが分かったときはちょっと感動してしまいました。

Python実装

せっかくなのでPythonでも実験してみました。 xy平面上にランダムに点をプロットした後、x軸に回転させて法線ベクトルを求めてます。

import numpy as np
from scipy.stats import norm
from scipy.spatial.transform import Rotation
from sklearn.decomposition import PCA
import matplotlib.pyplot as plt

# x, y平面上に点を作成
num = 20
xs = norm.rvs(size=num)
ys = norm.rvs(size=num)
zs = np.zeros(num)

# 回転行列を作る
# 45 degree
theta = np.pi / 4
r = np.array([[1, 0 , 0],
              [0, np.cos(theta), -np.sin(theta)],
              [0, np.sin(theta), np.cos(theta)]])
rot = Rotation.from_matrix(r)

# 回転させる
points = np.stack([xs, ys, zs]).T
points = rot.apply(points)

# plot
xs = points[:,0]
ys = points[:,1]
zs = points[:,2]

fig = plt.figure()
ax = fig.add_subplot(projection='3d')
ax.scatter(xs, ys, zs)
ax.set_zlim(-3, 3)
plt.xlim(-3, 3)
plt.ylim(-3, 3)
plt.show()
plt.savefig('point0.png')

# PCA
pca = PCA()
pca.fit(points)

# 固有ベクトルの表示
print(pca.components_)

一応プロットはしているが、傾いているかどうかはわかりにくいです。 三次元の散布図を書くのにzlim()が無いことに地味に苦労しました。

実行するとこのようになります。

[[ 0.96289737  0.19082538  0.19082538]
 [-0.26986784  0.68087126  0.68087126]
 [-0.          0.70710678 -0.70710678]]

最後の(0.0, 0.7, -0.7)が法線ベクトルになります。x軸に対してちょうど45°回転していますね。

上のソースのthetaを0にすると点はxy平面のみになります。 実行結果はこうなります。

[[-0.89910619 -0.43773059 -0.        ]
 [ 0.43773059 -0.89910619 -0.        ]
 [ 0.          0.          1.        ]]

(0, 0, 1)なので、ちゃんと法線ベクトルになってますね。

点群に対してPCAを行うと法線ベクトルが求められます。そのとき、次元圧縮では使わない固有値が一番少ない固有ベクトルが点群ではとても役に立つ所が面白かったです。