Estimating Surface Normals in a PointCloud

ポイントクラウドにおける表面法線の推定

出典: http://pointclouds.org/documentation/tutorials/normal_estimation.php#normal-estimation

表面法線は、幾何学的表面の重要な特性であり、CGアプリケーションなどの多くの領域で、シェーディングやその他の視覚効果を生成する正しい光源を適用するために頻繁に使用される。

幾何的表面が与えられれば、表面上のある点での法線の方向を、その点を起点とする表面に垂直なベクトルとして推測するのは通常易しいことである。 しかし、ポイントクラウドのデータセットは実際の表面上の点のサンプル集合であるため、次の2つの可能性がある。

* ポイントクラウドのデータセットからその下にある表面を表面メッシュを使用して取得し、メッシュから表面法線を計算する。
* 近似を使用してポイントクラウドのデータセットから表面法線を直接推論する。

このチュートリアルでは2番目を採用、つまりポイントクラウドのデータセットに対し、クラウド内の各点で表面法線を直接計算する。

https://www.youtube.com/embed/x1FSssJrfik

基本理論

さまざまな通常の推定方法が存在するが、このチュートリアルでは、最も簡単な方法の1つで、次のように定式化される方法に焦点をあてる。 表面上の点の法線を決定する問題は、表面に接する平面の法線を推定する問題として近似する、つまり、最小二乗平面あてはめ推定問題とする。

注意: 最小二乗問題の数式を含む詳細については、Rusu (2009) PhD Thesis http://mediatum.ub.tum.de/download/800632/800632.pdf を参照。

表面法線を推定する方法は、対象点の最近傍から作成される共分散行列の固有ベクトルと固有値(またはPCA - 主成分分析)の分析に還元される。 より具体的には、各点 $\boldsymbol{p}_i$ に対して、共分散行列 $\mathcal{C}$ を次のように組み立てる: $\displaystyle \mathcal{C} = \frac{1}{k} \sum^k_{i=1} (\boldsymbol{p}_i - \bar{\boldsymbol{p}_i})\cdot(\boldsymbol{p}_i - \bar{\boldsymbol{p}_i})^T$, $\mathcal{C}\cdot\vec{\bf v}_j = \lambda_j \cdot \vec{\bf v}_j$, $j \in \{0,1,2\}$

ここで k は$\boldsymbol{p}_i$の近傍で考慮される点近傍の数であり、$\overline {\boldsymbol{p}}$は最近傍の3D重心を表し、$\lambda_j$は共分散行列の$j$番目の固有値 、および$ \vec{{\mathsf v}_j}$は $j$番目の固有ベクトル。

点集合から共分散行列を推定するには、次のものを使用できる。

In [2]:
from py3d import *

print("Testing IO for point cloud ...")
pcd = read_point_cloud("./TestData/fragment.pcd")
Testing IO for point cloud ...
In [8]:
# まずnumpy配列にする
import numpy as np
print(pcd.has_points())
pcd_array=np.asarray(pcd.points)
True
In [30]:
pcd_array.shape
pcd_array_reshape= pcd_array.T
print(pcd_array[:5])
print(pcd_array_reshape[:,:5])
[[1.16796875 1.01803279 0.96484375]
 [1.16845131 1.01953125 0.96484375]
 [1.16796875 1.02158833 0.96484375]
 [1.16796875 1.01953125 0.96634495]
 [1.16796875 1.03289878 0.95703125]]
[[1.16796875 1.16845131 1.16796875 1.16796875 1.16796875]
 [1.01803279 1.01953125 1.02158833 1.01953125 1.03289878]
 [0.96484375 0.96484375 0.96484375 0.96634495 0.95703125]]
In [31]:
# 距離画像においてpcdの重心を求める
cx,cy,cz=[np.mean(x) for x in pcd_array_reshape]
print(cx, cy, cz)
2.2721989223271186 1.7389421708490782 1.3153000239180128
In [32]:
# 共分散行列
np.cov(pcd_array_reshape)
Out[32]:
array([[ 0.49565598, -0.0380351 ,  0.04354783],
       [-0.0380351 ,  0.33533837,  0.01938292],
       [ 0.04354783,  0.01938292,  0.01870923]])

一般に、法線の符号(向き)を決定するための数学的方法はないので、上記のような主成分分析(PCA)によって計算された向きは曖昧であり、ポイントクラウドデータセット全体にわたる一貫した方向は決定できない。 下の図は、キッチン環境の一部を表すより大きなデータセットの2つのセクションに対し、この効果を示すものである。 図の右側には、点球からのすべての法線の方向を表す、標準球とも呼ばれる拡張ガウス画像(EGI)が表示されている。 このデータセットは2.5次元であり、単一の視点から取得されているので、法線はこのEGIの球の半分にのみ存在しなければならない。 しかし、方向の非一貫性により、法線は球全体に広がっている。

視点${\mathsf v}_p$が実際にわかっていれば、この問題の解決は簡単である。視点に向かって一貫してすべての法線$\vec{\boldsymbol{n}}_ i$を向けるには、次の方程式を満たす必要がある:

