|TOP Page|


メッシュ生成のプログラミングTIPS
×

[PR]この広告は3ヶ月以上更新がないため表示されています。
ホームページを更新後24時間以内に表示されなくなります。

メッシュ生成のプログラミングTIPS

TIPS for programing mesh generation




Last Update:2008年02月17日


2次元の幾何情報取得プログラム

3角形の面積

p_1,p_2,p_3からなる三角形の面積Sは次のように表すことができる。

S(p_1,p_2,p_3)=\frac{1}{2}\|\begin{array} 1 & x_1 & y_1 \\ 1 & x_2 & y_2 \\ 1 & x_3 & y_3 \end{array}\|=\frac{1}{2}\|\begin{array} x_2-x_1 & y_2-y_1 \\ x_3-x_1 & y_3-y_1 \end{array}\|\\\;\;=\frac{1}{2}\{(x_2-x_1)*(y_3-y_1)-(x_3-x_1)*(y_2-y_1)\}

また、ヘロンの公式によれば、3角形の3辺の長さをa,b,cs=(a+b+c)/2としたとき

S=\sqrt{s(s-a)(s-b)(s-c)}

が成り立つ。

プログラミング例

double TriArea(const CVector2D& v1, const CVector2D& v2, const CVector2D& v3){
   return 0.5*( (v2.x-v1.x)*(v3.y-v1.y) - (v3.x-v1.x)*(v2.y-v1.y) );
}

内接円の半径

3角形の周の長さの半分をs、面積をSとすると、内接円の半径r

r=\frac{S}{s}

となる。

外接円の半径

三角形の面積をS、3辺の長さをa,b,cとすると、外接円の半径R

R=\frac{abc}{4S}

が成り立つ。

アスペクト比

3角形のアスペクト比Qは、外接円の半径R、内接円の半径rとすると次のように定義される。

Q=\frac{R}{r}

3角形の質を評価するためにアスペクト比がよく用いられる。三角形の三辺の長さをa,b,cs=(a+b+c)/2とおくと、ヘロンの公式を使って、アスペクト比Qは次のようになる。

Q=\frac{R}{r}=\frac{abc}{4S}\cdot\frac{s}{S}=\frac{abc}{4(s-a)(s-b)(s-c)}

外接円の中心

p_1,p_2,p_3で与えられる、3角形の外接円の中心q

a=|p_2-p_3|,\;b=|p_3-p_1|,\;c=|p_1-p_2|とすると、

q=\[a^2(b^2+c^2-a^2)p_1+b^2(c^2+a^2-b^2)p_2+c^2(a^2+b^2-c^2)p_3\]/(16S^2)

によって与えられる。但しSは3角形の面積

外接円の中心の求め方は、こちらを参考にした。

ベクトル的解法のヒント
http://robotics.me.es.osaka-u.ac.jp/~masutani/Jouki/Kyouzai/hint-vector2d.html

プログラミング例

bool CenterCircumcircle(
      const CVector2D& p0, 
      const CVector2D& p1, 
      const CVector2D& p2, 
      CVector2D& center){

   const double area = TriArea(p0,p1,p2);
   if( fabs(area) < 1.0e-10 ){ return false; }
   const double tmp_val = 1.0/(area*area*16.0);

   const double dtmp0 = SquareLength(p1,p2);
   const double dtmp1 = SquareLength(p0,p2);
   const double dtmp2 = SquareLength(p0,p1);

   const double etmp0 = tmp_val*dtmp0*(dtmp1+dtmp2-dtmp0);
   const double etmp1 = tmp_val*dtmp1*(dtmp0+dtmp2-dtmp1);
   const double etmp2 = tmp_val*dtmp2*(dtmp0+dtmp1-dtmp2);

   center.x = etmp0*p0.x + etmp1*p1.x + etmp2*p2.x;
   center.y = etmp0*p0.y + etmp1*p1.y + etmp2*p2.y;
   return true
}

3角形のドロネー条件判定

2次元3角形分割の場合はドロネー分割がある種の最適な分割である。ドロネー三角分割とは、3角形の外接円に他の三角形の頂点が入っていないことである。これをチェックするプログラムを考えよう。

