FreeCAD

TopoAlgorithm.cpp
1903 строки · 69.7 Кб
Перенос по словам
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 <boost/core/ignore_unused.hpp>
28
#include <queue>
29
#include <utility>
30
#endif
31

32
#include <Base/Console.h>
33
#include <Mod/Mesh/App/WildMagic4/Wm4MeshCurvature.h>
34

35
#include "Evaluation.h"
36
#include "Iterator.h"
37
#include "MeshKernel.h"
38
#include "TopoAlgorithm.h"
39
#include "Triangulation.h"
40

41

42
using namespace MeshCore;
43

44
MeshTopoAlgorithm::MeshTopoAlgorithm(MeshKernel& rclM)
45
    : _rclMesh(rclM)
46
{}
47

48
MeshTopoAlgorithm::~MeshTopoAlgorithm()
49
{
50
    if (_needsCleanup) {
51
        Cleanup();
52
    }
53
    EndCache();
54
}
55

56
bool MeshTopoAlgorithm::InsertVertex(FacetIndex ulFacetPos, const Base::Vector3f& rclPoint)
57
{
58
    MeshFacet& rclF = _rclMesh._aclFacetArray[ulFacetPos];
59
    MeshFacet clNewFacet1, clNewFacet2;
60

61
    // insert new point
62
    PointIndex ulPtCnt = _rclMesh._aclPointArray.size();
63
    PointIndex ulPtInd = this->GetOrAddIndex(rclPoint);
64
    FacetIndex ulSize = _rclMesh._aclFacetArray.size();
65

66
    if (ulPtInd < ulPtCnt) {
67
        return false;  // the given point is already part of the mesh => creating new facets would
68
                       // be an illegal operation
69
    }
70

71
    // adjust the facets
72
    //
73
    // first new facet
74
    clNewFacet1._aulPoints[0] = rclF._aulPoints[1];
75
    clNewFacet1._aulPoints[1] = rclF._aulPoints[2];
76
    clNewFacet1._aulPoints[2] = ulPtInd;
77
    clNewFacet1._aulNeighbours[0] = rclF._aulNeighbours[1];
78
    clNewFacet1._aulNeighbours[1] = ulSize + 1;
79
    clNewFacet1._aulNeighbours[2] = ulFacetPos;
80
    // second new facet
81
    clNewFacet2._aulPoints[0] = rclF._aulPoints[2];
82
    clNewFacet2._aulPoints[1] = rclF._aulPoints[0];
83
    clNewFacet2._aulPoints[2] = ulPtInd;
84
    clNewFacet2._aulNeighbours[0] = rclF._aulNeighbours[2];
85
    clNewFacet2._aulNeighbours[1] = ulFacetPos;
86
    clNewFacet2._aulNeighbours[2] = ulSize;
87
    // adjust the neighbour facet
88
    if (rclF._aulNeighbours[1] != FACET_INDEX_MAX) {
89
        _rclMesh._aclFacetArray[rclF._aulNeighbours[1]].ReplaceNeighbour(ulFacetPos, ulSize);
90
    }
91
    if (rclF._aulNeighbours[2] != FACET_INDEX_MAX) {
92
        _rclMesh._aclFacetArray[rclF._aulNeighbours[2]].ReplaceNeighbour(ulFacetPos, ulSize + 1);
93
    }
94
    // original facet
95
    rclF._aulPoints[2] = ulPtInd;
96
    rclF._aulNeighbours[1] = ulSize;
97
    rclF._aulNeighbours[2] = ulSize + 1;
98

99
    // insert new facets
100
    _rclMesh._aclFacetArray.push_back(clNewFacet1);
101
    _rclMesh._aclFacetArray.push_back(clNewFacet2);
102

103
    return true;
104
}
105

106
bool MeshTopoAlgorithm::SnapVertex(FacetIndex ulFacetPos, const Base::Vector3f& rP)
107
{
108
    MeshFacet& rFace = _rclMesh._aclFacetArray[ulFacetPos];
109
    if (!rFace.HasOpenEdge()) {
110
        return false;
111
    }
112
    Base::Vector3f cNo1 = _rclMesh.GetNormal(rFace);
113
    for (unsigned short i = 0; i < 3; i++) {
114
        if (rFace._aulNeighbours[i] == FACET_INDEX_MAX) {
115
            const Base::Vector3f& rPt1 = _rclMesh._aclPointArray[rFace._aulPoints[i]];
116
            const Base::Vector3f& rPt2 = _rclMesh._aclPointArray[rFace._aulPoints[(i + 1) % 3]];
117
            Base::Vector3f cNo2 = (rPt2 - rPt1) % cNo1;
118
            Base::Vector3f cNo3 = (rP - rPt1) % (rPt2 - rPt1);
119
            float fD2 = Base::DistanceP2(rPt1, rPt2);
120
            float fTV = (rP - rPt1) * (rPt2 - rPt1);
121

122
            // Point is on the edge
123
            if (cNo3.Length() < FLOAT_EPS) {
124
                return SplitOpenEdge(ulFacetPos, i, rP);
125
            }
126
            else if ((rP - rPt1) * cNo2 > 0.0f && fD2 >= fTV && fTV >= 0.0f) {
127
                MeshFacet cTria;
128
                cTria._aulPoints[0] = this->GetOrAddIndex(rP);
129
                cTria._aulPoints[1] = rFace._aulPoints[(i + 1) % 3];
130
                cTria._aulPoints[2] = rFace._aulPoints[i];
131
                cTria._aulNeighbours[1] = ulFacetPos;
132
                rFace._aulNeighbours[i] = _rclMesh.CountFacets();
133
                _rclMesh._aclFacetArray.push_back(cTria);
134
                return true;
135
            }
136
        }
137
    }
138

139
    return false;
140
}
141

142
void MeshTopoAlgorithm::OptimizeTopology(float fMaxAngle)
143
{
144
    // For each internal edge get the adjacent facets. When doing an edge swap we must update
145
    // this structure.
146
    std::map<std::pair<PointIndex, PointIndex>, std::vector<FacetIndex>> aEdge2Face;
147
    for (MeshFacetArray::_TIterator pI = _rclMesh._aclFacetArray.begin();
148
         pI != _rclMesh._aclFacetArray.end();
149
         ++pI) {
150
        for (int i = 0; i < 3; i++) {
151
            // ignore open edges
152
            if (pI->_aulNeighbours[i] != FACET_INDEX_MAX) {
153
                PointIndex ulPt0 =
154
                    std::min<PointIndex>(pI->_aulPoints[i], pI->_aulPoints[(i + 1) % 3]);
155
                PointIndex ulPt1 =
156
                    std::max<PointIndex>(pI->_aulPoints[i], pI->_aulPoints[(i + 1) % 3]);
157
                aEdge2Face[std::pair<PointIndex, PointIndex>(ulPt0, ulPt1)].push_back(
158
                    pI - _rclMesh._aclFacetArray.begin());
159
            }
160
        }
161
    }
162

163
    // fill up this list with all internal edges and perform swap edges until this list is empty
164
    std::list<std::pair<PointIndex, PointIndex>> aEdgeList;
165
    std::map<std::pair<PointIndex, PointIndex>, std::vector<FacetIndex>>::iterator pE;
166
    for (pE = aEdge2Face.begin(); pE != aEdge2Face.end(); ++pE) {
167
        if (pE->second.size() == 2) {  // make sure that we really have an internal edge
168
            aEdgeList.push_back(pE->first);
169
        }
170
    }
171

172
    // to be sure to avoid an endless loop
173
    unsigned long uMaxIter = 5 * aEdge2Face.size();
174

175
    // Perform a swap edge where needed
176
    while (!aEdgeList.empty() && uMaxIter > 0) {
177
        // get the first edge and remove it from the list
178
        std::pair<PointIndex, PointIndex> aEdge = aEdgeList.front();
179
        aEdgeList.pop_front();
180
        uMaxIter--;
181

182
        // get the adjacent facets to this edge
183
        pE = aEdge2Face.find(aEdge);
184

185
        // this edge has been removed some iterations before
186
        if (pE == aEdge2Face.end()) {
187
            continue;
188
        }
189

190
        // Is swap edge allowed and sensible?
191
        if (!ShouldSwapEdge(pE->second[0], pE->second[1], fMaxAngle)) {
192
            continue;
193
        }
194

195
        // ok, here we should perform a swap edge to minimize the maximum angle
196
        if (/*fMax12 > fMax34*/ true) {
197
            // swap the edge
198
            SwapEdge(pE->second[0], pE->second[1]);
199

200
            MeshFacet& rF1 = _rclMesh._aclFacetArray[pE->second[0]];
201
            MeshFacet& rF2 = _rclMesh._aclFacetArray[pE->second[1]];
202
            unsigned short side1 = rF1.Side(aEdge.first, aEdge.second);
203
            unsigned short side2 = rF2.Side(aEdge.first, aEdge.second);
204

205
            // adjust the edge list
206
            for (int i = 0; i < 3; i++) {
207
                std::map<std::pair<PointIndex, PointIndex>, std::vector<FacetIndex>>::iterator it;
208
                // first facet
209
                PointIndex ulPt0 =
210
                    std::min<PointIndex>(rF1._aulPoints[i], rF1._aulPoints[(i + 1) % 3]);
211
                PointIndex ulPt1 =
212
                    std::max<PointIndex>(rF1._aulPoints[i], rF1._aulPoints[(i + 1) % 3]);
213
                it = aEdge2Face.find(std::make_pair(ulPt0, ulPt1));
214
                if (it != aEdge2Face.end()) {
215
                    if (it->second[0] == pE->second[1]) {
216
                        it->second[0] = pE->second[0];
217
                    }
218
                    else if (it->second[1] == pE->second[1]) {
219
                        it->second[1] = pE->second[0];
220
                    }
221
                    aEdgeList.push_back(it->first);
222
                }
223

224
                // second facet
225
                ulPt0 = std::min<PointIndex>(rF2._aulPoints[i], rF2._aulPoints[(i + 1) % 3]);
226
                ulPt1 = std::max<PointIndex>(rF2._aulPoints[i], rF2._aulPoints[(i + 1) % 3]);
227
                it = aEdge2Face.find(std::make_pair(ulPt0, ulPt1));
228
                if (it != aEdge2Face.end()) {
229
                    if (it->second[0] == pE->second[0]) {
230
                        it->second[0] = pE->second[1];
231
                    }
232
                    else if (it->second[1] == pE->second[0]) {
233
                        it->second[1] = pE->second[1];
234
                    }
235
                    aEdgeList.push_back(it->first);
236
                }
237
            }
238

239
            // Now we must remove the edge and replace it through the new edge
240
            PointIndex ulPt0 = std::min<PointIndex>(rF1._aulPoints[(side1 + 1) % 3],
241
                                                    rF2._aulPoints[(side2 + 1) % 3]);
242
            PointIndex ulPt1 = std::max<PointIndex>(rF1._aulPoints[(side1 + 1) % 3],
243
                                                    rF2._aulPoints[(side2 + 1) % 3]);
244
            std::pair<PointIndex, PointIndex> aNewEdge = std::make_pair(ulPt0, ulPt1);
245
            aEdge2Face[aNewEdge] = pE->second;
246
            aEdge2Face.erase(pE);
247
        }
248
    }
249
}
250

