ポイントセット登録

コンピュータビジョン、パターン認識、およびロボット工学において、ポイントセット レジストレーション(ポイントクラウド レジストレーションまたはスキャン マッチングとも呼ばれる) は、2 つのポイントクラウドを揃える空間変換(例: スケーリング、回転、および変換)を見つけるプロセスです。このような変換を見つける目的には、複数のデータ セットをグローバルに一貫したモデル (または座標フレーム) にマージすることや、新しい測定値を既知のデータ セットにマッピングして特徴を識別したり姿勢を推定したりすることが含まれます。生の 3D ポイントクラウド データは通常、LiDARおよびRGB-D カメラから取得されます。3D ポイントクラウドは、三角測量、バンドル調整、さらに最近ではディープラーニングを使用した単眼画像の深度推定などのコンピュータビジョン アルゴリズムから生成することもできます。画像処理および特徴ベースの画像レジストレーションで使用される 2D ポイントセット レジストレーションの場合、ポイントセットは、コーナー検出など、画像からの特徴抽出によって取得された 2D ピクセル座標である場合があります。ポイントクラウドレジストレーションは、自動運転、^[1]動作推​​定と3D再構成、^[2]物体検出と姿勢推定、^{[3] [4]}ロボット操作、^[5]同時位置推定とマッピング（SLAM）、^{[6] [7]}パノラマスティッチング、^[8]仮想現実と拡張現実、^[9]医療用画像処理など、幅広い分野で応用されています。^[10]

特殊なケースとして、3D 回転のみが異なる (つまり、スケーリングと変換がない) 2 つの点集合の登録は、ワバ問題と呼ばれ、直交プロクラステス問題にも関連しています。

この問題は次のようにまとめられる。^[11]有限次元実ベクトル空間における有限サイズの 2 つの点集合を、それぞれ との点を含むものと する(例えば、と が3D 点集合である場合の典型的なケースを は回復する)。問題は、動いている「モデル」点集合に適用する変換を、 と静的な「シーン」集合との差 (通常は点ごとのユークリッド距離の意味で定義される)が最小になるように見つけることである。言い換えれば、 から へのマッピングが求められ、これにより、変換された「モデル」集合と「シーン」集合の間の位置合わせが最適になる。マッピングは、固定変換または非固定変換で構成できる。変換モデルは と記述でき、これを使用して、変換された登録済みモデル点集合は次のようになる。 ${\displaystyle \lbrace {\mathcal {M}},{\mathcal {S}}\rbrace }$${\displaystyle \mathbb {R} ^{d}}$${\displaystyle M}$${\displaystyle N}$ ${\displaystyle d=3}$${\displaystyle {\mathcal {M}}}$${\displaystyle {\mathcal {S}}}$${\displaystyle {\mathcal {M}}}$${\displaystyle {\mathcal {M}}}$${\displaystyle {\mathcal {S}}}$${\displaystyle \mathbb {R} ^{d}}$${\displaystyle \mathbb {R} ^{d}}$${\displaystyle T}$

したがって、点集合登録アルゴリズムの出力は、距離関数 の定義された概念に従って、がに最もよく一致するような最適な変換 です。 ${\displaystyle T^{\star}}$${\displaystyle {\mathcal {M}}}$${\displaystyle {\mathcal {S}}}$${\displaystyle \operatorname {dist} (\cdot ,\cdot )}$

ここで、は最適化が検索しようとするすべての可能な変換の集合を表すために使用されます。距離関数の最も一般的な選択は、すべての点の ペアのユークリッド距離の2乗を取ることです。${\displaystyle {\mathcal {T}}}$

ここで、 はベクトル 2 ノルムを表し、は変換後にセット内の特定の点までの最短距離を達成するセット内の対応する点です。このような関数をリジッドレジストレーションで最小化することは、最小二乗問題を解くことと同じです。 ${\displaystyle \|\cdot \|_{2}}$${\displaystyle s_{m}}$${\displaystyle {\mathcal {S}}}$${\displaystyle m}$${\displaystyle {\mathcal {M}}}$

対応関係（すなわち、 ）が最適化の前に、例えば特徴マッチング技術を使用して与えられている場合、最適化では変換を推定するだけで済みます。このタイプの登録は、対応関係に基づく登録と呼ばれます。一方、対応関係が不明な場合は、最適化では対応関係と変換を一緒に見つけ出す必要があります。このタイプの登録は、同時姿勢および対応関係登録と呼ばれます。 ${\displaystyle s_{m}\leftrightarrow m}$

2つの点集合が与えられた場合、剛体レジストレーションは、一方の点集合を他方の点集合にマッピングする剛体変換を生成します。剛体変換は、任意の2点間の距離を変更しない変換として定義されます。通常、このような変換は平行移動と回転で構成されます。^[12]まれに、点集合がミラーリングされることもあります。ロボット工学とコンピュータービジョンでは、剛体レジストレーションが最も多くの用途があります。

2 つの点集合が与えられた場合、非剛体レジストレーションは、一方の点集合を他方の点集合にマッピングする非剛体変換を生成する。非剛体変換には、スケーリングやシアーマッピングなどのアフィン変換が含まれる。しかし、点集合レジストレーションの文脈では、非剛体レジストレーションは通常、非線形変換を伴う。点集合の固有モードの変化が分かっている場合、非線形変換は固有値によってパラメータ化される可能性がある。^{[13]非線形変換は、}薄板スプラインとしてパラメータ化されることもある。^{[14] [13]}

点集合登録のアプローチには、より一般的なグラフマッチング問題を解決するアルゴリズムを使用するものもあります。^[11]しかし、このような方法は計算量が多くなる傾向があり、剛体登録に限定されます。この記事では、変換に3D回転と平行移動（均一なスケーリングも含まれる可能性があります）が含まれていると想定される剛体登録のアルゴリズムのみを検討します。

PCL （ポイントクラウドライブラリ）は、 n次元ポイントクラウドと3Dジオメトリ処理のためのオープンソースフレームワークです。いくつかのポイント登録アルゴリズムが含まれています。^[15]