$\vec{\boldsymbol{n}}_i \cdot({\mathsf v}_p - \boldsymbol{p}_i)> 0$

下図は、上図からのデータセット内のすべての法線が、視点に向かって一貫して向くようにした後の結果を示す。

方向が重要な場合は、orient_normals_to_align_with_directionorient_normals_towards_camera_locationなどの方向付け関数を呼び出す必要がある。

In [33]:
help(orient_normals_towards_camera_location)
Help on built-in function orient_normals_towards_camera_location in module py3d:

orient_normals_towards_camera_location(...) method of builtins.PyCapsule instance
    orient_normals_towards_camera_location(cloud: py3d.PointCloud, camera_location: numpy.ndarray[float64[3, 1]]=array([0., 0., 0.])) -> bool
    
    Function to orient the normals of a point cloud

In [34]:
help(orient_normals_to_align_with_direction)
Help on built-in function orient_normals_to_align_with_direction in module py3d:

orient_normals_to_align_with_direction(...) method of builtins.PyCapsule instance
    orient_normals_to_align_with_direction(cloud: py3d.PointCloud, orientation_reference: numpy.ndarray[float64[3, 1]]=array([0., 0., 1.])) -> bool
    
    Function to orient the normals of a point cloud

In [36]:
pcd.has_normals()
Out[36]:
True
In [40]:
orient_normals_towards_camera_location(pcd,np.array([0., 0., 0.]))
Out[40]:
True
In [42]:
orient_normals_to_align_with_direction(pcd,np.array([0., 0., 1.]))
Out[42]:
True

Selecting the right scale

正しいスケールの選択

前に説明したように、点の表面法線は、その点の周囲の点近傍サポート(k-近傍 ともいう)から推定する必要がある。

最近傍推定問題を詰めていくと、正しくスケールを推定する、つまりサンプリングされたポイントクラウドのデータセットが与えられた場合、対象の点の最近傍を決定するのに用いるkrの正しい値がどのようなものかという問題が起きる。

この問題は非常に重要であり、点特徴表現の自動推定(すなわち、ユーザが閾値を与えない場合)における制約要因となる。この問題をよりよく説明するために、以下の図を参照されたい。これは、より小さなスケール(すなわち、小さいrk)と大きなスケール(すなわち、大きなrk)を選択する効果を示すものである。それぞれの図の左部分は、妥当に選択された尺度係数を示しており、2つの平面に対してほぼ垂直に推定された表面法線と、テーブル全体に見える小さな辺を示している。しかし、スケールファクタが大きすぎる場合(右部分の図)、近傍の点集合の方が隣接する表面よりも範囲が大きくなり、推定された点特徴表現が、2つの平面のエッジにおいて回転表面法線により歪んでしまい、エッジも細部もぼやけてしまっている。

あまり詳細に踏み込むことなく、点の近傍の決定のスケールは、アプリケーションで要求されるレベルに基づき選択するものと考えれば十分である。 簡単に言えば、マグカップのハンドルと円筒形部分との間の縁部の曲率が重要である場合、スケールはそれらの細部を捕捉するのに十分小さくなければならないが、そうでなければ大きくしたほうが良い。

Estimating the normals

法線の推定

すでに http://lang.sist.chukyo-u.ac.jp/Classes/PCL/How3DFeaturesWork.html で例を見せているが、裏で何が起きているのかをよりよく説明するため、ここではその1つを修正する。

次のコードは、入力データセット内のすべての点の表面法線の集合を推定する。

In [ ]:
from py3d import *

print("Testing IO for point cloud ...")
pcd = read_point_cloud("./TestData/fragment.pcd")

print("Downsample the point cloud with a voxel of 0.05")
downpcd = voxel_down_sample(pcd, voxel_size = 0.05)

print("Compute the normal of the downsampled point cloud")
estimate_normals(downpcd, search_param = KDTreeSearchParamHybrid(
            radius = 0.1, max_nn = 30))

draw_geometries([downpcd])

ここでestimate_normals関数の呼び出しは内部的には次のことを行っている:

ポイントクラウドの各点pに対し

  1. pの最近傍点を取得し、
  2. pの表面法線$\boldsymbol{n}$を計算し
  3. $\boldsymbol{n}$が視点に一貫して向かうものであるかをチェックし、そうでなければ向きを変える

視点はデフォルトでは$(0,0,0)$であるが、PCLではsetViewPoint (float vpx, float vpy, float vpz);によって変更できる Open3Dでは orient_normals_towards_camera_location関数によって視点を変更したポイントクラウドを作成することができる。

PCLでは1点に対して法線を計算できるが、Open3Dではそのような関数やメソッドは見当たらない。Open3Dでやりたければ、Cropにより領域を絞るしかないだろう

Speeding Normal Estimation with OpenMP

OpenMPによる法線推定の高速化

Open3Dにはこれに相当するものは見当たらない