FreeCAD

1
/***************************************************************************
2
 *   Copyright (c) 2005 Imetric 3D GmbH                                    *
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

25
#ifndef _PreComp_
26
#include <algorithm>
27
#include <map>
28
#include <queue>
29
#endif
30

31
#include <boost/math/special_functions/fpclassify.hpp>
32

33
#include "Degeneration.h"
34
#include "Grid.h"
35
#include "Iterator.h"
36
#include "TopoAlgorithm.h"
37
#include "Triangulation.h"
38

39

40
using namespace MeshCore;
41

42
bool MeshEvalInvalids::Evaluate()
43
{
44
    const MeshFacetArray& rFaces = _rclMesh.GetFacets();
45
    for (const auto& it : rFaces) {
46
        if (!it.IsValid()) {
47
            return false;
48
        }
49
    }
50

51
    const MeshPointArray& rPoints = _rclMesh.GetPoints();
52
    for (const auto& it : rPoints) {
53
        if (!it.IsValid()) {
54
            return false;
55
        }
56
    }
57

58
    return true;
59
}
60

61
std::vector<FacetIndex> MeshEvalInvalids::GetIndices() const
62
{
63
    std::vector<FacetIndex> aInds;
64
    const MeshFacetArray& rFaces = _rclMesh.GetFacets();
65
    const MeshPointArray& rPoints = _rclMesh.GetPoints();
66
    FacetIndex ind = 0;
67
    for (MeshFacetArray::_TConstIterator it = rFaces.begin(); it != rFaces.end(); ++it, ind++) {
68
        if (!it->IsValid()) {
69
            aInds.push_back(ind);
70
        }
71
        else if (!rPoints[it->_aulPoints[0]].IsValid()) {
72
            aInds.push_back(ind);
73
        }
74
        else if (!rPoints[it->_aulPoints[1]].IsValid()) {
75
            aInds.push_back(ind);
76
        }
77
        else if (!rPoints[it->_aulPoints[2]].IsValid()) {
78
            aInds.push_back(ind);
79
        }
80
    }
81

82
    return aInds;
83
}
84

85
bool MeshFixInvalids::Fixup()
86
{
87
    _rclMesh.RemoveInvalids();
88
    return true;
89
}
90

91
// ----------------------------------------------------------------------
92

93
namespace MeshCore
94
{
95

96
using VertexIterator = MeshPointArray::_TConstIterator;
97
/*
98
 * When building up a mesh then usually the class MeshBuilder is used. This
99
 * class uses internally a std::set<MeshPoint> which uses the '<' operator of
100
 * MeshPoint to sort the points. Thus to be consistent (and avoid using the
101
 * '==' operator of MeshPoint) we use the same operator when comparing the
102
 * points in the function object.
103
 */
104
struct Vertex_EqualTo
105
{
106
    bool operator()(const VertexIterator& x, const VertexIterator& y) const
107
    {
108
        if ((*x) < (*y)) {
109
            return false;
110
        }
111
        else if ((*y) < (*x)) {
112
            return false;
113
        }
114
        return true;
115
    }
116
};
117

118
struct Vertex_Less
119
{
120
    bool operator()(const VertexIterator& x, const VertexIterator& y) const
121
    {
122
        return (*x) < (*y);
123
    }
124
};
125

126
}  // namespace MeshCore
127

128
bool MeshEvalDuplicatePoints::Evaluate()
129
{
130
    // get an const iterator to each vertex and sort them in ascending order by
131
    // their (x,y,z) coordinates
132
    const MeshPointArray& rPoints = _rclMesh.GetPoints();
133
    std::vector<VertexIterator> vertices;
134
    vertices.reserve(rPoints.size());
135
    for (MeshPointArray::_TConstIterator it = rPoints.begin(); it != rPoints.end(); ++it) {
136
        vertices.push_back(it);
137
    }
138

139
    // if there are two adjacent vertices which have the same coordinates
140
    std::sort(vertices.begin(), vertices.end(), Vertex_Less());
141
    if (std::adjacent_find(vertices.begin(), vertices.end(), Vertex_EqualTo()) < vertices.end()) {
142
        return false;
143
    }
144
    return true;
145
}
146

147
std::vector<PointIndex> MeshEvalDuplicatePoints::GetIndices() const
148
{
149
    // Note: We must neither use map or set to get duplicated indices because
150
    // the sort algorithms deliver different results compared to std::sort of
151
    // a vector.
152
    const MeshPointArray& rPoints = _rclMesh.GetPoints();
153
    std::vector<VertexIterator> vertices;
154
    vertices.reserve(rPoints.size());
155
    for (MeshPointArray::_TConstIterator it = rPoints.begin(); it != rPoints.end(); ++it) {
156
        vertices.push_back(it);
157
    }
158

159
    // if there are two adjacent vertices which have the same coordinates
160
    std::vector<PointIndex> aInds;
161
    Vertex_EqualTo pred;
162
    std::sort(vertices.begin(), vertices.end(), Vertex_Less());
163

164
    std::vector<VertexIterator>::iterator vt = vertices.begin();
165
    while (vt < vertices.end()) {
166
        // get first item which adjacent element has the same vertex
167
        vt = std::adjacent_find(vt, vertices.end(), pred);
168
        if (vt < vertices.end()) {
169
            ++vt;
170
            aInds.push_back(*vt - rPoints.begin());
171
        }
172
    }
173

174
    return aInds;
175
}
176

177
bool MeshFixDuplicatePoints::Fixup()
178
{
179
    // Note: We must neither use map or set to get duplicated indices because
180
    // the sort algorithms deliver different results compared to std::sort of
181
    // a vector.
182
    const MeshPointArray& rPoints = _rclMesh.GetPoints();
183
    std::vector<VertexIterator> vertices;
184
    vertices.reserve(rPoints.size());
185
    for (MeshPointArray::_TConstIterator it = rPoints.begin(); it != rPoints.end(); ++it) {
186
        vertices.push_back(it);
187
    }
188

189
    // get the indices of adjacent vertices which have the same coordinates
190
    std::sort(vertices.begin(), vertices.end(), Vertex_Less());
191

192
    Vertex_EqualTo pred;
193
    std::vector<VertexIterator>::iterator next = vertices.begin();
194
    std::map<PointIndex, PointIndex> mapPointIndex;
195
    std::vector<PointIndex> pointIndices;
196
    while (next < vertices.end()) {
197
        next = std::adjacent_find(next, vertices.end(), pred);
198
        if (next < vertices.end()) {
199
            std::vector<VertexIterator>::iterator first = next;
200
            PointIndex first_index = *first - rPoints.begin();
201
            ++next;
202
            while (next < vertices.end() && pred(*first, *next)) {
203
                PointIndex next_index = *next - rPoints.begin();
204
                mapPointIndex[next_index] = first_index;
205
                pointIndices.push_back(next_index);
206
                ++next;
207
            }
208
        }
209
    }
210

211
    // now set all facets to the correct index
212
    MeshFacetArray& rFacets = _rclMesh._aclFacetArray;
213
    for (auto& it : rFacets) {
214
        for (PointIndex& point : it._aulPoints) {
215
            std::map<PointIndex, PointIndex>::iterator pt = mapPointIndex.find(point);
216
            if (pt != mapPointIndex.end()) {
217
                point = pt->second;
218
            }
219
        }
220
    }
221

222
    // remove invalid indices
223
    _rclMesh.DeletePoints(pointIndices);
224
    _rclMesh.RebuildNeighbours();
225

226
    return true;
227
}
228

