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
#endif
28

29
#include <Base/Console.h>
30
#include <Base/Sequencer.h>
31

32
#include "Algorithm.h"
33
#include "Approximation.h"
34
#include "Elements.h"
35
#include "Grid.h"
36
#include "Iterator.h"
37
#include "Triangulation.h"
38

39

40
using namespace MeshCore;
41
using Base::BoundBox2d;
42
using Base::BoundBox3f;
43
using Base::Polygon2d;
44

45

46
bool MeshAlgorithm::IsVertexVisible(const Base::Vector3f& rcVertex,
47
                                    const Base::Vector3f& rcView,
48
                                    const MeshFacetGrid& rclGrid) const
49
{
50
    const float fMaxDistance = 0.001F;
51
    Base::Vector3f cDirection = rcVertex - rcView;
52
    float fDistance = cDirection.Length();
53
    Base::Vector3f cIntsct;
54
    FacetIndex uInd {};
55

56
    // search for the nearest facet to rcView in direction to rcVertex
57
    if (NearestFacetOnRay(rcView, cDirection, /*1.2f*fDistance,*/ rclGrid, cIntsct, uInd)) {
58
        // now check if the facet overlays the point
59
        float fLen = Base::Distance(rcView, cIntsct);
60
        if (fLen < fDistance) {
61
            // is it the same point?
62
            if (Base::Distance(rcVertex, cIntsct) > fMaxDistance) {
63
                // ok facet overlays the vertex
64
                return false;
65
            }
66
        }
67
    }
68

69
    return true;  // no facet between the two points
70
}
71

72
bool MeshAlgorithm::NearestFacetOnRay(const Base::Vector3f& rclPt,
73
                                      const Base::Vector3f& rclDir,
74
                                      Base::Vector3f& rclRes,
75
                                      FacetIndex& rulFacet) const
76
{
77
    return NearestFacetOnRay(rclPt, rclDir, Mathf::PI, rclRes, rulFacet);
78
}
79

80
bool MeshAlgorithm::NearestFacetOnRay(const Base::Vector3f& rclPt,
81
                                      const Base::Vector3f& rclDir,
82
                                      float fMaxAngle,
83
                                      Base::Vector3f& rclRes,
84
                                      FacetIndex& rulFacet) const
85
{
86
    Base::Vector3f clProj;
87
    Base::Vector3f clRes;
88
    bool bSol = false;
89
    FacetIndex ulInd = 0;
90

91
    // slow execution with no grid
92
    MeshFacetIterator clFIter(_rclMesh);
93
    for (clFIter.Init(); clFIter.More(); clFIter.Next()) {
94
        if (clFIter->Foraminate(rclPt, rclDir, clRes, fMaxAngle)) {
95
            if (!bSol) {
96
                // first solution
97
                bSol = true;
98
                clProj = clRes;
99
                ulInd = clFIter.Position();
100
            }
101
            else {
102
                // is closer to the point
103
                if ((clRes - rclPt).Length() < (clProj - rclPt).Length()) {
104
                    clProj = clRes;
105
                    ulInd = clFIter.Position();
106
                }
107
            }
108
        }
109
    }
110

111
    if (bSol) {
112
        rclRes = clProj;
113
        rulFacet = ulInd;
114
    }
115

116
    return bSol;
117
}
118

119
bool MeshAlgorithm::NearestFacetOnRay(const Base::Vector3f& rclPt,
120
                                      const Base::Vector3f& rclDir,
121
                                      const MeshFacetGrid& rclGrid,
122
                                      Base::Vector3f& rclRes,
123
                                      FacetIndex& rulFacet) const
124
{
125
    std::vector<FacetIndex> aulFacets;
126
    MeshGridIterator clGridIter(rclGrid);
127

128
    if (clGridIter.InitOnRay(rclPt, rclDir, aulFacets)) {
129
        if (!RayNearestField(rclPt, rclDir, aulFacets, rclRes, rulFacet)) {
130
            aulFacets.clear();
131
            while (clGridIter.NextOnRay(aulFacets)) {
132
                if (RayNearestField(rclPt, rclDir, aulFacets, rclRes, rulFacet)) {
133
                    return true;
134
                }
135
            }
136
        }
137
        else {
138
            return true;
139
        }
140
    }
141

142
    return false;
143
}
144

145
bool MeshAlgorithm::NearestFacetOnRay(const Base::Vector3f& rclPt,
146
                                      const Base::Vector3f& rclDir,
147
                                      float fMaxSearchArea,
148
                                      const MeshFacetGrid& rclGrid,
149
                                      Base::Vector3f& rclRes,
150
                                      FacetIndex& rulFacet) const
151
{
152
    const float fMaxAngle = 1.75F;
153
    std::vector<FacetIndex> aulFacets;
154
    MeshGridIterator clGridIter(rclGrid);
155

156
    if (clGridIter.InitOnRay(rclPt, rclDir, fMaxSearchArea, aulFacets)) {
157
        if (!RayNearestField(rclPt, rclDir, aulFacets, rclRes, rulFacet, fMaxAngle)) {
158
            aulFacets.clear();
159
            while (clGridIter.NextOnRay(aulFacets)) {
160
                if (RayNearestField(rclPt, rclDir, aulFacets, rclRes, rulFacet, fMaxAngle)) {
161
                    return true;
162
                }
163
            }
164
        }
165
        else {
166
            return true;
167
        }
168
    }
169

170
    return false;
171
}
172

173
bool MeshAlgorithm::NearestFacetOnRay(const Base::Vector3f& rclPt,
174
                                      const Base::Vector3f& rclDir,
175
                                      const std::vector<FacetIndex>& raulFacets,
176
                                      Base::Vector3f& rclRes,
177
                                      FacetIndex& rulFacet) const
178
{
179
    Base::Vector3f clProj;
180
    Base::Vector3f clRes;
181
    bool bSol = false;
182
    FacetIndex ulInd = 0;
183

184
    for (FacetIndex index : raulFacets) {
185
        MeshGeomFacet rclSFacet = _rclMesh.GetFacet(index);
186
        if (rclSFacet.Foraminate(rclPt, rclDir, clRes)) {
187
            if (!bSol) {  // first solution
188
                bSol = true;
189
                clProj = clRes;
190
                ulInd = index;
191
            }
192
            else {  // is closer to the point
193
                if ((clRes - rclPt).Length() < (clProj - rclPt).Length()) {
194
                    clProj = clRes;
195
                    ulInd = index;
196
                }
197
            }
198
        }
199
    }
200

201
    if (bSol) {
202
        rclRes = clProj;
203
        rulFacet = ulInd;
204
    }
205

206
    return bSol;
207
}
208

209
bool MeshAlgorithm::RayNearestField(const Base::Vector3f& rclPt,
210
                                    const Base::Vector3f& rclDir,
211
                                    const std::vector<FacetIndex>& raulFacets,
212
                                    Base::Vector3f& rclRes,
213
                                    FacetIndex& rulFacet,
214
                                    float /*fMaxAngle*/) const
215
{
216
    Base::Vector3f clProj, clRes;
217
    bool bSol = false;
218
    FacetIndex ulInd = 0;
219

220
    for (FacetIndex index : raulFacets) {
221
        if (_rclMesh.GetFacet(index).Foraminate(rclPt, rclDir, clRes /*, fMaxAngle*/)) {
222
            if (!bSol) {  // first solution
223
                bSol = true;
224
                clProj = clRes;
225
                ulInd = index;
226
            }
227
            else {  // is closer to the point
228
                if ((clRes - rclPt).Length() < (clProj - rclPt).Length()) {
229
                    clProj = clRes;
230
                    ulInd = index;
231
                }
232
            }
233
        }
234
    }
235

236
    if (bSol) {
237
        rclRes = clProj;
238
        rulFacet = ulInd;
239
    }
240

241
    return bSol;
242
}
243

244
bool MeshAlgorithm::FirstFacetToVertex(const Base::Vector3f& rPt,
245
                                       float fMaxDistance,
246
                                       const MeshFacetGrid& rGrid,
247
                                       FacetIndex& uIndex) const
248
{
249
    const float fEps = 0.001f;
250

251
    bool found = false;
252
    std::vector<FacetIndex> facets;
253

254
    // get the facets of the grid the point lies into
255
    rGrid.GetElements(rPt, facets);
256

257
    // Check all facets inside the grid if the point is part of it
258
    for (FacetIndex facet : facets) {
259
        MeshGeomFacet cFacet = this->_rclMesh.GetFacet(facet);
260
        if (cFacet.IsPointOfFace(rPt, fMaxDistance)) {
261
            found = true;
262
            uIndex = facet;
263
            break;
264
        }
265
        else {
266
            // if not then check the distance to the border of the triangle
267
            Base::Vector3f res;
268
            float fDist {};
269
            unsigned short uSide {};
270
            cFacet.ProjectPointToPlane(rPt, res);
271
            cFacet.NearestEdgeToPoint(res, fDist, uSide);
272
            if (fDist < fEps) {
273
                found = true;
274
                uIndex = facet;
275
                break;
276
            }
277
        }
278
    }
279

280
    return found;
281
}
282

283
float MeshAlgorithm::GetAverageEdgeLength() const
284
{
285
    float fLen = 0.0f;
286
    MeshFacetIterator cF(_rclMesh);
287
    for (cF.Init(); cF.More(); cF.Next()) {
288
        for (int i = 0; i < 3; i++) {
289
            fLen += Base::Distance(cF->_aclPoints[i], cF->_aclPoints[(i + 1) % 3]);
290
        }
291
    }
292

293
    fLen = fLen / (3.0f * _rclMesh.CountFacets());
294
    return fLen;
295
}
296

297
float MeshAlgorithm::GetMinimumEdgeLength() const
298
{
299
    float fLen = FLOAT_MAX;
300
    MeshFacetIterator cF(_rclMesh);
301
    for (cF.Init(); cF.More(); cF.Next()) {
302
        for (int i = 0; i < 3; i++) {
303
            fLen = std::min(fLen, Base::Distance(cF->_aclPoints[i], cF->_aclPoints[(i + 1) % 3]));
304
        }
305
    }
306

307
    return fLen;
308
}
309

310
float MeshAlgorithm::GetMaximumEdgeLength() const
311
{
312
    float fLen = 0.0f;
313
    MeshFacetIterator cF(_rclMesh);
314
    for (cF.Init(); cF.More(); cF.Next()) {
315
        for (int i = 0; i < 3; i++) {
316
            fLen = std::max(fLen, Base::Distance(cF->_aclPoints[i], cF->_aclPoints[(i + 1) % 3]));
317
        }
318
    }
319

320
    return fLen;
321
}
322

323
Base::Vector3f MeshAlgorithm::GetGravityPoint() const
324
{
325
    Base::Vector3f center;
326
    MeshPointIterator cP(_rclMesh);
327
    for (cP.Init(); cP.More(); cP.Next()) {
328
        center += *cP;
329
    }
330

331
    return center / static_cast<float>(_rclMesh.CountPoints());
332
}
333

334
void MeshAlgorithm::GetMeshBorders(std::list<std::vector<Base::Vector3f>>& rclBorders) const
335
{
336
    std::vector<FacetIndex> aulAllFacets(_rclMesh.CountFacets());
337
    FacetIndex k = 0;
338
    for (FacetIndex& index : aulAllFacets) {
339
        index = k++;
340
    }
341

342
    GetFacetBorders(aulAllFacets, rclBorders);
343
}
344

