(57)【要約】【課題】  レベルセット法の基本となる界面からの符号付距離の値を、従来例よりも高速に計算する為の技術を提供すること。【解決手段】  カットセル法に基づいて計算する事でＮ次元ベクトルとして得られる頂点の集合と頂点の順序付Ｎ個組によって向き付けられた（Ｎ−１）−単体の集合より構成される（Ｎ−１）−単体的複体Ｇにおいて、与えられた１つの頂点から１個の最近接頂点を検索し（Ｓ１０２）、最近接頂点を含む（Ｎ−１）−単体を検索し（Ｓ１０３）、検索した（Ｎ−１）−単体と与えられた１つの頂点との符号付距離を計算し（Ｓ１０４）、それらの絶対値の最小値を選択する（Ｓ１０５）。【選択図】  図１
2006285614.tif 000002
【特許請求の範囲】
【請求項１】
カットセル法に基づいて計算する事でＮ次元ベクトルとして得られる頂点の集合と頂点の順序付Ｎ個組によって向き付けられた（Ｎ−１）−単体の集合より構成される（Ｎ−１）−単体的複体Ｇにおいて、与えられた１つの頂点から１個の最近接頂点を検索する第１の検索手段と、
当該最近接頂点を含む（Ｎ−１）−単体を検索する第２の検索手段と、
前記第２の検索手段が検索した（Ｎ−１）−単体と、前記与えられた１つの頂点との符号付距離を計算し、それらの絶対値の最小値を選択する計算手段と
を備えることを特徴とする情報処理装置。
【請求項２】
カットセル法に基づいて計算する事でＮ次元ベクトルとして得られる頂点の集合と頂点の順序付Ｎ個組によって向き付けられた（Ｎ−１）−単体の集合より構成される（Ｎ−１）−単体的複体Ｇにおいて、与えられた１つの頂点から１以上の近接頂点を検索する第１の検索手段と、
当該近接頂点を含む（Ｎ−１）−単体を検索する第２の検索手段と、
前記第２の検索手段が検索した（Ｎ−１）−単体と、前記与えられた１つの頂点との符号付距離を計算し、それらの絶対値の最小値を選択する計算手段と
を備えることを特徴とする情報処理装置。
【請求項３】
前記第１の検索手段は、
前記与えられた１つの頂点から最近接頂点を検索する第１の最近接頂点検索手段と、
当該最近接頂点への符号付距離の絶対値から上限値を決定する第１の上限値決定手段と、
距離の絶対値が当該上限値より小さい節点を検索する第１の節点検索手段と
を備えることを特徴とする請求項２に記載の情報処理装置。
【請求項４】
前記第１の検索手段は、（Ｎ−１）−単体的複体が凸集合か否かを示す２値変数に基づいて、前記最近接頂点への符号付距離の絶対値を上限値とし、前記第１の節点検索手段による処理を省略することを特徴とする請求項２又は３に記載の情報処理装置。
【請求項５】
前記第１の検索手段は、
前記与えられた１つの頂点から最近接頂点を検索し、当該最近接頂点を含む（Ｎ−１）−単体を検索し、当該（Ｎ−１）−単体への符号付距離を計算し、当該符号付距離のうち絶対値が最小のものを選択する第２の最近接頂点検索手段と、
当該絶対値最小の符号付距離の絶対値から上限値を決定する第２の上限値決定手段と、
距離の絶対値が当該上限値より小さい節点を検索する節点検索手段と
を備えることを特徴とする請求項２に記載の情報処理装置。
【請求項６】
前記第１の上限値決定手段は、前記上限値を計算する複数の第１の上限値計算手段を備え、
当該第１の上限値計算手段は、
前記最近接節点への距離ｄと格子の分割幅Δ、及び空間の次元Ｎに対して、
Ｎ−１／２ｄ＝ｍΔ＋ｒ
なるｍ及びｒを計算し、Ｎ個の要素全てがｍΔなる値をもつベクトルＴＱと、Ｎ個の要素全てがｒなる値を持つベクトルＴＲを得、ＴＲ１を
ＴＲ１＝（Ｎ３／２ｒ，０，０）    ｉｆ  ｒ≦Ｎ−３／２Δ
＝（（１−Ｎ３／２）Δ＋Ｎ３／２ｒ，Δ，Δ）  ｅｌｓｅ  ｉｆ  （１−Ｎ−３／２）Δ＜ｒ
＝（Δ，Ｎ３／２ｒ−Δ，０）    ｏｔｈｅｒｗｉｓｅ
によって計算し、ベクトルＴＡを
ＴＡ＝ＴＱ＋ＴＲ１
によって計算し、上限値をＴＡのノルム‖ＴＡ‖２として計算することを特徴とする請求項３に記載の情報処理装置。
【請求項７】
前記第２の上限値決定手段は、前記上限値を計算する複数の第２の上限値計算手段を備え、
当該第２の上限値計算手段は、
前記最近接節点への距離ｄと格子の分割幅Δ、及び空間の次元Ｎに対して、
Ｎ−１／２ｄ＝ｍΔ＋ｒ
なるｍ及びｒを計算し、Ｎ個の要素全てがｍΔなる値をもつベクトルＴＱと、Ｎ個の要素全てがｒなる値を持つベクトルＴＲを得、ＴＲ１を
ＴＲ１＝（Ｎ３／２ｒ，０，０）    ｉｆ  ｒ≦Ｎ−３／２Δ
＝（（１−Ｎ３／２）Δ＋Ｎ３／２ｒ，Δ，Δ）  ｅｌｓｅ  ｉｆ  （１−Ｎ−３／２）Δ＜ｒ
＝（Δ，Ｎ３／２ｒ−Δ，０）    ｏｔｈｅｒｗｉｓｅ
によって計算し、ベクトルＴＡを
ＴＡ＝ＴＱ＋ＴＲ１
によって計算し、上限値をＴＡのノルム‖ＴＡ‖２として計算することを特徴とする請求項５に記載の情報処理装置。
【請求項８】
カットセル法に基づいて計算する事でＮ次元ベクトルとして得られる頂点の集合と頂点の順序付Ｎ個組によって向き付けられた（Ｎ−１）−単体の集合より構成される（Ｎ−１）−単体的複体Ｇにおいて、与えられた１つの頂点から１個の最近接頂点を検索する第１の検索工程と、
当該最近接頂点を含む（Ｎ−１）−単体を検索する第２の検索工程と、
前記第２の検索工程で検索した（Ｎ−１）−単体と、前記与えられた１つの頂点との符号付距離を計算し、それらの絶対値の最小値を選択する計算工程と
を備えることを特徴とする情報処理方法。
【請求項９】
カットセル法に基づいて計算する事でＮ次元ベクトルとして得られる頂点の集合と頂点の順序付Ｎ個組によって向き付けられた（Ｎ−１）−単体の集合より構成される（Ｎ−１）−単体的複体Ｇにおいて、与えられた１つの頂点から１以上の近接頂点を検索する第１の検索工程と、
当該近接頂点を含む（Ｎ−１）−単体を検索する第２の検索工程と、
前記第２の検索工程で検索した（Ｎ−１）−単体と、前記与えられた１つの頂点との符号付距離を計算し、それらの絶対値の最小値を選択する計算工程と
を備えることを特徴とする情報処理方法。
【請求項１０】
コンピュータに請求項８又は９に記載の情報処理方法を実行させることを特徴とするプログラム。
【請求項１１】
請求項１０に記載のプログラムを格納したことを特徴とする、コンピュータ読み取り可能な記憶媒体。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、多面体への符号付き距離の計算技術に関するものである。
【背景技術】
【０００２】
自由界面の運動は波動、混相流、気泡や液滴の運動・界面不安定問題、およびMarangoni対流や結晶成長等の工学的問題に関連して、現在幅広い分野で注目されている。こらの問題を数値的に解く手法は大きく３つに分類される：
１。  固定格子法（ＶＯＦ、ＬＳＭ、ＣＩＰ）
２。  移動格子法（ＡＬＥ、ＢＥＭ）
３。  粒子法
最初の２つは格子を用いるものであり、最後の粒子法は格子を用いない方法である。固定格子法は自由界面の運動をEuler的に観察し、その界面を捕獲する方法で、VOF(volume of fraction)法とLSM(level set method)、CIP(cubic interpolation pseudo-particle)法、格子ボルツマン法等がある。移動格子法は自由界面をLagrange的に観察し、その界面を追跡してゆく方法である。ALE(arbitrary Lagrangian-Eulerian)法および境界要素法がある。更に格子を用いない方法としてLagrange的に粒子を追跡する粒子法がある。これらの手法の中で現在最も期待されているのがLSMである。
【０００３】
LSMは領域Ωの界面∂Ωをレベルセット関数φで表現する：
φ（ｘ、ｙ、ｚ、ｔ）＜０  （ｘ、ｙ、ｚ）∈Ω−∂Ω    （１）
φ（ｘ、ｙ、ｚ、ｔ）＝０  （ｘ、ｙ、ｚ）∈∂Ω        （２）
φ（ｘ、ｙ、ｚ、ｔ）＞０  （ｘ、ｙ、ｚ）∈Ω          （３）
レベルセット関数φは一般に、任意の点における界面からの距離を与える符号付距離関数ψに単調非減少関数ξを作用させたものとして与えられる：
φ（ｘ、ｙ、ｚ、ｔ）＝ξ（ψ（ｘ、ｙ、ｚ、ｔ））    （４）
距離関数ψは、界面上での距離を０として、以下のEikonal方程式の解として与えられる：
|∇ψ|＝１      （５）
ここでψ（ｘ、ｙ、ｚ、ｔ）は界面∂Ωからの距離であり、界面∂Ω上でψ（ｘ、ｙ、ｚ、ｔ）＝０である。
【０００４】
界面∂Ωは関数として与えられる場合もあるが、複雑な形状の場合は、非構造メッシュ(三角分割によって得られる)や、カットセル法(適合格子を用いた区分線形近似表現)によって与えられる。
【０００５】
レベルセット関数の最も簡単でよく用いられるものは、ξが恒等写像の場合即ち、
φ（ｘ、ｙ、ｚ、ｔ）＝ψ（ｘ、ｙ、ｚ、ｔ）      （６）
である。
【０００６】
このとき、式（１）は式（５）を外側に向かって、また式（３）は式（５）を内側に向かって時間発展させることによって得られる。
【０００７】
式（１）及び式（３）の最も直接的な解法は、非特許文献１等に開示されているように、界面の時間発展を計算する方法である。この計算方法では、時間発展の回数が概ね一方向の節点の数（Ｎ個）だけ必要であると見積もることが出来るので、Ｎ３の節点の計算をＮ回繰り返すことから、Ｏ（Ｎ４）の計算量が必要である。
【０００８】
非特許文献２等で開示されている “narrow band level set method”では、界面の近傍に対してのみ時間発展を計算してゆく。この近傍幅をｋとすると、Ｏ（ｋＮ３）の計算量でレベルセット関数を計算することができる。しかし、界面の近傍における値には誤差が含まれ、この誤差が時間発展によって累積されることから、精度の点で問題が残る。
【０００９】
非特許文献３、特許文献１、２等で開示されている “fast marching method”は、Ｏ（Ｎ３ｌｏｇＮ）の計算法を提供している。
【００１０】
上記全ての方法が式（５）のEikonal方程式に基づいており、“fast marching method”の計算では、以下の離散近似式を用いて計算する：
｜∂ｕ｜≒（ｍａｘ［Ｄ−ｘｉｊｋｕ、０］２
＋ｍｉｎ［Ｄ＋ｘｉｊｋｕ、０］２
＋ｍａｘ［Ｄ−ｙｉｊｋｕ、０］２
＋ｍｉｎ［Ｄ＋ｙｉｊｋｕ、０］２
＋ｍａｘ［Ｄ−ｚｉｊｋｕ、０］２
＋ｍｉｎ［Ｄ＋ｚｉｊｋｕ、０］２
）１／２＝１      （７）
ここで、Ｄ−ｘｉｊｋおよびＤ＋ｘｉｊｋは風下、および風上差分演算子であり、ｈｘをｘ方向のグリッド幅とすると次式で定義される：
Ｄ−ｘｉｊｋｕ＝
（ｕｉｊｋ−ｕｉ−１ｊｋ）／ｈｘ、Ｄ＋ｘｉｊｋｕ＝（ｕｉ＋１ｊｋ−ｕｉｊｋ）／ｈｘ      （８）
Ｄ−ｙｉｊｋ、Ｄ＋ｙｉｊｋおよびＤ−ｚｉｊｋ、Ｄ＋ｚｉｊｋも同様に定義される。
【００１１】
また以下の離散近似式を用いても良い：
｜∂ｕ｜≒（ｍａｘ［Ｄ−ｘｉｊｋｕ、−Ｄ＋ｘｉｊｋｕ、０］２
＋ｍａｘ［Ｄ−ｙｉｊｋｕ、−Ｄ＋ｙｉｊｋｕ、０］２
＋ｍａｘ［Ｄ−ｚｉｊｋｕ、−Ｄ＋ｚｉｊｋｕ、０］２）１／２＝１      （９）
図１２を用いて非特許文献３、特許文献１、２等で開示されている“fast marching method”の処理について説明する。図１２は、“fast marching method”の処理のフローチャートである。
【００１２】
ステップＳ２００１では、境界上に位置する節点の距離を０に設定し、これらの節点にALIVEとタグを付ける（以下ALIVE節点と呼ぶ）。また、これらの点から１グリッド幅だけ離れた節点にCLOSEというタグを付ける（以下CLOSE節点と呼ぶ。他も同様）。それ以外のすべての節点にFARというタグを付ける。
【００１３】
ステップＳ２００２では、CLOSE節点の距離を式（７）、或いは式（９）を用いて計算する。
【００１４】
ステップＳ２００３では、距離の最小値を持つCLOSE節点を探索し、TRIALとタグ付けする。
【００１５】
ステップＳ２００４では、TRIAL節点の近傍のうち、ALIVE節点以外の節点をCLOSEとタグ付けする。
【００１６】
ステップＳ２００５では、TRIAL節点の近傍節点のうち全てのCLOSE節点に対して、距離を再計算する。
【００１７】
ステップＳ２００６では、TRIAL節点をALIVE節点とする。
【００１８】
ステップＳ２００７では、CLOSE節点およびFAR節点がなくなったかどうかを検査し、そうであれば処理を終了。そうでなければステップＳ２００３へ進む。
【００１９】
以上の反復処理はＯ（Ｎ３）回繰り返される。またステップＳ２００３で実行される距離の最小値をもつCLOSE節点を探索する処理は、高速アルゴリズムを利用することによって、前処理にＯ（ＮｌｏｇＮ）、実際の探索処理にＯ（ｌｏｇＮ）の計算量で出来るので、全体としての計算量はＯ（Ｎ３ｌｏｇＮ）となる。
【非特許文献１】S. Osher, JA. Sethian: “Fronts propagating with curvature dependent speed: Algorithms based on Hamilton-Jacobi formulation,” Journal of Computational Physics, 79, pp.12-49 (1988).
【非特許文献２】D. Adalsteinsson and JA. Sethian: “A fast level set method for propagating interfaces,” Journal of Computational Physics, 188, pp.269-277 (1995).
【非特許文献３】JA. Sethian: “Fast Marching Methods,” SIAM review, vol.41, no.2, pp.199-235 (1999).
【特許文献１】JA. Sethian and AM. Popovici: “Three-dimensional seismic imaging of complex velocity structures,” USP 6,018,499 (Jan. 25. 2000).
【特許文献２】AM. Popovici and JA. Sethian: “Second- and higher-order traveltimes for seismic imaging,” USP (US 6,324,478 B1) (Nov. 27. 2001).
【発明の開示】
【発明が解決しようとする課題】
【００２０】
しかしながら従来技術には以下に述べるような課題があった。
【００２１】
格子の節点における界面からの距離を、その近傍節点の距離の値を用いて、式（７）或いは式（９）で与えられる離散近似式を用いて計算するので、界面からの距離が大きくなるにつれて、離散誤差が累積するという欠点があった。更に従来技術は格子の節点全ての符号付距離を計算する方法であり、任意の１点から界面までの距離を計算する場合でも、任意の１点より界面より近い全ての頂点の符号付距離を計算しなければならないため、計算量が大きい。
【００２２】
本発明は以上の問題に鑑みて成されたものであり、レベルセット法の基本となる界面からの符号付距離の値を、従来例よりも高速に計算する為の技術を提供することを目的とする。
【課題を解決するための手段】
【００２３】
本発明の目的を達成するために、例えば本発明の情報処理装置は以下の構成を備える。
【００２４】
即ち、カットセル法に基づいて計算する事でＮ次元ベクトルとして得られる頂点の集合と頂点の順序付Ｎ個組によって向き付けられた（Ｎ−１）−単体の集合より構成される（Ｎ−１）−単体的複体Ｇにおいて、与えられた１つの頂点から１個の最近接頂点を検索する第１の検索手段と、
当該最近接頂点を含む（Ｎ−１）−単体を検索する第２の検索手段と、
前記第２の検索手段が検索した（Ｎ−１）−単体と、前記与えられた１つの頂点との符号付距離を計算し、それらの絶対値の最小値を選択する計算手段と
を備えることを特徴とする。
【００２５】
本発明の目的を達成するために、例えば本発明の情報処理装置は以下の構成を備える。
【００２６】
即ち、カットセル法に基づいて計算する事でＮ次元ベクトルとして得られる頂点の集合と頂点の順序付Ｎ個組によって向き付けられた（Ｎ−１）−単体の集合より構成される（Ｎ−１）−単体的複体Ｇにおいて、与えられた１つの頂点から１以上の近接頂点を検索する第１の検索手段と、
当該近接頂点を含む（Ｎ−１）−単体を検索する第２の検索手段と、
前記第２の検索手段が検索した（Ｎ−１）−単体と、前記与えられた１つの頂点との符号付距離を計算し、それらの絶対値の最小値を選択する計算手段と
を備えることを特徴とする。
【００２７】
本発明の目的を達成するために、例えば本発明の情報処理方法は以下の構成を備える。
【００２８】
即ち、カットセル法に基づいて計算する事でＮ次元ベクトルとして得られる頂点の集合と頂点の順序付Ｎ個組によって向き付けられた（Ｎ−１）−単体の集合より構成される（Ｎ−１）−単体的複体Ｇにおいて、与えられた１つの頂点から１個の最近接頂点を検索する第１の検索工程と、
当該最近接頂点を含む（Ｎ−１）−単体を検索する第２の検索工程と、
前記第２の検索工程で検索した（Ｎ−１）−単体と、前記与えられた１つの頂点との符号付距離を計算し、それらの絶対値の最小値を選択する計算工程と
を備えることを特徴とする。
【００２９】
本発明の目的を達成するために、例えば本発明の情報処理方法は以下の構成を備える。
【００３０】
即ち、カットセル法に基づいて計算する事でＮ次元ベクトルとして得られる頂点の集合と頂点の順序付Ｎ個組によって向き付けられた（Ｎ−１）−単体の集合より構成される（Ｎ−１）−単体的複体Ｇにおいて、与えられた１つの頂点から１以上の近接頂点を検索する第１の検索工程と、
当該近接頂点を含む（Ｎ−１）−単体を検索する第２の検索工程と、
前記第２の検索工程で検索した（Ｎ−１）−単体と、前記与えられた１つの頂点との符号付距離を計算し、それらの絶対値の最小値を選択する計算工程と
を備えることを特徴とする。
【発明の効果】
【００３１】
本発明の構成により、レベルセット法の基本となる境界からの符号付距離の値を、従来例よりも高速に計算することができる。
【発明を実施するための最良の形態】
【００３２】
以下添付図面を参照して、本発明を好適な実施形態に従って詳細に説明する。
【００３３】
［第１の実施形態］
図２は、本実施形態に係る情報処理装置として機能するコンピュータの基本構成を示すブロック図である。なお、このようなコンピュータとしては、一般のＰＣやＷＳなどが適用可能である。
【００３４】
同図において２０１はＣＰＵで、ＲＡＭ２０４に記憶されているプログラムやデータを用いて本コンピュータ全体の制御を行うと共に、後述する各処理を実行する。
【００３５】
２０２は表示装置で、ＣＲＴや液晶画面等により構成されており、ＣＰＵ２０１による処理結果を画像や文字などでもって表示することができる。
【００３６】
２０３は入力装置で、キーボードやマウスなどのデバイスにより構成されており、各種の指示をＣＰＵ２０１に対して入力することができる。
【００３７】
２０４はＲＡＭで、外部記憶装置２０５からロードされたプログラムやデータを一時的に記憶する為のエリアを備えると共に、ＣＰＵ２０１が各種の処理を実行する際に使用するワークエリアを備える。
【００３８】
２０５は外部記憶装置で、ハードディスクドライブ装置などの大容量情報記憶装置により構成されており、ここにＯＳ（オペレーティングシステム）やＣＰＵ２０１に後述する各処理を実行させるためのプログラムやデータが保存されており、これらの一部もしくは全部はＣＰＵ２０１の制御に従ってＲＡＭ２０４にロードされ、ＣＰＵ２０１による処理対象となる。
【００３９】
２０６は通信装置で、外部装置からデータを受信したり、逆に送信したりする為のインターフェース装置として機能するものである。
【００４０】
まず、陰関数表現された曲面をカットセル法によって三角分割する方法について説明する。図３は、適合分割格子によるカットセル法について説明する図である。同図において３０１は適合分割格子、実線３０２は界面、破線３０３は区分線形近似された界面（３次元の場合は三角分割）、３０４は再分割が必要な領域である。同図から３次元の場合を想像して頂きたい。
【００４１】
図４はカットセル法による三角分割の処理のフローチャートである。なお、同図のフローチャートに従った処理をＣＰＵ２０１に実行させるためのプログラムやデータはＲＡＭ２０４に記憶されており、これをＣＰＵ２０１が用いて処理を行うことで、本コンピュータは以下説明する各処理を実行することになる。
【００４２】
なお、本フローチャートに従った処理を行う前段で、先ず、計算対象の各格子点（各節点）の位置データが予め外部記憶装置２０５に保存されているので、これをＲＡＭ２０４にロードしておく。また、未処理の格子の集合を以下「Λ」とする。
【００４３】
先ず、集合Λに含まれている格子群のうち、１つを以下の処理対象とすべく、この１つの格子を構成する４つの格子点の位置データをＲＡＭ２０４の所定のエリア（例えば計算対象のデータを保持するエリア）に読み出す（ステップＳ４０１）。
【００４４】
次に、この格子の各稜線について、界面との交点を求める（ステップＳ４０２）。ステップＳ４０２における処理、即ち、交点を求める処理については図５を用いて後述する。
【００４５】
次に、交点が複数ある稜線が存在するかどうかをチェックし、存在する場合には処理をステップＳ４０４に進め、格子を分割する（ステップＳ４０４）。通常の適合格子では各軸ごとに２分割する。従って３次元空間における格子の場合は８分割されることになる。次に、分割された各格子について、稜線と界面の交点を求める（ステップＳ４０５）。
【００４６】
この処理は、分割によって生成されたそれぞれの格子（３次元の場合は８個）について図４のフローチャートに従った処理（即ちステップＳ４０１〜Ｓ４０６）を実行することにより成されるものである。
【００４７】
そして次に、まだ未処理の格子があるか否か、即ち集合Λが空集合であるのか否かをチェックし（ステップＳ４０６）、空集合ではない場合には処理をステップＳ４０１に戻し、未処理の格子について以下の処理を行う。
【００４８】
次に、上記ステップＳ４０２における処理、即ち、格子の各稜線について界面との交点を求める処理について説明する。
【００４９】
図５は、ステップＳ４０２における処理を説明するための図である。点５１６の座標値を(Xi、Yi、Zi)とする。また、格子間隔５０１、５０２、５０３をそれぞれ、dXi、dYi、dZiとする。
【００５０】
界面は陰関数表現、即ち
ｆ（ｘ、ｙ、ｚ）＝０      （１０）
という関数で表現されているとする。稜線５０４と界面との交点は、
ｆ（ｘ、Ｙｉ、Ｚｉ）＝０      （１１）
を満たす実数ｘのうち、［Ｘｉ、Ｘｉ＋ｄＸｉ］を満たすものをＸとすれば、（Ｘ、Ｙｉ、Ｚｉ）で与えられる。一般にＸは１個とは限らない。同様に稜線５０５〜５１１まで同じ処理を繰り返す。
【００５１】
１個の稜線と界面との交点が２個以上あった場合は格子を分割して、上記処理を実行する。従って、最終的には、１個の稜線と界面との交点は高々１個になっている。
【００５２】
以上の処理によって得られた交点を結ぶと多角形が得られる。得られる多角形については、図６に示すように６通りに分類することが出来るが、これらの交点より成る多角形は必ずしも三角形ではない。そこで図７に示すように三角形を超える多角形は三角形に分割する。以上のようにして、式（１０）によって陰関数表現された界面の三角形分割が得られる。
【００５３】
本実施形態に係る符号付距離計算方法は、上記処理によって得られる三角形分割で表現されている界面を対象としている。
【００５４】
次に、与えられた点Ａから界面までの符号付距離の計算方法について図１を用いて説明する。図１は、与えられた点Ａから界面までの符号付距離を計算する処理のフローチャートである。なお、同図のフローチャートに従った処理をＣＰＵ２０１に実行させるためのプログラムやデータはＲＡＭ２０４に記憶されており、これをＣＰＵ２０１が用いて処理を行うことで、本コンピュータは以下説明する各処理を実行することになる。
【００５５】
なお、本フローチャートに従った処理を行う前段で、先ず、計算対象の各格子点（各節点）の位置データが予め外部記憶装置２０５に保存されているので、これをＲＡＭ２０４にロードしておく。
【００５６】
また以下の説明では、対象とする空間の次元をＮ、界面は（Ｎ−１）−単体のみによって構成される（Ｎ−１）−単体的複体とする。即ち、Ｎが２の場合は界面は２次元空間における曲線であり（Ｎ−１）−単体的複体は折線、Ｎが３の場合は界面は３次元空間における曲面であり（Ｎ−１）−単体的複体は三角分割である。
【００５７】
先ず、界面形状のデータを外部記憶装置２０５からＲＡＭ２０４にロードし、これを陰関数の組み合わせ、或いは通常の３次元ＣＡＤデータに変換し、カットセル法によって表現する（ステップＳ１０１）。なお、予め変換済みのものを用いるようにしても良い。この時点で、各頂点の位置データが与えられることになる。
【００５８】
次に、１つの頂点（点Ａ）に最も近い頂点を上記位置データを参照して検索する（ステップＳ１０２）。検索した頂点を要素とする集合をΛｖとする。なお、三次元空間での点Ａの最近接点を探索する方法としては、例えば特開２０００−１１２９４３に提示されている方法が有効である。
【００５９】
図１に戻って次に、上記集合Λｖに含まれる頂点を含む全ての（Ｎ−１）−単体を探索する（ステップＳ１０３）。検索した（Ｎ−１）−単体の集合をΛＴとする。
【００６０】
次に、この集合ΛＴに含まれる（Ｎ−１）−単体への符号付距離ηを計算する（ステップＳ１０４）。
【００６１】
次に、この集合ΛＴに含まれる（Ｎ−１）−単体への符号付距離のうち、絶対値が最も小さいものを、点Ａから界面への符号付距離（最短距離）ηｍｉｎとして選択する（ステップＳ１０５）。
【００６２】
次に、ステップＳ１０４における処理の詳細について、図８から１１を用いて説明する。
【００６３】
まず、Ｎが２以上の（Ｎ−１）−単体は、その構成要素の並びによって向き付けられていることを前提とする。即ち、２次元空間に存在する点Ｖ１と点Ｖ２によって構成される線分（１−単体）は、［Ｖ１Ｖ２］と記述することによって、点Ｖ１から点Ｖ２へ向かうとき、右を外側、左を内側と定義する。図８において１２０１を点Ｖ１、１２０２を点Ｖ２とすると、上側が内側で下側が外側ということになる。線分［Ｖ１Ｖ２］が３次元空間に存在するときは、点Ｖ１及び点Ｖ２と共線関係にない点Ｖ３を与えることによってこれら３点を含む平面を構成し、その平面上で上記定義を採用すればよい。また３次元空間に存在する点Ｖ１、点Ｖ２と点Ｖ３によって構成される三角形（２−単体）は、［Ｖ１Ｖ２Ｖ３］と記述することによって、点Ｖ１→点Ｖ２→点Ｖ３→点Ｖ１の回転に対して右ねじの法則に従って、内側と外側を定義する。図１０において、１４０２を点Ｖ１、１４０３を点Ｖ２、１４０４を点Ｖ３とすると、１４０５が［Ｖ１Ｖ２Ｖ３］の単位法線ベクトルであり、その向かう方向が外側ということになる。
【００６４】
まず、図８及び９を用いて２次元空間における線分（１−単体）の符号付距離ηの計算方法について説明する。
【００６５】
ここで、内外判定の基準となる線分は図８の点１２０１（以下Ｖ１）と点１２０２（以下Ｖ２）より［Ｖ１Ｖ２］と記述されるものとし、内外判定の対象となる点は図８の１２０５（以下Ａ）、１２０７（以下Ｂ）、及び１２０９（以下Ｃ）とする。
【００６６】
先ず、点Ｖ１を原点、点Ｖ１から点Ｖ２へ向かうベクトルを単位ベクトルｅとして、座標軸Ｘを定義する。座標系Ｘでは‖ｅ‖２は１であるが、元の座標系では‖ｅ‖２は１でないことに注意されたい。
【００６７】
このような前提の元で、先ず、線分［Ｖ１Ｖ２］の単位法線ベクトルｎを以下の式（１２）に従って求める。
【００６８】
ｎ＝（ｅｙ、−ｅｘ）Ｔ／‖ｅ‖２      （１２）
次に、原点Ｖ１から与えられた点（Ａ、Ｂ、或いはＣ）へ向かうベクトル（それぞれａ、ｂ、ｃ）のＸ成分を計算する（ステップＳ１３０２）。例えばベクトルａのＸ成分は以下の式（１３）で計算される。
【００６９】
ａＸ＝ｅＴａ／‖ｅ‖２      （１３）
ベクトルａ、ｂ、ｃのＸ成分はそれぞれ、図８において１２０６、１２０８、１２１０に相当する。
【００７０】
次に、求めたＸ成分の値をチェックし、これが０以下の場合は処理をステップＳ１３０４に進め、０を超え１未満の場合は処理をステップＳ１３０５に進め、１以上の場合は処理をステップＳ１３０６に進める。
【００７１】
求めたＸ成分の値が０以下の場合、与えられた点Ｂから点Ｖ１への距離ｄを以下の式に従って求める（ステップＳ１３０４）。
【００７２】
ｄ＝‖ｂ‖２      （１４）
そして処理をステップＳ１３０７に進める。
【００７３】
一方、求めたＸ成分の値が０を超え１未満の場合、点から線分［Ｖ１Ｖ２］へ引いた垂線（線分［Ｖ１Ｖ２］の単位法線１２０４に平行）の長さｄを以下の式に従って求める（ステップＳ１３０５）。
【００７４】
ｄ＝‖ａ−ａＸｅ‖２      （１５）
この距離ｄは単位法線ｎより次式で計算することもできる。
【００７５】
ｄ＝‖ａＴｎ‖２      （１６）
そして処理をステップＳ１３０７に進める。
【００７６】
一方、求めたＸ成分の値が１以上の場合、与えられた点Ｃから点Ｖ２への距離ｄを以下の式に従って求める（ステップＳ１３０６）。
【００７７】
ｄ＝‖ｃ−ｅ‖２      （１７）
そして処理をステップＳ１３０７に進め、与えられた点Ａに対する距離の符号ｓｉｇｎ（Ａ）を以下の式に従って求める（ステップＳ１３０７）。
【００７８】
ｓｉｇｎ（Ａ）＝＋１    ｉｆ  ａＴｎ＞０
−１    ｏｔｈｅｒｗｉｓｅ                    （１８）
図８の場合には、点Ａと点Ｂの距離の符号は正、点Ｃの符号は負になる。
【００７９】
また、ステップＳ１３０７では、符号付距離ηを以下の式に従って求める。
【００８０】
η＝ｓｉｇｎ（Ａ）×ｄ      （１９）
次に、図１０及び１１を用いて３次元空間における三角形（２−単体）の符号付距離の計算方法について説明する。ここで、内外判定の基準となる三角形は図１０の点１４０２（以下Ｖ１）、点１４０３（以下Ｖ２）及び点１４０４（以下Ｖ３）より［Ｖ１Ｖ２Ｖ３］と記述されるものとし、内外判定の対象となる点は図１０の１４０６（以下Ａ）と１４０８（以下Ｂ）とする。また１４０１は三角形［Ｖ１Ｖ２Ｖ３］を含む平面Ωである。尚、三角形［Ｖ１Ｖ２Ｖ３］の単位法線ベクトル１４０５(以下ｎと書く)、及び平面Ωの方程式ｘＴｎ＝ｕは予め計算されているものとする。ここで、ベクトルｘに対して関数ｆ（ｘ）を次式で定義する。
【００８１】
ｆ（ｘ）＝ｘＴｎ−ｕ      （２０）
先ず、与えられた点Ａの位置ベクトルａに対し、次式を用いて点Ａの平面Ωへの正射影１４０７（以下Ａ’）の位置ベクトルａ’を求める（ステップＳ１５０１）。
【００８２】
ａ’＝ａ−ｆ（ａ）×ｎ      （２１）
また、ステップＳ１５１では距離の符号ｓｉｇｎ（Ａ）を以下の式に従って求める。
【００８３】
ｓｉｇｎ（Ａ）＝＋１    ｉｆ  ｆ（ａ）＞０
−１    ｏｔｈｅｒｗｉｓｅ      （２２）
次に、点Ａ’が三角形［Ｖ１Ｖ２Ｖ３］の内部（辺上を含む）にあるか、外部にあるかを以下の式を用いてチェックする（ステップＳ１５０２）。
【００８４】
例えば点Ａ’から三角形の各頂点へ向かう単位ベクトルをそれぞれｖ１、ｖ２、及びｖ３とすると、点Ａ’が三角形の内部にあるときは次式を満足する。
【００８５】
Ｃｏｓ−１（ｖ１、ｖ２）＋Ｃｏｓ−１（ｖ２、ｖ３）＋Ｃｏｓ−１（ｖ３、ｖ１）＝２π      （２３）
この式（２３）を満足すれば処理をステップＳ１５０３に進め、満足していなければ処理をステップＳ１５０４に進める。
【００８６】
式（２３）を満足している場合には、以下の式に従って符号付距離を計算し（ステップＳ１５０３）、処理を終了する。
【００８７】
ηｍｉｎ＝ｆ（ａ）      （２４）
一方、式（２３）を満足していない場合には先ず、三角形［Ｖ１Ｖ２Ｖ３］を構成する各辺、［Ｖ１Ｖ２］、［Ｖ２Ｖ３］、及び［Ｖ３Ｖ１］と点Ａとの距離を計算し、各辺、［Ｖ１Ｖ２］、［Ｖ２Ｖ３］、及び［Ｖ３Ｖ１］と点Ａとの距離の絶対値を、それぞれ、ｄ１、ｄ２、及びｄ３とすると界面から点Ａへの符号付距離ηｍｉｎを以下の式に従って求める（ステップＳ１５０４）。
【００８８】
ηｍｉｎ＝ｓｉｇｎ（Ａ）×ｍｉｎ［ｄ１、ｄ２、ｄ３］    （２５）
ここで、以下の３種類の計算方式の計算時間のオーダーを見積もる。
【００８９】
方式１  全ての（Ｎ−１）−単体に対する符号付距離を計算し、そのうち絶対値が最小のものを選択する。
【００９０】
方式２  ｆａｓｔ  ｍａｒｃｈｉｎｇ法。
【００９１】
方式３  上記本実施形態に係る方法。
【００９２】
任意の点から（Ｎ−１）−単体への距離の計算時間のオーダーをＫとすると、方式１はＯ（ＮＶＫ＋ｌｏｇＮ）である。
【００９３】
方式２は差分演算の計算時間のオーダーをＳとすると、Ｏ（ＳＮＶ３ｌｏｇＮＶ）である。
【００９４】
方式３は界面形状の複雑さに依存する。界面形状が凸集合の場合は、最短頂点探索のための計算時間Ｏ（ＮＶ）（或いはＯ（ＮＶｌｏｇＮＶ）の前処理をすればＯ（ｌｏｇＮＶ）で探索できる）と最短頂点を含む（Ｎ−１）−単体の符号付距離計算時間Ｏ（ＮＴＫ）の和Ｏ（ＮＴＫ＋ＮＶ）である。局所凸性が崩れるに従って処理すべき頂点数が増加し、それに伴って処理すべき（Ｎ−１）−単体の数も増加する。
【００９５】
計算精度の点から比較すると、方式１と方式３はカットセル法の近似精度と同等の精度が得られる。一方、方式２は微分を差分で近似し、それをもとに距離を和分によって算出するので、界面から遠ざかるにつれて精度は低下する。
【００９６】
以上説明したように、本実施形態によれば、レベルセット関数の基本となる界面からの距離関数を、従来例よりも高精度で高速に計算することが可能になった。
【００９７】
［第２の実施形態］
本実施形態では、第１の実施形態とは異なる方法でもって、与えられた点Ａから界面までの符号付距離の計算を行う。なお、方法は異なるが、本実施形態に係る符号付距離計算方法は、上記処理によって得られる三角形分割で表現されている界面を対象としている。
【００９８】
図１３は、与えられた点Ａから界面までの符号付距離を計算する処理のフローチャートである。なお、同図のフローチャートに従った処理をＣＰＵ２０１に実行させるためのプログラムやデータはＲＡＭ２０４に記憶されており、これをＣＰＵ２０１が用いて処理を行うことで、本コンピュータは以下説明する各処理を実行することになる。
【００９９】
以下の説明では、対象とする空間の次元をＮ、界面は（Ｎ−１）−単体のみによって構成される（Ｎ−１）−単体的複体とする。即ち、Ｎが２の場合は界面は２次元空間における曲線であり（Ｎ−１）−単体的複体は折線、Ｎが３の場合は界面は３次元空間における曲面であり（Ｎ−１）−単体的複体は三角分割である。
【０１００】
先ず、点Ａに最も近い頂点を探索する（ステップＳ６３００）。検索したこの頂点を要素とする集合をΛｖとする。
【０１０１】
次に、対象とする界面が凸集合（本実施形態においては、２次元の場合は凸多角形、３次元の場合は凸多面体）かどうかを示すフラグが１かどうかをチェックし（ステップＳ６３０１）、このフラグ値が１であれば処理をステップＳ１０５に進め、１でなければ処理をステップＳ６３０２に進める。
【０１０２】
このフラグ値が２でなければ、際近接頂点までの距離ｄｍｉｎを元に、近接頂点探索のための上限値ｄｓｕｐを計算する（ステップＳ６３０２）。この計算方法は図１４を用いて後述する。
【０１０３】
次に、距離が上限ｄｓｕｐ以下の頂点を探索する（ステップＳ６３０３）。検索した頂点を上記集合Λｖに追加する。
【０１０４】
次に、この集合Λｖに含まれる頂点を含む全ての（Ｎ−１）−単体を探索する（ステップＳ６３０４）。その集合をΛＴとする。
【０１０５】
次に、この集合ΛＴに含まれる（Ｎ−１）−単体への符号付距離ηを計算する（ステップＳ６３０５）。（Ｎ−１）−単体への符号付距離計算法については後述する。
【０１０６】
次に、この集合ΛＴに含まれる（Ｎ−１）−単体への符号付距離のうち、絶対値が最も小さいものを、点Ａから界面への符号付距離（最短距離）ηｍｉｎとして選択する（ステップＳ６３０６）。
【０１０７】
次に、上記ステップＳ６３０２における処理の詳細について図１４を用いて説明する。図１４は、２次元空間における、近接頂点探索のための上限値ｄｓｕｐを計算する処理を説明する図である。
【０１０８】
同図において８０１は適応格子である。８０２は与えられた点Ａである。界面への符号付距離ηｍｉｎを知りたい点は格子空間における節点であるので、点Ａは節点上にあることが保証される。
【０１０９】
今、点Ａ８０２から最短距離にある頂点（最近接頂点）が８０３であったとする。図中、確かに他の頂点は８０３より遠い位置にある。これは、点Ａ８０２から頂点８０３までの距離ｄｍｉｎを半径とする円８０４より明らかである。
【０１１０】
しかし、頂点８０５と８０６より構成される１−単体（線分）上の点８０７は、最近接頂点８０３より近くなることが可能である。即ち、頂点Ａ、Ｂの距離をｄ（Ａ，Ｂ）と書くと、
ｄｍｉｎ＝ｄ（８０２，８０３）＜ｄ（８０２，８０５）
かつｄｍｉｎ＝ｄ（８０２，８０３）＜ｄ（８０２，８０６）  （２６）
のとき
ｄｍｉｎ＝ｄ（８０２，８０３）＞ｄ（８０２，８０７）  （２７）
なる点８０７が存在し得る。このことを以下で厳密に説明する。
【０１１１】
与えられた点Ａはグリッドの節点上にあり、カットセル法によって得られた頂点は稜線上にあるので、節点Ａから頂点Ｖへ向かうベクトルＴＡ＝（ＴＡｘ，ＴＡｙ，ＴＡｚ）は以下のように記述することができる。
【０１１２】
ＴＡｘ＝ΣｖＡｘｊ２−ｊ＋ＴＡｘＲ
ＴＡｙ＝ΣｖＡｙｊ２−ｊ＋ＴＡｙＲ
ＴＡｚ＝ΣｖＡｚｊ２−ｊ＋ＴＡｚＲ      （２８）
ここでΣはｊ＝０〜Ｌ−１について加算することを意味する。
【０１１３】
また、ＴＱ＝（ＴＱｘ，ＴＱｙ，ＴＱｚ）を
ＴＱｘ＝ΣｖＱｘｊ２−ｊ
ＴＱｙ＝ΣｖＱｙｊ２−ｊ
ＴＱｚ＝ΣｖＱｚｊ２−ｊ    （２９）
とすると、ＴＱは節点Ａから節点ＶＱへ向かうベクトルであり、節点ＶＱは、頂点Ｖ及びベクトルＴＡに最も近い節点である。
【０１１４】
ここでＴＲ＝（ＴＲｘ，ＴＲｙ，ＴＲｚ）とすると、ＴＲは節点ＶＱから頂点Ｖへ向かうベクトルであり、TＡ＝ＴＱ＋ＴＲと書ける。
【０１１５】
さて点Ａの最近接頂点として、頂点Ｖが選ばれたとすると、その距離ｄｍｉｎ＝ｄ（Ａ，Ｖ）は次式で与えられる。
【０１１６】
ｄｍｉｎ＝ｄ（Ａ，Ｖ）＝‖ＴＡ‖２
＝｛（ＴＱｘ＋ＴＲｘ）２＋（ＴＱｙ＋ＴＲｙ）２＋（ＴＱｚ＋ＴＲｚ）２｝１／２    （３０）
ただし、‖Ｘ‖２はベクトルＸの２−ノルムである。
【０１１７】
一方、界面を構成する(Ｎ−1)−単体Ｔに対して、任意の点ＡからＴを構成する頂点ＶＴｊ，（ｊ＝１，２，・・・）へ向かうベクトルをＴＡｊ、ＡからＴへの最短距離に対応するベクトルをＴＡとし、ＴＡｊのノルムの最小値とＴＡのノルムとの比、即ち
ｍｉｎ［‖ＴＡ１‖２，‖ＴＡ２‖２，…］／‖ＴＡ‖２    （３１）
が最大になる場合は（図１４には２次元の場合を示しているが、３次元の場合も同様）
‖ＴＡ１‖２＝‖ＴＡ２‖２＝・・・＝ＴＡＭＩＮ      （３２）
のときである。ただしＴＡＭＩＮは実数値である。
【０１１８】
ここで上記式２９と同様に点Ｑ及びベクトルＴＱを考えると、ＴＡｊ＝ＴＱ＋ＴＲｊ、及びＴＡ＝ＴＱ＋ＴＲと書ける。これらの式によれば、上記式３１は、
ｍｉｎ［‖ＴＲ１‖２，‖ＴＲ２‖２，…］／‖ＴＲ‖２    （３３）
と書くことが出来、上式が最大になる条件は
‖ＴＲ１‖２＝‖ＴＲ２‖２＝・・・＝ＴＲＭＩＮ      （３４）
となる。ここでＴＲＭＩＮは実数値である。
【０１１９】
上記式３２、もしくは式３４が成り立つときの（Ｎ−１）−単体Ｔとそれを構成する頂点との関係は、Δを格子の分割幅とすると、‖ＴＲ‖２の値に基づいて以下の３通りに分類される：
‖ＴＲ‖２≦（１／Ｎ）Δのとき（図１５）
（１／Ｎ）Δ＜‖ＴＲ‖２≦（Ｎ１／２−１／Ｎ）Δのとき（図１６）
（Ｎ１／２−１／Ｎ）Δ＜‖ＴＲ‖２のとき（図１７）      （３５）
以下、それぞれの場合について、ＴＲＭＩＮの具体的な計算方法について説明する。
【０１２０】
先ず、上記式３２、もしくは式３４が成り立つときの（Ｎ−１）−単体Ｔとそれを構成する頂点との関係が、‖ＴＲ‖２≦（１／Ｎ）Δを満たす場合について図１５を用いて説明する。
【０１２１】
図１５において９０１は節点Ｑ、９０２は頂点Ｖ１、９０３は頂点Ｖ２、９０４は頂点Ｖ３である。９０１から９０２、９０３、及び９０４へ向かうベクトルが、それぞれ、ＴＲ１、ＴＲ２、及びＴＲ３である。９０２、９０３、及び９０４で構成される三角形が（Ｎ−１）−単体Ｔであり、９０５は点Ａから（Ｎ−１）−単体Ｔへの垂線の足である。従って９０１から９０５へ向かうベクトルがＴＲである。上記式３３の値は、簡単な計算からＮ（空間の次元数）であることが分かる。即ち、２次元の場合は２、３次元の場合は３である。これよりＴＲ１は次式で与えられる。
【０１２２】
ＴＲ１＝（Ｎ‖ＴＲ‖２，０，０）  （３６）
これはＴＲ２、ＴＲ３についても同様である。
【０１２３】
次に、上記式３２、もしくは式３４が成り立つときの（Ｎ−１）−単体Ｔとそれを構成する頂点との関係が、（１／Ｎ）Δ＜‖ＴＲ‖２≦（Ｎ１／２−１／Ｎ）Δを満たす場合について図１６を用いて説明する。
【０１２４】
９１１は節点Ｑ、９１２は頂点Ｖ１、９１３は頂点Ｖ２、９１４は頂点Ｖ３、９１５は頂点Ｖ４、９１６は頂点Ｖ５、９１７は頂点Ｖ６である。点Ａから頂点Ｖ１〜Ｖ６へ向かうベクトルが、それぞれ、ＴＲ１、ＴＲ２、、、、、及びＴＲ６である。頂点Ｖ１〜Ｖ６で構成される六角形が（Ｎ−１）−単体Ｔであり、９１８は点Ａから（Ｎ−１）−単体Ｔへの垂線の足である。従って９１１から９１８へ向かうベクトルがＴＲである。頂点９１９、９２０、及び９２１で構成される三角形を考えると、節点Ｑから頂点９１９の距離はＮＴＲであり（図では３次元だから３ＴＲ）、立方体の頂点９２２と９１２の距離は９２２と９１９の距離に等しいので、ＮＴＲ−Δである。従って、節点Ｑから頂点Ｖ１へ向かうベクトルＴＲ１は
ＴＲ１＝（Δ，Ｎ‖ＴＲ‖２−Ｄ，０）  （３７）
これはＴＲ２〜ＴＲ６についても同様である。
【０１２５】
次に、上記式３２、もしくは式３４が成り立つときの（Ｎ−１）−単体Ｔとそれを構成する頂点との関係が、（Ｎ１／２−１／Ｎ）Δ＜‖ＴＲ‖２を満たす場合について図１７を用いて説明する。
【０１２６】
９２１は節点Ｑ、９２２は頂点Ｖ１、９２３は頂点Ｖ２、９２４は頂点Ｖ３である。点Ａから頂点Ｖ１〜Ｖ３へ向かうベクトルを、それぞれ、ＴＲ１、ＴＲ２、及びＴＲ３とする。頂点Ｖ１〜Ｖ３で構成される三角形が（Ｎ−１）−単体Ｔであり、９２５は点Ａから（Ｎ−１）−単体Ｔへの垂線の足である。従って９２１から９２５へ向かうベクトルがＴＲである。すると、頂点９２６と点９２５の距離はＮ１／２Δ−‖ＴＲ‖２であるので、頂点９２６と頂点９２２の距離（頂点９２６と頂点９２３、頂点９２６と頂点９２４も同じ）はＮ（Ｎ１／２Δ−‖ＴＲ‖２）であり（図では３次元だから３（３１／２Δ−‖ＴＲ‖２））、従って頂点９２７と９２２の距離（頂点９２８と頂点９２３、頂点９２９と頂点９２４も同じ）は、
Δ−Ｎ（Ｎ１／２Δ−‖ＴＲ‖２）＝（１−Ｎ３／２）Δ＋Ｎ‖ＴＲ‖２  （３８）
と計算できる。これより、点Ｑから頂点Ｖ１へのベクトルＴＲ１は次式で与えられる。
【０１２７】
ＴＲ１＝（（１−Ｎ３／２）Δ＋Ｎ‖ＴＲ‖２，Δ，Δ）    （３９）
これはＴＲ２、ＴＲ３についても同様である。
【０１２８】
以上のようにＴＲが与えられればＴＲ１がわかり、ＴＲ１がわかれば点Ａから頂点Ｖ１へ向かうベクトルＴＡは次式で計算できる。
【０１２９】
ＴＡ＝ＴＱ＋ＴＲ１    （４０）
近接頂点探索のための上限値ｄｓｕｐは、ＴＡのノルム‖ＴＡ‖２として与えられる。
【０１３０】
次に、最近接頂点への距離ｄｍｉｎが与えられたときのＴＱ、及びＴＲの計算方法について説明する。
【０１３１】
式３１を最大化する（Ｎ−１）−単体への距離を与えるベクトルの各成分値は式３２より相等しいので、Ｎ−１／２ｄｍｉｎとなる。そこで、Ｎ−１／２ｄｍｉｎをΔで割った商をｍ、余りをｒとすると以下の式が成り立つ。
【０１３２】
Ｎ−１／２ｄｍｉｎ＝ｍΔ＋ｒ    （４１）
すると、上述の説明において、点Ａから点Ｑへ向かうベクトルＴＱは
ＴＱ＝（ｍΔ，ｍΔ，ｍΔ）    （４２）
また点Ｑから（Ｎ−１）−単体への垂線の足へ向かうベクトルＴＲは
ＴＲ＝（ｒ，ｒ，ｒ）    （４３）
ここで、‖ＴＲ‖２はＮ１／２ｒである。
【０１３３】
以上をまとめると、
１．  式（４１）〜式（４３）を用いてＴＱ及びＴＲを計算する。
【０１３４】
２．  ＴＲの値に基づいて次式よりＴＲ１を計算する。
【０１３５】
ＴＲ１＝（Ｎ３／２ｒ，０，０）  ｉｆ  ｒ≦Ｎ−３／２Δ
＝（（１−Ｎ３／２）Δ＋Ｎ３／２ｒ，Δ，Δ）  ｅｌｓｅ  ｉｆ  （１−Ｎ−３／２）Δ＜ｒ
＝（Δ，Ｎ３／２ｒ−Δ，０）  ｏｔｈｅｒｗｉｓｅ
３  ＴＲ１とＴＱより式（４０）を用いてＴＡを計算する。
【０１３６】
４  上限値ｄｓｕｐは、ＴＡのノルム‖ＴＡ‖２として与えられる。
【０１３７】
また、図１３のステップＳ６３０６における処理については第１の実施形態と同じであるので、これに関する説明は省略する。
【０１３８】
［その他の実施形態］
また、本発明の目的は、前述した実施形態の機能を実現するソフトウェアのプログラムコードを記録した記録媒体（または記憶媒体）を、システムあるいは装置に供給し、そのシステムあるいは装置のコンピュータ（またはＣＰＵやＭＰＵ）が記録媒体に格納されたプログラムコードを読み出し実行することによっても、達成されることは言うまでもない。この場合、記録媒体から読み出されたプログラムコード自体が前述した実施形態の機能を実現することになり、そのプログラムコードを記録した記録媒体は本発明を構成することになる。
【０１３９】
また、コンピュータが読み出したプログラムコードを実行することにより、前述した実施形態の機能が実現されるだけでなく、そのプログラムコードの指示に基づき、コンピュータ上で稼働しているオペレーティングシステム（ＯＳ）などが実際の処理の一部または全部を行い、その処理によって前述した実施形態の機能が実現される場合も含まれることは言うまでもない。
【０１４０】
さらに、記録媒体から読み出されたプログラムコードが、コンピュータに挿入された機能拡張カードやコンピュータに接続された機能拡張ユニットに備わるメモリに書込まれた後、そのプログラムコードの指示に基づき、その機能拡張カードや機能拡張ユニットに備わるＣＰＵなどが実際の処理の一部または全部を行い、その処理によって前述した実施形態の機能が実現される場合も含まれることは言うまでもない。
【０１４１】
本発明を上記記録媒体に適用する場合、その記録媒体には、先に説明したフローチャートに対応するプログラムコードが格納されることになる。
【図面の簡単な説明】
【０１４２】
【図１】与えられた点Ａから界面までの符号付距離の計算処理のフローチャートである。
【図２】本発明の第１の実施形態に係る情報処理装置として機能するコンピュータの基本構成を示すブロック図である。
【図３】適合分割格子によるカットセル法について説明する図である。
【図４】カットセル法による三角分割の処理のフローチャートである。
【図５】ステップＳ４０２における処理を説明するための図である。
【図６】得られる多角形の種類を示す図である。
【図７】三角形を越える多角形の分割を説明する図である。
【図８】ステップＳ１０４における処理の詳細を示す図である。
【図９】ステップＳ１０４における処理の詳細を示すフローチャートである。
【図１０】３次元空間における三角形（２−単体）の符号付距離の計算処理の詳細を示す図である。
【図１１】３次元空間における三角形（２−単体）の符号付距離の計算処理のフローチャートである。
【図１２】“fast marching method”の処理のフローチャートである。
【図１３】与えられた点Ａから界面までの符号付距離を計算する処理のフローチャートである。
【図１４】２次元空間における、近接頂点探索のための上限値ｄｓｕｐを計算する処理を説明する図である。
【図１５】ＴＲＭＩＮの具体的な計算方法を説明する図である。
【図１６】ＴＲＭＩＮの具体的な計算方法を説明する図である。
【図１７】ＴＲＭＩＮの具体的な計算方法を説明する図である。
【図１】
2006285614.tif 000003
【図２】
2006285614.tif 000004
【図３】
2006285614.tif 000005
【図４】
2006285614.tif 000006
【図５】
2006285614.tif 000007
【図６】
2006285614.tif 000008
【図７】
2006285614.tif 000009
【図８】
2006285614.tif 000010
【図９】
2006285614.tif 000011
【図１０】
2006285614.tif 000012
【図１１】
2006285614.tif 000013
【図１２】
2006285614.tif 000014
【図１３】
2006285614.tif 000015
【図１４】
2006285614.tif 000016
【図１５】
2006285614.tif 000017
【図１６】
2006285614.tif 000018
【図１７】
2006285614.tif 000019