开通博客赚积分
发布资源赚积分
分类

源码开发语言/平台
当前位置：首页 > 源码/资料 > 图形图像 > GIS编程 > 查看源码
DelaunayDoc.cpp
资源名称：VC++Delaunay.zip　[点击查看]
上传用户：azhong891
上传日期：2013-06-04
资源大小：197k
文件大小：32k
源码类别：

GIS编程

开发平台：

Visual C++

DelaunayDoc.cpp：源码内容
  1. // DelaunayDoc.cpp : implementation of the CDelaunayDoc class
  2. //
  4. #include "stdafx.h"
  5. #include "Delaunay.h"
  7. #include "DelaunayDoc.h"
  8. #include "pointview.h"
  10. #ifdef _DEBUG
  11. #define new DEBUG_NEW
  12. #undef THIS_FILE
  13. static char THIS_FILE[] = __FILE__;
  14. #endif
  16. /////////////////////////////////////////////////////////////////////////////
  17. // CDelaunayDoc
  19. IMPLEMENT_DYNCREATE(CDelaunayDoc, CDocument)
  21. BEGIN_MESSAGE_MAP(CDelaunayDoc, CDocument)
  22. //{{AFX_MSG_MAP(CDelaunayDoc)
  23. ON_COMMAND(ID_BUTTON_ADD, OnButtonAdd)
  24. ON_UPDATE_COMMAND_UI(ID_BUTTON_ADD, OnUpdateButtonAdd)
  25. ON_COMMAND(ID_BUTTON_BB, OnButtonBB)
  26. ON_UPDATE_COMMAND_UI(ID_BUTTON_BB, OnUpdateButtonBB)
  27. //}}AFX_MSG_MAP
  28. END_MESSAGE_MAP()
  30. /////////////////////////////////////////////////////////////////////////////
  31. // CDelaunayDoc construction/destruction
  33. CDelaunayDoc::CDelaunayDoc()
  34. {
  35. // TODO: add one-time construction code here
  38. }
  40. CDelaunayDoc::~CDelaunayDoc()
  41. {
  42. }
  44. BOOL CDelaunayDoc::OnNewDocument()
  45. {
  46. if (!CDocument::OnNewDocument())
  47. return FALSE;
  49. // TODO: add reinitialization code here
  50. // (SDI documents will reuse this document)
  51. //CPointPos* temp=new CPointPos;
  52. //m_point.Add(temp);
  53. return TRUE;
  54. }
  58. /////////////////////////////////////////////////////////////////////////////
  59. // CDelaunayDoc serialization
  61. void CDelaunayDoc::Serialize(CArchive& ar)
  62. {
  63. if (ar.IsStoring())
  64. {
  66. }
  67. else
  68. {
  69. // TODO: add loading code here
  70. }
  71. m_con.Serialize(ar);
  72. m_point.Serialize(ar);////////
  73. m_tri.Serialize(ar);
  74. // m_n.Serialize(ar);
  75. }
  77. /////////////////////////////////////////////////////////////////////////////
  78. // CDelaunayDoc diagnostics
  80. #ifdef _DEBUG
  81. void CDelaunayDoc::AssertValid() const
  82. {
  83. CDocument::AssertValid();
  84. }
  86. void CDelaunayDoc::Dump(CDumpContext& dc) const
  87. {
  88. CDocument::Dump(dc);
  89. }
  90. #endif //_DEBUG
  92. /////////////////////////////////////////////////////////////////////////////
  93. // CDelaunayDoc commands
  95. void CDelaunayDoc::AddPoint(double x, double y)
  96. {
  97. double z;
  98. int m_plen;//当前节点个数
  99. CPointPos *p;
  100. p=new CPointPos(x,y);
  101. if(x<=0.4)
  102. {
  103. z=0.03-(x-0.2)*(x-0.2)-(y-0.5)*(y-0.5);
  104.     if(z<0) z=0.1;
  105. p->m_z=sqrt(z)*double(3);
  106. }
  107. else{
  108.     if(x>=0.6)
  109. {
  110. z=0.05-(x-0.8)*(x-0.8)-(y-0.5)*(y-0.5);
  111.         if(z<0) z=0.1;
  112.     p->m_z=sqrt(z)*double(1.5);
  113. }
  114. else
  115.             p->m_z=0.5;
  116. }
  117. //p->m_z=sin((x-0.5)*(y-0.5));
  118. //z=0.25-(x-0.5)*(x-0.5)-(y-0.5)*(y-0.5);
  119. //if(z<0) z=0;
  120. //p->m_z=sqrt(z)*double(4);//;
  121. //p->m_z = (float)sin(3.0*x)*cos(3.0*y)/3.0f;
  122. ASSERT(p!=NULL);
  123. m_point.Add(p);
  124. POSITION pos;
  125. pos=GetFirstViewPosition();
  126. m_plen=m_point.GetSize();
  127. while(pos!=NULL)
  128. {
  129. CView *pp=GetNextView(pos);
  130. if(pp->IsKindOf(RUNTIME_CLASS(CPointView)))
  131. {
  132. CListCtrl &list=((CPointView*)pp)->GetListCtrl();
  133. LV_ITEM li;
  134. CString str;
  135. str.Format("%d",m_plen-1);
  136. li.mask=LVIF_TEXT;
  137. li.pszText=str.GetBuffer(10);;
  138. li.iItem=m_plen-1;
  139. li.iSubItem=0;
  140. list.InsertItem(&li);
  141. str.ReleaseBuffer();
  142. str.Format("%f",p->m_x);
  143. li.pszText=str.GetBuffer(10);
  144. li.iItem=m_plen-1;
  145. li.iSubItem=1;
  146. list.SetItem(&li);
  147. str.ReleaseBuffer();
  148. str.Format("%f",p->m_y);
  149. li.pszText=str.GetBuffer(10);
  150. li.iItem=m_plen-1;
  151. li.iSubItem=2;
  152. list.SetItem(&li);
  153. str.ReleaseBuffer();
  154. }
  155. }
  156. };
  157. ////////////////////////////////////////////
  158. void CDelaunayDoc::DeleteContents() 
  159. {
  160. int i,max;
  161. max=m_point.GetSize();
  162. for(i=0;i<max;i++)
  163. {
  164. if(m_point[i]!=NULL)
  165. delete m_point[i];
  166. //if(m_n[i]!=NULL)//????????
  167. // delete m_n[i];
  168. }
  169. m_point.RemoveAll();
  170. // m_n.RemoveAll();
  172. POSITION pos = m_tri.GetHeadPosition();
  173. while( pos != NULL )
  174. {
  175. delete m_tri.GetNext( pos );
  176. }
  177. m_tri.RemoveAll();
  178.     m_con.RemoveAll();
  179. m_index.RemoveAll();
  181. CDocument::DeleteContents();
  182. }
  183. ///////////////////////////////////////////////////
  184. void CDelaunayDoc::Center(CTriangle *temp)//外接圆心和半径
  185. {
  186. ASSERT(temp!=NULL);
  187. int p1,p2,p3;
  188. p1=temp->m_p1;
  189. p2=temp->m_p2;
  190. p3=temp->m_p3;
  191.     double k21,k31;
  192. CPointPos* p11=(CPointPos*)m_point[p1];
  193. CPointPos* p22=(CPointPos*)m_point[p2];
  194. CPointPos* p33=(CPointPos*)m_point[p3];
  195. if((p22->m_x-p11->m_x)==0) 
  196. {
  197. temp->m_yc=(p11->m_y+p22->m_y)/2;//circle 纵坐标
  198.         k31=(p33->m_y-p11->m_y)/(p33->m_x-p11->m_x);
  199. if (k31==0)
  200. {
  201. temp->m_xc=(p11->m_x+p33->m_x)/2;
  202. }
  203. else
  204. {
  205.    temp->m_xc=(p11->m_x+p33->m_x)/2+
  206.    ((p11->m_y+p33->m_y)/2-temp->m_yc)*k31;
  207. }
  208. }
  209. else 
  210. {
  211. k21=(p22->m_y-p11->m_y)/(p22->m_x-p11->m_x);
  212. if (k21==0)
  213. {
  214. temp->m_xc=(p11->m_x+p22->m_x)/2;
  215. if((p22->m_x-p33->m_x)==0) 
  216. {
  217. temp->m_yc=(p11->m_y+p33->m_y)/2;
  218. }
  219.     else
  220. {
  221.                 k31=(p33->m_y-p11->m_y)/
  222.            (p33->m_x-p11->m_x);
  223.         temp->m_yc=((p33->m_x+p11->m_x)/2-
  224.          temp->m_xc)/k31+
  225.  (p11->m_y+p33->m_y)/2;
  226. }
  227. }
  228. else
  229. {
  230.     if((p33->m_x-p11->m_x)==0)
  231. {
  232. temp->m_yc=(p33->m_y+p11->m_y)/2;
  233. temp->m_xc=k21*(p11->m_y/2+p22->m_y/2-temp->m_yc)+
  234. (p11->m_x+p22->m_x)/2;
  235. }
  236. else
  237. {
  238.  k31=(p33->m_y-p11->m_y)/(p33->m_x-p11->m_x);
  239.  if(k31==0)
  240.  {
  241.                       temp->m_xc=(p11->m_x+p33->m_x)/2;
  242.       temp->m_yc=((p11->m_x+p22->m_x)/2-
  243.          temp->m_xc)/k21+
  244.  (p11->m_y+p22->m_y)/2;
  245.  }
  246.  else
  247.  {
  248.                           temp->m_xc=(k21*k31*(p22->m_y-p33->m_y)+
  249. (p11->m_x+p22->m_x)*k31-(p11->m_x+p33->m_x)*k21)/(k31-k21);
  250.   temp->m_xc=temp->m_xc/2;
  251.                           temp->m_yc=((p11->m_x+p22->m_x)/2-
  252.          temp->m_xc)/k21+
  253.  (p11->m_y+p22->m_y)/2;
  254.  }
  255. }
  256. }
  257. }
  259. temp->m_rad=sqrt((p33->m_x-temp->m_xc)*(p33->m_x-temp->m_xc)
  260.               + (p33->m_y-temp->m_yc)*(p33->m_y-temp->m_yc));
  262. }
  263. ///////////////////////////////////////////////////
  264. void CDelaunayDoc::OnButtonAdd() 
  265. {
  266. m_DoWhat=DO_ADD;
  267. }
  269. void CDelaunayDoc::OnUpdateButtonAdd(CCmdUI* pCmdUI) 
  270. {
  271. if(m_DoWhat==DO_ADD)
  272. pCmdUI->SetCheck(1);
  273. else 
  274. pCmdUI->SetCheck(0);
  275. }
  279. void CDelaunayDoc::AddTriangle(int p)
  280. {
  281.    //add new triangles
  282. CTriangle* pTriangle;
  283. int max=m_edge.GetSize();
  284. for(int i=0;i<max;i++)
  285. {
  286.    //Every triangle is stored anticlockwise,too.
  287.    if(S(m_edge[i]->m_p1,m_edge[i]->m_p2,p)>0){
  288. pTriangle=new CTriangle(m_edge[i]->m_p1,m_edge[i]->m_p2,p);
  289.    }
  290.    else{
  291.               pTriangle=new CTriangle(m_edge[i]->m_p2,m_edge[i]->m_p1,p);
  292.    }
  293.    Center(pTriangle);
  294.        BaryCenter(pTriangle);
  295.    m_tri.AddTail(pTriangle);
  296. }
  297. }
  301. double CDelaunayDoc::S(int p1, int p2, int p3)
  302. {
  303. //求三角形面积,以右下角为原点时,s>0为逆时针(左转),
  304. //s=0为三点重合
  305. double s;
  306. //POI m_p1,m_p2,m_p3;
  307. CPointPos *m_p1=new CPointPos(m_point[p1]->m_x,m_point[p1]->m_y);
  308.     CPointPos *m_p2=new CPointPos(m_point[p2]->m_x,m_point[p2]->m_y);
  309.     CPointPos *m_p3=new CPointPos(m_point[p3]->m_x,m_point[p3]->m_y);
  310. s=m_p1->m_x*m_p2->m_y+m_p2->m_x*m_p3->m_y+m_p1->m_y*
  311. 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;
  312. s=s/2;
  313. delete m_p1;
  314. delete m_p2;
  315. delete m_p3;
  317. return s;
  318. }
  320. int CDelaunayDoc::TwoEdgeSuperposition(CBorder *b1, CBorder *b2)
  321. {
  322. if(b1==b2 || (b1->m_p1==b2->m_p2 && b1->m_p2==b2->m_p1) ) return 1;
  323. return 0;
  324. }
  326. int CDelaunayDoc::GetInitEdges(double x, double y, int p)
  327. {
  328. // GetInitEdges : ransack each trangle to record it's edges
  329. //when the piont belong it's circle,record the triangle's position in m_tri to m_index
  330. //'x' and 'y' are the coordinates of the insert point
  331. //'p' is the mark or th insert point in 'm_point'
  332. // the returned value 'k/2' is the number of cirecle the point belonged
  333. POSITION POS;
  334. CBorder *m_border;
  335. CBorder *m_border1;
  336. CTriangle* pTriangle;
  337. int j=0,i=0,k=0,max=0;
  338. CPointPos *point=new CPointPos(x,y);
  339.     POS = m_tri.GetHeadPosition();
  340. while(POS != NULL ){
  341.        i= m_tri.GetAt( POS )->Where(point);
  342.    
  343.    if(i==POS_ERROR){
  344.    return POS_ERROR;
  345.    }
  346.    if(i==POS_ON)//in circle POS_ON=2
  347.    {
  348.   k=k+i;
  349.   m_index.Add(POS);
  350.   pTriangle=m_tri.GetAt(POS);
  351.   //ransack three edge of a triangle,if an edge of a triangle
  352.   //have belong to the array 'm_pDoc->m_edge' delete both of them to form 
  353.   // the border of the inserted polygon which is stored in the array
  354.   //'m_pDoc->m_edge'
  355.   int array[4];
  356.   array[0]=pTriangle->m_p1;
  357.   array[1]=pTriangle->m_p2;
  358.   array[2]=pTriangle->m_p3;
  359.   array[3]=pTriangle->m_p1;
  360.           for(int h=0;h<3;h++)
  361.   {   
  362.  m_border=new CBorder(array[h],array[h+1]);
  363.  max=m_edge.GetSize();
  364.  if(max==0)
  365.  {
  366.                 m_edge.Add(m_border);
  367.  }
  368.  else
  369.  {
  370.     j=0;
  371.     for(int m=0;m<max;m++)
  372. {             
  373.    m_border1=m_edge[m];
  374.    if(TwoEdgeSuperposition(m_border,m_border1)==1)//?//???
  375.    {
  376.       m_edge.RemoveAt(m,1);
  377.   max=max-1;
  378.   j++;
  379.    }
  380. }
  381.     if(j==0)
  382. {
  383.                    m_edge.Add(m_border);
  384. }     
  385. }  
  386.   }/////for
  387.    }/////////////if(POS_ON)
  388.      pTriangle=m_tri.GetNext(POS);
  389. }//////////////while 
  390. //if(k==0)//the point don't belong any circle,
  391. //i.e.don't the convexity
  392. delete point;
  393. return k;
  394. }
  396. void CDelaunayDoc::DelTriMarked()
  397. {
  398.    int max=m_index.GetSize();
  399.    for(int i=0;i<max;i++)
  400.    {
  401.   m_tri.RemoveAt(m_index[i]);
  402.    }
  403. }
  405. void CDelaunayDoc::EditCon(int r, int l,int p)
  406. {
  407. // r : save the value of i as 右支撑点 or 应去凸包边界节点右边界   
  408. // l : save the value of i as 左支撑点 or 应去凸包边界节点左边界
  409. //'p' is the mark or th insert point in 'm_point'
  410.  CBorder *m_border;
  411.      int hr=0;//record the number of points on edge of convexity that should be deleted before r(include r) 
  412.  int hl=0;//record the number of points on edge of convexity that should be deleted after l(include l)
  413.  int i=r;
  414.  int i1;
  415.  int h=hr+hl;
  416.  int max=m_con.GetSize();
  417.  //search rightward, first
  418.  if((i+1)>max-1)  i1=i+1-max;  
  419.  else             i1=i+1;
  420.  while(S(p,m_con[i],m_con[i1])<0){                  
  421. m_border=new CBorder(m_con[i],m_con[i1]);
  422. m_edge.Add(m_border);
  423. r=i;
  424.             i++;
  425. hr++;
  426. if(i+1>max-1) i1=i+1-max;
  427. else         i1=i+1;
  428. if(i>max-1)     i=i-max;
  429. }
  430. i=l;
  431.     if((i-1)<0) i1=i-1+max;
  432. else        i1=i-1;
  433. while(S(m_con[i1],m_con[i],p)<0){
  434.     m_border=new CBorder(m_con[i],m_con[i1]);
  435. m_edge.Add(m_border);
  436. l=i;
  437. i--;
  438. hl++;
  439. if(i-1<0) i1=i-1+max;
  440. else      i1=i-1;
  441. if(i<0)   i=i+max;
  442. }
  443. h=hr+hl;
  444. if(h>0){
  445.        if(r<l){
  446.        if(hl>0){
  447.   if(hr>0){//both>0
  448.       m_con.RemoveAt(l,max-l);
  449.                   m_con.RemoveAt(0,hr);
  450.   }                            
  451.   else{//hl>0,hr=0
  452.       m_con.RemoveAt(l,max-l);
  453.       m_con.RemoveAt(0,hl-max+l);
  454.   }
  455.           m_con.InsertAt(0,p);                           
  456.   }//end if hl>0
  457.   else{//hl=0,hr>0         
  458.    m_con.RemoveAt(0,hr);
  459.        m_con.InsertAt(0,p);
  460.   }
  461. }//end if r<l
  462. else{//r>l
  463. if(hl>0){
  464.     m_con.RemoveAt(l,h);
  465.         m_con.InsertAt(l,p);                           
  466. }//end if hl>0
  467. else{//hl=0,hr>0         
  468.  m_con.RemoveAt(l+1,hr);
  469.      m_con.InsertAt(l+1,p);
  470. }
  471. }//end if r>l
  472. }//end if h>0
  473. else
  474. {
  475. if((r-l)>0){
  476.      if((r-l)==1) m_con.InsertAt(r,p);
  477.      else{
  478.  m_con.RemoveAt(l+1,r-l-1);
  479.  m_con.InsertAt(l+1,p);                  
  480.  }
  481. }
  482. else{//r<l,i.e. r and l 在 m_con[0] 两侧
  483. if(l<(max-1)) m_con.RemoveAt(l+1,max-l-1);
  484.             m_con.RemoveAt(0,r);
  485. m_con.InsertAt(0,p);       
  486. }
  487. }//h<0             
  488. }
  491. CPointPos* CDelaunayDoc::IntersectionPoint(CPointPos *point1, CPointPos *point2,CPointPos *point3, CPointPos *point4)
  492. {
  493. ASSERT(point1!=NULL);
  494.     ASSERT(point2!=NULL);
  495. ASSERT(point3!=NULL);
  496.     ASSERT(point4!=NULL);
  497.     CPointPos *pos=new CPointPos(0,0);
  498. double k21,k43;
  499. if((point2->m_x-point1->m_x)==0 || (point4->m_x-point3->m_x)==0)
  500. {
  501. if((point2->m_x-point1->m_x)==0 && (point4->m_x-point3->m_x)==0) return NULL;
  502. if((point2->m_x-point1->m_x)==0){
  503. k43=(point4->m_y-point3->m_y)/(point4->m_x-point3->m_x);
  504.             pos->m_x=point1->m_x;// x of the intersection point
  505.             pos->m_y=k43*(pos->m_x-point3->m_x)+point3->m_y;// y of the intersection point
  506.     return pos;
  507. }
  508. if((point4->m_x-point3->m_x)==0){
  509. k21=(point2->m_y-point1->m_y)/(point2->m_x-point1->m_x);
  510.             pos->m_x=point3->m_x;// x of the intersection point
  511.             pos->m_y=k21*(pos->m_x-point1->m_x)+point1->m_y;// y of the intersection point
  512.     return pos;
  513. }
  514. }
  515. else
  516. {
  517. k43=(point4->m_y-point3->m_y)/(point4->m_x-point3->m_x);
  518.         k21=(point2->m_y-point1->m_y)/(point2->m_x-point1->m_x);
  519. if(fabs(k43-k21)<0.0000000001) 
  520. return NULL;
  521. else{
  522. pos->m_x=(point3->m_y-point1->m_y+k21*point1->m_x-k43*point3->m_x)/(k21-k43);
  523.             pos->m_y=k43*(pos->m_x-point3->m_x)+point3->m_y;
  524. return pos;
  525. }
  526. }       
  527. }
  529. int CDelaunayDoc::TheOtherPoint(int p1, int p2, CTriangle *temp)
  530. {//p1,p2 are the mark recorded by m_tri in m_point
  531. int array[3];
  532. double s;
  533. array[0]=temp->m_p1;
  534. array[1]=temp->m_p2;
  535. array[2]=temp->m_p3;
  536. for(int i=0;i<3;i++)
  537. {
  538. s=S(array[i],p1,p2);
  539. if(fabs(s)>0.00000000001) return array[i];
  540. }
  541. }
  543. bool CDelaunayDoc::DelEdgeOrNot(int p1, int p2,int p)
  544. {
  545.    CTriangle* pTriangle;
  546.    //p1=m_pDoc->m_edge[j]->m_p1,record the point's mark in m_point
  547.    //p2=m_pDoc->m_edge[j]->m_p2,record the point's mark in m_point
  548.    int i,j;
  549.    for(i=0;i<m_index.GetSize();i++)
  550.    {
  551.    pTriangle=m_tri.GetAt(m_index[i]);
  552.    int array[3];
  553.    array[0]=pTriangle->m_p1;
  554.    array[1]=pTriangle->m_p2;
  555.    array[2]=pTriangle->m_p3;   
  556.    int a=-1,b=-1;
  557.    for(j=0;j<3;j++)
  558.    {
  559.    if(p1==array[j]) a=j;
  560.    }
  561.     for(j=0;j<3;j++)
  562.    {
  563.    if(p2==array[j]) b=j;
  564.    }
  565.        if((a>=0) & (b>=0))
  566.    {
  567.    int other_point=TheOtherPoint(p1,p2,pTriangle);
  568.    CPointPos* pPoint=IntersectionPoint(m_point[p1],m_point[p2],
  569.                       m_point[p],m_point[other_point]);
  570.    if(pTriangle->Where(pPoint)==POS_ON){//整数截断使交点在边外
  571.    return true;
  572.    }
  573.    }
  574.    }
  575.    return false;
  576. }
  578.  double CDelaunayDoc::S(CPointPos *p1,CPointPos *p2,CPointPos *p3)
  579.  {
  580.   //求三角形面积,以右下角为原点时,s>0为逆时针(左转),
  581.    //s=0为三点重合
  582.    double s;
  583.    s=p1->m_x*p2->m_y+p2->m_x*p3->m_y+p1->m_y*p3->m_x-
  584. p2->m_y*p3->m_x-p1->m_y*p2->m_x-p1->m_x*p3->m_y;
  585.    s=s/2.0;
  586.    return s;
  587.  }
  589. void CDelaunayDoc::OnButtonBB() 
  590. {
  591. m_DoWhat=DO_INTERPOLATION;
  592. int max=m_point.GetSize();
  593. for(int i=0;i<max;i++)
  594. {
  595. FindRelativeTri(i);//get Fx,Fy,N of a point in m_point
  596. }
  597. POSITION POS;CTriangle* temp;
  598. POS = m_tri.GetHeadPosition();
  599. temp=m_tri.GetAt(POS);
  600. while(POS != NULL ){
  601. temp=m_tri.GetNext(POS);
  602. double d;
  603. int p[4];p[0]=0;
  604.     p[1]=temp->m_p1;p[2]=temp->m_p2;p[3]=temp->m_p3;
  605. //////////////////////////////////////////////////////////////////////////
  606. //Get 6个 dij 计算无错 但沿边接的方向倒数的估算可能有问题,
  607. //1)可能不应单位化.2)?
  608. d=D(temp,2,3);
  609. temp->d12=m_point[p[2]]->m_z+d/double(3);
  611. d=D(temp,3,2);
  612. temp->d13=m_point[p[3]]->m_z+d/double(3);
  613.         
  614. d=D(temp,1,3);
  615. temp->d21=m_point[p[1]]->m_z+d/double(3);
  616.         
  617. d=D(temp,3,1);
  618. temp->d23=m_point[p[3]]->m_z+d/double(3);
  620. d=D(temp,1,2);
  621. temp->d31=m_point[p[1]]->m_z+d/double(3);
  623.         d=D(temp,2,1);
  624. temp->d32=m_point[p[2]]->m_z+d/double(3);
  625. //////////////////////////////////////////////////////////////////////////
  626. //计算无错
  627.         temp->c1=(temp->d21+temp->d31+m_point[p[1]]->m_z)/double(3);
  628. temp->c2=(temp->d12+temp->d32+m_point[p[2]]->m_z)/double(3);
  629. temp->c3=(temp->d13+temp->d23+m_point[p[3]]->m_z)/double(3);
  630. //////////////////////////////////////////////////////////////////////////
  631. /*Get 每边中点处f的法向导数值(XY平面上的) 
  632. 在这边的两顶点间对顶点的法向导数做插值
  633.  顶点的法向导数由该点的Fx,Fy取得,*/     
  634. temp->n1=F(temp,3,2);//计算无错
  635. temp->n2=F(temp,1,3);
  636. temp->n3=F(temp,2,1);
  638. //////////////////////////////////////////////////////////////////////////
  639. /*Get bi and b0(o), 220页公式(5.48) 应再推几遍.HCT 和 三角形面积 都用的是
  640. XY平面上的?????3维的又如何????*///计算无错
  641. CPointPos* h1;
  642. CPointPos* h2;
  643. CPointPos* h3;
  644. //三角形o,f2,f3
  645.         h1=GetChuiZu(temp->m_x,temp->m_y,m_point[p[2]],m_point[p[3]]);
  646. //三角形o,f3,f1
  647. h2=GetChuiZu(temp->m_x,temp->m_y,m_point[p[3]],m_point[p[1]]);
  648. //三角形o,f1,f2
  649. h3=GetChuiZu(temp->m_x,temp->m_y,m_point[p[1]],m_point[p[2]]);
  650. CPointPos* o=new CPointPos();//重心
  651. o->m_x=temp->m_x;
  652. o->m_y=temp->m_y;
  653. o->m_z=0;
  654. double s1=S(o,m_point[p[2]],m_point[p[3]]);
  655. double s2=S(o,m_point[p[3]],m_point[p[1]]);
  656. double s3=S(o,m_point[p[1]],m_point[p[2]]);
  657. temp->b1=S(o,m_point[p[3]],h2)/s2
  658.         *(m_point[p[1]]->m_z+temp->d21-temp->d23-m_point[p[3]]->m_z)+
  659.         S(o,m_point[p[1]],h3)/s3
  660.         *(m_point[p[2]]->m_z+temp->d32-temp->d31-m_point[p[1]]->m_z)-
  661. (temp->n3+temp->n2)/double(0.75)-temp->c3-temp->c2+
  662. m_point[p[1]]->m_z+m_point[p[3]]->m_z+temp->d21+temp->d32+(temp->d31+temp->d23)*
  663. double(2);
  664. temp->b1=temp->b1/double(6);
  666. temp->b2=S(o,m_point[p[2]],h1)/s1
  667.         *(m_point[p[3]]->m_z+temp->d13-temp->d12-m_point[p[2]]->m_z)+
  668.         S(o,m_point[p[1]],h3)/s3
  669.         *(m_point[p[2]]->m_z+temp->d32-temp->d31-m_point[p[1]]->m_z)-
  670. (temp->n3+temp->n1)/double(0.75)-temp->c3-temp->c1+
  671. m_point[p[1]]->m_z+m_point[p[2]]->m_z+temp->d32+temp->d13+(temp->d31+temp->d12)*
  672. double(2);
  673. temp->b2=temp->b2/double(6);
  675. temp->b3=S(o,m_point[p[2]],h1)/s1
  676.         *(m_point[p[3]]->m_z+temp->d13-temp->d12-m_point[p[2]]->m_z)+
  677.         S(o,m_point[p[3]],h2)/s2
  678.         *(m_point[p[1]]->m_z+temp->d21-temp->d23-m_point[p[3]]->m_z)-
  679. (temp->n1+temp->n2)/double(0.75)-temp->c1-temp->c2+
  680. m_point[p[2]]->m_z+m_point[p[3]]->m_z+temp->d13+temp->d21+(temp->d12+temp->d23)*
  681. double(2);
  682. temp->b3=temp->b3/double(6);
  683. delete o;
  684. temp->o=(temp->b1+temp->b2+temp->b3)/double(3);
  685. //////////////////////////////////////////////////////////////////////////
  686.        //Get ei//计算无错
  687.        temp->e1=double(3)*(-temp->b1+temp->b2+temp->b3)-(-temp->c1+temp->c2+temp->c3);
  688.    temp->e1=temp->e1/double(2);
  689.    temp->e2=double(3)*(temp->b1-temp->b2+temp->b3)-(temp->c1-temp->c2+temp->c3);
  690.    temp->e2=temp->e2/double(2);
  691.    temp->e3=double(3)*(temp->b1+temp->b2-temp->b3)-(temp->c1+temp->c2-temp->c3);
  692.    temp->e3=temp->e3/double(2);
  693. //////////////////////////////////////////////////////////////////////////
  694.    //NOW we have got the 19 control points in a triangle
  695. //       temp=m_tri.GetNext(POS);
  697.    delete h1;
  698.    delete h2;
  699.    delete h3;
  700. }   
  701. }
  703. void CDelaunayDoc::OnUpdateButtonBB(CCmdUI* pCmdUI) 
  704. {
  705. if(m_DoWhat==DO_INTERPOLATION)
  706. pCmdUI->SetCheck(1);
  707. else 
  708. pCmdUI->SetCheck(0);
  709. }
  711. void CDelaunayDoc::FindRelativeTri(int p)
  712. {
  713. CPointPos *point;
  714.     point=m_point[p];
  715. POSITION POS;
  716.     POS = m_tri.GetHeadPosition();
  717. while(POS != NULL ){
  718.   int array[3];
  719.   array[0]= m_tri.GetAt(POS)->m_p1;
  720.   array[1]= m_tri.GetAt(POS)->m_p2;
  721.   array[2]= m_tri.GetAt(POS)->m_p3;
  722.           for(int h=0;h<3;h++)
  723.   {     
  724.  if(point==m_point[array[h]])
  725.  {
  726.                 m_index.Add(POS);
  727.  }
  728.   }/////for
  729.       m_tri.GetNext(POS);
  730. }//////////////while 
  731. Get_Fx_Fy_N(p);
  732. }
  734. POI CDelaunayDoc::GetTriNormal(CTriangle *temp)
  735. {//得到三角片的单位法向
  736. POI a,b,n;
  737. int p1,p2,p3;
  738.     p1=temp->m_p1;
  739.     p2=temp->m_p2;
  740.     p3=temp->m_p3;
  741. a.x=m_point[p2]->m_x-m_point[p1]->m_x;
  742. a.y=m_point[p2]->m_y-m_point[p1]->m_y;
  743. a.z=m_point[p2]->m_z-m_point[p1]->m_z;
  744. b.x=m_point[p3]->m_x-m_point[p1]->m_x;
  745. b.y=m_point[p3]->m_y-m_point[p1]->m_y;
  746. b.z=m_point[p3]->m_z-m_point[p1]->m_z;
  747.     //get the normal n=a*b(外积)
  748. n=VectorProduct(a.x,a.y,a.z,b.x,b.y,b.z);
  749. //单位化 normal
  750. n=Unitization(n);
  751.     return n;
  752. }
  754. double CDelaunayDoc::GetDistance(CPointPos* p1, CPointPos* p2)
  755. {
  756. double distance;
  757. distance=(p1->m_x-p2->m_x)*(p1->m_x-p2->m_x)+
  758.      (p1->m_y-p2->m_y)*(p1->m_y-p2->m_y)+
  759.              (p1->m_z-p2->m_z)*(p1->m_z-p2->m_z);
  760. distance=sqrt(distance);
  761. return distance;
  762. }
  764. POI CDelaunayDoc::VectorProduct(double x1,double y1,double z1,double x2,double y2,double z2)
  765. {//get p1*p2(外积) 对!!!!
  766. POI product;
  767. product.x=y1*z2-z1*y2;
  768. product.y=z1*x2-x1*z2;
  769. product.z=x1*y2-y1*x2;
  770. return product;
  771. }
  773. POI CDelaunayDoc::Unitization(POI p)
  774. { //单位化 对!!!!
  775. double m;
  776. m=p.x*p.x+p.y*p.y+p.z*p.z;
  777. m=sqrt(m);
  778. p.x=p.x/m;
  779. p.y=p.y/m;
  780. p.z=p.z/m;
  781.     return p;
  782. }
  784. double CDelaunayDoc::DotProduction(double x1,double y1,double z1,double x2,double y2,double z2)
  785. {//(x1,y1,z1)(内积)(x2,y2,z2)
  786. double product;
  787.     product=x1*x2+y1*y2+z1*z2;
  788.     return product;
  789. }
  791. void CDelaunayDoc::Get_Fx_Fy_N(int p)
  792. {/*本处计算无错,但型值点法向没用,因算得是Fx and Fy,
  793. 然后用 Fx and Fy 在 xy 平面上沿边的方向计算方向倒数
  794. 没用型值点法向做沿空间三角形边的方向倒数:
  795.       D(i,i+1)=(Pi+1-Pi)-[(Pi+1-Pi)Ni]Ni
  796.   D(i,i+2)=(Pi+2-Pi)-[(Pi+2-Pi)Ni]Ni
  797. */
  798.    CPointPos *point;
  799.    point=m_point[p];
  800.    int max=m_index.GetSize();
  801.    double Fx=0;//点的x偏导
  802.    double Fy=0;
  803.   // POI n;//点的法向量
  804.   //n.x=0;n.y=0;n.z=0;
  805.    double w=0.0;//Little法
  806.    for(int i=0;i<max;i++)
  807.    {
  808.   CArray<double,double> a;//save 三角片所在平面的方程ax+by+cz+d=0的a 
  809.   CArray<double,double> b;//save 三角片所在平面的方程ax+by+cz+d=0的b
  810.   CArray<double,double> c;//save 三角片所在平面的方程ax+by+cz+d=0的c
  811.   int array[3];
  812.   CTriangle* Tri;
  813.   Tri=m_tri.GetAt(m_index[i]);
  814.   array[0]= Tri->m_p1;
  815.   array[1]= Tri->m_p2;
  816.   array[2]= Tri->m_p3;
  817.   //a,b,c 计算无错 对!!!!
  818.       a.SetAtGrow(i,(m_point[array[1]]->m_y-m_point[array[0]]->m_y)*
  819.    (m_point[array[2]]->m_z-m_point[array[0]]->m_z)-
  820.    (m_point[array[1]]->m_z-m_point[array[0]]->m_z)*
  821.    (m_point[array[2]]->m_y-m_point[array[0]]->m_y));
  822.   b.SetAtGrow(i,(m_point[array[1]]->m_z-m_point[array[0]]->m_z)*
  823.    (m_point[array[2]]->m_x-m_point[array[0]]->m_x)-
  824.    (m_point[array[1]]->m_x-m_point[array[0]]->m_x)*
  825.    (m_point[array[2]]->m_z-m_point[array[0]]->m_z));
  826.   c.SetAtGrow(i,(m_point[array[1]]->m_x-m_point[array[0]]->m_x)*
  827.    (m_point[array[2]]->m_y-m_point[array[0]]->m_y)-
  828.    (m_point[array[1]]->m_y-m_point[array[0]]->m_y)*
  829.    (m_point[array[2]]->m_x-m_point[array[0]]->m_x));
  830.   
  831.   CArray<CPointPos*,CPointPos*> P;
  832.   double l1,l2;   
  833.       for(int h=0;h<3;h++)
  834.   {     
  835.  if(point!=m_point[array[h]])
  836.  {
  837.             P.Add(m_point[array[h]]);
  838.  }
  839.   }/////for
  840.   l1=GetDistance(point,P[0]);
  841.   l2=GetDistance(point,P[1]);
  842.   P.RemoveAll();////????
  843.   l1=l1*l1;
  844.   l2=l2*l2;
  845.       w=w+double(1)/(l1*l2);//w ,Fx,Fy and n have been 初始化了 =0  
  846.   /*n.x=n.x+GetTriNormal(Tri).x/(l1*l2);
  847.   n.y=n.y+GetTriNormal(Tri).y/(l1*l2);
  848.   n.z=n.z+GetTriNormal(Tri).z/(l1*l2);*/
  849.   Fx=(-a[i]/c[i])/(l1*l2)+Fx;Fy=(-b[i]/c[i])/(l1*l2)+Fy;
  850.    }
  851.    /*m_point[p]->nx=n.x/w;///?????nx,ny,nz 有用吗
  852.    m_point[p]->ny=n.y/w;
  853.    m_point[p]->nz=n.z/w;*/
  854.    m_point[p]->Fx=Fx/w;m_point[p]->Fy=Fy/w;
  855.    m_index.RemoveAll();
  856. }
  858. void CDelaunayDoc::BaryCenter(CTriangle *temp)
  859. {//三角形重心
  860. ASSERT(temp!=NULL);
  861. int p1,p2,p3;
  862. p1=temp->m_p1;
  863. p2=temp->m_p2;
  864. p3=temp->m_p3;   
  865. CPointPos* p11=m_point[p1];
  866. CPointPos* p22=m_point[p2];
  867. CPointPos* p33=m_point[p3];
  868. temp->m_x=(p11->m_x+p22->m_x +p33->m_x )/double(3);
  869. temp->m_y=(p11->m_y+p22->m_y +p33->m_y )/double(3);
  870.     //temp->m_z=(p11->m_z+p22->m_z +p33->m_z )/double(3);
  871. }
  873. double CDelaunayDoc::GetMold(CPointPos *p)//取模
  874. {
  875. double m=0;
  876. m=p->m_x*p->m_x+p->m_y*p->m_y+p->m_z*p->m_z;
  877. m=sqrt(m);
  878. return m;
  880. }
  882. CPointPos* CDelaunayDoc::GetChuiZu(double x,double y, CPointPos* p2, CPointPos* p3)
  883. {//求 p1 对于p2和p3 的垂足.
  884. CPointPos* h=new CPointPos();
  885. double k23;
  886. if (p3->m_x-p2->m_x==0){
  887. h->m_x=p2->m_x;
  888. h->m_y=y;
  889. }
  890. else{
  891. k23=(p3->m_y-p2->m_y)/(p3->m_x-p2->m_x);
  892. if (k23==0){
  893. h->m_x=x;
  894. h->m_y=p2->m_y;
  895. }
  896. else{
  897.             h->m_x=(k23*k23*p2->m_x+x+(y-p2->m_y)*k23)/(1.0+k23*k23);
  898.             h->m_y=p2->m_y+k23*(h->m_x-p2->m_x);
  899. }
  900. }
  901. return h;
  902. }
  904. double CDelaunayDoc::D(CTriangle *temp, int i, int j)
  905. {//沿边方向导:i to j,用i点处的两个偏导Fx,Fy 和i to j 的边的方向上的单位矢量作内积
  906. int p[4];p[0]=0;
  907. p[1]=temp->m_p1;p[2]=temp->m_p2;p[3]=temp->m_p3;
  908. POI vector;
  909. double d;
  910. vector.x=m_point[p[j]]->m_x-m_point[p[i]]->m_x;
  911. vector.y=m_point[p[j]]->m_y-m_point[p[i]]->m_y;
  912.     vector=Unitization(vector);//要是单位化就错了
  913.     d=DotProduction(vector.x,vector.y,0,m_point[p[i]]->Fx,m_point[p[i]]->Fy,0);
  914. return d;
  915. }
  917. double CDelaunayDoc::F(CTriangle *temp, int i, int j)
  918. {//i to j,ij边上的法向导数
  919. int p[4];p[0]=0;
  920. p[1]=temp->m_p1;p[2]=temp->m_p2;p[3]=temp->m_p3;
  921. int h=TheOtherPoint(p[i],p[j],temp);
  922.     CPointPos* H=GetChuiZu(m_point[h]->m_x,m_point[h]->m_y,m_point[p[i]],m_point[p[j]]);
  923. POI n;
  924. n.x=H->m_x-m_point[h]->m_x;
  925. n.y=H->m_y-m_point[h]->m_y;
  926. n.z=0;
  927. //n=Unitization(n);//要是单位化就错了
  928. double n1,n2;
  929. n1=DotProduction(m_point[p[j]]->Fx,m_point[p[j]]->Fy,0,n.x,n.y,0);
  930. n2=DotProduction(m_point[p[i]]->Fx,m_point[p[i]]->Fy,0,n.x,n.y,0);
  931. delete H;
  932. return (n1+n2)/double(2);
  933. }
  935. void CDelaunayDoc::DrawTri(int m_p1,int m_p2,CTriangle *tri)
  936. {
  937. int i,j;
  938. POI p1,p2,o;
  939. p1.x=m_point[m_p1]->m_x;
  940. p1.y=m_point[m_p1]->m_y;
  941. p2.x=m_point[m_p2]->m_x;
  942. p2.y=m_point[m_p2]->m_y;
  943. o.x=tri->m_x;o.y=tri->m_y;
  944. int n;
  945. n=8;
  946. double x[8][8];
  947. double y[8][8];
  949. for(j=0;j<n;j++)
  950. {
  951. for(i=0;i<n-j;i++)
  952. {
  953. x[j][i]=p1.x+j*(p2.x-p1.x)/double(n-1)+
  954. i*(o.x+j*(p2.x-o.x)/double(n-1)-p1.x-j*(p2.x-p1.x)/double(n-1))
  955. /double(n-1-j);
  956. y[j][i]=p1.y+j*(p2.y-p1.y)/double(n-1)+
  957. i*(o.y+j*(p2.y-o.y)/double(n-1)-p1.y-j*(p2.y-p1.y)/double(n-1))
  958. /double(n-1-j);
  959. }
  960. }
  961. x[7][0]=p2.x;
  962. y[7][0]=p2.y;
  964.     for(j=0;j<n-1;j++)
  965. {
  966. for(i=0;i<n-j-1;i++)
  967. {
  968. //normal=GetTriNormal(m_hct[0],m_hct[1],m_hct[2]);
  969. //glNormal3d(normal.x,normal.y,normal.z);
  970. glBegin(GL_TRIANGLES);
  971.   glVertex3d(x[j][i],y[j][i],Bezier(x[j][i],y[j][i],m_p1,m_p2,tri));
  972.   glVertex3d(x[j+1][i],y[j+1][i],Bezier(x[j+1][i],y[j+1][i],m_p1,m_p2,tri));
  973.   glVertex3d(x[j][i+1],y[j][i+1],Bezier(x[j][i+1],y[j][i+1],m_p1,m_p2,tri));
  974. glEnd();
  975. }
  976. }
  977. for(j=1;j<n-1;j++)
  978. {
  979. for(i=0;i<n-j-1;i++)
  980. {
  981. //normal=GetTriNormal(m_hct[0],m_hct[1],m_hct[2]);
  982. //glNormal3d(normal.x,normal.y,normal.z);
  983. glBegin(GL_TRIANGLES);
  984.   glVertex3d(x[j][i],y[j][i],Bezier(x[j][i],y[j][i],m_p1,m_p2,tri));
  985.   glVertex3d(x[j][i+1],y[j][i+1],Bezier(x[j][i+1],y[j][i+1],m_p1,m_p2,tri));
  986.   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)); 
  987. glEnd();
  988. }
  989. }
  990. /*
  991. for(j=1;j<n;j++)
  992. {
  993. for(i=0;i<n-j;i++)
  994. {
  995. glBegin(GL_LINES);
  996.   glVertex3d(x[j][i],y[j][i],Bezier(x[j][i],y[j][i],m_p1,m_p2,tri));
  997.   glVertex3d(x[j-1][i],y[j-1][i],Bezier(x[j-1][i],y[j-1][i],m_p1,m_p2,tri));
  998.   glVertex3d(x[j][i],y[j][i],Bezier(x[j][i],y[j][i],m_p1,m_p2,tri));
  999.   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));
  1000. glEnd();
  1001. }
  1002. }
  1004. for(j=0;j<n-1;j++)
  1005. {
  1006. for(i=0;i<n-j-1;i++)
  1007. {
  1008. glBegin(GL_LINES);
  1009.   glVertex3d(x[j][i],y[j][i],Bezier(x[j][i],y[j][i],m_p1,m_p2,tri));
  1010.   glVertex3d(x[j][i+1],y[j][i+1],Bezier(x[j][i+1],y[j][i+1],m_p1,m_p2,tri));
  1011. glEnd();
  1012. }
  1013. }
  1014. */
  1016. }
  1018. double CDelaunayDoc::Bezier(double x, double y,int m_p1,int m_p2,CTriangle* tri)
  1019. {
  1020. //s1,s2,s3 : 面积坐标
  1021. double s1,s2,s3,s;
  1022. double b[4][4][4];
  1023. POI p1,p2,o,p;
  1024. p1.x=m_point[m_p1]->m_x;
  1025. p1.y=m_point[m_p1]->m_y;
  1026. p2.x=m_point[m_p2]->m_x;
  1027. p2.y=m_point[m_p2]->m_y;
  1028. o.x=tri->m_x;o.y=tri->m_y;
  1029. p.x=x;p.y=y;
  1031. s=S(o,p1,p2);
  1032. s1=S(p,p1,p2)/s;
  1033. s2=S(o,p,p2)/s;
  1034. s3=S(o,p1,p)/s;
  1035. if(fabs(s1)<0.000000000001) s1=0;
  1036. if(fabs(s2)<0.000000000001) s2=0;
  1037. if(fabs(s3)<0.000000000001) s3=0;
  1038. int i,j,k;
  1039. double B=0.0;
  1040. 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;
  1041. if(m_p1==tri->m_p1)
  1042. {
  1043. b[1][1][1]=tri->e3;
  1044. b[2][1][0]=tri->b1;b[1][2][0]=tri->c1;
  1045. b[0][2][1]=tri->d31;b[0][1][2]=tri->d32;
  1046. b[2][0][1]=tri->b2;b[1][0][2]=tri->c2;
  1048. for(i=0;i<4;i++){
  1049. for(j=0;j<=(3-i);j++){
  1050. for(k=0;k<=(3-i-j);k++){
  1051. if((i+j+k)==3){
  1052. B=B+b[i][j][k]*double(6)/(Factorial(i)*Factorial(j)*Factorial(k))*
  1053. Power(s1,i)*Power(s2,j)*Power(s3,k);
  1054. }
  1055. }
  1056. }
  1057. }
  1058. }
  1059. if(m_p1==tri->m_p2)
  1060. {
  1061. b[1][1][1]=tri->e1;
  1062. b[2][1][0]=tri->b2;b[1][2][0]=tri->c2;
  1063. b[0][2][1]=tri->d12;b[0][1][2]=tri->d13;
  1064. b[2][0][1]=tri->b3;b[1][0][2]=tri->c3;
  1066. for(i=0;i<4;i++){
  1067. for(j=0;j<=(3-i);j++){
  1068. for(k=0;k<=(3-i-j);k++){
  1069. if((i+j+k)==3){
  1070. B=B+b[i][j][k]*double(6)/(Factorial(i)*Factorial(j)*Factorial(k))*
  1071. Power(s1,i)*Power(s2,j)*Power(s3,k);
  1072. }
  1073. }
  1074. }
  1075. }
  1076. }
  1077. if(m_p1==tri->m_p3)
  1078. {
  1079. b[1][1][1]=tri->e2;
  1080. b[2][1][0]=tri->b3;b[1][2][0]=tri->c3;
  1081. b[0][2][1]=tri->d23;b[0][1][2]=tri->d21;
  1082. b[2][0][1]=tri->b1;b[1][0][2]=tri->c1;
  1084. for(i=0;i<4;i++){
  1085. for(j=0;j<=(3-i);j++){
  1086. for(k=0;k<=(3-i-j);k++){
  1087. if((i+j+k)==3){
  1088. B=B+b[i][j][k]*double(6)/(Factorial(i)*Factorial(j)*Factorial(k))*
  1089. Power(s1,i)*Power(s2,j)*Power(s3,k);
  1090. }
  1091. }
  1092. }
  1093. }
  1094. }
  1095. return B;
  1097. }
  1099. double CDelaunayDoc::S(POI p1, POI p2, POI p3)
  1100. {
  1101. //求三角形面积,以右下角为原点时,s>0为逆时针(左转),
  1102.    //s=0为三点重合
  1103.    double s;
  1104.    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;
  1105.    s=s/2.0;
  1106.    return s;
  1107. }
  1109. int CDelaunayDoc::Factorial(int n)
  1110. {//阶乘
  1111. int x=n;
  1112. if(n==0 ||n==1) 
  1113. return 1;
  1114. else{
  1115. for(int i=1;i<n;i++){
  1116. x=x*(n-i);
  1117. }
  1118. }
  1119.     return x;
  1120. }
  1122. double CDelaunayDoc::Power(double a, int e)
  1123. {
  1124. double x=1.0;
  1125. if(e==0){
  1126. return 1.0;
  1127. }
  1128. else{
  1129. if(a==0.0){
  1130. return 0.0;
  1131. }
  1132.         for(int i=1;i<=e;i++){
  1133. x=x*a;
  1134. }
  1135. }
  1136.     return x;
  1137. }
  1139. CTriangle* CDelaunayDoc::Belong(double x, double y)
  1140. {
  1141. CTriangle* tri;
  1142. POI p1,p2,p3,p;
  1143. p.x=x;p.y=y;
  1144. POSITION pos = m_tri.GetHeadPosition();
  1145. while( pos != NULL )
  1146. {
  1147.         tri=m_tri.GetAt( pos );
  1148. double s1,s2,s3;
  1149.     p1.x=m_point[tri->m_p1]->m_x;
  1150. p1.y=m_point[tri->m_p1]->m_y;
  1151. p2.x=m_point[tri->m_p2]->m_x;
  1152. p2.y=m_point[tri->m_p2]->m_y;
  1153.     p3.x=m_point[tri->m_p3]->m_x;
  1154. p3.y=m_point[tri->m_p3]->m_y;
  1155.     
  1156. s1=S(p,p1,p2);
  1157. s2=S(p,p2,p3);
  1158. s3=S(p,p3,p1);
  1159. if(s1>0 && s2>0 && s3>0)
  1160. return tri;
  1162. m_tri.GetNext( pos );
  1163. }
  1164.     return NULL;
  1165. }
  1167. int CDelaunayDoc::Belong(double x, double y, CTriangle *tri)
  1168. {
  1169. POI p1,p2,o,p;
  1170. p.x=x;p.y=y;
  1171. double s1,s2,s3;
  1172. o.x=tri->m_x;
  1173. o.y=tri->m_y;
  1174. int a[4];
  1175. a[0]=tri->m_p1;
  1176. a[1]=tri->m_p2;
  1177. a[2]=tri->m_p3;
  1178. a[3]=tri->m_p1;
  1179. for(int i=0;i<3;i++)
  1180. {
  1181. p1.x=m_point[a[i]]->m_x;
  1182. p1.y=m_point[a[i]]->m_y;
  1183. p2.x=m_point[a[i+1]]->m_x;
  1184. p2.y=m_point[a[i+1]]->m_y;
  1185. s1=S(p,o,p1);
  1186. s2=S(p,p1,p2);
  1187. s3=S(p,p2,o);
  1188. if(s1>0 && s2>0 && s3>0)
  1189.    return i;
  1190. }
  1191. }
  1193. void CDelaunayDoc::Wang()
  1194. {
  1195. int n=50;
  1196.     CTriangle* tri;
  1197. double x[60];
  1198. double y[60];
  1199. double z[60][60];
  1200. int i,j,what;
  1201. for(i=0;i<n;i++)
  1202. {
  1203. x[i]=double(i)/double(n);
  1204.     y[i]=double(i)/double(n);
  1205. }
  1207. for(i=0;i<n;i++)
  1208. {
  1209. for(j=0;j<n;j++)
  1210. {
  1211. tri=Belong(x[i],y[j]);
  1212. if(tri!=NULL)
  1213. {
  1214.    what=Belong(x[i],y[j],tri);
  1215.    switch(what)
  1216.    {
  1217.    case 0:
  1218.    z[i][j]=Bezier(x[i],y[j],tri->m_p1,tri->m_p2,tri);
  1219.    break;
  1220.    case 1:
  1221.    z[i][j]=Bezier(x[i],y[j],tri->m_p2,tri->m_p3,tri);
  1222.    break;
  1223.    default:
  1224.    z[i][j]=Bezier(x[i],y[j],tri->m_p3,tri->m_p1,tri);
  1225.    break;
  1226.    }
  1227.    //z[i][j]=z[i][j]-0.2;
  1228. }
  1229. else
  1230. z[i][j]=0;
  1231. }
  1232. }
  1233. for(i=0;i<n;i++)
  1234. {
  1235. for(j=0;j<n-1;j++)
  1236. {
  1237.     if((z[i][j]!=0) && (z[i][j+1]!=0))
  1238. {
  1239.    glBegin(GL_LINES);
  1240.  glVertex3d(x[i],y[j],z[i][j]);
  1241.  glVertex3d(x[i],y[j+1],z[i][j+1]);
  1242.    glEnd();
  1243. }
  1244. }
  1245. }
  1246. for(j=0;j<n;j++)
  1247. {
  1248. for(i=0;i<n-1;i++)
  1249. {
  1250. if((z[i][j]!=0) && (z[i+1][j]!=0))
  1251. {
  1252.    glBegin(GL_LINES);
  1253.     glVertex3d(x[i],y[j],z[i][j]);
  1254. glVertex3d(x[i+1],y[j],z[i+1][j]);
  1255.        glEnd();
  1256. }
  1258. }
  1259. }
  1261. }