delaunaydoc.cpp

三角网剖分 将平面数据点进行三角剖分 快速构建delaunay三角网

	  l1=l1*l1;
	  l2=l2*l2;
      w=w+double(1)/(l1*l2);//w ,Fx,Fy and n have been 初始化了	=0  
	  /*n.x=n.x+GetTriNormal(Tri).x/(l1*l2);
	  n.y=n.y+GetTriNormal(Tri).y/(l1*l2);
	  n.z=n.z+GetTriNormal(Tri).z/(l1*l2);*/
	  Fx=(-a[i]/c[i])/(l1*l2)+Fx;Fy=(-b[i]/c[i])/(l1*l2)+Fy;
   }
   /*m_point[p]->nx=n.x/w;///?????nx,ny,nz 有用吗
   m_point[p]->ny=n.y/w;
   m_point[p]->nz=n.z/w;*/
   m_point[p]->Fx=Fx/w;m_point[p]->Fy=Fy/w;
   m_index.RemoveAll();
}

void CDelaunayDoc::BaryCenter(CTriangle *temp)//MY WORK
{//三角形重心，用于向空间拓展
	ASSERT(temp!=NULL);
	int p1,p2,p3;
	p1=temp->m_p1;
	p2=temp->m_p2;
	p3=temp->m_p3;   
	CPointPos* p11=m_point[p1];
	CPointPos* p22=m_point[p2];
	CPointPos* p33=m_point[p3];
	temp->m_x=(p11->m_x+p22->m_x +p33->m_x )/double(3);
	temp->m_y=(p11->m_y+p22->m_y +p33->m_y )/double(3);
    //temp->m_z=(p11->m_z+p22->m_z +p33->m_z )/double(3);
}

double CDelaunayDoc::GetMold(CPointPos *p)//取模
{
	double m=0;
	m=p->m_x*p->m_x+p->m_y*p->m_y+p->m_z*p->m_z;
	m=sqrt(m);
	return m;

}

CPointPos* CDelaunayDoc::GetChuiZu(double x,double y, CPointPos* p2, CPointPos* p3)
{//求 p1 对于p2和p3 的垂足.
	CPointPos* h=new CPointPos();
	double k23;	
	if (p3->m_x-p2->m_x==0){
		h->m_x=p2->m_x;
		h->m_y=y;
	}
	else{
		k23=(p3->m_y-p2->m_y)/(p3->m_x-p2->m_x);
		if (k23==0){
			h->m_x=x;
			h->m_y=p2->m_y;
		}
		else{
            h->m_x=(k23*k23*p2->m_x+x+(y-p2->m_y)*k23)/(1.0+k23*k23);
            h->m_y=p2->m_y+k23*(h->m_x-p2->m_x);
		}
	}
	return h;
}

double CDelaunayDoc::D(CTriangle *temp, int i, int j)
{//沿边方向导:i to j,用i点处的两个偏导Fx,Fy 和i to j 的边的方向上的单位矢量作内积
	int p[4];p[0]=0;
	p[1]=temp->m_p1;p[2]=temp->m_p2;p[3]=temp->m_p3;
	POI vector;
	double d;
	vector.x=m_point[p[j]]->m_x-m_point[p[i]]->m_x;
	vector.y=m_point[p[j]]->m_y-m_point[p[i]]->m_y;
    vector=Unitization(vector);//要是单位化就错了
    d=DotProduction(vector.x,vector.y,0,m_point[p[i]]->Fx,m_point[p[i]]->Fy,0);
	return d;
}

double CDelaunayDoc::F(CTriangle *temp, int i, int j)
{//i to j,ij边上的法向导数
	int p[4];p[0]=0;
	p[1]=temp->m_p1;p[2]=temp->m_p2;p[3]=temp->m_p3;
	int h=TheOtherPoint(p[i],p[j],temp);
    CPointPos* H=GetChuiZu(m_point[h]->m_x,m_point[h]->m_y,m_point[p[i]],m_point[p[j]]);
	POI n;
	n.x=H->m_x-m_point[h]->m_x;
	n.y=H->m_y-m_point[h]->m_y;
	n.z=0;
	//n=Unitization(n);//要是单位化就错了
	double n1,n2;
	n1=DotProduction(m_point[p[j]]->Fx,m_point[p[j]]->Fy,0,n.x,n.y,0);
	n2=DotProduction(m_point[p[i]]->Fx,m_point[p[i]]->Fy,0,n.x,n.y,0);
	delete H;
	return (n1+n2)/double(2);
}

