FreeCAD

1
/***************************************************************************
2
 *   Copyright (c) 2012 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
#ifndef _PreComp_
25
#include <algorithm>
26
#endif
27

28
#include "Algorithm.h"
29
#include "Approximation.h"
30
#include "Segmentation.h"
31

32
using namespace MeshCore;
33

34
void MeshSurfaceSegment::Initialize(FacetIndex)
35
{}
36

37
bool MeshSurfaceSegment::TestInitialFacet(FacetIndex) const
38
{
39
    return true;
40
}
41

42
void MeshSurfaceSegment::AddFacet(const MeshFacet&)
43
{}
44

45
void MeshSurfaceSegment::AddSegment(const std::vector<FacetIndex>& segm)
46
{
47
    if (segm.size() >= minFacets) {
48
        segments.push_back(segm);
49
    }
50
}
51

52
MeshSegment MeshSurfaceSegment::FindSegment(FacetIndex index) const
53
{
54
    for (const auto& segment : segments) {
55
        if (std::find(segment.begin(), segment.end(), index) != segment.end()) {
56
            return segment;
57
        }
58
    }
59

60
    return {};
61
}
62

63
// --------------------------------------------------------
64

65
MeshDistancePlanarSegment::MeshDistancePlanarSegment(const MeshKernel& mesh,
66
                                                     unsigned long minFacets,
67
                                                     float tol)
68
    : MeshDistanceSurfaceSegment(mesh, minFacets, tol)
69
    , fitter(new PlaneFit)
70
{}
71

72
MeshDistancePlanarSegment::~MeshDistancePlanarSegment()
73
{
74
    delete fitter;
75
}
76

77
void MeshDistancePlanarSegment::Initialize(FacetIndex index)
78
{
79
    fitter->Clear();
80

81
    MeshGeomFacet triangle = kernel.GetFacet(index);
82
    basepoint = triangle.GetGravityPoint();
83
    normal = triangle.GetNormal();
84
    fitter->AddPoint(triangle._aclPoints[0]);
85
    fitter->AddPoint(triangle._aclPoints[1]);
86
    fitter->AddPoint(triangle._aclPoints[2]);
87
}
88

89
bool MeshDistancePlanarSegment::TestFacet(const MeshFacet& face) const
90
{
91
    if (!fitter->Done()) {
92
        fitter->Fit();
93
    }
94
    MeshGeomFacet triangle = kernel.GetFacet(face);
95
    for (auto pnt : triangle._aclPoints) {
96
        if (fabs(fitter->GetDistanceToPlane(pnt)) > tolerance) {
97
            return false;
98
        }
99
    }
100

101
    return true;
102
}
103

104
void MeshDistancePlanarSegment::AddFacet(const MeshFacet& face)
105
{
106
    MeshGeomFacet triangle = kernel.GetFacet(face);
107
    fitter->AddPoint(triangle.GetGravityPoint());
108
}
109

110
// --------------------------------------------------------
111

112
PlaneSurfaceFit::PlaneSurfaceFit()
113
    : fitter(new PlaneFit)
114
{}
115

116
PlaneSurfaceFit::PlaneSurfaceFit(const Base::Vector3f& b, const Base::Vector3f& n)
117
    : basepoint(b)
118
    , normal(n)
119
    , fitter(nullptr)
120
{}
121

122
PlaneSurfaceFit::~PlaneSurfaceFit()
123
{
124
    delete fitter;
125
}
126

127
void PlaneSurfaceFit::Initialize(const MeshCore::MeshGeomFacet& tria)
128
{
129
    if (fitter) {
130
        basepoint = tria.GetGravityPoint();
131
        normal = tria.GetNormal();
132

133
        fitter->Clear();
134

135
        fitter->AddPoint(tria._aclPoints[0]);
136
        fitter->AddPoint(tria._aclPoints[1]);
137
        fitter->AddPoint(tria._aclPoints[2]);
138
        fitter->Fit();
139
    }
140
}
141

142
bool PlaneSurfaceFit::TestTriangle(const MeshGeomFacet&) const
143
{
144
    return true;
145
}
146

147
void PlaneSurfaceFit::AddTriangle(const MeshCore::MeshGeomFacet& tria)
148
{
149
    if (fitter) {
150
        fitter->AddPoint(tria.GetGravityPoint());
151
    }
152
}
153

154
bool PlaneSurfaceFit::Done() const
155
{
156
    if (!fitter) {
157
        return true;
158
    }
159
    else {
160
        return fitter->Done();
161
    }
162
}
163

164
float PlaneSurfaceFit::Fit()
165
{
166
    if (!fitter) {
167
        return 0;
168
    }
169
    else {
170
        return fitter->Fit();
171
    }
172
}
173

174
float PlaneSurfaceFit::GetDistanceToSurface(const Base::Vector3f& pnt) const
175
{
176
    if (!fitter) {
177
        return pnt.DistanceToPlane(basepoint, normal);
178
    }
179
    else {
180
        return fitter->GetDistanceToPlane(pnt);
181
    }
182
}
183

184
std::vector<float> PlaneSurfaceFit::Parameters() const
185
{
186
    Base::Vector3f base = basepoint;
187
    Base::Vector3f norm = normal;
188
    if (fitter) {
189
        base = fitter->GetBase();
190
        norm = fitter->GetNormal();
191
    }
192

193
    std::vector<float> c;
194
    c.push_back(base.x);
195
    c.push_back(base.y);
196
    c.push_back(base.z);
197
    c.push_back(norm.x);
198
    c.push_back(norm.y);
199
    c.push_back(norm.z);
200
    return c;
201
}
202

203
// --------------------------------------------------------
204

205
CylinderSurfaceFit::CylinderSurfaceFit()
206
    : radius(FLOAT_MAX)
207
    , fitter(new CylinderFit)
208
{
209
    axis.Set(0, 0, 0);
210
}
211

212
/*!
213
 * \brief CylinderSurfaceFit::CylinderSurfaceFit
214
 * Set a predefined cylinder. Internal cylinder fits are not done, then.
215
 */