対応に基づく方法では、すべての点に対して推定上の対応が与えられていると仮定します。したがって、点集合と点集合の両方に点があり、対応が与えられている設定に到達します。 ${\displaystyle m\leftrightarrow s_{m}}$${\displaystyle m\in {\mathcal {M}}}$${\displaystyle {\mathcal {M}}}$${\displaystyle {\mathcal {S}}}$${\displaystyle N}$${\displaystyle m_{i}\leftrightarrow s_{i},i=1,\dots ,N}$

最も単純なケースでは、すべての対応が正しいと想定できます。つまり、ポイントは次のように生成されます。${\displaystyle m_{i},s_{i}\in \mathbb {R} ^{3}}$

ここで、は均一なスケーリング係数（多くの場合と仮定）、は適切な 3D 回転行列（次数の特殊直交群）、は 3D 変換ベクトルであり、 は未知の加法ノイズ（たとえば、ガウスノイズ）をモデル化します。具体的には、ノイズが標準偏差 のゼロ平均等方性ガウス分布に従うと仮定した場合、つまりの場合、次の最適化によって未知のスケール、回転、変換の最大尤度推定値が得られることがわかります。${\displaystyle l>0}$${\displaystyle l=1}$${\displaystyle R\in {\text{SO}}(3)}$${\displaystyle {\text{SO}}(d)}$${\displaystyle d}$${\displaystyle t\in \mathbb {R} ^{3}}$${\displaystyle \epsilon _{i}\in \mathbb {R} ^{3}}$ ${\displaystyle \epsilon _{i}}$${\displaystyle \sigma _{i}}$ ${\displaystyle \epsilon _{i}\sim {\mathcal {N}}(0,\sigma _{i}^{2}I_{3})}$

スケーリング係数が 1 で変換ベクトルが 0 の場合、最適化によってWahba 問題の定式化が回復されることに注意してください。集合の非凸性による最適化 ( cb.2 ) の非凸性にもかかわらず、 Berthold KP Hornによる独創的な研究では、スケール、回転、変換の推定を切り離すことで、( cb.2 ) は実際には閉じた形式のソリューションを許容することが示されました。 ^{[16]同様の結果が Arun}らによって発見されました。^[17]さらに、一意の変換を見つけるためには、各点集合に少なくとも非共線点が必要です。 ${\displaystyle {\text{SO}}(3)}$${\displaystyle (l,R,t)}$${\displaystyle N=3}$

さらに最近では、ブリアルズとゴンザレス・ヒメネスは、モデルセットに点、線、平面などのさまざまな3Dプリミティブが含まれる場合（モデルが3Dメッシュの場合）に、ラグランジュ双対性を使用した半正定値緩和法を開発しました。 ^[18]興味深いことに、半正定値緩和法は経験的に厳密であり、つまり、半正定値緩和法の解から 証明可能なグローバル最適解を抽出できます。${\displaystyle {\mathcal {M}}}$${\displaystyle {\mathcal {M}}}$

最小二乗法（cb.2 ）は外れ値が存在すると性能が著しく低下することが知られている。外れ値対応とは、生成モデル（cb.1 ）から外れた測定値のペアである。この場合、次のように異なる生成モデルを検討することができる：^[19]${\displaystyle s_{i}\leftrightarrow m_{i}}$

ここで、番目のペアがインライアであれば、外れ値のないモデル ( cb.1 )に従います。つまり、 は空間変換に小さなノイズを加えてから得られます。ただし、番目のペアが外れ値であれば、 は任意のベクトル になります。どの対応が外れ値であるかは事前にわからないため、生成モデル ( cb.3 ) による堅牢な登録は、現実の世界で展開されるコンピュータービジョンやロボット工学にとって極めて重要です。現在の特徴マッチング技術では、対応の 100 % 以上が外れ値になる可能性がある、非常に破損した対応を出力する傾向があるためです。^[20]${\displaystyle i-}$${\displaystyle s_{i}\leftrightarrow m_{i}}$ ${\displaystyle s_{i}}$${\displaystyle m_{i}}$${\displaystyle i-}$${\displaystyle s_{i}\leftrightarrow m_{i}}$${\displaystyle s_{i}}$${\displaystyle o_{i}}$${\displaystyle 95\%}$

次に、堅牢な登録のためのいくつかの一般的なパラダイムについて説明します。

最大合意は、空間変換の選択に対して生成モデル（ cb.1 ）と一致する対応の最大セットを見つけようとします。正式に言えば、最大合意は次の最適化を解決します。${\displaystyle (l,R,t)}$

ここで、 は集合 の濃度を表します。( cb.4 )の制約は、インライア集合内のすべての測定値のペアの残差が事前に定義されたしきい値 より小さくなければならないことを強制します。残念ながら、最近の分析では、問題 (cb.4) を大域的に解くのはNP 困難であり、大域的アルゴリズムでは通常、最悪の場合指数時間の計算量を要する分岐限定法 (BnB) に頼らざるを得ないことが示されています。 ^{[21] [22] [23] [24] [25]}${\displaystyle \vert {\mathcal {I}}\vert }$${\displaystyle {\mathcal {I}}}$${\displaystyle {\mathcal {I}}}$${\displaystyle \xi }$

コンセンサス最大化を正確に解くことは困難ですが、実際には非常に優れたパフォーマンスを発揮する効率的なヒューリスティックが存在します。最も人気のあるヒューリスティックの 1 つは、ランダム サンプル コンセンサス (RANSAC)スキームです。^[26] RANSAC は、反復的な仮説と検証の方法です。各反復で、この方法では、まず対応関係の総数から 3 つをランダムにサンプリングし、 Horn の方法を使用して仮説を計算します。 ^{[16]次に、(} cb.4 )の制約を評価して、実際にそのような仮説と一致する対応関係の数を数えます (つまり、残差を計算し、各測定ペアのしきい値と比較します)。アルゴリズムは、十分な対応関係を持つコンセンサス セットが見つかった後、または許可された反復の総数に達した後に終了します。各反復の主な計算は Horn の方法の閉形式のソリューションを実行することであるため、RANSAC は非常に効率的です。しかし、RANSACは非決定論的であり、外れ値比率が低い領域（例えば、下記）でのみうまく機能します。これは、外れ値比率に応じて実行時間が指数関数的に増加するためです。^[20]${\displaystyle N}$${\displaystyle (l,R,t)}$${\displaystyle \Vert s_{i}-lRm_{i}-t\Vert _{2}^{2}/\sigma _{i}^{2}}$${\displaystyle \xi }$${\displaystyle 50\%}$