229
// ----------------------------------------------------------------------
230

231
bool MeshEvalNaNPoints::Evaluate()
232
{
233
    const MeshPointArray& rPoints = _rclMesh.GetPoints();
234
    for (const auto& it : rPoints) {
235
        if (boost::math::isnan(it.x) || boost::math::isnan(it.y) || boost::math::isnan(it.z)) {
236
            return false;
237
        }
238
    }
239

240
    return true;
241
}
242

243
std::vector<PointIndex> MeshEvalNaNPoints::GetIndices() const
244
{
245
    std::vector<PointIndex> aInds;
246
    const MeshPointArray& rPoints = _rclMesh.GetPoints();
247
    for (MeshPointArray::_TConstIterator it = rPoints.begin(); it != rPoints.end(); ++it) {
248
        if (boost::math::isnan(it->x) || boost::math::isnan(it->y) || boost::math::isnan(it->z)) {
249
            aInds.push_back(it - rPoints.begin());
250
        }
251
    }
252

253
    return aInds;
254
}
255

256
bool MeshFixNaNPoints::Fixup()
257
{
258
    std::vector<PointIndex> aInds;
259
    const MeshPointArray& rPoints = _rclMesh.GetPoints();
260
    for (MeshPointArray::_TConstIterator it = rPoints.begin(); it != rPoints.end(); ++it) {
261
        if (boost::math::isnan(it->x) || boost::math::isnan(it->y) || boost::math::isnan(it->z)) {
262
            aInds.push_back(it - rPoints.begin());
263
        }
264
    }
265

266
    // remove invalid indices
267
    _rclMesh.DeletePoints(aInds);
268
    _rclMesh.RebuildNeighbours();
269

270
    return true;
271
}
272

273
// ----------------------------------------------------------------------
274

275
namespace MeshCore
276
{
277

278
using FaceIterator = MeshFacetArray::_TConstIterator;
279
/*
280
 * The facet with the lowset index is regarded as 'less'.
281
 */
282
struct MeshFacet_Less
283
{
284
    bool operator()(const FaceIterator& x, const FaceIterator& y) const
285
    {
286
        PointIndex tmp {};
287
        PointIndex x0 = x->_aulPoints[0];
288
        PointIndex x1 = x->_aulPoints[1];
289
        PointIndex x2 = x->_aulPoints[2];
290
        PointIndex y0 = y->_aulPoints[0];
291
        PointIndex y1 = y->_aulPoints[1];
292
        PointIndex y2 = y->_aulPoints[2];
293

294
        if (x0 > x1) {
295
            tmp = x0;
296
            x0 = x1;
297
            x1 = tmp;
298
        }
299
        if (x0 > x2) {
300
            tmp = x0;
301
            x0 = x2;
302
            x2 = tmp;
303
        }
304
        if (x1 > x2) {
305
            tmp = x1;
306
            x1 = x2;
307
            x2 = tmp;
308
        }
309
        if (y0 > y1) {
310
            tmp = y0;
311
            y0 = y1;
312
            y1 = tmp;
313
        }
314
        if (y0 > y2) {
315
            tmp = y0;
316
            y0 = y2;
317
            y2 = tmp;
318
        }
319
        if (y1 > y2) {
320
            tmp = y1;
321
            y1 = y2;
322
            y2 = tmp;
323
        }
324

325
        if (x0 < y0) {
326
            return true;
327
        }
328
        else if (x0 > y0) {
329
            return false;
330
        }
331
        else if (x1 < y1) {
332
            return true;
333
        }
334
        else if (x1 > y1) {
335
            return false;
336
        }
337
        else if (x2 < y2) {
338
            return true;
339
        }
340
        else {
341
            return false;
342
        }
343
    }
344
};
345

346
}  // namespace MeshCore
347

348
/*
349
 * Two facets are equal if all its three point indices refer to the same
350
 * location in the point array of the mesh kernel they belong to.
351
 */
352
struct MeshFacet_EqualTo
353
{
354
    bool operator()(const FaceIterator& x, const FaceIterator& y) const
355
    {
356
        for (int i = 0; i < 3; i++) {
357
            if (x->_aulPoints[0] == y->_aulPoints[i]) {
358
                if (x->_aulPoints[1] == y->_aulPoints[(i + 1) % 3]
359
                    && x->_aulPoints[2] == y->_aulPoints[(i + 2) % 3]) {
360
                    return true;
361
                }
362
                else if (x->_aulPoints[1] == y->_aulPoints[(i + 2) % 3]
363
                         && x->_aulPoints[2] == y->_aulPoints[(i + 1) % 3]) {
364
                    return true;
365
                }
366
            }
367
        }
368

369
        return false;
370
    }
371
};
372

373
bool MeshEvalDuplicateFacets::Evaluate()
374
{
375
    std::set<FaceIterator, MeshFacet_Less> aFaces;
376
    const MeshFacetArray& rFaces = _rclMesh.GetFacets();
377
    for (MeshFacetArray::_TConstIterator it = rFaces.begin(); it != rFaces.end(); ++it) {
378
        std::pair<std::set<FaceIterator, MeshFacet_Less>::iterator, bool> pI = aFaces.insert(it);
379
        if (!pI.second) {
380
            return false;
381
        }
382
    }
383

384
    return true;
385
}
386

