www.pudn.com > VC++Delaunay.rar > DelaunayDoc.cpp
// DelaunayDoc.cpp : implementation of the CDelaunayDoc class
//
#include "stdafx.h"
#include "Delaunay.h"
#include "DelaunayDoc.h"
#include "pointview.h"
#ifdef _DEBUG
#define new DEBUG_NEW
#undef THIS_FILE
static char THIS_FILE[] = __FILE__;
#endif
/////////////////////////////////////////////////////////////////////////////
// CDelaunayDoc
IMPLEMENT_DYNCREATE(CDelaunayDoc, CDocument)
BEGIN_MESSAGE_MAP(CDelaunayDoc, CDocument)
//{{AFX_MSG_MAP(CDelaunayDoc)
ON_COMMAND(ID_BUTTON_ADD, OnButtonAdd)
ON_UPDATE_COMMAND_UI(ID_BUTTON_ADD, OnUpdateButtonAdd)
ON_COMMAND(ID_BUTTON_BB, OnButtonBB)
ON_UPDATE_COMMAND_UI(ID_BUTTON_BB, OnUpdateButtonBB)
//}}AFX_MSG_MAP
END_MESSAGE_MAP()
/////////////////////////////////////////////////////////////////////////////
// CDelaunayDoc construction/destruction
CDelaunayDoc::CDelaunayDoc()
{
// TODO: add one-time construction code here
}
CDelaunayDoc::~CDelaunayDoc()
{
}
BOOL CDelaunayDoc::OnNewDocument()
{
if (!CDocument::OnNewDocument())
return FALSE;
// TODO: add reinitialization code here
// (SDI documents will reuse this document)
//CPointPos* temp=new CPointPos;
//m_point.Add(temp);
return TRUE;
}
/////////////////////////////////////////////////////////////////////////////
// CDelaunayDoc serialization
void CDelaunayDoc::Serialize(CArchive& ar)
{
if (ar.IsStoring())
{
}
else
{
// TODO: add loading code here
}
m_con.Serialize(ar);
m_point.Serialize(ar);////////
m_tri.Serialize(ar);
// m_n.Serialize(ar);
}
/////////////////////////////////////////////////////////////////////////////
// CDelaunayDoc diagnostics
#ifdef _DEBUG
void CDelaunayDoc::AssertValid() const
{
CDocument::AssertValid();
}
void CDelaunayDoc::Dump(CDumpContext& dc) const
{
CDocument::Dump(dc);
}
#endif //_DEBUG
/////////////////////////////////////////////////////////////////////////////
// CDelaunayDoc commands
void CDelaunayDoc::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;
}
//p->m_z=sin((x-0.5)*(y-0.5));
//z=0.25-(x-0.5)*(x-0.5)-(y-0.5)*(y-0.5);
//if(z<0) z=0;
//p->m_z=sqrt(z)*double(4);//;
//p->m_z = (float)sin(3.0*x)*cos(3.0*y)/3.0f;
ASSERT(p!=NULL);
m_point.Add(p);
POSITION pos;
pos=GetFirstViewPosition();
m_plen=m_point.GetSize();
while(pos!=NULL)
{
CView *pp=GetNextView(pos);
if(pp->IsKindOf(RUNTIME_CLASS(CPointView)))
{
CListCtrl &list=((CPointView*)pp)->GetListCtrl();
LV_ITEM li;
CString str;
str.Format("%d",m_plen-1);
li.mask=LVIF_TEXT;
li.pszText=str.GetBuffer(10);;
li.iItem=m_plen-1;
li.iSubItem=0;
list.InsertItem(&li);
str.ReleaseBuffer();
str.Format("%f",p->m_x);
li.pszText=str.GetBuffer(10);
li.iItem=m_plen-1;
li.iSubItem=1;
list.SetItem(&li);
str.ReleaseBuffer();
str.Format("%f",p->m_y);
li.pszText=str.GetBuffer(10);
li.iItem=m_plen-1;
li.iSubItem=2;
list.SetItem(&li);
str.ReleaseBuffer();
}
}
};
////////////////////////////////////////////
void CDelaunayDoc::DeleteContents()
{
int i,max;
max=m_point.GetSize();
for(i=0;im_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 CDelaunayDoc::OnButtonAdd()
{
m_DoWhat=DO_ADD;
}
void CDelaunayDoc::OnUpdateButtonAdd(CCmdUI* pCmdUI)
{
if(m_DoWhat==DO_ADD)
pCmdUI->SetCheck(1);
else
pCmdUI->SetCheck(0);
}
void CDelaunayDoc::AddTriangle(int p)
{
//add new triangles
CTriangle* pTriangle;
int max=m_edge.GetSize();
for(int i=0;im_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 CDelaunayDoc::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;
}
int CDelaunayDoc::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 CDelaunayDoc::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;mmax-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(r0){
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 rl
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{//rm_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;
}
}
}
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;im_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)
{
//求三角形面积,以右下角为原点时,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;
}
void CDelaunayDoc::OnButtonBB()
{
m_DoWhat=DO_INTERPOLATION;
int max=m_point.GetSize();
for(int i=0;im_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 a;//save 三角片所在平面的方程ax+by+cz+d=0的a
CArray b;//save 三角片所在平面的方程ax+by+cz+d=0的b
CArray 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 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();////????
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)
{//三角形重心
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;jm_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;im_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;im_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