JP4686762B2 - 三次元形状の位置あわせ方法及びプログラム - Google Patents

Description

本発明は、実物からのモデリングを必要とする、工業生産分野、デジタルコンテンツ生成分野、文化財保護、シミュレーション、ロボットナビゲーション等、また、既存の形状モデルに対するマッチングを行うことによる三次元物体形状認識、三次元形状検索等において用いることのできる三次元形状の位置あわせ方法及びプログラムに関する。

三次元形状を計測するために、ステレオ、光切断法、焦点、ぼけ、空間符号化、光レーダーなどさまざまな原理に基づくセンサーが提案されているが、光学的原理に基づくセンサーは受光素子として画像センサーを利用しているため、一般に三次元形状情報が画素に対応した距離画像としてデータが獲得される。距離画像はCV、CG、CAD、VR、ロボット等の分野で、対象の認識や形状モデル生成のための重要な三次元形状情報の入力形式となっている。

多くの場合、距離画像は対象物体の表面の一部しか計測できない。物体の全表面について形状モデルを生成するためには、対象物をいろいろな方向から計測する必要がある。それぞれの距離画像はセンサー中心の座標系で表されており、一つの形状モデルにまとめるには、それぞれの距離画像の位置関係を求める必要がある。このような位置関係をデータから推定するのが位置あわせである。

距離画像の位置あわせは、粗い位置あわせ(coarse registration)と細かい位置あわせ(fine registration)の二つの段階に分けることができる。粗い位置あわせの段階では、全く異なった視点から計測された距離画像間をある程度の誤差範囲内に収めるような位置あわせを行う。一旦粗く位置あわせできれば、それを初期値として細かい位置あわせを適用することができる。ICP(Iterative Closest Point) アルゴリズム（非特許文献１参照）、直接法（非特許文献２，３参照）、同時位置あわせによるモデリング手法（非特許文献４参照）は細かい段階に属する位置あわせ手法である。本発明は、粗い段階の位置あわせ手法を提案し、それを初期値として使用して細かい位置あわせ手法を適用することによってモデル生成を行う。

粗い位置あわせは、対象物にマーカーを装着したり、GPS などの電磁気的ポジショニングデバイスを利用したり、回転台やロボットアームのような能動的デバイスを利用したりして得ることもできるが、対象の大きさ、材質、計測する環境によっては、このようなデバイスが必ずしも使えるとは限らない。また、必ずしもそのような情報がデータと一緒に提供されるとも限らない。GUI を用いて人手により粗い位置あわせを与えることも可能であるが、データ数が増えると作業の労力は大きく、計算機によって補助されることが望ましい。

粗い位置あわせ手法は、対象認識として行われてきた研究と深い関連がある。最も代表的な手法は、位置あわせする幾何変換(剛体の場合ユークリッド変換) に対して不変な特徴量を使用する手法である。曲率は微分幾何学によって定義された不変量であり、数学的には極限的に微小な領域から求められる。実際の距離画像から曲率を計算するには、有限な範囲のデータを用いて安定に求めるために、さまざまな工夫が提案されてきた（非特許文献５参照）が、そのほとんどは離散的なデータに局所的に何らかの連続的な関数を当てはめることに基づいている。Feldmar とAyache （非特許文献６参照）は、曲率を計算する際に求められた主方向と法線によって構成される局所直交基底の組にRANSAC（ランダムに選択した最低限のサンプルを元に生成した仮説を全データで検証することを繰り返すことによる頑強なモデル推定アルゴリズム 非特許文献３０参照） を適用することによって、粗い位置あわせを行っている。曲率は二つの成分によって表され(最大−最小、ガウシアン−平均、curvedness−shape index) （非特許文献７参照）、対応づけを行うには必ずしも充分な情報ではなく、ランダムサンプリングや投票などの手法により多くの候補から選択することが必要になる。曲率に付帯するより大局的な情報も利用することで候補の検索を減らし安定した位置あわせを行うことができる。Kehtarnavaz and Mohan （非特許文献８参照）は曲率の符号によって分割した領域をグラフで表現し、部分グラフマッチングにより位置あわせを行っている。Higuchi ら（非特許文献９参照）は単位球面上にマッピングした曲率に類似した量のマッチングにより回転成分の推定を行っている。Krsek ら（非特許文献１０参照） は平均曲率が0 となるような曲線のマッチングをランダムサンプリングによって行っている。

入力形状を部分形状のセグメントに分割し、類似したパラメータで表されるセグメント間同士でマッチングを行う手法は、比較的大きな特徴を利用していると言える。Faugeras とHebertは平面や二次曲面で曲面を分割し、分割されたセグメント間の対応関係の仮説を検証していくことによって位置あわせを行っている（非特許文献１１参照）。Kawai ら（非特許文献１２参照） は三角パッチを利用している。これらのセグメントに基づく手法は、セグメントへの分割(分節) が安定に行えることが仮定されているが、ノイズやオクルージョンのある自由曲面は必ずしも対応するセグメントが同様に分割されるとは限らない。Wyngaerd とVan Gool （非特許文献１３参照）は二箇所で曲面に接する直線の接点が描くbitangent curve を利用して位置あわせを行っている。