387
std::vector<FacetIndex> MeshEvalDuplicateFacets::GetIndices() const
388
{
389
#if 1
390
    const MeshFacetArray& rFacets = _rclMesh.GetFacets();
391
    std::vector<FaceIterator> faces;
392
    faces.reserve(rFacets.size());
393
    for (MeshFacetArray::_TConstIterator it = rFacets.begin(); it != rFacets.end(); ++it) {
394
        faces.push_back(it);
395
    }
396

397
    // if there are two adjacent faces which references the same vertices
398
    std::vector<FacetIndex> aInds;
399
    MeshFacet_EqualTo pred;
400
    std::sort(faces.begin(), faces.end(), MeshFacet_Less());
401

402
    std::vector<FaceIterator>::iterator ft = faces.begin();
403
    while (ft < faces.end()) {
404
        // get first item which adjacent element has the same face
405
        ft = std::adjacent_find(ft, faces.end(), pred);
406
        if (ft < faces.end()) {
407
            ++ft;
408
            aInds.push_back(*ft - rFacets.begin());
409
        }
410
    }
411

412
    return aInds;
413
#else
414
    std::vector<FacetIndex> aInds;
415
    const MeshFacetArray& rFaces = _rclMesh.GetFacets();
416
    FacetIndex uIndex = 0;
417

418
    // get all facets
419
    std::set<FaceIterator, MeshFacet_Less> aFaceSet;
420
    for (MeshFacetArray::_TConstIterator it = rFaces.begin(); it != rFaces.end(); ++it, uIndex++) {
421
        std::pair<std::set<FaceIterator, MeshFacet_Less>::iterator, bool> pI = aFaceSet.insert(it);
422
        if (!pI.second) {
423
            aInds.push_back(uIndex);
424
        }
425
    }
426

427
    return aInds;
428
#endif
429
}
430

431
bool MeshFixDuplicateFacets::Fixup()
432
{
433
    FacetIndex uIndex = 0;
434
    std::vector<FacetIndex> aRemoveFaces;
435
    const MeshFacetArray& rFaces = _rclMesh.GetFacets();
436

437
    // get all facets
438
    std::set<FaceIterator, MeshFacet_Less> aFaceSet;
439
    for (MeshFacetArray::_TConstIterator it = rFaces.begin(); it != rFaces.end(); ++it, uIndex++) {
440
        std::pair<std::set<FaceIterator, MeshFacet_Less>::iterator, bool> pI = aFaceSet.insert(it);
441
        if (!pI.second) {
442
            aRemoveFaces.push_back(uIndex);
443
        }
444
    }
445

446
    _rclMesh.DeleteFacets(aRemoveFaces);
447
    _rclMesh.RebuildNeighbours();  // needs to be done here
448

449
    return true;
450
}
451

452
// ----------------------------------------------------------------------
453

454
bool MeshEvalInternalFacets::Evaluate()
455
{
456
    _indices.clear();
457
    FacetIndex uIndex = 0;
458
    const MeshFacetArray& rFaces = _rclMesh.GetFacets();
459

460
    // get all facets
461
    std::set<FaceIterator, MeshFacet_Less> aFaceSet;
462
    MeshFacetArray::_TConstIterator first = rFaces.begin();
463
    for (MeshFacetArray::_TConstIterator it = rFaces.begin(); it != rFaces.end(); ++it, uIndex++) {
464
        std::pair<std::set<FaceIterator, MeshFacet_Less>::iterator, bool> pI = aFaceSet.insert(it);
465
        if (!pI.second) {
466
            // collect both elements
467
            _indices.push_back(*pI.first - first);
468
            _indices.push_back(uIndex);
469
        }
470
    }
471

472
    return _indices.empty();
473
}
474

475
// ----------------------------------------------------------------------
476

477
bool MeshEvalDegeneratedFacets::Evaluate()
478
{
479
    MeshFacetIterator it(_rclMesh);
480
    for (it.Init(); it.More(); it.Next()) {
481
        if (it->IsDegenerated(fEpsilon)) {
482
            return false;
483
        }
484
    }
485

486
    return true;
487
}
488

489
unsigned long MeshEvalDegeneratedFacets::CountEdgeTooSmall(float fMinEdgeLength) const
490
{
491
    MeshFacetIterator clFIter(_rclMesh);
492
    unsigned long k = 0;
493

494
    while (!clFIter.EndReached()) {
495
        for (int i = 0; i < 3; i++) {
496
            if (Base::Distance(clFIter->_aclPoints[i], clFIter->_aclPoints[(i + 1) % 3])
497
                < fMinEdgeLength) {
498
                k++;
499
            }
500
        }
501
        ++clFIter;
502
    }
503

504
    return k;
505
}
506

507
std::vector<FacetIndex> MeshEvalDegeneratedFacets::GetIndices() const
508
{
509
    std::vector<FacetIndex> aInds;
510
    MeshFacetIterator it(_rclMesh);
511
    for (it.Init(); it.More(); it.Next()) {
512
        if (it->IsDegenerated(fEpsilon)) {
513
            aInds.push_back(it.Position());
514
        }
515
    }
516

517
    return aInds;
518
}
519

520
bool MeshFixDegeneratedFacets::Fixup()
521
{
522
    MeshTopoAlgorithm cTopAlg(_rclMesh);
523

524
    MeshFacetIterator it(_rclMesh);
525
    for (it.Init(); it.More(); it.Next()) {
526
        if (it->IsDegenerated(fEpsilon)) {
527
            FacetIndex uId = it.Position();
528
            bool removed = cTopAlg.RemoveDegeneratedFacet(uId);
529
            if (removed) {
530
                // due to a modification of the array the iterator became invalid
531
                it.Set(uId - 1);
532
            }
533
        }
534
    }
535

536
    return true;
537
}
538