251
// Cosine of the maximum angle in triangle (v1,v2,v3)
252
static float
253
cos_maxangle(const Base::Vector3f& v1, const Base::Vector3f& v2, const Base::Vector3f& v3)
254
{
255
    float a = Base::Distance(v2, v3);
256
    float b = Base::Distance(v3, v1);
257
    float c = Base::Distance(v1, v2);
258
    float A = a * (b * b + c * c - a * a);
259
    float B = b * (c * c + a * a - b * b);
260
    float C = c * (a * a + b * b - c * c);
261
    return 0.5f * std::min<float>(std::min<float>(A, B), C)
262
        / (a * b * c);  // min cosine == max angle
263
}
264

265
static float swap_benefit(const Base::Vector3f& v1,
266
                          const Base::Vector3f& v2,
267
                          const Base::Vector3f& v3,
268
                          const Base::Vector3f& v4)
269
{
270
    Base::Vector3f n124 = (v4 - v2) % (v1 - v2);
271
    Base::Vector3f n234 = (v3 - v2) % (v4 - v2);
272
    if ((n124 * n234) <= 0.0f) {
273
        return 0.0f;  // avoid normal flip
274
    }
275

276
    return std::max<float>(-cos_maxangle(v1, v2, v3), -cos_maxangle(v1, v3, v4))
277
        - std::max<float>(-cos_maxangle(v1, v2, v4), -cos_maxangle(v2, v3, v4));
278
}
279

280
float MeshTopoAlgorithm::SwapEdgeBenefit(FacetIndex f, int e) const
281
{
282
    const MeshFacetArray& faces = _rclMesh.GetFacets();
283
    const MeshPointArray& vertices = _rclMesh.GetPoints();
284

285
    FacetIndex n = faces[f]._aulNeighbours[e];
286
    if (n == FACET_INDEX_MAX) {
287
        return 0.0f;  // border edge
288
    }
289

290
    PointIndex v1 = faces[f]._aulPoints[e];
291
    PointIndex v2 = faces[f]._aulPoints[(e + 1) % 3];
292
    PointIndex v3 = faces[f]._aulPoints[(e + 2) % 3];
293
    unsigned short s = faces[n].Side(faces[f]);
294
    if (s == USHRT_MAX) {
295
        std::cerr << "MeshTopoAlgorithm::SwapEdgeBenefit: error in neighbourhood "
296
                  << "of faces " << f << " and " << n << std::endl;
297
        return 0.0f;  // topological error
298
    }
299
    PointIndex v4 = faces[n]._aulPoints[(s + 2) % 3];
300
    if (v3 == v4) {
301
        std::cerr << "MeshTopoAlgorithm::SwapEdgeBenefit: duplicate faces " << f << " and " << n
302
                  << std::endl;
303
        return 0.0f;  // duplicate faces
304
    }
305
    return swap_benefit(vertices[v2], vertices[v3], vertices[v1], vertices[v4]);
306
}
307

308
using FaceEdge = std::pair<FacetIndex, int>;  // (face, edge) pair
309
using FaceEdgePriority = std::pair<float, FaceEdge>;
310

311
void MeshTopoAlgorithm::OptimizeTopology()
312
{
313
    // Find all edges that can be swapped and insert them into a
314
    // priority queue
315
    const MeshFacetArray& faces = _rclMesh.GetFacets();
316
    FacetIndex nf = _rclMesh.CountFacets();
317
    std::priority_queue<FaceEdgePriority> todo;
318
    for (FacetIndex i = 0; i < nf; i++) {
319
        for (int j = 0; j < 3; j++) {
320
            float b = SwapEdgeBenefit(i, j);
321
            if (b > 0.0f) {
322
                todo.push(std::make_pair(b, std::make_pair(i, j)));
323
            }
324
        }
325
    }
326

327
    // Edges are sorted in decreasing order with respect to their benefit
328
    while (!todo.empty()) {
329
        FacetIndex f = todo.top().second.first;
330
        int e = todo.top().second.second;
331
        todo.pop();
332
        // Check again if the swap should still be done
333
        if (SwapEdgeBenefit(f, e) <= 0.0f) {
334
            continue;
335
        }
336
        // OK, swap the edge
337
        FacetIndex f2 = faces[f]._aulNeighbours[e];
338
        SwapEdge(f, f2);
339
        // Insert new edges into queue, if necessary
340
        for (int j = 0; j < 3; j++) {
341
            float b = SwapEdgeBenefit(f, j);
342
            if (b > 0.0f) {
343
                todo.push(std::make_pair(b, std::make_pair(f, j)));
344
            }
345
        }
346
        for (int j = 0; j < 3; j++) {
347
            float b = SwapEdgeBenefit(f2, j);
348
            if (b > 0.0f) {
349
                todo.push(std::make_pair(b, std::make_pair(f2, j)));
350
            }
351
        }
352
    }
353
}
354

355
void MeshTopoAlgorithm::DelaunayFlip(float fMaxAngle)
356
{
357
    // For each internal edge get the adjacent facets.
358
    std::set<std::pair<FacetIndex, FacetIndex>> aEdge2Face;
359
    FacetIndex index = 0;
360
    for (MeshFacetArray::_TIterator pI = _rclMesh._aclFacetArray.begin();
361
         pI != _rclMesh._aclFacetArray.end();
362
         ++pI, index++) {
363
        for (FacetIndex nbIndex : pI->_aulNeighbours) {
364
            // ignore open edges
365
            if (nbIndex != FACET_INDEX_MAX) {
366
                FacetIndex ulFt0 = std::min<FacetIndex>(index, nbIndex);
367
                FacetIndex ulFt1 = std::max<FacetIndex>(index, nbIndex);
368
                aEdge2Face.insert(std::pair<FacetIndex, FacetIndex>(ulFt0, ulFt1));
369
            }
370
        }
371
    }
372

373
    Base::Vector3f center;
374
    while (!aEdge2Face.empty()) {
375
        std::set<std::pair<FacetIndex, FacetIndex>>::iterator it = aEdge2Face.begin();
376
        std::pair<FacetIndex, FacetIndex> edge = *it;
377
        aEdge2Face.erase(it);
378
        if (ShouldSwapEdge(edge.first, edge.second, fMaxAngle)) {
379
            float radius = _rclMesh.GetFacet(edge.first).CenterOfCircumCircle(center);
380
            radius *= radius;
381
            const MeshFacet& face_1 = _rclMesh._aclFacetArray[edge.first];
382
            const MeshFacet& face_2 = _rclMesh._aclFacetArray[edge.second];
383
            unsigned short side = face_2.Side(edge.first);
384
            MeshPoint vertex = _rclMesh.GetPoint(face_2._aulPoints[(side + 1) % 3]);
385
            if (Base::DistanceP2(center, vertex) < radius) {
386
                SwapEdge(edge.first, edge.second);
387
                for (int i = 0; i < 3; i++) {
388
                    if (face_1._aulNeighbours[i] != FACET_INDEX_MAX
389
                        && face_1._aulNeighbours[i] != edge.second) {
390
                        FacetIndex ulFt0 =
391
                            std::min<FacetIndex>(edge.first, face_1._aulNeighbours[i]);
392
                        FacetIndex ulFt1 =
393
                            std::max<FacetIndex>(edge.first, face_1._aulNeighbours[i]);
394
                        aEdge2Face.insert(std::pair<FacetIndex, FacetIndex>(ulFt0, ulFt1));
395
                    }
396
                    if (face_2._aulNeighbours[i] != FACET_INDEX_MAX
397
                        && face_2._aulNeighbours[i] != edge.first) {
398
                        FacetIndex ulFt0 =
399
                            std::min<FacetIndex>(edge.second, face_2._aulNeighbours[i]);
400
                        FacetIndex ulFt1 =
401
                            std::max<FacetIndex>(edge.second, face_2._aulNeighbours[i]);
402
                        aEdge2Face.insert(std::pair<FacetIndex, FacetIndex>(ulFt0, ulFt1));
403
                    }
404
                }
405
            }
406
        }
407
    }
408
}
409

410
int MeshTopoAlgorithm::DelaunayFlip()
411
{
412
    int cnt_swap = 0;
413
    _rclMesh._aclFacetArray.ResetFlag(MeshFacet::TMP0);
414
    size_t cnt_facets = _rclMesh._aclFacetArray.size();
415
    for (size_t i = 0; i < cnt_facets; i++) {
416
        const MeshFacet& f_face = _rclMesh._aclFacetArray[i];
417
        if (f_face.IsFlag(MeshFacet::TMP0)) {
418
            continue;
419
        }
420
        for (int j = 0; j < 3; j++) {
421
            FacetIndex n = f_face._aulNeighbours[j];
422
            if (n != FACET_INDEX_MAX) {
423
                const MeshFacet& n_face = _rclMesh._aclFacetArray[n];
424
                if (n_face.IsFlag(MeshFacet::TMP0)) {
425
                    continue;
426
                }
427
                unsigned short k = n_face.Side(f_face);
428
                MeshGeomFacet f1 = _rclMesh.GetFacet(f_face);
429
                MeshGeomFacet f2 = _rclMesh.GetFacet(n_face);
430
                Base::Vector3f c1, c2, p1, p2;
431
                p1 = f1._aclPoints[(j + 2) % 3];
432
                p2 = f2._aclPoints[(k + 2) % 3];
433
                float r1 = f1.CenterOfCircumCircle(c1);
434
                r1 = r1 * r1;
435
                float r2 = f2.CenterOfCircumCircle(c2);
436
                r2 = r2 * r2;
437
                float d1 = Base::DistanceP2(c1, p2);
438
                float d2 = Base::DistanceP2(c2, p1);
439
                if (d1 < r1 || d2 < r2) {
440
                    SwapEdge(i, n);
441
                    cnt_swap++;
442
                    f_face.SetFlag(MeshFacet::TMP0);
443
                    n_face.SetFlag(MeshFacet::TMP0);
444
                }
445
            }
446
        }
447
    }
448

449
    return cnt_swap;
450
}
451

