Delaunay.cpp
上传用户:shlanyl88
上传日期:2013-03-14
资源大小:147k
文件大小:10k
源码类别:

3D图形编程

开发平台:

Visual C++

  1. /********************************************************************************
  2. Copyright (C) 2004-2005 Sjaak Priester
  3. This is free software; you can redistribute it and/or modify
  4. it under the terms of the GNU General Public License as published by
  5. the Free Software Foundation; either version 2 of the License, or
  6. (at your option) any later version.
  7. This file is distributed in the hope that it will be useful,
  8. but WITHOUT ANY WARRANTY; without even the implied warranty of
  9. MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  10. GNU General Public License for more details.
  11. You should have received a copy of the GNU General Public License
  12. along with Tinter; if not, write to the Free Software
  13. Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
  14. ********************************************************************************/
  15. // Delaunay
  16. // Class to perform Delaunay triangulation on a set of vertices
  17. //
  18. // Version 1.1 (C) 2005, Sjaak Priester, Amsterdam.
  19. // - Removed bug which gave incorrect results for co-circular vertices.
  20. //
  21. // Version 1.0 (C) 2004, Sjaak Priester, Amsterdam.
  22. // mailto:sjaak@sjaakpriester.nl
  23. #include "StdAfx.h"
  24. #include "Delaunay.h"
  25. const REAL sqrt3 = 1.732050808F;
  26. void triangle::SetCircumCircle()
  27. {
  28. REAL x0 = m_Vertices[0]->GetX();
  29. REAL y0 = m_Vertices[0]->GetY();
  30. REAL x1 = m_Vertices[1]->GetX();
  31. REAL y1 = m_Vertices[1]->GetY();
  32. REAL x2 = m_Vertices[2]->GetX();
  33. REAL y2 = m_Vertices[2]->GetY();
  34. REAL y10 = y1 - y0;
  35. REAL y21 = y2 - y1;
  36. bool b21zero = y21 > -REAL_EPSILON && y21 < REAL_EPSILON;
  37. if (y10 > -REAL_EPSILON && y10 < REAL_EPSILON)
  38. {
  39. if (b21zero) // All three vertices are on one horizontal line.
  40. {
  41. if (x1 > x0)
  42. {
  43. if (x2 > x1) x1 = x2;
  44. }
  45. else
  46. {
  47. if (x2 < x0) x0 = x2;
  48. }
  49. m_Center.X = (x0 + x1) * .5F;
  50. m_Center.Y = y0;
  51. }
  52. else // m_Vertices[0] and m_Vertices[1] are on one horizontal line.
  53. {
  54. REAL m1 = - (x2 - x1) / y21;
  55. REAL mx1 = (x1 + x2) * .5F;
  56. REAL my1 = (y1 + y2) * .5F;
  57. m_Center.X = (x0 + x1) * .5F;
  58. m_Center.Y = m1 * (m_Center.X - mx1) + my1;
  59. }
  60. }
  61. else if (b21zero) // m_Vertices[1] and m_Vertices[2] are on one horizontal line.
  62. {
  63. REAL m0 = - (x1 - x0) / y10;
  64. REAL mx0 = (x0 + x1) * .5F;
  65. REAL my0 = (y0 + y1) * .5F;
  66. m_Center.X = (x1 + x2) * .5F;
  67. m_Center.Y = m0 * (m_Center.X - mx0) + my0;
  68. }
  69. else // 'Common' cases, no multiple vertices are on one horizontal line.
  70. {
  71. REAL m0 = - (x1 - x0) / y10;
  72. REAL m1 = - (x2 - x1) / y21;
  73. REAL mx0 = (x0 + x1) * .5F;
  74. REAL my0 = (y0 + y1) * .5F;
  75. REAL mx1 = (x1 + x2) * .5F;
  76. REAL my1 = (y1 + y2) * .5F;
  77. m_Center.X = (m0 * mx0 - m1 * mx1 + my1 - my0) / (m0 - m1);
  78. m_Center.Y = m0 * (m_Center.X - mx0) + my0;
  79. }
  80. REAL dx = x0 - m_Center.X;
  81. REAL dy = y0 - m_Center.Y;
  82. m_R2 = dx * dx + dy * dy; // the radius of the circumcircle, squared
  83. m_R = (REAL) sqrt(m_R2); // the proper radius
  84. // Version 1.1: make m_R2 slightly higher to ensure that all edges
  85. // of co-circular vertices will be caught.
  86. // Note that this is a compromise. In fact, the algorithm isn't really
  87. // suited for very many co-circular vertices.
  88. m_R2 *= 1.000001f;
  89. }
  90. // Function object to check whether a triangle has one of the vertices in SuperTriangle.
  91. // operator() returns true if it does.
  92. class triangleHasVertex
  93. {
  94. public:
  95. triangleHasVertex(const vertex SuperTriangle[3]) : m_pSuperTriangle(SuperTriangle) {}
  96. bool operator()(const triangle& tri) const
  97. {
  98. for (int i = 0; i < 3; i++)
  99. {
  100. const vertex * p = tri.GetVertex(i);
  101. if (p >= m_pSuperTriangle && p < (m_pSuperTriangle + 3)) return true;
  102. }
  103. return false;
  104. }
  105. protected:
  106. const vertex * m_pSuperTriangle;
  107. };
  108. // Function object to check whether a triangle is 'completed', i.e. doesn't need to be checked
  109. // again in the algorithm, i.e. it won't be changed anymore.
  110. // Therefore it can be removed from the workset.
  111. // A triangle is completed if the circumcircle is completely to the left of the current vertex.
  112. // If a triangle is completed, it will be inserted in the output set, unless one or more of it's vertices
  113. // belong to the 'super triangle'.
  114. class triangleIsCompleted
  115. {
  116. public:
  117. triangleIsCompleted(cvIterator itVertex, triangleSet& output, const vertex SuperTriangle[3])
  118. : m_itVertex(itVertex)
  119. , m_Output(output)
  120. , m_pSuperTriangle(SuperTriangle)
  121. {}
  122. bool operator()(const triangle& tri) const
  123. {
  124. bool b = tri.IsLeftOf(m_itVertex);
  125. if (b)
  126. {
  127. triangleHasVertex thv(m_pSuperTriangle);
  128. if (! thv(tri)) m_Output.insert(tri);
  129. }
  130. return b;
  131. }
  132. protected:
  133. cvIterator m_itVertex;
  134. triangleSet& m_Output;
  135. const vertex * m_pSuperTriangle;
  136. };
  137. // Function object to check whether vertex is in circumcircle of triangle.
  138. // operator() returns true if it does.
  139. // The edges of a 'hot' triangle are stored in the edgeSet edges.
  140. class vertexIsInCircumCircle
  141. {
  142. public:
  143. vertexIsInCircumCircle(cvIterator itVertex, edgeSet& edges) : m_itVertex(itVertex), m_Edges(edges) {}
  144. bool operator()(const triangle& tri) const
  145. {
  146. bool b = tri.CCEncompasses(m_itVertex);
  147. if (b)
  148. {
  149. HandleEdge(tri.GetVertex(0), tri.GetVertex(1));
  150. HandleEdge(tri.GetVertex(1), tri.GetVertex(2));
  151. HandleEdge(tri.GetVertex(2), tri.GetVertex(0));
  152. }
  153. return b;
  154. }
  155. protected:
  156. void HandleEdge(const vertex * p0, const vertex * p1) const
  157. {
  158. const vertex * pVertex0(NULL);
  159. const vertex * pVertex1(NULL);
  160. // Create a normalized edge, in which the smallest vertex comes first.
  161. if (* p0 < * p1)
  162. {
  163. pVertex0 = p0;
  164. pVertex1 = p1;
  165. }
  166. else
  167. {
  168. pVertex0 = p1;
  169. pVertex1 = p0;
  170. }
  171. edge e(pVertex0, pVertex1);
  172. // Check if this edge is already in the buffer
  173. edgeIterator found = m_Edges.find(e);
  174. if (found == m_Edges.end()) m_Edges.insert(e); // no, it isn't, so insert
  175. else m_Edges.erase(found); // yes, it is, so erase it to eliminate double edges
  176. }
  177. cvIterator m_itVertex;
  178. edgeSet& m_Edges;
  179. };
  180. void Delaunay::Triangulate(const vertexSet& vertices, triangleSet& output)
  181. {
  182. if (vertices.size() < 3) return; // nothing to handle
  183. // Determine the bounding box.
  184. cvIterator itVertex = vertices.begin();
  185. REAL xMin = itVertex->GetX();
  186. REAL yMin = itVertex->GetY();
  187. REAL xMax = xMin;
  188. REAL yMax = yMin;
  189. ++itVertex; // If we're here, we know that vertices is not empty.
  190. for (; itVertex != vertices.end(); itVertex++)
  191. {
  192. xMax = itVertex->GetX(); // Vertices are sorted along the x-axis, so the last one stored will be the biggest.
  193. REAL y = itVertex->GetY();
  194. if (y < yMin) yMin = y;
  195. if (y > yMax) yMax = y;
  196. }
  197. REAL dx = xMax - xMin;
  198. REAL dy = yMax - yMin;
  199. // Make the bounding box slightly bigger, just to feel safe.
  200. REAL ddx = dx * 0.01F;
  201. REAL ddy = dy * 0.01F;
  202. xMin -= ddx;
  203. xMax += ddx;
  204. dx += 2 * ddx;
  205. yMin -= ddy;
  206. yMax += ddy;
  207. dy += 2 * ddy;
  208. // Create a 'super triangle', encompassing all the vertices. We choose an equilateral triangle with horizontal base.
  209. // We could have made the 'super triangle' simply very big. However, the algorithm is quite sensitive to
  210. // rounding errors, so it's better to make the 'super triangle' just big enough, like we do here.
  211. vertex vSuper[3];
  212. vSuper[0] = vertex(xMin - dy * sqrt3 / 3.0F, yMin); // Simple highschool geometry, believe me.
  213. vSuper[1] = vertex(xMax + dy * sqrt3 / 3.0F, yMin);
  214. vSuper[2] = vertex((xMin + xMax) * 0.5F, yMax + dx * sqrt3 * 0.5F);
  215. triangleSet workset;
  216. workset.insert(triangle(vSuper));
  217. for (itVertex = vertices.begin(); itVertex != vertices.end(); itVertex++)
  218. {
  219. // First, remove all 'completed' triangles from the workset.
  220. // A triangle is 'completed' if its circumcircle is entirely to the left of the current vertex.
  221. // Vertices are sorted in x-direction (the set container does this automagically).
  222. // Unless they are part of the 'super triangle', copy the 'completed' triangles to the output.
  223. // The algorithm also works without this step, but it is an important optimalization for bigger numbers of vertices.
  224. // It makes the algorithm about five times faster for 2000 vertices, and for 10000 vertices,
  225. // it's thirty times faster. For smaller numbers, the difference is negligible.
  226. tIterator itEnd = remove_if(workset.begin(), workset.end(), triangleIsCompleted(itVertex, output, vSuper));
  227. edgeSet edges;
  228. // A triangle is 'hot' if the current vertex v is inside the circumcircle.
  229. // Remove all hot triangles, but keep their edges.
  230. itEnd = remove_if(workset.begin(), itEnd, vertexIsInCircumCircle(itVertex, edges));
  231. workset.erase(itEnd, workset.end()); // remove_if doesn't actually remove; we have to do this explicitly.
  232. // Create new triangles from the edges and the current vertex.
  233. for (edgeIterator it = edges.begin(); it != edges.end(); it++)
  234. workset.insert(triangle(it->m_pV0, it->m_pV1, & (* itVertex)));
  235. }
  236. // Finally, remove all the triangles belonging to the 'super triangle' and move the remaining
  237. // triangles tot the output; remove_copy_if lets us do that in one go.
  238. tIterator where = output.begin();
  239. remove_copy_if(workset.begin(), workset.end(), inserter(output, where), triangleHasVertex(vSuper));
  240. }
  241. void Delaunay::TrianglesToEdges(const triangleSet& triangles, edgeSet& edges)
  242. {
  243. for (ctIterator it = triangles.begin(); it != triangles.end(); ++it)
  244. {
  245. HandleEdge(it->GetVertex(0), it->GetVertex(1), edges);
  246. HandleEdge(it->GetVertex(1), it->GetVertex(2), edges);
  247. HandleEdge(it->GetVertex(2), it->GetVertex(0), edges);
  248. }
  249. }
  250. void Delaunay::HandleEdge(const vertex * p0, const vertex * p1, edgeSet& edges)
  251. {
  252. const vertex * pV0(NULL);
  253. const vertex * pV1(NULL);
  254. if (* p0 < * p1)
  255. {
  256. pV0 = p0;
  257. pV1 = p1;
  258. }
  259. else
  260. {
  261. pV0 = p1;
  262. pV1 = p0;
  263. }
  264. // Insert a normalized edge. If it's already in edges, insertion will fail,
  265. // thus leaving only unique edges.
  266. edges.insert(edge(pV0, pV1));
  267. }