• 検索結果がありません。

陰関数曲面の非多様体への拡張,及びその高速な可視化

N/A
N/A
Protected

Academic year: 2021

シェア "陰関数曲面の非多様体への拡張,及びその高速な可視化"

Copied!
6
0
0

読み込み中.... (全文を見る)

全文

(1)グラフィクスとCAD 109−11 (2002. 11. 22). 陰関数曲面の非多様体への拡張,及びその高速な可視化 山崎 俊太郎∗. 加瀬 究†. 池内 克史‡. 概要 ボクセルサンプリングされた関数による陰関数曲面は多様体形状しか扱うことが出来ない. 本報では,関数 の不連続性を考慮することで非多様体形状を陰関数表現する手法を提案する.また,グラフィクスハードウェア を利用して陰関数表現された表面形状を効率的に描画方法を提案する.. Representation and visualization of non-manifold implicit surfaces Shuntaro Yamazaki∗. Kiwamu Kase†. Katsushi Ikeuchi‡. Abstract Implicit surfaces defined by a voxel-sampled implicit field cannot represent non-manifold shapes. We present the method of representing non-manifold surfaces in this form taking account of the discontinuity of both implicit functions and interpolating functions. We also present an effective method of visualizing implicit surfaces making the best of a modern graphics hardware.. 1. はじめに. コンピュータグラフィクスで物体の表面形状を記述する際に用いられる表現形式は,パラメトリック表現と陰 関数表現に分類することが出来る.本報では, 非多様体曲面をパラメトリック表現から陰関数表現へ変換する手 法を提案する. またそのための具体的なアルゴリズムとして, 三角形メッシュを入力としてボリューム化された陰 関数場を出力する手法について説明する. また,このようにして得られた陰関数場に対してボリュームレンダリ ングの手法を適用し,非多様体曲面を直接描画する手法に関しても提案する. この方法はグラフィクスカードの プログラム可能なテクスチャ合成機能を利用して効率的に実装することができ, 対話的な速度で高画質の曲面を 表示できる.. 2. 非多様体曲面の陰関数表現. 2.1 陰関数曲面 3 次元空間内の曲面 S ⊂ R3 に対して, 次を満たす実数値関数 f (p) を考える. p∈S. ⇔. f (p) = 0.. (1). ただし p ∈ R3 は 3 次元空間内の点である. f は S への距離に応じた値を返す関数と考えることができ, S は f の 値が 0 の等値面である. 曲面 S の式 (1) 右辺の形式による定義を陰関数表現と呼び, このとき S を陰関数曲面と 呼ぶ. 実数値関数 f で定義される関数場を計算機で処理する時には, f をサンプリングする必要がある.そこで,あ らかじめ f 格子点でサンプリングしてボリュームデータの形で保持することにより,陰関数場に対する処理を効 率的に行うと都合が良い. 本報ではこのボリュームデータを陰関数ボリュームと呼び, 陰関数場はボリュームで定 義されるものとする. 陰関数場がボリュームで与えられる場合, 格子点の間の点における陰関数の値はサンプル点 での値を補間して決定する. ∗ 東京大学. (University of Tokyo) ,理化学研究所 (RIKEN) (RIKEN) ‡ 東京大学 (University of Tokyo) † 理化学研究所. 1 −59−.

