FreeCAD

1
/***************************************************************************
2
 *   Copyright (c) 2005 Werner Mayer <wmayer[at]users.sourceforge.net>     *
3
 *                                                                         *
4
 *   This file is part of the FreeCAD CAx development system.              *
5
 *                                                                         *
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.      *
10
 *                                                                         *
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.                  *
15
 *                                                                         *
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                                *
20
 *                                                                         *
21
 ***************************************************************************/
22

23
#include "PreCompiled.h"
24
#ifndef _PreComp_
25
#include <queue>
26
#endif
27

28
#include <Base/Console.h>
29
#include <Base/Exception.h>
30
#include <Mod/Mesh/App/WildMagic4/Wm4Delaunay2.h>
31

32
#include "Approximation.h"
33
#include "MeshKernel.h"
34
#include "Triangulation.h"
35

36

37
using namespace MeshCore;
38

39

40
bool TriangulationVerifier::Accept(const Base::Vector3f& n,
41
                                   const Base::Vector3f& p1,
42
                                   const Base::Vector3f& p2,
43
                                   const Base::Vector3f& p3) const
44
{
45
    float ref_dist = (p2 - p1) * n;
46
    float tri_dist = (p3 - p1) * n;
47
    return (ref_dist * tri_dist <= 0.0f);
48
}
49

50
bool TriangulationVerifier::MustFlip(const Base::Vector3f& n1, const Base::Vector3f& n2) const
51
{
52
    return n1.Dot(n2) <= 0.0f;
53
}
54

55
bool TriangulationVerifierV2::Accept(const Base::Vector3f& n,
56
                                     const Base::Vector3f& p1,
57
                                     const Base::Vector3f& p2,
58
                                     const Base::Vector3f& p3) const
59
{
60
    float ref_dist = (p2 - p1) * n;
61
    float tri_dist = (p3 - p1) * n;
62
    float prod = ref_dist * tri_dist;
63
    (void)prod;
64
    return true;
65
}
66

67
bool TriangulationVerifierV2::MustFlip(const Base::Vector3f& n1, const Base::Vector3f& n2) const
68
{
69
    float dot = n1.Dot(n2);
70
    (void)dot;
71
    return false;
72
}
73

74
// ----------------------------------------------------------------------------
75

76
AbstractPolygonTriangulator::AbstractPolygonTriangulator()
77
    : _discard {false}
78
    , _verifier {new TriangulationVerifier()}
79
{}
80

81
AbstractPolygonTriangulator::~AbstractPolygonTriangulator()
82
{
83
    delete _verifier;
84
}
85

86
TriangulationVerifier* AbstractPolygonTriangulator::GetVerifier() const
87
{
88
    return _verifier;
89
}
90

91
void AbstractPolygonTriangulator::SetVerifier(TriangulationVerifier* v)
92
{
93
    delete _verifier;
94
    _verifier = v;
95
}
96

97
void AbstractPolygonTriangulator::SetPolygon(const std::vector<Base::Vector3f>& raclPoints)
98
{
99
    this->_points = raclPoints;
100
    if (!this->_points.empty()) {
101
        if (this->_points.front() == this->_points.back()) {
102
            this->_points.pop_back();
103
        }
104
    }
105
}
106

107
std::vector<Base::Vector3f> AbstractPolygonTriangulator::GetPolygon() const
108
{
109
    return _points;
110
}
111

112
float AbstractPolygonTriangulator::GetLength() const
113
{
114
    float len = 0.0f;
115
    if (_points.size() > 2) {
116
        for (std::vector<Base::Vector3f>::const_iterator it = _points.begin(); it != _points.end();
117
             ++it) {
118
            std::vector<Base::Vector3f>::const_iterator jt = it + 1;
119
            if (jt == _points.end()) {
120
                jt = _points.begin();
121
            }
122
            len += Base::Distance(*it, *jt);
123
        }
124
    }
125

126
    return len;
127
}
128

129
std::vector<Base::Vector3f> AbstractPolygonTriangulator::AddedPoints() const
130
{
131
    // Apply the inverse transformation to project back to world coordinates
132
    std::vector<Base::Vector3f> added;
133
    added.reserve(_newpoints.size());
134
    for (auto point : _newpoints) {
135
        added.push_back(_inverse * point);
136
    }
137
    return added;
138
}
139