345
void MeshAlgorithm::GetMeshBorders(std::list<std::vector<PointIndex>>& rclBorders) const
346
{
347
    std::vector<FacetIndex> aulAllFacets(_rclMesh.CountFacets());
348
    FacetIndex k = 0;
349
    for (FacetIndex& index : aulAllFacets) {
350
        index = k++;
351
    }
352

353
    GetFacetBorders(aulAllFacets, rclBorders, true);
354
}
355

356
void MeshAlgorithm::GetFacetBorders(const std::vector<FacetIndex>& raulInd,
357
                                    std::list<std::vector<Base::Vector3f>>& rclBorders) const
358
{
359
    const MeshPointArray& rclPAry = _rclMesh._aclPointArray;
360
    std::list<std::vector<PointIndex>> aulBorders;
361

362
    GetFacetBorders(raulInd, aulBorders, true);
363
    for (const auto& border : aulBorders) {
364
        std::vector<Base::Vector3f> boundary;
365
        boundary.reserve(border.size());
366

367
        for (PointIndex jt : border) {
368
            boundary.push_back(rclPAry[jt]);
369
        }
370

371
        rclBorders.push_back(boundary);
372
    }
373
}
374

375
void MeshAlgorithm::GetFacetBorders(const std::vector<FacetIndex>& raulInd,
376
                                    std::list<std::vector<PointIndex>>& rclBorders,
377
                                    bool ignoreOrientation) const
378
{
379
    const MeshFacetArray& rclFAry = _rclMesh._aclFacetArray;
380

381
    // mark all facets that are in the indices list
382
    ResetFacetFlag(MeshFacet::VISIT);
383
    for (FacetIndex it : raulInd) {
384
        rclFAry[it].SetFlag(MeshFacet::VISIT);
385
    }
386

387
    // collect all boundary edges (unsorted)
388
    std::list<std::pair<PointIndex, PointIndex>> aclEdges;
389
    for (FacetIndex it : raulInd) {
390
        const MeshFacet& rclFacet = rclFAry[it];
391
        for (unsigned short i = 0; i < 3; i++) {
392
            FacetIndex ulNB = rclFacet._aulNeighbours[i];
393
            if (ulNB != FACET_INDEX_MAX) {
394
                if (rclFAry[ulNB].IsFlag(MeshFacet::VISIT)) {
395
                    continue;
396
                }
397
            }
398

399
            aclEdges.push_back(rclFacet.GetEdge(i));
400
        }
401
    }
402

403
    if (aclEdges.empty()) {
404
        return;  // no borders found (=> solid)
405
    }
406

407
    // search for edges in the unsorted list
408
    PointIndex ulFirst {}, ulLast {};
409
    std::list<PointIndex> clBorder;
410
    ulFirst = aclEdges.begin()->first;
411
    ulLast = aclEdges.begin()->second;
412

413
    aclEdges.erase(aclEdges.begin());
414
    clBorder.push_back(ulFirst);
415
    clBorder.push_back(ulLast);
416

417
    while (!aclEdges.empty()) {
418
        // get adjacent edge
419
        std::list<std::pair<PointIndex, PointIndex>>::iterator pEI;
420
        for (pEI = aclEdges.begin(); pEI != aclEdges.end(); ++pEI) {
421
            if (pEI->first == ulLast) {
422
                ulLast = pEI->second;
423
                clBorder.push_back(ulLast);
424
                aclEdges.erase(pEI);
425
                pEI = aclEdges.begin();
426
                break;
427
            }
428
            else if (pEI->second == ulFirst) {
429
                ulFirst = pEI->first;
430
                clBorder.push_front(ulFirst);
431
                aclEdges.erase(pEI);
432
                pEI = aclEdges.begin();
433
                break;
434
            }
435
            // Note: Using this might result into boundaries with wrong orientation.
436
            // But if the mesh has some facets with wrong orientation we might get
437
            // broken boundary curves.
438
            else if (pEI->second == ulLast && ignoreOrientation) {
439
                ulLast = pEI->first;
440
                clBorder.push_back(ulLast);
441
                aclEdges.erase(pEI);
442
                pEI = aclEdges.begin();
443
                break;
444
            }
445
            else if (pEI->first == ulFirst && ignoreOrientation) {
446
                ulFirst = pEI->second;
447
                clBorder.push_front(ulFirst);
448
                aclEdges.erase(pEI);
449
                pEI = aclEdges.begin();
450
                break;
451
            }
452
        }
453

454
        // Note: Calling erase on list iterators doesn't force a re-allocation and
455
        // thus doesn't invalidate the iterator itself, only the referenced object
456
        if ((pEI == aclEdges.end()) || aclEdges.empty() || (ulLast == ulFirst)) {
457
            // no further edge found or closed polyline, respectively
458
            rclBorders.emplace_back(clBorder.begin(), clBorder.end());
459
            clBorder.clear();
460

461
            if (!aclEdges.empty()) {
462
                // start new boundary
463
                ulFirst = aclEdges.begin()->first;
464
                ulLast = aclEdges.begin()->second;
465
                aclEdges.erase(aclEdges.begin());
466
                clBorder.push_back(ulFirst);
467
                clBorder.push_back(ulLast);
468
            }
469
        }
470
    }
471
}
472

473
void MeshAlgorithm::GetFacetBorder(FacetIndex uFacet, std::list<PointIndex>& rBorder) const
474
{
475
    const MeshFacetArray& rFAry = _rclMesh._aclFacetArray;
476
    std::list<std::pair<PointIndex, PointIndex>> openEdges;
477
    if (uFacet >= rFAry.size()) {
478
        return;
479
    }
480
    // add the open edge to the beginning of the list
481
    MeshFacetArray::_TConstIterator face = rFAry.begin() + uFacet;
482
    for (unsigned short i = 0; i < 3; i++) {
483
        if (face->_aulNeighbours[i] == FACET_INDEX_MAX) {
484
            openEdges.push_back(face->GetEdge(i));
485
        }
486
    }
487

488
    if (openEdges.empty()) {
489
        return;  // facet is not a border facet
490
    }
491

492
    for (MeshFacetArray::_TConstIterator it = rFAry.begin(); it != rFAry.end(); ++it) {
493
        if (it == face) {
494
            continue;
495
        }
496
        for (unsigned short i = 0; i < 3; i++) {
497
            if (it->_aulNeighbours[i] == FACET_INDEX_MAX) {
498
                openEdges.push_back(it->GetEdge(i));
499
            }
500
        }
501
    }
502

503
    SplitBoundaryFromOpenEdges(openEdges, rBorder);
504
}
505

506
void MeshAlgorithm::GetFacetsBorders(const std::vector<FacetIndex>& uFacets,
507
                                     std::list<std::vector<PointIndex>>& rBorders) const
508
{
509
    ResetFacetFlag(MeshFacet::TMP0);
510
    SetFacetsFlag(uFacets, MeshFacet::TMP0);
511
    ResetPointFlag(MeshPoint::TMP0);
512

513
    const MeshFacetArray& rFAry = _rclMesh._aclFacetArray;
514
    const MeshPointArray& rPAry = _rclMesh._aclPointArray;
515
    std::list<std::pair<PointIndex, PointIndex>> openEdges;
516

517
    // add the open edge to the beginning of the list
518
    for (auto it : uFacets) {
519
        const MeshFacet& face = rFAry[it];
520
        for (int i = 0; i < 3; i++) {
521
            if (face._aulNeighbours[i] == FACET_INDEX_MAX) {
522
                std::pair<PointIndex, PointIndex> openEdge = face.GetEdge(i);
523
                openEdges.push_back(openEdge);
524
                // mark all points of open edges of the given facets
525
                rPAry[openEdge.first].SetFlag(MeshPoint::TMP0);
526
                rPAry[openEdge.second].SetFlag(MeshPoint::TMP0);
527
            }
528
        }
529
    }
530

531
    if (openEdges.empty()) {
532
        return;  // none of the facets are border facets
533
    }
534

535
    for (const auto& it : rFAry) {
536
        if (it.IsFlag(MeshFacet::TMP0)) {
537
            continue;
538
        }
539
        for (int i = 0; i < 3; i++) {
540
            if (it._aulNeighbours[i] == FACET_INDEX_MAX) {
541
                openEdges.push_back(it.GetEdge(i));
542
            }
543
        }
544
    }
545

546
    // if the first element is not an edge of "uFacets" then give up
547
    while (!openEdges.empty()) {
548
        PointIndex first = openEdges.begin()->first;
549
        PointIndex second = openEdges.begin()->second;
550
        if (!rPAry[first].IsFlag(MeshPoint::TMP0)) {
551
            break;
552
        }
553
        if (!rPAry[second].IsFlag(MeshPoint::TMP0)) {
554
            break;
555
        }
556

557
        std::list<PointIndex> boundary;
558
        SplitBoundaryFromOpenEdges(openEdges, boundary);
559
        rBorders.emplace_back(boundary.begin(), boundary.end());
560
    }
561
}
562

563
void MeshAlgorithm::SplitBoundaryFromOpenEdges(
564
    std::list<std::pair<PointIndex, PointIndex>>& openEdges,
565
    std::list<PointIndex>& boundary) const
566
{
567
    // Start with the edge that is associated to uFacet
568
    if (openEdges.empty()) {
569
        return;
570
    }
571

572
    PointIndex ulFirst = openEdges.begin()->first;
573
    PointIndex ulLast = openEdges.begin()->second;
574

575
    openEdges.erase(openEdges.begin());
576
    boundary.push_back(ulFirst);
577
    boundary.push_back(ulLast);
578

579
    while (ulLast != ulFirst) {
580
        // find adjacent edge
581
        std::list<std::pair<PointIndex, PointIndex>>::iterator pEI;
582
        for (pEI = openEdges.begin(); pEI != openEdges.end(); ++pEI) {
583
            if (pEI->first == ulLast) {
584
                ulLast = pEI->second;
585
                boundary.push_back(ulLast);
586
                openEdges.erase(pEI);
587
                pEI = openEdges.begin();
588
                break;
589
            }
590
            else if (pEI->second == ulFirst) {
591
                ulFirst = pEI->first;
592
                boundary.push_front(ulFirst);
593
                openEdges.erase(pEI);
594
                pEI = openEdges.begin();
595
                break;
596
            }
597
        }
598

599
        // cannot close the border
600
        if (pEI == openEdges.end()) {
601
            break;
602
        }
603
    }
604
}
605

606
void MeshAlgorithm::SplitBoundaryLoops(std::list<std::vector<PointIndex>>& aBorders)
607
{
608
    // Count the number of open edges for each point
609
    std::map<PointIndex, int> openPointDegree;
610
    for (const auto& jt : _rclMesh._aclFacetArray) {
611
        for (int i = 0; i < 3; i++) {
612
            if (jt._aulNeighbours[i] == FACET_INDEX_MAX) {
613
                openPointDegree[jt._aulPoints[i]]++;
614
                openPointDegree[jt._aulPoints[(i + 1) % 3]]++;
615
            }
616
        }
617
    }
618

619
    // go through all boundaries and split them if needed
620
    std::list<std::vector<PointIndex>> aSplitBorders;
621
    for (const auto& aBorder : aBorders) {
622
        bool split = false;
623
        for (auto jt : aBorder) {
624
            // two (or more) boundaries meet in one non-manifold point
625
            if (openPointDegree[jt] > 2) {
626
                split = true;
627
                break;
628
            }
629
        }
630

631
        if (!split) {
632
            aSplitBorders.push_back(aBorder);
633
        }
634
        else {
635
            SplitBoundaryLoops(aBorder, aSplitBorders);
636
        }
637
    }
638

639
    aBorders = aSplitBorders;
640
}
641