539
bool MeshRemoveNeedles::Fixup()
540
{
541
    using FaceEdge = std::pair<unsigned long, int>;  // (face, edge) pair
542
    using FaceEdgePriority = std::pair<float, FaceEdge>;
543

544
    MeshTopoAlgorithm topAlg(_rclMesh);
545
    MeshRefPointToFacets vf_it(_rclMesh);
546
    const MeshFacetArray& rclFAry = _rclMesh.GetFacets();
547
    const MeshPointArray& rclPAry = _rclMesh.GetPoints();
548
    rclFAry.ResetInvalid();
549
    rclPAry.ResetInvalid();
550
    rclPAry.ResetFlag(MeshPoint::VISIT);
551
    std::size_t facetCount = rclFAry.size();
552

553
    std::priority_queue<FaceEdgePriority, std::vector<FaceEdgePriority>, std::greater<>> todo;
554
    for (std::size_t index = 0; index < facetCount; index++) {
555
        const MeshFacet& facet = rclFAry[index];
556
        MeshGeomFacet tria(_rclMesh.GetFacet(facet));
557
        float perimeter = tria.Perimeter();
558
        float fMinLen = perimeter * fMinEdgeLength;
559
        for (int i = 0; i < 3; i++) {
560
            const Base::Vector3f& p1 = rclPAry[facet._aulPoints[i]];
561
            const Base::Vector3f& p2 = rclPAry[facet._aulPoints[(i + 1) % 3]];
562

563
            float distance = Base::Distance(p1, p2);
564
            if (distance < fMinLen) {
565
                unsigned long facetIndex = static_cast<unsigned long>(index);
566
                todo.push(std::make_pair(distance, std::make_pair(facetIndex, i)));
567
            }
568
        }
569
    }
570

571
    bool removedEdge = false;
572
    while (!todo.empty()) {
573
        FaceEdge faceedge = todo.top().second;
574
        todo.pop();
575

576
        // check if one of the face pairs was already processed
577
        if (!rclFAry[faceedge.first].IsValid()) {
578
            continue;
579
        }
580

581
        // the facet points may have changed, so check the current distance again
582
        const MeshFacet& facet = rclFAry[faceedge.first];
583
        MeshGeomFacet tria(_rclMesh.GetFacet(facet));
584
        float perimeter = tria.Perimeter();
585
        float fMinLen = perimeter * fMinEdgeLength;
586
        const Base::Vector3f& p1 = rclPAry[facet._aulPoints[faceedge.second]];
587
        const Base::Vector3f& p2 = rclPAry[facet._aulPoints[(faceedge.second + 1) % 3]];
588
        float distance = Base::Distance(p1, p2);
589
        if (distance >= fMinLen) {
590
            continue;
591
        }
592

593
        // collect the collapse-edge information
594
        EdgeCollapse ce;
595
        ce._fromPoint = rclFAry[faceedge.first]._aulPoints[faceedge.second];
596
        ce._toPoint = rclFAry[faceedge.first]._aulPoints[(faceedge.second + 1) % 3];
597

598
        ce._removeFacets.push_back(faceedge.first);
599
        FacetIndex neighbour = rclFAry[faceedge.first]._aulNeighbours[faceedge.second];
600
        if (neighbour != FACET_INDEX_MAX) {
601
            ce._removeFacets.push_back(neighbour);
602
        }
603

604
        std::set<FacetIndex> vf = vf_it[ce._fromPoint];
605
        vf.erase(faceedge.first);
606
        if (neighbour != FACET_INDEX_MAX) {
607
            vf.erase(neighbour);
608
        }
609
        ce._changeFacets.insert(ce._changeFacets.begin(), vf.begin(), vf.end());
610

611
        // get adjacent points
612
        std::set<PointIndex> vv;
613
        vv = vf_it.NeighbourPoints(ce._fromPoint);
614
        ce._adjacentFrom.insert(ce._adjacentFrom.begin(), vv.begin(), vv.end());
615
        vv = vf_it.NeighbourPoints(ce._toPoint);
616
        ce._adjacentTo.insert(ce._adjacentTo.begin(), vv.begin(), vv.end());
617

618
        if (topAlg.IsCollapseEdgeLegal(ce)) {
619
            topAlg.CollapseEdge(ce);
620
            for (auto it : ce._removeFacets) {
621
                vf_it.RemoveFacet(it);
622
            }
623
            for (auto it : ce._changeFacets) {
624
                vf_it.RemoveNeighbour(ce._fromPoint, it);
625
                vf_it.AddNeighbour(ce._toPoint, it);
626
            }
627
            removedEdge = true;
628
        }
629
    }
630

631
    if (removedEdge) {
632
        topAlg.Cleanup();
633
        _rclMesh.RebuildNeighbours();
634
    }
635

636
    return true;
637
}
638

639
// ----------------------------------------------------------------------
640

641
bool MeshFixCaps::Fixup()
642
{
643
    using FaceVertex = std::pair<unsigned long, int>;  // (face, vertex) pair
644
    using FaceVertexPriority = std::pair<float, FaceVertex>;
645

646
    MeshTopoAlgorithm topAlg(_rclMesh);
647
    const MeshFacetArray& rclFAry = _rclMesh.GetFacets();
648
    const MeshPointArray& rclPAry = _rclMesh.GetPoints();
649
    std::size_t facetCount = rclFAry.size();
650

651
    float fCosMaxAngle = static_cast<float>(cos(fMaxAngle));
652

653
    std::priority_queue<FaceVertexPriority, std::vector<FaceVertexPriority>, std::greater<>> todo;
654
    for (std::size_t index = 0; index < facetCount; index++) {
655
        for (int i = 0; i < 3; i++) {
656
            const MeshFacet& facet = rclFAry[index];
657
            const Base::Vector3f& p1 = rclPAry[facet._aulPoints[i]];
658
            const Base::Vector3f& p2 = rclPAry[facet._aulPoints[(i + 1) % 3]];
659
            const Base::Vector3f& p3 = rclPAry[facet._aulPoints[(i + 2) % 3]];
660
            Base::Vector3f dir1(p2 - p1);
661
            dir1.Normalize();
662
            Base::Vector3f dir2(p3 - p1);
663
            dir2.Normalize();
664

665
            float fCosAngle = dir1.Dot(dir2);
666
            if (fCosAngle < fCosMaxAngle) {
667
                unsigned long facetIndex = static_cast<unsigned long>(index);
668
                todo.push(std::make_pair(fCosAngle, std::make_pair(facetIndex, i)));
669
            }
670
        }
671
    }
672

673
    while (!todo.empty()) {
674
        FaceVertex facevertex = todo.top().second;
675
        todo.pop();
676

677
        // the facet points may have changed, so check the current distance again
678
        const MeshFacet& facet = rclFAry[facevertex.first];
679
        const Base::Vector3f& p1 = rclPAry[facet._aulPoints[facevertex.second]];
680
        const Base::Vector3f& p2 = rclPAry[facet._aulPoints[(facevertex.second + 1) % 3]];
681
        const Base::Vector3f& p3 = rclPAry[facet._aulPoints[(facevertex.second + 2) % 3]];
682
        Base::Vector3f dir1(p2 - p1);
683
        dir1.Normalize();
684
        Base::Vector3f dir2(p3 - p1);
685
        dir2.Normalize();
686

687
        // check that the criterion is still OK in case
688
        // an earlier edge-swap has an impact
689
        float fCosAngle = dir1.Dot(dir2);
690
        if (fCosAngle >= fCosMaxAngle) {
691
            continue;
692
        }
693

694
        // the triangle shouldn't be a needle, therefore the projection of the point with
695
        // the maximum angle must have a clear distance to the other corner points
696
        // as factor we choose a default value of 25% of the corresponding edge length
697
        Base::Vector3f p4 = p1.Perpendicular(p2, p3 - p2);
698
        float distP2P3 = Base::Distance(p2, p3);
699
        float distP2P4 = Base::Distance(p2, p4);
700
        float distP3P4 = Base::Distance(p3, p4);
701
        if (distP2P4 / distP2P3 < fSplitFactor || distP3P4 / distP2P3 < fSplitFactor) {
702
            continue;
703
        }
704

705
        FacetIndex facetpos = facevertex.first;
706
        FacetIndex neighbour = rclFAry[facetpos]._aulNeighbours[(facevertex.second + 1) % 3];
707
        if (neighbour != FACET_INDEX_MAX) {
708
            topAlg.SwapEdge(facetpos, neighbour);
709
        }
710
    }
711

712
    return true;
713
}
714

