JP2005202633A

JP2005202633A - 界面推定方法

Info

Publication number: JP2005202633A
Application number: JP2004007632A
Authority: JP
Inventors: Takahiko Tanahashi; ▲隆▼彦棚橋; Junya Imamura; 純也今村
Original assignee: Keio University
Current assignee: Keio University
Priority date: 2004-01-15
Filing date: 2004-01-15
Publication date: 2005-07-28

Abstract

【課題】設計寸法を押さえながら、ドット情報を取り込んで、境界形状を表現できるようにする。
【解決手段】実境界で区切られた３次元空間の状態値を格子点の値を未知パラメータとしてボクセルにより有限要素表現する。実境界上の所望の点について、状態値が０となるようにパッチを当てて状態値を補綴する。状態値をｆ(x,y,z)とし、内部状態が観測される確率Ｐ(x,y,z)を、Ｐ(x,y,z)＝１／{１＋exp(−ｆ(x,y,z))}とする。外部状態が観測される確率を(１−Ｐ(x,y,z))とし、実境界外部の複数のサンプリング点における各確率(１−Ｐ(x,y,z))と実境界内部の複数のサンプリング点における各確率Ｐ(x,y,z)との積が極大となる未知パラメータの値を求める。求めたパラメータにより状態値ｆ(x,y,z)を求める。状態値ｆ(x,y,z)から確率Ｐ(x,y,z)を求める。確率Ｐ(x,y,z)が１／２となる曲面を求める。その曲面を界面と推定する。
【選択図】図１

Description

本発明は、界面推定方法に関し、特に、確率密度関数の回帰モデルを有限要素表現して、ドット情報の分布から界面を推定する界面推定方法に関する。

界面の表現法は、有限要素法および境界要素法による形状表現が基本となっている。格子をベースとする有限要素法系列の第１の方法（多面体法）と第２の方法（ボクセル法）と、格子を用いない境界要素法系列の第３の方法（デルタ関数法）がある。第１の方法は、多面体の頂点座標値による表現方法であり、ポリゴンデータやマップトメッシュなどの表現方法が対応する。第２の方法は、ボクセル表現方法であり、ボクセルの属性を、物質なしか物質ありのいずれかとする方法である。第３の方法は、デルタ関数などの和として表現する方法である。

観測データから回帰関数を求める手法の１つにロジットモデルがある。ロジットモデルについては、多変量解析の教科書を参照されたい。ロジットモデルには、２項ロジット（binary logit：BL）モデルと、多項ロジット（multinominal logit：ML）モデルと、ネステッドロジット（nested logit：NL）モデルがある。単純なBLモデルでは項の数（事象の数）は２であり、モデルで予測したい変数の１次結合式を効用関数と呼び、ｉを事象番号（１と２）としてＶ_iで表わす。Ｖ_iは回帰式に相当する。それぞれの事象の確率Ｐ_iは、Ｖ_iの指数比で表わされると仮定する。
Ｐ_i＝exp(Ｖ_i)／Σ_i{exp(Ｖ_i)}
確率密度関数はロジスティック曲線となる。サンプリング番号をｎとし、サンプリング点での事象をδ_inとする。具体的には物体内なら（δ_1n＝１,δ_2n＝０）とし、物体外なら（δ_1n＝０,δ_2n＝１）とする。同時生起の確率をＬ^*として積の形で表わす。
Ｌ^*＝Π_nΠ_iＰ_in ^δin＝Π_n(Ｐ_1n ^δ1nＰ_2n ^δ2n)
δ_1n＝１のときδ_2n＝０、またδ_1n＝０のときδ_2n＝１であるから結局サンプリング点の事象の確率をＰ_inとして次の式になる。
Ｌ^*＝Π_nＰ_in
Ｖ_iの回帰係数をパラメータとして、Ｌ^*を尤度関数とし、同時生起の確率が最大となるようなパラメータを最尤推定する。

事象１と事象２が排他事象であるとすると、
Ｐ₁＝exp(Ｖ₁)／{exp(Ｖ₁)＋exp(Ｖ₂)}
Ｐ₁＋Ｐ₂＝１
Ｐ₁＝１／{１＋exp(Ｖ₂−Ｖ₁)}
＝１／{１＋exp(−ΔＶ)}
となる。

本発明者は、ロジットモデルを利用して境界を推定する方法をすでに提案した。この方法についての詳細は、非特許文献１〜３を参照されたい。ロジットモデルにボクセル解析法を組み合わせる方法である。２次元の場合は、空間をピクセルで表現する。数値計算の対象は、ボクセル群のうち、境界を含むボクセルのみとなるので、アダプティブ手法の考え方で、境界要素のみ細かい格子として精密計算できる。ロジットモデルによる境界形状表現は、入力データが、スキャナーやデジタルカメラなどから取得したドット情報の場合に適用して便利である。