140
Base::Matrix4D AbstractPolygonTriangulator::GetTransformToFitPlane() const
141
{
142
    PlaneFit planeFit;
143
    for (auto point : _points) {
144
        planeFit.AddPoint(point);
145
    }
146

147
    if (planeFit.Fit() >= FLOAT_MAX) {
148
        throw Base::RuntimeError("Plane fit failed");
149
    }
150

151
    Base::Vector3f bs = planeFit.GetBase();
152
    Base::Vector3f ex = planeFit.GetDirU();
153
    Base::Vector3f ey = planeFit.GetDirV();
154
    Base::Vector3f ez = planeFit.GetNormal();
155

156
    // build the matrix for the inverse transformation
157
    Base::Matrix4D rInverse;
158
    rInverse.setToUnity();
159
    rInverse[0][0] = static_cast<double>(ex.x);
160
    rInverse[0][1] = static_cast<double>(ey.x);
161
    rInverse[0][2] = static_cast<double>(ez.x);
162
    rInverse[0][3] = static_cast<double>(bs.x);
163

164
    rInverse[1][0] = static_cast<double>(ex.y);
165
    rInverse[1][1] = static_cast<double>(ey.y);
166
    rInverse[1][2] = static_cast<double>(ez.y);
167
    rInverse[1][3] = static_cast<double>(bs.y);
168

169
    rInverse[2][0] = static_cast<double>(ex.z);
170
    rInverse[2][1] = static_cast<double>(ey.z);
171
    rInverse[2][2] = static_cast<double>(ez.z);
172
    rInverse[2][3] = static_cast<double>(bs.z);
173

174
    return rInverse;
175
}
176

177
std::vector<Base::Vector3f> AbstractPolygonTriangulator::ProjectToFitPlane()
178
{
179
    std::vector<Base::Vector3f> proj = _points;
180
    _inverse = GetTransformToFitPlane();
181
    Base::Vector3f bs(static_cast<float>(_inverse[0][3]),
182
                      static_cast<float>(_inverse[1][3]),
183
                      static_cast<float>(_inverse[2][3]));
184
    Base::Vector3f ex(static_cast<float>(_inverse[0][0]),
185
                      static_cast<float>(_inverse[1][0]),
186
                      static_cast<float>(_inverse[2][0]));
187
    Base::Vector3f ey(static_cast<float>(_inverse[0][1]),
188
                      static_cast<float>(_inverse[1][1]),
189
                      static_cast<float>(_inverse[2][1]));
190
    for (auto& jt : proj) {
191
        jt.TransformToCoordinateSystem(bs, ex, ey);
192
    }
193
    return proj;
194
}
195

196
void AbstractPolygonTriangulator::PostProcessing(const std::vector<Base::Vector3f>& points)
197
{
198
    // For a good approximation we should have enough points, i.e. for 9 parameters
199
    // for the fit function we should have at least 50 points.
200
    unsigned int uMinPts = 50;
201

202
    PolynomialFit polyFit;
203
    Base::Vector3f bs(static_cast<float>(_inverse[0][3]),
204
                      static_cast<float>(_inverse[1][3]),
205
                      static_cast<float>(_inverse[2][3]));
206
    Base::Vector3f ex(static_cast<float>(_inverse[0][0]),
207
                      static_cast<float>(_inverse[1][0]),
208
                      static_cast<float>(_inverse[2][0]));
209
    Base::Vector3f ey(static_cast<float>(_inverse[0][1]),
210
                      static_cast<float>(_inverse[1][1]),
211
                      static_cast<float>(_inverse[2][1]));
212

213
    for (auto pt : points) {
214
        pt.TransformToCoordinateSystem(bs, ex, ey);
215
        polyFit.AddPoint(pt);
216
    }
217

218
    if (polyFit.CountPoints() >= uMinPts && polyFit.Fit() < FLOAT_MAX) {
219
        for (auto& newpoint : _newpoints) {
220
            newpoint.z = static_cast<float>(polyFit.Value(newpoint.x, newpoint.y));
221
        }
222
    }
223
}
224