void CDelaunayDoc::DrawTri(int m_p1,int m_p2,CTriangle *tri)
{
	int i,j;
	POI p1,p2,o;	
	p1.x=m_point[m_p1]->m_x;
	p1.y=m_point[m_p1]->m_y;
	p2.x=m_point[m_p2]->m_x;
	p2.y=m_point[m_p2]->m_y;
	o.x=tri->m_x;o.y=tri->m_y;
	int n;
	n=8;
	double x[8][8];
	double y[8][8];

	for(j=0;j<n;j++)
	{
		for(i=0;i<n-j;i++)
		{
			x[j][i]=p1.x+j*(p2.x-p1.x)/double(n-1)+
				i*(o.x+j*(p2.x-o.x)/double(n-1)-p1.x-j*(p2.x-p1.x)/double(n-1))
				/double(n-1-j);
			y[j][i]=p1.y+j*(p2.y-p1.y)/double(n-1)+
				i*(o.y+j*(p2.y-o.y)/double(n-1)-p1.y-j*(p2.y-p1.y)/double(n-1))
				/double(n-1-j);
		}
	}
	x[7][0]=p2.x;
	y[7][0]=p2.y;

    for(j=0;j<n-1;j++)
	{
		for(i=0;i<n-j-1;i++)
		{
			//normal=GetTriNormal(m_hct[0],m_hct[1],m_hct[2]);
			//glNormal3d(normal.x,normal.y,normal.z);
			glBegin(GL_TRIANGLES);
			  glVertex3d(x[j][i],y[j][i],Bezier(x[j][i],y[j][i],m_p1,m_p2,tri));
			  glVertex3d(x[j+1][i],y[j+1][i],Bezier(x[j+1][i],y[j+1][i],m_p1,m_p2,tri));
			  glVertex3d(x[j][i+1],y[j][i+1],Bezier(x[j][i+1],y[j][i+1],m_p1,m_p2,tri));
			glEnd();		
		}
	}
	for(j=1;j<n-1;j++)
	{
		for(i=0;i<n-j-1;i++)
		{
			//normal=GetTriNormal(m_hct[0],m_hct[1],m_hct[2]);
			//glNormal3d(normal.x,normal.y,normal.z);
			glBegin(GL_TRIANGLES);
			  glVertex3d(x[j][i],y[j][i],Bezier(x[j][i],y[j][i],m_p1,m_p2,tri));
			  glVertex3d(x[j][i+1],y[j][i+1],Bezier(x[j][i+1],y[j][i+1],m_p1,m_p2,tri));
			  glVertex3d(x[j-1][i+1],y[j-1][i+1],Bezier(x[j-1][i+1],y[j-1][i+1],m_p1,m_p2,tri)); 
			glEnd();		
		}
	}
/*
	for(j=1;j<n;j++)
	{
		for(i=0;i<n-j;i++)
		{
			glBegin(GL_LINES);
			  glVertex3d(x[j][i],y[j][i],Bezier(x[j][i],y[j][i],m_p1,m_p2,tri));
			  glVertex3d(x[j-1][i],y[j-1][i],Bezier(x[j-1][i],y[j-1][i],m_p1,m_p2,tri));			
			  glVertex3d(x[j][i],y[j][i],Bezier(x[j][i],y[j][i],m_p1,m_p2,tri));
			  glVertex3d(x[j-1][i+1],y[j-1][i+1],Bezier(x[j-1][i+1],y[j-1][i+1],m_p1,m_p2,tri));
			glEnd();		
		}
	}

	for(j=0;j<n-1;j++)
	{
		for(i=0;i<n-j-1;i++)
		{
			glBegin(GL_LINES);
			  glVertex3d(x[j][i],y[j][i],Bezier(x[j][i],y[j][i],m_p1,m_p2,tri));	
			  glVertex3d(x[j][i+1],y[j][i+1],Bezier(x[j][i+1],y[j][i+1],m_p1,m_p2,tri));
			glEnd();		
		}
	}
*/

}