高速だが不正確なRANSAC方式と正確だが網羅的なBnB最適化との間のギャップを埋めるために、最近の研究ではコンセンサス最大化を解決するための決定論的近似法が開発されている。^{[21] [22] [27] [23]}

外れ値除去法は、空間変換を推定する前に、非常に破損した対応関係のセットを前処理しようとします。外れ値除去の目的は、外れ値対応関係を維持しながら外れ値対応関係の数を大幅に減らすことです。これにより、変換の最適化がより容易かつ効率的になります (たとえば、外れ値比率が を超えると RANSAC はうまく機能しませんが、外れ値比率が を下回ると非常にうまく機能します)。 ${\displaystyle 95\%}$${\displaystyle 50\%}$

Parraらは、外れ値の対応を刈り込みながら、内値の対応を保証するために幾何学的制約を使用する、保証付き外れ値除去 (GORE) と呼ばれる方法を提案しました。^[20] GORE は外れ値比率を大幅に削減できることが示されており、RANSAC または BnB を使用したコンセンサス最大化のパフォーマンスを大幅に向上させることができます。 Yang と Carlone は、元の測定値セットからペアワイズ変換および回転不変測定値 (TRIM) を構築し、 3D ポイントをノードとするグラフのエッジとして TRIM を埋め込むことを提案しました。内値はスケールに関してペアワイズで一貫しているため、グラフ内でクリークを形成する必要があります。したがって、グラフの最大クリークを計算するための効率的なアルゴリズムを使用すると、内値を見つけて外れ値を効果的に刈り込むことができます。 ^[4]最大クリークに基づく外れ値除去方法は、実際のポイント セット レジストレーション問題でも非常に有用であることが示されています。^[19]同様の外れ値除去のアイデアは、Parraらによっても提案されている。 [ ^28]

M推定は、 （ cb.2 ）の最小二乗目的関数を外れ値の影響を受けにくい堅牢なコスト関数に置き換えます。正式には、M推定は次の問題を解決しようとします。

ここで、はロバストなコスト関数の選択を表します。 を選択すると、 ( cb.2 )の最小二乗推定が回復されることに注意してください。 一般的なロバストなコスト関数には、 - ノルム損失、Huber 損失、^[29] Geman-McClure 損失^[30]および切り捨て最小二乗損失があります。^{[19] [8] [4]} M 推定は、ロボット工学とコンピューター ビジョンにおけるロバスト推定の最も一般的なパラダイムの 1 つです。^{[31] [32]}ロバストな目的関数は通常非凸であるため (たとえば、切り捨て最小二乗損失と最小二乗損失)、非凸 M 推定を解決するためのアルゴリズムは通常、局所最適化 に基づいています。局所最適化では、最初に初期推定値が提供され、次に目的関数を減少させ続けるために変換が反復的に調整されます。 初期推定値が大域的最小値に近い場合は局所最適化がうまく機能する傾向がありますが、初期化が不十分な場合は局所最小値で止まってしまう傾向もあります。 ${\displaystyle \rho (\cdot )}$${\displaystyle \rho (x)=x^{2}}$${\displaystyle \ell _{1}}$

段階的非凸性（GNC）は、初期化なしで非凸最適化問題を解くための汎用フレームワークです。これは、初期の視覚および機械学習アプリケーションで成功を収めています。^{[33] [34]} GNCの背後にある重要なアイデアは、簡単な凸問題から始めて、難しい非凸問題を解決することです。具体的には、与えられた堅牢なコスト関数に対して、ハイパーパラメータを持つ代理関数を構築し、それを調整することで、代理関数の非凸性を徐々に増加させて、ターゲット関数に収束させることができます。^{[34] [35]}したがって、ハイパーパラメータの各レベルで、次の最適化が解決されます。${\displaystyle \rho (\cdot )}$${\displaystyle \rho _{\mu }(\cdot )}$${\displaystyle \mu }$${\displaystyle \rho _{\mu }(\cdot )}$${\displaystyle \rho (\cdot )}$${\displaystyle \mu }$

ブラックとランガラジャンは、各最適化の目的関数（cb.6 ）が、重み付き最小二乗法の合計と、各測定ペアにおける最適化の信頼性を決定する重みに関するいわゆる外れ値プロセス関数に双対化できることを証明しました。^[33]ブラック-ランガラジャン双対性とGeman-McClure関数に合わせて調整されたGNCを使用して、Zhouらは、対応関係の外れ値に対して堅牢な高速グローバル登録アルゴリズムを開発しました。 ^[30]最近では、ヤンらは、 GNC（Geman-McClure関数と切り捨て最小二乗関数に合わせて調整）とブラック-ランガラジャン双対性の併用により、点群やメッシュ登録などの堅牢な登録問題に対する汎用ソルバーを実現できることを示しました。^[35]${\displaystyle 80\%}$

上で述べた堅牢な登録アルゴリズムのほとんどは（最悪の場合でも指数時間で実行される BnB アルゴリズムを除く） 、パフォーマンス保証が付いていません。つまり、これらのアルゴリズムは予告なしに完全に誤った推定値を返す可能性があります。したがって、これらのアルゴリズムは、自動運転などの安全性が重要なアプリケーションには適していません。

ごく最近、ヤンらは、最初の証明可能な堅牢なレジストレーションアルゴリズムである、切り捨て最小二乗推定と半定値緩和（TEASER）を開発しました。^[19]点群レジストレーションの場合、TEASERは変換の推定値を出力するだけでなく、与えられた推定値の最適性を定量化します。TEASERは、次の切り捨て最小二乗（TLS）推定値を採用しています。