225
MeshGeomFacet AbstractPolygonTriangulator::GetTriangle(const MeshPointArray& points,
226
                                                       const MeshFacet& facet) const
227
{
228
    MeshGeomFacet triangle;
229
    triangle._aclPoints[0] = points[facet._aulPoints[0]];
230
    triangle._aclPoints[1] = points[facet._aulPoints[1]];
231
    triangle._aclPoints[2] = points[facet._aulPoints[2]];
232
    return triangle;
233
}
234

235
bool AbstractPolygonTriangulator::TriangulatePolygon()
236
{
237
    try {
238
        if (!this->_indices.empty() && this->_points.size() != this->_indices.size()) {
239
            Base::Console().Log("Triangulation: %d points <> %d indices\n",
240
                                _points.size(),
241
                                _indices.size());
242
            return false;
243
        }
244
        bool ok = Triangulate();
245
        if (ok) {
246
            Done();
247
        }
248
        return ok;
249
    }
250
    catch (const Base::Exception& e) {
251
        Base::Console().Log("Triangulation: %s\n", e.what());
252
        return false;
253
    }
254
    catch (const std::exception& e) {
255
        Base::Console().Log("Triangulation: %s\n", e.what());
256
        return false;
257
    }
258
    catch (...) {
259
        return false;
260
    }
261
}
262

263
std::vector<PointIndex> AbstractPolygonTriangulator::GetInfo() const
264
{
265
    return _info;
266
}
267

268
void AbstractPolygonTriangulator::Discard()
269
{
270
    if (!_discard) {
271
        _discard = true;
272
        _info.pop_back();
273
    }
274
}
275

276
void AbstractPolygonTriangulator::Reset()
277
{}
278

279
void AbstractPolygonTriangulator::Done()
280
{
281
    _info.push_back(_points.size());
282
    _discard = false;
283
}
284

285
// -------------------------------------------------------------
286

287
EarClippingTriangulator::EarClippingTriangulator() = default;
288

289
bool EarClippingTriangulator::Triangulate()
290
{
291
    _facets.clear();
292
    _triangles.clear();
293

294
    std::vector<Base::Vector3f> pts = ProjectToFitPlane();
295
    std::vector<PointIndex> result;
296

297
    //  Invoke the triangulator to triangulate this polygon.
298
    Triangulate::Process(pts, result);
299

300
    // print out the results.
301
    size_t tcount = result.size() / 3;
302

303
    bool ok = tcount + 2 == _points.size();
304
    if (tcount > _points.size()) {
305
        return false;  // no valid triangulation
306
    }
307

308
    MeshGeomFacet clFacet;
309
    MeshFacet clTopFacet;
310
    for (size_t i = 0; i < tcount; i++) {
311
        if (Triangulate::_invert) {
312
            clFacet._aclPoints[0] = _points[result[i * 3 + 0]];
313
            clFacet._aclPoints[2] = _points[result[i * 3 + 1]];
314
            clFacet._aclPoints[1] = _points[result[i * 3 + 2]];
315
            clTopFacet._aulPoints[0] = result[i * 3 + 0];
316
            clTopFacet._aulPoints[2] = result[i * 3 + 1];
317
            clTopFacet._aulPoints[1] = result[i * 3 + 2];
318
        }
319
        else {
320
            clFacet._aclPoints[0] = _points[result[i * 3 + 0]];
321
            clFacet._aclPoints[1] = _points[result[i * 3 + 1]];
322
            clFacet._aclPoints[2] = _points[result[i * 3 + 2]];
323
            clTopFacet._aulPoints[0] = result[i * 3 + 0];
324
            clTopFacet._aulPoints[1] = result[i * 3 + 1];
325
            clTopFacet._aulPoints[2] = result[i * 3 + 2];
326
        }
327

328
        _triangles.push_back(clFacet);
329
        _facets.push_back(clTopFacet);
330
    }
331

332
    return ok;
333
}
334

335
float EarClippingTriangulator::Triangulate::Area(const std::vector<Base::Vector3f>& contour)
336
{
337
    int n = contour.size();
338

339
    float A = 0.0f;
340

341
    for (int p = n - 1, q = 0; q < n; p = q++) {
342
        A += contour[p].x * contour[q].y - contour[q].x * contour[p].y;
343
    }
344
    return A * 0.5f;
345
}
346

