Preconditioned Single-step Transforms for Non-rigid ICP (Eurographics 2025)

Non-rigid(非剛体)ICPは対応探索と形状更新を反復する手法であり、形状更新では一般に平滑化正則化を含む線形最小二乗問題を解いてターゲット形状へ近づける。
しかし初期反復では対応が不正確であるため、各反復で線形最小二乗を高精度に解くことは計算的に非効率になる。
本研究は、各反復で厳密解を求める代わりに、共役勾配法（CG）を1ステップのみ実行して変形を更新するsingle-step transform枠組みを提案する。

Transcript

1. Preconditioned Single-step Transforms for Non-rigid ICP Eurographics 2025 Yucheol Jung, Hyomin Kim, Hyejeong Yoon, Seungyong Lee エクサウィザーズ ⾦⽥綾乃
2. Non-Rigid ICPって︖ Non-Rigid Iterative Closest Points Source を⾮剛体(Non-Rigid)に変形して Target 形状に⼀致させる Target Source Non-Rigid(⾮剛体)変形
3. Source を⾮剛体(Non-Rigid)に変形して Target 形状に⼀致させる Target Source Global alignment Non-Rigid(⾮剛体)変形 対応点を更新 (Cを更新) Source形状を⾮剛体変形で更新 (Ax = b) 各々の更新を 繰り返して最適化 Non-Rigid ICPって︖ Non-Rigid Iterative Closest Points Iterative Closest Points
4. Source を⾮剛体(Non-Rigid)に変形して Target 形状に⼀致させる Target Source Global alignment Non-Rigid(⾮剛体)変形 対応点を更新 (Cを更新) Source形状を⾮剛体変形で更新 (Ax = b) 各々の更新を 繰り返して最適化 今回の論⽂はこちらの最適化を * 多数の CG ステップで厳密に解く代わりに * Sobolev–Jacobi 前処理付きの 1step PCG （Single-step transform） ⾼速化かつ安定化を実現 Non-Rigid ICPって︖ Non-Rigid Iterative Closest Points Iterative Closest Points
5. Non-Rigid ICPって︖ Source を⾮剛体(Non-Rigid)に変形して Target 形状に⼀致させる Target Source Global alignment 対応点を更新 (Cを更新) Source形状を⾮剛体変形で更新 (Ax = b) 各々の更新を 繰り返して最適化 今回の論⽂はこちらの最適化を * 多数の CG ステップで厳密に解く代わりに * Sobolev–Jacobi 前処理付きの 1step PCG （Single-step transform） ⾼速化かつ安定化を実現 Iterative Closest Points
6. Source形状を⾮剛体変形で更新 Ax=b 𝐸!"! = 𝐸#$#% + 𝜆%&' 𝐸%&' + 𝜆()* 𝐸()* 𝐸!"!# 𝑆$ , 𝑇 = 1 |𝑐$ | ) (&,#)∈*! ( 𝑣& 𝑆$ − 𝑣# ′) / 𝑛# ′ " 𝐸#+& 𝑆$ , 𝑇 = 1 |𝑐#+& | ) (&,#)∈*"#$ 𝑣& (𝑆$ ) − 𝑣# ′ " 𝐸,-. 𝑆$ = 1 𝑁 ) & / ) +∈/(&) 𝑥& (𝑆$ ) − 𝑥+ (𝑆$ ) " x に関して線形関数の⼆乗の⾜し合わせになっている 𝐸!"! 𝑥 = . # 𝑟# 𝑥 $ = 𝑟 𝑥 $ = 𝐽𝑥 − 𝑑 $ ①Point-to-plane 誤差 ②Landmark 誤差 ③Laplacian 誤差 𝐴𝑥 = 𝑏 𝐽3𝐽𝑥 = 𝐽3𝑑 * 𝑣& 𝑆$ , 𝑥& (𝑆$ )は頂点座標ベクトル x の線形結合で書ける 𝐴 ∈ ℝ0×0 SPD 隣り合う頂点の位置（あるいは変位）が あまりバラバラにならないようにする“平滑化項” 今の変形後の Source 頂点 𝑣! 𝑆" が、 Target の局所平⾯から法線⽅向にどれだけ離れているか
7. Non-Rigid ICP における Ax = b の役割と課題 対応点を更新 (Cを更新) Source形状を⾮剛体変形で更新 (Ax = b) 各々の更新を 繰り返して最適化 通常であれば • 直接法（Cholesky） • PCG（前処理付き共役勾配…. 等で厳密解を求めていく 1. 対応点が更新されるごとに𝑨, 𝒃も更新される 2. 特に早い段階では対応が定まっておらず、 そのイテレーションだけ厳密解を求めても意味が薄い 今の対応に対して“そこそこ良い”⽅向に動いてくれれば⼗分 iteration
8. PCGのプロセス A𝑥 = b 𝐴 ∈ ℝ0×0 SPD 𝑟# = 𝑓 − 𝐴𝑥 𝑧# = 𝑀%&𝑟# 𝛽# = 𝑟# '𝑧# 𝑟#%& ' 𝑧#%& 𝑝# = 𝑧# + 𝛽# 𝑝#%& 𝛼# = 𝑟# '𝑧# 𝑝# ' 𝐻𝑝# 𝑥#(& = 𝑥# + 𝛼# 𝑝# 1. 残差 2. 前処理 4. 探索⽅向 5. ステップ⻑ 6. 位置更新 3.共役⽅向の係数 通常は k=0..K 収束まで何回も繰り返す 𝐵 = 𝑀%&𝐴 の条件数 𝜅 𝐵 = 𝜆+23 𝜆+$0 が⼩さいほど早く収束する → 良い前処理⾏列 M 対応点を更新 (Cを更新) Source形状を⾮剛体変形で更新 (Ax = b)
9. PCGのプロセス A𝑥 = b 𝐴 ∈ ℝ0×0 SPD 𝑟# = 𝑓 − 𝐴𝑥 𝑧# = 𝑀%&𝑟# 𝛽# = 𝑟# '𝑧# 𝑟#%& ' 𝑧#%& 𝑝# = 𝑧# + 𝛽# 𝑝#%& 𝛼# = 𝑟# '𝑧# 𝑝# ' 𝐻𝑝# 𝑥#(& = 𝑥# + 𝛼# 𝑝# 1. 残差 2. 前処理 4. 探索⽅向 5. ステップ⻑ 6. 位置更新 3.共役⽅向の係数 通常は k=0..K 収束まで何回も繰り返す 𝐵 = 𝑀%&𝐴 の条件数 𝜅 𝐵 = 𝜆+23 𝜆+$0 が⼩さいほど早く収束する → 良い前処理⾏列 M 対応点を更新 (Cを更新) Source形状を⾮剛体変形で更新 (Ax = b) 今このステップだけに対して “厳密解” を取っても 次のステップ(Cを更新)で捨てられる
10. PCGのプロセス A𝑥 = b 𝐴 ∈ ℝ0×0 SPD 𝑟# = 𝑓 − 𝐴𝑥 𝑧# = 𝑀%&𝑟# 𝛽# = 𝑟# '𝑧# 𝑟#%& ' 𝑧#%& 𝑝# = 𝑧# + 𝛽# 𝑝#%& 𝛼# = 𝑟# '𝑧# 𝑝# ' 𝐻𝑝# 𝑥#(& = 𝑥# + 𝛼# 𝑝# 1. 残差 2. 前処理 4. 探索⽅向 5. ステップ⻑ 6. 位置更新 3.共役⽅向の係数 通常は k=0..K 収束まで何回も繰り返す 𝐵 = 𝑀%&𝐴 の条件数 𝜅 𝐵 = 𝜆+23 𝜆+$0 が⼩さいほど早く収束する → 良い前処理⾏列 M 対応点を更新 (Cを更新) Source形状を⾮剛体変形で更新 (Ax = b) 今このステップだけに対して “厳密解” を取っても 次のステップ(Cを更新)で捨てられる NR-ICPの構造に特化して1回のPCGステップでまあまあ良い更新が出 せればいいのでは
11. PCGのプロセス A𝑥 = b 𝐴 ∈ ℝ0×0 SPD 𝑟# = 𝑓 − 𝐴𝑥 𝑧# = 𝑀%&𝑟# 𝛽# = 𝑟# '𝑧# 𝑟#%& ' 𝑧#%& 𝑝# = 𝑧# + 𝛽# 𝑝#%& 𝛼# = 𝑟# '𝑧# 𝑝# ' 𝐻𝑝# 𝑥#(& = 𝑥# + 𝛼# 𝑝# 1. 残差 2. 前処理 4. 探索⽅向 5. ステップ⻑ 6. 位置更新 3.共役⽅向の係数 通常は k=0..K 収束まで何回も繰り返す 𝐵 = 𝑀%&𝐴 の条件数 𝜅 𝐵 = 𝜆+23 𝜆+$0 が⼩さいほど早く収束する → 良い前処理⾏列 M 対応点を更新 (Cを更新) Source形状を⾮剛体変形で更新 (Ax = b) 今このステップだけに対して “厳密解” を取っても 次のステップ(Cを更新)で捨てられる NR-ICPの構造に特化して1回のPCGステップでまあまあ良い更新が出 せればいいのでは →そのために適切な前処理を設計していこう
12. ⼀回の更新で満たしたい条件 条件 前処理M で実現したいこと 1. エネルギーを⼗分減らす Bの条件数を下げて,1ステップでも residual が落ちるようにする 2. 滑らかで⼤域的な変形 Laplacian を使って残差を拡散し，⾼周波ノイズを抑える 3. coarse-to-fine な挙動 iteration 番号 k に応じて Sobolev の強さを変える 4. コストが⼩さい メッシュのトポロジは固定 → Laplacian 部分を事前に分解して再利⽤ Sobolev 前処理︓ ② 滑らか かつ ③coarse-to-fineのベース Jacobi 前処理︓ ① エネルギーを⼗分減らす Adaptive Sobolev–Jacobi前処理 : ③coarse-to-fine かつ ④コスト (前処理の再利⽤)
13. Sobolev 前処理/Jacobi前処理 𝑀="> = 𝐼 + 𝑡 𝐿 # 𝐿 ∶ラプラシアン 𝑡 > 0 ∶ 拡散のスケール 𝑝 ∶ 拡散の強さ(iteration) 𝑧# = 𝑀%&𝑟# 局所的な𝑟# を、⾯全体の変形⽅向 𝑧# に変換する 𝑀+,- ./ 𝑢 𝑢 ∈ ℝ) , 𝑢* = 𝛿*+ 各頂点 𝑖に定義されたスカラー場を可視化 悪い点 : 条件数の改善は限定的で,fineへの切り替えは単体ではできない Sobolev 前処理により1点のスパイク(左)が 近傍へ拡散し，球⾯上で滑らかな分布 条件数はあまり下がらな い Sobolev 前処理︓Laplacian による平滑化 ⼀気に下がる 良い点 : 滑らか かつ coarse な⽅向へのバイアス
14. Sobolev 前処理/Jacobi前処理 𝑀="> = 𝐼 + 𝑡 𝐿 # 𝐿 ∶ラプラシアン 𝑡 > 0 ∶ 拡散のスケール 𝑝 ∶ 拡散の強さ(iteration) 𝑧# = 𝑀%&𝑟# 局所的な𝑟# を、⾯全体の変形⽅向 𝑧# に変換する 𝑀+,- ./ 𝑢 𝑢 ∈ ℝ) , 𝑢* = 𝛿*+ Sobolev 前処理により1点のスパイク(左)が 近傍へ拡散し，球⾯上で滑らかな分布 条件数はあまり下がらな い Sobolev 前処理︓Laplacian による平滑化 𝑀DEF = 𝑑𝑖𝑎𝑔 𝐴 Jacobi 前処理︓対⾓成分による逆⾏列近似 良い点 : 対数成分を取るのみなので計算コストが低い 条件数を下げてくれる(収束を早める) 悪い点 : 各成分をバラバラにスケールするのでメッシュの滑らかさは考慮しない ⼀気に下がる 各頂点 𝑖に定義されたスカラー場を可視化 良い点 : 滑らか かつ coarse な⽅向へのバイアス 悪い点 : 条件数の改善は限定的で,fineへの切り替えは単体ではできない
15. Sobolev 前処理/Jacobi前処理 𝑀="> = 𝐼 + 𝑡 𝐿 # 𝐿 ∶ラプラシアン 𝑡 > 0 ∶ 拡散のスケール 𝑝 ∶ 拡散の強さ(iteration) 𝑧# = 𝑀%&𝑟# 局所的な𝑟# を、⾯全体の変形⽅向 𝑧# に変換する 𝑀+,- ./ 𝑢 𝑢 ∈ ℝ) , 𝑢* = 𝛿*+ 良い点 : 滑らか かつ coarse な⽅向へのバイアス Sobolev 前処理により1点のスパイク(左)が 近傍へ拡散し，球⾯上で滑らかな分布 条件数はあまり下がらな い Sobolev 前処理︓Laplacian による平滑化 𝑀DEF = 𝑑𝑖𝑎𝑔 𝐴 Jacobi 前処理︓対⾓成分による逆⾏列近似 良い点 : 対数成分を取るのみなので計算コストが低い 条件数を下げてくれる(収束を早める) 悪い点 : 各成分をバラバラにスケールするのでメッシュの滑らかさは考慮しない 𝑀=D = 𝑀="> + 𝛽𝑀DEF Smoothかつ coarse 制約 条件数を下げる ⼀気に下がる Sobolev + Jacobi 前処理 (単純な⾜しあわせ) 各頂点 𝑖に定義されたスカラー場を可視化 悪い点 : 条件数の改善は限定的で,fineへの切り替えは単体ではできない
16. Adaptive Sobolev–Jacobi 前処理 𝑀E=D (𝑘) = 5 𝑀="> (𝑘) + 𝛽′ 5 𝑀DEF (k) ' 𝑀+,-(𝑘) = (𝐼 + 𝛾/ 0𝑡/𝐿)(𝐼 + 𝛾1 0𝑡1𝐿) 𝑘: Non-rigid ICPのiteration番号 𝛾4 , 𝛾" (0 < 𝛾 ≤ 1): 減衰率 k が⼩さい（序盤）ときは, 𝛾4 &, 𝛾" & ≈ 1 で : 𝑀567 = (𝐼 + 𝑡4 𝐿)(𝐼 + 𝑡" 𝐿) （強い拡散） k が進むと, ⽚⽅の項がほぼ I に近づき， : 𝑀567 = (𝐼 + ̃ 𝑡𝐿) （弱い拡散） ' 𝑀ABC(𝑘) = 𝟏 𝛾/ 0 < 𝜔 3 𝑀ABC あるkまでは Jacobi を完全に切っておき, 途中からonになる 𝟏 ⋅ ︓条件を満たすとき 1 になるインジケータ 𝑘によって求められる性質が違う 序盤(kが⼩さい) … Sobolev が強い (coarse) 後半(kが⼩さい) … Sobolev が弱い (fine) + 途中からJacobi で収束を早めていく まとめると ラプラシアン𝐿は⼀度計算すれば再利⽤可能 (コストが⼩さい) Coarse to fineが実現
17. 𝐶!(& = 𝐶 𝑥! 1. 対応点を更新する 各Non-rigid ICP イテレーションtでの⽐較 標準 Non-Rigid ICP（1回の反復 t） Single-step Transform(本論⽂) 2. 対応𝐶894から𝐴894, 𝑏894を構成 3. 𝐴894𝑥 = 𝑏894を収束するまで解く 4. 得られた𝑥894を次の形状とする 𝐶!(& = 𝐶 𝑥! 1. 対応点を更新する 2. 対応𝐶894から𝐴894, 𝑏894を構成 3. 𝐴894𝑥 = 𝑏894を１ステップだけ解く 4. 得られた𝑥894を次の形状とする 対応点を更新 (Cを更新) Source形状を⾮剛体変形で更新 (Ax = b) 前処理はAdaptive Sobolev–Jacobiを使⽤ 𝑴𝒂𝒔𝒋 (𝑡) : iteration t によって変化
18. 3DCaricShop での精度・速度⽐較 [table 1] Quasi-Newton NR-registration 3DCaricShopのベースライン Neural Deformation Pyramid Geodesic CPD correspondence-free point-set registration 収束判定の閾値ε Balance ε = 5e-6 Quality ε = 5e-7 Sourse : 顔メッシュ Sourse: 球メッシュ 古典的なNR-ICP Sphere⽤のベースライン 特に向きと滑らかさ(法 線⽅向)がとても良い 圧倒的に⾼速
19. 前処理ごとの single-step NR-ICP の性能⽐較 [table 2] 3DCaricShop のアーティストメッシュ 10 個をターゲットに使⽤ Face / Sphere テンプレートの 2 ケースについて 平均値 を算出 前処理⾏列作成 single-step NR-ICP の外側ループ Chamfer距離 (⼩さいほどいい) surface normal 誤差 (⼩さいほどいい) Jacobiは少ないiterationで収束す るが、Normal誤差が⼤きい Sob/AdaSobは収束に必要な iterationがかなり⼤きくなる Sobolev 動的 Sobolev Jacobi Jacobi+ Sobolev 動的 Jacobi+ Sobolev Ours は Jac と同じくらい速く, Sob 系と同じくらい滑らか 前処理
20. 各前処理付き single-step NR-ICP の収束曲線 [Fig 4] Jacobi単体は 最初下がるが頭打ち 序盤からスムーズに残差が 落ちる 停⽌条件 𝜀 = 0 Max ite=400 Sob/AdaSob/Sob-Jacは ゆっくり下がる
21. ⽋損メッシュをtargetにした場合(sourceは同じ姿勢のテンプレートメッシュ* [Fig 6/7]) Hole付近の三⾓形は境界に引っ張られる⼒が残るので ちょっと歪むがおおむねいい感じに ⽋損していても、もともとそこには連続な⾯があった うまくいかないケース うまくいくケース ⾜の部分が覆われることはなく崩れている ⽳の境界が形状の本来の境界／鋭いエッジに沿っている *Sourseのメッシュに関してテンプレートのみで明記なしなので別姿勢ではないと判断しました 境界に沿うような ⽋損 point-to-plane は境界の法線に沿って Source を押し付けようとする → ⽳の中⾝ ではなく ⽳の縁 に合わせに⾏く point-to-plane + Laplacian だけでは部分⼀致（partial-to-partial）問題は扱えない
22. ガウシアンノイズを与えた結果[Fig 5] 𝜎 = 10% 𝜎 = 30% Targetにガウシアンノイズを付与した結果 (Sourseは球メッシュ) results 平均エッジ⻑の ノイズ⼊り⼊⼒ → なめらかな表⾯への復元がそこそこできている ノイズにもフィットしていってしまい、 形状は破綻してしまっている ⼤きなノイズや外れ値を積極的に無視して復元する能⼒はない デノイズ器としては不⼗分
23. Face template をソースメッシュとして変形[Fig 11] 顎が消えてる ディティールが消えてる 位置はいい感じだけど滑らかさがない *参考 こちらは平均値です
24. Sphere meshをソースとして変形 [Fig 12] 多くのケースで失敗 Sphere startに特化した⼿法と同等の結果だがoursが早い *参考 こちらは平均値です Sphere startに特化した⼿法
25. まとめ 𝐶!(& = 𝐶 𝑥! 1. 対応点を更新する 2. 対応𝐶894から𝐴894, 𝑏894を構成 3. 𝐴894𝑥 = 𝑏894を収束するまで解く 4. 得られた𝑥894を次の形状とする ③で Ax=b を Cholesky や多ステップ PCG でほぼ厳密に解く しかし対応が毎回変わるため，(Aもbも更新) 各ステップで厳密解を求めるのはコストの割にリターンが⼩さい Non-Rigid ICP 課題 提案 : Single-step Transform 各 Non-rigid ICP イテレーションで PCG を 1 ステップだけ実⾏ その 1 歩の質を上げるために 前処理 M を設計 𝑀E=D (𝑘) = 5 𝑀="> (𝑘) + 𝛽′ 5 𝑀DEF (k) ⾼周波ノイズを抑えた 滑らかな更新 対⾓スケーリングで条 件数を下げる Iteration k で両者を使い分けるadaptive設計 Target Source Global alignment 個⼈的メモ:アイデアが⾯⽩いと思った。終盤で1 iterationでresidual落としていけるのすごい（ほぼ対応が更新されなくなるのかな︖） Source テンプレート側はトポロジー固定じゃないとダメ ハイパラは固定解像度でablationにて決めている。テストに使ったものはほぼ10K〜11K vertsで解像度が⼤きく違うもの が少ない。 もしかしたら…😨 Source を⾮剛体(Non-Rigid)に変形して Target 形状に⼀致させる