452
void MeshTopoAlgorithm::AdjustEdgesToCurvatureDirection()
453
{
454
    std::vector<Wm4::Vector3<float>> aPnts;
455
    MeshPointIterator cPIt(_rclMesh);
456
    aPnts.reserve(_rclMesh.CountPoints());
457
    for (cPIt.Init(); cPIt.More(); cPIt.Next()) {
458
        aPnts.emplace_back(cPIt->x, cPIt->y, cPIt->z);
459
    }
460

461
    // get all point connections
462
    std::vector<int> aIdx;
463
    const MeshFacetArray& raFts = _rclMesh.GetFacets();
464
    aIdx.reserve(3 * raFts.size());
465

466
    // Build map of edges to the referencing facets
467
    FacetIndex k = 0;
468
    std::map<std::pair<PointIndex, PointIndex>, std::list<FacetIndex>> aclEdgeMap;
469
    for (std::vector<MeshFacet>::const_iterator jt = raFts.begin(); jt != raFts.end(); ++jt, k++) {
470
        for (int i = 0; i < 3; i++) {
471
            PointIndex ulT0 = jt->_aulPoints[i];
472
            PointIndex ulT1 = jt->_aulPoints[(i + 1) % 3];
473
            PointIndex ulP0 = std::min<PointIndex>(ulT0, ulT1);
474
            PointIndex ulP1 = std::max<PointIndex>(ulT0, ulT1);
475
            aclEdgeMap[std::make_pair(ulP0, ulP1)].push_front(k);
476
            aIdx.push_back(static_cast<int>(jt->_aulPoints[i]));
477
        }
478
    }
479

480
    // compute vertex based curvatures
481
    Wm4::MeshCurvature<float> meshCurv(static_cast<int>(_rclMesh.CountPoints()),
482
                                       &(aPnts[0]),
483
                                       static_cast<int>(_rclMesh.CountFacets()),
484
                                       &(aIdx[0]));
485

486
    // get curvature information now
487
    const Wm4::Vector3<float>* aMaxCurvDir = meshCurv.GetMaxDirections();
488
    const Wm4::Vector3<float>* aMinCurvDir = meshCurv.GetMinDirections();
489
    const float* aMaxCurv = meshCurv.GetMaxCurvatures();
490
    const float* aMinCurv = meshCurv.GetMinCurvatures();
491

492
    raFts.ResetFlag(MeshFacet::VISIT);
493
    const MeshPointArray& raPts = _rclMesh.GetPoints();
494
    for (auto& kt : aclEdgeMap) {
495
        if (kt.second.size() == 2) {
496
            PointIndex uPt1 = kt.first.first;
497
            PointIndex uPt2 = kt.first.second;
498
            FacetIndex uFt1 = kt.second.front();
499
            FacetIndex uFt2 = kt.second.back();
500

501
            const MeshFacet& rFace1 = raFts[uFt1];
502
            const MeshFacet& rFace2 = raFts[uFt2];
503
            if (rFace1.IsFlag(MeshFacet::VISIT) || rFace2.IsFlag(MeshFacet::VISIT)) {
504
                continue;
505
            }
506

507
            PointIndex uPt3 {}, uPt4 {};
508
            unsigned short side = rFace1.Side(uPt1, uPt2);
509
            uPt3 = rFace1._aulPoints[(side + 2) % 3];
510
            side = rFace2.Side(uPt1, uPt2);
511
            uPt4 = rFace2._aulPoints[(side + 2) % 3];
512

513
            Wm4::Vector3<float> dir;
514
            float fActCurvature {};
515
            if (fabs(aMinCurv[uPt1]) > fabs(aMaxCurv[uPt1])) {
516
                fActCurvature = aMinCurv[uPt1];
517
                dir = aMaxCurvDir[uPt1];
518
            }
519
            else {
520
                fActCurvature = aMaxCurv[uPt1];
521
                dir = aMinCurvDir[uPt1];
522
            }
523

524
            Base::Vector3f cMinDir(dir.X(), dir.Y(), dir.Z());
525
            Base::Vector3f cEdgeDir1 = raPts[uPt1] - raPts[uPt2];
526
            Base::Vector3f cEdgeDir2 = raPts[uPt3] - raPts[uPt4];
527
            cMinDir.Normalize();
528
            cEdgeDir1.Normalize();
529
            cEdgeDir2.Normalize();
530

531
            // get the plane and calculate the distance to the fourth point
532
            MeshGeomFacet cPlane(raPts[uPt1], raPts[uPt2], raPts[uPt3]);
533
            // positive or negative distance
534
            float fDist = raPts[uPt4].DistanceToPlane(cPlane._aclPoints[0], cPlane.GetNormal());
535

536
            float fLength12 = Base::Distance(raPts[uPt1], raPts[uPt2]);
537
            float fLength34 = Base::Distance(raPts[uPt3], raPts[uPt4]);
538
            if (fabs(cEdgeDir1 * cMinDir) < fabs(cEdgeDir2 * cMinDir)) {
539
                if (IsSwapEdgeLegal(uFt1, uFt2) && fLength34 < 1.05f * fLength12
540
                    && fActCurvature * fDist > 0.0f) {
541
                    SwapEdge(uFt1, uFt2);
542
                    rFace1.SetFlag(MeshFacet::VISIT);
543
                    rFace2.SetFlag(MeshFacet::VISIT);
544
                }
545
            }
546
        }
547
    }
548
}
549

550
bool MeshTopoAlgorithm::InsertVertexAndSwapEdge(FacetIndex ulFacetPos,
551
                                                const Base::Vector3f& rclPoint,
552
                                                float fMaxAngle)
553
{
554
    if (!InsertVertex(ulFacetPos, rclPoint)) {
555
        return false;
556
    }
557

558
    // get the created elements
559
    FacetIndex ulF1Ind = _rclMesh._aclFacetArray.size() - 2;
560
    FacetIndex ulF2Ind = _rclMesh._aclFacetArray.size() - 1;
561
    MeshFacet& rclF1 = _rclMesh._aclFacetArray[ulFacetPos];
562
    MeshFacet& rclF2 = _rclMesh._aclFacetArray[ulF1Ind];
563
    MeshFacet& rclF3 = _rclMesh._aclFacetArray[ulF2Ind];
564

565
    // first facet
566
    for (FacetIndex uNeighbour : rclF1._aulNeighbours) {
567
        if (uNeighbour != FACET_INDEX_MAX && uNeighbour != ulF1Ind && uNeighbour != ulF2Ind) {
568
            if (ShouldSwapEdge(ulFacetPos, uNeighbour, fMaxAngle)) {
569
                SwapEdge(ulFacetPos, uNeighbour);
570
                break;
571
            }
572
        }
573
    }
574
    for (FacetIndex uNeighbour : rclF2._aulNeighbours) {
575
        // second facet
576
        if (uNeighbour != FACET_INDEX_MAX && uNeighbour != ulFacetPos && uNeighbour != ulF2Ind) {
577
            if (ShouldSwapEdge(ulF1Ind, uNeighbour, fMaxAngle)) {
578
                SwapEdge(ulF1Ind, uNeighbour);
579
                break;
580
            }
581
        }
582
    }
583

584
    // third facet
585
    for (FacetIndex uNeighbour : rclF3._aulNeighbours) {
586
        if (uNeighbour != FACET_INDEX_MAX && uNeighbour != ulFacetPos && uNeighbour != ulF1Ind) {
587
            if (ShouldSwapEdge(ulF2Ind, uNeighbour, fMaxAngle)) {
588
                SwapEdge(ulF2Ind, uNeighbour);
589
                break;
590
            }
591
        }
592
    }
593

594
    return true;
595
}
596

597
bool MeshTopoAlgorithm::IsSwapEdgeLegal(FacetIndex ulFacetPos, FacetIndex ulNeighbour) const
598
{
599
    MeshFacet& rclF = _rclMesh._aclFacetArray[ulFacetPos];
600
    MeshFacet& rclN = _rclMesh._aclFacetArray[ulNeighbour];
601

602
    unsigned short uFSide = rclF.Side(rclN);
603
    unsigned short uNSide = rclN.Side(rclF);
604

605
    if (uFSide == USHRT_MAX || uNSide == USHRT_MAX) {
606
        return false;  // not neighbours
607
    }
608

609
    Base::Vector3f cP1 = _rclMesh._aclPointArray[rclF._aulPoints[uFSide]];
610
    Base::Vector3f cP2 = _rclMesh._aclPointArray[rclF._aulPoints[(uFSide + 1) % 3]];
611
    Base::Vector3f cP3 = _rclMesh._aclPointArray[rclF._aulPoints[(uFSide + 2) % 3]];
612
    Base::Vector3f cP4 = _rclMesh._aclPointArray[rclN._aulPoints[(uNSide + 2) % 3]];
613

614
    // do not allow to create degenerated triangles
615
    MeshGeomFacet cT3(cP4, cP3, cP1);
616
    if (cT3.IsDegenerated(MeshDefinitions::_fMinPointDistanceP2)) {
617
        return false;
618
    }
619
    MeshGeomFacet cT4(cP3, cP4, cP2);
620
    if (cT4.IsDegenerated(MeshDefinitions::_fMinPointDistanceP2)) {
621
        return false;
622
    }
623

624
    // We must make sure that the two adjacent triangles builds a convex polygon, otherwise
625
    // the swap edge operation is illegal
626
    Base::Vector3f cU = cP2 - cP1;
627
    Base::Vector3f cV = cP4 - cP3;
628
    // build a helper plane through cP1 that must separate cP3 and cP4
629
    Base::Vector3f cN1 = (cU % cV) % cU;
630
    if (((cP3 - cP1) * cN1) * ((cP4 - cP1) * cN1) >= 0.0f) {
631
        return false;  // not convex
632
    }
633
    // build a helper plane through cP3 that must separate cP1 and cP2
634
    Base::Vector3f cN2 = (cU % cV) % cV;
635
    if (((cP1 - cP3) * cN2) * ((cP2 - cP3) * cN2) >= 0.0f) {
636
        return false;  // not convex
637
    }
638

639
    return true;
640
}
641

642
bool MeshTopoAlgorithm::ShouldSwapEdge(FacetIndex ulFacetPos,
643
                                       FacetIndex ulNeighbour,
644
                                       float fMaxAngle) const
645
{
646
    if (!IsSwapEdgeLegal(ulFacetPos, ulNeighbour)) {
647
        return false;
648
    }
649

650
    MeshFacet& rclF = _rclMesh._aclFacetArray[ulFacetPos];
651
    MeshFacet& rclN = _rclMesh._aclFacetArray[ulNeighbour];
652

653
    unsigned short uFSide = rclF.Side(rclN);
654
    unsigned short uNSide = rclN.Side(rclF);
655

656
    Base::Vector3f cP1 = _rclMesh._aclPointArray[rclF._aulPoints[uFSide]];
657
    Base::Vector3f cP2 = _rclMesh._aclPointArray[rclF._aulPoints[(uFSide + 1) % 3]];
658
    Base::Vector3f cP3 = _rclMesh._aclPointArray[rclF._aulPoints[(uFSide + 2) % 3]];
659
    Base::Vector3f cP4 = _rclMesh._aclPointArray[rclN._aulPoints[(uNSide + 2) % 3]];
660

661
    MeshGeomFacet cT1(cP1, cP2, cP3);
662
    float fMax1 = cT1.MaximumAngle();
663
    MeshGeomFacet cT2(cP2, cP1, cP4);
664
    float fMax2 = cT2.MaximumAngle();
665
    MeshGeomFacet cT3(cP4, cP3, cP1);
666
    float fMax3 = cT3.MaximumAngle();
667
    MeshGeomFacet cT4(cP3, cP4, cP2);
668
    float fMax4 = cT4.MaximumAngle();
669

670
    // get the angle between the triangles
671
    Base::Vector3f cN1 = cT1.GetNormal();
672
    Base::Vector3f cN2 = cT2.GetNormal();
673
    if (cN1.GetAngle(cN2) > fMaxAngle) {
674
        return false;
675
    }
676

677
    float fMax12 = std::max<float>(fMax1, fMax2);
678
    float fMax34 = std::max<float>(fMax3, fMax4);
679

680
    return fMax12 > fMax34;
681
}
682

683
void MeshTopoAlgorithm::SwapEdge(FacetIndex ulFacetPos, FacetIndex ulNeighbour)
684
{
685
    MeshFacet& rclF = _rclMesh._aclFacetArray[ulFacetPos];
686
    MeshFacet& rclN = _rclMesh._aclFacetArray[ulNeighbour];
687

688
    unsigned short uFSide = rclF.Side(rclN);
689
    unsigned short uNSide = rclN.Side(rclF);
690

691
    if (uFSide == USHRT_MAX || uNSide == USHRT_MAX) {
692
        return;  // not neighbours
693
    }
694

695
    // adjust the neighbourhood
696
    if (rclF._aulNeighbours[(uFSide + 1) % 3] != FACET_INDEX_MAX) {
697
        _rclMesh._aclFacetArray[rclF._aulNeighbours[(uFSide + 1) % 3]].ReplaceNeighbour(
698
            ulFacetPos,
699
            ulNeighbour);
700
    }
701
    if (rclN._aulNeighbours[(uNSide + 1) % 3] != FACET_INDEX_MAX) {
702
        _rclMesh._aclFacetArray[rclN._aulNeighbours[(uNSide + 1) % 3]].ReplaceNeighbour(ulNeighbour,
703
                                                                                        ulFacetPos);
704
    }
705

706
    // swap the point and neighbour indices
707
    rclF._aulPoints[(uFSide + 1) % 3] = rclN._aulPoints[(uNSide + 2) % 3];
708
    rclN._aulPoints[(uNSide + 1) % 3] = rclF._aulPoints[(uFSide + 2) % 3];
709
    rclF._aulNeighbours[uFSide] = rclN._aulNeighbours[(uNSide + 1) % 3];
710
    rclN._aulNeighbours[uNSide] = rclF._aulNeighbours[(uFSide + 1) % 3];
711
    rclF._aulNeighbours[(uFSide + 1) % 3] = ulNeighbour;
712
    rclN._aulNeighbours[(uNSide + 1) % 3] = ulFacetPos;
713
}
714