715
// ----------------------------------------------------------------------
716

717
bool MeshEvalDeformedFacets::Evaluate()
718
{
719
    float fCosMinAngle = cos(fMinAngle);
720
    float fCosMaxAngle = cos(fMaxAngle);
721

722
    MeshFacetIterator it(_rclMesh);
723
    for (it.Init(); it.More(); it.Next()) {
724
        if (it->IsDeformed(fCosMinAngle, fCosMaxAngle)) {
725
            return false;
726
        }
727
    }
728

729
    return true;
730
}
731

732
std::vector<FacetIndex> MeshEvalDeformedFacets::GetIndices() const
733
{
734
    float fCosMinAngle = cos(fMinAngle);
735
    float fCosMaxAngle = cos(fMaxAngle);
736

737
    std::vector<FacetIndex> aInds;
738
    MeshFacetIterator it(_rclMesh);
739
    for (it.Init(); it.More(); it.Next()) {
740
        if (it->IsDeformed(fCosMinAngle, fCosMaxAngle)) {
741
            aInds.push_back(it.Position());
742
        }
743
    }
744

745
    return aInds;
746
}
747

748
bool MeshFixDeformedFacets::Fixup()
749
{
750
    float fCosMinAngle = cos(fMinAngle);
751
    float fCosMaxAngle = cos(fMaxAngle);
752

753
    Base::Vector3f u, v;
754
    MeshTopoAlgorithm cTopAlg(_rclMesh);
755

756
    MeshFacetIterator it(_rclMesh);
757
    for (it.Init(); it.More(); it.Next()) {
758
        // possibly deformed but not degenerated
759
        if (!it->IsDegenerated(fEpsilon)) {
760
            // store the angles to avoid to compute twice
761
            float fCosAngles[3] = {0, 0, 0};
762
            bool done = false;
763

764
            for (int i = 0; i < 3; i++) {
765
                u = it->_aclPoints[(i + 1) % 3] - it->_aclPoints[i];
766
                v = it->_aclPoints[(i + 2) % 3] - it->_aclPoints[i];
767
                u.Normalize();
768
                v.Normalize();
769

770
                float fCosAngle = u * v;
771
                fCosAngles[i] = fCosAngle;
772
            }
773

774
            // first check for angle > 120 deg: in this case we swap with the opposite edge
775
            for (int i = 0; i < 3; i++) {
776
                float fCosAngle = fCosAngles[i];
777
                if (fCosAngle < fCosMaxAngle) {
778
                    const MeshFacet& face = it.GetReference();
779
                    FacetIndex uNeighbour = face._aulNeighbours[(i + 1) % 3];
780
                    if (uNeighbour != FACET_INDEX_MAX
781
                        && cTopAlg.ShouldSwapEdge(it.Position(), uNeighbour, fMaxSwapAngle)) {
782
                        cTopAlg.SwapEdge(it.Position(), uNeighbour);
783
                        done = true;
784
                    }
785
                    break;
786
                }
787
            }
788

789
            // we have swapped already
790
            if (done) {
791
                continue;
792
            }
793

794
            // now check for angle < 30 deg: in this case we swap with one of the edges the corner
795
            // is part of
796
            for (int j = 0; j < 3; j++) {
797
                float fCosAngle = fCosAngles[j];
798
                if (fCosAngle > fCosMinAngle) {
799
                    const MeshFacet& face = it.GetReference();
800

801
                    FacetIndex uNeighbour = face._aulNeighbours[j];
802
                    if (uNeighbour != FACET_INDEX_MAX
803
                        && cTopAlg.ShouldSwapEdge(it.Position(), uNeighbour, fMaxSwapAngle)) {
804
                        cTopAlg.SwapEdge(it.Position(), uNeighbour);
805
                        break;
806
                    }
807

808
                    uNeighbour = face._aulNeighbours[(j + 2) % 3];
809
                    if (uNeighbour != FACET_INDEX_MAX
810
                        && cTopAlg.ShouldSwapEdge(it.Position(), uNeighbour, fMaxSwapAngle)) {
811
                        cTopAlg.SwapEdge(it.Position(), uNeighbour);
812
                        break;
813
                    }
814
                }
815
            }
816
        }
817
    }
818

819
    return true;
820
}
821

822
// ----------------------------------------------------------------------
823

824
bool MeshFixMergeFacets::Fixup()
825
{
826
    MeshCore::MeshRefPointToPoints vv_it(_rclMesh);
827
    MeshCore::MeshRefPointToFacets vf_it(_rclMesh);
828
    unsigned long countPoints = _rclMesh.CountPoints();
829

830
    std::vector<MeshFacet> newFacets;
831
    newFacets.reserve(countPoints / 20);  // 5% should be sufficient
832

833
    MeshTopoAlgorithm topAlg(_rclMesh);
834
    for (unsigned long i = 0; i < countPoints; i++) {
835
        if (vv_it[i].size() == 3 && vf_it[i].size() == 3) {
836
            VertexCollapse vc;
837
            vc._point = i;
838
            const std::set<PointIndex>& adjPts = vv_it[i];
839
            vc._circumPoints.insert(vc._circumPoints.begin(), adjPts.begin(), adjPts.end());
840
            const std::set<FacetIndex>& adjFts = vf_it[i];
841
            vc._circumFacets.insert(vc._circumFacets.begin(), adjFts.begin(), adjFts.end());
842
            topAlg.CollapseVertex(vc);
843
        }
844
    }
845

846
    topAlg.Cleanup();
847
    return true;
848
}
849