微小領域に基づく特徴は安定性と記述能力の問題があり、大域的な特徴は分節の安定性とオクルージョンへの頑強性の問題がある。Stein とMedioni は、中心頂点を囲む円周上の曲面の法線から求められる“splash” と呼ばれる特徴量を提案し（非特許文献１４参照）、Chua とJarvis は曲面と球面の交差を接平面上に投影した曲線から求められるpoint signature を提案している（非特許文献１５参照）。Johnson とHebert は中心頂点で接平面に垂直な画像平面を法線周りに回転して集めた計測点のヒストグラム画像をspin image とし（非特許文献１６参照）、対称認識に応用している。Spin image は固有画像展開により次元圧縮された空間で対応点探索を行っている。Huber とHebertはspin image を利用し、見えに基づく検証と組み合わせることにより、複数距離画像の自動位置あわせを行う手法を提案している（非特許文献１７参照）。これらの手法は曲率より広い領域からより豊富な情報の記述を得ることにより、対応探索を容易にしようとする試みである。

特徴による対応付けを行うのではなく、位置あわせの状態を記述する誤差関数を最適化することにより解くアプローチもある。粗い位置あわせには良い初期値が提供されているわけではなく、探索するパラメータの空間が広いので、確率的な最適化手法が用いられる。Blais とLevine はsimulated annealing(SA)を用い（非特許文献１８参照）、Brunnstrom とStoddart は最適な対応関係をgenetic algorithm (GA) で探索している（非特許文献１９参照）。Silva らは複数距離画像の位置あわせを、誤差関数が最小になるように位置あわせパラメータをGA で最適化することによって解いている（非特許文献２０参照）。このような確率的最適化に基づく手法は非常に多くの回数誤差関数を評価しなくてはならない。

P. J. Besl and N. D. McKay: "A method forregistration of 3-D shapes", IEEE Transactions on Pattern Analysis and MachineIntelligence, 14, 2, pp. 239−256(1992). B. K. P. Horn and J. G. Harris: "Rigid bodymotion from range image sequences", CVGIP: Image Understanding, 53, 1, pp. 1−13 (1991). M. Yamamoto, P. Boulanger, J.-A. Beraldin and M.Rioux: "Direct estimation of range flow on deformable shape from a video raterange camera", IEEE Transactions on Pattern Analysis and Machine Intelligence, 15, 1, pp. 82−89 (1993). T. Masuda: "Registration and integration ofmultiple range images by matching signed distance fields for object shapemodeling", Computer Vision and Image Understanding, 87, pp. 51−65 (2003). P. J. Besl and R. C. Jain: "Invariant surfacecharacteristics for 3D object recognition in range images", Computer Vision,Graphics, and Image Processing, 33,pp. 33−80 (1986). J. Feldmar and N. Ayache: "Rigid, affine andlocally affine rewgistration of free-form surfaces", International Journal ofComputer Vision, 18, 2, pp. 99−119(1996). J. J. Koenderink: "Solid Shapes", The MIT Press(1990). N. Kehtarnavaz and S. Mohan: "A framework forestimation of motion parameters from range images", Computer Vision, Graphics,and Image Processing, 45, 1, pp.88−105 (1989). K. Higuchi, M. Hebert and K. Ikeuchi: "Building3-D models from unregistered range images", Graphical Models and ImageProcessing, 57, 4, pp. 315−333(1995). P. Krsek, T. Pajdla and V. Hlavac:"Differential invariants as the base of triangulated surface registration",Computer Vision and Image Understanding, 87, pp. 27−38 (2002). O. D. Faugeras and M. Hebert: "Therepresentation, recognition, and locating of 3-D objects", The InternationalJournal of Robotics Research, 5,3 (1986). Y. Kawai, T. Ueshiba, T. Yoshimi and M. Oshima:"Reconstruction of 3D objects by integration of multiple range data",Proceedings of the 11th International Conference on Pattern Recognition, pp.154−157 (1992). J. V. Wyngaerd and L. V. Gool: "Automatic crudepatch registration: Toward automatic 3d model building", Computer Vision andImage Understanding, 87, pp. 8−26(2002). F. Stein and G. Medioni: "Structural indexing:efficient 3-d object recognition", IEEE Transactions on Pattern Analysis andMachine Intelligence, 14, 2, pp.125−145 (1992). C. Chua and R. Jarvis: "3-d free-form surfaceregistration and object recognition", International Journal of Computer Vision,17, 1, pp.77−99 (1996). A. E. Johnson and M. Hebert: "Using spin imagesfor efficient object recognition in cluttered 3d scenes", IEEE Transactions onPattern Analysis and Machine Intelligence, 21, 5, pp. 433−449 (1999). D. Huber and M. Hebert: "Fully automaticregistration of multiple3d data sets", Image and Vision Computing, Vol. 21, pp.637-650(2003). G. Blais and M. D. Levine: "Registeringmultiview range data tocreate 3D computer objects", IEEE Transactions onPattern Analysisand Machine Intelligence, 17, 8, pp. 820−824 (1995). K. Brunnstr¨om and A. Stoddart: "Geneticalgorithms for free-form surface matching", Proceedings of the 13thInternational Conferenceon Pattern Recognition, Vol. 4, pp. 689−693 (1996). L. Silva, O. R. Bellon and K. L. Boyer:"Enhanced, robust genetic algorithms for multiview range image registration",Proceedings of the4th International Conference on 3-D Digital Imaging andModeling(2003). G. Sandini and V. Tagliasco: "Ananthropomorphic retina-like structure for scene analysis", Computer Vision,Graphics, and Image Processing,14,pp. 365−372 (1980). K. Hotta, T. Mishima and T. Kurita: "Scaleinvariant face detection and classification method using shift invariantfeatures extracted from log-polar image", Transactions of IEICE, E84-D, pp. 867−878(2001). G.Wolberg and S. Zokai: "Robust imageregistration using log-polar transform", Proceedings of the IEEE InternationalConference on Image Processing (ICIP2000) (2000). T. Masuda: "Surface curvature estimation fromthe signed distance field", Proceedings of the 4th International Conference on3-D Digital Imaging and Modeling, pp. 361−368 (2003). H. Murase and S. Nayer: "Visual learning andrecognition of 3-Dobjects from appearance", International Journal of ComputerVision,14, pp. 5−24 (1995). S. A. Nene and S. K. Nayar: "A simple algorithmfor nearest neighbor search in high dimensions", IEEE Transactions on PatternAnalysis and Machine Intelligence, 19,9, pp. 989−1003 (1997). P. J. Rousseeuw and A. M. Leroy: "RobustRegression and Outlier Detection", Wiley (1987). D. W. Eggert, A. Lorusse and R. B. Fisher:"Estimating 3-D rigid body transformations: a comparison of four majoralgorithms", Machine Vision and Applications, 9, pp. 272−290 (1997). "The stanford 3d scanning repository", http://graphics.stanford.edu/data/3Dscanrep/. M. A. Fischler andR. C. Bolles, Random Sample Consensus: A Paradigm for Model Fitting withApplications to Image Analysis and Automated Cartography, Communications of theACM, Vol. 24, pp. 381−395, 1981