これは、TLS ロバスト コスト関数 を選択することで得られます。ここで、 は、許容される残差の最大値を決定する定義済み定数です。TLS 目的関数には、インライア対応 ( ) に対して通常の最小二乗ペナルティが適用され、アウトライア対応 ( ) に対してはペナルティは適用されず、アウトライアが破棄されるという特性があります。TLS 最適化 ( cb.7 ) がグローバル最適に解決された場合、インライア対応のみで Horn 法を実行するのと同等になります。 ${\displaystyle \rho (x)=\min(x^{2},{\bar {c}}^{2})}$${\displaystyle {\bar {c}}^{2}}$${\displaystyle \Vert s_{i}-lRm_{i}-t\Vert _{2}^{2}/\sigma _{i}^{2}<{\bar {c}}^{2}}$${\displaystyle \Vert s_{i}-lRm_{i}-t\Vert _{2}^{2}/\sigma _{i}^{2}>{\bar {c}}^{2}}$

しかし、( cb.7 )を解くのは、その組み合わせ的な性質上、かなり困難です。 TEASERは( cb.7 )を次のように解きます。(i) スケール、回転、および変換の推定を切り離して別々に解くことができるように不変の測定を構築します。これは、元のHornの方法に触発された戦略です。(ii) 3つのサブ問題のそれぞれに同じTLS推定が適用されます。スケールTLS問題は、適応投票と呼ばれるアルゴリズムを使用して正確に解くことができ、回転TLS問題は、半正定値計画法(SDP)に緩和することができ、その緩和は実際には正確であり、^[8]大量の外れ値があっても正確です。変換TLS問題は、コンポーネントごとの適応投票を使用して解決できます。GNCを活用した高速な実装は、こちらでオープンソース化されています。実際には、TEASERは外れ値対応以上を許容し、数ミリ秒で実行されます。 ${\displaystyle 99\%}$

TEASERの開発に加えて、ヤンらは、点群データに対するいくつかの穏やかな条件下では、TEASERの推定変換が真の変換からの誤差を制限していることも証明しました。^[19]

反復最近点（ICP）アルゴリズムは、BeslとMcKayによって導入されました。^[36] このアルゴリズムは、（i）変換が与えられた場合、内のすべての点に対して、内の最近点を見つける；（ii）対応が与えられた場合、最小二乗問題（cb.2 ）を解くことによって最適な剛体変換を見つける、という動作を交互に行うことで、反復的に剛体レジストレーションを実行します。したがって、の初期姿勢がに十分に近い場合に、このアルゴリズムは最適に機能します。疑似コードでは、基本的なアルゴリズムは次のように実装されています。 ${\displaystyle {\mathcal {S}}}$${\displaystyle {\mathcal {M}}}$${\displaystyle {\mathcal {M}}}$${\displaystyle {\mathcal {S}}}$

アルゴリズム ICP( M , S ) 
    θ  := θ 0登録されてい
    ない場合:
         X  := ∅ 
        m i ∊ T ( M , θ )の場合:  
            ŝ i := S内のm i
            に 最も近い点X  := X + ⟨ m i , ŝ i ⟩ 
        θ  := least_squares( X )
     θを返す 

ここで、関数はホーン^[16]とアラン^[17]による閉形式の解を使用して、各ペアの距離を最小化する最小二乗least_squares最適化を実行します。${\displaystyle \langle m_{i},{\hat {s}}_{i}\rangle }$

レジストレーションのコスト関数は、内のすべての点に最も近い 内の点を見つけることに依存するため、アルゴリズムの実行中に変化する可能性があります。そのため、ICP が実際に局所最適値に正確に収束することを証明することは困難です。^[37]実際、経験的には、ICP と EM-ICP はコスト関数の局所最小値に収束しません。^[37]それでも、ICP は直感的に理解でき、実装が簡単なため、最も一般的に使用されているポイント セット レジストレーション アルゴリズムです。^[37] ICP の多くのバリエーションが提案されており、ポイントの選択とマッチングから最小化戦略まで、アルゴリズムのすべてのフェーズに影響します。^{[13] [38]} たとえば、期待最大化アルゴリズムを ICP アルゴリズムに適用して EM-ICP 法を形成し、Levenberg-Marquardt アルゴリズムをICP アルゴリズムに適用して LM-ICP 法を形成します。^[12]${\displaystyle {\mathcal {S}}}$${\displaystyle {\mathcal {M}}}$

ロバストポイントマッチング（RPM）は、Goldらによって導入されました。^[39]この方法は、決定論的アニーリングと点集合間の対応のソフト割り当てを使用して位置合わせを実行します。ICPでは最近傍ヒューリスティックによって生成される対応はバイナリですが、RPMはソフト対応を使用します。ソフト対応では、任意の2点間の対応は0から1の範囲になりますが、最終的には0または1に収束します。RPMで見つかった対応は常に1対1ですが、ICPでは必ずしもそうではありません。^[14]をの番目の点、をの 番目の点とします。マッチ行列は次のように定義されます。 ${\displaystyle m_{i}}$${\displaystyle i}$${\displaystyle {\mathcal {M}}}$${\displaystyle s_{j}}$${\displaystyle j}$${\displaystyle {\mathcal {S}}}$ ${\displaystyle \mathbf {\mu } }$

問題は次のように定義されます。2つの点集合が与えられ、それらを最もよく関連付けるアフィン変換とマッチング行列を見つけます。 ^[39]最適な変換がわかればマッチング行列を簡単に決定できます。逆もまた同様です。ただし、RPMアルゴリズムは両方を同時に決定します。変換は、変換ベクトルと変換行列に分解できます。 ${\displaystyle {\mathcal {M}}}$${\displaystyle {\mathcal {S}}}$ ${\displaystyle T}$${\displaystyle \mathbf {\mu } }$

${\displaystyle T(m)=\mathbf {A} m+\mathbf {t} }$

2D のマトリックスは、スケール、回転、垂直および水平せん断成分の 4 つの個別のパラメータで構成されます。コスト関数は次のようになります。${\displaystyle \mathbf {A} }$${\displaystyle \lbrace a,\theta ,b,c\rbrace }$