Delaunay条件を満たす三角形 Delaunay条件を満たさない三角形

p_1,p_2,p_3を頂点とする三角形について、

H(x,y,z)=\|\begin{array} 1 & x_1 & y_1 & x_1^2+y_1^2 \\ 1 & x_2 & y_2 & x_2^2+y_2^2 \\ 1 & x_3 & y_3 & x_3^2+y_3^2 \\ 1 & x & y & x^2+y^2 \end{array}\|

となる。さて、点p_4がこの外接円に入っているかどうか調べよう。これにはH(p_4)の符号を調べればよい。

余因子展開を用いてこれは次のように書ける。

\|\begin{array} 1 & x_1 & y_1 & x_1^2+y_1^2 \\ 1 & x_2 & y_2 & x_2^2+y_2^2 \\ 1 & x_3 & y_3 & x_3^2+y_3^2 \\ 1 & x_4 & y_4 & x_4^2+y_4^2 \end{array}\|   =   \|\begin{array} x_1-x_4 & y_1-y_4 & (x_1-x_4)(x_1+x_4)+(y_1-y_4)(y_1+y_4) \\ x_2-x_4 & y_2-y_4 & (x_2-x_4)(x_2+x_4)+(y_2-y_4)(y_2+y_4) \\ x_3-x_4 & y_3-y_4 & (x_3-x_4)(x_3+x_4)+(y_3-y_4)(y_3+y_4) \end{array}\|\\\;\;    =    +\{\bar{x}_1(x_1+x_4)+\bar{y}_1(y_1+y_4)\}\|\begin{array}\bar{x}_2 & \bar{y}_2\\ \bar{x}_3 & \bar{y}_3 \end{array}\|\\\;\;\;\;-\{\bar{x}_2(x_2+x_4)+\bar{y}_2(y_2+y_4)\}\|\begin{array}\bar{x}_1 & \bar{y}_1\\ \bar{x}_3 & \bar{y}_3 \end{array}\|\\\;\;\;\; +\{\bar{x}_3(x_3+x_4)+\bar{y}_3(y_3+y_4)\}\|\begin{array}\bar{x}_1 & \bar{y}_1\\ \bar{x}_2 & \bar{y}_2 \end{array}\|

但し\bar{x}_i=x_i-x_4,\;\;\bar{y}_i=y_i-y_4である。

プログラミング例

double DetDelaunay(
      const CVector2D& p0, 
      const CVector2D& p1, 
      const CVector2D& p2, 
      const CVector2D& p3){

   a[6] = {
      p0.x-p3.x // 0
      p1.x-p3.x // 1
      p2.x-p3.x // 2
      p0.y-p3.y // 3
      p1.y-p3.y // 4
      p2.y-p3.y // 5
   };
   return +(a[0]*(p0.x+p3.x)+a[3]*(p0.y+p3.y))*(a[1]*a[5]-a[2]*a[4])
          -(a[1]*(p1.x+p3.x)+a[4]*(p1.y+p3.y))*(a[0]*a[5]-a[2]*a[3])
          +(a[2]*(p2.x+p3.x)+a[5]*(p2.y+p3.y))*(a[0]*a[4]-a[1]*a[3]);
}

3次元の幾何情報取得プログラム

3角形の面積

外積を利用して、面積を求める。

S=\frac{1}{2}|(p_2-p_1)\times (p_3-p_1)|

プログラミング例

double TriArea(const CVector3D& v1, const CVector3D& v2, const CVector3D& v3)
{
   double x, y, z;
   x = ( v2.y - v1.y )*( v3.z - v1.z ) - ( v3.y - v1.y )*( v2.z - v1.z );
   y = ( v2.z - v1.z )*( v3.x - v1.x ) - ( v3.z - v1.z )*( v2.x - v1.x );
   z = ( v2.x - v1.x )*( v3.y - v1.y ) - ( v3.x - v1.x )*( v2.y - v1.y );
   return 0.5*sqrt( x*x + y*y + z*z );
}

四面体の体積

p_1,p_2,p_3,p_4からなる四面体の体積Vは次のとおり。

