delaunaydoc.cpp

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

			r=i;
            i++;
			hr++;
			if(i+1>max-1)	i1=i+1-max;			
			else	        i1=i+1;
			if(i>max-1)     i=i-max;
	}	
	i=l;
    if((i-1)<0) i1=i-1+max;
	else        i1=i-1;
	while(S(m_con[i1],m_con[i],p)<0)
	{
		    m_border=new CBorder(m_con[i],m_con[i1]);
			m_edge.Add(m_border);
			l=i;
			i--;
			hl++;
			if(i-1<0) i1=i-1+max;
			else      i1=i-1;
			if(i<0)   i=i+max;		
	}
	h=hr+hl;
	if(h>0)
	{
	       if(r<l){
	    		  if(hl>0){
						  if(hr>0){//both>0
							      m_con.RemoveAt(l,max-l);
				                  m_con.RemoveAt(0,hr);
						  }                            						
						  else{//hl>0,hr=0
							      m_con.RemoveAt(l,max-l);
							      m_con.RemoveAt(0,hl-max+l);
						  }
				          m_con.InsertAt(0,p);                           
				  }//end if hl>0
				  else{//hl=0,hr>0							        
					   m_con.RemoveAt(0,hr);							
				       m_con.InsertAt(0,p);						
				  }
			}//end if r<l
			else{//r>l
				if(hl>0){
					    m_con.RemoveAt(l,h);
				        m_con.InsertAt(l,p);                           
				}//end if hl>0
				else{//hl=0,hr>0							        
					 m_con.RemoveAt(l+1,hr);							
				     m_con.InsertAt(l+1,p);						
				}
			}//end if r>l
	}//end if h>0
	else
	{
		if((r-l)>0){
		     if((r-l)==1) m_con.InsertAt(r,p);
		     else{
				 m_con.RemoveAt(l+1,r-l-1);							
				 m_con.InsertAt(l+1,p);                  
			 }
		}
		else{//r<l,i.e. r and l 在 m_con[0] 两侧
			if(l<(max-1)) m_con.RemoveAt(l+1,max-l-1);
            m_con.RemoveAt(0,r);							
			m_con.InsertAt(0,p);       
		}
	}//h<0             	
}								
			

CPointPos* CDelaunayDoc::IntersectionPoint(CPointPos *point1, CPointPos *point2,CPointPos *point3, CPointPos *point4)//四点形成的凸包进行判断和优化
{
	ASSERT(point1!=NULL);
    ASSERT(point2!=NULL);
	ASSERT(point3!=NULL);
    ASSERT(point4!=NULL);
    CPointPos *pos=new CPointPos(0,0);
	double k21,k43;
	if((point2->m_x-point1->m_x)==0 || (point4->m_x-point3->m_x)==0)
	{
		if((point2->m_x-point1->m_x)==0 && (point4->m_x-point3->m_x)==0) 
		if((point2->m_x-point1->m_x)==0)
		{
			k43=(point4->m_y-point3->m_y)/(point4->m_x-point3->m_x);
            pos->m_x=point1->m_x;// x of the intersection point
            pos->m_y=k43*(pos->m_x-point3->m_x)+point3->m_y;// y of the intersection point
		    return pos;
		}
		if((point4->m_x-point3->m_x)==0){
			k21=(point2->m_y-point1->m_y)/(point2->m_x-point1->m_x);
            pos->m_x=point3->m_x;// x of the intersection point
            pos->m_y=k21*(pos->m_x-point1->m_x)+point1->m_y;// y of the intersection point
		    return pos;
		}
	}
	else
	{
		k43=(point4->m_y-point3->m_y)/(point4->m_x-point3->m_x);
        k21=(point2->m_y-point1->m_y)/(point2->m_x-point1->m_x);
		if(fabs(k43-k21)<0.0000000001) 
			return NULL;
		else{		
			pos->m_x=(point3->m_y-point1->m_y+k21*point1->m_x-k43*point3->m_x)/(k21-k43);
            pos->m_y=k43*(pos->m_x-point3->m_x)+point3->m_y;
			return pos;
		}
	}       
}

int CDelaunayDoc::TheOtherPoint(int p1, int p2, CTriangle *temp)
{//p1,p2 are the mark recorded by m_tri in m_point
	int array[3];
	double s;
	array[0]=temp->m_p1;
	array[1]=temp->m_p2;
	array[2]=temp->m_p3;
	for(int i=0;i<3;i++)
	{
		s=S(array[i],p1,p2);
		if(fabs(s)>0.00000000001) return array[i];
	}
}