850
// ----------------------------------------------------------------------
851

852
bool MeshEvalDentsOnSurface::Evaluate()
853
{
854
    this->indices.clear();
855
    MeshRefPointToFacets clPt2Facets(_rclMesh);
856
    const MeshPointArray& rPntAry = _rclMesh.GetPoints();
857
    MeshFacetArray::_TConstIterator f_beg = _rclMesh.GetFacets().begin();
858

859
    MeshGeomFacet rTriangle;
860
    Base::Vector3f tmp;
861
    unsigned long ctPoints = _rclMesh.CountPoints();
862
    for (unsigned long index = 0; index < ctPoints; index++) {
863
        std::vector<PointIndex> point;
864
        point.push_back(index);
865

866
        // get the local neighbourhood of the point
867
        std::set<PointIndex> nb = clPt2Facets.NeighbourPoints(point, 1);
868
        const std::set<FacetIndex>& faces = clPt2Facets[index];
869

870
        for (PointIndex pt : nb) {
871
            const MeshPoint& mp = rPntAry[pt];
872
            for (FacetIndex ft : faces) {
873
                // the point must not be part of the facet we test
874
                if (f_beg[ft]._aulPoints[0] == pt) {
875
                    continue;
876
                }
877
                if (f_beg[ft]._aulPoints[1] == pt) {
878
                    continue;
879
                }
880
                if (f_beg[ft]._aulPoints[2] == pt) {
881
                    continue;
882
                }
883
                // is the point projectable onto the facet?
884
                rTriangle = _rclMesh.GetFacet(f_beg[ft]);
885
                if (rTriangle.IntersectWithLine(mp, rTriangle.GetNormal(), tmp)) {
886
                    const std::set<FacetIndex>& f = clPt2Facets[pt];
887
                    this->indices.insert(this->indices.end(), f.begin(), f.end());
888
                    break;
889
                }
890
            }
891
        }
892
    }
893

894
    // remove duplicates
895
    std::sort(this->indices.begin(), this->indices.end());
896
    this->indices.erase(std::unique(this->indices.begin(), this->indices.end()),
897
                        this->indices.end());
898

899
    return this->indices.empty();
900
}
901

902
std::vector<FacetIndex> MeshEvalDentsOnSurface::GetIndices() const
903
{
904
    return this->indices;
905
}
906

907
/*
908
Forbidden is:
909
 + two facets share a common point but not a common edge
910

911
 Repair:
912
 + store the point indices which can be projected on a face
913
 + store the face indices on which a point can be projected
914
 + remove faces with an edge length smaller than a certain threshold (e.g. 0.01) from the stored
915
triangles or that reference one of the stored points
916
 + for this edge merge the two points
917
 + if a point of a face can be projected onto another face and they have a common point then split
918
the second face if the distance is under a certain threshold
919
 */
920
bool MeshFixDentsOnSurface::Fixup()
921
{
922
    MeshEvalDentsOnSurface eval(_rclMesh);
923
    if (!eval.Evaluate()) {
924
        std::vector<FacetIndex> inds = eval.GetIndices();
925
        _rclMesh.DeleteFacets(inds);
926
    }
927

928
    return true;
929
}
930

931
// ----------------------------------------------------------------------
932

933
bool MeshEvalFoldsOnSurface::Evaluate()
934
{
935
    this->indices.clear();
936
    const MeshFacetArray& rFAry = _rclMesh.GetFacets();
937
    unsigned long ct = 0;
938
    for (MeshFacetArray::const_iterator it = rFAry.begin(); it != rFAry.end(); ++it, ct++) {
939
        for (int i = 0; i < 3; i++) {
940
            FacetIndex n1 = it->_aulNeighbours[i];
941
            FacetIndex n2 = it->_aulNeighbours[(i + 1) % 3];
942
            Base::Vector3f v1 = _rclMesh.GetFacet(*it).GetNormal();
943
            if (n1 != FACET_INDEX_MAX && n2 != FACET_INDEX_MAX) {
944
                Base::Vector3f v2 = _rclMesh.GetFacet(n1).GetNormal();
945
                Base::Vector3f v3 = _rclMesh.GetFacet(n2).GetNormal();
946
                if (v2 * v3 > 0.0f) {
947
                    if (v1 * v2 < -0.1f && v1 * v3 < -0.1f) {
948
                        indices.push_back(n1);
949
                        indices.push_back(n2);
950
                        indices.push_back(ct);
951
                    }
952
                }
953
            }
954
        }
955
    }
956

957
    // remove duplicates
958
    std::sort(this->indices.begin(), this->indices.end());
959
    this->indices.erase(std::unique(this->indices.begin(), this->indices.end()),
960
                        this->indices.end());
961
    return this->indices.empty();
962
}
963

964
std::vector<FacetIndex> MeshEvalFoldsOnSurface::GetIndices() const
965
{
966
    return this->indices;
967
}
968

969
// ----------------------------------------------------------------------
970

971
bool MeshEvalFoldsOnBoundary::Evaluate()
972
{
973
    // remove all boundary facets with two open edges and where
974
    // the angle to the neighbour is more than 60 degree
975
    this->indices.clear();
976
    const MeshFacetArray& rFacAry = _rclMesh.GetFacets();
977
    for (MeshFacetArray::_TConstIterator it = rFacAry.begin(); it != rFacAry.end(); ++it) {
978
        if (it->CountOpenEdges() == 2) {
979
            for (FacetIndex nbIndex : it->_aulNeighbours) {
980
                if (nbIndex != FACET_INDEX_MAX) {
981
                    MeshGeomFacet f1 = _rclMesh.GetFacet(*it);
982
                    MeshGeomFacet f2 = _rclMesh.GetFacet(nbIndex);
983
                    float cos_angle = f1.GetNormal() * f2.GetNormal();
984
                    if (cos_angle <= 0.5f) {  // ~ 60 degree
985
                        indices.push_back(it - rFacAry.begin());
986
                    }
987
                }
988
            }
989
        }
990
    }
991

992
    return this->indices.empty();
993
}
994

995
std::vector<FacetIndex> MeshEvalFoldsOnBoundary::GetIndices() const
996
{
997
    return this->indices;
998
}
999