(2) +. -. +. -. +. -. (a-1). +. -. +. -. -. +. +. +. -. -. (a-2). +. -. +. -. (b-1). +. -. *. *. original surface nearest point interpolated surface. (b-2). 図 1: 陰関数表現された非多様体形状. (a-1) (a-2): 従来法. (b-1) (b-2): 提案手法.. 2.2 非多様体曲面 3 次元空間で, 集合 S が 2-多様体であるとは, S 上の任意の点における無限に小さい近傍が 2 次元の円盤と位相 同相であることをいう. 以降では, 多様体という言葉を用いた場合, 2- 多様体を指すものとする. 陰関数場が,サ ンプリングされた値の線形補間で与えられる場合には,空間内の任意の点で関数値が微分可能である.従って等 値面上の任意の点 p の近傍は法線が ∇ f の微小平面近似することができ, このとき陰関数曲面は多様体である. 一方, S が非多様体であるとは, 表面上のある点の近傍領域において, 座標軸が 1 つ以下しかとれないか, もしく は 2 次元局所座標系のとり方が 3 通り以上あることをいう. 本報で扱う形状はパラメトリック曲面, 特に三角形 メッシュで表現できる曲面であり, 個々の三角形は自己交差を考えない限り,境界以外は多様体である. 従って曲 面上で多様体の性質を満たさなくなる点が存在するのは三角形の境界のみであり,これは次の 2 つの場合に分類 できる. 表面に境界がある場合: ある三角形の境界上の点が, 別のどの三角形の境界上にもない場合. 2 次元で模式化され た例を図 1 (a-1) に示す. 表面に分岐がある場合: ある三角形の境界上の点が, 別の 2 つ以上の異なる三角形の境界上にある場合. 2 次元で 模式化された例を図 1(a-2) に示す. ただし入力の三角形の集合内の異なる 2 つの三角形が, その境界以外の点を共有する場合には, 共有している点を 含む直線に沿ってこれらの三角形を分割することにより, 三角形同士の交差は生じないものとする. 図 1(a-1)(a-2) にある通り, 従来手法を用いてこれらの非多様体形状を陰関数曲面に変換すると形状が大きく変 わってしまう. そこで,非多様体曲面を正しく陰関数表現するためには, これらの特徴を陰関数場の中で保持し, また表面として再構成するアルゴリズムを構築する必要がある.. 2.3 距離関数と補間関数の設計 2.1 節で述べたように, 陰関数ボリュームを用いて曲面を陰関数表現するためには, 陰関数場を構成するための 距離関数と, 表面の位置を決定するための補間関数を定義する必要がある. そこでこれらをうまく定めることによ り, 2.2 節で挙げた 2 つの非多様体の特長を陰関数曲面の枠組みの中で表現する方法を考える. まず曲面に境界がある場合に関して, 従来の陰関数表現では境界部分が延長されて, もともと存在しなかった表 面が過剰に生成されてしまう (図 1(a-1)). これは, 従来の陰関数を用いて定義された陰関数ボリューム中にも正し い表面の位置が定義されているが, 補間によって不要な面が生成されていることを意味する. そこで元の曲面が存 在しない位置には表面を生成しないように補間関数を拡張することで, 表面の境界を表現することが可能である (図 ??(b-1)). 距離関数として符号付距離を用いた場合には, 補間関数として線形補間を用いて表面の位置を決定していた. そ こでこれを拡張し, 補間の対象となる 2 つの互いに隣接するボクセルでの符号付距離の差が一定幅より大きい場 合には, このボクセル間には表面を生成しないようにする. 即ち任意の 2 つのサンプル点の間隔を w, 値をそれぞ れ u, v (u ≤ v) とするとき,これら 2 つの点の間に表面が存在するのは u ∈ (−∞, t]. (2). v ∈ [t, ∞)   0 < (−u) − (−t) + (v − t) < αw. (3) (4). が成立する場合のみとする. ただし α (≥ 1) は曲面の陰関数場への変換が誤差を含む場合にも表面を正しく発生さ せるために導入する係数である.陰関数場が格子点サンプリングされた値を線形補間した値で与えられる場合,. 2 −60−.