715
bool MeshTopoAlgorithm::SplitEdge(FacetIndex ulFacetPos,
716
                                  FacetIndex ulNeighbour,
717
                                  const Base::Vector3f& rP)
718
{
719
    MeshFacet& rclF = _rclMesh._aclFacetArray[ulFacetPos];
720
    MeshFacet& rclN = _rclMesh._aclFacetArray[ulNeighbour];
721

722
    unsigned short uFSide = rclF.Side(rclN);
723
    unsigned short uNSide = rclN.Side(rclF);
724

725
    if (uFSide == USHRT_MAX || uNSide == USHRT_MAX) {
726
        return false;  // not neighbours
727
    }
728

729
    PointIndex uPtCnt = _rclMesh._aclPointArray.size();
730
    PointIndex uPtInd = this->GetOrAddIndex(rP);
731
    FacetIndex ulSize = _rclMesh._aclFacetArray.size();
732

733
    // the given point is already part of the mesh => creating new facets would
734
    // be an illegal operation
735
    if (uPtInd < uPtCnt) {
736
        return false;
737
    }
738

739
    // adjust the neighbourhood
740
    if (rclF._aulNeighbours[(uFSide + 1) % 3] != FACET_INDEX_MAX) {
741
        _rclMesh._aclFacetArray[rclF._aulNeighbours[(uFSide + 1) % 3]].ReplaceNeighbour(ulFacetPos,
742
                                                                                        ulSize);
743
    }
744
    if (rclN._aulNeighbours[(uNSide + 2) % 3] != FACET_INDEX_MAX) {
745
        _rclMesh._aclFacetArray[rclN._aulNeighbours[(uNSide + 2) % 3]].ReplaceNeighbour(ulNeighbour,
746
                                                                                        ulSize + 1);
747
    }
748

749
    MeshFacet cNew1, cNew2;
750
    cNew1._aulPoints[0] = uPtInd;
751
    cNew1._aulPoints[1] = rclF._aulPoints[(uFSide + 1) % 3];
752
    cNew1._aulPoints[2] = rclF._aulPoints[(uFSide + 2) % 3];
753
    cNew1._aulNeighbours[0] = ulSize + 1;
754
    cNew1._aulNeighbours[1] = rclF._aulNeighbours[(uFSide + 1) % 3];
755
    cNew1._aulNeighbours[2] = ulFacetPos;
756

757
    cNew2._aulPoints[0] = rclN._aulPoints[uNSide];
758
    cNew2._aulPoints[1] = uPtInd;
759
    cNew2._aulPoints[2] = rclN._aulPoints[(uNSide + 2) % 3];
760
    cNew2._aulNeighbours[0] = ulSize;
761
    cNew2._aulNeighbours[1] = ulNeighbour;
762
    cNew2._aulNeighbours[2] = rclN._aulNeighbours[(uNSide + 2) % 3];
763

764
    // adjust the facets
765
    rclF._aulPoints[(uFSide + 1) % 3] = uPtInd;
766
    rclF._aulNeighbours[(uFSide + 1) % 3] = ulSize;
767
    rclN._aulPoints[uNSide] = uPtInd;
768
    rclN._aulNeighbours[(uNSide + 2) % 3] = ulSize + 1;
769

770
    // insert new facets
771
    _rclMesh._aclFacetArray.push_back(cNew1);
772
    _rclMesh._aclFacetArray.push_back(cNew2);
773

774
    return true;
775
}
776

777
bool MeshTopoAlgorithm::SplitOpenEdge(FacetIndex ulFacetPos,
778
                                      unsigned short uSide,
779
                                      const Base::Vector3f& rP)
780
{
781
    MeshFacet& rclF = _rclMesh._aclFacetArray[ulFacetPos];
782
    if (rclF._aulNeighbours[uSide] != FACET_INDEX_MAX) {
783
        return false;  // not open
784
    }
785

786
    PointIndex uPtCnt = _rclMesh._aclPointArray.size();
787
    PointIndex uPtInd = this->GetOrAddIndex(rP);
788
    FacetIndex ulSize = _rclMesh._aclFacetArray.size();
789

790
    if (uPtInd < uPtCnt) {
791
        return false;  // the given point is already part of the mesh => creating new facets would
792
                       // be an illegal operation
793
    }
794

795
    // adjust the neighbourhood
796
    if (rclF._aulNeighbours[(uSide + 1) % 3] != FACET_INDEX_MAX) {
797
        _rclMesh._aclFacetArray[rclF._aulNeighbours[(uSide + 1) % 3]].ReplaceNeighbour(ulFacetPos,
798
                                                                                       ulSize);
799
    }
800

801
    MeshFacet cNew;
802
    cNew._aulPoints[0] = uPtInd;
803
    cNew._aulPoints[1] = rclF._aulPoints[(uSide + 1) % 3];
804
    cNew._aulPoints[2] = rclF._aulPoints[(uSide + 2) % 3];
805
    cNew._aulNeighbours[0] = FACET_INDEX_MAX;
806
    cNew._aulNeighbours[1] = rclF._aulNeighbours[(uSide + 1) % 3];
807
    cNew._aulNeighbours[2] = ulFacetPos;
808

809
    // adjust the facets
810
    rclF._aulPoints[(uSide + 1) % 3] = uPtInd;
811
    rclF._aulNeighbours[(uSide + 1) % 3] = ulSize;
812

813
    // insert new facets
814
    _rclMesh._aclFacetArray.push_back(cNew);
815
    return true;
816
}
817

818
bool MeshTopoAlgorithm::Vertex_Less::operator()(const Base::Vector3f& u,
819
                                                const Base::Vector3f& v) const
820
{
821
    if (fabs(u.x - v.x) > FLOAT_EPS) {
822
        return u.x < v.x;
823
    }
824
    if (fabs(u.y - v.y) > FLOAT_EPS) {
825
        return u.y < v.y;
826
    }
827
    if (fabs(u.z - v.z) > FLOAT_EPS) {
828
        return u.z < v.z;
829
    }
830
    return false;
831
}
832

833
void MeshTopoAlgorithm::BeginCache()
834
{
835
    if (_cache) {
836
        delete _cache;
837
    }
838
    _cache = new tCache();
839
    PointIndex nbPoints = _rclMesh._aclPointArray.size();
840
    for (unsigned int pntCpt = 0; pntCpt < nbPoints; ++pntCpt) {
841
        _cache->insert(std::make_pair(_rclMesh._aclPointArray[pntCpt], pntCpt));
842
    }
843
}
844

845
void MeshTopoAlgorithm::EndCache()
846
{
847
    if (_cache) {
848
        _cache->clear();
849
        delete _cache;
850
        _cache = nullptr;
851
    }
852
}
853

854
PointIndex MeshTopoAlgorithm::GetOrAddIndex(const MeshPoint& rclPoint)
855
{
856
    if (!_cache) {
857
        return _rclMesh._aclPointArray.GetOrAddIndex(rclPoint);
858
    }
859

860
    unsigned long sz = _rclMesh._aclPointArray.size();
861
    std::pair<tCache::iterator, bool> retval = _cache->insert(std::make_pair(rclPoint, sz));
862
    if (retval.second) {
863
        _rclMesh._aclPointArray.push_back(rclPoint);
864
    }
865
    return retval.first->second;
866
}
867

868
std::vector<FacetIndex> MeshTopoAlgorithm::GetFacetsToPoint(FacetIndex uFacetPos,
869
                                                            PointIndex uPointPos) const
870
{
871
    // get all facets this point is referenced by
872
    std::list<FacetIndex> aReference;
873
    aReference.push_back(uFacetPos);
874
    std::set<FacetIndex> aRefFacet;
875
    while (!aReference.empty()) {
876
        FacetIndex uIndex = aReference.front();
877
        aReference.pop_front();
878
        aRefFacet.insert(uIndex);
879
        MeshFacet& rFace = _rclMesh._aclFacetArray[uIndex];
880
        for (int i = 0; i < 3; i++) {
881
            if (rFace._aulPoints[i] == uPointPos) {
882
                if (rFace._aulNeighbours[i] != FACET_INDEX_MAX) {
883
                    if (aRefFacet.find(rFace._aulNeighbours[i]) == aRefFacet.end()) {
884
                        aReference.push_back(rFace._aulNeighbours[i]);
885
                    }
886
                }
887
                if (rFace._aulNeighbours[(i + 2) % 3] != FACET_INDEX_MAX) {
888
                    if (aRefFacet.find(rFace._aulNeighbours[(i + 2) % 3]) == aRefFacet.end()) {
889
                        aReference.push_back(rFace._aulNeighbours[(i + 2) % 3]);
890
                    }
891
                }
892
                break;
893
            }
894
        }
895
    }
896

897
    // copy the items
898
    std::vector<FacetIndex> aRefs;
899
    aRefs.insert(aRefs.end(), aRefFacet.begin(), aRefFacet.end());
900
    return aRefs;
901
}
902

903
void MeshTopoAlgorithm::Cleanup()
904
{
905
    _rclMesh.RemoveInvalids();
906
    _needsCleanup = false;
907
}
908

909
bool MeshTopoAlgorithm::CollapseVertex(const VertexCollapse& vc)
910
{
911
    if (vc._circumFacets.size() != vc._circumPoints.size()) {
912
        return false;
913
    }
914

915
    if (vc._circumFacets.size() != 3) {
916
        return false;
917
    }
918

919
    if (!_rclMesh._aclPointArray[vc._point].IsValid()) {
920
        return false;  // the point is marked invalid from a previous run
921
    }
922

923
    MeshFacet& rFace1 = _rclMesh._aclFacetArray[vc._circumFacets[0]];
924
    MeshFacet& rFace2 = _rclMesh._aclFacetArray[vc._circumFacets[1]];
925
    MeshFacet& rFace3 = _rclMesh._aclFacetArray[vc._circumFacets[2]];
926

927
    // get the point that is not shared by rFace1
928
    PointIndex ptIndex = POINT_INDEX_MAX;
929
    std::vector<PointIndex>::const_iterator it;
930
    for (it = vc._circumPoints.begin(); it != vc._circumPoints.end(); ++it) {
931
        if (!rFace1.HasPoint(*it)) {
932
            ptIndex = *it;
933
            break;
934
        }
935
    }
936

937
    if (ptIndex == POINT_INDEX_MAX) {
938
        return false;
939
    }
940

941
    FacetIndex neighbour1 = FACET_INDEX_MAX;
942
    FacetIndex neighbour2 = FACET_INDEX_MAX;
943

944
    const std::vector<FacetIndex>& faces = vc._circumFacets;
945
    // get neighbours that are not part of the faces to be removed
946
    for (int i = 0; i < 3; i++) {
947
        if (std::find(faces.begin(), faces.end(), rFace2._aulNeighbours[i]) == faces.end()) {
948
            neighbour1 = rFace2._aulNeighbours[i];
949
        }
950
        if (std::find(faces.begin(), faces.end(), rFace3._aulNeighbours[i]) == faces.end()) {
951
            neighbour2 = rFace3._aulNeighbours[i];
952
        }
953
    }
954

955
    // adjust point and neighbour indices
956
    rFace1.Transpose(vc._point, ptIndex);
957
    rFace1.ReplaceNeighbour(vc._circumFacets[1], neighbour1);
958
    rFace1.ReplaceNeighbour(vc._circumFacets[2], neighbour2);
959

960
    if (neighbour1 != FACET_INDEX_MAX) {
961
        MeshFacet& rFace4 = _rclMesh._aclFacetArray[neighbour1];
962
        rFace4.ReplaceNeighbour(vc._circumFacets[1], vc._circumFacets[0]);
963
    }
964
    if (neighbour2 != FACET_INDEX_MAX) {
965
        MeshFacet& rFace5 = _rclMesh._aclFacetArray[neighbour2];
966
        rFace5.ReplaceNeighbour(vc._circumFacets[2], vc._circumFacets[0]);
967
    }
968

969
    // the two facets and the point can be marked for removal
970
    rFace2.SetInvalid();
971
    rFace3.SetInvalid();
972
    _rclMesh._aclPointArray[vc._point].SetInvalid();
973

974
    _needsCleanup = true;
975

976
    return true;
977
}
978