V(p_1,p_2,p_3,p_4)=\frac{1}{6}\|\begin{array} 1 & x_1 & y_1 & z_1 \\ 1 & x_2 & y_2 & z_2 \\ 1 & x_3 & y_3 & z_3 \\ 1 & x_4 & y_4 & z_4 \end{array}\| = \frac{1}{6}\|\begin{array}x_2-x_1 & y_2-y_1 & z_2-z_1\\ x_3-x_1 & y_3-y_1 & z_3-z_1 \\ x_4-x_1 & y_4-y_1 & z_4-z_1 \end{array}\|\\\;\; = \frac{1}{6}\left\{(x_2-x_1)\|\begin{array}y_3-y_1 & z_3-z_1\\ y_4-y_1 & z_4-z_1\end{array}\|-(y_2-y_1)\|\begin{array}x_3-x_1 & z_3-z_1\\ x_4-x_1 & z_4-z_1\end{array}\|+(z_2-z_1)\|\begin{array}x_3-x_1 & y_3-y_1\\ x_4-x_1 & y_4-y_1\end{array}\|\right\}

プログラミング例

double TetVolume(
      const CVector3D& v1,
      const CVector3D& v2, 
      const CVector3D& v3, 
      const CVector3D& v4 )
{
   return	
      (   ( v2.x - v1.x )*( ( v3.y - v1.y )*( v4.z - v1.z ) - ( v4.y - v1.y )*( v3.z - v1.z ) )		
        - ( v2.y - v1.y )*( ( v3.x - v1.x )*( v4.z - v1.z ) - ( v4.x - v1.x )*( v3.z - v1.z ) )		
        + ( v2.z - v1.z )*( ( v3.x - v1.x )*( v4.y - v1.y ) - ( v4.x - v1.x )*( v3.y - v1.y ) )
      ) * 0.16666666666666666666666666666667;		
}

四面体の内接球の半径

四面体の体積をV、表面積をSとする。

内接球の半径をr

r=\frac{V}{3S}

となる。

プログラミング例

double InscribedRadius(
      const CVector3D& ipo0, 
      const CVector3D& ipo1, 
      const CVector3D& ipo2, 
      const CVector3D& ipo3)
   double surface = TriArea( ipo1, ipo3, ipo2 )
      + TriArea( ipo0, ipo2, ipo3 )
      + TriArea( ipo0, ipo3, ipo1 )
      + TriArea( ipo0, ipo1, ipo2 );
   return TetVolume(ipo0, ipo1, ipo2, ipo3) * 3.0 / surface;
}

四面体の外接球の半径

外心の座標を求めるために、外接円の中心の座標を求めて、頂点からの距離を計算する。

外心の座標の求め方についてはこちらを参考にした。

三角形、四面体との相似
http://homepage2.nifty.com/PAF00305/math/triangle/node7.html

プログラミング例

inline double Circumradius(
      const CVector3D& ipo0, 
      const CVector3D& ipo1, 
      const CVector3D& ipo2, 
      const CVector3D& ipo3)
{
   double base[3][3] = {
		{ ipo1.x-ipo0.x, ipo1.y-ipo0.y, ipo1.z-ipo0.z },
		{ ipo2.x-ipo0.x, ipo2.y-ipo0.y, ipo2.z-ipo0.z },
		{ ipo3.x-ipo0.x, ipo3.y-ipo0.y, ipo3.z-ipo0.z }
   };
   double s[6] = {
      base[0][0]*base[0][0]+base[0][1]*base[0][1]+base[0][2]*base[0][2],
      base[1][0]*base[1][0]+base[1][1]*base[1][1]+base[1][2]*base[1][2],
      base[2][0]*base[2][0]+base[2][1]*base[2][1]+base[2][2]*base[2][2],
      base[1][0]*base[2][0]+base[1][1]*base[2][1]+base[1][2]*base[2][2],
      base[2][0]*base[0][0]+base[2][1]*base[0][1]+base[2][2]*base[0][2],
      base[0][0]*base[1][0]+base[0][1]*base[1][1]+base[0][2]*base[1][2],
   };
   const double vol = TetVolume(ipo0,ipo1,ipo2,ipo3)*6.0;
   if( vol < 1.0e-20 ){
      std::cout << "Error!-->Too Small " << std::endl;
      assert(0);
   }
   const double inv_det = 1.0 / (vol*vol);
   double t[6] = {
      (s[1]*s[2]-s[3]*s[3])*0.5*inv_det,
      (s[2]*s[0]-s[4]*s[4])*0.5*inv_det,
      (s[0]*s[1]-s[5]*s[5])*0.5*inv_det,
      (s[4]*s[5]-s[0]*s[3])*0.5*inv_det,
      (s[5]*s[3]-s[1]*s[4])*0.5*inv_det,
      (s[3]*s[4]-s[2]*s[5])*0.5*inv_det,
   };
   double u[3] = {
      t[0]*s[0]+t[5]*s[1]+t[4]*s[2],
      t[5]*s[0]+t[1]*s[1]+t[3]*s[2],
      t[4]*s[0]+t[3]*s[1]+t[2]*s[2],
   };
   return sqrt(0.5*(u[0]*s[0]+u[1]*s[1]+u[2]*s[2]));
}

