Developer Documentation

- Main Page
- Related Pages
- Modules
- Namespaces
- Classes
- Files

- File List
- File Members

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Modules Pages
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
180 Triangulator::~Triangulator()
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
197 bool Triangulator::isReflexVertex(int i) const
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 }
ACG
Namespace providing different geometric functions concerning angles.
Definition: DBSCANT.cc:51
ACG::Triangulator::~Triangulator
virtual ~Triangulator()
Destructor.
Definition: Triangulator.cc:180
ACG::Vec3f
VectorT< float, 3 > Vec3f
Definition: VectorT.hh:125
OpenMesh::VectorT
Definition: Vector11T.hh:83
ACG::Triangulator::isReflexVertex
bool isReflexVertex(int _i) const
Check if a vertex is reflex.
Definition: Triangulator.cc:197
ACG::Triangulator::Triangulator
Triangulator(const std::vector< Vec3f > &_pos)
Execute the triangulation algorithm on a polygon.
Definition: Triangulator.cc:61