double CDelaunayDoc::Bezier(double x, double y,int m_p1,int m_p2,CTriangle* tri)//点的定位
{
	//s1,s2,s3 : 面积坐标
	double s1,s2,s3,s;
	double b[4][4][4];
	POI p1,p2,o,p;	
	p1.x=m_point[m_p1]->m_x;
	p1.y=m_point[m_p1]->m_y;
	p2.x=m_point[m_p2]->m_x;
	p2.y=m_point[m_p2]->m_y;
	o.x=tri->m_x;o.y=tri->m_y;
	p.x=x;p.y=y;

	s=S(o,p1,p2);
	s1=S(p,p1,p2)/s;
	s2=S(o,p,p2)/s;
	s3=S(o,p1,p)/s;
	if(fabs(s1)<0.000000000001) s1=0;
	if(fabs(s2)<0.000000000001) s2=0;
	if(fabs(s3)<0.000000000001) s3=0;
	int i,j,k;
	double B=0.0;
	b[0][0][3]=m_point[m_p2]->m_z;b[0][3][0]=m_point[m_p1]->m_z;b[3][0][0]=tri->o;
	if(m_p1==tri->m_p1)
	{
		b[1][1][1]=tri->e3;
		b[2][1][0]=tri->b1;b[1][2][0]=tri->c1;
		b[0][2][1]=tri->d31;b[0][1][2]=tri->d32;
		b[2][0][1]=tri->b2;b[1][0][2]=tri->c2;		
	
		for(i=0;i<4;i++){
			for(j=0;j<=(3-i);j++){
				for(k=0;k<=(3-i-j);k++){
					if((i+j+k)==3){
					B=B+b[i][j][k]*double(6)/(Factorial(i)*Factorial(j)*Factorial(k))*
						Power(s1,i)*Power(s2,j)*Power(s3,k);
					}
				}
			}
		}
	}
	if(m_p1==tri->m_p2)
	{
		b[1][1][1]=tri->e1;
		b[2][1][0]=tri->b2;b[1][2][0]=tri->c2;
		b[0][2][1]=tri->d12;b[0][1][2]=tri->d13;
		b[2][0][1]=tri->b3;b[1][0][2]=tri->c3;		
	
		for(i=0;i<4;i++){
			for(j=0;j<=(3-i);j++){
				for(k=0;k<=(3-i-j);k++){
					if((i+j+k)==3){
					B=B+b[i][j][k]*double(6)/(Factorial(i)*Factorial(j)*Factorial(k))*
						Power(s1,i)*Power(s2,j)*Power(s3,k);
					}
				}
			}
		}
	}
	if(m_p1==tri->m_p3)
	{
		b[1][1][1]=tri->e2;
		b[2][1][0]=tri->b3;b[1][2][0]=tri->c3;
		b[0][2][1]=tri->d23;b[0][1][2]=tri->d21;
		b[2][0][1]=tri->b1;b[1][0][2]=tri->c1;		
	
		for(i=0;i<4;i++){
			for(j=0;j<=(3-i);j++){
				for(k=0;k<=(3-i-j);k++){
					if((i+j+k)==3){
					B=B+b[i][j][k]*double(6)/(Factorial(i)*Factorial(j)*Factorial(k))*
						Power(s1,i)*Power(s2,j)*Power(s3,k);
					}
				}
			}
		}
	}
return B;

}

double CDelaunayDoc::S(POI p1, POI p2, POI p3)
{
	//求三角形面积,以右下角为原点时,s>0为逆时针(左转),
  	//s=0为三点重合
  	double s;
  	s=p1.x*p2.y+p2.x*p3.y+p1.y*p3.x-p2.y*p3.x-p1.y*p2.x-p1.x*p3.y;
  	s=s/2.0;
  	return s;
}

int CDelaunayDoc::Factorial(int n)
{//阶乘
	int x=n;
	if(n==0 ||n==1) 
		return 1;
	else{		
		for(int i=1;i<n;i++){
			x=x*(n-i);
		}
	}
    return x;
}