２次元の場合、サンプリング点情報（ｘ,ｙ座標）を独立変数とし、回帰式を有限要素関数とし、場をあらわす有限要素関数の値をΔＶとする。図６に示すように、界面を含むピクセル要素のみ取り出し、要素内サンプリング点が境界内にあれば確率をＰとし、
Ｐ＝１／{１＋exp(−ΔＶ)}
境界外ならば確率を(１−Ｐ)とする。要素内で、それらの同時確率式を表わし、最尤推定式を連立させて解く。図７に示すように、ΔＶ＝0.0を境界面として表わす。

３次元の場合、変数は座標値(ｘ,ｙ,ｚ)となる。ΔＶを有限要素表現する。効用関数差ΔＶを、最も簡単なボクセルで表現する。パラメータは有限要素法に基づきノード値で表現する。ΔＶのパラメータをαとし、同時生起の確率の尤度関数Ｌ^*とする。サンプリング点における確率は、系外では(１−Ｐ)となり、系内ではＰとなる。Ｌ^*を最大とするようなαを最尤推定する。Ｌ^*を最大にすることと、その対数Ｌを最大にすることは等価である。すなわち、δＬ＝０とするパラメータを推定すればよい。非線形式であるから、Hesse行列で表わして、反復計算して求める。

ΔＶに要求されるのは、Ｃ₀-連続性である。Ｃ₀-連続な最少自由度のボクセルは、よく知られた３重１次要素である。境界形状は、Ｃ₀-連続よりＣ₁-連続な要素で表現する方が、視覚的により滑らかである。Ｃ₀-要素を球体に適用して、数値計算した結果を図８に示す。１×１×１の３次元空間に、直径0.96の球体を置き、等分割ボクセルで近似したものである。分割数は、(a)２×２×２，(b)４×４×４，（c）10×10×10，(d)20×20×20，(e)60×60×60とし、(f)には、オリジナルの球断面を円弧でプロットした。いずれも、中心からピッチ0.05で９つの断面を重ね描きしたものである。このスケールの図で(d)(e)の分割では、オリジナルとの差はほとんど識別できない。
今村、"自由表面問題におけるVOF関数へのロジットモデルの適用"、第47回応用力学連合講演会予稿集, (1998), p.37. 今村、"プロビットモデルへの有限要素法の適用"、土木学会大会梗概集（講演No.41380）,(1999). 今村純也、″非集計行動モデルと計算工学の考察″、計算工学講演会論文集(計算工学会)、Vol.8, No.1 (2003年5月)、pp.363-364.

しかし、従来の方法では、境界上の特定点データがわかっていても利用することができないので、精度の高い界面の推定ができないという問題があった。本発明は、上記従来の問題を解決して、特定点の設計寸法を押さえながら、ドット情報を取り込んで、境界形状を表現できるようにすることを目的とする。

上記の課題を解決するために、本発明では、境界推定方法を、実境界で区切られた平面の状態を格子点の値を未知パラメータとして正方格子により有限要素表現し、実境界上の所望の点について、状態値が０となるようにパッチを当てて状態値を補綴し、状態値をｆ(x,y)とし、内部状態が観測される確率Ｐ(x,y)を、
Ｐ(x,y)＝１／{１＋exp(−ｆ(x,y))}
とし、外部状態が観測される確率を(１−Ｐ(x,y))とし、実境界外部の複数のサンプリング点における各確率(１−Ｐ(x,y))と実境界内部の複数のサンプリング点における各確率Ｐ(x,y)との積が極大となる未知パラメータの値を求め、求めたパラメータにより状態値ｆ(x,y)を求め、状態値ｆ(x,y)から確率Ｐ(x,y)を求め、確率Ｐ(x,y)が１／２となる曲線を求め、その曲線を境界と推定する方法とした。

本発明では、上記のように構成したことにより、ボクセル解法の利点であるデータ入力の容易さを損なわずに、精度高い界面の推定が可能となる。境界が時系列に移動する問題のオイラー型解析も可能となる。時間軸を変数に含めれば、接触問題の解法にも応用できる。設計情報などのように、界面の一部が確定した数値（押さえ寸法）として得られる場合には、六面体パッチ要素法と組み合わせて、押さえ寸法を制約条件とすることで、物体のより正確な表現が可能となる。

以下、本発明を実施するための最良の形態について、図１〜図５を参照しながら詳細に説明する。

本発明の実施例は、ロジットモデルにボクセル解析法と有限パッチ要素法を組み合わせて、最尤推定法で界面を推定する界面推定方法である。