979
bool MeshTopoAlgorithm::CollapseEdge(FacetIndex ulFacetPos, FacetIndex ulNeighbour)
980
{
981
    MeshFacet& rclF = _rclMesh._aclFacetArray[ulFacetPos];
982
    MeshFacet& rclN = _rclMesh._aclFacetArray[ulNeighbour];
983

984
    unsigned short uFSide = rclF.Side(rclN);
985
    unsigned short uNSide = rclN.Side(rclF);
986

987
    if (uFSide == USHRT_MAX || uNSide == USHRT_MAX) {
988
        return false;  // not neighbours
989
    }
990

991
    if (!rclF.IsValid() || !rclN.IsValid()) {
992
        return false;  // the facets are marked invalid from a previous run
993
    }
994

995
    // get the point index we want to remove
996
    PointIndex ulPointPos = rclF._aulPoints[uFSide];
997
    PointIndex ulPointNew = rclN._aulPoints[uNSide];
998

999
    // get all facets this point is referenced by
1000
    std::vector<FacetIndex> aRefs = GetFacetsToPoint(ulFacetPos, ulPointPos);
1001
    for (FacetIndex it : aRefs) {
1002
        MeshFacet& rFace = _rclMesh._aclFacetArray[it];
1003
        rFace.Transpose(ulPointPos, ulPointNew);
1004
    }
1005

1006
    // set the new neighbourhood
1007
    if (rclF._aulNeighbours[(uFSide + 1) % 3] != FACET_INDEX_MAX) {
1008
        _rclMesh._aclFacetArray[rclF._aulNeighbours[(uFSide + 1) % 3]].ReplaceNeighbour(
1009
            ulFacetPos,
1010
            rclF._aulNeighbours[(uFSide + 2) % 3]);
1011
    }
1012
    if (rclF._aulNeighbours[(uFSide + 2) % 3] != FACET_INDEX_MAX) {
1013
        _rclMesh._aclFacetArray[rclF._aulNeighbours[(uFSide + 2) % 3]].ReplaceNeighbour(
1014
            ulFacetPos,
1015
            rclF._aulNeighbours[(uFSide + 1) % 3]);
1016
    }
1017
    if (rclN._aulNeighbours[(uNSide + 1) % 3] != FACET_INDEX_MAX) {
1018
        _rclMesh._aclFacetArray[rclN._aulNeighbours[(uNSide + 1) % 3]].ReplaceNeighbour(
1019
            ulNeighbour,
1020
            rclN._aulNeighbours[(uNSide + 2) % 3]);
1021
    }
1022
    if (rclN._aulNeighbours[(uNSide + 2) % 3] != FACET_INDEX_MAX) {
1023
        _rclMesh._aclFacetArray[rclN._aulNeighbours[(uNSide + 2) % 3]].ReplaceNeighbour(
1024
            ulNeighbour,
1025
            rclN._aulNeighbours[(uNSide + 1) % 3]);
1026
    }
1027

1028
    // isolate the both facets and the point
1029
    rclF._aulNeighbours[0] = FACET_INDEX_MAX;
1030
    rclF._aulNeighbours[1] = FACET_INDEX_MAX;
1031
    rclF._aulNeighbours[2] = FACET_INDEX_MAX;
1032
    rclF.SetInvalid();
1033
    rclN._aulNeighbours[0] = FACET_INDEX_MAX;
1034
    rclN._aulNeighbours[1] = FACET_INDEX_MAX;
1035
    rclN._aulNeighbours[2] = FACET_INDEX_MAX;
1036
    rclN.SetInvalid();
1037
    _rclMesh._aclPointArray[ulPointPos].SetInvalid();
1038

1039
    _needsCleanup = true;
1040

1041
    return true;
1042
}
1043

1044
bool MeshTopoAlgorithm::IsCollapseEdgeLegal(const EdgeCollapse& ec) const
1045
{
1046
    // http://stackoverflow.com/a/27049418/148668
1047
    // Check connectivity
1048
    //
1049
    std::vector<PointIndex> commonPoints;
1050
    std::set_intersection(ec._adjacentFrom.begin(),
1051
                          ec._adjacentFrom.end(),
1052
                          ec._adjacentTo.begin(),
1053
                          ec._adjacentTo.end(),
1054
                          std::back_insert_iterator<std::vector<PointIndex>>(commonPoints));
1055
    if (commonPoints.size() > 2) {
1056
        return false;
1057
    }
1058

1059
    // Check geometry
1060
    std::vector<FacetIndex>::const_iterator it;
1061
    for (it = ec._changeFacets.begin(); it != ec._changeFacets.end(); ++it) {
1062
        MeshFacet f = _rclMesh._aclFacetArray[*it];
1063
        if (!f.IsValid()) {
1064
            return false;
1065
        }
1066

1067
        // ignore the facet(s) at this edge
1068
        if (f.HasPoint(ec._fromPoint) && f.HasPoint(ec._toPoint)) {
1069
            continue;
1070
        }
1071

1072
        MeshGeomFacet tria1 = _rclMesh.GetFacet(f);
1073
        f.Transpose(ec._fromPoint, ec._toPoint);
1074
        MeshGeomFacet tria2 = _rclMesh.GetFacet(f);
1075

1076
        if (tria1.GetNormal() * tria2.GetNormal() < 0.0f) {
1077
            return false;
1078
        }
1079
    }
1080

1081
    // If the data structure is valid and the algorithm works as expected
1082
    // it should never happen to reject the edge-collapse here!
1083
    for (it = ec._removeFacets.begin(); it != ec._removeFacets.end(); ++it) {
1084
        MeshFacet f = _rclMesh._aclFacetArray[*it];
1085
        if (!f.IsValid()) {
1086
            return false;
1087
        }
1088
    }
1089

1090
    if (!_rclMesh._aclPointArray[ec._fromPoint].IsValid()) {
1091
        return false;
1092
    }
1093

1094
    if (!_rclMesh._aclPointArray[ec._toPoint].IsValid()) {
1095
        return false;
1096
    }
1097

1098
    return true;
1099
}
1100

1101
bool MeshTopoAlgorithm::CollapseEdge(const EdgeCollapse& ec)
1102
{
1103
    std::vector<FacetIndex>::const_iterator it;
1104
    for (it = ec._removeFacets.begin(); it != ec._removeFacets.end(); ++it) {
1105
        MeshFacet& f = _rclMesh._aclFacetArray[*it];
1106
        f.SetInvalid();
1107

1108
        // adjust the neighbourhood
1109
        std::vector<FacetIndex> neighbours;
1110
        for (FacetIndex nbIndex : f._aulNeighbours) {
1111
            // get the neighbours of the facet that won't be invalidated
1112
            if (nbIndex != FACET_INDEX_MAX) {
1113
                if (std::find(ec._removeFacets.begin(), ec._removeFacets.end(), nbIndex)
1114
                    == ec._removeFacets.end()) {
1115
                    neighbours.push_back(nbIndex);
1116
                }
1117
            }
1118
        }
1119

1120
        if (neighbours.size() == 2) {
1121
            MeshFacet& n1 = _rclMesh._aclFacetArray[neighbours[0]];
1122
            n1.ReplaceNeighbour(*it, neighbours[1]);
1123
            MeshFacet& n2 = _rclMesh._aclFacetArray[neighbours[1]];
1124
            n2.ReplaceNeighbour(*it, neighbours[0]);
1125
        }
1126
        else if (neighbours.size() == 1) {
1127
            MeshFacet& n1 = _rclMesh._aclFacetArray[neighbours[0]];
1128
            n1.ReplaceNeighbour(*it, FACET_INDEX_MAX);
1129
        }
1130
    }
1131

1132
    for (it = ec._changeFacets.begin(); it != ec._changeFacets.end(); ++it) {
1133
        MeshFacet& f = _rclMesh._aclFacetArray[*it];
1134
        f.Transpose(ec._fromPoint, ec._toPoint);
1135
    }
1136

1137
    _rclMesh._aclPointArray[ec._fromPoint].SetInvalid();
1138

1139
    _needsCleanup = true;
1140
    return true;
1141
}
1142

1143
bool MeshTopoAlgorithm::CollapseFacet(FacetIndex ulFacetPos)
1144
{
1145
    MeshFacet& rclF = _rclMesh._aclFacetArray[ulFacetPos];
1146
    if (!rclF.IsValid()) {
1147
        return false;  // the facet is marked invalid from a previous run
1148
    }
1149

1150
    // get the point index we want to remove
1151
    PointIndex ulPointInd0 = rclF._aulPoints[0];
1152
    PointIndex ulPointInd1 = rclF._aulPoints[1];
1153
    PointIndex ulPointInd2 = rclF._aulPoints[2];
1154

1155
    // move the vertex to the gravity center
1156
    Base::Vector3f cCenter = _rclMesh.GetGravityPoint(rclF);
1157
    _rclMesh._aclPointArray[ulPointInd0] = cCenter;
1158

1159
    // set the new point indices for all facets that share one of the points to be deleted
1160
    std::vector<FacetIndex> aRefs = GetFacetsToPoint(ulFacetPos, ulPointInd1);
1161
    for (FacetIndex it : aRefs) {
1162
        MeshFacet& rFace = _rclMesh._aclFacetArray[it];
1163
        rFace.Transpose(ulPointInd1, ulPointInd0);
1164
    }
1165

1166
    aRefs = GetFacetsToPoint(ulFacetPos, ulPointInd2);
1167
    for (FacetIndex it : aRefs) {
1168
        MeshFacet& rFace = _rclMesh._aclFacetArray[it];
1169
        rFace.Transpose(ulPointInd2, ulPointInd0);
1170
    }
1171

1172
    // set the neighbourhood of the circumjacent facets
1173
    for (FacetIndex nbIndex : rclF._aulNeighbours) {
1174
        if (nbIndex == FACET_INDEX_MAX) {
1175
            continue;
1176
        }
1177
        MeshFacet& rclN = _rclMesh._aclFacetArray[nbIndex];
1178
        unsigned short uNSide = rclN.Side(rclF);
1179

1180
        if (rclN._aulNeighbours[(uNSide + 1) % 3] != FACET_INDEX_MAX) {
1181
            _rclMesh._aclFacetArray[rclN._aulNeighbours[(uNSide + 1) % 3]].ReplaceNeighbour(
1182
                nbIndex,
1183
                rclN._aulNeighbours[(uNSide + 2) % 3]);
1184
        }
1185
        if (rclN._aulNeighbours[(uNSide + 2) % 3] != FACET_INDEX_MAX) {
1186
            _rclMesh._aclFacetArray[rclN._aulNeighbours[(uNSide + 2) % 3]].ReplaceNeighbour(
1187
                nbIndex,
1188
                rclN._aulNeighbours[(uNSide + 1) % 3]);
1189
        }
1190

1191
        // Isolate the neighbours from the topology
1192
        rclN._aulNeighbours[0] = FACET_INDEX_MAX;
1193
        rclN._aulNeighbours[1] = FACET_INDEX_MAX;
1194
        rclN._aulNeighbours[2] = FACET_INDEX_MAX;
1195
        rclN.SetInvalid();
1196
    }
1197

1198
    // Isolate this facet and make two of its points invalid
1199
    rclF._aulNeighbours[0] = FACET_INDEX_MAX;
1200
    rclF._aulNeighbours[1] = FACET_INDEX_MAX;
1201
    rclF._aulNeighbours[2] = FACET_INDEX_MAX;
1202
    rclF.SetInvalid();
1203
    _rclMesh._aclPointArray[ulPointInd1].SetInvalid();
1204
    _rclMesh._aclPointArray[ulPointInd2].SetInvalid();
1205

1206
    _needsCleanup = true;
1207

1208
    return true;
1209
}
1210