、、の制約を受けます。この項は、マッチ行列に1がもっと多く含まれる場合にコストを下げることで、目的関数を相関が強くなるようにバイアスします。この関数は、スケール成分とシアー成分の大きな値をペナルティとしてアフィン変換を正規化します。 ${\textstyle \forall j~\sum _{i=1}^{M}\mu _{ij}\leq 1}$${\textstyle \forall i~\sum _{j=1}^{N}\mu _{ij}\leq 1}$${\textstyle \forall ij~\mu _{ij}\in \lbrace 0,1\rbrace }$${\displaystyle \alpha }$${\displaystyle g(\mathbf {A} )}$

${\displaystyle g(\mathbf {A} (a,\theta ,b,c))=\gamma (a^{2}+b^{2}+c^{2})}$

ある正規化パラメータに対して。 ${\displaystyle \gamma }$

RPM 法は、Softassign アルゴリズムを使用してコスト関数を最適化します。ここでは 1D の場合を導出します。 の変数のセットが与えられます。変数はそれぞれに関連付けられ、 となります。 目標はを最大化する を見つけることです。 これは、制御パラメータ を導入することで連続問題として定式化できます。決定論的アニーリング法では、アルゴリズムの実行中に制御パラメータが徐々に増加します。 を次のようにします。 ${\displaystyle \lbrace Q_{j}\rbrace }$${\displaystyle Q_{j}\in \mathbb {R} ^{1}}$${\displaystyle \mu _{j}}$${\displaystyle Q_{j}}$${\textstyle \sum _{j=1}^{J}\mu _{j}=1}$${\displaystyle \mathbf {\mu } }$${\textstyle \sum _{j=1}^{J}\mu _{j}Q_{j}}$${\displaystyle \beta >0}$${\displaystyle \beta }$${\displaystyle \mathbf {\mu } }$

これはソフトマックス関数として知られています。 が増加すると、方程式 ( rpm.1 )で必要なバイナリ値に近づきます。 問題は 2D の場合に一般化され、 を最大化する代わりに、以下が最大化されます。 ${\displaystyle \beta }$${\textstyle \sum _{j=1}^{J}\mu _{j}Q_{j}}$

どこ

${\displaystyle Q_{ij}=-(\lVert s_{j}-\mathbf {t} -\mathbf {A} m_{i}\rVert ^{2}-\alpha )=-{\frac {\partial \operatorname {cost} }{\partial \mu _{ij}}}}$

これは簡単ですが、制約が二重確率行列制約、および になります。そのため、式 ( rpm.3 )の分母は、 2D の場合に単純に表現することはできません。制約を満たすには、Sinkhorn の結果^[39]を使用できます。この結果によると、すべての正の要素を持つ任意の正方行列から、行と列を交互に正規化する反復プロセスによって二重確率行列が得られます。したがって、アルゴリズムは次のように記述されます。^[39]${\displaystyle \mu }$${\textstyle \forall j~\sum _{i=1}^{M}\mu _{ij}=1}$${\textstyle \forall i~\sum _{j=1}^{N}\mu _{ij}=1}$

アルゴリズム RPM2D ${\displaystyle ({\mathcal {M}},{\mathcal {S}})}$
    t  := 0 
    a , θ  b , c  := 0 
    β  := β 0 
    while β < β f : while μが収束していない間:
             // ソフト割り当てによって対応パラメータを更新// シンクホーン法を適用whileが収束していない間: //すべての行にわたって正規化することによって更新: //すべての列にわたって正規化することによって更新: // 座標降下によって姿勢パラメータを更新解析解を使用してθ を
            更新${\displaystyle {\hat {\mu }}_{ij}:=1+\epsilon }$
     
         
            ${\displaystyle Q_{ij}:=-{\frac {\partial \operatorname {cost} }{\partial \mu _{ij}}}}$
            ${\displaystyle \mu _{ij}^{0}:=\exp(\beta Q_{ij})}$
            
             ${\displaystyle {\hat {\mu }}}$
                ${\displaystyle {\hat {\mu }}}$
                ${\displaystyle {\hat {\mu }}_{ij}^{1}:={\frac {{\hat {\mu }}_{ij}^{0}}{\sum _{i=1}^{M+1}{\hat {\mu }}_{ij}^{0}}}}$
                ${\displaystyle {\hat {\mu }}}$
                ${\displaystyle {\hat {\mu }}_{ij}^{0}:={\frac {{\hat {\mu }}_{ij}^{1}}{\sum _{j=1}^{N+1}{\hat {\mu }}_{ij}^{1}}}}$
            解析解を使用してtを
            更新するニュートン法を使用してa、b、cを
            更新し、a、b、c、θ、tを返します。
        ${\displaystyle \beta :=\beta _{r}\beta }$
        ${\displaystyle \gamma :={\frac {\gamma }{\beta _{r}}}}$
     

ここで、決定論的アニーリング制御パラメータは、最初は に設定され、最大値 に達するまで係数 ずつ増加します。 の制約が不等式であるため、正規化ステップでの合計は、 と ではなく とになります。そのため、番目の要素と番目の要素はスラック変数です。 ${\displaystyle \beta }$${\displaystyle \beta _{0}}$${\displaystyle \beta _{r}}$${\displaystyle \beta _{f}}$${\displaystyle M+1}$${\displaystyle N+1}$${\displaystyle M}$${\displaystyle N}$${\displaystyle \mu }$${\displaystyle M+1}$${\displaystyle N+1}$

このアルゴリズムは、3次元以上の点集合にも拡張できる。対応行列の制約は、3次元の場合も2次元の場合も同じである。したがって、アルゴリズムの構造は変わらないが、主な違いは回転行列と並進行列の解き方である。^[39]${\displaystyle \mathbf {\mu } }$

ChuiとRangarajanによる薄板スプラインロバストポイントマッチング（TPS-RPM）アルゴリズムは、RPM法を拡張して、変換を薄板スプラインとしてパラメータ化することで非剛体レジストレーションを実行します。^[14] しかし、薄板スプラインのパラメータ化は3次元でのみ存在するため、この方法は4次元以上の問題に拡張することはできません。