642
void MeshAlgorithm::SplitBoundaryLoops(const std::vector<PointIndex>& rBound,
643
                                       std::list<std::vector<PointIndex>>& aBorders)
644
{
645
    std::map<PointIndex, int> aPtDegree;
646
    std::vector<PointIndex> cBound;
647
    for (PointIndex it : rBound) {
648
        int deg = (aPtDegree[it]++);
649
        if (deg > 0) {
650
            for (std::vector<PointIndex>::iterator jt = cBound.begin(); jt != cBound.end(); ++jt) {
651
                if (*jt == it) {
652
                    std::vector<PointIndex> cBoundLoop;
653
                    cBoundLoop.insert(cBoundLoop.end(), jt, cBound.end());
654
                    cBoundLoop.push_back(it);
655
                    cBound.erase(jt, cBound.end());
656
                    aBorders.push_back(cBoundLoop);
657
                    (aPtDegree[it]--);
658
                    break;
659
                }
660
            }
661
        }
662

663
        cBound.push_back(it);
664
    }
665
}
666

667
bool MeshAlgorithm::FillupHole(const std::vector<PointIndex>& boundary,
668
                               AbstractPolygonTriangulator& cTria,
669
                               MeshFacetArray& rFaces,
670
                               MeshPointArray& rPoints,
671
                               int level,
672
                               const MeshRefPointToFacets* pP2FStructure) const
673
{
674
    if (boundary.front() == boundary.back()) {
675
        // first and last vertex are identical
676
        if (boundary.size() < 4) {
677
            return false;  // something strange
678
        }
679
    }
680
    else if (boundary.size() < 3) {
681
        return false;  // something strange
682
    }
683

684
    // Get a facet as reference coordinate system
685
    MeshGeomFacet rTriangle;
686
    MeshFacet rFace;
687
    PointIndex refPoint0 = *(boundary.begin());
688
    PointIndex refPoint1 = *(boundary.begin() + 1);
689
    if (pP2FStructure) {
690
        const std::set<FacetIndex>& ring1 = (*pP2FStructure)[refPoint0];
691
        const std::set<FacetIndex>& ring2 = (*pP2FStructure)[refPoint1];
692
        std::vector<FacetIndex> f_int;
693
        std::set_intersection(ring1.begin(),
694
                              ring1.end(),
695
                              ring2.begin(),
696
                              ring2.end(),
697
                              std::back_insert_iterator<std::vector<FacetIndex>>(f_int));
698
        if (f_int.size() != 1) {
699
            return false;  // error, this must be an open edge!
700
        }
701

702
        rFace = _rclMesh._aclFacetArray[f_int.front()];
703
        rTriangle = _rclMesh.GetFacet(rFace);
704
    }
705
    else {
706
        bool ready = false;
707
        for (MeshFacetArray::_TConstIterator it = _rclMesh._aclFacetArray.begin();
708
             it != _rclMesh._aclFacetArray.end();
709
             ++it) {
710
            for (int i = 0; i < 3; i++) {
711
                if (((it->_aulPoints[i] == refPoint0) && (it->_aulPoints[(i + 1) % 3] == refPoint1))
712
                    || ((it->_aulPoints[i] == refPoint1)
713
                        && (it->_aulPoints[(i + 1) % 3] == refPoint0))) {
714
                    rFace = *it;
715
                    rTriangle = _rclMesh.GetFacet(*it);
716
                    ready = true;
717
                    break;
718
                }
719
            }
720

721
            if (ready) {
722
                break;
723
            }
724
        }
725
    }
726

727
    // add points to the polygon
728
    std::vector<Base::Vector3f> polygon;
729
    for (PointIndex jt : boundary) {
730
        polygon.push_back(_rclMesh._aclPointArray[jt]);
731
        rPoints.push_back(_rclMesh._aclPointArray[jt]);
732
    }
733

734
    // remove the last added point if it is duplicated
735
    std::vector<PointIndex> bounds = boundary;
736
    if (boundary.front() == boundary.back()) {
737
        bounds.pop_back();
738
        polygon.pop_back();
739
        rPoints.pop_back();
740
    }
741

742
    // There is no easy way to check whether the boundary is interior (a hole) or exterior before
743
    // performing the triangulation. Afterwards we can compare the normals of the created triangles
744
    // with the z-direction of our local coordinate system. If the scalar product is positive it was
745
    // a hole, otherwise not.
746
    cTria.SetPolygon(polygon);
747
    cTria.SetIndices(bounds);
748

749
    std::vector<Base::Vector3f> surf_pts = cTria.GetPolygon();
750
    if (pP2FStructure && level > 0) {
751
        std::set<PointIndex> index = pP2FStructure->NeighbourPoints(boundary, level);
752
        for (PointIndex it : index) {
753
            Base::Vector3f pt(_rclMesh._aclPointArray[it]);
754
            surf_pts.push_back(pt);
755
        }
756
    }
757

758
    if (cTria.TriangulatePolygon()) {
759
        // if we have enough points then we fit a surface through the points and project
760
        // the added points onto this surface
761
        cTria.PostProcessing(surf_pts);
762
        // get the facets and add the additional points to the array
763
        rFaces.insert(rFaces.end(), cTria.GetFacets().begin(), cTria.GetFacets().end());
764
        std::vector<Base::Vector3f> newVertices = cTria.AddedPoints();
765
        for (const auto& vertex : newVertices) {
766
            rPoints.push_back(vertex);
767
        }
768

769
        // Unfortunately, some algorithms do not care about the orientation of the polygon so we
770
        // cannot rely on the normal criterion to decide whether it's a hole or not.
771
        //
772
        std::vector<MeshFacet> faces = cTria.GetFacets();
773

774
        // Special case handling for a hole with three edges: the resulting facet might be
775
        // coincident with the reference facet
776
        if (faces.size() == 1) {
777
            MeshFacet first = faces.front();
778
            if (cTria.NeedsReindexing()) {
779
                first._aulPoints[0] = boundary[first._aulPoints[0]];
780
                first._aulPoints[1] = boundary[first._aulPoints[1]];
781
                first._aulPoints[2] = boundary[first._aulPoints[2]];
782
            }
783
            if (first.IsEqual(rFace)) {
784
                rFaces.clear();
785
                rPoints.clear();
786
                cTria.Discard();
787
                return false;
788
            }
789
        }
790

791
        // Get the new neighbour to our reference facet
792
        MeshFacet facet;
793
        unsigned short ref_side = rFace.Side(refPoint0, refPoint1);
794
        unsigned short tri_side = USHRT_MAX;
795
        if (cTria.NeedsReindexing()) {
796
            // the referenced indices of the polyline
797
            refPoint0 = 0;
798
            refPoint1 = 1;
799
        }
800
        if (ref_side < USHRT_MAX) {
801
            for (const auto& face : faces) {
802
                tri_side = face.Side(refPoint0, refPoint1);
803
                if (tri_side < USHRT_MAX) {
804
                    facet = face;
805
                    break;
806
                }
807
            }
808
        }
809

810
        // in case the reference facet has not an open edge print a log message
811
        if (ref_side == USHRT_MAX || tri_side == USHRT_MAX) {
812
            Base::Console().Log(
813
                "MeshAlgorithm::FillupHole: Expected open edge for facet <%d, %d, %d>\n",
814
                rFace._aulPoints[0],
815
                rFace._aulPoints[1],
816
                rFace._aulPoints[2]);
817
            rFaces.clear();
818
            rPoints.clear();
819
            cTria.Discard();
820
            return false;
821
        }
822

823
#if 1
824
        MeshGeomFacet triangle;
825
        triangle = cTria.GetTriangle(rPoints, facet);
826

827
        TriangulationVerifier* verifier = cTria.GetVerifier();
828
        if (!verifier) {
829
            return true;
830
        }
831

832
        // Now we have two adjacent triangles which we check for overlaps.
833
        // Therefore we build a separation plane that must separate the two diametrically opposed
834
        // points.
835
        Base::Vector3f planeNormal = rTriangle.GetNormal()
836
            % (rTriangle._aclPoints[(ref_side + 1) % 3] - rTriangle._aclPoints[ref_side]);
837
        Base::Vector3f planeBase = rTriangle._aclPoints[ref_side % 3];
838
        Base::Vector3f ref_point = rTriangle._aclPoints[(ref_side + 2) % 3];
839
        Base::Vector3f tri_point = triangle._aclPoints[(tri_side + 2) % 3];
840

841
        if (!verifier->Accept(planeNormal, planeBase, ref_point, tri_point)) {
842
            rFaces.clear();
843
            rPoints.clear();
844
            cTria.Discard();
845
            return false;
846
        }
847

848
        // we know to have filled a polygon, now check for the orientation
849
        if (verifier->MustFlip(triangle.GetNormal(), rTriangle.GetNormal())) {
850
            for (auto& rFace : rFaces) {
851
                rFace.FlipNormal();
852
            }
853
        }
854
#endif
855

856
        return true;
857
    }
858

859
    return false;
860
}
861

862
void MeshAlgorithm::SetFacetsProperty(const std::vector<FacetIndex>& raulInds,
863
                                      const std::vector<unsigned long>& raulProps) const
864
{
865
    if (raulInds.size() != raulProps.size()) {
866
        return;
867
    }
868

869
    std::vector<unsigned long>::const_iterator iP = raulProps.begin();
870
    for (std::vector<FacetIndex>::const_iterator i = raulInds.begin(); i != raulInds.end();
871
         ++i, ++iP) {
872
        _rclMesh._aclFacetArray[*i].SetProperty(*iP);
873
    }
874
}
875

876
void MeshAlgorithm::SetFacetsFlag(const std::vector<FacetIndex>& raulInds,
877
                                  MeshFacet::TFlagType tF) const
878
{
879
    for (FacetIndex it : raulInds) {
880
        _rclMesh._aclFacetArray[it].SetFlag(tF);
881
    }
882
}
883

884
void MeshAlgorithm::SetPointsFlag(const std::vector<FacetIndex>& raulInds,
885
                                  MeshPoint::TFlagType tF) const
886
{
887
    for (PointIndex it : raulInds) {
888
        _rclMesh._aclPointArray[it].SetFlag(tF);
889
    }
890
}
891

892
void MeshAlgorithm::GetFacetsFlag(std::vector<FacetIndex>& raulInds, MeshFacet::TFlagType tF) const
893
{
894
    raulInds.reserve(raulInds.size() + CountFacetFlag(tF));
895
    MeshFacetArray::_TConstIterator beg = _rclMesh._aclFacetArray.begin();
896
    MeshFacetArray::_TConstIterator end = _rclMesh._aclFacetArray.end();
897
    for (MeshFacetArray::_TConstIterator it = beg; it != end; ++it) {
898
        if (it->IsFlag(tF)) {
899
            raulInds.push_back(it - beg);
900
        }
901
    }
902
}
903

904
void MeshAlgorithm::GetPointsFlag(std::vector<PointIndex>& raulInds, MeshPoint::TFlagType tF) const
905
{
906
    raulInds.reserve(raulInds.size() + CountPointFlag(tF));
907
    MeshPointArray::_TConstIterator beg = _rclMesh._aclPointArray.begin();
908
    MeshPointArray::_TConstIterator end = _rclMesh._aclPointArray.end();
909
    for (MeshPointArray::_TConstIterator it = beg; it != end; ++it) {
910
        if (it->IsFlag(tF)) {
911
            raulInds.push_back(it - beg);
912
        }
913
    }
914
}
915