次の論文のように、直接外接球の半径を求める方法もあるらしい。ただ、2次元の場合と違ってあまり単純にはならない。

四面体の外接球の半径について   五十嵐貫 佐藤義隆 杉江道男 高遠節夫 豊成敏隆 向山一男
http://ci.nii.ac.jp/naid/110000198608/

四面体のドロネー条件判定

p_1,p_2,p_3,p_4を頂点とする四面体について、

H(x,y,z)=\|\begin{array} 1 & x_1 & y_1 & z_1 & x_1^2+y_1^2+z_1^2 \\ 1 & x_2 & y_2 & z_2 & x_2^2+y_2^2+z_2^2 \\ 1 & x_3 & y_3 & z_3 & x_3^2+y_3^2+z_3^2 \\ 1 & x_4 & y_4 & z_4 & x_4^2+y_4^2+z_4^2 \\ 1 & x & y & z & x^2+y^2+z^2 \end{array}\|

  1. H(x,y,z)<0のとき四面体の外接球の内側の方程式
  2. H(x,y,z)=0のとき四面体の外接球の方程式
  3. H(x,y,z)>0のとき四面体の外接球の外側の方程式

を与える。

p_5がこの四面体の中に入っているかどうか調べるためには次の符号を調べればよい。

H(p_5)    =    \|\begin{array} 1 & x_1 & y_1 & z_1 & x_1^2+y_1^2+z_1^2 \\ 1 & x_2 & y_2 & z_2 & x_2^2+y_2^2+z_2^2 \\ 1 & x_3 & y_3 & z_3 & x_3^2+y_3^2+z_3^2 \\ 1 & x_4 & y_4 & z_4 & x_4^2+y_4^2+z_4^2 \\ 1 & x_5 & y_5 & z_5 & x_5^2+y_5^2+z_5^2 \end{array}\|\\\;\;\;    =   \|\begin{array} x_1-x_5 & y_1-y_5 & z_1-z_5 & (x_1-x_5)(x_1+x_5)+(y_1-y_5)(y_1+y_5)+(z_1-z_5)(z_1+z_5) \\ x_2-x_5 & y_2-y_5 & z_2-z_5 & (x_2-x_5)(x_2+x_5)+(y_2-y_5)(y_2+y_5)+(z_2-z_5)(z_2+z_5) \\ x_3-x_5 & y_3-y_5 & z_3-z_5 & (x_3-x_5)(x_3+x_5)+(y_3-y_5)(y_3+y_5)+(z_3-z_5)(z_3+z_5) \\ x_4-x_5 & y_4-y_5 & z_4-z_5 & (x_4-x_5)(x_4+x_5)+(y_4-y_5)(y_4+y_5)+(z_4-z_5)(z_4+z_5) \end{array}\|\\\;\;\;    =     -\{\bar{x}_1(x_1+x_5)+\bar{y}_1(y_1+y_5)+\bar{z}_1(z_1+z_5)\}\|\begin{array} \bar{x}_2 & \bar{y}_2 & \bar{z}_2 & \\ \bar{x}_3 & \bar{y}_3 & \bar{z}_3 \\ \bar{x}_4 & \bar{y}_4 & \bar{z}_4\end{array}\|\\\;\;\;\;\;\;   +\{\bar{x}_2(x_2+x_5)+\bar{y}_2(y_2+y_5)+\bar{z}_2(z_2+z_5)\}\|\begin{array} \bar{x}_1 & \bar{y}_1 & \bar{z}_1 & \\ \bar{x}_3 & \bar{y}_3 & \bar{z}_3 \\ \bar{x}_4 & \bar{y}_4 & \bar{z}_4\end{array}\|\\\;\;\;\;\;\;   -\{\bar{x}_3(x_3+x_5)+\bar{y}_3(y_3+y_5)+\bar{z}_3(z_3+z_5)\}\|\begin{array} \bar{x}_1 & \bar{y}_1 & \bar{z}_1 & \\ \bar{x}_2 & \bar{y}_2 & \bar{z}_2 \\ \bar{x}_4 & \bar{y}_4 & \bar{z}_4\end{array}\|\\\;\;\;\;\;    -\{\bar{x}_4(x_4+x_5)+\bar{y}_4(y_4+y_5)+\bar{z}_4(z_4+z_5)\}\|\begin{array} \bar{x}_1 & \bar{y}_1 & \bar{z}_1 & \\ \bar{x}_2 & \bar{y}_2 & \bar{z}_2 \\ \bar{x}_3 & \bar{y}_3 & \bar{z}_3\end{array}\|