ポイントセットレジストレーションのカーネル相関（KC）アプローチは、TsinとKanadeによって導入されました。^[37] ICPと比較すると、KCアルゴリズムはノイズの多いデータに対してより堅牢です。すべてのモデルポイントに対して最も近いシーンポイントのみが考慮されるICPとは異なり、ここではすべてのシーンポイントがすべてのモデルポイントに影響を与えます。^[37]そのため、これは多重リンクレジストレーションアルゴリズムです。あるカーネル関数 に対して、 2つのポイントのカーネル相関は次のように定義されます。^[37]${\displaystyle K}$${\displaystyle KC}$${\displaystyle x_{i},x_{j}}$

点集合の登録に選択されるカーネル関数は 、通常、対称かつ非負のカーネルであり、パルゼンウィンドウ密度推定で使用されるものと同様である。ガウスカーネルは、その単純さから一般的に使用されるが、エパネチニコフカーネルやトライキューブカーネルなどの他のカーネルを代用することもできる。^[37]点集合全体のカーネル相関は、集合内のすべての点と集合内の他のすべての点とのカーネル相関の合計として定義される。^[37]${\displaystyle K}$${\displaystyle {\mathcal {\chi }}}$

ポイント セットの KC の対数は、定数倍の範囲内で情報エントロピーに比例します。KC はポイント セットの「コンパクトさ」の尺度であることに注意してください。ポイント セット内のすべてのポイントが同じ場所にある場合、KC は大きな値に評価されます。変換パラメータに対するポイント セット登録アルゴリズムのコスト関数は、次のように定義されます。 ${\displaystyle \theta }$

いくつかの代数操作により次の結果が得られます。

この式は、が に依存しないことを観察することで簡略化されます。さらに、固定位置合わせを想定すると、が変更されても は不変です。これは、すべての点のペア間のユークリッド距離が固定変換 の下で同じままであるためです。したがって、上記の式は次のように書き直すことができます。 ${\displaystyle KC({\mathcal {S}})}$${\displaystyle \theta }$${\displaystyle KC(T({\mathcal {M}},\theta ))}$${\displaystyle \theta }$

カーネル密度推定値は次のように定義されます。

${\displaystyle P_{\mathcal {M}}(x,\theta )={\frac {1}{M}}\sum _{m\in {\mathcal {M}}}K(x,T(m,\theta ))}$

${\displaystyle P_{\mathcal {S}}(x)={\frac {1}{N}}\sum _{s\in {\mathcal {S}}}K(x,s)}$

コスト関数は、2 つのカーネル密度推定値の相関として示されます。

コスト関数を確立したら、アルゴリズムは単純に勾配降下法を使用して最適な変換を見つけます。すべての反復でコスト関数を最初から計算するのは計算コストが高いため、コスト関数方程式の離散バージョン ( kc.6 ) が使用されます。カーネル密度推定値はグリッド ポイントで評価し、ルックアップ テーブルに格納できます。ICP や関連方法とは異なり、最近傍を見つける必要がないため、KC アルゴリズムの実装は比較的簡単です。 ${\displaystyle P_{\mathcal {M}},P_{\mathcal {S}}}$

ノイズの多い2Dおよび3Dポイントセットに対するICPおよびEM-ICPと比較すると、KCアルゴリズムはノイズの影響を受けにくく、より頻繁に正しい登録が得られます。^[37]

カーネル密度推定値はガウス分布の合計であるため、ガウス混合モデル（GMM）として表すことができます。^{[40] JianとVemuriは、KCレジストレーションアルゴリズムのGMMバージョンを使用して、}薄板スプラインによってパラメータ化された非剛性レジストレーションを実行します。

コヒーレントポイントドリフト（CPD）は、ミロネンコとソングによって導入されました。^{[13] [41]このアルゴリズムは、GMM KC法と同様に、ポイントセットの位置合わせに確率的アプローチを採用しています。}薄板スプライン変換モデル を前提とする非剛体レジストレーションの以前のアプローチとは異なり、CPDは使用される変換モデルに関して不可知です。ポイントセットは、ガウス混合モデル（GMM）の重心を表します。2つのポイントセットが最適に位置合わせされている場合、対応は、特定のデータポイントのGMM事後確率の最大値です。ポイントセットの位相構造を維持するために、GMMの重心はグループとしてコヒーレントに移動するように強制されます。期待値最大化アルゴリズムを使用して、コスト関数を最適化します。^[13]${\displaystyle {\mathcal {M}}}$

にM個の点があり、にN個の点があるとします。点sのGMM確率密度関数は次のようになります。 ${\displaystyle {\mathcal {M}}}$${\displaystyle {\mathcal {S}}}$

ここで、D次元では、点を中心とするガウス分布です。 ${\displaystyle p(s|i)}$${\displaystyle m_{i}\in {\mathcal {M}}}$

${\displaystyle p(s|i)={\frac {1}{(2\pi \sigma ^{2})^{D/2}}}\exp {\left(-{\frac {\lVert s-m_{i}\rVert ^{2}}{2\sigma ^{2}}}\right)}}$

メンバーシップ確率はすべてのGMMコンポーネントで等しい。一様分布の重みは と表される。混合モデルは次のようになる。 ${\displaystyle P(i)={\frac {1}{M}}}$${\displaystyle w\in [0,1]}$

GMMの重心は、尤度を最大化することによって推定されたパラメータのセットによって再パラメータ化されます。これは、負の対数尤度関数を最小化することと同じです。 ${\displaystyle \theta }$

ここで、データは独立かつ同一に分布していると仮定します。2つの点と間の対応確率は、データ点が与えられたときのGMM重心の 事後確率として定義されます。${\displaystyle m_{i}}$${\displaystyle s_{j}}$

${\displaystyle P(i|s_{j})={\frac {P(i)p(s_{j}|i)}{p(s_{j})}}}$

期待値最大化(EM) アルゴリズムは、と を見つけるために使用されます。EM アルゴリズムは 2 つのステップで構成されます。まず、E ステップまたは推定ステップで、パラメータの値 (「古い」パラメータ値) を推測し、ベイズの定理を使用して混合成分の事後確率分布を計算します。次に、M ステップまたは最大化ステップで、完全な負の対数尤度関数、つまりコスト関数の期待値を最小化することによって、「新しい」パラメータ値を見つけます。 ${\displaystyle \theta }$${\displaystyle \sigma ^{2}}$${\displaystyle P^{\text{old}}(i,s_{j})}$