916
void MeshAlgorithm::ResetFacetsFlag(const std::vector<FacetIndex>& raulInds,
917
                                    MeshFacet::TFlagType tF) const
918
{
919
    for (FacetIndex it : raulInds) {
920
        _rclMesh._aclFacetArray[it].ResetFlag(tF);
921
    }
922
}
923

924
void MeshAlgorithm::ResetPointsFlag(const std::vector<FacetIndex>& raulInds,
925
                                    MeshPoint::TFlagType tF) const
926
{
927
    for (PointIndex it : raulInds) {
928
        _rclMesh._aclPointArray[it].ResetFlag(tF);
929
    }
930
}
931

932
void MeshAlgorithm::SetFacetFlag(MeshFacet::TFlagType tF) const
933
{
934
    _rclMesh._aclFacetArray.SetFlag(tF);
935
}
936

937
void MeshAlgorithm::SetPointFlag(MeshPoint::TFlagType tF) const
938
{
939
    _rclMesh._aclPointArray.SetFlag(tF);
940
}
941

942
void MeshAlgorithm::ResetFacetFlag(MeshFacet::TFlagType tF) const
943
{
944
    _rclMesh._aclFacetArray.ResetFlag(tF);
945
}
946

947
void MeshAlgorithm::ResetPointFlag(MeshPoint::TFlagType tF) const
948
{
949
    _rclMesh._aclPointArray.ResetFlag(tF);
950
}
951

952
unsigned long MeshAlgorithm::CountFacetFlag(MeshFacet::TFlagType tF) const
953
{
954
    MeshIsFlag<MeshFacet> flag;
955
    return std::count_if(_rclMesh._aclFacetArray.begin(),
956
                         _rclMesh._aclFacetArray.end(),
957
                         [flag, tF](const MeshFacet& f) {
958
                             return flag(f, tF);
959
                         });
960
}
961

962
unsigned long MeshAlgorithm::CountPointFlag(MeshPoint::TFlagType tF) const
963
{
964
    MeshIsFlag<MeshPoint> flag;
965
    return std::count_if(_rclMesh._aclPointArray.begin(),
966
                         _rclMesh._aclPointArray.end(),
967
                         [flag, tF](const MeshPoint& f) {
968
                             return flag(f, tF);
969
                         });
970
}
971

972
void MeshAlgorithm::GetFacetsFromToolMesh(const MeshKernel& rToolMesh,
973
                                          const Base::Vector3f& rcDir,
974
                                          std::vector<FacetIndex>& raclCutted) const
975
{
976
    MeshFacetIterator cFIt(_rclMesh);
977
    MeshFacetIterator cTIt(rToolMesh);
978

979
    BoundBox3f cBB = rToolMesh.GetBoundBox();
980

981
    Base::SequencerLauncher seq("Check facets...", _rclMesh.CountFacets());
982

983
    // check all facets
984
    Base::Vector3f tmp;
985
    for (cFIt.Init(); cFIt.More(); cFIt.Next()) {
986
        // check each point of each facet
987
        for (const auto& pnt : cFIt->_aclPoints) {
988
            // at least the point must be inside the bounding box of the tool mesh
989
            if (cBB.IsInBox(pnt)) {
990
                // should not cause performance problems since the tool mesh is usually rather
991
                // lightweight
992
                int ct = 0;
993
                for (cTIt.Init(); cTIt.More(); cTIt.Next()) {
994
                    if (cTIt->IsPointOfFace(pnt, MeshPoint::epsilon())) {
995
                        ct = 1;
996
                        break;  // the point lies on the tool mesh
997
                    }
998
                    else if (cTIt->Foraminate(pnt, rcDir, tmp)) {
999
                        // check if the intersection point lies in direction rcDir of the considered
1000
                        // point
1001
                        if ((tmp - pnt) * rcDir > 0) {
1002
                            ct++;
1003
                        }
1004
                    }
1005
                }
1006

1007
                // odd number => point is inside the tool mesh
1008
                if (ct % 2 == 1) {
1009
                    raclCutted.push_back(cFIt.Position());
1010
                    break;
1011
                }
1012
            }
1013
        }
1014

1015
        seq.next();
1016
    }
1017
}
1018

1019
void MeshAlgorithm::GetFacetsFromToolMesh(const MeshKernel& rToolMesh,
1020
                                          const Base::Vector3f& rcDir,
1021
                                          const MeshFacetGrid& rGrid,
1022
                                          std::vector<FacetIndex>& raclCutted) const
1023
{
1024
    // iterator over grid structure
1025
    MeshGridIterator clGridIter(rGrid);
1026
    BoundBox3f cBB = rToolMesh.GetBoundBox();
1027
    Base::Vector3f tmp;
1028

1029
    MeshFacetIterator cFIt(_rclMesh);
1030
    MeshFacetIterator cTIt(rToolMesh);
1031
    MeshAlgorithm cToolAlg(rToolMesh);
1032

1033
    // To speed up the algorithm we use the grid built up from the associated mesh. For each grid
1034
    // element we check whether it lies completely inside or outside the toolmesh or even intersects
1035
    // with the toolmesh. So we can reduce the number of facets with further tests dramatically.
1036
    // If the grid box is outside the toolmesh all the facets inside can be skipped. If the grid
1037
    // box is inside the toolmesh all facets are stored with no further tests because they must
1038
    // also lie inside the toolmesh. Finally, if the grid box intersects with the toolmesh we must
1039
    // also check for each whether it intersects with the toolmesh as well.
1040
    std::vector<FacetIndex> aulInds;
1041
    for (clGridIter.Init(); clGridIter.More(); clGridIter.Next()) {
1042
        int ret = cToolAlg.Surround(clGridIter.GetBoundBox(), rcDir);
1043

1044
        // the box is completely inside the toolmesh
1045
        if (ret == 1) {
1046
            // these facets can be removed without more checks
1047
            clGridIter.GetElements(raclCutted);
1048
        }
1049
        // the box intersects with toolmesh
1050
        else if (ret == 0) {
1051
            // these facets must be tested for intersections with the toolmesh
1052
            clGridIter.GetElements(aulInds);
1053
        }
1054
        // the box is outside the toolmesh but this could still mean that the triangles
1055
        // inside the grid intersect with the toolmesh
1056
        else if (ret == -1) {
1057
            // these facets must be tested for intersections with the toolmesh
1058
            clGridIter.GetElements(aulInds);
1059
        }
1060
    }
1061

1062
    // remove duplicates
1063
    std::sort(aulInds.begin(), aulInds.end());
1064
    aulInds.erase(std::unique(aulInds.begin(), aulInds.end()), aulInds.end());
1065
    std::sort(raclCutted.begin(), raclCutted.end());
1066
    raclCutted.erase(std::unique(raclCutted.begin(), raclCutted.end()), raclCutted.end());
1067

1068
    Base::SequencerLauncher seq("Check facets...", aulInds.size());
1069

1070
    // check all facets
1071
    for (FacetIndex it : aulInds) {
1072
        cFIt.Set(it);
1073

1074
        // check each point of each facet
1075
        for (auto point : cFIt->_aclPoints) {
1076
            // at least the point must be inside the bounding box of the tool mesh
1077
            if (cBB.IsInBox(point)) {
1078
                // should not cause performance problems since the tool mesh is usually rather
1079
                // lightweight
1080
                int ct = 0;
1081
                for (cTIt.Init(); cTIt.More(); cTIt.Next()) {
1082
                    if (cTIt->IsPointOfFace(point, MeshPoint::epsilon())) {
1083
                        ct = 1;
1084
                        break;  // the point lies on the tool mesh
1085
                    }
1086
                    else if (cTIt->Foraminate(point, rcDir, tmp)) {
1087
                        // check if the intersection point lies in direction rcDir of the considered
1088
                        // point
1089
                        if ((tmp - point) * rcDir > 0) {
1090
                            ct++;
1091
                        }
1092
                    }
1093
                }
1094

1095
                // odd number => point is inside the tool mesh
1096
                if (ct % 2 == 1) {
1097
                    raclCutted.push_back(cFIt.Position());
1098
                    break;
1099
                }
1100
            }
1101
        }
1102

1103
        seq.next();
1104
    }
1105

1106
    // remove duplicates
1107
    std::sort(raclCutted.begin(), raclCutted.end());
1108
    raclCutted.erase(std::unique(raclCutted.begin(), raclCutted.end()), raclCutted.end());
1109
}
1110

1111
int MeshAlgorithm::Surround(const Base::BoundBox3f& rBox, const Base::Vector3f& rcDir)
1112
{
1113
    Base::Vector3f pt1, pt2, tmp;
1114
    const BoundBox3f& cBB = _rclMesh.GetBoundBox();
1115

1116
    // at least both boxes intersect
1117
    if (cBB && rBox) {
1118
        // check for intersections with the actual mesh
1119
        Base::Vector3f cCorner[8] = {Base::Vector3f(rBox.MinX, rBox.MinY, rBox.MinZ),
1120
                                     Base::Vector3f(rBox.MaxX, rBox.MinY, rBox.MinZ),
1121
                                     Base::Vector3f(rBox.MaxX, rBox.MaxY, rBox.MinZ),
1122
                                     Base::Vector3f(rBox.MinX, rBox.MaxY, rBox.MinZ),
1123
                                     Base::Vector3f(rBox.MinX, rBox.MinY, rBox.MaxZ),
1124
                                     Base::Vector3f(rBox.MaxX, rBox.MinY, rBox.MaxZ),
1125
                                     Base::Vector3f(rBox.MaxX, rBox.MaxY, rBox.MaxZ),
1126
                                     Base::Vector3f(rBox.MinX, rBox.MaxY, rBox.MaxZ)};
1127

1128
        MeshFacetIterator cTIt(_rclMesh);
1129

1130
        // triangulation of the box
1131
        int triangles[36] = {0, 1, 2, 0, 2, 3, 0, 1, 5, 0, 5, 4, 0, 4, 7, 0, 7, 3,
1132
                             6, 7, 4, 6, 4, 5, 6, 2, 3, 6, 3, 7, 6, 1, 2, 6, 5, 1};
1133

1134
        std::vector<MeshGeomFacet> cFacet(12);
1135
        int id = 0;
1136
        for (size_t ii = 0; ii < 12; ii++) {
1137
            cFacet[ii]._aclPoints[0] = cCorner[triangles[id++]];
1138
            cFacet[ii]._aclPoints[1] = cCorner[triangles[id++]];
1139
            cFacet[ii]._aclPoints[2] = cCorner[triangles[id++]];
1140
        }
1141

1142
        // check for intersections of the box with the mesh
1143
        for (const auto& it : cFacet) {
1144
            for (cTIt.Init(); cTIt.More(); cTIt.Next()) {
1145
                int ret = cTIt->IntersectWithFacet(it, pt1, pt2);
1146

1147
                // the box intersects the mesh?
1148
                if (ret != 0) {
1149
                    return 0;  // => no more investigations required
1150
                }
1151
            }
1152
        }
1153

1154
        // Now we know that the box doesn't intersect with the mesh. This means that either the box
1155
        // is completely inside or outside the mesh. To check this we test one point of the box
1156
        // whether it is inside or outside.
1157
        int ct = 0;
1158
        for (cTIt.Init(); cTIt.More(); cTIt.Next()) {
1159
            if (cTIt->IsPointOfFace(cCorner[0], MeshPoint::epsilon())) {
1160
                ct = 1;
1161
                break;  // the point lies on the tool mesh
1162
            }
1163
            else if (cTIt->Foraminate(cCorner[0], rcDir, tmp)) {
1164
                // check if the intersection point lies in direction rcDir of the considered point
1165
                if ((tmp - cCorner[0]) * rcDir > 0) {
1166
                    ct++;
1167
                }
1168
            }
1169
        }
1170

1171
        // odd number => point (i.e. the box) is inside the mesh, even number => point is outside
1172
        // the mesh
1173
        return (ct % 2 == 1) ? 1 : -1;
1174
    }
1175

1176
    // no intersection the box is outside the mesh
1177
    return -1;
1178
}
1179

