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

源码开发语言/平台
当前位置：首页 > 源码/资料 > 中间件编程 > 查看源码
Triangulate.cpp
资源名称：08.zip　[点击查看]
上传用户：ynjin1970
上传日期：2014-10-13
资源大小：6438k
文件大小：16k
源码类别：

中间件编程

开发平台：

Visual C++

Triangulate.cpp：源码内容
  1. // GenerateTri1.cpp: implementation of the CTriangulate class.
  2. //
  3. //////////////////////////////////////////////////////////////////////
  5. #include "stdafx.h"
  6. //#include "ct.h"
  7. #include "Triangulate.h"
  8. #include "globalfunctions.h"
  9. #include <math.h>
  10. #include <strstrea.h>
  12. #ifdef _DEBUG
  13. #undef THIS_FILE
  14. static char THIS_FILE[]=__FILE__;
  15. #define new DEBUG_NEW
  16. #endif
  18. //////////////////////////////////////////////////////////////////////
  19. // Construction/Destruction
  20. //////////////////////////////////////////////////////////////////////
  22. CTriangulate::CTriangulate()
  23. {
  24. nv_ = 0;
  25. ntri_ = 0;
  26. }
  28. CTriangulate::~CTriangulate()
  29. {
  30. delete []pv_;
  31. delete ptri_;
  32. }
  36. /*******************************************************************
  37.    Triangulation subroutine
  38.    Takes as input NV vertices in array pxyz
  39.    Returned is a list of ntri triangular faces in the array v
  40.    These triangles are arranged in a consistent clockwise order.
  41.    The triangle array 'v' should be malloced to 3 * nv
  42.    The vertex array pxyz must be big enough to hold 3 more points
  43.    The vertex array must be sorted in increasing x values say
  45.    qsort(p,nv,sizeof(XYZ),XYZCompare);
  46.       :
  47.    int XYZCompare(void *v1,void *v2)
  48.    {
  49.       XYZ *p1,*p2;
  50.       p1 = v1;
  51.       p2 = v2;
  52.       if (p1->x < p2->x)
  53.          return(-1);
  54.       else if (p1->x > p2->x)
  55.          return(1);
  56.       else
  57.          return(0);
  58.    }
  59. ******************************************************************/
  61. int CTriangulate::Triangulate(int nv, XYZ *pxyz, int *ntri, TRIANGLE *v)
  62.   // nv为存储在顶点数组pxyz中的顶点的个数
  63.   // 返回存储在V中的三角形，其个数由ntri指定
  64. {
  65.    int *complete = NULL; //在点已按x轴排序的情况下表示三角形是否需要处理
  66.    IEDGE *edges = NULL;
  67.    int nedge = 0;
  68.    int trimax; // 三角形的最多个数
  69.    int emax = 200; //边缓冲区的最大长度，会自然增长。与edges关联
  70.    int status = 0;  //非0表示错误
  72.    int inside;
  73.    int i,j,k;
  74.    double xp,yp,x1,y1,x2,y2,x3,y3,xc,yc,r; //判断点是否在外接圆内
  75.    double xmin,xmax,ymin,ymax,xmid,ymid,dx, dy, dmax; //超级三角形
  77.    /* Allocate memory for the completeness list, flag for each triangle */
  78.    trimax = 4 * nv;
  79.    if ((complete = (int *) malloc(trimax*sizeof(int))) == NULL) {
  80.       status = 1;
  81.       goto skip;
  82.    }
  84.    /* Allocate memory for the edge list */
  85.    if ((edges = (IEDGE *) malloc(emax*(long)sizeof(IEDGE))) == NULL) {
  86.       status = 2;
  87.       goto skip;
  88.    }
  90.    /*
  91.       Find the maximum and minimum vertex bounds.
  92.       This is to allow calculation of the bounding triangle
  93.    */
  94.    xmin = pxyz[0].x;
  95.    ymin = pxyz[0].y;
  96.    xmax = xmin;
  97.    ymax = ymin;
  98.    for (i=1;i<nv;i++) {
  99.       if (pxyz[i].x < xmin) xmin = pxyz[i].x;
  100.       if (pxyz[i].x > xmax) xmax = pxyz[i].x;
  101.       if (pxyz[i].y < ymin) ymin = pxyz[i].y;
  102.       if (pxyz[i].y > ymax) ymax = pxyz[i].y;
  103.    }
  104.    dx = xmax - xmin;
  105.    dy = ymax - ymin;
  106.    dmax = (dx > dy) ? dx : dy;
  107.    xmid = (xmax + xmin) / 2.0;
  108.    ymid = (ymax + ymin) / 2.0;
  110.    /*
  111.       Set up the supertriangle
  112.       This is a triangle which encompasses all the sample points.
  113.       The supertriangle coordinates are added to the end of the
  114.       vertex list. The supertriangle is the first triangle in
  115.       the triangle list.
  116.    */
  117.    pxyz[nv+0].x = xmid - 20 * dmax;
  118.    pxyz[nv+0].y = ymid - dmax;
  119.    pxyz[nv+0].z = 0.0;
  120.    pxyz[nv+1].x = xmid;
  121.    pxyz[nv+1].y = ymid + 20 * dmax;
  122.    pxyz[nv+1].z = 0.0;
  123.    pxyz[nv+2].x = xmid + 20 * dmax;
  124.    pxyz[nv+2].y = ymid - dmax;
  125.    pxyz[nv+2].z = 0.0;
  126.    v[0].p1 = nv;
  127.    v[0].p2 = nv+1;
  128.    v[0].p3 = nv+2;
  129.    complete[0] = FALSE;
  130.    *ntri = 1;
  132.    /*
  133.       Include each point one at a time into the existing mesh
  134.    */
  135.    for (i=0;i<nv;i++) {
  137.       xp = pxyz[i].x;
  138.       yp = pxyz[i].y;
  139.       nedge = 0;
  141.       /*
  142.          Set up the edge buffer.
  143.          If the point (xp,yp) lies inside the circumcircle then the
  144.          three edges of that triangle are added to the edge buffer
  145.          and that triangle is removed.
  146.       */
  147.       for (j=0;j<(*ntri);j++) {
  148.          if (complete[j])
  149.             continue;
  150.          x1 = pxyz[v[j].p1].x;
  151.          y1 = pxyz[v[j].p1].y;
  152.          x2 = pxyz[v[j].p2].x;
  153.          y2 = pxyz[v[j].p2].y;
  154.          x3 = pxyz[v[j].p3].x;
  155.          y3 = pxyz[v[j].p3].y;
  156.          inside = CircumCircle(xp,yp,x1,y1,x2,y2,x3,y3,&xc,&yc,&r);
  157.          if (xc + r < xp)
  158.             complete[j] = TRUE;
  159.          if (inside) {
  160.             /* Check that we haven't exceeded the edge list size */
  161.             if (nedge+3 >= emax) {
  162.                emax += 100;
  163.                if ((edges = (IEDGE *) realloc(edges,emax*(long)sizeof(IEDGE))) == NULL) {
  164.                   status = 3;
  165.                   goto skip;
  166.                }
  167.             }
  168.             edges[nedge+0].p1 = v[j].p1;
  169.             edges[nedge+0].p2 = v[j].p2;
  170.             edges[nedge+1].p1 = v[j].p2;
  171.             edges[nedge+1].p2 = v[j].p3;
  172.             edges[nedge+2].p1 = v[j].p3;
  173.             edges[nedge+2].p2 = v[j].p1;
  174.             nedge += 3;
  175.             v[j] = v[(*ntri)-1];
  176.             complete[j] = complete[(*ntri)-1];
  177.             (*ntri)--;
  178.             j--;
  179.          }
  180.       }
  182.       /*
  183.          Tag multiple edges
  184.          Note: if all triangles are specified anticlockwise then all
  185.                interior edges are opposite pointing in direction.
  186.       */
  187.       for (j=0;j<nedge-1;j++) {
  188.          for (k=j+1;k<nedge;k++) {
  189.             if ((edges[j].p1 == edges[k].p2) && (edges[j].p2 == edges[k].p1)) {
  190.                edges[j].p1 = -1;
  191.                edges[j].p2 = -1;
  192.                edges[k].p1 = -1;
  193.                edges[k].p2 = -1;
  194.             }
  195.             /* Shouldn't need the following, see note above */
  196.             if ((edges[j].p1 == edges[k].p1) && (edges[j].p2 == edges[k].p2)) {
  197.                edges[j].p1 = -1;
  198.                edges[j].p2 = -1;
  199.                edges[k].p1 = -1;
  200.                edges[k].p2 = -1;
  201.             }
  202.          }
  203.       }
  205.       /*
  206.          Form new triangles for the current point
  207.          Skipping over any tagged edges.
  208.          All edges are arranged in clockwise order.
  209.       */
  210.       for (j=0;j<nedge;j++) {
  211.          if (edges[j].p1 < 0 || edges[j].p2 < 0)
  212.             continue;
  213.          if ((*ntri) >= trimax) {
  214.             status = 4;
  215.             goto skip;
  216.          }
  217.          v[*ntri].p1 = edges[j].p1;
  218.          v[*ntri].p2 = edges[j].p2;
  219.          v[*ntri].p3 = i;
  220.          complete[*ntri] = FALSE;
  221.          (*ntri)++;
  222.       }
  223.    }
  225.    /*
  226.       Remove triangles with supertriangle vertices
  227.       These are triangles which have a vertex number greater than nv
  228.    */
  229.    for (i=0;i<(*ntri);i++) {
  230.       if (v[i].p1 >= nv || v[i].p2 >= nv || v[i].p3 >= nv) {
  231.          v[i] = v[(*ntri)-1];
  232.          (*ntri)--;
  233.          i--;
  234.       }
  235.    }
  237. skip:
  238.    free(edges);
  239.    free(complete);
  240.    return(status);
  241. }
  244. /////////////////////////////////////////////////////////////////////
  245. /*
  246.    Return TRUE if a point (xp,yp) is inside the circumcircle made up
  247.    of the points (x1,y1), (x2,y2), (x3,y3)
  248.    The circumcircle centre is returned in (xc,yc) and the radius r
  249.    NOTE: A point on the edge is inside the circumcircle
  250. */
  251. int CTriangulate::CircumCircle(double xp,double yp,
  252.    double x1,double y1,double x2,double y2,double x3,double y3,
  253.    double *xc,double *yc,double *r)
  254. {
  255.    double m1,m2,mx1,mx2,my1,my2;
  256.    double dx,dy,rsqr,drsqr;
  258.    double EPSILON = 0.000001;
  260.    /* Check for coincident points */
  261.    if (fabs(y1-y2) < EPSILON && fabs(y2-y3) < EPSILON)
  262.        return(FALSE);
  264.    if (fabs(y2-y1) < EPSILON) {
  265.       m2 = - (x3-x2) / (y3-y2);
  266.       mx2 = (x2 + x3) / 2.0;
  267.       my2 = (y2 + y3) / 2.0;
  268.       *xc = (x2 + x1) / 2.0;
  269.       *yc = m2 * (*xc - mx2) + my2;
  270.    } else if (fabs(y3-y2) < EPSILON) {
  271.       m1 = - (x2-x1) / (y2-y1);
  272.       mx1 = (x1 + x2) / 2.0;
  273.       my1 = (y1 + y2) / 2.0;
  274.       *xc = (x3 + x2) / 2.0;
  275.       *yc = m1 * (*xc - mx1) + my1;
  276.    } else {
  277.       m1 = - (x2-x1) / (y2-y1);
  278.       m2 = - (x3-x2) / (y3-y2);
  279.       mx1 = (x1 + x2) / 2.0;
  280.       mx2 = (x2 + x3) / 2.0;
  281.       my1 = (y1 + y2) / 2.0;
  282.       my2 = (y2 + y3) / 2.0;
  283.       *xc = (m1 * mx1 - m2 * mx2 + my2 - my1) / (m1 - m2);
  284.       *yc = m1 * (*xc - mx1) + my1;
  285.    }
  287.    dx = x2 - *xc;
  288.    dy = y2 - *yc;
  289.    rsqr = dx*dx + dy*dy;
  290.    *r = sqrt(rsqr);
  292.    dx = xp - *xc;
  293.    dy = yp - *yc;
  294.    drsqr = dx*dx + dy*dy;
  296.    return((drsqr <= rsqr) ? TRUE : FALSE);
  297. }
  301. void CTriangulate::QSortVertexByX() //X增大排序
  302. {
  303. qsort((void *)pv_, (size_t)nv_, sizeof(XYZ), compare);
  304. }
  307. ///////////////////////////////////////////////////////////////////////
  308. /*
  309.    用成员变量m_szFileName指定的文件初始化pv_和nv_.
  310.    成功返回true，否则返回false
  311. */
  312. bool CTriangulate::Initial(char *szFileName) 
  313. {
  314. strcpy(szFileName_, szFileName); // 初始化成员变量，该文件保存了顶点信息
  316. char *buf = NULL;
  318. if ( this->szFileName_ == NULL ) return false;
  320. // 使用内存映像文件来读取文件数据
  321. HANDLE hFile = CreateFile( szFileName_,
  322.    GENERIC_READ,
  323.    FILE_SHARE_READ,
  324.    NULL, //security attributes
  325.    OPEN_EXISTING,
  326.    FILE_FLAG_SEQUENTIAL_SCAN,
  327.    NULL); //template file
  328. if ( hFile == INVALID_HANDLE_VALUE) return false; //因错误返回
  329. HANDLE hFileMapping = CreateFileMapping(hFile, 
  330. NULL,  // Security attributes
  331. PAGE_READONLY, 
  332. 0,   //内存映射的最大大小，这位64位中的高位
  333. 0, //64位中的低位，0为匹配文件的实际长度
  334. NULL); //给内存映射命名
  335. if ( hFileMapping == NULL ) return false;
  337. LPVOID pMem =  MapViewOfFile(hFileMapping,
  338.   FILE_MAP_READ,
  339.   0, //视图在映射文件的起始位置,64位中的高位
  340.   0, // 64位中的低位，这是从文件头开始映射
  341.   0); //映射长度
  342. if (pMem == NULL) return false;
  344. buf = (char *)pMem;
  347. ////////////////////////////////////////////////////////////////
  348. // TODO: Add your process in here
  349. //
  350. char linebuf[1024];
  351. char *pos;
  352. XYZ *pTemp;
  353. long MAXVERTEX = 1000;
  354. //pv_  = (XYZ *)malloc( MAXVERTEX * sizeof(XYZ) );
  355. if ((pv_ = new XYZ[MAXVERTEX]) == NULL )  return false;
  357. istrstream is(buf, GetFileSize(hFile, NULL));
  358. while ( ! is.eof() )
  359. {
  360. is.getline(linebuf, 1024);
  361. pos = strchr(linebuf, '#'); // 符号#后为注释
  362. if ( pos != NULL ) *pos = ''; //remove comment
  363. if ( sscanf(linebuf, "%le %le %le", &pv_[nv_].x, &pv_[nv_].y, &pv_[nv_].z ) != 3)
  364. continue; // error! continue to read next line
  365. else
  366. {
  367. nv_++;
  368. if ( nv_ >= MAXVERTEX)
  369. {
  370. if ((pTemp = new XYZ[MAXVERTEX + 1000]) == NULL ) return false;
  371. memcpy(pTemp, pv_, sizeof(XYZ) * MAXVERTEX);
  372. delete []pv_;
  373. pv_ = pTemp;
  374. MAXVERTEX += 1000;
  375. }
  376. }
  377. }
  378. //
  379. // end your process
  380. //////////////////////////////////////////////////////////////////
  383. CloseHandle(hFile);
  384. CloseHandle(hFileMapping);
  385. UnmapViewOfFile(pMem);
  388. // set RBvertex_ and LTvertex_
  389. LTvertex_ = pv_[0];
  390. RBvertex_ = LTvertex_;
  391. for ( int i = 1; i < nv_; i++ )
  392. {
  393. if ( pv_[i].x > LTvertex_.x ) LTvertex_.x = pv_[i].x;
  394. if ( pv_[i].x < RBvertex_.x ) RBvertex_.x = pv_[i].x;
  395. if ( pv_[i].y > LTvertex_.y ) LTvertex_.y = pv_[i].y;
  396. if ( pv_[i].y < RBvertex_.y ) RBvertex_.y = pv_[i].y;
  397. if ( pv_[i].z > LTvertex_.z ) LTvertex_.z = pv_[i].z;
  398. if ( pv_[i].z < RBvertex_.z ) RBvertex_.z = pv_[i].z;
  399. }
  400. return true; //正常返回为0
  401. }
  404. //////////////////////////////////////////////////////////////////////
  405. /*
  406.    根据nv_和pv_创建三角形，创建的三角形保存在ptri_中
  407.    三角形的个数保存在ntri_中；
  408.    其中需要调用三角形的剖分算法Triangulate()
  410.    注意：
  411.        有很多种三角剖分算法，Triangulate()是一种增量算法。
  412. */
  413. bool CTriangulate::CreateTriangle()
  414. {
  415. SortAndRemoveRepeatPoints(); // 根据x的大小对pv_进行排序,并且删掉重复的点
  416. // Triangulate()函数需要已经排序过的pv_
  417. AfxMessageBox("排序和删除重复点结束。",0,0);
  419. if ((ptri_ = new TRIANGLE[4 * nv_]) == NULL)
  420. return false; //三角形的最多个数是顶点个数的4倍
  421. Triangulate(nv_, pv_, &ntri_, ptri_);
  423. AfxMessageBox("三角剖分结束。",0,0);
  425. MeshedTris();
  426. AfxMessageBox("建立三角形之间的关系结束。",0,0);
  428. return true;
  429. }
  435. XYZ CTriangulate::GetLTvertex() // 获得散点的xyz的最大值
  436. {
  437. return LTvertex_;
  438. }
  440. XYZ CTriangulate::GetRBvertex() // 获得散点的xyz的最小值
  441. {
  442. return RBvertex_;
  443. }
  447. // 为TRIANGLE的T1，T2，T3设定值。即建立三角形的邻接关系
  448. void CTriangulate::MeshedTris() //对各个独立的三角形建立网络联系
  449. {
  450. int i, j;
  452. for (i = 0; i < ntri_; i++) // remove 不规范的点, need it?
  453. {
  454. if (ptri_[i].p1 == ptri_[i].p2 || ptri_[i].p1 == ptri_[i].p3 || ptri_[i].p2 == ptri_[i].p3 )
  455. {
  456. ptri_[i] = ptri_[ntri_ - 1];
  457. ntri_--;
  458. }
  459. }
  461. for (i = 0; i < ntri_; i++)
  462. {
  463. if (ptri_[i].t3 == NEIGHBOR_UNKNOWN)
  464. for ( j = i + 1; j < ntri_; j++)
  465. {
  466. if (   (ptri_[i].p1 == ptri_[j].p1 || ptri_[i].p1 == ptri_[j].p2 || ptri_[i].p1 == ptri_[j].p3 )
  467. && (ptri_[i].p2 == ptri_[j].p1 || ptri_[i].p2 == ptri_[j].p2 || ptri_[i].p2 == ptri_[j].p3 ) )
  468. {
  469. ptri_[i].t3 = j; // 设置第i个三角形的边p1p2的邻接三角形
  471. // 由于两个邻接三角形是相互邻接的，所以同时设置“邻接三角形j”的邻接三角形为i
  472. if ( ptri_[j].p1 + ptri_[j].p2 == ptri_[i].p1 + ptri_[i].p2 ) 
  473. ptri_[j].t3 = i; 
  474. else if ( ptri_[j].p1 + ptri_[j].p3 == ptri_[i].p1 + ptri_[i].p2 ) 
  475. ptri_[j].t2 = i;
  476. else //( ptri_[j].p3 + ptri_[j].p2 == ptri_[i].p1 + ptri_[i].p2 ) 
  477. ptri_[j].t1 = i;
  479. break;
  480. }
  481. }
  482. if (ptri_[i].t2 == NEIGHBOR_UNKNOWN)
  483. for ( j = i + 1; j < ntri_; j++)
  484. {
  485. if (   (ptri_[i].p1 == ptri_[j].p1 || ptri_[i].p1 == ptri_[j].p2 || ptri_[i].p1 == ptri_[j].p3 )
  486. && (ptri_[i].p3 == ptri_[j].p1 || ptri_[i].p3 == ptri_[j].p2 || ptri_[i].p3 == ptri_[j].p3 ) )
  487. {
  488. ptri_[i].t2 = j; // 设置第i个三角形的边p1p2的邻接三角形
  490. // 由于两个邻接三角形是相互邻接的，所以同时设置“邻接三角形j”的邻接三角形为i
  491. if ( ptri_[j].p1 + ptri_[j].p2 == ptri_[i].p1 + ptri_[i].p3 ) 
  492. ptri_[j].t3 = i; 
  493. else if ( ptri_[j].p1 + ptri_[j].p3 == ptri_[i].p1 + ptri_[i].p3 ) 
  494. ptri_[j].t2 = i;
  495. else //if ( ptri_[j].p3 + ptri_[j].p2 == ptri_[i].p1 + ptri_[i].p3 ) 
  496. ptri_[j].t1 = i;
  498. break;
  499. }
  500. }
  501. if (ptri_[i].t1 == NEIGHBOR_UNKNOWN)
  502. for ( j = i + 1; j < ntri_; j++)
  503. {
  504. if (   (ptri_[i].p3 == ptri_[j].p1 || ptri_[i].p3 == ptri_[j].p2 || ptri_[i].p3 == ptri_[j].p3 )
  505. && (ptri_[i].p2 == ptri_[j].p1 || ptri_[i].p2 == ptri_[j].p2 || ptri_[i].p2 == ptri_[j].p3 ) )
  506. {
  507. ptri_[i].t1 = j; // 设置第i个三角形的边p1p2的邻接三角形
  509. // 由于两个邻接三角形是相互邻接的，所以同时设置“邻接三角形j”的邻接三角形为i
  510. if ( ptri_[j].p1 + ptri_[j].p2 == ptri_[i].p3 + ptri_[i].p2 ) 
  511. ptri_[j].t3 = i; 
  512. else if ( ptri_[j].p1 + ptri_[j].p3 == ptri_[i].p3 + ptri_[i].p2 ) 
  513. ptri_[j].t2 = i;
  514. else //if ( ptri_[j].p3 + ptri_[j].p2 == ptri_[i].p3 + ptri_[i].p2 ) 
  515. ptri_[j].t1 = i;
  517. break;
  518. }
  519. }
  520. } // end i
  521. }
  523. void CTriangulate::SortAndRemoveRepeatPoints()
  524. {
  525. int i, j;
  526. int pos(0), pos2(0);
  527. int readpos(0), writepos(0);
  529. QSortVertexByX(); // 按X增大排序
  531. bool *flagkilled = new bool[nv_];
  532. for ( i = 0; i < nv_; i++)
  533. flagkilled[i] = false;
  535. while( pos < nv_)
  536. {
  537. while( pv_[++pos2].x == pv_[pos].x );
  539. for(i = pos; i < pos2; i++)
  540. {
  541. if(flagkilled[i]) continue;
  542. for (j = i + 1; j < pos2; j++)
  543. if(pv_[j].y == pv_[i].y) 
  544. flagkilled[j] = true;
  545. }
  546. pos = pos2;
  547. }
  549. while(readpos < nv_)
  550. {
  551. while(flagkilled[readpos]) // skip killed points
  552. readpos++;
  553. if (writepos != readpos)
  554. pv_[writepos] = pv_[readpos];
  555. writepos++;
  556. readpos++;
  557. }
  559. nv_ = writepos;
  560. delete []flagkilled;
  562. //TRACE("output norepeat points: nv=%dn", nv_);
  563. //for (i=0; i<nv_; i++)
  564. // TRACE("%lf  %lf  %lfn", pv_[i].x, pv_[i].y, pv_[i].z);
  565. //TRACE("End outputn");
  566. }