gendelaunaydoc.cpp

一个二维Delaunay三角网格剖分算法

// GenDelaunayDoc.cpp : implementation of the CGenDelaunayDoc class
//

#include "stdafx.h"
#include "GenDelaunay.h"

#include "GenDelaunayDoc.h"

#ifdef _DEBUG
#define new DEBUG_NEW
#undef THIS_FILE
static char THIS_FILE[] = __FILE__;
#endif





/////////////////////////////////////////////////////////////////////////////
// CGenDelaunayDoc

IMPLEMENT_DYNCREATE(CGenDelaunayDoc, CDocument)

BEGIN_MESSAGE_MAP(CGenDelaunayDoc, CDocument)
	//{{AFX_MSG_MAP(CGenDelaunayDoc)
		// NOTE - the ClassWizard will add and remove mapping macros here.
		//    DO NOT EDIT what you see in these blocks of generated code!
	//}}AFX_MSG_MAP
END_MESSAGE_MAP()

/////////////////////////////////////////////////////////////////////////////
// CGenDelaunayDoc construction/destruction

CGenDelaunayDoc::CGenDelaunayDoc()
{
	// TODO: add one-time construction code here

}

CGenDelaunayDoc::~CGenDelaunayDoc()
{
}

BOOL CGenDelaunayDoc::OnNewDocument()
{
	if (!CDocument::OnNewDocument())
		return FALSE;

	// TODO: add reinitialization code here
	// (SDI documents will reuse this document)

	return TRUE;
}



/////////////////////////////////////////////////////////////////////////////
// CGenDelaunayDoc serialization

void CGenDelaunayDoc::Serialize(CArchive& ar)
{
	if (ar.IsStoring())
	{
		// TODO: add storing code here
	}
	else
	{
		// TODO: add loading code here
	}
}

/////////////////////////////////////////////////////////////////////////////
// CGenDelaunayDoc diagnostics

#ifdef _DEBUG
void CGenDelaunayDoc::AssertValid() const
{
	CDocument::AssertValid();
}

void CGenDelaunayDoc::Dump(CDumpContext& dc) const
{
	CDocument::Dump(dc);
}
#endif //_DEBUG

/////////////////////////////////////////////////////////////////////////////
// CGenDelaunayDoc commands


void CGenDelaunayDoc::AddPoint(double x, double y)
{
	double z;
	int m_plen;//当前节点个数
	CPointPos *p;
	p=new CPointPos(x,y);
	if(x<=0.4)
	{
		z=0.03-(x-0.2)*(x-0.2)-(y-0.5)*(y-0.5);
		if(z<0) z=0.1;
		p->m_z=sqrt(z)*double(3);
	}
	else{
		if(x>=0.6)
		{
			z=0.05-(x-0.8)*(x-0.8)-(y-0.5)*(y-0.5);
			if(z<0) z=0.1;
			p->m_z=sqrt(z)*double(1.5);
		}
		else
            p->m_z=0.5;
	}
	ASSERT(p!=NULL);	
	m_point.Add(p);
	POSITION pos;
	pos=GetFirstViewPosition();
	m_plen=m_point.GetSize();
};


double CGenDelaunayDoc::S(int p1, int p2, int p3)
{
	//求三角形面积,以右下角为原点时,s>0为逆时针(左转),
	//s=0为三点重合
	double s;
	//POI m_p1,m_p2,m_p3;
	CPointPos *m_p1=new CPointPos(m_point[p1]->m_x,m_point[p1]->m_y);
    CPointPos *m_p2=new CPointPos(m_point[p2]->m_x,m_point[p2]->m_y);
    CPointPos *m_p3=new CPointPos(m_point[p3]->m_x,m_point[p3]->m_y);
	s=m_p1->m_x*m_p2->m_y+m_p2->m_x*m_p3->m_y+m_p1->m_y*
		m_p3->m_x-m_p2->m_y*m_p3->m_x-m_p1->m_y*m_p2->m_x-m_p1->m_x*m_p3->m_y;
	s=s/2;
	delete m_p1;
	delete m_p2;
	delete m_p3;
	
	return s;
}