216
CylinderSurfaceFit::CylinderSurfaceFit(const Base::Vector3f& b, const Base::Vector3f& a, float r)
217
    : basepoint(b)
218
    , axis(a)
219
    , radius(r)
220
    , fitter(nullptr)
221
{}
222

223
CylinderSurfaceFit::~CylinderSurfaceFit()
224
{
225
    delete fitter;
226
}
227

228
void CylinderSurfaceFit::Initialize(const MeshCore::MeshGeomFacet& tria)
229
{
230
    if (fitter) {
231
        fitter->Clear();
232
        fitter->AddPoint(tria._aclPoints[0]);
233
        fitter->AddPoint(tria._aclPoints[1]);
234
        fitter->AddPoint(tria._aclPoints[2]);
235
    }
236
}
237

238
void CylinderSurfaceFit::AddTriangle(const MeshCore::MeshGeomFacet& tria)
239
{
240
    if (fitter) {
241
        fitter->AddPoint(tria._aclPoints[0]);
242
        fitter->AddPoint(tria._aclPoints[1]);
243
        fitter->AddPoint(tria._aclPoints[2]);
244
    }
245
}
246

247
bool CylinderSurfaceFit::TestTriangle(const MeshGeomFacet& tria) const
248
{
249
    // This is to filter out triangles whose points lie on the cylinder and
250
    // that whose normals are more or less parallel to the cylinder axis
251
    float dot = axis.Dot(tria.GetNormal());
252
    return fabs(dot) < 0.5f;
253
}
254

255
bool CylinderSurfaceFit::Done() const
256
{
257
    if (fitter) {
258
        return fitter->Done();
259
    }
260

261
    return true;
262
}
263

264
float CylinderSurfaceFit::Fit()
265
{
266
    if (!fitter) {
267
        return 0;
268
    }
269

270
    float fit = fitter->Fit();
271
    if (fit < FLOAT_MAX) {
272
        basepoint = fitter->GetBase();
273
        axis = fitter->GetAxis();
274
        radius = fitter->GetRadius();
275
    }
276
    return fit;
277
}
278