1211
void MeshTopoAlgorithm::SplitFacet(FacetIndex ulFacetPos,
1212
                                   const Base::Vector3f& rP1,
1213
                                   const Base::Vector3f& rP2)
1214
{
1215
    float fEps = MESH_MIN_EDGE_LEN;
1216
    MeshFacet& rFace = _rclMesh._aclFacetArray[ulFacetPos];
1217
    MeshPoint& rVertex0 = _rclMesh._aclPointArray[rFace._aulPoints[0]];
1218
    MeshPoint& rVertex1 = _rclMesh._aclPointArray[rFace._aulPoints[1]];
1219
    MeshPoint& rVertex2 = _rclMesh._aclPointArray[rFace._aulPoints[2]];
1220

1221
    auto pointIndex = [=](const Base::Vector3f& rP) {
1222
        unsigned short equalP = USHRT_MAX;
1223
        if (Base::Distance(rVertex0, rP) < fEps) {
1224
            equalP = 0;
1225
        }
1226
        else if (Base::Distance(rVertex1, rP) < fEps) {
1227
            equalP = 1;
1228
        }
1229
        else if (Base::Distance(rVertex2, rP) < fEps) {
1230
            equalP = 2;
1231
        }
1232
        return equalP;
1233
    };
1234

1235
    unsigned short equalP1 = pointIndex(rP1);
1236
    unsigned short equalP2 = pointIndex(rP2);
1237

1238
    // both points are coincident with the corner points
1239
    if (equalP1 != USHRT_MAX && equalP2 != USHRT_MAX) {
1240
        return;  // must not split the facet
1241
    }
1242

1243
    if (equalP1 != USHRT_MAX) {
1244
        // get the edge to the second given point and perform a split edge operation
1245
        SplitFacetOnOneEdge(ulFacetPos, rP2);
1246
    }
1247
    else if (equalP2 != USHRT_MAX) {
1248
        // get the edge to the first given point and perform a split edge operation
1249
        SplitFacetOnOneEdge(ulFacetPos, rP1);
1250
    }
1251
    else {
1252
        SplitFacetOnTwoEdges(ulFacetPos, rP1, rP2);
1253
    }
1254
}
1255

1256
void MeshTopoAlgorithm::SplitFacetOnOneEdge(FacetIndex ulFacetPos, const Base::Vector3f& rP)
1257
{
1258
    float fMinDist = FLOAT_MAX;
1259
    unsigned short iEdgeNo = USHRT_MAX;
1260
    MeshFacet& rFace = _rclMesh._aclFacetArray[ulFacetPos];
1261

1262
    for (unsigned short i = 0; i < 3; i++) {
1263
        Base::Vector3f cBase(_rclMesh._aclPointArray[rFace._aulPoints[i]]);
1264
        Base::Vector3f cEnd(_rclMesh._aclPointArray[rFace._aulPoints[(i + 1) % 3]]);
1265
        Base::Vector3f cDir = cEnd - cBase;
1266

1267
        float fDist = rP.DistanceToLine(cBase, cDir);
1268
        if (fDist < fMinDist) {
1269
            fMinDist = fDist;
1270
            iEdgeNo = i;
1271
        }
1272
    }
1273

1274
    if (fMinDist < 0.05f) {
1275
        if (rFace._aulNeighbours[iEdgeNo] != FACET_INDEX_MAX) {
1276
            SplitEdge(ulFacetPos, rFace._aulNeighbours[iEdgeNo], rP);
1277
        }
1278
        else {
1279
            SplitOpenEdge(ulFacetPos, iEdgeNo, rP);
1280
        }
1281
    }
1282
}
1283

1284
void MeshTopoAlgorithm::SplitFacetOnTwoEdges(FacetIndex ulFacetPos,
1285
                                             const Base::Vector3f& rP1,
1286
                                             const Base::Vector3f& rP2)
1287
{
1288
    // search for the matching edges
1289
    unsigned short iEdgeNo1 = USHRT_MAX, iEdgeNo2 = USHRT_MAX;
1290
    float fMinDist1 = FLOAT_MAX, fMinDist2 = FLOAT_MAX;
1291
    MeshFacet& rFace = _rclMesh._aclFacetArray[ulFacetPos];
1292

1293
    for (unsigned short i = 0; i < 3; i++) {
1294
        Base::Vector3f cBase(_rclMesh._aclPointArray[rFace._aulPoints[i]]);
1295
        Base::Vector3f cEnd(_rclMesh._aclPointArray[rFace._aulPoints[(i + 1) % 3]]);
1296
        Base::Vector3f cDir = cEnd - cBase;
1297

1298
        float fDist = rP1.DistanceToLine(cBase, cDir);
1299
        if (fDist < fMinDist1) {
1300
            fMinDist1 = fDist;
1301
            iEdgeNo1 = i;
1302
        }
1303
        fDist = rP2.DistanceToLine(cBase, cDir);
1304
        if (fDist < fMinDist2) {
1305
            fMinDist2 = fDist;
1306
            iEdgeNo2 = i;
1307
        }
1308
    }
1309

1310
    if (iEdgeNo1 == iEdgeNo2 || fMinDist1 >= 0.05f || fMinDist2 >= 0.05f) {
1311
        return;  // no valid configuration
1312
    }
1313

1314
    // make first point lying on the previous edge
1315
    Base::Vector3f cP1 = rP1;
1316
    Base::Vector3f cP2 = rP2;
1317
    if ((iEdgeNo2 + 1) % 3 == iEdgeNo1) {
1318
        std::swap(iEdgeNo1, iEdgeNo2);
1319
        std::swap(cP1, cP2);
1320
    }
1321

1322
    // insert new points
1323
    PointIndex cntPts1 = this->GetOrAddIndex(cP1);
1324
    PointIndex cntPts2 = this->GetOrAddIndex(cP2);
1325
    FacetIndex cntFts = _rclMesh.CountFacets();
1326

1327
    unsigned short v0 = (iEdgeNo2 + 1) % 3;
1328
    unsigned short v1 = iEdgeNo1;
1329
    unsigned short v2 = iEdgeNo2;
1330

1331
    PointIndex p0 = rFace._aulPoints[v0];
1332
    PointIndex p1 = rFace._aulPoints[v1];
1333
    PointIndex p2 = rFace._aulPoints[v2];
1334

1335
    FacetIndex n0 = rFace._aulNeighbours[v0];
1336
    FacetIndex n1 = rFace._aulNeighbours[v1];
1337
    FacetIndex n2 = rFace._aulNeighbours[v2];
1338

1339
    // Modify and add facets
1340
    //
1341
    rFace._aulPoints[v0] = cntPts2;
1342
    rFace._aulPoints[v1] = cntPts1;
1343
    rFace._aulNeighbours[v0] = cntFts + 1;
1344

1345
    float dist1 = Base::DistanceP2(_rclMesh._aclPointArray[p0], cP1);
1346
    float dist2 = Base::DistanceP2(_rclMesh._aclPointArray[p1], cP2);
1347

1348
    if (dist1 > dist2) {
1349
        AddFacet(p0, p1, cntPts2, n0, cntFts + 1, n2);
1350
        AddFacet(p1, cntPts1, cntPts2, n1, ulFacetPos, cntFts);
1351
    }
1352
    else {
1353
        AddFacet(p0, p1, cntPts1, n0, n1, cntFts + 1);
1354
        AddFacet(p0, cntPts1, cntPts2, cntFts, ulFacetPos, n2);
1355
    }
1356

1357
    std::vector<FacetIndex> fixIndices;
1358
    fixIndices.push_back(ulFacetPos);
1359

1360
    if (n0 != FACET_INDEX_MAX) {
1361
        fixIndices.push_back(n0);
1362
    }
1363

1364
    // split up the neighbour facets
1365
    if (n1 != FACET_INDEX_MAX) {
1366
        fixIndices.push_back(n1);
1367
        MeshFacet& rN = _rclMesh._aclFacetArray[n1];
1368
        for (FacetIndex nbIndex : rN._aulNeighbours) {
1369
            fixIndices.push_back(nbIndex);
1370
        }
1371
        SplitFacet(n1, p1, p2, cntPts1);
1372
    }
1373

1374
    if (n2 != FACET_INDEX_MAX) {
1375
        fixIndices.push_back(n2);
1376
        MeshFacet& rN = _rclMesh._aclFacetArray[n2];
1377
        for (FacetIndex nbIndex : rN._aulNeighbours) {
1378
            fixIndices.push_back(nbIndex);
1379
        }
1380
        SplitFacet(n2, p2, p0, cntPts2);
1381
    }
1382

1383
    FacetIndex cntFts2 = _rclMesh.CountFacets();
1384
    for (FacetIndex i = cntFts; i < cntFts2; i++) {
1385
        fixIndices.push_back(i);
1386
    }
1387

1388
    std::sort(fixIndices.begin(), fixIndices.end());
1389
    fixIndices.erase(std::unique(fixIndices.begin(), fixIndices.end()), fixIndices.end());
1390
    HarmonizeNeighbours(fixIndices);
1391
}
1392

1393
void MeshTopoAlgorithm::SplitFacet(FacetIndex ulFacetPos,
1394
                                   PointIndex P1,
1395
                                   PointIndex P2,
1396
                                   PointIndex Pn)
1397
{
1398
    MeshFacet& rFace = _rclMesh._aclFacetArray[ulFacetPos];
1399
    unsigned short side = rFace.Side(P1, P2);
1400
    if (side != USHRT_MAX) {
1401
        PointIndex V1 = rFace._aulPoints[(side + 1) % 3];
1402
        PointIndex V2 = rFace._aulPoints[(side + 2) % 3];
1403
        FacetIndex size = _rclMesh._aclFacetArray.size();
1404

1405
        rFace._aulPoints[(side + 1) % 3] = Pn;
1406
        FacetIndex N1 = rFace._aulNeighbours[(side + 1) % 3];
1407
        if (N1 != FACET_INDEX_MAX) {
1408
            _rclMesh._aclFacetArray[N1].ReplaceNeighbour(ulFacetPos, size);
1409
        }
1410

1411
        rFace._aulNeighbours[(side + 1) % 3] = ulFacetPos;
1412
        AddFacet(Pn, V1, V2);
1413
    }
1414
}
1415