およびに依存しない定数を無視すると、式( cpd.4 )は次のように表すことができます。 ${\displaystyle \theta }$${\displaystyle \sigma }$

どこ

${\displaystyle N_{\mathbf {P} }=\sum _{j=0}^{N}\sum _{i=0}^{M}P^{\text{old}}(i|s_{j})\leq N}$

の場合にのみ成り立ちます。以前のパラメータ値を使用して計算されたGMMコンポーネントの事後確率は次のようになります。 ${\displaystyle N=N_{\mathbf {P} }}$${\displaystyle w=0}$${\displaystyle P^{\text{old}}}$

式（ cpd.5 ）のコスト関数を最小化すると、式（cpd.3 ）の負の対数尤度関数Eが、すでに局所最小値に達していない限り、必然的に減少する。^[13]したがって、アルゴリズムは次の疑似コードを使用して表現することができ、点集合とがそれぞれ行列ととして表される：^[13]${\displaystyle {\mathcal {M}}}$${\displaystyle {\mathcal {S}}}$${\displaystyle M\times D}$${\displaystyle N\times D}$${\displaystyle \mathbf {M} }$${\displaystyle \mathbf {S} }$

アルゴリズム CPD ${\displaystyle ({\mathcal {M}},{\mathcal {S}})}$
    θ  := θ 0
    初期化 0 ≤ w ≤ 1登録されていない
    間:
         // Eステップ、i ∊ [1, M ]およびj ∊ [1, N ]に対してPを計算: // Mステップ、最適変換を解く{ θ , σ 2 } := solve ( S , M , P ) return θ${\displaystyle \sigma ^{2}:={\frac {1}{DNM}}\sum _{j=1}^{N}\sum _{i=1}^{M}\lVert s_{j}-m_{i}\rVert ^{2}}$
    
         
            ${\displaystyle p_{ij}:={\frac {\exp \left(-{\frac {1}{2\sigma ^{2}}}\lVert s_{j}-T(m_{i},\theta )\rVert ^{2}\right)}{\sum _{k=1}^{M}\exp \left(-{\frac {1}{2\sigma ^{2}}}\lVert s_{j}-T(m_{k},\theta )\rVert ^{2}\right)+(2\pi \sigma ^{2})^{\frac {D}{2}}{\frac {w}{1-w}}{\frac {M}{N}}}}}$
        
        
     

ここで、ベクトルは 1 の列ベクトルです。関数は、実行される位置合わせの種類によって異なります。たとえば、リジッド位置合わせでは、出力はスケールa、回転行列、および変換ベクトルです。パラメータは、次のタプルとして記述できます。 ${\displaystyle \mathbf {1} }$solve${\displaystyle \mathbf {R} }$${\displaystyle \mathbf {t} }$${\displaystyle \theta }$

${\displaystyle \theta =\lbrace a,\mathbf {R} ,\mathbf {t} \rbrace }$

これは1に初期化され、単位行列とゼロの列ベクトルになります。

${\displaystyle \theta _{0}=\lbrace 1,\mathbf {I} ,\mathbf {0} \rbrace }$

整列されたポイント セットは次のとおりです。

${\displaystyle T(\mathbf {M} )=a\mathbf {M} \mathbf {R} ^{T}+\mathbf {1} \mathbf {t} ^{T}}$

剛体レジストレーションの関数solve_rigidは次のように記述でき、代数の導出はMyronenkoの2010年の論文で説明されている。^[13]

solve_rigid ( S , M , P ) 
    N P  := 1 T P1 
    U , V  := svd ( A ) // Aの特異値分解= UΣV T C  := diag(1, …, 1, det( UV T )) // diag( ξ )はベクトル ξ から形成された対角行列ですR  := UCV T // trは行列のトレースですt  := μ s − a R μ m return { a , R , t }, σ 2${\displaystyle \mu _{s}:={\frac {1}{N_{\mathbf {P} }}}\mathbf {S} ^{T}\mathbf {P} ^{T}\mathbf {1} }$
    ${\displaystyle \mu _{m}:={\frac {1}{N_{\mathbf {P} }}}\mathbf {M} ^{T}\mathbf {P} \mathbf {1} }$
    ${\displaystyle {\hat {\mathbf {S} }}:=\mathbf {S} -\mathbf {1} \mu _{s}^{T}}$
    ${\displaystyle {\hat {\mathbf {M} }}:=\mathbf {M} -\mathbf {1} \mu _{m}^{T}}$
    ${\displaystyle \mathbf {A} :={\hat {\mathbf {S} ^{T}}}\mathbf {P} ^{T}{\hat {\mathbf {M} }}}$
      
     
    
    ${\displaystyle a:={\frac {\operatorname {tr} (\mathbf {A} ^{T}\mathbf {R} )}{\operatorname {tr} (\mathbf {{\hat {\mathbf {M} }}^{T}\operatorname {diag} (\mathbf {P} \mathbf {1} ){\hat {\mathbf {M} }}} )}}}$   
    
    ${\displaystyle \sigma ^{2}:={\frac {1}{N_{\mathbf {P} }D}}(\operatorname {tr} (\mathbf {{\hat {\mathbf {S} }}^{T}\operatorname {diag} (\mathbf {P} ^{T}\mathbf {1} ){\hat {\mathbf {S} }}} )-a\operatorname {tr} (\mathbf {A} ^{T}\mathbf {R} ))}$
    

アフィン登録の場合、目標は固定変換ではなくアフィン変換を見つけることであり、出力はアフィン変換行列と、位置合わせされた点セットが次のようになるような 変換です。${\displaystyle \mathbf {B} }$${\displaystyle \mathbf {t} }$

${\displaystyle T(\mathbf {M} )=\mathbf {M} \mathbf {B} ^{T}+\mathbf {1} \mathbf {t} ^{T}}$

剛体レジストレーションの関数solve_affineは次のように記述でき、代数の導出はMyronenkoの2010年の論文で説明されている。^[13]