279
float CylinderSurfaceFit::GetDistanceToSurface(const Base::Vector3f& pnt) const
280
{
281
    if (fitter && !fitter->Done()) {
282
        // collect some points
283
        return 0;
284
    }
285
    float dist = pnt.DistanceToLine(basepoint, axis);
286
    return (dist - radius);
287
}
288

289
std::vector<float> CylinderSurfaceFit::Parameters() const
290
{
291
    Base::Vector3f base = basepoint;
292
    Base::Vector3f norm = axis;
293
    float radval = radius;
294
    if (fitter) {
295
        base = fitter->GetBase();
296
        norm = fitter->GetAxis();
297
        radval = fitter->GetRadius();
298
    }
299

300
    std::vector<float> c;
301
    c.push_back(base.x);
302
    c.push_back(base.y);
303
    c.push_back(base.z);
304
    c.push_back(norm.x);
305
    c.push_back(norm.y);
306
    c.push_back(norm.z);
307
    c.push_back(radval);
308
    return c;
309
}
310

311
// --------------------------------------------------------
312

313
SphereSurfaceFit::SphereSurfaceFit()
314
    : radius(FLOAT_MAX)
315
    , fitter(new SphereFit)
316
{
317
    center.Set(0, 0, 0);
318
}
319

320
SphereSurfaceFit::SphereSurfaceFit(const Base::Vector3f& c, float r)
321
    : center(c)
322
    , radius(r)
323
    , fitter(nullptr)
324
{}
325

326
SphereSurfaceFit::~SphereSurfaceFit()
327
{
328
    delete fitter;
329
}
330

331
void SphereSurfaceFit::Initialize(const MeshCore::MeshGeomFacet& tria)
332
{
333
    if (fitter) {
334
        fitter->Clear();
335
        fitter->AddPoint(tria._aclPoints[0]);
336
        fitter->AddPoint(tria._aclPoints[1]);
337
        fitter->AddPoint(tria._aclPoints[2]);
338
    }
339
}
340

341
void SphereSurfaceFit::AddTriangle(const MeshCore::MeshGeomFacet& tria)
342
{
343
    if (fitter) {
344
        fitter->AddPoint(tria._aclPoints[0]);
345
        fitter->AddPoint(tria._aclPoints[1]);
346
        fitter->AddPoint(tria._aclPoints[2]);
347
    }
348
}
349

350
bool SphereSurfaceFit::TestTriangle(const MeshGeomFacet&) const
351
{
352
    // Already handled by GetDistanceToSurface
353
    return true;
354
}
355

356
bool SphereSurfaceFit::Done() const
357
{
358
    if (fitter) {
359
        return fitter->Done();
360
    }
361

362
    return true;
363
}
364

365
float SphereSurfaceFit::Fit()
366
{
367
    if (!fitter) {
368
        return 0;
369
    }
370

371
    float fit = fitter->Fit();
372
    if (fit < FLOAT_MAX) {
373
        center = fitter->GetCenter();
374
        radius = fitter->GetRadius();
375
    }
376
    return fit;
377
}
378

379
float SphereSurfaceFit::GetDistanceToSurface(const Base::Vector3f& pnt) const
380
{
381
    float dist = Base::Distance(pnt, center);
382
    return (dist - radius);
383
}
384

385
std::vector<float> SphereSurfaceFit::Parameters() const
386
{
387
    Base::Vector3f base = center;
388
    float radval = radius;
389
    if (fitter) {
390
        base = fitter->GetCenter();
391
        radval = fitter->GetRadius();
392
    }
393

394
    std::vector<float> c;
395
    c.push_back(base.x);
396
    c.push_back(base.y);
397
    c.push_back(base.z);
398
    c.push_back(radval);
399
    return c;
400
}
401

402
// --------------------------------------------------------
403

404
MeshDistanceGenericSurfaceFitSegment::MeshDistanceGenericSurfaceFitSegment(AbstractSurfaceFit* fit,
405
                                                                           const MeshKernel& mesh,
406
                                                                           unsigned long minFacets,
407
                                                                           float tol)