347
/*
348
  InsideTriangle decides if a point P is Inside of the triangle
349
  defined by A, B, C.
350
*/
351
bool EarClippingTriangulator::Triangulate::InsideTriangle(float Ax,
352
                                                          float Ay,
353
                                                          float Bx,
354
                                                          float By,
355
                                                          float Cx,
356
                                                          float Cy,
357
                                                          float Px,
358
                                                          float Py)
359
{
360
    float ax {}, ay {}, bx {}, by {}, cx {}, cy {}, apx {}, apy {}, bpx {}, bpy {}, cpx {}, cpy {};
361
    float cCROSSap {}, bCROSScp {}, aCROSSbp {};
362

363
    ax = Cx - Bx;
364
    ay = Cy - By;
365
    bx = Ax - Cx;
366
    by = Ay - Cy;
367
    cx = Bx - Ax;
368
    cy = By - Ay;
369
    apx = Px - Ax;
370
    apy = Py - Ay;
371
    bpx = Px - Bx;
372
    bpy = Py - By;
373
    cpx = Px - Cx;
374
    cpy = Py - Cy;
375

376
    aCROSSbp = ax * bpy - ay * bpx;
377
    cCROSSap = cx * apy - cy * apx;
378
    bCROSScp = bx * cpy - by * cpx;
379

380
    return ((aCROSSbp >= FLOAT_EPS) && (bCROSScp >= FLOAT_EPS) && (cCROSSap >= FLOAT_EPS));
381
}
382

383
bool EarClippingTriangulator::Triangulate::Snip(const std::vector<Base::Vector3f>& contour,
384
                                                int u,
385
                                                int v,
386
                                                int w,
387
                                                int n,
388
                                                int* V)
389
{
390
    int p {};
391
    float Ax {}, Ay {}, Bx {}, By {}, Cx {}, Cy {}, Px {}, Py {};
392

393
    Ax = contour[V[u]].x;
394
    Ay = contour[V[u]].y;
395

396
    Bx = contour[V[v]].x;
397
    By = contour[V[v]].y;
398

399
    Cx = contour[V[w]].x;
400
    Cy = contour[V[w]].y;
401

402
    if (FLOAT_EPS > (((Bx - Ax) * (Cy - Ay)) - ((By - Ay) * (Cx - Ax)))) {
403
        return false;
404
    }
405

406
    for (p = 0; p < n; p++) {
407
        if ((p == u) || (p == v) || (p == w)) {
408
            continue;
409
        }
410
        Px = contour[V[p]].x;
411
        Py = contour[V[p]].y;
412
        if (InsideTriangle(Ax, Ay, Bx, By, Cx, Cy, Px, Py)) {
413
            return false;
414
        }
415
    }
416

417
    return true;
418
}
419

420
bool EarClippingTriangulator::Triangulate::_invert = false;
421

422
bool EarClippingTriangulator::Triangulate::Process(const std::vector<Base::Vector3f>& contour,
423
                                                   std::vector<PointIndex>& result)
424
{
425
    /* allocate and initialize list of Vertices in polygon */
426

427
    int n = contour.size();
428
    if (n < 3) {
429
        return false;
430
    }
431

432
    int* V = new int[n];
433

434
    /* we want a counter-clockwise polygon in V */
435

436
    if (0.0f < Area(contour)) {
437
        for (int v = 0; v < n; v++) {
438
            V[v] = v;
439
        }
440
        _invert = true;
441
    }
442
    //    for(int v=0; v<n; v++) V[v] = (n-1)-v;
443
    else {
444
        for (int v = 0; v < n; v++) {
445
            V[v] = (n - 1) - v;
446
        }
447
        _invert = false;
448
    }
449

450
    int nv = n;
451

452
    /*  remove nv-2 Vertices, creating 1 triangle every time */
453
    int count = 2 * nv; /* error detection */
454

455
    for (int v = nv - 1; nv > 2;) {
456
        /* if we loop, it is probably a non-simple polygon */
457
        if (0 >= (count--)) {
458
            //** Triangulate: ERROR - probable bad polygon!
459
            delete[] V;
460
            return false;
461
        }
462

463
        /* three consecutive vertices in current polygon, <u,v,w> */
464
        int u = v;
465
        if (nv <= u) {
466
            u = 0; /* previous */
467
        }
468
        v = u + 1;
469
        if (nv <= v) {
470
            v = 0; /* new v    */
471
        }
472
        int w = v + 1;
473
        if (nv <= w) {
474
            w = 0; /* next     */
475
        }
476

477
        if (Snip(contour, u, v, w, nv, V)) {
478
            int a {}, b {}, c {}, s {}, t {};
479

480
            /* true names of the vertices */
481
            a = V[u];
482
            b = V[v];
483
            c = V[w];
484

485
            /* output Triangle */
486
            result.push_back(a);
487
            result.push_back(b);
488
            result.push_back(c);
489

490
            /* remove v from remaining polygon */
491
            for (s = v, t = v + 1; t < nv; s++, t++) {
492
                V[s] = V[t];
493
            }
494

495
            nv--;
496

497
            /* reset error detection counter */
498
            count = 2 * nv;
499
        }
500
    }
501

502
    delete[] V;
503

504
    return true;
505
}
506