本発明は、精密な位置あわせを行うための初期値として、ばらばらの座標系で記述されたデータの大まかな位置関係を求めることを目的としている。これによって、類似画像検索に基づく距離画像の粗い位置あわせ手法を提供する。

図１は、本発明による三次元形状の位置あわせ方法の概略説明図である。
ステップＳ１：特徴としては接平面を画像平面としてlog-polar座標系で記述された局所距離画像を使用する。高速に生成するために符号付距離場により入力距離画像を粗く表現し、高速に局所距離画像を生成するように実装している。なお、本発明はこの実装方法に限定されず、距離画像から局所log-polar距離画像を生成する際に、例えば接平面に距離画像を直接正射影するような生成手法も用いることができる。
ステップＳ２：生成された画像は角度方向に周期的であり、対応する画像間では角度方向に不定性がある。この不定性を正規化して回転に対して不変な量にするために、法線周りの角度軸についてFourier 級数展開を用いてパワースペクトラムに変換している。本発明はこの正規化方法に限定されず、他の直交展開手法やモーメント等の手法も用いることができる。
ステップＳ３：法線周りの回転に対して不変とした画像である特徴量について、最も近い対を探索することにより対応点を求める。画像間の距離を求める際に、大きな固有値に対応する固有画像のみ使用する次元圧縮を行うことにより高速に対応関係を求められるように実装することができるが、本発明は、使用している次元圧縮の方法や最近傍点の探索手法には限定されない。
ステップＳ４：対応点に基づき距離画像間の位置関係を決定する。
ステップＳ５：全距離画像の位置関係を決定し、位置あわせ結果とする。

本発明は、局所Log-Polar 距離画像を特徴として用い、複数の距離画像の粗い位置あわせを行う。入力距離画像を符号付距離場によって粗くサンプリングすることにより、局所Log-Polar 距離画像を高速に生成できる。局所Log-Polar 距離画像をフーリエ変換によって圧縮し、さらに固有画像によって次元圧縮し、最近傍探索によって対応点対を生成する。検証によって誤対応を取り除き、各入力距離画像間のユークリッド変換をランダムサンプリングによって推定し、その変換に属するinlierが最大になるような木構造としてView Tree を構築することにより、入力距離画像の粗い位置あわせが行われる。粗い位置あわせ結果を、同時位置あわせ手法の初期状態として使用することにより、統合形状を生成することができる。

本発明による提案手法は、Point Signature（非特許文献１４，１５参照）やSpin Image （非特許文献１６，１７参照）などの手法と類似点がある。Point Signature に類する特徴量は、注目点を中心とした円周状の曲線上での特性から求められるが、提案手法では円盤状領域内の接平面からの奥行き値をそのまま使用する。提案手法は特徴量が画像である点でSpin Image にも類似しているが、Spin Image は接平面に垂直な画像平面で切った断面の集積で生成される。Spin Image が画像平面を法線の回りに回転して集積することにより回転不変性を持たせているのに対して、提案手法では回転方向にパワースペクトルをとることで回転不変性を持たせている。同一の半径であっても、SpinImage は原点から線状に伸びるパターンが生成されるのに対して、局所Log-Polar 距離画像では画素が埋まっており、フーリエ変換で圧縮したことによる誤対応も相関をとることによって検証することができる。

