第4回 モデリング

コンピュータグラフィクス論
– モデリング (3) –
2015年4月30日
高山 健志
陰関数表現による形状モデリング
• 3D空間上のスカラー場 𝐹 𝑥, 𝑦, 𝑧 に対し、
𝐹 𝑥, 𝑦, 𝑧 = 0 で定義されるサーフェス
• 別名:ゼロ等値面、陰関数 (ボリューム) 表現
•  サーフェス表現
• 𝐹 𝒑 が正 (負) のとき、𝒑 は物体の外部 (内部)
• 勾配 𝛻𝐹 は法線方向に一致
半径 2 の球
𝐹 𝑥, 𝑦, 𝑧 = 𝑥 2 + 𝑦 2 + 𝑧 2 − 4 = 0
• 𝛻𝐹(𝒑) = 1 ∀𝒑 のとき、F は物体表面からの距離
2
陰関数の例
大半径 R, 小半径 a のトーラス
3
陰関数の例:等電位面 (Metaball)
𝑞𝑖
𝐹𝑖 𝒑 =
𝒑 − 𝒑𝑖
𝐹 𝒑 = 𝐹1 𝒑 + 𝐹2 𝒑 + 𝐹3 𝒑 + 𝐹4 𝒑
𝐹 𝒑 = 𝐹1 𝒑 + 𝐹2 (𝒑)
𝐹 𝒑 = 𝐹1 𝒑 − 𝐹2 (𝒑)
4
陰関数の線形補間によるモーフィング
𝐹 𝒑 = 𝐹1 (𝒑)
1
2
𝐹 𝒑 = 𝐹1 𝒑 + 𝐹2 𝒑
3
3
2
1
𝐹 𝒑 = 𝐹1 𝒑 + 𝐹2 𝒑
3
3
𝐹 𝒑 = 𝐹2 (𝒑)
5
複数の陰関数を組み合わせたモデリング
6
Constructive Solid Geometry (Boolean演算)
和:𝐹𝐴∪𝐵 𝒑 = max 𝐹𝐴 𝒑 , 𝐹𝐵 𝒑
積:𝐹𝐴∩𝐵 𝒑 = min 𝐹𝐴 𝒑 , 𝐹𝐵 𝒑
差:𝐹𝐴−𝐵 𝒑 = min 𝐹𝐴 𝒑 , −𝐹𝐵 𝒑
7
陰関数の表示方法:Marching Cubes
• 等値面を三角形メッシュとして抽出
• 立方体格子の各セルに対し、
(1) 立方体の 8 頂点で関数値を計算
(2) その正負のパターンから、
生成する面のタイプを決定
• 対称性から 15 通りに分類
(3) 関数値の線形補間から
面の位置を決定
• 最も有名 (特許問題でも )
Marching Cubes: A High Resolution 3D Surface Construction Algorithm [Lorensen SIGGRAPH87]
8
Marching Cubes の曖昧性
隣接するセルの間で面が整合しない
不整合を解決する新たなルール
The asymptotic decider: resolving the ambiguity in marching cubes [Nielson VIS91]
9
Marching Tetrahedra
• 立方体の代わりに四面体を使う
• パターンが少なく、曖昧性が無い
 実装が簡単
