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 #include "Triangulator.hh"
43 
44 
45 #include <iostream>
46 
47 
48 
49 namespace ACG {
50 
51 
52 Triangulator::Triangulator(const std::vector<Vec3f>& _pos)
53  : polySize_(_pos.size()), numRemaningVertices_(_pos.size()), numTris_(0),
54  numReflexVertices_(0),
55  ok_(false), convex_(false)
56 {
57  if (polySize_ < 3)
58  return;
59 
60 
61  if (polySize_ == 3)
62  {
63  numTris_ = 1;
64  tris_.resize(3);
65  tris_[0] = 0;
66  tris_[1] = 1;
67  tris_[2] = 2;
68 
69  numRemaningVertices_ = 0;
70  convex_ = true;
71  ok_ = true;
72  }
73  else
74  {
75  // project vertices onto the 2d plane of the polygon.
76  // the projection plane is chosen orthogonal to the polygon surface normal
77 
78  // use Newell's Method to compute the surface normal
79  // http://www.opengl.org/wiki/Calculating_a_Surface_Normal
80 
81  Vec3f n(0.0f, 0.0f, 0.0f);
82 
83  for (size_t i = 0; i < polySize_; ++i)
84  {
85  const size_t next = (i + 1) % polySize_;
86 
87  Vec3f a = _pos[i] - _pos[next];
88  Vec3f b = _pos[i] + _pos[next];
89 
90  n[0] += a[1] * b[2];
91  n[1] += a[2] * b[0];
92  n[2] += a[0] * b[1];
93  }
94 
95  // project to 2d
96  pos_.resize(polySize_);
97 
98  Vec3f axis[3] = { Vec3f(1.0f, 0.0f, 0.0f), Vec3f(0.0f, 1.0f, 0.0f), n };
99 
100  // orthonormalize projection axes
101  axis[2].normalize();
102 
103  // make sure first axis is linearly independent from the normal
104  while (std::abs(axis[0] | axis[2]) > 0.95f || (axis[0].sqrnorm() < 0.001f))
105  {
106  for (int i = 0; i < 3; ++i)
107  axis[0][i] = float(rand()) / float(RAND_MAX) * 2.0f - 1.0f;
108 
109  axis[0].normalize();
110  }
111 
112  // make axis[0] orthogonal to normal
113  axis[0] = axis[0] - axis[2] * (axis[0] | axis[2]);
114  axis[0].normalize();
115  axis[1] = axis[2] % axis[0];
116 
117 
118  for (size_t i = 0; i < polySize_; ++i)
119  {
120  // project onto polygon plane
121  pos_[i][0] = axis[0] | _pos[i];
122  pos_[i][1] = axis[1] | _pos[i];
123  }
124 
125 
126  // create triangle fans if there is at most one concave vertex
127  int reflexVertexID = 0;
128 
129  for (size_t i = 0; i < polySize_; ++i)
130  {
131  // test vertex (i+1)
132  if (isReflexVertex(pos_[i], pos_[(i + 1) % polySize_], pos_[(i + 2) % polySize_]))
133  {
134  ++numReflexVertices_;
135  reflexVertexID = (i + 1) % polySize_;
136  }
137  }
138 
139  convex_ = !numReflexVertices_;
140 
141 
142  if (numReflexVertices_ <= 1)
143  {
144  // create triangle fans
145  numTris_ = polySize_ - 2;
146  tris_.resize(numTris_ * 3);
147  numRemaningVertices_ = 0;
148  ok_ = true;
149 
150  for (size_t i = 0; i < numTris_; ++i)
151  {
152  tris_[i * 3] = reflexVertexID;
153  tris_[i * 3 + 1] = (reflexVertexID + i + 1) % polySize_;
154  tris_[i * 3 + 2] = (reflexVertexID + i + 2) % polySize_;
155  }
156 
157  }
158  else
159  {
160  // use the ear clipping algorithm
161 
162  earClippingN2();
163 // earClippingN3();
164 
165 // triangulateExternal();
166  }
167  }
168 }
169 
170 
172 {
173 }
174 
175 
176 bool Triangulator::isReflexVertex(const Vec2f& v0, const Vec2f& v1, const Vec2f& v2) const
177 {
178  // compute the sign of the cross product of the edges sharing vertex v1
179  // <0 : inner angle greater than 180 deg (reflex)
180  // >0 : inner angle less than 180 deg (convex)
181 
182  Vec2f u = v2 - v1;
183  Vec2f v = v0 - v1;
184 
185  return u[0] * v[1] - u[1] * v[0] < 0.0f;
186 }
187 
189 {
190  int p = (i + polySize_ - 1) % polySize_;
191  int n = (i + 1) % polySize_;
192 
193  return isReflexVertex(pos_[p], pos_[i], pos_[n]);
194 }
195 
196 float Triangulator::triangleAreaSign(const Vec2f& v0, const Vec2f& v1, const Vec2f& v2) const
197 {
198  // cross product
199  return (v0[0] - v2[0]) * (v1[1] - v2[1]) - (v1[0] - v2[0]) * (v0[1] - v2[1]);
200 }
201 
202 float Triangulator::distancePointToSegmentSq(const Vec2f& v0, const Vec2f& v1, const Vec2f& pt) const
203 {
204  Vec2f segment = v1 - v0;
205  Vec2f vec = pt - v0;
206  float dp = vec | segment;
207 
208  if (dp < 0.0f)
209  return vec.sqrnorm();
210 
211  float segSq = segment.sqrnorm();
212  dp /= segSq;
213 
214  if (dp < 1.0f)
215  return vec.sqrnorm() - dp * dp * segSq;
216 
217  vec = pt - v1;
218  return vec.sqrnorm();
219 }
220 
221 bool Triangulator::pointInTriangle(const Vec2f& v0, const Vec2f& v1, const Vec2f& v2, const Vec2f& pt) const
222 {
223  // ACG implementation based on barycentric coordinates (slow)
224 // return Geometry::isInTriangle(pt, v0, v1, v2);
225 
226 
227  // fast implementation based on triangle areas
228  // http://www.gamedev.net/topic/295943-is-this-a-better-point-in-triangle-test-2d/
229 
230  return triangleAreaSign(pt, v0, v1) >= 0.0f &&
231  triangleAreaSign(pt, v1, v2) >= 0.0f &&
232  triangleAreaSign(pt, v2, v0) >= 0.0f;
233 
234 
235 // // more accurate algorithm:
236 // // http://totologic.blogspot.de/2014/01/accurate-point-in-triangle-test.html
237 // // note: didn't improve accuracy at all for problematic polygons
238 //
239 // static const float eps = 1e-4f;
240 // static const float eps2 = eps*eps;
241 //
242 //
243 // // point in aabb of triangle
244 // Vec2f aabbMin = v0;
245 //
246 // aabbMin.minimize(v1);
247 // aabbMin.minimize(v2);
248 //
249 // if (pt[0] + eps < aabbMin[0] || pt[1] + eps < aabbMin[1])
250 // return false;
251 //
252 // Vec2f aabbMax = v0;
253 // aabbMax.maximize(v1);
254 // aabbMax.maximize(v2);
255 //
256 // if (pt[0] > aabbMax[0] + eps || pt[1] > aabbMax[1] + eps)
257 // return false;
258 //
259 //
260 // if (triangleAreaSign(pt, v0, v1) >= 0.0f &&
261 // triangleAreaSign(pt, v1, v2) >= 0.0f &&
262 // triangleAreaSign(pt, v2, v0) >= 0.0f)
263 // return true;
264 //
265 //
266 // return (distancePointToSegmentSq(v0, v1, pt) <= eps2 ||
267 // distancePointToSegmentSq(v1, v2, pt) <= eps2 ||
268 // distancePointToSegmentSq(v2, v0, pt) <= eps2);
269 }
270 
271 
272 void Triangulator::initVertexList()
273 {
274  vertices_.resize(polySize_);
275  reflexVertices_.clear();
276  for (size_t i = 0; i < polySize_; ++i)
277  {
278  size_t p = (i + polySize_ - 1) % polySize_;
279  size_t n = (i + 1) % polySize_;
280 
281  vertices_[i] = RingVertex(i, isReflexVertex(pos_[p], pos_[i], pos_[n]), pos_[i], &vertices_[p], &vertices_[n]);
282 
283  if (vertices_[i].reflex)
284  reflexVertices_.push_back(&vertices_[i]);
285  }
286 }
287 
288 
289 
290 int Triangulator::earClippingN3()
291 {
292  // http://www.geometrictools.com/Documentation/TriangulationByEarClipping.pdf
293  // O(n^3)
294 
295  numTris_ = 0;
296  ok_ = true;
297 
298  initVertexList();
299 
300  tris_.resize((polySize_ - 2) * 3);
301 
302  int numIterations = 0;
303  RingVertex* firstVertex = &vertices_[0];
304 
305  while (numTris_ < polySize_ - 2)
306  {
307  // find an ear in the remaining polygon
308  bool hasEars = false;
309 
310  RingVertex* curVertex = firstVertex;
311  do
312  {
313  curVertex->reflex = isReflexVertex(curVertex->prev->pos, curVertex->pos, curVertex->next->pos);
314 
315  if (!curVertex->reflex)
316  {
317  // triangle containment test
318  bool isEar = true;
319 
320  // check all remaining vertices r for containment
321  for (RingVertex* r = curVertex->next->next; r != curVertex->prev; r = r->next)
322  {
323  if (pointInTriangle(curVertex->prev->pos, curVertex->pos, curVertex->next->pos, r->pos))
324  {
325  isEar = false;
326  break;
327  }
328  }
329 
330  // found an ear
331  if (isEar)
332  {
333  // triangulate ear
334  hasEars = true;
335  tris_[numTris_ * 3] = curVertex->prev->id;
336  tris_[numTris_ * 3 + 1] = curVertex->id;
337  tris_[numTris_ * 3 + 2] = curVertex->next->id;
338  ++numTris_;
339 
340  // remove vertex from linked list
341  curVertex->prev->next = curVertex->next;
342  curVertex->next->prev = curVertex->prev;
343 
344  break;
345  }
346  }
347 
348  curVertex = curVertex->next;
349  ++numIterations;
350  } while (curVertex != firstVertex);
351 
352  firstVertex = firstVertex->next;
353 
354  // create triangle fans and hope for good result if there are no more ears
355  if (!hasEars && (numTris_ + 2 < polySize_))
356  {
357  for (RingVertex* iteratorVertex = firstVertex->next; iteratorVertex != firstVertex->prev; iteratorVertex = iteratorVertex->next)
358  {
359  tris_[numTris_ * 3] = firstVertex->id;
360  tris_[numTris_ * 3 + 1] = iteratorVertex->id;
361  tris_[numTris_ * 3 + 2] = iteratorVertex->next->id;
362  ++numTris_;
363  }
364 
365  assert(numTris_ == polySize_ - 2);
366  }
367  }
368 
369  return numTris_;
370 }
371 
372 
373 
374 int Triangulator::earClippingN2()
375 {
376  // http://www.geometrictools.com/Documentation/TriangulationByEarClipping.pdf
377  // O(n^2)
378 
379  numTris_ = 0;
380  ok_ = true;
381 
382  initVertexList();
383 
384  // triangulate
385 
386  size_t numTries = 0; // # checked vertices per iteration that aren't ears
387  numRemaningVertices_ = polySize_; // size of currently remaining polygon
388  int numIterations = 0;
389 
390  tris_.resize((polySize_ - 2) * 3);
391 
392  RingVertex* curVertex = &vertices_[0];
393  while (numRemaningVertices_ > 3)
394  {
395  // check if the current vertex is an ear tip
396  bool isEar = false;
397 
398  if (!curVertex->reflex)
399  {
400  // test current vertex for ear property
401  isEar = true;
402  for (std::list<RingVertex*>::iterator it = reflexVertices_.begin(); isEar && it != reflexVertices_.end(); ++it)
403  {
404  // skip direct neighbors
405  if (*it == curVertex->prev || *it == curVertex->next)
406  continue;
407 
408  // if any remaining vertex is inside the triangle, the current vertex is not an ear
409  if (pointInTriangle(curVertex->prev->pos, curVertex->pos, curVertex->next->pos, (*it)->pos))
410  isEar = false;
411  }
412 
413 
414  // found an ear
415  if (isEar)
416  {
417  addEar(curVertex);
418  numTries = 0;
419  }
420  }
421 
422  if (!isEar)
423  ++numTries;
424 
425  if (numTries > numRemaningVertices_)
426  {
427  // something went wrong
428  // create a triangle anyway and hope the result is ok
429 
430  addEar(curVertex);
431  numTries = 0;
432 
433  ok_ = false;
434  }
435 
436  curVertex = curVertex->next;
437  ++numIterations;
438  }
439 
440 
441  // add the last remaining triangle
442  if (numRemaningVertices_ == 3)
443  {
444  tris_[numTris_ * 3 + 0] = curVertex->prev->id;
445  tris_[numTris_ * 3 + 1] = curVertex->id;
446  tris_[numTris_ * 3 + 2] = curVertex->next->id;
447  ++numTris_;
448  }
449 
450 
451  return numTris_;
452 }
453 
454 bool Triangulator::updateReflexVertex(RingVertex* v)
455 {
456  if (v->reflex)
457  {
458  // check reflex property
459  v->reflex = isReflexVertex(v->prev->pos, v->pos, v->next->pos);
460 
461  // update list of reflex vertices
462  if (!v->reflex)
463  reflexVertices_.remove(v);
464  }
465 
466  return v->reflex;
467 }
468 
469 void Triangulator::addEar(RingVertex* _earTip)
470 {
471  // add ear triangle
472  tris_[numTris_ * 3 + 0] = _earTip->prev->id;
473  tris_[numTris_ * 3 + 1] = _earTip->id;
474  tris_[numTris_ * 3 + 2] = _earTip->next->id;
475 
476  // remove ear vertex from the linked list
477  _earTip->prev->next = _earTip->next;
478  _earTip->next->prev = _earTip->prev;
479 
480  // update reflex vertices list by checking the neighboring vertices
481  updateReflexVertex(_earTip->prev);
482  updateReflexVertex(_earTip->next);
483 
484  --numRemaningVertices_;
485  ++numTris_;
486 }
487 
488 
489 
490 }
virtual ~Triangulator()
Destructor.
Namespace providing different geometric functions concerning angles.
VectorT< float, 3 > Vec3f
Definition: VectorT.hh:119
Triangulator(const std::vector< Vec3f > &_pos)
Execute the triangulation algorithm on a polygon.
Definition: Triangulator.cc:52
decltype(std::declval< S >() *std::declval< S >()) sqrnorm() const
compute squared euclidean norm
Definition: Vector11T.hh:421
bool isReflexVertex(int _i) const
Check if a vertex is reflex.
auto normalize() -> decltype(*this/=std::declval< VectorT< S, DIM >>().norm())
Definition: Vector11T.hh:453