原点回りでの回転が平行移動で表すことができるのがLog-Polar 座標系を使用する利点であるが、もう一つの利点は、原点からの距離を対数圧縮している点である。原点から遠い部分は粗くサンプリングされ、画素数は中心部分より相対的に少なくなるので、計測されなかった部分を含むことによる影響を少なくすることができる。また、Log-Polar 画像の半径方向の画素を増やすことによる被覆半径の増加は指数的なので、距離画像全体を含むようなLog-Polar 距離画像も比較的少ない画素数で生成することができる。

局所Log-Polar 距離画像を位置あわせに適用する場合、その被覆が入力距離画像間の重なり幅より小さいことが必要となる。現在、符号付距離のサンプリング間隔δ は対象物の1/64 程度の大きさにしている。これは実用的な時間で局所Log-Polar 距離画像が生成され対応点対が探索できることから経験的に使用している設定であるが、R が4 以上は必要であることが実験的に得られているので、この場合、大まかには距離画像の1/8 程度の幅は重なっていることが必要になる。

一方で、形状が非常に単純な場合には、被覆が局所的過ぎると、局所的に類似した形状が多すぎて対応点対が正しい確率が減る。その場合には、ある程度被覆を大きくする必要があるが、重なりが少ない距離画像間で位置あわせすることができなくなる。非常に単純な形状は、微細な部分に特徴の情報が局在しているとも考えられ、そのような距離画像については不要な対応点対を生じさせないように特徴点を抽出したほうが効率的ではあるが、特徴点のある箇所が重ならないと位置あわせできないことになる。また、当然ながら対称な形状は形状だけで位置あわせするのは困難であり、テクスチャや事前知識などの補助的な情報が必要になる。

提案手法で最も計算時間を要しているのは、最近傍検索である。最近傍検索は、学習やデータマイニングの分野で研究が進んでおり、より高速なアルゴリズムを利用すれば、大幅な高速化を図ることができる。現在、固有画像を用いて圧縮しているが、経験的には入力距離画像によらず類似した固有画像が得られる傾向があり、最適性にこだわらなければ固有画像を共通化できる可能性もある。
提案手法は、誤った位置あわせを対話的に修正できるようなGUI や同時位置あわせアルゴリズムと組み合わせることにより、実用的な位置あわせソフトウェアを構成することができる。また、ロボットのための対象認識・環境地図生成・経路探索などにも応用できる。

本発明によれば、局所Log-Polar 距離画像を特徴として用い、複数の距離画像の粗い位置あわせを行う手法を提供することができ、従来用いられている特徴量に比べ情報が豊富であり、誤対応の検証もできる。

（特徴画像の生成）
図３は、特徴画像の生成を説明する図である。図示したように、入力距離画像はまず符号付距離場(SDF) で粗くサンプリングしなおす。表面上の各最近点を中心として局所Log-Polar 距離画像(LR) を生成する。図で黒く示されている未計測領域は奥行き値0で埋める。θ 軸方向にフーリエ変換をかけパワースペクトルを求め(FLR)、固有画像展開により圧縮する。FLR は少数の固有画像の線形結合として近似でき、その係数をCFLRとする。図ではCFLR として、再構成されたFLR を示している。元のFLR が近似されていることがわかる。以下、これらについて順次説明する。

（局所Log-Polar 距離画像）
対象物体について計測した入力距離画像をSα(1 <=α <=NS)とする。この入力距離画像から局所Log-Polar（「対数極座標」） 距離画像を生成する(図２、図３)。図３に示すように、入力距離画像はまず符号付距離場(SDF) で粗くサンプリングしなおす。表面上の各最近点を中心として局所Log-Polar 距離画像(LR) を生成する。図で黒く示されている未計測領域は奥行き値0 で埋める。

局所Log-Polar 距離画像は、曲面上のある点の接平面を画像平面として、中心点の周りに張ったLog-Polar 座標系（非特許文献21〜非特許文献23参照）で表した距離画像である。画像平面上の直交座標系での位置が(u、 v) であるとき、Log-Polar 座標系による位置は(ζ、 θ) = (log ｒ、 arctan(v/u)) となる。ここでr = √（u² + v²）である。局所Log-Polar 距離画像(LR) はLog-Polar 座標(ζ、 θ)に正射影される曲面上の点の法線方向の高さを値として持つ。

入力距離画像をδ の精度で位置あわせすることを目指す場合、入力距離画像の点群の密度は非常に高い場合があり、δ より細かい形状情報は不要である。粗い形状表現を得るために、Masuda が使用している曲面の符号付距離場(SDF: Signed Distance Field) を利用する（非特許文献４参照）。サンプル点p での符号付距離場SDF のサンプルは、p から曲面へ最も近い点：c、そこでの法線：n、p への符号付距離：s の組で与えられ、これらは式p = c + sn を満たす。（図９参照）。法線方向n はc からp の方向の単位ベクトルとして、局所的な平面を当てはめて法線を求める手法よりも安定して求められる。符号付距離場SDF から各サンプルについて距離画像LR を求めるには、各距離画像を符号付距離場SDF によって標本化し、曲面上の最近点c を中心点、n を接平面の法線ベクトルとして利用する。曲面上でδ 四方の領域辺りに一点程度の密度で距離画像LR を生成すれば充分なので、|s| < δ/2 を満たすサンプルを使用する。