1416
void MeshTopoAlgorithm::AddFacet(PointIndex P1, PointIndex P2, PointIndex P3)
1417
{
1418
    MeshFacet facet;
1419
    facet._aulPoints[0] = P1;
1420
    facet._aulPoints[1] = P2;
1421
    facet._aulPoints[2] = P3;
1422

1423
    _rclMesh._aclFacetArray.push_back(facet);
1424
}
1425

1426
void MeshTopoAlgorithm::AddFacet(PointIndex P1,
1427
                                 PointIndex P2,
1428
                                 PointIndex P3,
1429
                                 FacetIndex N1,
1430
                                 FacetIndex N2,
1431
                                 FacetIndex N3)
1432
{
1433
    MeshFacet facet;
1434
    facet._aulPoints[0] = P1;
1435
    facet._aulPoints[1] = P2;
1436
    facet._aulPoints[2] = P3;
1437
    facet._aulNeighbours[0] = N1;
1438
    facet._aulNeighbours[1] = N2;
1439
    facet._aulNeighbours[2] = N3;
1440

1441
    _rclMesh._aclFacetArray.push_back(facet);
1442
}
1443

1444
void MeshTopoAlgorithm::HarmonizeNeighbours(const std::vector<FacetIndex>& ulFacets)
1445
{
1446
    for (FacetIndex it : ulFacets) {
1447
        for (FacetIndex jt : ulFacets) {
1448
            HarmonizeNeighbours(it, jt);
1449
        }
1450
    }
1451
}
1452

1453
void MeshTopoAlgorithm::HarmonizeNeighbours(FacetIndex facet1, FacetIndex facet2)
1454
{
1455
    if (facet1 == facet2) {
1456
        return;
1457
    }
1458

1459
    MeshFacet& rFace1 = _rclMesh._aclFacetArray[facet1];
1460
    MeshFacet& rFace2 = _rclMesh._aclFacetArray[facet2];
1461

1462
    unsigned short side = rFace1.Side(rFace2);
1463
    if (side != USHRT_MAX) {
1464
        rFace1._aulNeighbours[side] = facet2;
1465
    }
1466

1467
    side = rFace2.Side(rFace1);
1468
    if (side != USHRT_MAX) {
1469
        rFace2._aulNeighbours[side] = facet1;
1470
    }
1471
}
1472

1473
void MeshTopoAlgorithm::SplitNeighbourFacet(FacetIndex ulFacetPos,
1474
                                            unsigned short uFSide,
1475
                                            const Base::Vector3f rPoint)
1476
{
1477
    MeshFacet& rclF = _rclMesh._aclFacetArray[ulFacetPos];
1478

1479
    FacetIndex ulNeighbour = rclF._aulNeighbours[uFSide];
1480
    MeshFacet& rclN = _rclMesh._aclFacetArray[ulNeighbour];
1481

1482
    unsigned short uNSide = rclN.Side(rclF);
1483

1484
    PointIndex uPtInd = this->GetOrAddIndex(rPoint);
1485
    FacetIndex ulSize = _rclMesh._aclFacetArray.size();
1486

1487
    // adjust the neighbourhood
1488
    if (rclN._aulNeighbours[(uNSide + 1) % 3] != FACET_INDEX_MAX) {
1489
        _rclMesh._aclFacetArray[rclN._aulNeighbours[(uNSide + 1) % 3]].ReplaceNeighbour(ulNeighbour,
1490
                                                                                        ulSize);
1491
    }
1492

1493
    MeshFacet cNew;
1494
    cNew._aulPoints[0] = uPtInd;
1495
    cNew._aulPoints[1] = rclN._aulPoints[(uNSide + 1) % 3];
1496
    cNew._aulPoints[2] = rclN._aulPoints[(uNSide + 2) % 3];
1497
    cNew._aulNeighbours[0] = ulFacetPos;
1498
    cNew._aulNeighbours[1] = rclN._aulNeighbours[(uNSide + 1) % 3];
1499
    cNew._aulNeighbours[2] = ulNeighbour;
1500

1501
    // adjust the facet
1502
    rclN._aulPoints[(uNSide + 1) % 3] = uPtInd;
1503
    rclN._aulNeighbours[(uNSide + 1) % 3] = ulSize;
1504

1505
    // insert new facet
1506
    _rclMesh._aclFacetArray.push_back(cNew);
1507
}
1508

1509
#if 0
1510
  // create 3 new facets
1511
  MeshGeomFacet clFacet;
1512

1513
  // facet [P1, Ei+1, P2]
1514
  clFacet._aclPoints[0] = cP1;
1515
  clFacet._aclPoints[1] = _rclMesh._aclPointArray[rFace._aulPoints[(iEdgeNo1+1)%3]];
1516
  clFacet._aclPoints[2] = cP2;
1517
  clFacet.CalcNormal();
1518
  _aclNewFacets.push_back(clFacet);
1519
  // facet [P2, Ei+2, Ei]
1520
  clFacet._aclPoints[0] = cP2;
1521
  clFacet._aclPoints[1] = _rclMesh._aclPointArray[rFace._aulPoints[(iEdgeNo1+2)%3]];
1522
  clFacet._aclPoints[2] = _rclMesh._aclPointArray[rFace._aulPoints[iEdgeNo1]];
1523
  clFacet.CalcNormal();
1524
  _aclNewFacets.push_back(clFacet);
1525
  // facet [P2, Ei, P1]
1526
  clFacet._aclPoints[0] = cP2;
1527
  clFacet._aclPoints[1] = _rclMesh._aclPointArray[rFace._aulPoints[iEdgeNo1]];
1528
  clFacet._aclPoints[2] = cP1;
1529
  clFacet.CalcNormal();
1530
  _aclNewFacets.push_back(clFacet);
1531
#endif
1532

1533
bool MeshTopoAlgorithm::RemoveDegeneratedFacet(FacetIndex index)
1534
{
1535
    if (index >= _rclMesh._aclFacetArray.size()) {
1536
        return false;
1537
    }
1538
    MeshFacet& rFace = _rclMesh._aclFacetArray[index];
1539

1540
    // coincident corners (either topological or geometrical)
1541
    for (int i = 0; i < 3; i++) {
1542
        const MeshPoint& rE0 = _rclMesh._aclPointArray[rFace._aulPoints[i]];
1543
        const MeshPoint& rE1 = _rclMesh._aclPointArray[rFace._aulPoints[(i + 1) % 3]];
1544
        if (rE0 == rE1) {
1545
            FacetIndex uN1 = rFace._aulNeighbours[(i + 1) % 3];
1546
            FacetIndex uN2 = rFace._aulNeighbours[(i + 2) % 3];
1547
            if (uN2 != FACET_INDEX_MAX) {
1548
                _rclMesh._aclFacetArray[uN2].ReplaceNeighbour(index, uN1);
1549
            }
1550
            if (uN1 != FACET_INDEX_MAX) {
1551
                _rclMesh._aclFacetArray[uN1].ReplaceNeighbour(index, uN2);
1552
            }
1553

1554
            // isolate the face and remove it
1555
            rFace._aulNeighbours[0] = FACET_INDEX_MAX;
1556
            rFace._aulNeighbours[1] = FACET_INDEX_MAX;
1557
            rFace._aulNeighbours[2] = FACET_INDEX_MAX;
1558
            _rclMesh.DeleteFacet(index);
1559
            return true;
1560
        }
1561
    }
1562

1563
    // We have a facet of the form
1564
    // P0 +----+------+P2
1565
    //         P1
1566
    for (int j = 0; j < 3; j++) {
1567
        Base::Vector3f cVec1 = _rclMesh._aclPointArray[rFace._aulPoints[(j + 1) % 3]]
1568
            - _rclMesh._aclPointArray[rFace._aulPoints[j]];
1569
        Base::Vector3f cVec2 = _rclMesh._aclPointArray[rFace._aulPoints[(j + 2) % 3]]
1570
            - _rclMesh._aclPointArray[rFace._aulPoints[j]];
1571

1572
        // adjust the neighbourhoods and point indices
1573
        if (cVec1 * cVec2 < 0.0f) {
1574
            FacetIndex uN1 = rFace._aulNeighbours[(j + 1) % 3];
1575
            if (uN1 != FACET_INDEX_MAX) {
1576
                // get the neighbour and common edge side
1577
                MeshFacet& rNb = _rclMesh._aclFacetArray[uN1];
1578
                unsigned short side = rNb.Side(index);
1579

1580
                // bend the point indices
1581
                rFace._aulPoints[(j + 2) % 3] = rNb._aulPoints[(side + 2) % 3];
1582
                rNb._aulPoints[(side + 1) % 3] = rFace._aulPoints[j];
1583

1584
                // set correct neighbourhood
1585
                FacetIndex uN2 = rFace._aulNeighbours[(j + 2) % 3];
1586
                rNb._aulNeighbours[side] = uN2;
1587
                if (uN2 != FACET_INDEX_MAX) {
1588
                    _rclMesh._aclFacetArray[uN2].ReplaceNeighbour(index, uN1);
1589
                }
1590
                FacetIndex uN3 = rNb._aulNeighbours[(side + 1) % 3];
1591
                rFace._aulNeighbours[(j + 1) % 3] = uN3;
1592
                if (uN3 != FACET_INDEX_MAX) {
1593
                    _rclMesh._aclFacetArray[uN3].ReplaceNeighbour(uN1, index);
1594
                }
1595
                rNb._aulNeighbours[(side + 1) % 3] = index;
1596
                rFace._aulNeighbours[(j + 2) % 3] = uN1;
1597
            }
1598
            else {
1599
                _rclMesh.DeleteFacet(index);
1600
            }
1601

1602
            return true;
1603
        }
1604
    }
1605

1606
    return false;
1607
}
1608

1609
bool MeshTopoAlgorithm::RemoveCorruptedFacet(FacetIndex index)
1610
{
1611
    if (index >= _rclMesh._aclFacetArray.size()) {
1612
        return false;
1613
    }
1614
    MeshFacet& rFace = _rclMesh._aclFacetArray[index];
1615

1616
    // coincident corners (topological)
1617
    for (int i = 0; i < 3; i++) {
1618
        if (rFace._aulPoints[i] == rFace._aulPoints[(i + 1) % 3]) {
1619
            FacetIndex uN1 = rFace._aulNeighbours[(i + 1) % 3];
1620
            FacetIndex uN2 = rFace._aulNeighbours[(i + 2) % 3];
1621
            if (uN2 != FACET_INDEX_MAX) {
1622
                _rclMesh._aclFacetArray[uN2].ReplaceNeighbour(index, uN1);
1623
            }
1624
            if (uN1 != FACET_INDEX_MAX) {
1625
                _rclMesh._aclFacetArray[uN1].ReplaceNeighbour(index, uN2);
1626
            }
1627

1628
            // isolate the face and remove it
1629
            rFace._aulNeighbours[0] = FACET_INDEX_MAX;
1630
            rFace._aulNeighbours[1] = FACET_INDEX_MAX;
1631
            rFace._aulNeighbours[2] = FACET_INDEX_MAX;
1632
            _rclMesh.DeleteFacet(index);
1633
            return true;
1634
        }
1635
    }
1636

1637
    return false;
1638
}
1639