• 各立方体セルを、6 個の四面体に分割
• (隣接セル間で分割の向きを合わせることに注意)
• きれいな三角形メッシュを取り出す工夫
http://paulbourke.net/geometry/polygonise/
Regularised marching tetrahedra: improved iso-surface extraction [Treece C&G99]
10
シャープなエッジを保持した等値面抽出
格子サイズ:65×65×65
Marching Cubes
改善版
Feature Sensitive Surface Extraction from Volume Data [Kobbelt SIGGRAPH01]
http://www.graphics.rwth-aachen.de/IsoEx/
改善版 (勾配方向も考慮)
Marching Cubes (値のみ考慮)
11
陰関数の表示方法:レイトレーシング
Modelling with implicit surfaces that interpolate [Turk TOG02]
GPU-based Fast Ray Casting for a Large Number of Metaballs [Kanamori EG08]
12
応用例:スケッチベース 3D モデリング
https://www.youtube.com/watch?v=GyJUG2VSvqw
A sketching interface for modeling the internal structures of 3d shapes [Owada SmartGraphics03]
ShapeShop: Sketch-Based Solid Modeling with BlobTrees [Schmidt SBIM05]
13
応用例:メッシュの穴埋め
Filling holes in complex surfaces using volumetric diffusion [Davis 3DPVT02]
14
点群からのサーフェス再構成
15
3D 形状の計測
Range Scanner
(LIDAR)
• 得られるデータ:点群
• 3D座標
• 法線 (面の向き)
Depth Camera
Structured Light
Multi-View Stereo
• 得られない場合もある
• ノイズが多すぎる場合もある
16
点群からのサーフェス形状再構成
• 入力:N 個の点群データ
y
• 座標 𝐱𝑖 = 𝑥𝑖 , 𝑦𝑖 , 𝑧𝑖 と法線 𝐧𝑖 = 𝑛𝑖x , 𝑛𝑖 , 𝑛𝑖z , 𝑖 ∈ 1, … , 𝑁
• 出力:関数 𝑓(𝐱) で、値と勾配の制約を満たすもの
• 𝑓 𝐱𝑖 = 0𝑓𝑖
• 𝛁𝑓 𝐱𝑖 = 𝐧𝑖
• 等値面 𝑓 𝐱 = 0 が出力サーフェス形状
• “Scattered Data Interpolation” と呼ばれる問題
• Moving Least Squares
• Radial Basis Function
CG以外の分野 (e.g. 機械学習) でも重要
17
勾配を制約する二通りの方法
• 法線方向にオフセットした位置に値の制約を追加
• 簡単
• 数学表現そのものに勾配制約を取り入れる (エルミート補間)
• 高品質
値と勾配の制約
エルミート補間
Modelling with implicit surfaces that interpolate [Turk TOG02]
Hermite Radial Basis Functions Implicits [Macedo CGF10]
オフセット法
18
Moving Least Squares による補間
(移動最小二乗)
19
出発点:Least SQuares (最小二乗)
• 求めたい関数が線形だと仮定する:𝑓 𝐱 = 𝑎𝑥 + 𝑏𝑦 + 𝑐𝑧 + 𝑑
• 𝑎, 𝑏, 𝑐, 𝑑 が未知係数
𝐱 ≔ (𝑥, 𝑦, 𝑧)
• データ点における値の制約
・・・
𝐴
𝑓 𝐱𝑁 = 𝑎𝑥𝑁 + 𝑏𝑦𝑁 + 𝑐𝑧𝑁 + 𝑑 = 𝑓𝑁
𝑥𝑁 𝑦𝑁 𝑧𝑁 1
𝑎
𝑏
𝐜𝑐
𝑑
𝑓1
𝑓2
=
𝐟
・・・
𝑥1 𝑦1 𝑧1 1
𝑥2 𝑦2 𝑧2 1
・・・
𝑓 𝐱1 = 𝑎𝑥1 + 𝑏𝑦1 + 𝑐𝑧1 + 𝑑 = 𝑓1
𝑓 𝐱2 = 𝑎𝑥2 + 𝑏𝑦2 + 𝑐𝑧2 + 𝑑 = 𝑓2
𝑓𝑁
• (勾配制約は今は考えない)
20
Overconstrained System
• #未知数 < #制約 (i.e. 縦長の行列)  全ての制約を同時に満たせない
𝐴
𝐴⊺
𝐜 = 𝐟
𝐴
𝐜 =
𝐴⊺
𝐟
𝐴⊺ 𝐴
𝐴⊺
𝐟
• fitting の誤差を最小化:
𝑨𝐜−𝐟
2
=
𝑁
𝑖=1
𝑓 𝐱𝑖 − 𝑓𝑖
2
𝐜 =
−1
21
LSQ の幾何的な解釈
𝑝x 𝑞x
𝑝y 𝑞y
𝑝z 𝑞z
𝑟x
𝛼
= 𝑟y
𝛽
𝑟z
𝐪
𝐫
𝑑
𝛽
𝐩
𝛼
• 𝐩 と 𝐪 が張る空間中で 𝐫 に最も近い点を求める (投影する) ことに相当
• fitting 誤差は投影距離に相当:
𝑑2 = 𝛼𝐩 + 𝛽𝐪 − 𝐫
2
22
Weighted Least Squares (重み付き最小二乗)
• 各データ点ごとの誤差に、重み 𝑤𝑖 をつける
• 重要度、確信度
• 以下の誤差を最小化:
𝑁
𝑖=1
𝑤1
𝑤𝑖 𝑓 𝐱𝑖 − 𝑓𝑖
𝑥1 𝑦1 𝑧1 1
𝑥2 𝑦2 𝑧2 1
𝑤2
𝑤𝑁
𝑥𝑁 𝑦𝑁 𝑧𝑁 1
𝑤1
=
𝑓1
𝑓2
𝑤2
𝑊
𝐟
・・・
𝐴
・・・
𝑊
𝑎
𝑏
𝐜𝑐
𝑑
2
𝑤𝑁
𝑓𝑁
23
Weighted Least Squares (重み付き最小二乗)
𝑊
𝐴
𝐜 =
𝐜 =
𝑊
𝐴⊺ 𝑊 2 𝐴
−1
𝐟
𝐴⊺ 𝑊 2
𝐟
24
Moving Least Squares (移動最小二乗)
• 重み 𝑤𝑖 が、評価位置 𝐱 に依存:
𝑤𝑖 𝐱 = 𝑤( 𝐱 − 𝐱 𝑖 )
• よく使われる関数 (Kernel):
• 𝑤 𝑟 =
2 /𝜎 2
−𝑟
𝑒
• 𝑤 𝑟 =
1
𝑟 2 +𝜖2
評価位置に近いほど
大きな重み
• 重み行列 𝑊 が 𝐱 に依存
 係数 𝑎, 𝑏, 𝑐, 𝑑 が 𝐱 に依存