void CGenDelaunayDoc::Center(CTriangle *temp)//外接圆心和半径
{
	ASSERT(temp!=NULL);
	int p1,p2,p3;
	p1=temp->m_p1;
	p2=temp->m_p2;
	p3=temp->m_p3;
    double k21,k31;
	CPointPos* p11=(CPointPos*)m_point[p1];
	CPointPos* p22=(CPointPos*)m_point[p2];
	CPointPos* p33=(CPointPos*)m_point[p3];
	if((p22->m_x-p11->m_x)==0) 
	{
		temp->m_yc=(p11->m_y+p22->m_y)/2;//circle 纵坐标
        k31=(p33->m_y-p11->m_y)/(p33->m_x-p11->m_x);
		if (k31==0)
		{
			temp->m_xc=(p11->m_x+p33->m_x)/2;
		}
		else
		{
			temp->m_xc=(p11->m_x+p33->m_x)/2+
				((p11->m_y+p33->m_y)/2-temp->m_yc)*k31;
		}
	}
	else 
	{
		k21=(p22->m_y-p11->m_y)/(p22->m_x-p11->m_x);
		if (k21==0)
		{
			temp->m_xc=(p11->m_x+p22->m_x)/2;
			if((p22->m_x-p33->m_x)==0) 
			{
				temp->m_yc=(p11->m_y+p33->m_y)/2;
			}
			else
			{
                k31=(p33->m_y-p11->m_y)/
					(p33->m_x-p11->m_x);
				temp->m_yc=((p33->m_x+p11->m_x)/2-
					temp->m_xc)/k31+
					(p11->m_y+p33->m_y)/2;
			}
		}
		else
		{
			if((p33->m_x-p11->m_x)==0)
			{
				temp->m_yc=(p33->m_y+p11->m_y)/2;
				temp->m_xc=k21*(p11->m_y/2+p22->m_y/2-temp->m_yc)+
					(p11->m_x+p22->m_x)/2;
			}
			else
			{
				k31=(p33->m_y-p11->m_y)/(p33->m_x-p11->m_x);
				if(k31==0)
				{
					temp->m_xc=(p11->m_x+p33->m_x)/2;
					temp->m_yc=((p11->m_x+p22->m_x)/2-
						temp->m_xc)/k21+
						(p11->m_y+p22->m_y)/2;
				}
				else
				{
					temp->m_xc=(k21*k31*(p22->m_y-p33->m_y)+
						(p11->m_x+p22->m_x)*k31-(p11->m_x+p33->m_x)*k21)/(k31-k21);
					temp->m_xc=temp->m_xc/2;
					temp->m_yc=((p11->m_x+p22->m_x)/2-
						temp->m_xc)/k21+
						(p11->m_y+p22->m_y)/2;
				}
			}
		}
	}
	
	temp->m_rad=sqrt((p33->m_x-temp->m_xc)*(p33->m_x-temp->m_xc)
		+ (p33->m_y-temp->m_yc)*(p33->m_y-temp->m_yc));
	
}
void CGenDelaunayDoc::BaryCenter(CTriangle *temp)
{//三角形重心
	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);
}