本発明の実施例における界面推定方法の手順を説明する。図１は、本発明の実施例における界面推定方法の処理手順を示す図である。図２と図３と図４は、２次元でパッチを当てる方法を示す図である。図５は、円柱にパッチを当てる方法を示す図である。

最初に、図１を参照しながら、界面推定方法の処理手順の概略を説明する。（１）実境界で区切られた３次元空間の状態値を格子点の値を未知パラメータとしてボクセルにより有限要素表現する。
（２）実境界上の所望の点について、状態値が０となるようにパッチを当てて状態値を補綴する。３次元では、８要素１組のパッチ系とする。稜線の境界面貫通点にパッチノードが位置するよう当てる。境界ボクセルの隅を三角錐の形に切取る場合、ボクセル稜線の境界面貫通点は３点あるので、パッチ系は３組必要となる。
（３）サンプリング点における確率を状態値から計算する。
（４）同時生起の確率として、各確率の積を計算する。
（５）実境界外部の複数のサンプリング点における各確率と実境界内部の複数のサンプリング点における各確率との積が極大となる未知パラメータの値を求める。
（６）求めたパラメータにより状態値を求める。状態値から確率を求める。確率が１／２となる曲面を求める。その曲面を境界と推定する。

ここで、用語の説明をする。境界に接するか境界を含む要素一般を、境界要素と呼ぶ。実際の境界として与えられた境界を実境界（real boundary）と呼ぶ。例えば、測定点や設計点を結んで境界としたものなどである。ボクセルやピクセルの外側境界を仮想境界（virtual boundary）と呼ぶ。実境界内の要素領域を実領域（real domain）と呼ぶ。実境界外の領域を仮想領域（virtual domain）と呼ぶ。有限パッチ要素法では、基本系(base system)の基底関数ｆaを、パッチ要素ｆbで、ｆ＝ｆa＋ｆbの和形式の形に補綴する。重合メッシュ法の考えと同じで、ｆb群の外周ノード値はゼロとする。Ｃ₁-要素では、外周勾配はゼロとする。基本系を補綴する系を、パッチ系（patch-system）と呼ぶ。

最も簡単なパッチ系は、２次元では田の字モデルである。３次元では８個の立方体からなるパッチ系８要素である。中心ノードのみ未知数とし、周辺ノードはすべてゼロとする関数を用いる。これらを１つまたは複数個組み合わせて、パッチ系を構成する。基本系の格子に対しては、パッチ系の位置する格子のみに影響し、他の格子には何ら影響を及ぼさない。実境界と基本系の要素辺との交点を選び、その点における基本系の関数値をゼロとするように、パッチ系の中心ノード値を与えて、パッチ系の中心ノードを、その点に重ねれば、少なくともこの点では境界条件を満足する。同様にして、他の辺と実境界との交点でも補綴する。

図２は、和形式の簡易モデルである。境界ボクセルをパッチ系で補綴して、境界条件を近似的に満足させることを、２次元で説明するための図である。実境界が通る基本系のボクセルに細かく区切った形のパッチ系を当てる。網掛けの部分がパッチ系である。パッチ要素のノード値は、境界ボクセルと内部ボクセル間の格子上ノード（白抜きの四角で示すノード）でゼロとする。仮想領域内ノード（黒塗りの四角で示すノード）には、基本系の値の符号を逆にして与えて相殺する。実領域内ノード（黒塗りの丸で示すノード）は、ガレルキン法で決定する。

これを変形した方法を、図３に示す。実境界を含むパッチ要素の実領域内ノード（白抜きの四角で示すノード）には、ゼロを与える。仮想領域内ノード（白抜きの丸で示すノード）には、基本系の値を逆符号で与えて相殺する。薄い網掛けの部分は、相殺された結果としてゼロ値であるから、積分不要となる。パッチ系のノードに基本系の値を逆符号で与えるときの値は、基本系のノードパラメータ（未知数）で表現できる。近似モデルの連立方程式は、結局基本系のみで記述される。

図４は、押さえ点にパッチを当てる方法を示す図である。この方法では、設計寸法を押さえて、曲線や曲面で滑らかに補間することができる。境界要素内をいずれかの軸が押さえ点を通るようにして、サブ境界要素に分割する。ボクセル内に境界の極値が現れないよう、境界の各空間軸極値点にノード設定する。各サブ境界要素にパッチを当てる。図５に、直径Ｄの円柱周りをＤ/6×Ｄ/6のピクセル格子で切った基本系に、和形式で補綴する方法を示す。