(3) 空間内の点における関数値は最大格子点間隔の半分程度の誤差を含むと考えられるため,実験では α = 1.5 を使 用した. 上の 3 つの式を満たすとき, 表面が存在する位置は. q=. t−u v−u. (5). で求められる. ただし q (∈ [0, 1]) は表面の位置が, q = 0 の時に値が u の格子点上に, q = 1 の時に値が v の格子点 上になるように正規化されている. 次に表面に分岐がある場合に関して, 従来の陰関数表現では分岐部分において接続している表面の一部が除去 されてしまう (図 1(a-2)). これは, 連続な実数値関数を用いて定義された陰関数ボリューム中では分岐を持つ表面 の位置を正しく表現できないことが原因である. 曲面が多様体の場合, 曲面上の任意の点の近傍において面の表裏を定義することが出来る. したがって表面上の 任意の点は, 最近傍点がその表側にある領域と裏側にある領域の 2 つの領域に挟まれる部分に存在し, このような 領域分割を符号付距離関数を用いて表現することが出来る. ところが表面に分岐があり非多様体になっている場 合には面の裏と表を定義することは出来ないため, 特に分岐線上の点の周辺領域を最近傍点の面の方向から 2 種 類に分類することはできない. 入力の曲面を表裏の判別が出来る曲面パッチに分解しその表裏にそれぞれ番号を振ってこれらを区別すること により, 最近傍点の面の番号を用いて空間を複数領域に分類すると, 分岐のある曲面を陰関数表現することが可能 になる (図 ??(b-2)). 以下にその領域分割のアルゴリズムを示す.. 1. 入力の非多様体曲面を分岐線上に沿って分割し, それぞれ分岐を持たない曲面パッチに分解する. 2. 得られたパッチに順に番号を振る.パッチの表裏を区別し, 表面の番号を i+ , 裏面の番号を i− とする. 3. 空間を格子点 p でサンプリングし,格子点に曲面への Euclid 距離 dE (p) と最近傍点の面の番号 i(p) を割り当 てる.   4. 各格子点 p に対して, 6 隣接点 pn で i(pn ) を調べ, i(p)  i(p n) であるような i(p), i(p n ) の組を列挙する. 5. 上で作られた番号の組を別の新たな番号で置換する.ただし置換の結果,最初に i+ と i− であった番号が同じ 番号になる場合にはその組み合わせに対する置換を行わない. 最終的に番号が 0 から順に並ぶようにする. 6. 置換表に従い各格子点 p での領域番号 i(p) を書き換える. このようにして得られたボリュームの領域番号 i(p) と, 各ボクセルにおける表面への Euclid 距離 d E (p) から, 実 数値の陰関数ボリュームを構成する. そのために, 符号付距離を拡張した, 領域付距離を定義する. 領域付距離にお ける「領域」とは符号付距離における正と負の 2 種類の「符号」を任意の数に拡張したものであり, 距離を表す 実数空間を領域分割することにより実現する. 分割された領域の幅を Ds とする時,領域 i に含まれる距離 d is は   dis ∈ Ds i, D s (i + 1) (6) である. 以上のような領域付距離は以下の手順で計算できる. まず Ds = 2 B と置き,B = Bmax − log2 n から領域の大き さを決定する.ただし n は領域の数である.Bmax は十分に大きい数であり, 計算機上で表現する場合には領域付 距離を表す変数のビット数として与えることができる. このとき各ボクセルの位置 p で, dE (p) と i(p) から領域付 距離 f s (p) を次の式で計算する. f s (p) = min(d E , 2 B − ) + 2 B i(p) (7) ただし (> 0) は f s (p) が式 (6) の半開区間に含まれるように d E (p) を切り捨てるための微小な正の実数である.実 験では  を格子点間隔,Bmax = 8 とした. 領域付距離関数は連続関数ではないため, その値を線形に補間することは出来ない.そこで,与えられた領域 付距離場中の領域数を n,距離場の 2 つのサンプル点の間隔を w, 値をそれぞれ u, v (u ≤ v) とするとき,これら 2 つの点の間に表面が存在するのは   u ∈ 2B i, 2 B (i + 1) , (8)   v ∈ 2 B j, 2 B ( j + 1) , (9). 0 < (u − 2 B i) + (v − 2 B j) < αw, 3 −61−. (10).