408
    : MeshDistanceSurfaceSegment(mesh, minFacets, tol)
409
    , fitter(fit)
410
{}
411

412
MeshDistanceGenericSurfaceFitSegment::~MeshDistanceGenericSurfaceFitSegment()
413
{
414
    delete fitter;
415
}
416

417
void MeshDistanceGenericSurfaceFitSegment::Initialize(FacetIndex index)
418
{
419
    MeshGeomFacet triangle = kernel.GetFacet(index);
420
    fitter->Initialize(triangle);
421
}
422

423
bool MeshDistanceGenericSurfaceFitSegment::TestInitialFacet(FacetIndex index) const
424
{
425
    MeshGeomFacet triangle = kernel.GetFacet(index);
426
    for (auto pnt : triangle._aclPoints) {
427
        if (fabs(fitter->GetDistanceToSurface(pnt)) > tolerance) {
428
            return false;
429
        }
430
    }
431
    return fitter->TestTriangle(triangle);
432
}
433

434
bool MeshDistanceGenericSurfaceFitSegment::TestFacet(const MeshFacet& face) const
435
{
436
    if (!fitter->Done()) {
437
        fitter->Fit();
438
    }
439
    MeshGeomFacet triangle = kernel.GetFacet(face);
440
    for (auto ptIndex : triangle._aclPoints) {
441
        if (fabs(fitter->GetDistanceToSurface(ptIndex)) > tolerance) {
442
            return false;
443
        }
444
    }
445

446
    return fitter->TestTriangle(triangle);
447
}
448

449
void MeshDistanceGenericSurfaceFitSegment::AddFacet(const MeshFacet& face)
450
{
451
    MeshGeomFacet triangle = kernel.GetFacet(face);
452
    fitter->AddTriangle(triangle);
453
}
454

455
std::vector<float> MeshDistanceGenericSurfaceFitSegment::Parameters() const
456
{
457
    return fitter->Parameters();
458
}
459

460
// --------------------------------------------------------
461

462
bool MeshCurvaturePlanarSegment::TestFacet(const MeshFacet& rclFacet) const
463
{
464
    for (PointIndex ptIndex : rclFacet._aulPoints) {
465
        const CurvatureInfo& ci = GetInfo(ptIndex);
466
        if (fabs(ci.fMinCurvature) > tolerance) {
467
            return false;
468
        }
469
        if (fabs(ci.fMaxCurvature) > tolerance) {
470
            return false;
471
        }
472
    }
473

474
    return true;
475
}
476

477
bool MeshCurvatureCylindricalSegment::TestFacet(const MeshFacet& rclFacet) const
478
{
479
    for (PointIndex ptIndex : rclFacet._aulPoints) {
480
        const CurvatureInfo& ci = GetInfo(ptIndex);
481
        float fMax = std::max<float>(fabs(ci.fMaxCurvature), fabs(ci.fMinCurvature));
482
        float fMin = std::min<float>(fabs(ci.fMaxCurvature), fabs(ci.fMinCurvature));
483
        if (fMin > toleranceMin) {
484
            return false;
485
        }
486
        if (fabs(fMax - curvature) > toleranceMax) {
487
            return false;
488
        }
489
    }
490

491
    return true;
492
}
493

494
bool MeshCurvatureSphericalSegment::TestFacet(const MeshFacet& rclFacet) const
495
{
496
    for (PointIndex ptIndex : rclFacet._aulPoints) {
497
        const CurvatureInfo& ci = GetInfo(ptIndex);
498
        if (ci.fMaxCurvature * ci.fMinCurvature < 0) {
499
            return false;
500
        }
501
        float diff {};
502
        diff = fabs(ci.fMinCurvature) - curvature;
503
        if (fabs(diff) > tolerance) {
504
            return false;
505
        }
506
        diff = fabs(ci.fMaxCurvature) - curvature;
507
        if (fabs(diff) > tolerance) {
508
            return false;
509
        }
510
    }
511

512
    return true;
513
}
514