各最近点c_i を中心とした距離画像LR を生成するには、同一距離画像の符号付距離場SDF の異なる最近点c_j のうち、半径Rδ 内にあり(||c_i 。 c_j || < Rδ)、法線が同じ向き(n_i ・ n_j > 0) のものを使用する。各c_j を接平面である画像平面に正射影し、対応するLog-Polar 座標と、高さn_i ・ (c_j − c_i) を求める(図２)。同一Log-Polar 座標の画素に射影される高さのうち、最大値を画素値とする。距離画像には未計測領域があるが、そのような部分は高さ0とする。第α番距離画像S^α から生成された符号付距離場SDF の最近点c^α _iにおける距離画像LR をLR[c^α _i](ζ、 θ) と表す。

実際に画像を生成するには、座標値(ζ、 θ) を量子化しなくてはならない。θ 軸のとりうる範囲[−π、 π) を2Nθ 画素に分割する場合、(ζ、 θ) 空間の両方の軸を均等に分割するには、i番目のζ 座標値はζ_i = logr_i = (π/Nθ)i を満たす必要があり、ζ 軸に必要な画素数はN_ζ = (Nθ/π) log(R) となる。ここで、R は、δ を単位とした、局所画像が被覆する領域の半径である。中心に近い領域ζ < 0 は使用しない。例えばN_θ = 16でR = 8 の場合、N_ζ = 11 なので生成される距離画像LR の大きさは[0、 10] × [−16、 15] = 352 画素となる。Log-Polar 座標のθ を決めるための局所直交座標としては、主方向によって構成される基底を利用しているが（非特許文献２４参照）、そのとり方に結果が依存しないように以下の処理を構成する。

（パワースペクトラム）
図３に示すように、生成された局所距離画像に対してθ 軸方向にフーリエ変換をかけパワースペクトルを求め(FLR)、固有画像展開により圧縮する。FLR は少数の固有画像の線形結合として近似でき、その係数をCFLR とする。図３ではCFLR として、再構成されたFLR を示している。元のFLR が近似されていることがわかる。

距離画像LR はθ 軸方向に周期的なので(LR(ζ、 θ) = LR(ζ、 θ + 2π))、フーリエ級数で展開することができる。フーリエ展開して求めたθ 軸方向のパワースペクトラムは、θ の位相に不変な量となる：

ここでk は周波数とする。ここで定義したFLR を、各中心点に付随する特徴量とし、対応付けに使用する。フーリエ変換はFFT を用いて高速に計算できる。距離画像LR は実関数なのでFLR は偶関数となる(FLR(ζ、 k) = FLR(ζ、−k)) ので、FLR の大きさは距離画像LR の半分となる。例えば前節で示した例の場合、FLR の大きさは[0、 10] × [0、 15] = 176 画素となる。FLR の画素数をDFLRと表す。

（高次局所自己相関特徴）
不定性の正規化を、上述のようにフーリエ級数展開してパワースペクトラムとすることに代えて、高次局所自己相関特徴とすることによって行うこともできる。N次の変位a1,...,aNに対応する画像f(r)の高次相関関数は、以下の式によって定義される。

変位や次数の組み合わせは無限に存在し得るが，画像データに適用する場合は、実装上は、N=0,1,2でa1,...,aNを3x3の範囲内にとり、平行移動による冗長性を取り除くと、図１０に示すように、35個のパターンに集約され、一つの画像について35次元のベクトルを特徴量として取り出すことができる。図１０は、高次局所自己相関のマスクパターンを示しており，3x3マトリックスに示す数字の回数の画素値の積をとり，各パターンについて画像全体で和をとることによって特徴ベクトルを計算する。

（固有画像による次元圧縮）
図３に示すように、フーリエ変換によりパワースペクトル(FLR)を求めた後、固有画像展開により圧縮する。FLR の情報は、ζ が大きく周波数が低い部分に偏在している。FLR の集合の固有画像を求めることにより、低い次元の空間に情報を圧縮することができる。FLR は少数の固有画像の線形結合として近似でき、その係数をCFLR とする。図３ではCFLR として、再構成されたFLR を示している。元のFLR が近似されていることがわかる。
全ての入力距離画像について、NFLR 個のFLR が生成されているとする。これらのFLR を一行として積み重ねることにより、N_FLR × D_FLR 行列M_FLR を生成する。FLR の共分散行列(M_FLR)^TM_FLR の固有ベクトルは、M_FLR を特異値分解(SVD)することで求めることができる。固有ベクトルをEFLR[c^α _i、 l]とし、FLR の展開係数をCFLR[c^α _i](l) = EFLR[c^α _i、 l] ・ FLR[c^α _i]とするとき、大きいほうからD 個の固有値に対応する固有ベクトルを用いて、元のFLR は次のように近似できる。

FRL の代わりにこの展開係数CFLR[c^α _i](l)、 (1 <=l <=D) を用いて、D 次元の空間で対応点を検索する。
類似した画像を検索するために固有画像を利用する場合、平均値を引いたり輝度のノルムを1 にそろえたりするような正規化を適用するが（非特許文献２５参照）、提案手法ではFLR の絶対値が無視できないため、そのような正規化は行わない。