507
// -------------------------------------------------------------
508

509
QuasiDelaunayTriangulator::QuasiDelaunayTriangulator() = default;
510

511
bool QuasiDelaunayTriangulator::Triangulate()
512
{
513
    if (!EarClippingTriangulator::Triangulate()) {
514
        return false;  // no valid triangulation
515
    }
516

517
    // For each internal edge get the adjacent facets. When doing an edge swap we must update
518
    // this structure.
519
    std::map<std::pair<PointIndex, PointIndex>, std::vector<FacetIndex>> aEdge2Face;
520
    for (std::vector<MeshFacet>::iterator pI = _facets.begin(); pI != _facets.end(); ++pI) {
521
        for (int i = 0; i < 3; i++) {
522
            PointIndex ulPt0 = std::min<PointIndex>(pI->_aulPoints[i], pI->_aulPoints[(i + 1) % 3]);
523
            PointIndex ulPt1 = std::max<PointIndex>(pI->_aulPoints[i], pI->_aulPoints[(i + 1) % 3]);
524
            // ignore borderlines of the polygon
525
            if ((ulPt1 - ulPt0) % (_points.size() - 1) > 1) {
526
                aEdge2Face[std::pair<PointIndex, PointIndex>(ulPt0, ulPt1)].push_back(
527
                    pI - _facets.begin());
528
            }
529
        }
530
    }
531

532
    // fill up this list with all internal edges and perform swap edges until this list is empty
533
    std::list<std::pair<PointIndex, PointIndex>> aEdgeList;
534
    std::map<std::pair<PointIndex, PointIndex>, std::vector<FacetIndex>>::iterator pE;
535
    for (pE = aEdge2Face.begin(); pE != aEdge2Face.end(); ++pE) {
536
        aEdgeList.push_back(pE->first);
537
    }
538

539
    // to be sure to avoid an endless loop
540
    size_t uMaxIter = 5 * aEdge2Face.size();
541

542
    // Perform a swap edge where needed
543
    while (!aEdgeList.empty() && uMaxIter > 0) {
544
        // get the first edge and remove it from the list
545
        std::pair<PointIndex, PointIndex> aEdge = aEdgeList.front();
546
        aEdgeList.pop_front();
547
        uMaxIter--;
548

549
        // get the adjacent facets to this edge
550
        pE = aEdge2Face.find(aEdge);
551

552
        // this edge has been removed some iterations before
553
        if (pE == aEdge2Face.end()) {
554
            continue;
555
        }
556

557
        MeshFacet& rF1 = _facets[pE->second[0]];
558
        MeshFacet& rF2 = _facets[pE->second[1]];
559
        unsigned short side1 = rF1.Side(aEdge.first, aEdge.second);
560

561
        Base::Vector3f cP1 = _points[rF1._aulPoints[side1]];
562
        Base::Vector3f cP2 = _points[rF1._aulPoints[(side1 + 1) % 3]];
563
        Base::Vector3f cP3 = _points[rF1._aulPoints[(side1 + 2) % 3]];
564

565
        unsigned short side2 = rF2.Side(aEdge.first, aEdge.second);
566
        Base::Vector3f cP4 = _points[rF2._aulPoints[(side2 + 2) % 3]];
567

568
        MeshGeomFacet cT1(cP1, cP2, cP3);
569
        float fMax1 = cT1.MaximumAngle();
570
        MeshGeomFacet cT2(cP2, cP1, cP4);
571
        float fMax2 = cT2.MaximumAngle();
572
        MeshGeomFacet cT3(cP4, cP3, cP1);
573
        float fMax3 = cT3.MaximumAngle();
574
        MeshGeomFacet cT4(cP3, cP4, cP2);
575
        float fMax4 = cT4.MaximumAngle();
576

577
        float fMax12 = std::max<float>(fMax1, fMax2);
578
        float fMax34 = std::max<float>(fMax3, fMax4);
579

580
        // We must make sure that the two adjacent triangles builds a convex polygon, otherwise
581
        // the swap edge operation is illegal
582
        Base::Vector3f cU = cP2 - cP1;
583
        Base::Vector3f cV = cP4 - cP3;
584
        // build a helper plane through cP1 that must separate cP3 and cP4
585
        Base::Vector3f cN1 = (cU % cV) % cU;
586
        if (((cP3 - cP1) * cN1) * ((cP4 - cP1) * cN1) >= 0.0f) {
587
            continue;  // not convex
588
        }
589
        // build a helper plane through cP3 that must separate cP1 and cP2
590
        Base::Vector3f cN2 = (cU % cV) % cV;
591
        if (((cP1 - cP3) * cN2) * ((cP2 - cP3) * cN2) >= 0.0f) {
592
            continue;  // not convex
593
        }
594

595
        // ok, here we should perform a swap edge to minimize the maximum angle
596
        if (fMax12 > fMax34) {
597
            rF1._aulPoints[(side1 + 1) % 3] = rF2._aulPoints[(side2 + 2) % 3];
598
            rF2._aulPoints[(side2 + 1) % 3] = rF1._aulPoints[(side1 + 2) % 3];
599

600
            // adjust the edge list
601
            for (int i = 0; i < 3; i++) {
602
                std::map<std::pair<PointIndex, PointIndex>, std::vector<FacetIndex>>::iterator it;
603
                // first facet
604
                PointIndex ulPt0 =
605
                    std::min<PointIndex>(rF1._aulPoints[i], rF1._aulPoints[(i + 1) % 3]);
606
                PointIndex ulPt1 =
607
                    std::max<PointIndex>(rF1._aulPoints[i], rF1._aulPoints[(i + 1) % 3]);
608
                it = aEdge2Face.find(std::make_pair(ulPt0, ulPt1));
609
                if (it != aEdge2Face.end()) {
610
                    if (it->second[0] == pE->second[1]) {
611
                        it->second[0] = pE->second[0];
612
                    }
613
                    else if (it->second[1] == pE->second[1]) {
614
                        it->second[1] = pE->second[0];
615
                    }
616
                    aEdgeList.push_back(it->first);
617
                }
618

619
                // second facet
620
                ulPt0 = std::min<PointIndex>(rF2._aulPoints[i], rF2._aulPoints[(i + 1) % 3]);
621
                ulPt1 = std::max<PointIndex>(rF2._aulPoints[i], rF2._aulPoints[(i + 1) % 3]);
622
                it = aEdge2Face.find(std::make_pair(ulPt0, ulPt1));
623
                if (it != aEdge2Face.end()) {
624
                    if (it->second[0] == pE->second[0]) {
625
                        it->second[0] = pE->second[1];
626
                    }
627
                    else if (it->second[1] == pE->second[0]) {
628
                        it->second[1] = pE->second[1];
629
                    }
630
                    aEdgeList.push_back(it->first);
631
                }
632
            }
633

634
            // Now we must remove the edge and replace it through the new edge
635
            PointIndex ulPt0 = std::min<PointIndex>(rF1._aulPoints[(side1 + 1) % 3],
636
                                                    rF2._aulPoints[(side2 + 1) % 3]);
637
            PointIndex ulPt1 = std::max<PointIndex>(rF1._aulPoints[(side1 + 1) % 3],
638
                                                    rF2._aulPoints[(side2 + 1) % 3]);
639
            std::pair<PointIndex, PointIndex> aNewEdge = std::make_pair(ulPt0, ulPt1);
640
            aEdge2Face[aNewEdge] = pE->second;
641
            aEdge2Face.erase(pE);
642
        }
643
    }
644

645
    return true;
646
}
647