図６で２次元の場合の境界推定手順を説明する。実境界で区切られた平面の状態を格子点の値を未知パラメータとして正方格子により有限要素表現する。実境界上の所望の点について、状態値が０となるようにパッチを当てて状態値を補綴する。状態値をｆ(x,y)とし、内部状態が観測される確率Ｐ(x,y)を、
Ｐ(x,y)＝１／{１＋exp(−ｆ(x,y))}
とする。外部状態が観測される確率を(１−Ｐ(x,y))とする。実境界外部の複数のサンプリング点における各確率(１−Ｐ(x,y))と実境界内部の複数のサンプリング点における各確率Ｐ(x,y)との積が極大となる未知パラメータの値を求める。求めたパラメータにより状態値ｆ(x,y)を求める。状態値ｆ(x,y)から確率Ｐ(x,y)を求める。確率Ｐ(x,y)が１／２となる曲線を求める。その曲線を境界と推定する。

３次元の場合の界面推定手順を説明する。実境界で区切られた３次元空間の状態を格子点の値を未知パラメータとしてボクセルにより有限要素表現する。実境界上の所望の点について、状態値が０となるようにパッチを当てて状態値を補綴する。状態値をｆ(x,y,z)とし、内部状態が観測される確率Ｐ(x,y,z)を、
Ｐ(x,y,z)＝１／{１＋exp(−ｆ(x,y,z))}
とする。外部状態が観測される確率を(１−Ｐ(x,y,z))とする。実境界外部の複数のサンプリング点における各確率(１−Ｐ(x,y,z))と実境界内部の複数のサンプリング点における各確率Ｐ(x,y,z)との積が極大となる未知パラメータの値を求める。求めたパラメータにより状態値ｆ(x,y,z)を求める。状態値ｆ(x,y,z)から確率Ｐ(x,y,z)を求める。確率Ｐ(x,y,z)が１／２となる曲面を求める。その曲面を界面と推定する。

上記のように、本発明の実施例では、境界推定方法を、ロジットモデルにボクセル解析法と有限パッチ要素法を組み合わせて、最尤推定法で境界を推定する方法としたので、物体のより正確な表現が可能となる。

本発明の境界推定方法は、押さえ寸法を制約条件とする方法と組み合わせているので、設計情報の精密な補間ができ、CAD／CAMで利用する境界推定方法として最適である。クレイモデル等から得られるデジタル情報を画面化する技法としても応用でき、且つ設計情報としての特定点の値を与えることもできる。

本発明の実施例における境界推定方法の処理手順を示す図、本発明の実施例における境界推定方法で、パッチを当てる方法を示す図、本発明の実施例における境界推定方法で、パッチを当てる方法を示す図、本発明の実施例における境界推定方法で、パッチを当てる方法を示す図、本発明の実施例における境界推定方法で、円柱にパッチを当てる方法を示す図、従来の境界推定方法による計算方法を説明する図、従来の境界推定方法に境界推定を説明する図、従来の境界推定方法による計算例を示す図である。

Claims

実境界で区切られた平面の状態値を格子点の値を未知パラメータとして正方格子により有限要素表現し、実境界上の所望の点について、状態値が０となるようにパッチを当てて状態値を補綴し、状態値をｆ(x,y)とし、内部状態が観測される確率Ｐ(x,y)を、
Ｐ(x,y)＝１／{１＋exp(−ｆ(x,y))}
とし、外部状態が観測される確率を(１−Ｐ(x,y))とし、実境界外部の複数のサンプリング点における各確率(１−Ｐ(x,y))と実境界内部の複数のサンプリング点における各確率Ｐ(x,y)との積が極大となる未知パラメータの値を求め、求めたパラメータにより状態値ｆ(x,y)を求め、状態値ｆ(x,y)から確率Ｐ(x,y)を求め、確率Ｐ(x,y)が１／２となる曲線を求め、その曲線を境界と推定することを特徴とする境界推定方法。
実境界で区切られた３次元空間の状態値を格子点の値を未知パラメータとしてボクセルにより有限要素表現し、実境界上の所望の点について、状態値が０となるようにパッチを当てて状態値を補綴し、状態値をｆ(x,y,z)とし、内部状態が観測される確率Ｐ(x,y,z)を、
Ｐ(x,y,z)＝１／{１＋exp(−ｆ(x,y,z))}
とし、外部状態が観測される確率を(１−Ｐ(x,y,z))し、実境界外部の複数のサンプリング点における各確率(１−Ｐ(x,y,z))と実境界内部の複数のサンプリング点における各確率Ｐ(x,y,z)との積が極大となる未知パラメータの値を求め、求めたパラメータにより状態値ｆ(x,y,z)を求め、状態値ｆ(x,y,z)から確率Ｐ(x,y,z)を求め、確率Ｐ(x,y,z)が１／２となる曲面を求め、その曲面を界面と推定することを特徴とする界面推定方法。