（対応付け）
図４は、対応付けを説明する図である。図示したように、異なる入力距離画像のCFLR で互いに最も近いもの同士を対応点対とする(左)。FLR が似ていても不要な不変性のためにつく誤対応を距離画像LR の相互相関によって検証する(中)。対応点間にRANSAC を適用して、入力距離画像間のユークリッド変換と、inlier（該当する日本語訳は無いが、反対語のoutlierは「外れ値」と訳されている） となる対応点対を求める(右)。以下、これらについて順次説明する。

（対応点の探索）
全ての入力距離画像の各中心点についてCFLR を生成した後、異なる入力距離画像のCFLR 間で互いに最も近い組を探す。入力距離画像S^α の中心点c^α _iについて求められたCFLR[c^α _i] について、他の入力距離画像Sβ(β≠α) の中心点c^β _j について求められたCFLR[c^βj ] が最も近い場合、点の対応[c^α _i、 c^β _j ] を対応点のリストに加える。

CFLR はD 次元のベクトルなので、汎用の最近傍点の探索アルゴリズムを適用することができる。例えば、Nene とNayerにより提案された手法（非特許文献２６参照）を使用することができる。各軸への射影を利用して、最近傍を探す点を中心とする幅2εの超立方体内の候補点を抽出し、その候補の中で最も近い点を全探索するという単純なアルゴリズムである。

このアルゴリズムを適用するには、閾値_ の設定が重要である。閾値ε が大きすぎると全探索する候補が増え探索効率が悪くなり、小さすぎると最近点の探索に失敗する。今回の実装では

により設定している。ここで、V は点が分布する体積であり、各軸で一次元のLMS(LMedS:
LeastMedian of Squares) 推定する際に使用する半数の点を覆う最小幅から推定している（非特許文献２７参照）。
各CFLR について最近傍点に対応づけた後、相互に最近傍になっていないような対応点は除去する(図４ 左)。

（相互相関による検証）
FLR の性質として、次の式で表されるθ 軸方向への位相に対する不変性は望ましい性質である。

しかしながら、FLR には次の式で表されるような符号や向きに対する不変性も兼ね備えており、これらは望ましいとはいえない。

θ 軸方向のみのパワースペクトラムを用いているので、ζ の値によって位相が異なるようなずれ方に対しても不変性を持ってしまう。

このような不要な不変性による誤対応は、距離画像LR の相互相関をとることによって検証できる。向き± と位相Δθ を変えたLR[c^β _j ](ζ、±θ + Δθ) は、Log-Polar 座標系の性質から、ベクトル要素のシフトや置換などの簡易な操作で生成することができる。対応づけられているLR[c^α _i] とLR[c^β _j ] の対について、向きと位相を変え、次の正規化相関を計算する。

もし全ての± とΔθ の組み合わせについての最大値が閾値(実装では0．5) より小さい場合、または最大値になるのが逆向き(向きが’−’) の場合、この点対の対応関係は誤りであると判別され除去される(図４ 中)。

現在の実装では、FLR が生成された段階で距離画像LR はメモリから消去し、この検証を行うのに必要な距離画像LR のみ再生成しているが、実験結果では充分高速な実行速度が得られる(実施例参照)。

（剛体性の検証）
複数視点から計測された距離画像の位置あわせする際に必要となる幾何変換は剛体の運動を表すユークリッド変換である。各入力距離画像間には一つのユークリッド変換が対応しているはずである。
これまでの手続きで生成され検証によって除去されずに残っている点対を、対応する入力距離画像の対毎に分類する。例えば対応点対[c^α _i、 c^β _j ] は入力距離画像対[S^α、 S^β] の要素に分類される。それぞれの入力距離画像対に属する対応点対について、RANSAC を適用し、ユークリッド変換を推定するとともにinlier である対応点対を求める。

入力距離画像対[S^α、 S^β] に属する対応点対から三対[c^α _i1 、 c^β _j1 ]、[c^α _i2 、 c^β _j2 ]、 [c^α _i3 、 c^β _j3 ] を選ぶ。三対の対応点対があれば、ユークリッド変換T = {R、 t} を一意に求めることができる（非特許文献２８参照）。この変換を、同じ入力距離画像対に属する対応点に適用し、条件||Tc^α _i。 c^β _j || < δ、 Rn^α _i・ n^β _j >
τn(実装ではcos(π/8)) を満たす点対をinlier とみなし、その数を数える。三対の対応点対の選択方法として、対応点数が少ない場合は全ての組み合わせ、多い場合はランダムに選択し、以上の手続きを繰り返す。実装では対応点対が22 を超える入力距離画像対に対しては、10000 回ランダムに選択している。入力距離画像対[S^α、 S^β] 間で最も多くのinlier が属する変換をその入力距離画像対間の変換T^β _αとする(図４ 右)。

（View Tree ：木構造の生成）
図５は、それぞれの入力距離画像対についてRANSAC アルゴリズムにより選び出されたinlier の数の表(上) と、全体のinlier の数を最大化するような入力距離画像のspanning tree として生成されたView Tree(下)を示している。