double CDelaunayDoc::Power(double a, int e)
{
	double x=1.0;
	if(e==0){
		return 1.0;
	}
	else{
		if(a==0.0){
			return 0.0;
		}
        for(int i=1;i<=e;i++){
			x=x*a;
		}
	}
    return x;
}

CTriangle* CDelaunayDoc::Belong(double x, double y)//判断点的位置
{
	CTriangle* tri;
	POI p1,p2,p3,p;
	p.x=x;p.y=y;
	POSITION pos = m_tri.GetHeadPosition();	
	while( pos != NULL )
	{
        tri=m_tri.GetAt( pos );
		double s1,s2,s3;
	    p1.x=m_point[tri->m_p1]->m_x;
		p1.y=m_point[tri->m_p1]->m_y;
		p2.x=m_point[tri->m_p2]->m_x;
		p2.y=m_point[tri->m_p2]->m_y;
	    p3.x=m_point[tri->m_p3]->m_x;
		p3.y=m_point[tri->m_p3]->m_y;
	    
		s1=S(p,p1,p2);
		s2=S(p,p2,p3);
		s3=S(p,p3,p1);
		if(s1>0 && s2>0 && s3>0)
			return tri;
	
		m_tri.GetNext( pos );
	}
    return NULL;
}

int CDelaunayDoc::Belong(double x, double y, CTriangle *tri)
{
	POI p1,p2,o,p;
	p.x=x;p.y=y;
	double s1,s2,s3;
	o.x=tri->m_x;
	o.y=tri->m_y;
	int a[4];
	a[0]=tri->m_p1;
	a[1]=tri->m_p2;
	a[2]=tri->m_p3;
	a[3]=tri->m_p1;
	for(int i=0;i<3;i++)
	{
		p1.x=m_point[a[i]]->m_x;
		p1.y=m_point[a[i]]->m_y;
		p2.x=m_point[a[i+1]]->m_x;
		p2.y=m_point[a[i+1]]->m_y;
		s1=S(p,o,p1);
		s2=S(p,p1,p2);
		s3=S(p,p2,o);
		if(s1>0 && s2>0 && s3>0)
		   return i;
	}
}

void CDelaunayDoc::Wang()//空间的网格化
{
	int n=50;
    CTriangle* tri;
	double x[60];
	double y[60];
	double z[60][60];
	int i,j,what;
	for(i=0;i<n;i++)
	{
		x[i]=double(i)/double(n);
	    y[i]=double(i)/double(n);
	}
		
	for(i=0;i<n;i++)
	{
		for(j=0;j<n;j++)
		{
			tri=Belong(x[i],y[j]);
			if(tri!=NULL)
			{
			   what=Belong(x[i],y[j],tri);
			   switch(what)
			   {
			   case 0:
				   z[i][j]=Bezier(x[i],y[j],tri->m_p1,tri->m_p2,tri);
				   break;
			   case 1:
				   z[i][j]=Bezier(x[i],y[j],tri->m_p2,tri->m_p3,tri);
				   break;
			   default:
				   z[i][j]=Bezier(x[i],y[j],tri->m_p3,tri->m_p1,tri);
				   break;
			   }
			   //z[i][j]=z[i][j]-0.2;
			}
			else
				z[i][j]=0;
		}
	}
	for(i=0;i<n;i++)
	{
		for(j=0;j<n-1;j++)
		{
		    if((z[i][j]!=0) && (z[i][j+1]!=0))
			{
			   glBegin(GL_LINES);
				 glVertex3d(x[i],y[j],z[i][j]);
				 glVertex3d(x[i],y[j+1],z[i][j+1]);
			   glEnd();
			}
		}
	}
	for(j=0;j<n;j++)
	{
		for(i=0;i<n-1;i++)
		{
			if((z[i][j]!=0) && (z[i+1][j]!=0))
			{
			   glBegin(GL_LINES);
			    glVertex3d(x[i],y[j],z[i][j]);
				glVertex3d(x[i+1],y[j],z[i+1][j]);
		       glEnd();
			}
			
		}
	}

}