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
171 Triangulator::~Triangulator()
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
188 bool Triangulator::isReflexVertex(int i) const
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 }
ACG
Namespace providing different geometric functions concerning angles.
Definition: BaseObjectData.hh:68
OpenMesh::VectorT::sqrnorm
decltype(std::declval< S >() *std::declval< S >()) sqrnorm() const
compute squared euclidean norm
Definition: Vector11T.hh:397
ACG::Vec3f
VectorT< float, 3 > Vec3f
Definition: VectorT.hh:119
OpenMesh::VectorT::normalize
auto normalize() -> decltype(*this/=std::declval< VectorT< S, DIM >>().norm())
Definition: Vector11T.hh:429
ACG::Triangulator::~Triangulator
virtual ~Triangulator()
Destructor.
Definition: Triangulator.cc:171
OpenMesh::VectorT
Definition: Vector11T.hh:83
ACG::Triangulator::Triangulator
Triangulator(const std::vector< Vec3f > &_pos)
Execute the triangulation algorithm on a polygon.
Definition: Triangulator.cc:52
ACG::Triangulator::RingVertex
Definition: Triangulator.hh:143
ACG::Triangulator::isReflexVertex
bool isReflexVertex(int _i) const
Check if a vertex is reflex.
Definition: Triangulator.cc:188