1
/***************************************************************************
2
* Copyright (c) 2005 Werner Mayer <werner.wm.mayer@gmx.de> *
4
* This file is part of the FreeCAD CAx development system. *
6
* This library is free software; you can redistribute it and/or *
7
* modify it under the terms of the GNU Library General Public *
8
* License as published by the Free Software Foundation; either *
9
* version 2 of the License, or (at your option) any later version. *
11
* This library is distributed in the hope that it will be useful, *
12
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
13
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
14
* GNU Library General Public License for more details. *
16
* You should have received a copy of the GNU Library General Public *
17
* License along with this library; see the file COPYING.LIB. If not, *
18
* write to the Free Software Foundation, Inc., 59 Temple Place, *
19
* Suite 330, Boston, MA 02111-1307, USA *
21
***************************************************************************/
24
#include "PreCompiled.h"
28
#include "Triangulation.h"
29
#include "Approximation.h"
30
#include "MeshKernel.h"
32
#include <Mod/Mesh/App/WildMagic4/Wm4Delaunay2.h>
35
using namespace MeshCore;
38
AbstractPolygonTriangulator::AbstractPolygonTriangulator()
43
AbstractPolygonTriangulator::~AbstractPolygonTriangulator()
47
void AbstractPolygonTriangulator::SetPolygon(const std::vector<Base::Vector3f>& raclPoints)
49
this->_points = raclPoints;
50
if (this->_points.size() > 0) {
51
if (this->_points.front() == this->_points.back())
52
this->_points.pop_back();
56
std::vector<Base::Vector3f> AbstractPolygonTriangulator::GetPolygon() const
61
float AbstractPolygonTriangulator::GetLength() const
64
if (_points.size() > 2) {
65
for (std::vector<Base::Vector3f>::const_iterator it = _points.begin(); it != _points.end();++it) {
66
std::vector<Base::Vector3f>::const_iterator jt = it + 1;
67
if (jt == _points.end()) jt = _points.begin();
68
len += Base::Distance(*it, *jt);
75
std::vector<Base::Vector3f> AbstractPolygonTriangulator::AddedPoints() const
77
// Apply the inverse transformation to project back to world coordinates
78
std::vector<Base::Vector3f> added;
79
added.reserve(_newpoints.size());
80
for (std::vector<Base::Vector3f>::const_iterator pt = _newpoints.begin(); pt != _newpoints.end(); ++pt)
81
added.push_back(_inverse * *pt);
85
Base::Matrix4D AbstractPolygonTriangulator::GetTransformToFitPlane() const
88
for (std::vector<Base::Vector3f>::const_iterator it = _points.begin(); it!=_points.end(); ++it)
89
planeFit.AddPoint(*it);
93
Base::Vector3f bs = planeFit.GetBase();
94
Base::Vector3f ex = planeFit.GetDirU();
95
Base::Vector3f ey = planeFit.GetDirV();
96
Base::Vector3f ez = planeFit.GetNormal();
98
// build the matrix for the inverse transformation
99
Base::Matrix4D rInverse;
100
rInverse.setToUnity();
101
rInverse[0][0] = ex.x; rInverse[0][1] = ey.x; rInverse[0][2] = ez.x; rInverse[0][3] = bs.x;
102
rInverse[1][0] = ex.y; rInverse[1][1] = ey.y; rInverse[1][2] = ez.y; rInverse[1][3] = bs.y;
103
rInverse[2][0] = ex.z; rInverse[2][1] = ey.z; rInverse[2][2] = ez.z; rInverse[2][3] = bs.z;
108
std::vector<Base::Vector3f> AbstractPolygonTriangulator::ProjectToFitPlane()
110
std::vector<Base::Vector3f> proj = _points;
111
_inverse = GetTransformToFitPlane();
112
Base::Vector3f bs((float)_inverse[0][3], (float)_inverse[1][3], (float)_inverse[2][3]);
113
Base::Vector3f ex((float)_inverse[0][0], (float)_inverse[1][0], (float)_inverse[2][0]);
114
Base::Vector3f ey((float)_inverse[0][1], (float)_inverse[1][1], (float)_inverse[2][1]);
115
for (std::vector<Base::Vector3f>::iterator jt = proj.begin(); jt!=proj.end(); ++jt)
116
jt->TransformToCoordinateSystem(bs, ex, ey);
120
void AbstractPolygonTriangulator::ProjectOntoSurface(const std::vector<Base::Vector3f>& points)
122
// For a good approximation we should have enough points, i.e. for 9 parameters
123
// for the fit function we should have at least 50 points.
124
unsigned int uMinPts = 50;
126
PolynomialFit polyFit;
127
Base::Vector3f bs((float)_inverse[0][3], (float)_inverse[1][3], (float)_inverse[2][3]);
128
Base::Vector3f ex((float)_inverse[0][0], (float)_inverse[1][0], (float)_inverse[2][0]);
129
Base::Vector3f ey((float)_inverse[0][1], (float)_inverse[1][1], (float)_inverse[2][1]);
131
for (std::vector<Base::Vector3f>::const_iterator it = points.begin(); it != points.end(); ++it) {
132
Base::Vector3f pt = *it;
133
pt.TransformToCoordinateSystem(bs, ex, ey);
134
polyFit.AddPoint(pt);
139
if (polyFit.CountPoints() >= uMinPts) {
140
for (std::vector<Base::Vector3f>::iterator pt = _newpoints.begin(); pt != _newpoints.end(); ++pt)
141
pt->z = (float)polyFit.Value(pt->x, pt->y);
145
bool AbstractPolygonTriangulator::TriangulatePolygon()
147
bool ok = Triangulate();
152
std::vector<unsigned long> AbstractPolygonTriangulator::GetInfo() const
157
void AbstractPolygonTriangulator::Discard()
165
void AbstractPolygonTriangulator::Done()
167
_info.push_back(_points.size());
171
// -------------------------------------------------------------
173
EarClippingTriangulator::EarClippingTriangulator()
177
EarClippingTriangulator::~EarClippingTriangulator()
181
bool EarClippingTriangulator::Triangulate()
186
std::vector<Base::Vector3f> pts = ProjectToFitPlane();
187
std::vector<unsigned long> result;
189
// Invoke the triangulator to triangulate this polygon.
190
Triangulate::Process(pts,result);
192
// print out the results.
193
unsigned long tcount = result.size()/3;
195
bool ok = tcount+2 == _points.size();
196
if (tcount > _points.size())
197
return false; // no valid triangulation
199
MeshGeomFacet clFacet;
200
MeshFacet clTopFacet;
201
for (unsigned long i=0; i<tcount; i++) {
202
if (Triangulate::_invert) {
203
clFacet._aclPoints[0] = _points[result[i*3+0]];
204
clFacet._aclPoints[2] = _points[result[i*3+1]];
205
clFacet._aclPoints[1] = _points[result[i*3+2]];
206
clTopFacet._aulPoints[0] = result[i*3+0];
207
clTopFacet._aulPoints[2] = result[i*3+1];
208
clTopFacet._aulPoints[1] = result[i*3+2];
211
clFacet._aclPoints[0] = _points[result[i*3+0]];
212
clFacet._aclPoints[1] = _points[result[i*3+1]];
213
clFacet._aclPoints[2] = _points[result[i*3+2]];
214
clTopFacet._aulPoints[0] = result[i*3+0];
215
clTopFacet._aulPoints[1] = result[i*3+1];
216
clTopFacet._aulPoints[2] = result[i*3+2];
219
_triangles.push_back(clFacet);
220
_facets.push_back(clTopFacet);
226
float EarClippingTriangulator::Triangulate::Area(const std::vector<Base::Vector3f> &contour)
228
int n = contour.size();
232
for(int p=n-1,q=0; q<n; p=q++) {
233
A+= contour[p].x*contour[q].y - contour[q].x*contour[p].y;
239
InsideTriangle decides if a point P is Inside of the triangle
242
bool EarClippingTriangulator::Triangulate::InsideTriangle(float Ax, float Ay, float Bx,
243
float By, float Cx, float Cy,
246
float ax, ay, bx, by, cx, cy, apx, apy, bpx, bpy, cpx, cpy;
247
float cCROSSap, bCROSScp, aCROSSbp;
249
ax = Cx - Bx; ay = Cy - By;
250
bx = Ax - Cx; by = Ay - Cy;
251
cx = Bx - Ax; cy = By - Ay;
252
apx= Px - Ax; apy= Py - Ay;
253
bpx= Px - Bx; bpy= Py - By;
254
cpx= Px - Cx; cpy= Py - Cy;
256
aCROSSbp = ax*bpy - ay*bpx;
257
cCROSSap = cx*apy - cy*apx;
258
bCROSScp = bx*cpy - by*cpx;
260
return ((aCROSSbp >= FLOAT_EPS) && (bCROSScp >= FLOAT_EPS) && (cCROSSap >= FLOAT_EPS));
263
bool EarClippingTriangulator::Triangulate::Snip(const std::vector<Base::Vector3f> &contour,
264
int u,int v,int w,int n,int *V)
267
float Ax, Ay, Bx, By, Cx, Cy, Px, Py;
269
Ax = contour[V[u]].x;
270
Ay = contour[V[u]].y;
272
Bx = contour[V[v]].x;
273
By = contour[V[v]].y;
275
Cx = contour[V[w]].x;
276
Cy = contour[V[w]].y;
278
if (FLOAT_EPS > (((Bx-Ax)*(Cy-Ay)) - ((By-Ay)*(Cx-Ax)))) return false;
281
if( (p == u) || (p == v) || (p == w) ) continue;
282
Px = contour[V[p]].x;
283
Py = contour[V[p]].y;
284
if (InsideTriangle(Ax,Ay,Bx,By,Cx,Cy,Px,Py)) return false;
290
bool EarClippingTriangulator::Triangulate::_invert = false;
292
bool EarClippingTriangulator::Triangulate::Process(const std::vector<Base::Vector3f> &contour,
293
std::vector<unsigned long> &result)
295
/* allocate and initialize list of Vertices in polygon */
297
int n = contour.size();
298
if ( n < 3 ) return false;
302
/* we want a counter-clockwise polygon in V */
304
if (0.0f < Area(contour)) {
305
for (int v=0; v<n; v++) V[v] = v;
308
// for(int v=0; v<n; v++) V[v] = (n-1)-v;
310
for(int v=0; v<n; v++) V[v] = (n-1)-v;
316
/* remove nv-2 Vertices, creating 1 triangle every time */
317
int count = 2*nv; /* error detection */
319
for(int m=0, v=nv-1; nv>2; ) {
320
/* if we loop, it is probably a non-simple polygon */
321
if (0 >= (count--)) {
322
//** Triangulate: ERROR - probable bad polygon!
326
/* three consecutive vertices in current polygon, <u,v,w> */
327
int u = v ; if (nv <= u) u = 0; /* previous */
328
v = u+1; if (nv <= v) v = 0; /* new v */
329
int w = v+1; if (nv <= w) w = 0; /* next */
331
if (Snip(contour,u,v,w,nv,V)) {
334
/* true names of the vertices */
335
a = V[u]; b = V[v]; c = V[w];
337
/* output Triangle */
338
result.push_back( a );
339
result.push_back( b );
340
result.push_back( c );
344
/* remove v from remaining polygon */
345
for(s=v,t=v+1;t<nv;s++,t++) V[s] = V[t]; nv--;
347
/* resest error detection counter */
357
// -------------------------------------------------------------
359
QuasiDelaunayTriangulator::QuasiDelaunayTriangulator()
363
QuasiDelaunayTriangulator::~QuasiDelaunayTriangulator()
367
bool QuasiDelaunayTriangulator::Triangulate()
369
if (EarClippingTriangulator::Triangulate() == false)
370
return false; // no valid triangulation
372
// For each internal edge get the adjacent facets. When doing an edge swap we must update
374
std::map<std::pair<unsigned long, unsigned long>, std::vector<unsigned long> > aEdge2Face;
375
for (std::vector<MeshFacet>::iterator pI = _facets.begin(); pI != _facets.end(); pI++) {
376
for (int i = 0; i < 3; i++) {
377
unsigned long ulPt0 = std::min<unsigned long>(pI->_aulPoints[i], pI->_aulPoints[(i+1)%3]);
378
unsigned long ulPt1 = std::max<unsigned long>(pI->_aulPoints[i], pI->_aulPoints[(i+1)%3]);
379
// ignore borderlines of the polygon
380
if ((ulPt1-ulPt0)%(_points.size()-1) > 1)
381
aEdge2Face[std::pair<unsigned long, unsigned long>(ulPt0, ulPt1)].push_back(pI - _facets.begin());
385
// fill up this list with all internal edges and perform swap edges until this list is empty
386
std::list<std::pair<unsigned long, unsigned long> > aEdgeList;
387
std::map<std::pair<unsigned long, unsigned long>, std::vector<unsigned long> >::iterator pE;
388
for (pE = aEdge2Face.begin(); pE != aEdge2Face.end(); ++pE)
389
aEdgeList.push_back(pE->first);
391
// to be sure to avoid an endless loop
392
unsigned long uMaxIter = 5 * aEdge2Face.size();
394
// Perform a swap edge where needed
395
while (!aEdgeList.empty() && uMaxIter > 0) {
396
// get the first edge and remove it from the list
397
std::pair<unsigned long, unsigned long> aEdge = aEdgeList.front();
398
aEdgeList.pop_front();
401
// get the adjacent facets to this edge
402
pE = aEdge2Face.find( aEdge );
404
// this edge has been removed some iterations before
405
if (pE == aEdge2Face.end())
408
MeshFacet& rF1 = _facets[pE->second[0]];
409
MeshFacet& rF2 = _facets[pE->second[1]];
410
unsigned short side1 = rF1.Side(aEdge.first, aEdge.second);
412
Base::Vector3f cP1 = _points[rF1._aulPoints[side1]];
413
Base::Vector3f cP2 = _points[rF1._aulPoints[(side1+1)%3]];
414
Base::Vector3f cP3 = _points[rF1._aulPoints[(side1+2)%3]];
416
unsigned short side2 = rF2.Side(aEdge.first, aEdge.second);
417
Base::Vector3f cP4 = _points[rF2._aulPoints[(side2+2)%3]];
419
MeshGeomFacet cT1(cP1, cP2, cP3); float fMax1 = cT1.MaximumAngle();
420
MeshGeomFacet cT2(cP2, cP1, cP4); float fMax2 = cT2.MaximumAngle();
421
MeshGeomFacet cT3(cP4, cP3, cP1); float fMax3 = cT3.MaximumAngle();
422
MeshGeomFacet cT4(cP3, cP4, cP2); float fMax4 = cT4.MaximumAngle();
424
float fMax12 = std::max<float>(fMax1, fMax2);
425
float fMax34 = std::max<float>(fMax3, fMax4);
427
// We must make sure that the two adjacent triangles builds a convex polygon, otherwise
428
// the swap edge operation is illegal
429
Base::Vector3f cU = cP2-cP1;
430
Base::Vector3f cV = cP4-cP3;
431
// build a helper plane through cP1 that must separate cP3 and cP4
432
Base::Vector3f cN1 = (cU % cV) % cU;
433
if (((cP3-cP1)*cN1)*((cP4-cP1)*cN1) >= 0.0f)
434
continue; // not convex
435
// build a helper plane through cP3 that must separate cP1 and cP2
436
Base::Vector3f cN2 = (cU % cV) % cV;
437
if (((cP1-cP3)*cN2)*((cP2-cP3)*cN2) >= 0.0f)
438
continue; // not convex
440
// ok, here we should perform a swap edge to minimize the maximum angle
441
if (fMax12 > fMax34) {
442
rF1._aulPoints[(side1+1)%3] = rF2._aulPoints[(side2+2)%3];
443
rF2._aulPoints[(side2+1)%3] = rF1._aulPoints[(side1+2)%3];
445
// adjust the edge list
446
for (int i=0; i<3; i++) {
447
std::map<std::pair<unsigned long, unsigned long>, std::vector<unsigned long> >::iterator it;
449
unsigned long ulPt0 = std::min<unsigned long>(rF1._aulPoints[i], rF1._aulPoints[(i+1)%3]);
450
unsigned long ulPt1 = std::max<unsigned long>(rF1._aulPoints[i], rF1._aulPoints[(i+1)%3]);
451
it = aEdge2Face.find( std::make_pair(ulPt0, ulPt1) );
452
if (it != aEdge2Face.end()) {
453
if (it->second[0] == pE->second[1])
454
it->second[0] = pE->second[0];
455
else if (it->second[1] == pE->second[1])
456
it->second[1] = pE->second[0];
457
aEdgeList.push_back(it->first);
461
ulPt0 = std::min<unsigned long>(rF2._aulPoints[i], rF2._aulPoints[(i+1)%3]);
462
ulPt1 = std::max<unsigned long>(rF2._aulPoints[i], rF2._aulPoints[(i+1)%3]);
463
it = aEdge2Face.find( std::make_pair(ulPt0, ulPt1) );
464
if (it != aEdge2Face.end()) {
465
if (it->second[0] == pE->second[0])
466
it->second[0] = pE->second[1];
467
else if (it->second[1] == pE->second[0])
468
it->second[1] = pE->second[1];
469
aEdgeList.push_back(it->first);
473
// Now we must remove the edge and replace it through the new edge
474
unsigned long ulPt0 = std::min<unsigned long>(rF1._aulPoints[(side1+1)%3], rF2._aulPoints[(side2+1)%3]);
475
unsigned long ulPt1 = std::max<unsigned long>(rF1._aulPoints[(side1+1)%3], rF2._aulPoints[(side2+1)%3]);
476
std::pair<unsigned long, unsigned long> aNewEdge = std::make_pair(ulPt0, ulPt1);
477
aEdge2Face[aNewEdge] = pE->second;
478
aEdge2Face.erase(pE);
485
// -------------------------------------------------------------
488
namespace Triangulation {
489
struct Vertex2d_Less : public std::binary_function<const Base::Vector3f&, const Base::Vector3f&, bool>
491
bool operator()(const Base::Vector3f& p, const Base::Vector3f& q) const
493
if (fabs(p.x - q.x) < MeshDefinitions::_fMinPointDistanceD1) {
494
if (fabs(p.y - q.y) < MeshDefinitions::_fMinPointDistanceD1) {
495
return false; } else return p.y < q.y;
496
} else return p.x < q.x; return true;
499
struct Vertex2d_EqualTo : public std::binary_function<const Base::Vector3f&, const Base::Vector3f&, bool>
501
bool operator()(const Base::Vector3f& p, const Base::Vector3f& q) const
503
if (fabs(p.x - q.x) < MeshDefinitions::_fMinPointDistanceD1
504
&& fabs(p.y - q.y) < MeshDefinitions::_fMinPointDistanceD1)
505
return true; return false;
511
DelaunayTriangulator::DelaunayTriangulator()
515
DelaunayTriangulator::~DelaunayTriangulator()
519
bool DelaunayTriangulator::Triangulate()
521
// before starting the triangulation we must make sure that all polygon
522
// points are different
523
std::vector<Base::Vector3f> aPoints = _points;
524
// sort the points ascending x,y coordinates
525
std::sort(aPoints.begin(), aPoints.end(), Triangulation::Vertex2d_Less());
526
// if there are two adjacent points whose distance is less then an epsilon
527
if (std::adjacent_find(aPoints.begin(), aPoints.end(),
528
Triangulation::Vertex2d_EqualTo()) < aPoints.end())
534
std::vector<Wm4::Vector2d> akVertex;
535
akVertex.reserve(_points.size());
536
for (std::vector<Base::Vector3f>::iterator it = _points.begin(); it != _points.end(); it++) {
537
akVertex.push_back(Wm4::Vector2d(it->x, it->y));
540
Wm4::Delaunay2d del(akVertex.size(), &(akVertex[0]), 0.001, false, Wm4::Query::QT_INT64);
541
int iTQuantity = del.GetSimplexQuantity();
542
std::vector<int> aiTVertex(3*iTQuantity);
543
size_t uiSize = 3*iTQuantity*sizeof(int);
544
Wm4::System::Memcpy(&(aiTVertex[0]),uiSize,del.GetIndices(),uiSize);
546
// If H is the number of hull edges and N is the number of vertices,
547
// then the triangulation must have 2*N-2-H triangles and 3*N-3-H
551
del.GetHull(iEQuantity,aiIndex);
552
int iUniqueVQuantity = del.GetUniqueVertexQuantity();
553
int iTVerify = 2*iUniqueVQuantity - 2 - iEQuantity;
554
(void)iTVerify; // avoid warning in release build
555
bool succeeded = (iTVerify == iTQuantity);
556
int iEVerify = 3*iUniqueVQuantity - 3 - iEQuantity;
557
(void)iEVerify; // avoid warning about unused variable
560
MeshGeomFacet triangle;
562
for (int i = 0; i < iTQuantity; i++) {
563
for (int j=0; j<3; j++) {
564
facet._aulPoints[j] = aiTVertex[3*i+j];
565
triangle._aclPoints[j].x = (float)akVertex[aiTVertex[3*i+j]].X();
566
triangle._aclPoints[j].y = (float)akVertex[aiTVertex[3*i+j]].Y();
569
_triangles.push_back(triangle);
570
_facets.push_back(facet);
576
// -------------------------------------------------------------
578
FlatTriangulator::FlatTriangulator()
582
FlatTriangulator::~FlatTriangulator()
586
bool FlatTriangulator::Triangulate()
589
// before starting the triangulation we must make sure that all polygon
590
// points are different
591
std::vector<Base::Vector3f> aPoints = ProjectToFitPlane();
592
std::vector<Base::Vector3f> tmp = aPoints;
593
// sort the points ascending x,y coordinates
594
std::sort(tmp.begin(), tmp.end(), Triangulation::Vertex2d_Less());
595
// if there are two adjacent points whose distance is less then an epsilon
596
if (std::adjacent_find(tmp.begin(), tmp.end(),
597
Triangulation::Vertex2d_EqualTo()) < tmp.end() )
603
// Todo: Implement algorithm for constraint delaunay triangulation
604
QuasiDelaunayTriangulator tria;
605
tria.SetPolygon(this->GetPolygon());
606
bool succeeded = tria.TriangulatePolygon();
607
this->_facets = tria.GetFacets();
608
this->_triangles = tria.GetTriangles();
613
void FlatTriangulator::ProjectOntoSurface(const std::vector<Base::Vector3f>&)
617
// -------------------------------------------------------------
619
ConstraintDelaunayTriangulator::ConstraintDelaunayTriangulator(float area)
624
ConstraintDelaunayTriangulator::~ConstraintDelaunayTriangulator()
628
bool ConstraintDelaunayTriangulator::Triangulate()
631
// before starting the triangulation we must make sure that all polygon
632
// points are different
633
std::vector<Base::Vector3f> aPoints = ProjectToFitPlane();
634
std::vector<Base::Vector3f> tmp = aPoints;
635
// sort the points ascending x,y coordinates
636
std::sort(tmp.begin(), tmp.end(), Triangulation::Vertex2d_Less());
637
// if there are two adjacent points whose distance is less then an epsilon
638
if (std::adjacent_find(tmp.begin(), tmp.end(),
639
Triangulation::Vertex2d_EqualTo()) < tmp.end() )
645
// Todo: Implement algorithm for constraint delaunay triangulation
646
QuasiDelaunayTriangulator tria;
647
tria.SetPolygon(this->GetPolygon());
648
bool succeeded = tria.TriangulatePolygon();
649
this->_facets = tria.GetFacets();
650
this->_triangles = tria.GetTriangles();
added
removed
src/Mod/Mesh/App/Core/Triangulation.cpp