648
// -------------------------------------------------------------
649

650
namespace MeshCore
651
{
652
namespace Triangulation
653
{
654
struct Vertex2d_Less
655
{
656
    bool operator()(const Base::Vector3f& p, const Base::Vector3f& q) const
657
    {
658
        if (fabs(p.x - q.x) < MeshDefinitions::_fMinPointDistanceD1) {
659
            if (fabs(p.y - q.y) < MeshDefinitions::_fMinPointDistanceD1) {
660
                return false;
661
            }
662
            else {
663
                return p.y < q.y;
664
            }
665
        }
666
        else {
667
            return p.x < q.x;
668
        }
669
    }
670
};
671
struct Vertex2d_EqualTo
672
{
673
    bool operator()(const Base::Vector3f& p, const Base::Vector3f& q) const
674
    {
675
        if (fabs(p.x - q.x) < MeshDefinitions::_fMinPointDistanceD1
676
            && fabs(p.y - q.y) < MeshDefinitions::_fMinPointDistanceD1) {
677
            return true;
678
        }
679

680
        return false;
681
    }
682
};
683
}  // namespace Triangulation
684
}  // namespace MeshCore
685

686
DelaunayTriangulator::DelaunayTriangulator() = default;
687

688
bool DelaunayTriangulator::Triangulate()
689
{
690
    // before starting the triangulation we must make sure that all polygon
691
    // points are different
692
    std::vector<Base::Vector3f> aPoints = _points;
693
    // sort the points ascending x,y coordinates
694
    std::sort(aPoints.begin(), aPoints.end(), Triangulation::Vertex2d_Less());
695
    // if there are two adjacent points whose distance is less then an epsilon
696
    if (std::adjacent_find(aPoints.begin(), aPoints.end(), Triangulation::Vertex2d_EqualTo())
697
        < aPoints.end()) {
698
        return false;
699
    }
700

701
    _facets.clear();
702
    _triangles.clear();
703

704
    std::vector<Wm4::Vector2d> akVertex;
705
    akVertex.reserve(_points.size());
706
    for (const auto& point : _points) {
707
        akVertex.emplace_back(static_cast<double>(point.x), static_cast<double>(point.y));
708
    }
709

710
    Wm4::Delaunay2d del(static_cast<int>(akVertex.size()),
711
                        &(akVertex[0]),
712
                        0.001,
713
                        false,
714
                        Wm4::Query::QT_INT64);
715
    int iTQuantity = del.GetSimplexQuantity();
716
    std::vector<int> aiTVertex(static_cast<size_t>(3 * iTQuantity));
717

718
    bool succeeded = false;
719
    if (iTQuantity > 0) {
720
        size_t uiSize = static_cast<size_t>(3 * iTQuantity) * sizeof(int);
721
        Wm4::System::Memcpy(&(aiTVertex[0]), uiSize, del.GetIndices(), uiSize);
722

723
        // If H is the number of hull edges and N is the number of vertices,
724
        // then the triangulation must have 2*N-2-H triangles and 3*N-3-H
725
        // edges.
726
        int iEQuantity = 0;
727
        int* aiIndex = nullptr;
728
        del.GetHull(iEQuantity, aiIndex);
729
        int iUniqueVQuantity = del.GetUniqueVertexQuantity();
730
        int iTVerify = 2 * iUniqueVQuantity - 2 - iEQuantity;
731
        (void)iTVerify;  // avoid warning in release build
732
        succeeded = (iTVerify == iTQuantity);
733
        int iEVerify = 3 * iUniqueVQuantity - 3 - iEQuantity;
734
        (void)iEVerify;  // avoid warning about unused variable
735
        delete[] aiIndex;
736
    }
737

738
    MeshGeomFacet triangle;
739
    MeshFacet facet;
740
    for (int i = 0; i < iTQuantity; i++) {
741
        for (int j = 0; j < 3; j++) {
742
            size_t index = static_cast<size_t>(aiTVertex[static_cast<size_t>(3 * i + j)]);
743
            facet._aulPoints[j] = static_cast<PointIndex>(index);
744
            triangle._aclPoints[j].x = static_cast<float>(akVertex[index].X());
745
            triangle._aclPoints[j].y = static_cast<float>(akVertex[index].Y());
746
        }
747

748
        _triangles.push_back(triangle);
749
        _facets.push_back(facet);
750
    }
751

752
    return succeeded;
753
}
754