1180
void MeshAlgorithm::CheckFacets(const MeshFacetGrid& rclGrid,
1181
                                const Base::ViewProjMethod* pclProj,
1182
                                const Base::Polygon2d& rclPoly,
1183
                                bool bInner,
1184
                                std::vector<FacetIndex>& raulFacets) const
1185
{
1186
    std::vector<FacetIndex>::iterator it;
1187
    MeshFacetIterator clIter(_rclMesh, 0);
1188
    Base::Vector3f clPt2d;
1189
    Base::Vector3f clGravityOfFacet;
1190
    bool bNoPointInside {};
1191
    // Cache current view projection matrix since calls to Coin's projection are expensive
1192
    Base::ViewProjMatrix fixedProj(pclProj->getComposedProjectionMatrix());
1193
    // Precompute the polygon's bounding box
1194
    Base::BoundBox2d clPolyBBox = rclPoly.CalcBoundBox();
1195

1196
    // if true use grid on mesh to speed up search
1197
    if (bInner) {
1198
        BoundBox3f clBBox3d;
1199
        BoundBox2d clViewBBox;
1200
        std::vector<FacetIndex> aulAllElements;
1201
        // iterator for the bounding box grids
1202
        MeshGridIterator clGridIter(rclGrid);
1203
        for (clGridIter.Init(); clGridIter.More(); clGridIter.Next()) {
1204
            clBBox3d = clGridIter.GetBoundBox();
1205
            clViewBBox = clBBox3d.ProjectBox(&fixedProj);
1206
            if (clViewBBox.Intersect(clPolyBBox)) {
1207
                // collect all elements in aulAllElements
1208
                clGridIter.GetElements(aulAllElements);
1209
            }
1210
        }
1211

1212
        // remove duplicates
1213
        std::sort(aulAllElements.begin(), aulAllElements.end());
1214
        aulAllElements.erase(std::unique(aulAllElements.begin(), aulAllElements.end()),
1215
                             aulAllElements.end());
1216

1217
        Base::SequencerLauncher seq("Check facets", aulAllElements.size());
1218

1219
        for (it = aulAllElements.begin(); it != aulAllElements.end(); ++it) {
1220
            bNoPointInside = true;
1221
            clGravityOfFacet.Set(0.0f, 0.0f, 0.0f);
1222
            MeshGeomFacet rclFacet = _rclMesh.GetFacet(*it);
1223
            for (const auto& pnt : rclFacet._aclPoints) {
1224
                clPt2d = fixedProj(pnt);
1225
                clGravityOfFacet += clPt2d;
1226
                if (clPolyBBox.Contains(Base::Vector2d(clPt2d.x, clPt2d.y))
1227
                    && rclPoly.Contains(Base::Vector2d(clPt2d.x, clPt2d.y))) {
1228
                    raulFacets.push_back(*it);
1229
                    bNoPointInside = false;
1230
                    break;
1231
                }
1232
            }
1233

1234
            // if no facet point is inside the polygon then check also the gravity
1235
            if (bNoPointInside) {
1236
                clGravityOfFacet *= 1.0f / 3.0f;
1237

1238
                if (clPolyBBox.Contains(Base::Vector2d(clGravityOfFacet.x, clGravityOfFacet.y))
1239
                    && rclPoly.Contains(Base::Vector2d(clGravityOfFacet.x, clGravityOfFacet.y))) {
1240
                    raulFacets.push_back(*it);
1241
                }
1242
            }
1243

1244
            seq.next();
1245
        }
1246
    }
1247
    // When cutting triangles outside then go through all elements
1248
    else {
1249
        Base::SequencerLauncher seq("Check facets", _rclMesh.CountFacets());
1250
        for (clIter.Init(); clIter.More(); clIter.Next()) {
1251
            for (const auto& pnt : clIter->_aclPoints) {
1252
                clPt2d = fixedProj(pnt);
1253
                if ((clPolyBBox.Contains(Base::Vector2d(clPt2d.x, clPt2d.y))
1254
                     && !rclPoly.Contains(Base::Vector2d(clPt2d.x, clPt2d.y)))) {
1255
                    raulFacets.push_back(clIter.Position());
1256
                    break;
1257
                }
1258
            }
1259
            seq.next();
1260
        }
1261
    }
1262
}
1263

1264
void MeshAlgorithm::CheckFacets(const Base::ViewProjMethod* pclProj,
1265
                                const Base::Polygon2d& rclPoly,
1266
                                bool bInner,
1267
                                std::vector<FacetIndex>& raulFacets) const
1268
{
1269
    const MeshPointArray& p = _rclMesh.GetPoints();
1270
    const MeshFacetArray& f = _rclMesh.GetFacets();
1271
    Base::Vector3f pt2d;
1272
    // Use a bounding box to reduce number of call to Polygon::Contains
1273
    Base::BoundBox2d bb = rclPoly.CalcBoundBox();
1274
    // Precompute the screen projection matrix as Coin's projection function is expensive
1275
    Base::ViewProjMatrix fixedProj(pclProj->getComposedProjectionMatrix());
1276

1277
    FacetIndex index = 0;
1278
    for (MeshFacetArray::_TConstIterator it = f.begin(); it != f.end(); ++it, ++index) {
1279
        for (PointIndex ptIndex : it->_aulPoints) {
1280
            pt2d = fixedProj(p[ptIndex]);
1281

1282
            // First check whether the point is in the bounding box of the polygon
1283
            if ((bb.Contains(Base::Vector2d(pt2d.x, pt2d.y))
1284
                 && rclPoly.Contains(Base::Vector2d(pt2d.x, pt2d.y)))
1285
                ^ !bInner) {
1286
                raulFacets.push_back(index);
1287
                break;
1288
            }
1289
        }
1290
    }
1291
}
1292

1293
float MeshAlgorithm::Surface() const
1294
{
1295
    float fTotal = 0.0f;
1296
    MeshFacetIterator clFIter(_rclMesh);
1297

1298
    for (clFIter.Init(); clFIter.More(); clFIter.Next()) {
1299
        fTotal += clFIter->Area();
1300
    }
1301

1302
    return fTotal;
1303
}
1304

1305
void MeshAlgorithm::SubSampleByDist(float fDist, std::vector<Base::Vector3f>& rclPoints) const
1306
{
1307
    rclPoints.clear();
1308
    MeshFacetIterator clFIter(_rclMesh);
1309
    for (clFIter.Init(); clFIter.More(); clFIter.Next()) {
1310
        size_t k = rclPoints.size();
1311
        clFIter->SubSample(fDist, rclPoints);
1312
        if (rclPoints.size() == k) {
1313
            rclPoints.push_back(clFIter->GetGravityPoint());  // min. add middle point
1314
        }
1315
    }
1316
}
1317

1318
void MeshAlgorithm::SubSampleAllPoints(std::vector<Base::Vector3f>& rclPoints) const
1319
{
1320
    rclPoints.clear();
1321

1322
    // Add all Points
1323
    //
1324
    MeshPointIterator clPIter(_rclMesh);
1325
    for (clPIter.Init(); clPIter.More(); clPIter.Next()) {
1326
        rclPoints.push_back(*clPIter);
1327
    }
1328
}
1329

1330
void MeshAlgorithm::SubSampleByCount(unsigned long ulCtPoints,
1331
                                     std::vector<Base::Vector3f>& rclPoints) const
1332
{
1333
    float fDist = float(sqrt(Surface() / float(ulCtPoints)));
1334
    SubSampleByDist(fDist, rclPoints);
1335
}
1336

1337
void MeshAlgorithm::SearchFacetsFromPolyline(const std::vector<Base::Vector3f>& rclPolyline,
1338
                                             float fRadius,
1339
                                             const MeshFacetGrid& rclGrid,
1340
                                             std::vector<FacetIndex>& rclResultFacetsIndices) const
1341
{
1342
    rclResultFacetsIndices.clear();
1343
    if (rclPolyline.size() < 3) {
1344
        return;  // no polygon defined
1345
    }
1346

1347
    std::set<FacetIndex> aclFacets;
1348
    for (std::vector<Base::Vector3f>::const_iterator pV = rclPolyline.begin();
1349
         pV < (rclPolyline.end() - 1);
1350
         ++pV) {
1351
        const Base::Vector3f &rclP0 = *pV, &rclP1 = *(pV + 1);
1352

1353
        // BB eines Polyline-Segments
1354
        BoundBox3f clSegmBB(rclP0.x, rclP0.y, rclP0.z, rclP0.x, rclP0.y, rclP0.z);
1355
        clSegmBB.Add(rclP1);
1356
        clSegmBB.Enlarge(fRadius);  // BB um Suchradius vergroessern
1357

1358
        std::vector<FacetIndex> aclBBFacets;
1359
        unsigned long k = rclGrid.Inside(clSegmBB, aclBBFacets, false);
1360
        for (unsigned long i = 0; i < k; i++) {
1361
            if (_rclMesh.GetFacet(aclBBFacets[i]).DistanceToLineSegment(rclP0, rclP1) < fRadius) {
1362
                aclFacets.insert(aclBBFacets[i]);
1363
            }
1364
        }
1365
    }
1366

1367
    rclResultFacetsIndices.insert(rclResultFacetsIndices.begin(),
1368
                                  aclFacets.begin(),
1369
                                  aclFacets.end());
1370
}
1371

1372
void MeshAlgorithm::CutBorderFacets(std::vector<FacetIndex>& raclFacetIndices,
1373
                                    unsigned short usLevel) const
1374
{
1375
    std::vector<FacetIndex> aclToDelete;
1376

1377
    CheckBorderFacets(raclFacetIndices, aclToDelete, usLevel);
1378

1379
    // alle gefunden "Rand"-Facetsindizes" aus dem Array loeschen
1380
    std::vector<FacetIndex> aclResult;
1381
    std::set<FacetIndex> aclTmp(aclToDelete.begin(), aclToDelete.end());
1382

1383
    for (FacetIndex facetIndex : raclFacetIndices) {
1384
        if (aclTmp.find(facetIndex) == aclTmp.end()) {
1385
            aclResult.push_back(facetIndex);
1386
        }
1387
    }
1388

1389
    raclFacetIndices = aclResult;
1390
}
1391

1392
unsigned long MeshAlgorithm::CountBorderEdges() const
1393
{
1394
    unsigned long cnt = 0;
1395
    const MeshFacetArray& rclFAry = _rclMesh._aclFacetArray;
1396
    MeshFacetArray::_TConstIterator end = rclFAry.end();
1397
    for (MeshFacetArray::_TConstIterator it = rclFAry.begin(); it != end; ++it) {
1398
        for (FacetIndex facetIndex : it->_aulNeighbours) {
1399
            if (facetIndex == FACET_INDEX_MAX) {
1400
                cnt++;
1401
            }
1402
        }
1403
    }
1404

1405
    return cnt;
1406
}
1407