(4) が全て成立する i, j (0 ≤ i ≤ j ≤ n − 1) が存在する場合のみとする. α ≥ 1 は式 (5) で用いた係数である. このとき 表面が存在する位置は u − 2Bi q= (11) (u − 2 B i) + (v − 2 B j) で求められる. q の定義は式 (5) と同様である.. 3. 陰関数曲面の可視化. 3.1 スライスを用いたボリュームレンダリング法 陰関数ボリュームに対して ray casting 法 [4] によるボリュー ムレンダリングの手法を適用することで陰関数曲面を直接, 描 画することが可能である. ray casting 法では高画質な描画が可 能である. ray casting は任意の視線上での関数値のサンプリングが必 要であるため比較的コストの高い計算であるが,Lacroute ら はボリュームを x, y, z 軸と垂直な 3 組のスライスの集合で表 現して (図 2 左), ボリュームのサンプリングをスライス平面上 のみで行うことで計算効率を上げて描画を高速化する手法を 提案した [3]. この方法はポリゴンとテクスチャマップを用い て効率的に実装することが可能である [1]. この手法では,ボ リュームのスライスを透明度付き 2D テクスチャとしてポリ 図 2: 左:ボリュームの軸に沿ったスライス生成. 任意の ゴンにマップし, 得られた半透明ポリゴンの集合を視線に対し 方向から表示する場合には x, y, z の 3 種類を持ち, 最も て奥の方から順に α ブレンドしながら表示する.この表示処 視線と垂直なものを選択する. 右:描画面と垂直なスラ 理ではグラフィクスカードの描画加速を利用できるため, 高速 イス生成. 視点変更のたびにスライスを再生成する. に描画できる. またグラフィクスカードを使う場合には,3D テクスチャ機能を利用することで, スライスの方向を固定せずに常に視線方向と垂直なポリゴンを発生させて描 画することで (図 2 右), 多少の速度の減少と引き換えにより高画質な描画を行うことも可能である. 本報告におけ る実験では,視線と垂直なポリゴンを用いる方法を利用した.. 3.2. Pre-integration 法 [2]. スライスを使ったボリュームレンダリングの大きな問題点 は, ボリュームの値の補間が, スライス上では任意の解像度で isosurface 行われるのに対し, スライスと垂直な方向では固定のスライス slice 間隔でしか行われないため, 補間精度の不一致によりモアレと 呼ばれる縞状の不具合や亀裂が観測される点である. この問題 back は, 発生させるスライスの枚数を, 必要な補間精度に応じて変 slab 更することにより回避可能であるが, 描画速度とスライス枚数 front は反比例の関係にあり, 多くの場合モアレが気にならない程度 にスライスを発生させることは現実的ではない. pre-integration この問題に対して Engel らは, 各スライスに, そのスライス (table lookup) 上でのボリュームの値だけでなく, 隣接するスライス上でのボ リュームの値も同時にマップし, スライス間に存在するボリュー ムの影響も考慮することでモアレを除去する, pre-integration 図 3: 隣接するスライス上のテクスチャを利用. 法を提案した. 図 3 左のようにボリュームをスライス集合で表現し, front で示されるスライス S f ront をポリゴンで描画する状 況を考える. 従来のスライス法ではこのポリゴンに対し, スライス上でのボリュームの値を表すテクスチャT f ront をマップして表示を行うため, スライス間のボリュームの値は描画結果に影響しない. 一方, pre-integration 法では T f ront と同時に back で示される隣接スライス S back 上のテクスチャT back もマップし, 2 つのテクスチャを用いて S f ront と S back で囲まれる領域のボリュームの値を補間して表示する. その結果, スライスが存在しない領域も含. }. −62− 4.