bool CDelaunayDoc::DelEdgeOrNot(int p1, int p2,int p)//是否删除边
{
   CTriangle* pTriangle;
   //p1=m_pDoc->m_edge[j]->m_p1,record the point's mark in m_point
   //p2=m_pDoc->m_edge[j]->m_p2,record the point's mark in m_point
   int i,j;
   for(i=0;i<m_index.GetSize();i++)
   {
	   pTriangle=m_tri.GetAt(m_index[i]);
	   int array[3];
	   array[0]=pTriangle->m_p1;
	   array[1]=pTriangle->m_p2;
	   array[2]=pTriangle->m_p3;	  
	   int a=-1,b=-1;
	   for(j=0;j<3;j++)
	   {
		   if(p1==array[j]) a=j;
	   }
	    for(j=0;j<3;j++)
	   {
		   if(p2==array[j]) b=j;
	   }
       if((a>=0) & (b>=0))
	   {
		   int other_point=TheOtherPoint(p1,p2,pTriangle);
		   CPointPos* pPoint=IntersectionPoint(m_point[p1],m_point[p2],
		                      m_point[p],m_point[other_point]);
		   if(pTriangle->Where(pPoint)==POS_ON){//整数截断使交点在边外
			   return true;
		   }
	   }	
   }
   return false;
}

 double CDelaunayDoc::S(CPointPos *p1,CPointPos *p2,CPointPos *p3)//求面积
 {
 	
  	double s;
  	s=p1->m_x*p2->m_y+p2->m_x*p3->m_y+p1->m_y*p3->m_x-
		p2->m_y*p3->m_x-p1->m_y*p2->m_x-p1->m_x*p3->m_y;
  	s=s/2.0;
  	return s;
 }

