平成 18 年度大学院数理情報学専攻修士論文審査

滑らかな境界をもつ領域のための等方的メッシュの生成法
56214 長井 超慧
指導教員 杉原 厚吉 教授
平成 18 年 2 月 5 日

1 概要
メッシュは流体力学やコンピュータグラフィックスなど多くの分野で用いられるデータ表現法である．2 次元三角形メッシュおよび 3 次元四面体メッシュに関して，滑らかな境界をもつ領域に対する等方的なメッシュ生成法の研究を行った．
2 次元三角形メッシュ: 等方性の向上を目的とする．Perssonand Strang による，領域境界の曲線部分の表現に優れた等方的メッシュ生成アルゴリズム DistMesh を改良した．
DistMesh で生成したメッシュの頂点密度が低い部分に頂点を追加し，スムージングを施してメッシュを再構成する手法を提案した．
3 次元四面体メッシュ: 等方性を保ちつつ，領域境界を忠実に表現することが目的である．四面体メッシュ表面上の頂点位置の決定に quadric error metric を用いて幾何誤差を減少した．頂点移動で等方性が低下するのを抑えるため，接続する要素の等方性を基に各頂点の移動距離を定めた．

22 次元三角形メッシュの等方性向上法の提案

本提案手法は Persson and Strang [3] による 2 次元三角形メッシュ生成アルゴリズム DistMesh の改良である．DistMeshで生成したメッシュに頂点を追加して，より等方的な 2 次元三角形メッシュを生成する．の隣接頂点同士の間の辺長の総和) − (c の周長) の値が閾値より大きいなら，vi 周辺は密度が低いとする．
以上の判定をメッシュ内部の全頂点に対して行う．頂点密度が低いと判断された頂点がなければ 4. へ．
2. 頂点の追加: 1. で周囲の頂点密度が低いと判断された
頂点 vi の近傍に新たに頂点を 1 つ追加する．vi に接続する辺のうち，lj /lj∗ が最大の辺 ej を特定する．この辺に接続する 2 つの三角形のうち，外接円の半径が大きい方の外心を，頂点追加位置とする．
3. 接続関係と頂点位置の更新: DistMesh に 従 っ て 全 頂点を移動し，Delaunay 三角形分割を行うことを収束するまで繰り返す．1. へ．
4. スムージング: 一様メッシュにはラプラシアンスムー
ジング を，非 一様メッシ ュには optimal Delaunay
triangulation スムージング [2] を施す．

DistMesh の辺長関数は辺長の変化を用いて辺長を定めるため，アルゴリズムの反復のたびに異なる関数値を返す．従ってこのままだと，アルゴリズム 頂点追加 1. の判定基準も反復のたびに変わってしまい，アルゴリズムが収束しないことがある．
これを防ぐため，本手法では新たに絶対的な辺長関数:
h(xei )lj2 /Hj

2.1 DistMesh
生成領域の境界の形状を陰関数で表現することで，境界の曲線部分を滑らかに表現でき，境界上の頂点位置決定も容易になる．頂点同士の接続関係は Delaunay 三角形分割で定める．接続関係を固定したとときの最適な頂点位置は，以下のエネルギー関数を最小化する位置とする．

vi ∈Vej ∈E(vi )kq(lj − lj∗ )(1)

V はメッシュの頂点集合，E(vi ) は頂点 vi に接続する辺集合，k は定数，lj は辺 ej の長さ，lj∗ は辺 ej の理想長を表す．また，
関数 q は以下で定義される関数である:
x (x > 0),q(x) =(2)0 (x ≤ 0).
辺の理想長は通常ユーザ指定の辺長関数で与えられる．しかしDistMesh では，ユーザは辺長変化を定める関数 h を入力し，それを基に独自の辺長関数:

h(xei )lj2 /
h2 (xej )jj(3)を用いて理想長を定める（辺 ei の指定は中点 xei を与えることで行う）．このため相対的な辺長関数ということができる．
頂点の初期位置はメッシュ生成領域内部の格子点である．また，頂点数が初期配置後に増減することはない．
2.2 頂点追加
DistMesh で生成した等方的なメッシュに対し，頂点を追加する以下のアルゴリズムを提案する．
アルゴリズム 頂点追加
入力: DistMesh による等方的な 2 次元三角形メッシュ，このメッシュを生成する際 DistMesh に与えた入力
出力: 等方的な 2 次元三角形メッシュ
1. 頂点密度判定: メッシュ内部の頂点 vi を中心とし，viにおける辺長関数値を半径とする円 c を考える．(vi(4)を定義し，頂点密度判定の際にはこれを用いる．H は入力メッ
シュの j h2 (xej ) の値である．
頂点 vi の次数を ni ，隣接する頂点の座標を vi,j とする（以下本稿中では，頂点 v の座標を v で表記する）．頂点 vi の移動
先はラプラシアンスムージングでは

1  nivi,j(5)j=1niとなる．また，|Ti,j | で vi に接続する要素（2 次元では三角形）
Ti,j の体積（2 次元では面積）を，ci,j で要素 Ti,j の外心の座標を表すと，ODT スムージングによる移動先は以下となる:
 ni
j=1 |Ti,j |ci,j ni.
(6)
j=1 |Ti,j |DistMesh の移動法は頂点を均一に分布させるよう働く．3.
でこの移動法を用い，頂点密度判定に式 4 の絶対的な辺長関数を採用したため，いずれ頂点密度が十分高まり，アルゴリズム頂点追加は収束する．
2.3 結果150 頂点の一様メッシュに対し提案アルゴリズムを施した結果，16 個の頂点が追加された (図 1(a))．非一様メッシュの例
では，687 頂点から 749 頂点になった（図 1(b)）．アルゴリズ
ム適用後の各メッシュと，それぞれとほぼ同数の頂点を持つメッシュとの品質を比較した（表 1）．比較対象のメッシュはDistMesh で生成し，メッシュの種類に応じたスムージングをかけたものである．厳密に同じ頂点数でないのは DistMesh では頂点数の厳密な指定ができないためである．等方性は，各三角形の面積と辺長の二乗和の比で定義されるアスペクト比（値域は [0, 1]，値が大きいほど等方的であることを示す）で測った．辺長関数との整合性を測る指標である辺長比は，辺 ei に対して li∗ /li (li∗ ≤ li ) で，または li /li∗ (li < li∗ ) で定義される．0以上 1 以下の値をとり，1 に近いほど良い．アスペクト比の平表 1. アルゴリズム 頂点追加を実行した結果と比較対象 (詳細は本文参照) のアスペクト比と辺長比の比較．
一様メッシュ
非一様メッシュ
比較対象
提案手法
比較対象
提案手法
頂点数169166751749
アスペクト比 平均
アスペクト比 分散

辺長比 平均

0.970
0.975
0.977
0.981
1.60 ×10−3 1.31 ×10−3 7.81×10−4 5.24×10−40.938
0.945

均値が大きく，分散が小さくなっていることから，頂点追加でより等方的なメッシュが得られた事が分かる．非一様メッシュ
では，辺長比の平均値も大きくなっている．これらの結果から提案アルゴリズムの有効性が確認できた．また，メッシュ内部の頂点の次数を調べた結果，2 次元三角形メッシュで最適な次数 6 のものの割合が増加していた．

3 表面形状を再現する四面体メッシュ生成
本節では，曲面からなる境界をもつ 3 次元領域に対して，等方的な一様四面体メッシュを生成することを目的とする*1 ．
第 2 節とは異なり，本節で提案するアルゴリズムは，等方性は多少低下しても境界との幾何誤差を減少する方針である．
3.1 四面体メッシュ生成
以下に提案するアルゴリズムを示す．頂点の移動や接続関係の変更は全て四面体メッシュに対してのみ行い，表面メッシュは一切変更しない．
アルゴリズム 四面体メッシュ生成
入力: 内部を四面体メッシュ生成領域とする表面メッシュ
Msurf ，四面体メッシュの頂点数を定めるパラメータ
出力: 等方的な四面体メッシュ M
1. 初期四面体メッシュの生成
2. 頂点移動と接続関係の更新: 収束するまで以下を繰り返す．
2-1. 領域内部に対する操作
2-2. 領域境界に対する操作
3. 四面体メッシュ M の抽出
初期四面体メッシュの生成: 入力で与えたパラメータに基づ
く細かさの格子で Msurf 内部を分割し，格子点を M の初期頂点とする．初期頂点に対して Delaunay 単体分割を行う．
領域内部に対する操作: ODT スムージングで生成領域の内部の頂点を動かし，全頂点に対し Delaunay 単体分割を行うことを収束するまで繰り返す．ODT スムージングを用いるのは，境界を除けば ODT スムージングは接続する単体を 2 次元なら正三角形，3 次元なら正四面体に近づけようとするので等方的なメッシュ生成に適するからである（但し正四面体では 3 次元空間を充填できない）．
領域境界に対する操作: まず領域境界を表現する頂点を決定する．次にそれらの頂点を後述の quadric error metric (QEM)を用いて決定した位置に移動し，Delaunay 単体分割を行う．
境界を表現する頂点は，以下のように選ぶ [1]: Msurf の各頂点が M の頂点で最も近い位置にあるものを選ぶとき，Msurfの少なくとも一つの頂点に選ばれた M の頂点を，境界を表現する頂点とする．境界を表現する頂点 vi （座標 (xi , yi , zi )）を選んだ Msurf の頂点集合 VMsurf (vi ) は，vi の Voronoi 領域とMsurf の共通部分の離散的な近似になっている．vi の，表面メッシュ Msurf に対する誤差 (QEM) を∆′ (vi ) = ṽiT

f|f |f · f T /Tf|f | ṽi(7)と定義する．ここで ṽi = (xi , yi , zi , 1) であり，f は Msurfの面 f （面積 |f |）を含む平面 ax + by + cz + d = 0 から定*1 本節のアルゴリズムは第 2 節とは全く独立なものである．

まるベクトル (a, b, c, d)T である．式 7 の和はどちらも，頂点vj ∈ VMsurf (vi ) に接続する面集合に関してとる．アルゴリズ
ム 四面体メッシュ生成 3-1. では，誤差 ∆′ (vi ) が大きい順に次の手順で頂点 vi を移動する．
アルゴリズム 領域境界を表現する頂点 vi の位置決定
1.QEM を最小にする位置 vi∗ の算出
2. 頂点 vi の移動: vi∗ を，vi の移動先候補 vi† の初期値と†する．vi に vi を移動するかどうかを決定する．以下†の基準を共に満たす時のみ vi に vi を移す．
†• vi に vi を移動後，接続する要素の等方性の平均値と最小値の低下が，移動前のものに対して共に一定の割合未満である• vi と vi† との距離が一定以下である
基準を満たさなかった場合には，現在の位置 vi と vi∗を結ぶ線分上で二分探索を行う．この線分の中点を改††めて移動先候補 vi とし判定を行うことを，vi が vi に十分近くなるか移動先が決定するまで繰り返す．
等方性は内接球の半径と外接球の半径の比で測る．
3. 接続関係の更新: Delaunay 単体分割をする
メッシュ抽出: Delaunay 四面体分割は点集合の凸包に対する分割であるため，四面体集合の中からメッシュ生成領域を表現するものを選ぶ必要がある．M の要素として抽出する四面体は，重心が Msurf 内部にあり，かつ，外心から Msurf までの距離の外接球の半径に対する比が十分小さいものとする [1]．
3.2 結果アルゴリズム 四面体メッシュ生成に対する計算機実験の結果が図 2 である．表面メッシュ Msurf （図 2 右半分に，M に重ねて表示）の頂点集合と M の表面上の頂点集合とのハウス
ドルフ距離は 8.40 × 10−2 ，QEM の平均値は 2.17 × 10−6 であった．

(a)

(b)

図 1. 頂点追加後のメッシュの例．大きい四角が新たに追加された頂点，小さい四角が周囲に頂点が少ないと判断された頂点を表す．(a) 一様メッシュ，(b) 非一様メッシュに提案手法を施した結果，それぞれ 16 頂点，62 頂点追加された．

図 2. 左 上: 四 面
体 メ ッ シ ュ M．
左 下: M の 断 面 ．
右: Msurf のワイヤフレームと M を重ねて表示したもの．

参考文献[1] Pierre Alliez, David Cohen-Steiner, Mariette Yvinec,
and Mathieu Desbrun. Variational tetrahedral meshing.
ACM Transactions on Graphics, Vol. 24, No. 3, pp. 617—625, 2005.
[2] Long Chen. Mesh smoothing schemes based on optimal
Delaunay triangulations. In Proceedings of the 13th International Meshing Roundtable, pp. 109—120, 2004.
[3] Per-Olof Persson and Gilbert Strang. A simple meshgeneration in MATLAB. SIAM review, Vol. 46, No. 2,
pp. 329—345, 2004.