1000
bool MeshFixFoldsOnBoundary::Fixup()
1001
{
1002
    MeshEvalFoldsOnBoundary eval(_rclMesh);
1003
    if (!eval.Evaluate()) {
1004
        std::vector<FacetIndex> inds = eval.GetIndices();
1005
        _rclMesh.DeleteFacets(inds);
1006
    }
1007

1008
    return true;
1009
}
1010

1011
// ----------------------------------------------------------------------
1012

1013
bool MeshEvalFoldOversOnSurface::Evaluate()
1014
{
1015
    this->indices.clear();
1016
    const MeshCore::MeshFacetArray& facets = _rclMesh.GetFacets();
1017
    MeshCore::MeshFacetArray::_TConstIterator f_it, f_beg = facets.begin(), f_end = facets.end();
1018

1019
    Base::Vector3f n1, n2;
1020
    for (f_it = facets.begin(); f_it != f_end; ++f_it) {
1021
        for (int i = 0; i < 3; i++) {
1022
            FacetIndex index1 = f_it->_aulNeighbours[i];
1023
            FacetIndex index2 = f_it->_aulNeighbours[(i + 1) % 3];
1024
            if (index1 != FACET_INDEX_MAX && index2 != FACET_INDEX_MAX) {
1025
                // if the topology is correct but the normals flip from
1026
                // two neighbours we have a fold
1027
                if (f_it->HasSameOrientation(f_beg[index1])
1028
                    && f_it->HasSameOrientation(f_beg[index2])) {
1029
                    n1 = _rclMesh.GetFacet(index1).GetNormal();
1030
                    n2 = _rclMesh.GetFacet(index2).GetNormal();
1031
                    if (n1 * n2 < -0.5f) {  // angle > 120 deg
1032
                        this->indices.push_back(f_it - f_beg);
1033
                        break;
1034
                    }
1035
                }
1036
            }
1037
        }
1038
    }
1039

1040
    return this->indices.empty();
1041
}
1042

1043
// ----------------------------------------------------------------
1044

1045
bool MeshEvalBorderFacet::Evaluate()
1046
{
1047
    const MeshCore::MeshFacetArray& facets = _rclMesh.GetFacets();
1048
    MeshCore::MeshFacetArray::_TConstIterator f_it, f_beg = facets.begin(), f_end = facets.end();
1049
    MeshCore::MeshRefPointToPoints vv_it(_rclMesh);
1050
    MeshCore::MeshRefPointToFacets vf_it(_rclMesh);
1051

1052
    for (f_it = facets.begin(); f_it != f_end; ++f_it) {
1053
        bool ok = true;
1054
        for (PointIndex index : f_it->_aulPoints) {
1055
            if (vv_it[index].size() == vf_it[index].size()) {
1056
                ok = false;
1057
                break;
1058
            }
1059
        }
1060

1061
        if (ok) {
1062
            _facets.push_back(f_it - f_beg);
1063
        }
1064
    }
1065

1066
    return _facets.empty();
1067
}
1068

1069
// ----------------------------------------------------------------------
1070

1071
bool MeshEvalRangeFacet::Evaluate()
1072
{
1073
    const MeshFacetArray& rFaces = _rclMesh.GetFacets();
1074
    FacetIndex ulCtFacets = rFaces.size();
1075

1076
    for (const auto& it : rFaces) {
1077
        for (FacetIndex nbFacet : it._aulNeighbours) {
1078
            if ((nbFacet >= ulCtFacets) && (nbFacet < FACET_INDEX_MAX)) {
1079
                return false;
1080
            }
1081
        }
1082
    }
1083

1084
    return true;
1085
}
1086

1087
std::vector<FacetIndex> MeshEvalRangeFacet::GetIndices() const
1088
{
1089
    std::vector<FacetIndex> aInds;
1090
    const MeshFacetArray& rFaces = _rclMesh.GetFacets();
1091
    FacetIndex ulCtFacets = rFaces.size();
1092

1093
    FacetIndex ind = 0;
1094
    for (MeshFacetArray::_TConstIterator it = rFaces.begin(); it != rFaces.end(); ++it, ind++) {
1095
        for (FacetIndex nbIndex : it->_aulNeighbours) {
1096
            if ((nbIndex >= ulCtFacets) && (nbIndex < FACET_INDEX_MAX)) {
1097
                aInds.push_back(ind);
1098
                break;
1099
            }
1100
        }
1101
    }
1102

1103
    return aInds;
1104
}
1105

1106
bool MeshFixRangeFacet::Fixup()
1107
{
1108
    _rclMesh.RebuildNeighbours();
1109
    return true;
1110
}
1111

1112
// ----------------------------------------------------------------------
1113

1114
bool MeshEvalRangePoint::Evaluate()
1115
{
1116
    const MeshFacetArray& rFaces = _rclMesh.GetFacets();
1117
    PointIndex ulCtPoints = _rclMesh.CountPoints();
1118

1119
    for (const auto& it : rFaces) {
1120
        if (std::find_if(it._aulPoints,
1121
                         it._aulPoints + 3,
1122
                         [ulCtPoints](PointIndex i) {
1123
                             return i >= ulCtPoints;
1124
                         })
1125
            < it._aulPoints + 3) {
1126
            return false;
1127
        }
1128
    }
1129

1130
    return true;
1131
}
1132

1133
std::vector<PointIndex> MeshEvalRangePoint::GetIndices() const
1134
{
1135
    std::vector<PointIndex> aInds;
1136
    const MeshFacetArray& rFaces = _rclMesh.GetFacets();
1137
    PointIndex ulCtPoints = _rclMesh.CountPoints();
1138

1139
    PointIndex ind = 0;
1140
    for (MeshFacetArray::_TConstIterator it = rFaces.begin(); it != rFaces.end(); ++it, ind++) {
1141
        if (std::find_if(it->_aulPoints,
1142
                         it->_aulPoints + 3,
1143
                         [ulCtPoints](PointIndex i) {
1144
                             return i >= ulCtPoints;
1145
                         })
1146
            < it->_aulPoints + 3) {
1147
            aInds.push_back(ind);
1148
        }
1149
    }
1150

1151
    return aInds;
1152
}
1153

1154
bool MeshFixRangePoint::Fixup()
1155
{
1156
    MeshEvalRangePoint eval(_rclMesh);
1157
    if (_rclMesh.CountPoints() == 0) {
1158
        // if no points are there but facets then the whole mesh can be cleared
1159
        _rclMesh.Clear();
1160
    }
1161
    else {
1162
        // facets with point indices out of range cannot be directly deleted because
1163
        // 'DeleteFacets' will segfault. But setting all point indices to 0 works.
1164
        std::vector<PointIndex> invalid = eval.GetIndices();
1165
        if (!invalid.empty()) {
1166
            for (PointIndex it : invalid) {
1167
                _rclMesh.SetFacetPoints(it, 0, 0, 0);
1168
            }
1169

1170
            _rclMesh.DeleteFacets(invalid);
1171
        }
1172
    }
1173
    return true;
1174
}
1175