アフィン解を解く( S , M , P ) 
    N P  := 1 T P1 
    t  := μ s − B μ mを返す{ B , t }, σ 2${\displaystyle \mu _{s}:={\frac {1}{N_{\mathbf {P} }}}\mathbf {S} ^{T}\mathbf {P} ^{T}\mathbf {1} }$
    ${\displaystyle \mu _{m}:={\frac {1}{N_{\mathbf {P} }}}\mathbf {M} ^{T}\mathbf {P} \mathbf {1} }$
    ${\displaystyle {\hat {\mathbf {S} }}:=\mathbf {S} -\mathbf {1} \mu _{s}^{T}}$
    ${\displaystyle {\hat {\mathbf {M} }}:=\mathbf {M} -\mathbf {1} \mu _{m}^{T}}$
    ${\displaystyle \mathbf {B} :=({\hat {\mathbf {S} }}^{T}\mathbf {P} ^{T}{\hat {\mathbf {M} }})({\hat {\mathbf {M} }}^{T}\operatorname {diag} (\mathbf {P} \mathbf {1} ){\hat {\mathbf {M} }})^{-1}}$
    
    ${\displaystyle \sigma ^{2}:={\frac {1}{N_{\mathbf {P} }D}}(\operatorname {tr} ({\hat {\mathbf {S} }}^{T}\operatorname {diag} (\mathbf {P} ^{T}\mathbf {1} ){\hat {\mathbf {S} }})-\operatorname {tr} ({\hat {\mathbf {S} }}^{T}\mathbf {P} ^{T}{\hat {\mathbf {M} }}\mathbf {B} ^{T}))}$
    

変分法を用いて導出されたパラメータ化を用いた非剛体レジストレーションでCPDを使用することも可能である。^[13]

ガウス分布の和は高速ガウス変換（FGT）を使って線形時間で計算できる。 ^[13]したがって、CPDの時間計算量は であり、漸近的には方法よりもはるかに高速である。^[13]${\displaystyle O(M+N)}$${\displaystyle O(MN)}$

ベイズコヒーレントポイントドリフト（BCPD）と呼ばれるコヒーレントポイントドリフトの変形は、ポイントセットレジストレーションのベイズ定式化によって導き出されました。 ^[42] BCPDはCPDに比べていくつかの利点があります。たとえば、（1）非剛性および剛性のレジストレーションを単一のアルゴリズムで実行できます。（2）動きのコヒーレンスを定義するグラム行列のガウス性に関係なくアルゴリズムを高速化できます。（3）外れ値分布の定義がより合理的であるため、アルゴリズムは外れ値に対してより堅牢です。さらに、ベイズ定式化では、動きのコヒーレンスは変位ベクトルの事前分布を通じて導入され、動きのコヒーレンスを制御するチューニングパラメータ間に明確な違いをもたらします。BCPDは、（1）ポイントセットのダウンサンプリング、（2）ダウンサンプリングされたポイントセットのレジストレーション、および（3）変形フィールドの補間からなる3段階の手順であるBCPD++と呼ばれる方法によってさらに高速化されました。 ^[43] この方法は、登録精度を維持しながら1000万点以上の点群を登録することができる。

コヒーレントポイントドリフトの変形として、剛性ポイントクラウドレジストレーション用のローカルサーフェスジオメトリ付きCPD（LSG-CPD）と呼ばれるものがあります。^[44]この方法は、ローカルサーフェスの平坦度に基づいて、ポイントツーポイントペナルティに加えて、ポイントツープレーンペナルティのさまざまなレベルを適応的に追加します。これにより、元のCPDの等方性共分散ではなく、異方性共分散を持つGMMコンポーネントが生成されます。^[13]異方性共分散行列は次のようにモデル化されます。

どこ

${\displaystyle \Sigma _{m}}$はターゲットセットのm番目の点の異方性共分散行列です。は同じ点に対応する法線ベクトルです。 は正規化子として機能し、問題を不適切性から引き離す単位行列です。はペナルティ係数（修正シグモイド関数）であり、局所表面の平坦度に応じて異なるレベルの点対平面ペナルティを追加するように適応的に設定されます。 これは、m番目のターゲットポイントの近傍内の表面変化^[45]を評価することによって実現されます。はペナルティの上限です。 ${\displaystyle \mathbf {n} _{m}}$${\displaystyle \mathbf {I} }$${\displaystyle \alpha _{m}}$${\displaystyle \kappa _{m}}$${\displaystyle \alpha _{max}}$

ポイントクラウドのレジストレーションは、最大尤度推定 (MLE) 問題として定式化され、期待値最大化 (EM) アルゴリズムで解決されます。E ステップでは、対応計算が単純な行列操作に作り直され、GPU で効率的に計算されます。M ステップでは、行列リー群の制約なしの最適化が、レジストレーションの剛体変換を効率的に更新するように設計されています。ローカル幾何学的共分散を利用することで、この方法は、ベースライン CPD と比較して、精度とノイズおよび外れ値に対する堅牢性の点で優れたパフォーマンスを示します。^[46] GPU による対応計算の高速化により、実行時パフォーマンスの向上が期待されます。LSG-CPD の実装は、こちらでオープンソース化されています。

このアルゴリズムは、ソナー画像のレジストレーションに対応するために、2013 年に H. Assalih によって導入されました。^[47]この種の画像はノイズが多くなりがちなので、マッチングするポイント セットには多くの外れ値が含まれることが予想されます。SCS は外れ値に対して高い堅牢性を提供し、外れ値がある場合でも ICP や CPD のパフォーマンスを上回ることができます。SCS は高次元空間での反復最適化を使用せず、確率的でもスペクトルでもありません。SCS は、剛体変換と非剛体変換をマッチングでき、ターゲット変換が 3 から 6 の自由度の間にあるときに最高のパフォーマンスを発揮します。

• ポイントフィーチャマッチング
• 点集合三角測量
• 正規分布の変換

• 薄板スプラインの堅牢なポイントマッチングのリファレンス実装
• カーネル相関ポイントセット登録のリファレンス実装
• コヒーレントポイントドリフトのリファレンス実装
• ICPバリアントのリファレンス実装
• ベイジアンコヒーレントポイントドリフトのリファレンス実装
• LSG-CPDのリファレンス実装