755
// -------------------------------------------------------------
756

757
FlatTriangulator::FlatTriangulator() = default;
758

759
bool FlatTriangulator::Triangulate()
760
{
761
    _newpoints.clear();
762
    // before starting the triangulation we must make sure that all polygon
763
    // points are different
764
    std::vector<Base::Vector3f> aPoints = ProjectToFitPlane();
765
    std::vector<Base::Vector3f> tmp = aPoints;
766
    // sort the points ascending x,y coordinates
767
    std::sort(tmp.begin(), tmp.end(), Triangulation::Vertex2d_Less());
768
    // if there are two adjacent points whose distance is less then an epsilon
769
    if (std::adjacent_find(tmp.begin(), tmp.end(), Triangulation::Vertex2d_EqualTo()) < tmp.end()) {
770
        return false;
771
    }
772

773
    _facets.clear();
774
    _triangles.clear();
775

776
    // Todo: Implement algorithm for constraint delaunay triangulation
777
    QuasiDelaunayTriangulator tria;
778
    tria.SetPolygon(this->GetPolygon());
779
    bool succeeded = tria.TriangulatePolygon();
780
    this->_facets = tria.GetFacets();
781
    this->_triangles = tria.GetTriangles();
782

783
    return succeeded;
784
}
785

786
void FlatTriangulator::PostProcessing(const std::vector<Base::Vector3f>&)
787
{}
788