(5) めてボリューム全体の影響を描画結果に反映させることが出来る. 実際にある視線での見え方を描画する時には, 視線と S f ront , S back の交点を計算して, front 側のテクセル値 t f ront と back 側のテクセル値 tback から,スライス間 で視線上にあるボリュームがこの位置で観測される色と透明度に与える影響を計算して,ポリゴン上に表示する. 例えば図 3 右にあるように, スライスを用いたボリュームレンダリングの手法で,表面を表示する状況を考え る. 図の中央の視線での見え方を考えるとき,従来のスライスを用いた手法ではスライスが表面上にはないため 表面が描画されない.一方, pre-integration 法ではスライス間のボクセル値の変化を考慮することが出来るため, 少ない枚数のスライスを用いた描画でも正しい結果が得られる. t f ront , tback の組から,観測される色と透明度を計算する処理を前もって行っておき, 結果を pre-integration table と呼ばれる早見表として保存しておくことで計算量を減らすことが可能である. またその表は 2 次元テクスチャ としてグラフィクスカード内に保存しておき,描画時にテクスチャ合成機能を利用して参照することにより高速 に描画処理できる. ボリュームはフォン・シェーディングによって照光処理される.そのため,ボリュームとして与えられるスカ ラー場の勾配をあらかじめ計算しておく必要がある.各ボクセルにおけるスカラー値 s と曲面の勾配 (u x , uy , uz ) が与えられるとき,ここからピクセル値 (s, uy , uz , u x ) を持つ 3D テクスチャを生成し,これを読み込む.描画時に はマルチテクスチャ機能を利用し,3D テクスチャをテクスチャユニット 0 と 1 に読み込み, ユニット 0 から tback が,ユニット 1 から t f ront が取り出せるようにテクスチャ座標を与える. pre-integration table は 2D テクスチャで 表現して t f ront ,tback で与えられる任意のスカラー値の組に対する色を RGB チャネルに,透明度を A チャネルに割 り当てる.ただし照光計算の際にスライス間で勾配を補間するために,B チャネルに勾配を補間する係数を割り 当て,色は実際には RG の 2 チャネルのみで表現する.このテクスチャをユニット 3 に読み込むと,このユニッ トに対するテクスチャ座標がグラフィクスハードウェアによって (t f ront , tback ) に設定され,pre-integration table が 参照される.. 3.3 陰関数曲面の可視化手法 領域付距離場ボリュームで表された陰関数曲面を pre-integration 法で描画するためには, 陰関数ボリュームを 3D テクスチャとして使える形に変換し, t f ront , tback からなる任意の領域付距離の組に対して, それらの間の描画色 を決定するための pre-integration table を作ればよい. まず陰関数ボリュームに関して, 領域付距離は一般に実数で与えられるが, 現行のグラフィクスカードでテクス チャとして使うために距離を 8bit 整数 (0–255) で量子化する. またシェーディングに必要な勾配は,陰関数場の 生成の際に元の曲面の勾配から計算し,同様に量子化しておく.pre-integration table は 256 × 256 のサイズの 2 次 元の RGBA テクスチャとして作成する. RG の 2 つのチャネルに割り当てる色は,領域付距離の距離の組ごとに 任意の色を設定でき, 表面を均一に描画する時には単一色を与える. A チャネルには, スライス間に面が存在しな い場合には A=0(透明), 存在する場合には A=1(不透明) とする. 最後に B チャネルに割り当てる補間係数として式 (11) から計算される, q ∈ [0, 1] を 8bit 整数で量子化した値を用いる. 例として図 ??左図に示す非多様体曲面を陰関数表現して可視化することを考える. このモデルは球面が平面に はめ込まれた形の非多様体形状であり, 空間は 3 つの領域に分類できる. これに対応する pre-integration table は右 図の通りである. 表の横軸が t f ront ∈ [0, 255], 縦軸が tback ∈ [0, 255] に対応し, 表の各点の色が描画される表面の色 を色がない部分は表面が存在せずに透明色で描画される組み合わせに対応する.色の濃さは勾配の補間のための 係数に対応している.. 4. 実験結果. 可視化の実験に利用した計算機は Pentium4 1.7GHz, 主記憶 1.0GB, GeForce4 Ti4600 搭載グラフィクスカード, VRAM 128MB の PC である. 図 4 は非多様体表面モデルを陰関数表現に変換し, 可視化した結果である. 左が入 力として与えたパラメトリック表現された表面モデルで,領域付距離場を用いて 2563 サイズの陰関数ボリュー ムに変換し, 可視化した結果が右である. この表面は境界線と表面の分岐をもつ非多様体三角形メッシュである が, 本手法では境界や分岐が正しく表現できることがわかる. 描画時の色は, 表面が存在する領域の対ごとに, 表面 と裏面を区別して設定することが可能である. この例では全ての組み合わせに異なる色を設定している.. 5 −63−.