各入力距離画像間でRANSAC によってinlier を選択した後、全体としてinlier 数を最大化するようなspanning tree を構成する。この入力距離画像の木構造をView Tree と呼ぶ。例えば、図５ において、入力距離画像対[S₁、 S₂]、[S₂、 S₃]、[S₁、 S₃] にはそれぞれ143、22、5 対のinlier 対応点対があるが、この3 つの入力距離画像を接続するには[S₁、 S₂] と[S₂、 S₃] で充分であり、inlier の数が少ない[S1、 S3] は使用しない。View Tree は次のようなアルゴリズムで構成することができる。

1．空のグラフを用意し、全ての入力距離画像対を対象とする。
2．対象となる入力距離画像対の最大inlier 数がτconnection 以下になるまで以下の操作を繰り返す。
2．1．最大のinlier 数の入力距離画像対をグラフに加える。
2．2．両端の距離画像がグラフに接続済みの入力距離画像対を対象から除く。

このアルゴリズムによって、最小のinlier 数がτconnection(実装では5) 以上で、構成するinlier 数の総和が最大化されたspanning tree としてのView Tree が構成できる(図５)。
基底距離画像Sbase(図５ ではS1) の変換をTbase = {I_3×3、 0₃}として、View Tree の辺に対応するそれぞれの入力距離画像対のユークリッド変換T^β _αを逐次的に適用していくことで、各入力距離画像S^α の基底距離画像S_base に対する相対的なユークリッド変換T_α を求めることができる。

（補正）
各入力距離画像S^α のユークリッド変換Tα が決まれば、全体の大まかな形状を組み立てることができる。上述のView Tree で使用されなかった対応点対も利用して、形状全体での補正を行う。上述の検証により残された点対全てのうち、条件||T_αc^α _i− T_βc^β _j || < δ、 R_αn^α _i・ R_βn^β _j > τ_n を満たす点対を選択し、各入力距離画像S^α についてc^α _iを対応点対の中点(T_αc^α _i+ T_βc^β _j )/2 に対応付けることにより、変換T_α を更新する。この対応点対の選択と変換の更新を、条件を満たす対応点対間の位置あわせ誤差の変動が充分小さくなるまで繰り返す。

対応点間での補正では、対応がついている部分しか利用していない上に、View Tree で離れた位置にある入力距離画像間では位置合わせ誤差が累積したままである。このような位置あわせ誤差は、全入力距離画像の同時位置あわせアルゴリズムを適用することで解消することができる。ここでは、Masuda により提唱されたアルゴリズム（非特許文献４参照）を利用した。提案手法で推定された、位置あわせ精度δ の範囲内での粗い位置あわせ結果を、まずサンプリング間隔をδ として同時位置あわせを行い、さらにサンプリング間隔を細かくすることにより精細な形状を得ることができる(図６)。

図６は、提案手法によりδ = 4mm の粗さで位置あわせされた入力距離画像(左)。これを初期値として、同時位置あわせ手法によって統合した形状モデル(中)。さらにその結果を初期値として、サンプリング間隔をδ = 0.5mm に細かくして生成した精細な形状モデル(右)を示している。

ここで使用した同時位置あわせアルゴリズムは符号付距離場(SDF) を形状の記述として使用している点で提案手法と共通点がある。位置あわせと同時に形状統合も行い、統合形状は符号付距離場SDFで記述されているので、新たな距離画像や他の部分を統合した符号付距離場SDF を提案手法の入力として利用することにより、大きなデータセットを分割して処理することができる。

上述した手法の説明に使用してきた距離画像は、Stanford 3D Scanning Repository （非特許文献２９参照） から入手した‘Bunny’というデータセットである。全部で10 枚の距離画像(図５) で構成されている。対象物は約25cm の大きさがあり、符号付距離場SDF はδ = 4mm のサンプリング間隔で生成した。局所Log-Polar 距離画像は、N_θ = 16、R = 8、D_FLR = 11 × 16 = 176 の設定により、全部でN_FLR = 19156 枚の画像が生成された。これらの画像はD = 8 枚の固有画像により圧縮され(図３)、そのときの累積寄与率は92.9% であった。上述した方法によって決定したεにより、17916 点についての最近点が探索され、そのうち3711 対が相互に最近傍になっている対応点対であった。相互相関による検証により1129 対が残り、剛体性の検証により462 対がinlier として残され、View Tree は366 対の対応点対によって構成された(図５)。計算時間として、入力距離画像の読込みと符号付距離場SDF のサンプリングには2 分、局所Log-Polar 距離画像の生成と圧縮には4 分、最近点の探索に10分、対応点対の検証に1 分を要した。計算時間はXeon 1.7GHzシングルプロセッサーで計測した。

図８は、同一のデータに対して、異なる設定を適用した結果を示す表である。この表は、さまざまな設定によるBunny データの位置あわせ結果を示している。ここで、D_FLR = N_θ ×N_ζ、 CP%: D 次元まで使用したときの累積寄与率、#inliers: View Tree を構成するinlier の総数である。