void CDelaunayDoc::OnButtonBB() 
{
	m_DoWhat=DO_INTERPOLATION;
	int max=m_point.GetSize();
	for(int i=0;i<max;i++)
	{
		FindRelativeTri(i);//get Fx,Fy,N of a point in m_point
	}
	POSITION POS;CTriangle* temp;
	POS = m_tri.GetHeadPosition();
	temp=m_tri.GetAt(POS);
	while(POS != NULL ){
		temp=m_tri.GetNext(POS);
		double d;
		int p[4];p[0]=0;
	    p[1]=temp->m_p1;p[2]=temp->m_p2;p[3]=temp->m_p3;
//////////////////////////////////////////////////////////////////////////
		//Get 6个 dij 计算无错 但沿边接的方向倒数的估算可能有问题,
		//1)可能不应单位化.2)?
		d=D(temp,2,3);
		temp->d12=m_point[p[2]]->m_z+d/double(3);

		d=D(temp,3,2);
		temp->d13=m_point[p[3]]->m_z+d/double(3);
        
		d=D(temp,1,3);
		temp->d21=m_point[p[1]]->m_z+d/double(3);
        
		d=D(temp,3,1);
		temp->d23=m_point[p[3]]->m_z+d/double(3);

		d=D(temp,1,2);
		temp->d31=m_point[p[1]]->m_z+d/double(3);

        d=D(temp,2,1);
		temp->d32=m_point[p[2]]->m_z+d/double(3);
//////////////////////////////////////////////////////////////////////////
		//计算无错
        temp->c1=(temp->d21+temp->d31+m_point[p[1]]->m_z)/double(3);
		temp->c2=(temp->d12+temp->d32+m_point[p[2]]->m_z)/double(3);
		temp->c3=(temp->d13+temp->d23+m_point[p[3]]->m_z)/double(3);
//////////////////////////////////////////////////////////////////////////
		/*Get 每边中点处f的法向导数值(XY平面上的) 
		在这边的两顶点间对顶点的法向导数做插值
		 顶点的法向导数由该点的Fx,Fy取得,*/	    
		temp->n1=F(temp,3,2);//计算无错
		temp->n2=F(temp,1,3);
		temp->n3=F(temp,2,1);
		
//////////////////////////////////////////////////////////////////////////
	/*Get bi and b0(o), 220页公式(5.48) 应再推几遍.HCT 和 三角形面积 都用的是
		XY平面上的?????3维的又如何????*///计算无错
		CPointPos* h1;
		CPointPos* h2;
		CPointPos* h3;
		//三角形o,f2,f3
        h1=GetChuiZu(temp->m_x,temp->m_y,m_point[p[2]],m_point[p[3]]);
		//三角形o,f3,f1
		h2=GetChuiZu(temp->m_x,temp->m_y,m_point[p[3]],m_point[p[1]]);
		//三角形o,f1,f2
		h3=GetChuiZu(temp->m_x,temp->m_y,m_point[p[1]],m_point[p[2]]);
		CPointPos* o=new CPointPos();//重心
		o->m_x=temp->m_x;
		o->m_y=temp->m_y;
		o->m_z=0;
		double s1=S(o,m_point[p[2]],m_point[p[3]]);
		double s2=S(o,m_point[p[3]],m_point[p[1]]);
		double s3=S(o,m_point[p[1]],m_point[p[2]]);
		temp->b1=S(o,m_point[p[3]],h2)/s2
        *(m_point[p[1]]->m_z+temp->d21-temp->d23-m_point[p[3]]->m_z)+
        S(o,m_point[p[1]],h3)/s3
        *(m_point[p[2]]->m_z+temp->d32-temp->d31-m_point[p[1]]->m_z)-
		(temp->n3+temp->n2)/double(0.75)-temp->c3-temp->c2+
		m_point[p[1]]->m_z+m_point[p[3]]->m_z+temp->d21+temp->d32+(temp->d31+temp->d23)*
			double(2);
		temp->b1=temp->b1/double(6);

temp->b2=S(o,m_point[p[2]],h1)/s1
        *(m_point[p[3]]->m_z+temp->d13-temp->d12-m_point[p[2]]->m_z)+
        S(o,m_point[p[1]],h3)/s3
        *(m_point[p[2]]->m_z+temp->d32-temp->d31-m_point[p[1]]->m_z)-
		(temp->n3+temp->n1)/double(0.75)-temp->c3-temp->c1+
		m_point[p[1]]->m_z+m_point[p[2]]->m_z+temp->d32+temp->d13+(temp->d31+temp->d12)*
		double(2);
temp->b2=temp->b2/double(6);

temp->b3=S(o,m_point[p[2]],h1)/s1
        *(m_point[p[3]]->m_z+temp->d13-temp->d12-m_point[p[2]]->m_z)+
        S(o,m_point[p[3]],h2)/s2
        *(m_point[p[1]]->m_z+temp->d21-temp->d23-m_point[p[3]]->m_z)-
		(temp->n1+temp->n2)/double(0.75)-temp->c1-temp->c2+
		m_point[p[2]]->m_z+m_point[p[3]]->m_z+temp->d13+temp->d21+(temp->d12+temp->d23)*
		double(2);
temp->b3=temp->b3/double(6);
delete o;
temp->o=(temp->b1+temp->b2+temp->b3)/double(3);
//////////////////////////////////////////////////////////////////////////
       //Get ei//计算无错
       temp->e1=double(3)*(-temp->b1+temp->b2+temp->b3)-(-temp->c1+temp->c2+temp->c3);
	   temp->e1=temp->e1/double(2);
	   temp->e2=double(3)*(temp->b1-temp->b2+temp->b3)-(temp->c1-temp->c2+temp->c3);
	   temp->e2=temp->e2/double(2);
	   temp->e3=double(3)*(temp->b1+temp->b2-temp->b3)-(temp->c1+temp->c2-temp->c3);
	   temp->e3=temp->e3/double(2);
//////////////////////////////////////////////////////////////////////////
	   //NOW we have got the 19 control points in a triangle
//       temp=m_tri.GetNext(POS);

	   delete h1;
	   delete h2;
	   delete h3;
	}   
}

void CDelaunayDoc::OnUpdateButtonBB(CCmdUI* pCmdUI) 
{
	if(m_DoWhat==DO_INTERPOLATION)
		pCmdUI->SetCheck(1);
	else 
		pCmdUI->SetCheck(0);	
}

void CDelaunayDoc::FindRelativeTri(int p)
{
	CPointPos *point;
    point=m_point[p];
	POSITION POS;
    POS = m_tri.GetHeadPosition();
	while(POS != NULL ){
		  int array[3];
		  array[0]= m_tri.GetAt(POS)->m_p1;
		  array[1]= m_tri.GetAt(POS)->m_p2;
		  array[2]= m_tri.GetAt(POS)->m_p3;	
          for(int h=0;h<3;h++)
		  {					  			 
			 if(point==m_point[array[h]])
			 {
                m_index.Add(POS);
			 }
		  }/////for
	      m_tri.GetNext(POS);
	}//////////////while 
	Get_Fx_Fy_N(p);
}

POI CDelaunayDoc::GetTriNormal(CTriangle *temp)
{//得到三角片的单位法向
	POI a,b,n;
	int p1,p2,p3;
    p1=temp->m_p1;
    p2=temp->m_p2;
    p3=temp->m_p3;
	a.x=m_point[p2]->m_x-m_point[p1]->m_x;
	a.y=m_point[p2]->m_y-m_point[p1]->m_y;
	a.z=m_point[p2]->m_z-m_point[p1]->m_z;
	b.x=m_point[p3]->m_x-m_point[p1]->m_x;
	b.y=m_point[p3]->m_y-m_point[p1]->m_y;
	b.z=m_point[p3]->m_z-m_point[p1]->m_z;
    //get the normal n=a*b(外积)
	n=VectorProduct(a.x,a.y,a.z,b.x,b.y,b.z);
	//单位化 normal
	n=Unitization(n);
    return n;
}