プログラミング例

double DetDelaunay3D(
      const CVector3D& v1, 
      const CVector3D& v2, 
      const CVector3D& v3, 
      const CVector3D& v4, 
      const CVector3D& v5)
{
   const double a[12] = {
      v1.x-v5.x,	//  0
      v2.x-v5.x,	//  1
      v3.x-v5.x,	//  2
      v4.x-v5.x,	//  3
      v1.y-v5.y,	//  4
      v2.y-v5.y,	//  5
      v3.y-v5.y,	//  6
      v4.y-v5.y,	//  7
      v1.z-v5.z,	//  8
      v2.z-v5.z,	//  9
      v3.z-v5.z,	// 10
      v4.z-v5.z,	// 11
   };
   const double b[6] = {
      a[ 6]*a[11]-a[ 7]*a[10],	// 0
      a[ 5]*a[11]-a[ 7]*a[ 9],	// 1
      a[ 5]*a[10]-a[ 6]*a[ 9],	// 2
      a[ 7]*a[ 8]-a[ 4]*a[11],	// 3
      a[ 6]*a[ 8]-a[ 4]*a[10],	// 4
      a[ 4]*a[ 9]-a[ 5]*a[ 8],	// 5 
   };
   return -( a[0]*(v1.x+v5.x)+a[4]*(v1.y+v5.y)+a[ 8]*(v1.z+v5.z) )*( a[ 1]*b[0]-a[ 2]*b[1]+a[ 3]*b[2] )
          +( a[1]*(v2.x+v5.x)+a[5]*(v2.y+v5.y)+a[ 9]*(v2.z+v5.z) )*( a[ 0]*b[0]+a[ 2]*b[3]-a[ 3]*b[4] )
          -( a[2]*(v3.x+v5.x)+a[6]*(v3.y+v5.y)+a[10]*(v3.z+v5.z) )*( a[ 0]*b[1]+a[ 1]*b[3]+a[ 3]*b[5] )
          +( a[3]*(v4.x+v5.x)+a[7]*(v4.y+v5.y)+a[11]*(v4.z+v5.z) )*( a[ 0]*b[2]+a[ 1]*b[4]+a[ 2]*b[5] );
}

参考にしたもの

LINK

四面体と3角形の類似、四面体の内心、外心、垂心の求め方が載っている。
http://homepage2.nifty.com/PAF00305/math/triangle/node7.html
ベクトル的解法のヒント、3角形の内接円、外接円、線に接する円の求め方が載っている。
http://robotics.me.es.osaka-u.ac.jp/~masutani/Jouki/Kyouzai/hint-vector2d.html
3角形の面積を紐解く
http://www.nikonet.or.jp/spring/menseki/menseki.pdf


Made by Nobuyuki UMETANI  梅谷 信行
n.umetani@gmail.com