Developer Documentation
Triangulator.cc
1 /*===========================================================================*\
2 * *
3 * OpenFlipper *
4 * Copyright (c) 2001-2015, RWTH-Aachen University *
5 * Department of Computer Graphics and Multimedia *
6 * All rights reserved. *
7 * www.openflipper.org *
8 * *
9 *---------------------------------------------------------------------------*
10 * This file is part of OpenFlipper. *
11 *---------------------------------------------------------------------------*
12 * *
13 * Redistribution and use in source and binary forms, with or without *
14 * modification, are permitted provided that the following conditions *
15 * are met: *
16 * *
17 * 1. Redistributions of source code must retain the above copyright notice, *
18 * this list of conditions and the following disclaimer. *
19 * *
20 * 2. Redistributions in binary form must reproduce the above copyright *
21 * notice, this list of conditions and the following disclaimer in the *
22 * documentation and/or other materials provided with the distribution. *
23 * *
24 * 3. Neither the name of the copyright holder nor the names of its *
25 * contributors may be used to endorse or promote products derived from *
26 * this software without specific prior written permission. *
27 * *
28 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS *
29 * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED *
30 * TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A *
31 * PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER *
32 * OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, *
33 * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, *
34 * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR *
35 * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF *
36 * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING *
37 * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS *
38 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. *
39 * *
40 \*===========================================================================*/
41 
42 /*===========================================================================*\
43 * *
44 * $Revision: 21022 $ *
45 * $Author: moebius $ *
46 * $Date: 2015-07-17 08:23:03 +0200 (Fri, 17 Jul 2015) $ *
47 * *
48 \*===========================================================================*/
49 
50 
51 #include "Triangulator.hh"
52 
53 
54 #include <iostream>
55 
56 
57 
58 namespace ACG {
59 
60 
61 Triangulator::Triangulator(const std::vector<Vec3f>& _pos)
62  : polySize_(_pos.size()), numRemaningVertices_(_pos.size()), numTris_(0),
63  numReflexVertices_(0),
64  status_(-1), convex_(false)
65 {
66  if (polySize_ < 3)
67  return;
68 
69 
70  if (polySize_ == 3)
71  {
72  numTris_ = 1;
73  tris_.resize(3);
74  tris_[0] = 0;
75  tris_[1] = 1;
76  tris_[2] = 2;
77 
78  numRemaningVertices_ = 0;
79  convex_ = true;
80  status_ = 0;
81  }
82  else
83  {
84  // project vertices onto the 2d plane of the polygon.
85  // the projection plane is chosen orthogonal to the polygon surface normal
86 
87  // use Newell's Method to compute the surface normal
88  // http://www.opengl.org/wiki/Calculating_a_Surface_Normal
89 
90  Vec3f n(0.0f, 0.0f, 0.0f);
91 
92  for (int i = 0; i < polySize_; ++i)
93  {
94  int next = (i + 1) % polySize_;
95 
96  Vec3f a = _pos[i] - _pos[next];
97  Vec3f b = _pos[i] + _pos[next];
98 
99  n[0] += a[1] * b[2];
100  n[1] += a[2] * b[0];
101  n[2] += a[0] * b[1];
102  }
103 
104  // project to 2d
105  pos_.resize(polySize_);
106 
107  Vec3f axis[3] = { Vec3f(1.0f, 0.0f, 0.0f), Vec3f(0.0f, 1.0f, 0.0f), n };
108 
109  // orthonormalize projection axes
110  axis[2].normalize();
111 
112  // make sure first axis is linearly independent from the normal
113  while (std::abs(axis[0] | axis[2]) > 0.95f || (axis[0].sqrnorm() < 0.001f))
114  {
115  for (int i = 0; i < 3; ++i)
116  axis[0][i] = float(rand()) / float(RAND_MAX) * 2.0f - 1.0f;
117 
118  axis[0].normalize();
119  }
120 
121  // make axis[0] orthogonal to normal
122  axis[0] = axis[0] - axis[2] * (axis[0] | axis[2]);
123  axis[0].normalize();
124  axis[1] = axis[2] % axis[0];
125 
126 
127  for (int i = 0; i < polySize_; ++i)
128  {
129  // project onto polygon plane
130  pos_[i][0] = axis[0] | _pos[i];
131  pos_[i][1] = axis[1] | _pos[i];
132  }
133 
134 
135  // create triangle fans if there is at most one concave vertex
136  int reflexVertexID = 0;
137 
138  for (int i = 0; i < polySize_; ++i)
139  {
140  // test vertex (i+1)
141  if (isReflexVertex(pos_[i], pos_[(i + 1) % polySize_], pos_[(i + 2) % polySize_]))
142  {
143  ++numReflexVertices_;
144  reflexVertexID = (i + 1) % polySize_;
145  }
146  }
147 
148  convex_ = !numReflexVertices_;
149 
150 
151  if (numReflexVertices_ <= 1)
152  {
153  // create triangle fans
154  numTris_ = polySize_ - 2;
155  tris_.resize(numTris_ * 3);
156  numRemaningVertices_ = 0;
157  status_ = 0;
158 
159  for (int i = 0; i < numTris_; ++i)
160  {
161  tris_[i * 3] = reflexVertexID;
162  tris_[i * 3 + 1] = (reflexVertexID + i + 1) % polySize_;
163  tris_[i * 3 + 2] = (reflexVertexID + i + 2) % polySize_;
164  }
165 
166  }
167  else
168  {
169  // use the ear clipping algorithm
170 
171  earClippingN2();
172 // earClippingN3();
173 
174 // triangulateExternal();
175  }
176  }
177 }
178 
179 
181 {
182 }
183 
184 
185 bool Triangulator::isReflexVertex(const Vec2f& v0, const Vec2f& v1, const Vec2f& v2) const
186 {
187  // compute the sign of the cross product of the edges sharing vertex v1
188  // <0 : inner angle greater than 180 deg (reflex)
189  // >0 : inner angle less than 180 deg (convex)
190 
191  Vec2f u = v2 - v1;
192  Vec2f v = v0 - v1;
193 
194  return u[0] * v[1] - u[1] * v[0] < 0.0f;
195 }
196 
198 {
199  int p = (i + polySize_ - 1) % polySize_;
200  int n = (i + 1) % polySize_;
201 
202  return isReflexVertex(pos_[p], pos_[i], pos_[n]);
203 }
204 
205 float Triangulator::triangleAreaSign(const Vec2f& v0, const Vec2f& v1, const Vec2f& v2) const
206 {
207  // cross product
208  return (v0[0] - v2[0]) * (v1[1] - v2[1]) - (v1[0] - v2[0]) * (v0[1] - v2[1]);
209 }
210 
211 float Triangulator::distancePointToSegmentSq(const Vec2f& v0, const Vec2f& v1, const Vec2f& pt) const
212 {
213  Vec2f segment = v1 - v0;
214  Vec2f vec = pt - v0;
215  float dp = vec | segment;
216 
217  if (dp < 0.0f)
218  return vec.sqrnorm();
219 
220  float segSq = segment.sqrnorm();
221  dp /= segSq;
222 
223  if (dp < 1.0f)
224  return vec.sqrnorm() - dp * dp * segSq;
225 
226  vec = pt - v1;
227  return vec.sqrnorm();
228 }
229 
230 bool Triangulator::pointInTriangle(const Vec2f& v0, const Vec2f& v1, const Vec2f& v2, const Vec2f& pt) const
231 {
232  // ACG implementation based on barycentric coordinates (slow)
233 // return Geometry::isInTriangle(pt, v0, v1, v2);
234 
235 
236  // fast implementation based on triangle areas
237  // http://www.gamedev.net/topic/295943-is-this-a-better-point-in-triangle-test-2d/
238 
239  return triangleAreaSign(pt, v0, v1) >= 0.0f &&
240  triangleAreaSign(pt, v1, v2) >= 0.0f &&
241  triangleAreaSign(pt, v2, v0) >= 0.0f;
242 
243 
244 // // more accurate algorithm:
245 // // http://totologic.blogspot.de/2014/01/accurate-point-in-triangle-test.html
246 // // note: didn't improve accuracy at all for problematic polygons
247 //
248 // static const float eps = 1e-4f;
249 // static const float eps2 = eps*eps;
250 //
251 //
252 // // point in aabb of triangle
253 // Vec2f aabbMin = v0;
254 //
255 // aabbMin.minimize(v1);
256 // aabbMin.minimize(v2);
257 //
258 // if (pt[0] + eps < aabbMin[0] || pt[1] + eps < aabbMin[1])
259 // return false;
260 //
261 // Vec2f aabbMax = v0;
262 // aabbMax.maximize(v1);
263 // aabbMax.maximize(v2);
264 //
265 // if (pt[0] > aabbMax[0] + eps || pt[1] > aabbMax[1] + eps)
266 // return false;
267 //
268 //
269 // if (triangleAreaSign(pt, v0, v1) >= 0.0f &&
270 // triangleAreaSign(pt, v1, v2) >= 0.0f &&
271 // triangleAreaSign(pt, v2, v0) >= 0.0f)
272 // return true;
273 //
274 //
275 // return (distancePointToSegmentSq(v0, v1, pt) <= eps2 ||
276 // distancePointToSegmentSq(v1, v2, pt) <= eps2 ||
277 // distancePointToSegmentSq(v2, v0, pt) <= eps2);
278 }
279 
280 
281 void Triangulator::initVertexList()
282 {
283  vertices_.resize(polySize_);
284  reflexVertices_.clear();
285  for (int i = 0; i < polySize_; ++i)
286  {
287  int p = (i + polySize_ - 1) % polySize_;
288  int n = (i + 1) % polySize_;
289 
290  vertices_[i] = RingVertex(i, isReflexVertex(pos_[p], pos_[i], pos_[n]), pos_[i], &vertices_[p], &vertices_[n]);
291 
292  if (vertices_[i].reflex)
293  reflexVertices_.push_back(&vertices_[i]);
294  }
295 }
296 
297 
298 
299 int Triangulator::earClippingN3()
300 {
301  // http://www.geometrictools.com/Documentation/TriangulationByEarClipping.pdf
302  // O(n^3)
303 
304  numTris_ = 0;
305  status_ = 0;
306 
307  initVertexList();
308 
309  tris_.resize((polySize_ - 2) * 3);
310 
311  int numIterations = 0;
312  RingVertex* firstVertex = &vertices_[0];
313 
314  while (numTris_ < polySize_ - 2)
315  {
316  // find an ear in the remaining polygon
317  bool hasEars = false;
318 
319  RingVertex* curVertex = firstVertex;
320  do
321  {
322  curVertex->reflex = isReflexVertex(curVertex->prev->pos, curVertex->pos, curVertex->next->pos);
323 
324  if (!curVertex->reflex)
325  {
326  // triangle containment test
327  bool isEar = true;
328 
329  // check all remaining vertices r for containment
330  for (RingVertex* r = curVertex->next->next; r != curVertex->prev; r = r->next)
331  {
332  if (pointInTriangle(curVertex->prev->pos, curVertex->pos, curVertex->next->pos, r->pos))
333  {
334  isEar = false;
335  break;
336  }
337  }
338 
339  // found an ear
340  if (isEar)
341  {
342  // triangulate ear
343  hasEars = true;
344  tris_[numTris_ * 3] = curVertex->prev->id;
345  tris_[numTris_ * 3 + 1] = curVertex->id;
346  tris_[numTris_ * 3 + 2] = curVertex->next->id;
347  ++numTris_;
348 
349  // remove vertex from linked list
350  curVertex->prev->next = curVertex->next;
351  curVertex->next->prev = curVertex->prev;
352 
353  break;
354  }
355  }
356 
357  curVertex = curVertex->next;
358  ++numIterations;
359  } while (curVertex != firstVertex);
360 
361  firstVertex = firstVertex->next;
362 
363  // create triangle fans and hope for good result if there are no more ears
364  if (!hasEars && (numTris_ + 2 < polySize_))
365  {
366  for (RingVertex* iteratorVertex = firstVertex->next; iteratorVertex != firstVertex->prev; iteratorVertex = iteratorVertex->next)
367  {
368  tris_[numTris_ * 3] = firstVertex->id;
369  tris_[numTris_ * 3 + 1] = iteratorVertex->id;
370  tris_[numTris_ * 3 + 2] = iteratorVertex->next->id;
371  ++numTris_;
372  }
373 
374  assert(numTris_ == polySize_ - 2);
375  }
376  }
377 
378  return numTris_;
379 }
380 
381 
382 
383 int Triangulator::earClippingN2()
384 {
385  // http://www.geometrictools.com/Documentation/TriangulationByEarClipping.pdf
386  // O(n^2)
387 
388  numTris_ = 0;
389  status_ = 0;
390 
391  initVertexList();
392 
393  // triangulate
394 
395  int numTries = 0; // # checked vertices per iteration that aren't ears
396  numRemaningVertices_ = polySize_; // size of currently remaining polygon
397  int numIterations = 0;
398 
399  tris_.resize((polySize_ - 2) * 3);
400 
401  RingVertex* curVertex = &vertices_[0];
402  while (numRemaningVertices_ > 3)
403  {
404  // check if the current vertex is an ear tip
405  bool isEar = false;
406 
407  if (!curVertex->reflex)
408  {
409  // test current vertex for ear property
410  isEar = true;
411  for (std::list<RingVertex*>::iterator it = reflexVertices_.begin(); isEar && it != reflexVertices_.end(); ++it)
412  {
413  // skip direct neighbors
414  if (*it == curVertex->prev || *it == curVertex->next)
415  continue;
416 
417  // if any remaining vertex is inside the triangle, the current vertex is not an ear
418  if (pointInTriangle(curVertex->prev->pos, curVertex->pos, curVertex->next->pos, (*it)->pos))
419  isEar = false;
420  }
421 
422 
423  // found an ear
424  if (isEar)
425  {
426  addEar(curVertex);
427  numTries = 0;
428  }
429  }
430 
431  if (!isEar)
432  ++numTries;
433 
434  if (numTries > numRemaningVertices_)
435  {
436  // something went wrong
437  // create a triangle anyway and hope the result is ok
438 
439  addEar(curVertex);
440  numTries = 0;
441 
442  status_ = 1;
443  }
444 
445  curVertex = curVertex->next;
446  ++numIterations;
447  }
448 
449 
450  // add the last remaining triangle
451  if (numRemaningVertices_ == 3)
452  {
453  tris_[numTris_ * 3 + 0] = curVertex->prev->id;
454  tris_[numTris_ * 3 + 1] = curVertex->id;
455  tris_[numTris_ * 3 + 2] = curVertex->next->id;
456  ++numTris_;
457  }
458 
459 
460  return numTris_;
461 }
462 
463 bool Triangulator::updateReflexVertex(RingVertex* v)
464 {
465  if (v->reflex)
466  {
467  // check reflex property
468  v->reflex = isReflexVertex(v->prev->pos, v->pos, v->next->pos);
469 
470  // update list of reflex vertices
471  if (!v->reflex)
472  reflexVertices_.remove(v);
473  }
474 
475  return v->reflex;
476 }
477 
478 void Triangulator::addEar(RingVertex* _earTip)
479 {
480  // add ear triangle
481  tris_[numTris_ * 3 + 0] = _earTip->prev->id;
482  tris_[numTris_ * 3 + 1] = _earTip->id;
483  tris_[numTris_ * 3 + 2] = _earTip->next->id;
484 
485  // remove ear vertex from the linked list
486  _earTip->prev->next = _earTip->next;
487  _earTip->next->prev = _earTip->prev;
488 
489  // update reflex vertices list by checking the neighboring vertices
490  updateReflexVertex(_earTip->prev);
491  updateReflexVertex(_earTip->next);
492 
493  --numRemaningVertices_;
494  ++numTris_;
495 }
496 
497 
498 
499 }
bool isReflexVertex(int _i) const
Check if a vertex is reflex.
VectorT< float, 3 > Vec3f
Definition: VectorT.hh:125
virtual ~Triangulator()
Destructor.
Namespace providing different geometric functions concerning angles.
Definition: DBSCANT.cc:51
auto normalize() -> decltype(*this/=std::declval< VectorT< S, DIM >>().norm())
Definition: Vector11T.hh:428
decltype(std::declval< S >()*std::declval< S >()) sqrnorm() const
compute squared euclidean norm
Definition: Vector11T.hh:396
Triangulator(const std::vector< Vec3f > &_pos)
Execute the triangulation algorithm on a polygon.
Definition: Triangulator.cc:61