double CDelaunayDoc::GetDistance(CPointPos* p1, CPointPos* p2)
{
	double distance;
	distance=(p1->m_x-p2->m_x)*(p1->m_x-p2->m_x)+
		     (p1->m_y-p2->m_y)*(p1->m_y-p2->m_y)+
             (p1->m_z-p2->m_z)*(p1->m_z-p2->m_z);
	distance=sqrt(distance);
	return distance;
}

POI CDelaunayDoc::VectorProduct(double x1,double y1,double z1,double x2,double y2,double z2)
{//get p1*p2(外积) 对!!!!
	POI product;
	product.x=y1*z2-z1*y2;
	product.y=z1*x2-x1*z2;
	product.z=x1*y2-y1*x2;
	return product;
}

POI CDelaunayDoc::Unitization(POI p)
{	//单位化 对!!!!
	double m;
	m=p.x*p.x+p.y*p.y+p.z*p.z;
	m=sqrt(m);
	p.x=p.x/m;
	p.y=p.y/m;
	p.z=p.z/m;
    return p;
}

double CDelaunayDoc::DotProduction(double x1,double y1,double z1,double x2,double y2,double z2)
{//(x1,y1,z1)(内积)(x2,y2,z2)
	double product;
    product=x1*x2+y1*y2+z1*z2;
    return product;
}

void CDelaunayDoc::Get_Fx_Fy_N(int p)
{/*本处计算无错,但型值点法向没用,因算得是Fx and Fy,
	然后用 Fx and Fy 在 xy 平面上沿边的方向计算方向倒数
	没用型值点法向做沿空间三角形边的方向倒数:
	      D(i,i+1)=(Pi+1-Pi)-[(Pi+1-Pi)Ni]Ni
		  D(i,i+2)=(Pi+2-Pi)-[(Pi+2-Pi)Ni]Ni
	*/
   CPointPos *point;
   point=m_point[p];
   int max=m_index.GetSize();
   double Fx=0;//点的x偏导
   double Fy=0;
  // POI n;//点的法向量
  //n.x=0;n.y=0;n.z=0;
   double w=0.0;//Little法
   for(int i=0;i<max;i++)
   {
	  CArray<double,double> a;//save 三角片所在平面的方程ax+by+cz+d=0的a 
	  CArray<double,double> b;//save 三角片所在平面的方程ax+by+cz+d=0的b
	  CArray<double,double> c;//save 三角片所在平面的方程ax+by+cz+d=0的c
	  int array[3];
	  CTriangle* Tri;
	  Tri=m_tri.GetAt(m_index[i]);
	  array[0]= Tri->m_p1;
	  array[1]= Tri->m_p2;
	  array[2]= Tri->m_p3;	
	  //a,b,c 计算无错 对!!!!
      a.SetAtGrow(i,(m_point[array[1]]->m_y-m_point[array[0]]->m_y)*
		   (m_point[array[2]]->m_z-m_point[array[0]]->m_z)-
		   (m_point[array[1]]->m_z-m_point[array[0]]->m_z)*
		   (m_point[array[2]]->m_y-m_point[array[0]]->m_y));
	  b.SetAtGrow(i,(m_point[array[1]]->m_z-m_point[array[0]]->m_z)*
		   (m_point[array[2]]->m_x-m_point[array[0]]->m_x)-
		   (m_point[array[1]]->m_x-m_point[array[0]]->m_x)*
		   (m_point[array[2]]->m_z-m_point[array[0]]->m_z));
	  c.SetAtGrow(i,(m_point[array[1]]->m_x-m_point[array[0]]->m_x)*
		   (m_point[array[2]]->m_y-m_point[array[0]]->m_y)-
		   (m_point[array[1]]->m_y-m_point[array[0]]->m_y)*
		   (m_point[array[2]]->m_x-m_point[array[0]]->m_x));
	  
	  CArray<CPointPos*,CPointPos*> P;
	  double l1,l2;	  
      for(int h=0;h<3;h++)
	  {					  			 
		 if(point!=m_point[array[h]])
		 {
            P.Add(m_point[array[h]]);
		 }
	  }/////for
	  l1=GetDistance(point,P[0]);
	  l2=GetDistance(point,P[1]);
	  P.RemoveAll();////????