1408
void MeshAlgorithm::CheckBorderFacets(const std::vector<FacetIndex>& raclFacetIndices,
1409
                                      std::vector<FacetIndex>& raclResultIndices,
1410
                                      unsigned short usLevel) const
1411
{
1412
    ResetFacetFlag(MeshFacet::TMP0);
1413
    SetFacetsFlag(raclFacetIndices, MeshFacet::TMP0);
1414

1415
    const MeshFacetArray& rclFAry = _rclMesh._aclFacetArray;
1416

1417
    for (unsigned short usL = 0; usL < usLevel; usL++) {
1418
        for (FacetIndex facetIndex : raclFacetIndices) {
1419
            for (FacetIndex ulNB : rclFAry[facetIndex]._aulNeighbours) {
1420
                if (ulNB == FACET_INDEX_MAX) {
1421
                    raclResultIndices.push_back(facetIndex);
1422
                    rclFAry[facetIndex].ResetFlag(MeshFacet::TMP0);
1423
                    continue;
1424
                }
1425
                if (!rclFAry[ulNB].IsFlag(MeshFacet::TMP0)) {
1426
                    raclResultIndices.push_back(facetIndex);
1427
                    rclFAry[facetIndex].ResetFlag(MeshFacet::TMP0);
1428
                    continue;
1429
                }
1430
            }
1431
        }
1432
    }
1433
}
1434

1435
void MeshAlgorithm::GetBorderPoints(const std::vector<FacetIndex>& raclFacetIndices,
1436
                                    std::set<PointIndex>& raclResultPointsIndices) const
1437
{
1438
    ResetFacetFlag(MeshFacet::TMP0);
1439
    SetFacetsFlag(raclFacetIndices, MeshFacet::TMP0);
1440

1441
    const MeshFacetArray& rclFAry = _rclMesh._aclFacetArray;
1442

1443
    for (FacetIndex facetIndex : raclFacetIndices) {
1444
        for (int i = 0; i < 3; i++) {
1445
            const MeshFacet& rclFacet = rclFAry[facetIndex];
1446
            FacetIndex ulNB = rclFacet._aulNeighbours[i];
1447
            if (ulNB == FACET_INDEX_MAX) {
1448
                raclResultPointsIndices.insert(rclFacet._aulPoints[i]);
1449
                raclResultPointsIndices.insert(rclFacet._aulPoints[(i + 1) % 3]);
1450
                continue;
1451
            }
1452
            if (!rclFAry[ulNB].IsFlag(MeshFacet::TMP0)) {
1453
                raclResultPointsIndices.insert(rclFacet._aulPoints[i]);
1454
                raclResultPointsIndices.insert(rclFacet._aulPoints[(i + 1) % 3]);
1455
                continue;
1456
            }
1457
        }
1458
    }
1459
}
1460

1461
bool MeshAlgorithm::NearestPointFromPoint(const Base::Vector3f& rclPt,
1462
                                          FacetIndex& rclResFacetIndex,
1463
                                          Base::Vector3f& rclResPoint) const
1464
{
1465
    if (_rclMesh.CountFacets() == 0) {
1466
        return false;
1467
    }
1468

1469
    // calc each facet
1470
    float fMinDist = FLOAT_MAX;
1471
    FacetIndex ulInd = FACET_INDEX_MAX;
1472
    MeshFacetIterator pF(_rclMesh);
1473
    for (pF.Init(); pF.More(); pF.Next()) {
1474
        float fDist = pF->DistanceToPoint(rclPt);
1475
        if (fDist < fMinDist) {
1476
            fMinDist = fDist;
1477
            ulInd = pF.Position();
1478
        }
1479
    }
1480

1481
    MeshGeomFacet rclSFacet = _rclMesh.GetFacet(ulInd);
1482
    rclSFacet.DistanceToPoint(rclPt, rclResPoint);
1483
    rclResFacetIndex = ulInd;
1484

1485
    return true;
1486
}
1487

1488
bool MeshAlgorithm::NearestPointFromPoint(const Base::Vector3f& rclPt,
1489
                                          const MeshFacetGrid& rclGrid,
1490
                                          FacetIndex& rclResFacetIndex,
1491
                                          Base::Vector3f& rclResPoint) const
1492
{
1493
    FacetIndex ulInd = rclGrid.SearchNearestFromPoint(rclPt);
1494

1495
    if (ulInd == FACET_INDEX_MAX) {
1496
        return false;
1497
    }
1498

1499
    MeshGeomFacet rclSFacet = _rclMesh.GetFacet(ulInd);
1500
    rclSFacet.DistanceToPoint(rclPt, rclResPoint);
1501
    rclResFacetIndex = ulInd;
1502

1503
    return true;
1504
}
1505

1506
bool MeshAlgorithm::NearestPointFromPoint(const Base::Vector3f& rclPt,
1507
                                          const MeshFacetGrid& rclGrid,
1508
                                          float fMaxSearchArea,
1509
                                          FacetIndex& rclResFacetIndex,
1510
                                          Base::Vector3f& rclResPoint) const
1511
{
1512
    FacetIndex ulInd = rclGrid.SearchNearestFromPoint(rclPt, fMaxSearchArea);
1513

1514
    if (ulInd == FACET_INDEX_MAX) {
1515
        return false;  // no facets inside BoundingBox
1516
    }
1517

1518
    MeshGeomFacet rclSFacet = _rclMesh.GetFacet(ulInd);
1519
    rclSFacet.DistanceToPoint(rclPt, rclResPoint);
1520
    rclResFacetIndex = ulInd;
1521

1522
    return true;
1523
}
1524

1525
bool MeshAlgorithm::CutWithPlane(const Base::Vector3f& clBase,
1526
                                 const Base::Vector3f& clNormal,
1527
                                 const MeshFacetGrid& rclGrid,
1528
                                 std::list<std::vector<Base::Vector3f>>& rclResult,
1529
                                 float fMinEps,
1530
                                 bool bConnectPolygons) const
1531
{
1532
    std::vector<FacetIndex> aulFacets;
1533

1534
    // Search grid
1535
    MeshGridIterator clGridIter(rclGrid);
1536
    for (clGridIter.Init(); clGridIter.More(); clGridIter.Next()) {
1537
        // if Gridvoxel intersects the plane: pick up all facets for cutting
1538
        if (clGridIter.GetBoundBox().IsCutPlane(clBase, clNormal)) {
1539
            clGridIter.GetElements(aulFacets);
1540
        }
1541
    }
1542

1543
    // remove multiple triangles
1544
    std::sort(aulFacets.begin(), aulFacets.end());
1545
    aulFacets.erase(std::unique(aulFacets.begin(), aulFacets.end()), aulFacets.end());
1546

1547
    // intersect all facets with plane
1548
    std::list<std::pair<Base::Vector3f, Base::Vector3f>>
1549
        clTempPoly;  // Field with intersection lines (unsorted, not chained)
1550

1551
    for (FacetIndex facetIndex : aulFacets) {
1552
        Base::Vector3f clE1, clE2;
1553
        const MeshGeomFacet clF(_rclMesh.GetFacet(facetIndex));
1554

1555
        // Cut the facet and store the cutting path
1556
        if (clF.IntersectWithPlane(clBase, clNormal, clE1, clE2)) {
1557
            clTempPoly.emplace_back(clE1, clE2);
1558
        }
1559
    }
1560

1561
    if (bConnectPolygons) {
1562
        // std::list<std::pair<Base::Vector3f, Base::Vector3f> > rclTempLines;
1563
        std::list<std::pair<Base::Vector3f, Base::Vector3f>> rclResultLines(clTempPoly.begin(),
1564
                                                                            clTempPoly.end());
1565
        std::list<std::vector<Base::Vector3f>> tempList;
1566
        ConnectLines(clTempPoly, tempList, fMinEps);
1567
        ConnectPolygons(tempList, clTempPoly);
1568

1569
        for (auto& iter : clTempPoly) {
1570
            rclResultLines.push_front(iter);
1571
        }
1572

1573
        return ConnectLines(rclResultLines, rclResult, fMinEps);
1574
    }
1575

1576
    return ConnectLines(clTempPoly, rclResult, fMinEps);
1577
}
1578

1579
bool MeshAlgorithm::ConnectLines(std::list<std::pair<Base::Vector3f, Base::Vector3f>>& rclLines,
1580
                                 std::list<std::vector<Base::Vector3f>>& rclPolylines,
1581
                                 float fMinEps) const
1582
{
1583
    using TCIter = std::list<std::pair<Base::Vector3f, Base::Vector3f>>::iterator;
1584

1585
    // square search radius
1586
    // const float fMinEps = 1.0e-2f; // := 10 micrometer distance
1587
    fMinEps = fMinEps * fMinEps;
1588

1589
    // remove all lines whose distance is smaller than epsilon
1590
    std::list<TCIter> _clToDelete;
1591
    float fToDelDist = fMinEps / 10.0f;
1592
    for (TCIter pF = rclLines.begin(); pF != rclLines.end(); ++pF) {
1593
        if (Base::DistanceP2(pF->first, pF->second) < fToDelDist) {
1594
            _clToDelete.push_back(pF);
1595
        }
1596
    }
1597

1598
    for (auto& pI : _clToDelete) {
1599
        rclLines.erase(pI);
1600
    }
1601

1602
    while (!rclLines.empty()) {
1603
        TCIter pF;
1604

1605
        // new polyline
1606
        std::list<Base::Vector3f> clPoly;
1607

1608
        // add first line and delete from the list
1609
        Base::Vector3f clFront = rclLines.begin()->first;  // current start point of the polyline
1610
        Base::Vector3f clEnd = rclLines.begin()->second;   // current end point of the polyline
1611
        clPoly.push_back(clFront);
1612
        clPoly.push_back(clEnd);
1613
        rclLines.erase(rclLines.begin());
1614

1615
        // search for the next line on the begin/end of the polyline and add it
1616
        TCIter pFront, pEnd;
1617
        bool bFoundLine {};
1618
        do {
1619
            float fFrontMin = fMinEps, fEndMin = fMinEps;
1620
            bool bFrontFirst = false, bEndFirst = false;
1621

1622
            pFront = rclLines.end();
1623
            pEnd = rclLines.end();
1624
            bFoundLine = false;
1625

1626
            for (pF = rclLines.begin(); pF != rclLines.end(); ++pF) {
1627
                if (Base::DistanceP2(clFront, pF->first) < fFrontMin) {
1628
                    fFrontMin = Base::DistanceP2(clFront, pF->first);
1629
                    pFront = pF;
1630
                    bFrontFirst = true;
1631
                }
1632
                else if (Base::DistanceP2(clEnd, pF->first) < fEndMin) {
1633
                    fEndMin = Base::DistanceP2(clEnd, pF->first);
1634
                    pEnd = pF;
1635
                    bEndFirst = true;
1636
                }
1637
                else if (Base::DistanceP2(clFront, pF->second) < fFrontMin) {
1638
                    fFrontMin = Base::DistanceP2(clFront, pF->second);
1639
                    pFront = pF;
1640
                    bFrontFirst = false;
1641
                }
1642
                else if (Base::DistanceP2(clEnd, pF->second) < fEndMin) {
1643
                    fEndMin = Base::DistanceP2(clEnd, pF->second);
1644
                    pEnd = pF;
1645
                    bEndFirst = false;
1646
                }
1647
            }
1648

1649
            if (pFront != rclLines.end()) {
1650
                bFoundLine = true;
1651
                if (bFrontFirst) {
1652
                    clPoly.push_front(pFront->second);
1653
                    clFront = pFront->second;
1654
                }
1655
                else {
1656
                    clPoly.push_front(pFront->first);
1657
                    clFront = pFront->first;
1658
                }
1659

1660
                rclLines.erase(pFront);
1661
            }
1662

1663
            if (pEnd != rclLines.end()) {
1664
                bFoundLine = true;
1665
                if (bEndFirst) {
1666
                    clPoly.push_back(pEnd->second);
1667
                    clEnd = pEnd->second;
1668
                }
1669
                else {
1670
                    clPoly.push_back(pEnd->first);
1671
                    clEnd = pEnd->first;
1672
                }
1673

1674
                rclLines.erase(pEnd);
1675
            }
1676
        } while (bFoundLine);
1677

1678
        rclPolylines.emplace_back(clPoly.begin(), clPoly.end());
1679
    }
1680

1681
    // remove all polylines with too few length
1682
    using TPIter = std::list<std::vector<Base::Vector3f>>::iterator;
1683
    std::list<TPIter> _clPolyToDelete;
1684
    for (TPIter pJ = rclPolylines.begin(); pJ != rclPolylines.end(); ++pJ) {
1685
        if (pJ->size() == 2) {  // only one line segment
1686
            if (Base::DistanceP2(*pJ->begin(), *(pJ->begin() + 1)) <= fMinEps) {
1687
                _clPolyToDelete.push_back(pJ);
1688
            }
1689
        }
1690
    }
1691

1692
    for (auto& pK : _clPolyToDelete) {
1693
        rclPolylines.erase(pK);
1694
    }
1695

1696
    return true;
1697
}
1698