515
bool MeshCurvatureFreeformSegment::TestFacet(const MeshFacet& rclFacet) const
516
{
517
    for (PointIndex ptIndex : rclFacet._aulPoints) {
518
        const CurvatureInfo& ci = GetInfo(ptIndex);
519
        if (fabs(ci.fMinCurvature - c2) > toleranceMin) {
520
            return false;
521
        }
522
        if (fabs(ci.fMaxCurvature - c1) > toleranceMax) {
523
            return false;
524
        }
525
    }
526

527
    return true;
528
}
529

530
// --------------------------------------------------------
531

532
MeshSurfaceVisitor::MeshSurfaceVisitor(MeshSurfaceSegment& segm, std::vector<FacetIndex>& indices)
533
    : indices(indices)
534
    , segm(segm)
535
{}
536

537
bool MeshSurfaceVisitor::AllowVisit(const MeshFacet& face,
538
                                    const MeshFacet&,
539
                                    FacetIndex,
540
                                    unsigned long,
541
                                    unsigned short)
542
{
543
    return segm.TestFacet(face);
544
}
545

546
bool MeshSurfaceVisitor::Visit(const MeshFacet& face,
547
                               const MeshFacet&,
548
                               FacetIndex ulFInd,
549
                               unsigned long)
550
{
551
    indices.push_back(ulFInd);
552
    segm.AddFacet(face);
553
    return true;
554
}
555

556
// --------------------------------------------------------
557

558
void MeshSegmentAlgorithm::FindSegments(std::vector<MeshSurfaceSegmentPtr>& segm)
559
{
560
    // reset VISIT flags
561
    FacetIndex startFacet {};
562
    MeshCore::MeshAlgorithm cAlgo(myKernel);
563
    cAlgo.ResetFacetFlag(MeshCore::MeshFacet::VISIT);
564

565
    const MeshCore::MeshFacetArray& rFAry = myKernel.GetFacets();
566
    MeshCore::MeshFacetArray::_TConstIterator iCur = rFAry.begin();
567
    MeshCore::MeshFacetArray::_TConstIterator iBeg = rFAry.begin();
568
    MeshCore::MeshFacetArray::_TConstIterator iEnd = rFAry.end();
569

570
    // start from the first not visited facet
571
    cAlgo.CountFacetFlag(MeshCore::MeshFacet::VISIT);
572
    std::vector<FacetIndex> resetVisited;
573

574
    for (auto& it : segm) {
575
        cAlgo.ResetFacetsFlag(resetVisited, MeshCore::MeshFacet::VISIT);
576
        resetVisited.clear();
577

578
        MeshCore::MeshIsNotFlag<MeshCore::MeshFacet> flag;
579
        iCur = std::find_if(iBeg, iEnd, [flag](const MeshFacet& f) {
580
            return flag(f, MeshFacet::VISIT);
581
        });
582
        if (iCur < iEnd) {
583
            startFacet = iCur - iBeg;
584
        }
585
        else {
586
            startFacet = FACET_INDEX_MAX;
587
        }
588
        while (startFacet != FACET_INDEX_MAX) {
589
            // collect all facets of the same geometry
590
            std::vector<FacetIndex> indices;
591
            it->Initialize(startFacet);
592
            if (it->TestInitialFacet(startFacet)) {
593
                indices.push_back(startFacet);
594
            }
595
            MeshSurfaceVisitor pv(*it, indices);
596
            myKernel.VisitNeighbourFacets(pv, startFacet);
597

598
            // add or discard the segment
599
            if (indices.size() <= 1) {
600
                resetVisited.push_back(startFacet);
601
            }
602
            else {
603
                it->AddSegment(indices);
604
            }
605

606
            // search for the next start facet
607
            iCur = std::find_if(iCur, iEnd, [flag](const MeshFacet& f) {
608
                return flag(f, MeshFacet::VISIT);
609
            });
610
            if (iCur < iEnd) {
611
                startFacet = iCur - iBeg;
612
            }
613
            else {
614
                startFacet = FACET_INDEX_MAX;
615
            }
616
        }
617
    }
618
}
619