𝑓 𝐱 = 𝑥 𝑦 𝑧 1
𝑎(𝐱)
𝑏(𝐱)
𝐴⊺ 𝑊 𝐱
𝑐(𝐱)
𝑑(𝐱)
2 𝐴 −1
𝐴⊺ 𝑊 𝐱
2
𝐟
25
法線制約の導入
• 各データ点が表す 1 次式を考える:
𝑔𝑖 𝐱 = 𝑓𝑖 + 𝐱 − 𝐱 𝑖 ⊺ 𝐧𝑖
• 各 𝑔𝑖 を現在位置で評価したときの誤差を最小化:
𝑁
𝑖=1
𝑤𝑖 𝐱 𝑓 𝐱 − 𝑔𝑖 𝐱
𝑥 𝑦 𝑧 1
𝑥 𝑦 𝑧 1
𝑤1 (𝐱)
𝑤2 (𝐱)
𝑥 𝑦 𝑧 1
𝑔1 𝐱
𝑔2 𝐱
𝑤1 (𝐱)
𝑤2 (𝐱)
=
・・・
・・・
𝑤𝑁 (𝐱)
𝑎
𝑏
𝑐
𝑑
2
𝑤𝑁 (𝐱)
Interpolating and Approximating Implicit Surfaces from Polygon Soup [Shen SIGGRAPH04]
𝑔𝑁 𝐱
26
法線制約の導入
法線制約を利用
Input : Polygon Soup
Interpolation
法線方向にオフセットして値を制約
Approximation 1
Approximation 2
Interpolating and Approximating Implicit Surfaces from Polygon Soup [Shen SIGGRAPH04]
Approximation 3
27
Radial Basis Function による補間
(放射基底関数)
28
基本的な考え方
• 関数 𝑓(𝐱) を、基底関数 𝜙(𝐱) の重み付き和として定義:
𝑓 𝐱 =
𝑁
𝑖=1
𝑤𝑖 𝜙(𝐱 − 𝐱𝑖 )
基底関数をデータ位置
𝐱 𝑖 に平行移動
• 放射基底関数 𝜙(𝐱):𝐱 の長さのみに依存
• 𝜙 𝐱 =
• 𝜙 𝐱 =
2 /𝜎 2
−
𝐱
𝑒
1
𝐱 2 +𝑐 2
(Gaussian)
(Inverse Multiquadric)
• 各データ点における制約 𝑓 𝐱𝑖 = 𝑓𝑖 から、重み係数 𝑤𝑖 を求める
29
基本的な考え方
𝜙𝑖,𝑗 = 𝜙 𝐱 𝑖 − 𝐱𝑗 と表記する
𝑓 𝐱1 = 𝑤1 𝜙1,1 + 𝑤2 𝜙1,2 + ⋯ + 𝑤𝑁 𝜙1,𝑁 = 𝑓1
𝑓 𝐱2 = 𝑤1 𝜙2,1 + 𝑤2 𝜙2,2 + ⋯ + 𝑤𝑁 𝜙2,𝑁 = 𝑓2
𝜙1,1 𝜙1,2
𝜙2,1 𝜙2,2
𝑤1
𝑤2
𝜙𝑁,𝑁
𝑓1
𝑓2
=
𝑤𝑁
𝐟
・・・
・・・
・・・
𝐰
Φ
𝜙𝑁,1 𝜙𝑁,2
𝑓 𝐱𝑁 = 𝑤1 𝜙𝑁,1 + 𝑤2 𝜙𝑁,2 + ⋯ + 𝑤𝑁 𝜙𝑁,𝑁 = 𝑓𝑁
𝜙1,𝑁
𝜙2,𝑁
𝑓𝑁
これを解けば良い
30
Gaussian 基底関数を使う場合
𝜙 𝐱 =𝑒
− 𝐱 2 /𝜎 2
• パラメタ 𝜎 の選び方によって、結果が大きく変わる!
𝜎
小
大
• なるべく滑らかな結果を得るには?
Scattered Data Interpolation for Computer Graphics [Anjyo SIGGRAPH14 Course]
31
関数の “曲がり具合” の尺度 (Thin-Plate Energy)
• 2 階微分 (=曲率) の大きさを空間全体で積分したもの:
𝐸2 𝑓 =
𝐱∈ℝ𝑑
Δ𝑓(𝐱) 2 𝑑𝐱
• 1 次元空間の場合:
𝐸2 𝑓 =
𝑓 ′′ 𝑥 2 𝑑𝑥
𝑥∈ℝ
• 2 次元空間の場合:
𝐸2 𝑓 =
𝐱∈ℝ2
𝑓xx 𝐱
2
+ 2𝑓xy 𝐱
2
+ 𝑓yy 𝐱
2
𝑑𝐱
• 3 次元空間の場合:
𝐸2 𝑓 =
𝐱∈ℝ3
𝑓xx 𝐱
2
+ 𝑓yy 𝐱
2
+ 𝑓zz 𝐱
2
+ 2𝑓xy 𝐱
2
+ 2𝑓yz 𝐱
2
+ 2𝑓zx 𝐱
2
𝑑𝐱
32
数学分野の知見
• 制約 𝑓 𝐱𝑖 = 𝑓𝑖 を満たす関数全体のうち、𝐸2 を最小化する関数は
以下の基底を使った RBF として表せる:
• 1 次元空間の場合:𝜙 𝑥 = 𝑥
• 2 次元空間の場合:𝜙 𝐱 = 𝐱
3
2 log
𝐱
• 3 次元空間の場合:𝜙 𝐱 = 𝐱
• 参考
• 有限要素法の場合:離散化した領域上で 𝐸2 を最小化する 𝑓 を近似的に求める
• RBF の場合:グリーン関数を使って 𝐸2 を最小化する 𝑓 を解析的に求める
Scattered Data Interpolation for Computer Graphics [Anjyo SIGGRAPH14 Course]
33
線形項の追加
• 𝐸2 [𝑓] は 2 階微分を使って定義される
 任意の線形項 𝑝 𝐱 = 𝑎𝑥 + 𝑏𝑦 + 𝑐𝑧 + 𝑑 を加えても不変:
𝐸2 𝑓 + 𝑝 = 𝐸2 [𝑓]
• 線形項を未知数に含めることで、関数を一意に定める:
𝑓 𝐱 =
𝑁
𝑖=1
𝑤𝑖 𝜙 𝐱 − 𝐱𝑖 + 𝑎𝑥 + 𝑏𝑦 + 𝑐𝑧 + 𝑑
34
線形項の追加
𝑓 𝐱1 = 𝑤1 𝜙1,1 + 𝑤2 𝜙1,2 + ⋯ + 𝑤𝑁 𝜙1,𝑁 + 𝑎𝑥1 + 𝑏𝑦1 + 𝑐𝑧1 + 𝑑 = 𝑓1
𝑓 𝐱 2 = 𝑤1 𝜙2,1 + 𝑤2 𝜙2,2 + ⋯ + 𝑤𝑁 𝜙2,𝑁 + 𝑎𝑥2 + 𝑏𝑦2 + 𝑐𝑧2 + 𝑑 = 𝑓2
・・・
𝑓 𝐱 𝑁 = 𝑤1 𝜙𝑁,1 + 𝑤2 𝜙𝑁,2 + ⋯ + 𝑤𝑁 𝜙𝑁,𝑁 + 𝑎𝑥𝑁 + 𝑏𝑦𝑁 + 𝑐𝑧𝑁 + 𝑑 = 𝑓𝑁
𝜙𝑁,𝑁 𝑥𝑁 𝑦𝑁 𝑧𝑁 1
𝑓1
𝑓2
𝐰
𝑤𝑁
𝑎
𝑏
𝐜𝑐
𝑑
=
・・・
𝜙𝑁,1 𝜙𝑁,2
P
・・・
Φ
𝜙1,𝑁 𝑥1 𝑦1 𝑧1 1
𝜙2,𝑁 𝑥2 𝑦2 𝑧2 1
・・・
𝜙1,1 𝜙1,2
𝜙2,1 𝜙2,2
𝑤1
𝑤2
𝐟
𝑓𝑁
4 個の未知数 𝑎, 𝑏, 𝑐, 𝑑
が追加されたので、
4 個の制約を追加する
必要がある
35
追加の制約条件:線形関数の再現性
• 「全てのデータ点の制約 𝐱𝑖 , 𝑓𝑖 がある線形関数からのサンプリング
であるとき、RBF による補間結果はその線形関数と一致する」
• これを満たすための条件:
𝑁
𝑖=1 𝑥𝑖 𝑤𝑖
•
𝑁
𝑖=1 𝑦𝑖 𝑤𝑖
=0
•
𝑁
𝑖=1 𝑧𝑖 𝑤𝑖
=0
=0
Φ
P
𝜙𝑁,1 𝜙𝑁,2 𝜙𝑁,𝑁 𝑥𝑁 𝑦𝑁 𝑧𝑁 1
𝑥1 𝑥2
𝑥𝑁 0 0 0 0
𝑦1 𝑦2 ・・・
𝑦𝑁 0 0 0 0
⊺
P
𝑧1 𝑧2
𝑦𝑁 0 0 0 0
1 1
1 0 0 0 0
Scattered Data Interpolation for Computer Graphics [Anjyo SIGGRAPH14 Course]
𝑓1
𝑓2
𝐰
𝐟
𝑤𝑁
𝑎
𝑏
𝐜𝑐
𝑑
=
・・・
•
=0
𝑤1
𝑤2
・・・
𝑁
𝑖=1 𝑤𝑖
𝜙1,𝑁 𝑥1 𝑦1 𝑧1 1
𝜙2,𝑁 𝑥2 𝑦2 𝑧2 1
・・・
•
𝜙1,1 𝜙1,2
𝜙2,1 𝜙2,2
𝑓𝑁
0
0
0
0
36
勾配制約の導入
• 基底関数の勾配 𝛁𝜙 の重み付き和を導入:
𝑓 𝐱 =
𝑁
𝑖=1
• 𝑓 の勾配:
𝛁𝑓 𝐱 =
𝑤𝑖 𝜙 𝐱 − 𝐱 𝑖 + 𝐯𝑖⊺ 𝛁𝜙 𝐱 − 𝐱 𝑖
𝑁
𝑖=1
+ 𝑎𝑥 + 𝑏𝑦 + 𝑐𝑧 + 𝑑
𝑎
𝑤𝑖 𝛁𝜙 𝐱 − 𝐱𝑖 + H𝜙 𝐱 − 𝐱𝑖 𝐯𝑖 + 𝑏
𝑐
• 勾配の制約 𝛁𝑓 𝐱𝑖 = 𝐧𝑖 を追加
Hermite Radial Basis Functions Implicits [Macedo CGF10]
𝜙xx
H𝜙 = 𝜙yx
𝜙zx
𝜙xy
𝜙yy
𝜙zy
Hessian 行列
𝜙xz
𝜙yz
𝜙zz
37
𝑤1
勾配制約の導入
𝐯1
• 1 番目のデータ点について:
𝑤2
値の制約:
𝑓 𝐱1 = 𝑤1 𝜙1,1 + 𝐯1⊺ 𝛁𝜙1,1 + 𝑤2 𝜙1,2 + 𝐯2⊺ 𝛁𝜙1,2 + ⋯ + 𝑤𝑁 𝜙1,𝑁 + 𝐯𝑁⊺ 𝛁𝜙1,𝑁 + 𝐯𝑎𝑥
2 1 + 𝑏𝑦1 + 𝑐𝑧1 + 𝑑 = 𝑓1
・・・
勾配の制約:
1,2
𝛁𝜙1,𝑁
1,1
𝛁𝜙1,2
𝛁𝜙1,1
𝑎
𝛁𝑓 𝐱1 = 𝑤1 𝛁𝜙1,1 + H𝜙1,1 𝐯1 + 𝑤2 𝛁𝜙1,2 + H𝜙1,2 𝐯2 + ⋯ + 𝑤𝑁 𝛁𝜙1,𝑁 + H𝜙1,𝑁 𝐯𝑁 + 𝑏 = 𝐧1
𝑐
𝑤𝑁
⊺
⊺
⊺
𝑥1 𝑦1 𝑧1 1
𝜙1,1 𝛁𝜙1,1
𝜙1,2 𝛁𝜙1,2
𝜙1,𝑁 𝛁𝜙1,𝑁
𝑓1
1 0 0 0
𝐯𝑁
・・・
=
ΦH𝜙
ΦH𝜙
ΦH𝜙
P
1,1
1,2
1,𝑁
1
0 1 0 0
𝐧
Hermite Radial Basis Functions Implicits [Macedo CGF10]
1,𝑁
1
0 0 1 0
𝑎
𝑏
𝑐
𝑑
38
Φ1,1
Φ1,2
Φ1,𝑁
P1
・・・
Φ2,1
Φ2,2
Φ2,𝑁
Φ𝑁,𝑁
・・・
𝑃1⊺
𝑃2⊺
𝑃𝑁⊺
𝐯1
𝐧1
𝑤2
𝑓2
𝐯2
𝐧2
P𝑁
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
=
・・・
Φ𝑁,2
𝑓1
・・・
・・・
・・・
Φ𝑁,1
P2
𝑤1
𝑤𝑁
𝑓𝑁
𝐯𝑁
𝐧𝑁
𝑎
𝑏
𝑐
𝑑
0
0
0
0
39
比較
勾配を制約
オフセットして値を制約
40
参考サーベイ等
• State of the Art in Surface Reconstruction from Point Clouds [Berger
EG14 STAR]
• A survey of methods for moving least squares surfaces [Cheng PBG08]
• Scattered Data Interpolation for Computer Graphics [Anjyo
SIGGRAPH14 Course]
• An as-short-as-possible introduction to the least squares, weighted
least squares and moving least squares for scattered data
approximation and interpolation [Nealen TechRep04]
41
参考ページ
• http://en.wikipedia.org/wiki/Implicit_surface
• http://en.wikipedia.org/wiki/Radial_basis_function
• http://en.wikipedia.org/wiki/Thin_plate_spline
• http://en.wikipedia.org/wiki/Polyharmonic_spline
42