1699
bool MeshAlgorithm::ConnectPolygons(
1700
    std::list<std::vector<Base::Vector3f>>& clPolyList,
1701
    std::list<std::pair<Base::Vector3f, Base::Vector3f>>& rclLines) const
1702
{
1703

1704
    for (std::list<std::vector<Base::Vector3f>>::iterator OutIter = clPolyList.begin();
1705
         OutIter != clPolyList.end();
1706
         ++OutIter) {
1707
        if (OutIter->empty()) {
1708
            continue;
1709
        }
1710
        std::pair<Base::Vector3f, Base::Vector3f> currentSort;
1711
        float fDist = Base::Distance(OutIter->front(), OutIter->back());
1712
        currentSort.first = OutIter->front();
1713
        currentSort.second = OutIter->back();
1714

1715
        for (std::list<std::vector<Base::Vector3f>>::iterator InnerIter = clPolyList.begin();
1716
             InnerIter != clPolyList.end();
1717
             ++InnerIter) {
1718
            if (OutIter == InnerIter) {
1719
                continue;
1720
            }
1721

1722
            if (Base::Distance(OutIter->front(), InnerIter->front()) < fDist) {
1723
                currentSort.second = InnerIter->front();
1724
                fDist = Base::Distance(OutIter->front(), InnerIter->front());
1725
            }
1726

1727
            if (Base::Distance(OutIter->front(), InnerIter->back()) < fDist) {
1728
                currentSort.second = InnerIter->back();
1729
                fDist = Base::Distance(OutIter->front(), InnerIter->back());
1730
            }
1731
        }
1732

1733
        rclLines.push_front(currentSort);
1734
    }
1735

1736
    return true;
1737
}
1738

1739
void MeshAlgorithm::GetFacetsFromPlane(const MeshFacetGrid& rclGrid,
1740
                                       const Base::Vector3f& clNormal,
1741
                                       float d,
1742
                                       const Base::Vector3f& rclLeft,
1743
                                       const Base::Vector3f& rclRight,
1744
                                       std::vector<FacetIndex>& rclRes) const
1745
{
1746
    std::vector<FacetIndex> aulFacets;
1747

1748
    Base::Vector3f clBase = d * clNormal;
1749

1750
    Base::Vector3f clPtNormal(rclLeft - rclRight);
1751
    clPtNormal.Normalize();
1752

1753
    // search grid
1754
    MeshGridIterator clGridIter(rclGrid);
1755
    for (clGridIter.Init(); clGridIter.More(); clGridIter.Next()) {
1756
        // add facets from grid if the plane if cut the grid-voxel
1757
        if (clGridIter.GetBoundBox().IsCutPlane(clBase, clNormal)) {
1758
            clGridIter.GetElements(aulFacets);
1759
        }
1760
    }
1761

1762
    // testing facet against planes
1763
    for (FacetIndex facetIndex : aulFacets) {
1764
        MeshGeomFacet clSFacet = _rclMesh.GetFacet(facetIndex);
1765
        if (clSFacet.IntersectWithPlane(clBase, clNormal)) {
1766
            bool bInner = false;
1767
            for (int i = 0; (i < 3) && !bInner; i++) {
1768
                Base::Vector3f clPt = clSFacet._aclPoints[i];
1769
                if ((clPt.DistanceToPlane(rclLeft, clPtNormal) <= 0.0f)
1770
                    && (clPt.DistanceToPlane(rclRight, clPtNormal) >= 0.0f)) {
1771
                    bInner = true;
1772
                }
1773
            }
1774

1775
            if (bInner) {
1776
                rclRes.push_back(facetIndex);
1777
            }
1778
        }
1779
    }
1780
}
1781

1782
void MeshAlgorithm::PointsFromFacetsIndices(const std::vector<FacetIndex>& rvecIndices,
1783
                                            std::vector<Base::Vector3f>& rvecPoints) const
1784
{
1785
    const MeshFacetArray& rclFAry = _rclMesh._aclFacetArray;
1786
    const MeshPointArray& rclPAry = _rclMesh._aclPointArray;
1787

1788
    std::set<PointIndex> setPoints;
1789

1790
    for (FacetIndex facetIndex : rvecIndices) {
1791
        for (PointIndex pointIndex : rclFAry[facetIndex]._aulPoints) {
1792
            setPoints.insert(pointIndex);
1793
        }
1794
    }
1795

1796
    rvecPoints.clear();
1797
    for (PointIndex pointIndex : setPoints) {
1798
        rvecPoints.push_back(rclPAry[pointIndex]);
1799
    }
1800
}
1801

1802
bool MeshAlgorithm::Distance(const Base::Vector3f& rclPt,
1803
                             FacetIndex ulFacetIdx,
1804
                             float fMaxDistance,
1805
                             float& rfDistance) const
1806
{
1807
    const MeshFacetArray& rclFAry = _rclMesh._aclFacetArray;
1808
    const MeshPointArray& rclPAry = _rclMesh._aclPointArray;
1809
    const PointIndex* pulIdx = rclFAry[ulFacetIdx]._aulPoints;
1810

1811
    BoundBox3f clBB;
1812
    clBB.Add(rclPAry[*(pulIdx++)]);
1813
    clBB.Add(rclPAry[*(pulIdx++)]);
1814
    clBB.Add(rclPAry[*pulIdx]);
1815
    clBB.Enlarge(fMaxDistance);
1816

1817
    if (!clBB.IsInBox(rclPt)) {
1818
        return false;
1819
    }
1820

1821
    rfDistance = _rclMesh.GetFacet(ulFacetIdx).DistanceToPoint(rclPt);
1822

1823
    return rfDistance < fMaxDistance;
1824
}
1825

1826
float MeshAlgorithm::CalculateMinimumGridLength(float fLength,
1827
                                                const Base::BoundBox3f& rBBox,
1828
                                                unsigned long maxElements) const
1829
{
1830
    // Max. limit of grid elements
1831
    float fMaxGridElements = static_cast<float>(maxElements);
1832

1833
    // estimate the minimum allowed grid length
1834
    float fMinGridLen = static_cast<float>(
1835
        pow((rBBox.LengthX() * rBBox.LengthY() * rBBox.LengthZ() / fMaxGridElements), 0.3333f));
1836
    return std::max<float>(fMinGridLen, fLength);
1837
}
1838

1839
// ----------------------------------------------------
1840

1841
void MeshRefPointToFacets::Rebuild()
1842
{
1843
    _map.clear();
1844

1845
    const MeshPointArray& rPoints = _rclMesh.GetPoints();
1846
    const MeshFacetArray& rFacets = _rclMesh.GetFacets();
1847
    _map.resize(rPoints.size());
1848

1849
    MeshFacetArray::_TConstIterator pFBegin = rFacets.begin();
1850
    for (MeshFacetArray::_TConstIterator pFIter = rFacets.begin(); pFIter != rFacets.end();
1851
         ++pFIter) {
1852
        _map[pFIter->_aulPoints[0]].insert(pFIter - pFBegin);
1853
        _map[pFIter->_aulPoints[1]].insert(pFIter - pFBegin);
1854
        _map[pFIter->_aulPoints[2]].insert(pFIter - pFBegin);
1855
    }
1856
}
1857

1858
Base::Vector3f MeshRefPointToFacets::GetNormal(PointIndex pos) const
1859
{
1860
    const std::set<FacetIndex>& n = _map[pos];
1861
    Base::Vector3f normal;
1862
    MeshGeomFacet f;
1863
    for (FacetIndex it : n) {
1864
        f = _rclMesh.GetFacet(it);
1865
        normal += f.Area() * f.GetNormal();
1866
    }
1867

1868
    normal.Normalize();
1869
    return normal;
1870
}
1871

1872
std::set<PointIndex> MeshRefPointToFacets::NeighbourPoints(const std::vector<PointIndex>& pt,
1873
                                                           int level) const
1874
{
1875
    std::set<PointIndex> cp, nb, lp;
1876
    cp.insert(pt.begin(), pt.end());
1877
    lp.insert(pt.begin(), pt.end());
1878
    MeshFacetArray::_TConstIterator f_it = _rclMesh.GetFacets().begin();
1879
    for (int i = 0; i < level; i++) {
1880
        std::set<PointIndex> cur;
1881
        for (PointIndex it : lp) {
1882
            const std::set<FacetIndex>& ft = (*this)[it];
1883
            for (FacetIndex jt : ft) {
1884
                for (PointIndex index : f_it[jt]._aulPoints) {
1885
                    if (cp.find(index) == cp.end() && nb.find(index) == nb.end()) {
1886
                        nb.insert(index);
1887
                        cur.insert(index);
1888
                    }
1889
                }
1890
            }
1891
        }
1892

1893
        lp = cur;
1894
        if (lp.empty()) {
1895
            break;
1896
        }
1897
    }
1898
    return nb;
1899
}
1900

1901
std::set<PointIndex> MeshRefPointToFacets::NeighbourPoints(PointIndex pos) const
1902
{
1903
    std::set<PointIndex> p;
1904
    const std::set<FacetIndex>& vf = _map[pos];
1905
    for (FacetIndex it : vf) {
1906
        PointIndex p1 {}, p2 {}, p3 {};
1907
        _rclMesh.GetFacetPoints(it, p1, p2, p3);
1908
        if (p1 != pos) {
1909
            p.insert(p1);
1910
        }
1911
        if (p2 != pos) {
1912
            p.insert(p2);
1913
        }
1914
        if (p3 != pos) {
1915
            p.insert(p3);
1916
        }
1917
    }
1918

1919
    return p;
1920
}
1921