(6) 図 4: 描画結果. 上段:入力として与えた三角形メッシュモデル. 下段:陰関数表現に変換した後, 可視化した結果.. 5. まとめと将来課題. 本報では, 陰関数曲面を定義する距離関数と補間関数として不連続関数を用いることにより, 陰関数表現で非多 様体曲面を扱えるように拡張した. またこのように定義された多様体・非多様体陰関数曲面を高速に描画する手 法を提案した. 実験で用いた入力のパラメトリック曲面は三角形メッシュのみであるが, 提案手法は自由曲面にも適用可能で ある.ただし自由曲面を多様体パッチに分割するアルゴリズムなどの細かい変更は必要であると予想される. ま た形状の変形や集合演算を行う際には曲面の陰関数表現が有効であるが, 提案手法で定義された非多様体陰関数 曲面表現に対してもこれらの処理が適用可能であるか考察する必要がある. また本手法ではグラフィクスハードウェアの制限により,陰関数を 8bit 整数で量子化したが,これが原因で形 状を正確に表示することができない.今後,32bit 浮動小数点テクスチャが利用可能なハードウェアが登場する ため,これらを用いた実験を行う予定である. 従来の陰関数曲面では取り扱うことの出来る形状に制限が多く, 表面の境界と表面の分岐という大きな問題に 関して, 提案手法で解決することが可能である. しかしこれ以外にも多くの制約があり, 例えば, 向きつけ不可能な 曲面の取り扱いなどにも拡張することが考えられるが, これらは今後の課題である.. 参考文献 [1] M. Brady, K. Jung, H.T. Nguyen, and T Nguyen. Two-phase perspective ray casting for interactive volume navigation. In Visualization 97, pp. 243–56, 1997. [2] K. Engel, M. Kraus, and T. Ertl. High-quality pre-integrated volume rendering using hardware-accelerated pixel shading. In Eurographics / SIGGRAPH Workshop on Graphics Hardware ’01, pp. 9–16, 2001. [3] Philippe Lacroute and Marc Levoy. Fast volume rendering using a shear-warp factorization of the viewing transformation. In SIGGRAPH ’94, pp. 451–458, 1994. [4] H. Tuy and L. Tuy. Direct 2d display of 3d objects. IEEE Mag. Computer Graphics and Applications, Vol. 4, No. 10, pp. 29–33, 1984.. 6 −64−.

(7)

図 4: 描画結果. 上段:入力として与えた三角形メッシュモデル. 下段:陰関数表現に変換した後, 可視化した結果. 5 まとめと将来課題 本報では, 陰関数曲面を定義する距離関数と補間関数として不連続関数を用いることにより, 陰関数表現で非多 様体曲面を扱えるように拡張した

参照

関連したドキュメント

K orman , Global solution branches and exact multiplicity of solutions for two point boundary value problems, Handbook of Differential Equations, Ordinary Differential Equa- tions,

This shows that by considering some other conditions of local existence and local unicity of the implicit function instead of the conditions from Theorem 1, we can produce

Our goal is to define and examine the “manifold” of all solutions of the system ( ∗ ) using a generalized notion of manifold which, in effect, allows for non-standard solutions..

For example, if we restrict to the class of closed, irreducible 3-manifolds, then as said above, each manifold has a bounded number of incompressible surfaces, but clearly there is

As we have anticipated, Theo- rem 4.1 of [11] ensures that each immersed minimal surface having properly embedded ends with finite total curvature that is in a neighbourhood of M k

In Section 3 the extended Rapcs´ ak system with curvature condition is considered in the n-dimensional generic case, when the eigenvalues of the Jacobi curvature tensor Φ are

Recently, Velin [44, 45], employing the fibering method, proved the existence of multiple positive solutions for a class of (p, q)-gradient elliptic systems including systems

In section 2 we present the model in its original form and establish an equivalent formulation using boundary integrals. This is then used to devise a semi-implicit algorithm