789
// -------------------------------------------------------------
790

791
ConstraintDelaunayTriangulator::ConstraintDelaunayTriangulator(float area)
792
    : fMaxArea(area)
793
{
794
    // silent warning: -Wunused-private-field
795
    (void)fMaxArea;
796
}
797

798
bool ConstraintDelaunayTriangulator::Triangulate()
799
{
800
    _newpoints.clear();
801
    // before starting the triangulation we must make sure that all polygon
802
    // points are different
803
    std::vector<Base::Vector3f> aPoints = ProjectToFitPlane();
804
    std::vector<Base::Vector3f> tmp = aPoints;
805
    // sort the points ascending x,y coordinates
806
    std::sort(tmp.begin(), tmp.end(), Triangulation::Vertex2d_Less());
807
    // if there are two adjacent points whose distance is less then an epsilon
808
    if (std::adjacent_find(tmp.begin(), tmp.end(), Triangulation::Vertex2d_EqualTo()) < tmp.end()) {
809
        return false;
810
    }
811

812
    _facets.clear();
813
    _triangles.clear();
814

815
    // Todo: Implement algorithm for constraint delaunay triangulation
816
    QuasiDelaunayTriangulator tria;
817
    tria.SetPolygon(this->GetPolygon());
818
    bool succeeded = tria.TriangulatePolygon();
819
    this->_facets = tria.GetFacets();
820
    this->_triangles = tria.GetTriangles();
821

822
    return succeeded;
823
}
824

825
// -------------------------------------------------------------
826

827
#if 0
828
Triangulator::Triangulator(const MeshKernel& k, bool flat) : _kernel(k)
829
{
830
}
831

832
Triangulator::~Triangulator()
833
{
834
}
835

836
bool Triangulator::Triangulate()
837
{
838
    return false;
839
}
840

841
MeshGeomFacet Triangulator::GetTriangle(const MeshPointArray&,
842
                                        const MeshFacet& facet) const
843
{
844
    return MeshGeomFacet();
845
}
846

847
void Triangulator::PostProcessing(const std::vector<Base::Vector3f>&)
848
{
849
}
850

851
void Triangulator::Discard()
852
{
853
    AbstractPolygonTriangulator::Discard();
854
}
855

856
void Triangulator::Reset()
857
{
858
}
859
#endif
860