1922
void MeshRefPointToFacets::Neighbours(FacetIndex ulFacetInd,
1923
                                      float fMaxDist,
1924
                                      MeshCollector& collect) const
1925
{
1926
    std::set<FacetIndex> visited;
1927
    Base::Vector3f clCenter = _rclMesh.GetFacet(ulFacetInd).GetGravityPoint();
1928

1929
    const MeshFacetArray& rFacets = _rclMesh.GetFacets();
1930
    SearchNeighbours(rFacets, ulFacetInd, clCenter, fMaxDist * fMaxDist, visited, collect);
1931
}
1932

1933
void MeshRefPointToFacets::SearchNeighbours(const MeshFacetArray& rFacets,
1934
                                            FacetIndex index,
1935
                                            const Base::Vector3f& rclCenter,
1936
                                            float fMaxDist2,
1937
                                            std::set<FacetIndex>& visited,
1938
                                            MeshCollector& collect) const
1939
{
1940
    if (visited.find(index) != visited.end()) {
1941
        return;
1942
    }
1943

1944
    const MeshFacet& face = rFacets[index];
1945
    if (Base::DistanceP2(rclCenter, _rclMesh.GetFacet(face).GetGravityPoint()) > fMaxDist2) {
1946
        return;
1947
    }
1948

1949
    visited.insert(index);
1950
    collect.Append(_rclMesh, index);
1951
    for (PointIndex ptIndex : face._aulPoints) {
1952
        const std::set<FacetIndex>& f = (*this)[ptIndex];
1953

1954
        for (FacetIndex j : f) {
1955
            SearchNeighbours(rFacets, j, rclCenter, fMaxDist2, visited, collect);
1956
        }
1957
    }
1958
}
1959

1960
MeshFacetArray::_TConstIterator MeshRefPointToFacets::GetFacet(FacetIndex index) const
1961
{
1962
    return _rclMesh.GetFacets().begin() + index;
1963
}
1964

1965
const std::set<FacetIndex>& MeshRefPointToFacets::operator[](PointIndex pos) const
1966
{
1967
    return _map[pos];
1968
}
1969

1970
std::vector<FacetIndex> MeshRefPointToFacets::GetIndices(PointIndex pos1, PointIndex pos2) const
1971
{
1972
    std::vector<FacetIndex> intersection;
1973
    std::back_insert_iterator<std::vector<FacetIndex>> result(intersection);
1974
    const std::set<FacetIndex>& set1 = _map[pos1];
1975
    const std::set<FacetIndex>& set2 = _map[pos2];
1976
    std::set_intersection(set1.begin(), set1.end(), set2.begin(), set2.end(), result);
1977
    return intersection;
1978
}
1979

1980
std::vector<FacetIndex>
1981
MeshRefPointToFacets::GetIndices(PointIndex pos1, PointIndex pos2, PointIndex pos3) const
1982
{
1983
    std::vector<FacetIndex> intersection;
1984
    std::back_insert_iterator<std::vector<FacetIndex>> result(intersection);
1985
    std::vector<FacetIndex> set1 = GetIndices(pos1, pos2);
1986
    const std::set<FacetIndex>& set2 = _map[pos3];
1987
    std::set_intersection(set1.begin(), set1.end(), set2.begin(), set2.end(), result);
1988
    return intersection;
1989
}
1990

1991
void MeshRefPointToFacets::AddNeighbour(PointIndex pos, FacetIndex facet)
1992
{
1993
    _map[pos].insert(facet);
1994
}
1995

1996
void MeshRefPointToFacets::RemoveNeighbour(PointIndex pos, FacetIndex facet)
1997
{
1998
    _map[pos].erase(facet);
1999
}
2000

2001
void MeshRefPointToFacets::RemoveFacet(FacetIndex facetIndex)
2002
{
2003
    PointIndex p0 {}, p1 {}, p2 {};
2004
    _rclMesh.GetFacetPoints(facetIndex, p0, p1, p2);
2005

2006
    _map[p0].erase(facetIndex);
2007
    _map[p1].erase(facetIndex);
2008
    _map[p2].erase(facetIndex);
2009
}
2010

2011
//----------------------------------------------------------------------------
2012

2013
void MeshRefFacetToFacets::Rebuild()
2014
{
2015
    _map.clear();
2016

2017
    const MeshFacetArray& rFacets = _rclMesh.GetFacets();
2018
    _map.resize(rFacets.size());
2019

2020
    MeshRefPointToFacets vertexFace(_rclMesh);
2021
    MeshFacetArray::_TConstIterator pFBegin = rFacets.begin();
2022
    for (MeshFacetArray::_TConstIterator pFIter = pFBegin; pFIter != rFacets.end(); ++pFIter) {
2023
        for (PointIndex ptIndex : pFIter->_aulPoints) {
2024
            const std::set<FacetIndex>& faces = vertexFace[ptIndex];
2025
            for (FacetIndex face : faces) {
2026
                _map[pFIter - pFBegin].insert(face);
2027
            }
2028
        }
2029
    }
2030
}
2031

2032
const std::set<FacetIndex>& MeshRefFacetToFacets::operator[](FacetIndex pos) const
2033
{
2034
    return _map[pos];
2035
}
2036

2037
std::vector<FacetIndex> MeshRefFacetToFacets::GetIndices(FacetIndex pos1, FacetIndex pos2) const
2038
{
2039
    std::vector<FacetIndex> intersection;
2040
    std::back_insert_iterator<std::vector<FacetIndex>> result(intersection);
2041
    const std::set<FacetIndex>& set1 = _map[pos1];
2042
    const std::set<FacetIndex>& set2 = _map[pos2];
2043
    std::set_intersection(set1.begin(), set1.end(), set2.begin(), set2.end(), result);
2044
    return intersection;
2045
}
2046

2047
//----------------------------------------------------------------------------
2048

2049
void MeshRefPointToPoints::Rebuild()
2050
{
2051
    _map.clear();
2052

2053
    const MeshPointArray& rPoints = _rclMesh.GetPoints();
2054
    _map.resize(rPoints.size());
2055

2056
    const MeshFacetArray& rFacets = _rclMesh.GetFacets();
2057
    for (const auto& rFacet : rFacets) {
2058
        PointIndex ulP0 = rFacet._aulPoints[0];
2059
        PointIndex ulP1 = rFacet._aulPoints[1];
2060
        PointIndex ulP2 = rFacet._aulPoints[2];
2061

2062
        _map[ulP0].insert(ulP1);
2063
        _map[ulP0].insert(ulP2);
2064
        _map[ulP1].insert(ulP0);
2065
        _map[ulP1].insert(ulP2);
2066
        _map[ulP2].insert(ulP0);
2067
        _map[ulP2].insert(ulP1);
2068
    }
2069
}
2070

2071
Base::Vector3f MeshRefPointToPoints::GetNormal(PointIndex pos) const
2072
{
2073
    const MeshPointArray& rPoints = _rclMesh.GetPoints();
2074
    MeshCore::PlaneFit pf;
2075
    pf.AddPoint(rPoints[pos]);
2076
    MeshCore::MeshPoint center = rPoints[pos];
2077
    const std::set<PointIndex>& cv = _map[pos];
2078
    for (PointIndex cv_it : cv) {
2079
        pf.AddPoint(rPoints[cv_it]);
2080
        center += rPoints[cv_it];
2081
    }
2082

2083
    pf.Fit();
2084

2085
    Base::Vector3f normal = pf.GetNormal();
2086
    normal.Normalize();
2087
    return normal;
2088
}
2089

2090
float MeshRefPointToPoints::GetAverageEdgeLength(PointIndex index) const
2091
{
2092
    const MeshPointArray& rPoints = _rclMesh.GetPoints();
2093
    float len = 0.0f;
2094
    const std::set<PointIndex>& n = (*this)[index];
2095
    const Base::Vector3f& p = rPoints[index];
2096
    for (PointIndex it : n) {
2097
        len += Base::Distance(p, rPoints[it]);
2098
    }
2099
    return (len / n.size());
2100
}
2101

2102
const std::set<PointIndex>& MeshRefPointToPoints::operator[](PointIndex pos) const
2103
{
2104
    return _map[pos];
2105
}
2106

2107
void MeshRefPointToPoints::AddNeighbour(PointIndex pos, PointIndex facet)
2108
{
2109
    _map[pos].insert(facet);
2110
}
2111

2112
void MeshRefPointToPoints::RemoveNeighbour(PointIndex pos, PointIndex facet)
2113
{
2114
    _map[pos].erase(facet);
2115
}
2116

2117
//----------------------------------------------------------------------------
2118

2119
void MeshRefEdgeToFacets::Rebuild()
2120
{
2121
    _map.clear();
2122

2123
    const MeshFacetArray& rFacets = _rclMesh.GetFacets();
2124
    FacetIndex index = 0;
2125
    for (MeshFacetArray::_TConstIterator it = rFacets.begin(); it != rFacets.end(); ++it, ++index) {
2126
        for (int i = 0; i < 3; i++) {
2127
            MeshEdge e;
2128
            e.first = it->_aulPoints[i];
2129
            e.second = it->_aulPoints[(i + 1) % 3];
2130
            std::map<MeshEdge, MeshFacetPair, EdgeOrder>::iterator jt = _map.find(e);
2131
            if (jt == _map.end()) {
2132
                _map[e].first = index;
2133
                _map[e].second = FACET_INDEX_MAX;
2134
            }
2135
            else {
2136
                _map[e].second = index;
2137
            }
2138
        }
2139
    }
2140
}
2141

2142
const std::pair<FacetIndex, FacetIndex>& MeshRefEdgeToFacets::operator[](const MeshEdge& edge) const
2143
{
2144
    return _map.find(edge)->second;
2145
}
2146

2147
//----------------------------------------------------------------------------
2148

2149
void MeshRefNormalToPoints::Rebuild()
2150
{
2151
    _norm.clear();
2152

2153
    const MeshPointArray& rPoints = _rclMesh.GetPoints();
2154
    _norm.resize(rPoints.size());
2155

2156
    const MeshFacetArray& rFacets = _rclMesh.GetFacets();
2157
    for (const auto& rFacet : rFacets) {
2158
        const MeshPoint& p0 = rPoints[rFacet._aulPoints[0]];
2159
        const MeshPoint& p1 = rPoints[rFacet._aulPoints[1]];
2160
        const MeshPoint& p2 = rPoints[rFacet._aulPoints[2]];
2161
        float l2p01 = Base::DistanceP2(p0, p1);
2162
        float l2p12 = Base::DistanceP2(p1, p2);
2163
        float l2p20 = Base::DistanceP2(p2, p0);
2164

2165
        Base::Vector3f facenormal = _rclMesh.GetFacet(rFacet).GetNormal();
2166
        _norm[rFacet._aulPoints[0]] += facenormal * (1.0f / (l2p01 * l2p20));
2167
        _norm[rFacet._aulPoints[1]] += facenormal * (1.0f / (l2p12 * l2p01));
2168
        _norm[rFacet._aulPoints[2]] += facenormal * (1.0f / (l2p20 * l2p12));
2169
    }
2170
    for (auto& it : _norm) {
2171
        it.Normalize();
2172
    }
2173
}
2174

2175
const Base::Vector3f& MeshRefNormalToPoints::operator[](PointIndex pos) const
2176
{
2177
    return _norm[pos];
2178
}
2179