1176
// ----------------------------------------------------------------------
1177

1178
bool MeshEvalCorruptedFacets::Evaluate()
1179
{
1180
    const MeshFacetArray& rFaces = _rclMesh.GetFacets();
1181

1182
    for (const auto& it : rFaces) {
1183
        // duplicated point indices
1184
        if (it.IsDegenerated()) {
1185
            return false;
1186
        }
1187
    }
1188

1189
    return true;
1190
}
1191

1192
std::vector<FacetIndex> MeshEvalCorruptedFacets::GetIndices() const
1193
{
1194
    std::vector<FacetIndex> aInds;
1195
    const MeshFacetArray& rFaces = _rclMesh.GetFacets();
1196
    FacetIndex ind = 0;
1197

1198
    for (MeshFacetArray::_TConstIterator it = rFaces.begin(); it != rFaces.end(); ++it, ind++) {
1199
        if (it->IsDegenerated()) {
1200
            aInds.push_back(ind);
1201
        }
1202
    }
1203

1204
    return aInds;
1205
}
1206

1207
bool MeshFixCorruptedFacets::Fixup()
1208
{
1209
    MeshTopoAlgorithm cTopAlg(_rclMesh);
1210

1211
    MeshFacetIterator it(_rclMesh);
1212
    for (it.Init(); it.More(); it.Next()) {
1213
        if (it.GetReference().IsDegenerated()) {
1214
            unsigned long uId = it.Position();
1215
            bool removed = cTopAlg.RemoveCorruptedFacet(uId);
1216
            if (removed) {
1217
                // due to a modification of the array the iterator became invalid
1218
                it.Set(uId - 1);
1219
            }
1220
        }
1221
    }
1222

1223
    return true;
1224
}
1225

1226
// ----------------------------------------------------------------------
1227

1228
bool MeshEvalPointOnEdge::Evaluate()
1229
{
1230
    MeshFacetGrid facetGrid(_rclMesh);
1231
    const MeshPointArray& points = _rclMesh.GetPoints();
1232
    const MeshFacetArray& facets = _rclMesh.GetFacets();
1233

1234
    auto IsPointOnEdge = [&points](PointIndex idx, const MeshFacet& facet) {
1235
        // point must not be a corner of the facet
1236
        if (!facet.HasPoint(idx)) {
1237
            for (int i = 0; i < 3; i++) {
1238
                MeshGeomEdge edge;
1239
                edge._aclPoints[0] = points[facet._aulPoints[i]];
1240
                edge._aclPoints[1] = points[facet._aulPoints[(i + 1) % 3]];
1241

1242
                if (edge.GetBoundBox().IsInBox(points[idx])) {
1243
                    if (edge.IsPointOf(points[idx], 0.001f)) {
1244
                        return true;
1245
                    }
1246
                }
1247
            }
1248
        }
1249
        return false;
1250
    };
1251

1252
    PointIndex maxPoints = _rclMesh.CountPoints();
1253
    for (PointIndex i = 0; i < maxPoints; i++) {
1254
        std::vector<FacetIndex> elements;
1255
        facetGrid.GetElements(points[i], elements);
1256

1257
        for (const auto& it : elements) {
1258
            const MeshFacet& face = facets[it];
1259
            if (IsPointOnEdge(i, face)) {
1260
                pointsIndices.push_back(i);
1261
                if (face.HasOpenEdge()) {
1262
                    facetsIndices.push_back(it);
1263
                }
1264
            }
1265
        }
1266
    }
1267
    return pointsIndices.empty();
1268
}
1269

1270
std::vector<PointIndex> MeshEvalPointOnEdge::GetPointIndices() const
1271
{
1272
    return pointsIndices;
1273
}
1274

1275
std::vector<FacetIndex> MeshEvalPointOnEdge::GetFacetIndices() const
1276
{
1277
    return facetsIndices;
1278
}
1279

1280
bool MeshFixPointOnEdge::Fixup()
1281
{
1282
    MeshEvalPointOnEdge eval(_rclMesh);
1283
    eval.Evaluate();
1284
    std::vector<PointIndex> pointsIndices = eval.GetPointIndices();
1285
    std::vector<FacetIndex> facetsIndices = eval.GetFacetIndices();
1286

1287
    if (!pointsIndices.empty()) {
1288
        if (fillBoundary) {
1289
            MarkBoundaries(facetsIndices);
1290
        }
1291

1292
        _rclMesh.DeletePoints(pointsIndices);
1293

1294
        if (fillBoundary) {
1295
            std::list<std::vector<PointIndex>> borderList;
1296
            FindBoundaries(borderList);
1297
            if (!borderList.empty()) {
1298
                FillBoundaries(borderList);
1299
            }
1300
        }
1301
    }
1302

1303
    return true;
1304
}
1305

1306
void MeshFixPointOnEdge::MarkBoundaries(const std::vector<FacetIndex>& facetsIndices)
1307
{
1308
    MeshAlgorithm meshalg(_rclMesh);
1309
    meshalg.ResetFacetFlag(MeshFacet::TMP0);
1310
    meshalg.SetFacetsFlag(facetsIndices, MeshFacet::TMP0);
1311
}
1312

1313
void MeshFixPointOnEdge::FindBoundaries(std::list<std::vector<PointIndex>>& borderList)
1314
{
1315
    std::vector<FacetIndex> tmp;
1316
    MeshAlgorithm meshalg(_rclMesh);
1317
    meshalg.GetFacetsFlag(tmp, MeshFacet::TMP0);
1318

1319
    if (!tmp.empty()) {
1320
        meshalg.GetFacetsBorders(tmp, borderList);
1321
    }
1322
}
1323

1324
void MeshFixPointOnEdge::FillBoundaries(const std::list<std::vector<PointIndex>>& borderList)
1325
{
1326
    FlatTriangulator tria;
1327
    tria.SetVerifier(new MeshCore::TriangulationVerifierV2);
1328
    MeshTopoAlgorithm topalg(_rclMesh);
1329
    std::list<std::vector<PointIndex>> failed;
1330
    topalg.FillupHoles(1, tria, borderList, failed);
1331
}
1332