1640
void MeshTopoAlgorithm::FillupHoles(unsigned long length,
1641
                                    int level,
1642
                                    AbstractPolygonTriangulator& cTria,
1643
                                    std::list<std::vector<PointIndex>>& aFailed)
1644
{
1645
    // get the mesh boundaries as an array of point indices
1646
    std::list<std::vector<PointIndex>> aBorders, aFillBorders;
1647
    MeshAlgorithm cAlgo(_rclMesh);
1648
    cAlgo.GetMeshBorders(aBorders);
1649

1650
    // split boundary loops if needed
1651
    cAlgo.SplitBoundaryLoops(aBorders);
1652

1653
    for (const auto& aBorder : aBorders) {
1654
        if (aBorder.size() - 1 <= length) {  // ignore boundary with too many edges
1655
            aFillBorders.push_back(aBorder);
1656
        }
1657
    }
1658

1659
    if (!aFillBorders.empty()) {
1660
        FillupHoles(level, cTria, aFillBorders, aFailed);
1661
    }
1662
}
1663

1664
void MeshTopoAlgorithm::FillupHoles(int level,
1665
                                    AbstractPolygonTriangulator& cTria,
1666
                                    const std::list<std::vector<PointIndex>>& aBorders,
1667
                                    std::list<std::vector<PointIndex>>& aFailed)
1668
{
1669
    // get the facets to a point
1670
    MeshRefPointToFacets cPt2Fac(_rclMesh);
1671
    MeshAlgorithm cAlgo(_rclMesh);
1672

1673
    MeshFacetArray newFacets;
1674
    MeshPointArray newPoints;
1675
    unsigned long numberOfOldPoints = _rclMesh._aclPointArray.size();
1676
    for (const auto& aBorder : aBorders) {
1677
        MeshFacetArray cFacets;
1678
        MeshPointArray cPoints;
1679
        std::vector<PointIndex> bound = aBorder;
1680
        if (cAlgo.FillupHole(bound, cTria, cFacets, cPoints, level, &cPt2Fac)) {
1681
            if (bound.front() == bound.back()) {
1682
                bound.pop_back();
1683
            }
1684
            // the triangulation may produce additional points which we must take into account when
1685
            // appending to the mesh
1686
            if (cPoints.size() > bound.size()) {
1687
                unsigned long countBoundaryPoints = bound.size();
1688
                unsigned long countDifference = cPoints.size() - countBoundaryPoints;
1689
                MeshPointArray::_TIterator pt = cPoints.begin() + countBoundaryPoints;
1690
                for (unsigned long i = 0; i < countDifference; i++, pt++) {
1691
                    bound.push_back(numberOfOldPoints++);
1692
                    newPoints.push_back(*pt);
1693
                }
1694
            }
1695
            if (cTria.NeedsReindexing()) {
1696
                for (auto& cFacet : cFacets) {
1697
                    cFacet._aulPoints[0] = bound[cFacet._aulPoints[0]];
1698
                    cFacet._aulPoints[1] = bound[cFacet._aulPoints[1]];
1699
                    cFacet._aulPoints[2] = bound[cFacet._aulPoints[2]];
1700
                    newFacets.push_back(cFacet);
1701
                }
1702
            }
1703
            else {
1704
                for (auto& cFacet : cFacets) {
1705
                    newFacets.push_back(cFacet);
1706
                }
1707
            }
1708
        }
1709
        else {
1710
            aFailed.push_back(aBorder);
1711
        }
1712
    }
1713

1714
    // insert new points and faces into the mesh structure
1715
    _rclMesh._aclPointArray.insert(_rclMesh._aclPointArray.end(),
1716
                                   newPoints.begin(),
1717
                                   newPoints.end());
1718
    for (const auto& newPoint : newPoints) {
1719
        _rclMesh._clBoundBox.Add(newPoint);
1720
    }
1721
    if (!newFacets.empty()) {
1722
        // Do some checks for invalid point indices
1723
        MeshFacetArray addFacets;
1724
        addFacets.reserve(newFacets.size());
1725
        unsigned long ctPoints = _rclMesh.CountPoints();
1726
        for (auto& newFacet : newFacets) {
1727
            if (newFacet._aulPoints[0] >= ctPoints || newFacet._aulPoints[1] >= ctPoints
1728
                || newFacet._aulPoints[2] >= ctPoints) {
1729
                Base::Console().Log("Ignore invalid face <%d, %d, %d> (%d vertices)\n",
1730
                                    newFacet._aulPoints[0],
1731
                                    newFacet._aulPoints[1],
1732
                                    newFacet._aulPoints[2],
1733
                                    ctPoints);
1734
            }
1735
            else {
1736
                addFacets.push_back(newFacet);
1737
            }
1738
        }
1739
        _rclMesh.AddFacets(addFacets, true);
1740
    }
1741
}
1742

1743
void MeshTopoAlgorithm::FindHoles(unsigned long length,
1744
                                  std::list<std::vector<PointIndex>>& aBorders) const
1745
{
1746
    std::list<std::vector<PointIndex>> border;
1747
    MeshAlgorithm cAlgo(_rclMesh);
1748
    cAlgo.GetMeshBorders(border);
1749
    for (const auto& it : border) {
1750
        if (it.size() <= length) {
1751
            aBorders.push_back(it);
1752
        }
1753
    }
1754
}
1755

1756
void MeshTopoAlgorithm::FindComponents(unsigned long count, std::vector<FacetIndex>& findIndices)
1757
{
1758
    std::vector<std::vector<FacetIndex>> segments;
1759
    MeshComponents comp(_rclMesh);
1760
    comp.SearchForComponents(MeshComponents::OverEdge, segments);
1761

1762
    for (const auto& segment : segments) {
1763
        if (segment.size() <= count) {
1764
            findIndices.insert(findIndices.end(), segment.begin(), segment.end());
1765
        }
1766
    }
1767
}
1768

1769
void MeshTopoAlgorithm::RemoveComponents(unsigned long count)
1770
{
1771
    std::vector<FacetIndex> removeFacets;
1772
    FindComponents(count, removeFacets);
1773
    if (!removeFacets.empty()) {
1774
        _rclMesh.DeleteFacets(removeFacets);
1775
    }
1776
}
1777

1778
void MeshTopoAlgorithm::HarmonizeNormals()
1779
{
1780
    std::vector<FacetIndex> uIndices = MeshEvalOrientation(_rclMesh).GetIndices();
1781
    for (FacetIndex index : uIndices) {
1782
        _rclMesh._aclFacetArray[index].FlipNormal();
1783
    }
1784
}
1785

1786
void MeshTopoAlgorithm::FlipNormals()
1787
{
1788
    for (MeshFacetArray::_TIterator i = _rclMesh._aclFacetArray.begin();
1789
         i < _rclMesh._aclFacetArray.end();
1790
         ++i) {
1791
        i->FlipNormal();
1792
    }
1793
}
1794

1795
// ---------------------------------------------------------------------------
1796

1797
/**
1798
 * Some important formulas:
1799
 *
1800
 * Ne = 3Nv - Nb + 3B + 6(G-R)
1801
 * Nt = 2Nv - Nb + 2B + 4(G-R)
1802
 *
1803
 * Ne <= 3Nv + 6(G-R)
1804
 * Nt <= 2Nv + 4(G-R)
1805
 *
1806
 * Ne ~ 3Nv, Nv >> G, Nv >> R
1807
 * Nt ~ 2Nv, Nv >> G, Nv >> R
1808
 *
1809
 * Ne = #Edges
1810
 * Nt = #Facets
1811
 * Nv = #Vertices
1812
 * Nb = #Boundary vertices
1813
 * B  = #Boundaries
1814
 * G  = Genus (Number of holes)
1815
 * R  = #components
1816
 *
1817
 * See also http://max-limper.de/publications/Euler/
1818
 */
1819

1820
MeshComponents::MeshComponents(const MeshKernel& rclMesh)
1821
    : _rclMesh(rclMesh)
1822
{}
1823

1824
void MeshComponents::SearchForComponents(TMode tMode,
1825
                                         std::vector<std::vector<FacetIndex>>& aclT) const
1826
{
1827
    // all facets
1828
    std::vector<FacetIndex> aulAllFacets(_rclMesh.CountFacets());
1829
    FacetIndex k = 0;
1830
    for (FacetIndex& aulAllFacet : aulAllFacets) {
1831
        aulAllFacet = k++;
1832
    }
1833

1834
    SearchForComponents(tMode, aulAllFacets, aclT);
1835
}
1836

1837
void MeshComponents::SearchForComponents(TMode tMode,
1838
                                         const std::vector<FacetIndex>& aSegment,
1839
                                         std::vector<std::vector<FacetIndex>>& aclT) const
1840
{
1841
    FacetIndex ulStartFacet {};
1842

1843
    if (_rclMesh.CountFacets() == 0) {
1844
        return;
1845
    }
1846

1847
    // reset VISIT flags
1848
    MeshAlgorithm cAlgo(_rclMesh);
1849
    cAlgo.SetFacetFlag(MeshFacet::VISIT);
1850
    cAlgo.ResetFacetsFlag(aSegment, MeshFacet::VISIT);
1851

1852
    const MeshFacetArray& rFAry = _rclMesh.GetFacets();
1853
    MeshFacetArray::_TConstIterator iTri = rFAry.begin();
1854
    MeshFacetArray::_TConstIterator iBeg = rFAry.begin();
1855
    MeshFacetArray::_TConstIterator iEnd = rFAry.end();
1856

1857
    // start from the first not visited facet
1858
    unsigned long ulVisited = cAlgo.CountFacetFlag(MeshFacet::VISIT);
1859
    MeshIsNotFlag<MeshFacet> flag;
1860
    iTri = std::find_if(iTri, iEnd, [flag](const MeshFacet& f) {
1861
        return flag(f, MeshFacet::VISIT);
1862
    });
1863
    ulStartFacet = iTri - iBeg;
1864

1865
    // visitor
1866
    std::vector<FacetIndex> aclComponent;
1867
    std::vector<std::vector<FacetIndex>> aclConnectComp;
1868
    MeshTopFacetVisitor clFVisitor(aclComponent);
1869

1870
    while (ulStartFacet != FACET_INDEX_MAX) {
1871
        // collect all facets of a component
1872
        aclComponent.clear();
1873
        if (tMode == OverEdge) {
1874
            ulVisited += _rclMesh.VisitNeighbourFacets(clFVisitor, ulStartFacet);
1875
        }
1876
        else if (tMode == OverPoint) {
1877
            ulVisited += _rclMesh.VisitNeighbourFacetsOverCorners(clFVisitor, ulStartFacet);
1878
        }
1879

1880
        // get also start facet
1881
        aclComponent.push_back(ulStartFacet);
1882
        aclConnectComp.push_back(aclComponent);
1883

1884
        // if the mesh consists of several topologic independent components
1885
        // We can search from position 'iTri' on because all elements _before_ are already visited
1886
        // what we know from the previous iteration.
1887
        iTri = std::find_if(iTri, iEnd, [flag](const MeshFacet& f) {
1888
            return flag(f, MeshFacet::VISIT);
1889
        });
1890

1891
        if (iTri < iEnd) {
1892
            ulStartFacet = iTri - iBeg;
1893
        }
1894
        else {
1895
            ulStartFacet = FACET_INDEX_MAX;
1896
        }
1897
    }
1898

1899
    // sort components by size (descending order)
1900
    std::sort(aclConnectComp.begin(), aclConnectComp.end(), CNofFacetsCompare());
1901
    aclT = aclConnectComp;
1902
    boost::ignore_unused(ulVisited);
1903
}
1904
FreeCAD

Использование cookies