int CGenDelaunayDoc::GetInitEdges(double x, double y, int p)
{
	// GetInitEdges : ransack each trangle to record it's edges
	//when the piont belong it's circle,record the triangle's position in m_tri to m_index
	//'x' and 'y' are the coordinates of the insert point
	//'p' is the mark or th insert point in 'm_point'
	// the returned value 'k/2' is the number of cirecle the point belonged
	POSITION POS;
	CBorder *m_border;
	CBorder *m_border1;
	CTriangle* pTriangle;
	int j=0,i=0,k=0,max=0;
	CPointPos *point=new CPointPos(x,y);
    POS = m_tri.GetHeadPosition();
	while(POS != NULL ){
		i= m_tri.GetAt( POS )->Where(point);
		
		if(i==POS_ERROR){
			return POS_ERROR;
		}
		if(i==POS_ON)//in circle POS_ON=2
		{
			k=k+i;
			m_index.Add(POS);
			pTriangle=m_tri.GetAt(POS);
			//ransack three edge of a triangle,if an edge of a triangle
			//have belong to the array 'm_pDoc->m_edge' delete both of them to form 
			// the border of the inserted polygon which is stored in the array
			//'m_pDoc->m_edge'
			int array[4];
			array[0]=pTriangle->m_p1;
			array[1]=pTriangle->m_p2;
			array[2]=pTriangle->m_p3;
			array[3]=pTriangle->m_p1;
			for(int h=0;h<3;h++)
			{					  
				m_border=new CBorder(array[h],array[h+1]);
				max=m_edge.GetSize();
				if(max==0)
				{
					m_edge.Add(m_border);
				}
				else
				{
					j=0;
					for(int m=0;m<max;m++)
					{						  				         
						m_border1=m_edge[m];
						if(TwoEdgeSuperposition(m_border,m_border1)==1)//?//???
						{
							m_edge.RemoveAt(m,1);
							max=max-1;
							j++;
						}
					}
					if(j==0)
					{
						m_edge.Add(m_border);
					}				    
				}		 
			}/////for
		}/////////////if(POS_ON)
		pTriangle=m_tri.GetNext(POS);
	}//////////////while 
	//if(k==0)//the point don't belong any circle,
	//i.e.don't the convexity
	delete point;
	return k;
}


bool CGenDelaunayDoc::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;
}

int CGenDelaunayDoc::TwoEdgeSuperposition(CBorder *b1, CBorder *b2)
{
	if(b1==b2 || (b1->m_p1==b2->m_p2 && b1->m_p2==b2->m_p1) ) return 1;
	return 0;
}
int CGenDelaunayDoc::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];
	}
}
CPointPos* CGenDelaunayDoc::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) return NULL;
		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;
		}
	}       
}

void CGenDelaunayDoc::EditCon(int r, int l,int p)
{
	// r : save the value of i as 右支撑点 or 应去凸包边界节点右边界		  
	// l : save the value of i as 左支撑点 or 应去凸包边界节点左边界
	//'p' is the mark or th insert point in 'm_point'
	CBorder *m_border;
	int hr=0;//record the number of points on edge of convexity that should be deleted before r(include r) 
	int hl=0;//record the number of points on edge of convexity that should be deleted after l(include l)
	int i=r;
	int i1;
	int h=hr+hl;
	int max=m_con.GetSize();
	//search rightward, first
	if((i+1)>max-1)  i1=i+1-max;	 
	else  	          i1=i+1;	
	while(S(p,m_con[i],m_con[i1])<0){                  
		m_border=new CBorder(m_con[i],m_con[i1]);
		m_edge.Add(m_border);	
		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             	
}		


void CGenDelaunayDoc::DelTriMarked()
{
	int max=m_index.GetSize();
	for(int i=0;i<max;i++)
	{
		m_tri.RemoveAt(m_index[i]);
	}
}


void CGenDelaunayDoc::AddTriangle(int p)
{
	//add new triangles
	CTriangle* pTriangle;
	int max=m_edge.GetSize();	
	for(int i=0;i<max;i++)
	{
		//Every triangle is stored anticlockwise,too.
		if(S(m_edge[i]->m_p1,m_edge[i]->m_p2,p)>0){
			pTriangle=new CTriangle(m_edge[i]->m_p1,m_edge[i]->m_p2,p);
		}
		else{
			pTriangle=new CTriangle(m_edge[i]->m_p2,m_edge[i]->m_p1,p);
		}
		Center(pTriangle);
		BaryCenter(pTriangle);
		m_tri.AddTail(pTriangle);						
	}
}


double CGenDelaunayDoc::S(CPointPos *p1,CPointPos *p2,CPointPos *p3)
{
	//求三角形面积,以右下角为原点时,s>0为逆时针(左转),
	//s=0为三点重合
	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;
}