R とD が小さすぎると位置あわせが失敗しており、情報量が少なすぎて正しい対応点対が見つからなかったと考えられる。N_θ は小さい値に設定しても位置あわせはできているが、最後に残されるinlier の数が少なく、性能は劣化している。R = 4、 D_FLR = 8、N_θ = 4 がおおよそ最低限必要とされる設定である。D を大きくすれば性能は向上するが、最近点を探索する空間の次元が高くなり多くの計算時間が必要となる。N_θ は大きくするとDFLR が大きくなりSVD の計算量が増えるが、あまり大きくしても性能が向上するわけではない。一方、R を大きくしすぎると、重複部分が少ない入力距離画像間での位置あわせに失敗するようになるため、ある程度の情報量が確保できる範囲で小さいほうが位置あわせには適している。

提案手法を、同じくStanford 3D Scanning Repository （非特許文献２９参照）から取得した‘Dragon’ データセットに適用した(図７)。図７は、粗く位置あわせした結果(左)、これを同時位置あわせすることにより得られた統合形状モデル(中)、及び生成された精細な統合形状モデル(右)を示す図である。

‘Dragon’データセットは全部で70 枚の距離画像を含むが、まず、そのうち30枚について、提案手法で粗い位置あわせを行った。対象物の大きさは約25cm であり、δ = 4mm、N_θ = 8、R = 4、D_FLR = 8と設定した。38745 枚の局所Log-Polar 距離画像が生成され、609 対の対応点対でView Tree が構成され、738 対を用いて対応点対による補正が行われた(図７ 左)。計算時間は約1 時間であった。

粗い位置あわせ結果を初期値として、同じδ = 4mm で同時位置あわせアルゴリズムを適用することにより、統合形状が生成された(図７ 中)。この統合形状に、残りの40 枚を4 回に分け、粗く位置あわせして統合形状を生成する手続きを繰り返して付加していき、全距離画像について形状統合を行った。その結果を元に、サンプリング間隔をδ = 1mm まで細分化することによって、全距離画像を用いた詳細な形状モデルを生成することができた(図７ 右)。

大きなデータセットを処理する際に、一度に全部処理すると対応点対の探索の負荷が大きく、データ数が多くなると関係が把握しにくくなるため、分割して処理することができることは有用である。分割の方法はいろいろ考えられ、ここで行った方法の他に、距離画像を部分毎に分類し、部分毎に統合形状を作り、統合形状同士で位置あわせし、最後に再び距離画像を使用して同時位置あわせをしなおすような方法も考えられる。

本発明による三次元形状の位置あわせ方法の概略説明図である。 局所log-polar距離画像の生成を説明する図である。 特徴画像の生成を説明する図である。 対応付けを説明する図である。 それぞれの入力距離画像対についてRANSAC アルゴリズムにより選び出されたinlier の数の表(上) と、全体のinlier の数を最大化するような入力距離画像のspanning tree として生成されたView Tree(下)を示す図である。 提案手法によりδ = 4mm の粗さで位置あわせされた入力距離画像(左)。これを初期値として、同時位置あわせ手法によって統合した形状モデル(中)。さらにその結果を初期値として、サンプリング間隔をδ = 0.5mm に細かくして生成した精細な形状モデル(右)を示す図である。 粗く位置あわせした結果(左)、これを同時位置あわせすることにより得られた統合形状モデル(中)、及び生成された精細な統合形状モデル(右)を示す図である。 同一のデータに対して、異なる設定を適用した結果を示す表である。 曲面の符号付距離場を説明する図である。 高次局所自己相関のマスクパターンを例示する図である。

Claims (4)

1. 計算機は、画像センサーにより対象物体の表面の三次元形状を計測した入力距離画像から、前記対象物体の表面上のある点の接平面を画像平面としてその点を中心に周りに張ったlog-polar座標系で記述された局所距離画像を生成し、
   前記局所距離画像に対して、前記対象物体の表面上のある点における法線まわりの角度軸についてフーリエ級数展開してパワースペクトラムを算出して特徴量とし、
   異なる入力距離画像間において最も近い特徴量を有する対応点を検索して点の対応関係を求め、
   求めた対応点に基づき前記入力距離画像間の位置関係を決定する
   ことから成る三次元形状の位置あわせ方法。
2. 前記入力距離画像を符号付距離場でサンプリングすることにより前記局所距離画像を生成する請求項１に記載の三次元形状の位置あわせ方法。
3. 前記特徴量の対応関係を求める際に、固有ベクトルを用いて固有空間で展開することにより次元圧縮し、次元圧縮した空間で最も近い値を有する対応点を検索する請求項１に記載の三次元形状の位置あわせ方法。
4. 画像センサーにより対象物体の表面の三次元形状を計測した入力距離画像から、前記対象物体の表面上のある点の接平面を画像平面としてその点を中心に周りに張ったlog-polar座標系で記述された局所距離画像を生成し、
   前記局所距離画像に対して、前記対象物体の表面上のある点における法線まわりの角度軸についてフーリエ級数展開してパワースペクトラムを算出して特徴量とし、
   異なる入力距離画像間において最も近い特徴量を有する対応点を検索して点の対応関係を求め、
   求めた対応点に基づき距離画像間の位置関係を決定し、
   決定された全距離画像間での位置関係を初期値として用いて、複数の入力距離画像を位置あわせすることにより、形状モデルを生成する、
   各手順を計算機に実行させるための三次元形状の位置あわせプログラム。
   

Images