FreeCAD

1
/***************************************************************************
2
 *   Copyright (c) 2011 Werner Mayer <wmayer[at]users.sourceforge.net>     *
3
 *                                                                         *
4
 *   This file is part of the FreeCAD CAx development system.              *
5
 *                                                                         *
6
 *   This library is free software; you can redistribute it and/or         *
7
 *   modify it under the terms of the GNU Library General Public           *
8
 *   License as published by the Free Software Foundation; either          *
9
 *   version 2 of the License, or (at your option) any later version.      *
10
 *                                                                         *
11
 *   This library  is distributed in the hope that it will be useful,      *
12
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
13
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
14
 *   GNU Library General Public License for more details.                  *
15
 *                                                                         *
16
 *   You should have received a copy of the GNU Library General Public     *
17
 *   License along with this library; see the file COPYING.LIB. If not,    *
18
 *   write to the Free Software Foundation, Inc., 59 Temple Place,         *
19
 *   Suite 330, Boston, MA  02111-1307, USA                                *
20
 *                                                                         *
21
 ***************************************************************************/
22

23
#include "PreCompiled.h"
24

25
#ifndef _PreComp_
26
#include <boost/core/ignore_unused.hpp>
27
#include <numeric>
28

29
#include <BRepBuilderAPI_MakeVertex.hxx>
30
#include <BRepClass3d_SolidClassifier.hxx>
31
#include <BRepExtrema_DistShapeShape.hxx>
32
#include <BRepGProp_Face.hxx>
33
#include <TopExp.hxx>
34
#include <TopoDS.hxx>
35
#include <gp_Pnt.hxx>
36

37
#include <QEventLoop>
38
#include <QFuture>
39
#include <QFutureWatcher>
40
#include <QtConcurrentMap>
41
#endif
42

43
#include <Base/Console.h>
44
#include <Base/FutureWatcherProgress.h>
45
#include <Base/Sequencer.h>
46
#include <Base/Stream.h>
47

48
#include <Mod/Mesh/App/Core/Algorithm.h>
49
#include <Mod/Mesh/App/Core/Grid.h>
50
#include <Mod/Mesh/App/Core/Iterator.h>
51
#include <Mod/Mesh/App/Core/MeshKernel.h>
52
#include <Mod/Mesh/App/MeshFeature.h>
53
#include <Mod/Part/App/PartFeature.h>
54
#include <Mod/Points/App/PointsFeature.h>
55
#include <Mod/Points/App/PointsGrid.h>
56

57
#include "InspectionFeature.h"
58

59

60
using namespace Inspection;
61
namespace sp = std::placeholders;
62

63
InspectActualMesh::InspectActualMesh(const Mesh::MeshObject& rMesh)
64
    : _mesh(rMesh.getKernel())
65
{
66
    Base::Matrix4D tmp;
67
    _clTrf = rMesh.getTransform();
68
    _bApply = _clTrf != tmp;
69
}
70

71
InspectActualMesh::~InspectActualMesh() = default;
72

73
unsigned long InspectActualMesh::countPoints() const
74
{
75
    return _mesh.CountPoints();
76
}
77

78
Base::Vector3f InspectActualMesh::getPoint(unsigned long index) const
79
{
80
    Base::Vector3f point = _mesh.GetPoint(index);
81
    if (_bApply) {
82
        _clTrf.multVec(point, point);
83
    }
84
    return point;
85
}
86

87
// ----------------------------------------------------------------
88

89
InspectActualPoints::InspectActualPoints(const Points::PointKernel& rPoints)
90
    : _rKernel(rPoints)
91
{}
92

93
unsigned long InspectActualPoints::countPoints() const
94
{
95
    return _rKernel.size();
96
}
97

98
Base::Vector3f InspectActualPoints::getPoint(unsigned long index) const
99
{
100
    Base::Vector3d pnt = _rKernel.getPoint(index);
101
    return Base::Vector3f(float(pnt.x), float(pnt.y), float(pnt.z));
102
}
103

104
// ----------------------------------------------------------------
105

106
InspectActualShape::InspectActualShape(const Part::TopoShape& shape)
107
    : _rShape(shape)
108
{
109
    Standard_Real deflection = _rShape.getAccuracy();
110
    fetchPoints(deflection);
111
}
112

113
void InspectActualShape::fetchPoints(double deflection)
114
{
115
    // get points from faces or sub-sampled edges
116
    TopTools_IndexedMapOfShape mapOfShapes;
117
    TopExp::MapShapes(_rShape.getShape(), TopAbs_FACE, mapOfShapes);
118
    if (!mapOfShapes.IsEmpty()) {
119
        std::vector<Data::ComplexGeoData::Facet> faces;
120
        _rShape.getFaces(points, faces, deflection);
121
    }
122
    else {
123
        TopExp::MapShapes(_rShape.getShape(), TopAbs_EDGE, mapOfShapes);
124
        if (!mapOfShapes.IsEmpty()) {
125
            std::vector<Data::ComplexGeoData::Line> lines;
126
            _rShape.getLines(points, lines, deflection);
127
        }
128
        else {
129
            std::vector<Base::Vector3d> normals;
130
            _rShape.getPoints(points, normals, deflection);
131
        }
132
    }
133
}
134

135
unsigned long InspectActualShape::countPoints() const
136
{
137
    return points.size();
138
}
139

140
Base::Vector3f InspectActualShape::getPoint(unsigned long index) const
141
{
142
    return Base::toVector<float>(points[index]);
143
}
144

145
// ----------------------------------------------------------------
146

147
namespace Inspection
148
{
149
class MeshInspectGrid: public MeshCore::MeshGrid
150
{
151
public:
152
    MeshInspectGrid(const MeshCore::MeshKernel& mesh, float fGridLen, const Base::Matrix4D& mat)
153
        : MeshCore::MeshGrid(mesh)
154
        , _transform(mat)
155
    {
156
        Base::BoundBox3f clBBMesh = _pclMesh->GetBoundBox().Transformed(mat);
157
        Rebuild(std::max<unsigned long>((unsigned long)(clBBMesh.LengthX() / fGridLen), 1),
158
                std::max<unsigned long>((unsigned long)(clBBMesh.LengthY() / fGridLen), 1),
159
                std::max<unsigned long>((unsigned long)(clBBMesh.LengthZ() / fGridLen), 1));
160
    }
161

162
    void Validate(const MeshCore::MeshKernel& kernel) override
163
    {
164
        // do nothing
165
        boost::ignore_unused(kernel);
166
    }
167

168
    void Validate()
169
    {
170
        // do nothing
171
    }
172

173
    bool Verify() const override
174
    {
175
        // do nothing
176
        return true;
177
    }
178

179
protected:
180
    void CalculateGridLength(unsigned long /*ulCtGrid*/, unsigned long /*ulMaxGrids*/) override
181
    {
182
        // do nothing
183
    }
184

185
    void CalculateGridLength(int /*iCtGridPerAxis*/) override
186
    {
187
        // do nothing
188
    }
189

190
    unsigned long HasElements() const override
191
    {
192
        return _pclMesh->CountFacets();
193
    }
194

195
    void Pos(const Base::Vector3f& rclPoint,
196
             unsigned long& rulX,
197
             unsigned long& rulY,
198
             unsigned long& rulZ) const
199
    {
200
        rulX = (unsigned long)((rclPoint.x - _fMinX) / _fGridLenX);
201
        rulY = (unsigned long)((rclPoint.y - _fMinY) / _fGridLenY);
202
        rulZ = (unsigned long)((rclPoint.z - _fMinZ) / _fGridLenZ);
203

204
        assert((rulX < _ulCtGridsX) && (rulY < _ulCtGridsY) && (rulZ < _ulCtGridsZ));
205
    }
206

207
    void AddFacet(const MeshCore::MeshGeomFacet& rclFacet, unsigned long ulFacetIndex)
208
    {
209
        unsigned long ulX1;
210
        unsigned long ulY1;
211
        unsigned long ulZ1;
212
        unsigned long ulX2;
213
        unsigned long ulY2;
214
        unsigned long ulZ2;
215

216
        Base::BoundBox3f clBB;
217
        clBB.Add(rclFacet._aclPoints[0]);
218
        clBB.Add(rclFacet._aclPoints[1]);
219
        clBB.Add(rclFacet._aclPoints[2]);
220

221
        Pos(Base::Vector3f(clBB.MinX, clBB.MinY, clBB.MinZ), ulX1, ulY1, ulZ1);
222
        Pos(Base::Vector3f(clBB.MaxX, clBB.MaxY, clBB.MaxZ), ulX2, ulY2, ulZ2);
223

224

225
        if ((ulX1 < ulX2) || (ulY1 < ulY2) || (ulZ1 < ulZ2)) {
226
            for (unsigned long ulX = ulX1; ulX <= ulX2; ulX++) {
227
                for (unsigned long ulY = ulY1; ulY <= ulY2; ulY++) {
228
                    for (unsigned long ulZ = ulZ1; ulZ <= ulZ2; ulZ++) {
229
                        if (rclFacet.IntersectBoundingBox(GetBoundBox(ulX, ulY, ulZ))) {
230
                            _aulGrid[ulX][ulY][ulZ].insert(ulFacetIndex);
231
                        }
232
                    }
233
                }
234
            }
235
        }
236
        else {
237
            _aulGrid[ulX1][ulY1][ulZ1].insert(ulFacetIndex);
238
        }
239
    }
240

241
    void InitGrid() override
242
    {
243
        unsigned long i, j;
244

245
        Base::BoundBox3f clBBMesh = _pclMesh->GetBoundBox().Transformed(_transform);
246

247
        float fLengthX = clBBMesh.LengthX();
248
        float fLengthY = clBBMesh.LengthY();
249
        float fLengthZ = clBBMesh.LengthZ();
250

251
        _fGridLenX = (1.0f + fLengthX) / float(_ulCtGridsX);
252
        _fMinX = clBBMesh.MinX - 0.5f;
253

254
        _fGridLenY = (1.0f + fLengthY) / float(_ulCtGridsY);
255
        _fMinY = clBBMesh.MinY - 0.5f;
256

257
        _fGridLenZ = (1.0f + fLengthZ) / float(_ulCtGridsZ);
258
        _fMinZ = clBBMesh.MinZ - 0.5f;
259

260
        _aulGrid.clear();
261
        _aulGrid.resize(_ulCtGridsX);
262
        for (i = 0; i < _ulCtGridsX; i++) {
263
            _aulGrid[i].resize(_ulCtGridsY);
264
            for (j = 0; j < _ulCtGridsY; j++) {
265
                _aulGrid[i][j].resize(_ulCtGridsZ);
266
            }
267
        }
268
    }
269

270
    void RebuildGrid() override
271
    {
272
        _ulCtElements = _pclMesh->CountFacets();
273
        InitGrid();
274

275
        unsigned long i = 0;
276
        MeshCore::MeshFacetIterator clFIter(*_pclMesh);
277
        clFIter.Transform(_transform);
278
        for (clFIter.Init(); clFIter.More(); clFIter.Next()) {
279
            AddFacet(*clFIter, i++);
280
        }
281
    }
282

283
private:
284
    Base::Matrix4D _transform;
285
};
286
}  // namespace Inspection
287

288
InspectNominalMesh::InspectNominalMesh(const Mesh::MeshObject& rMesh, float offset)
289
    : _mesh(rMesh.getKernel())
290
{
291
    Base::Matrix4D tmp;
292
    _clTrf = rMesh.getTransform();
293
    _bApply = _clTrf != tmp;
294

295
    // Max. limit of grid elements
296
    float fMaxGridElements = 8000000.0f;
297
    Base::BoundBox3f box = _mesh.GetBoundBox().Transformed(rMesh.getTransform());
298

299
    // estimate the minimum allowed grid length
300
    float fMinGridLen =
301
        (float)pow((box.LengthX() * box.LengthY() * box.LengthZ() / fMaxGridElements), 0.3333f);
302
    float fGridLen = 5.0f * MeshCore::MeshAlgorithm(_mesh).GetAverageEdgeLength();
303

304
    // We want to avoid to get too small grid elements otherwise building up the grid structure
305
    // would take too much time and memory. Having quite a dense grid speeds up more the following
306
    // algorithms extremely. Due to the issue above it's always a compromise between speed and
307
    // memory usage.
308
    fGridLen = std::max<float>(fMinGridLen, fGridLen);
309

310
    // build up grid structure to speed up algorithms
311
    _pGrid = new MeshInspectGrid(_mesh, fGridLen, rMesh.getTransform());
312
    _box = box;
313
    _box.Enlarge(offset);
314
}
315

316
InspectNominalMesh::~InspectNominalMesh()
317
{
318
    delete this->_pGrid;
319
}
320

321
float InspectNominalMesh::getDistance(const Base::Vector3f& point) const
322
{
323
    if (!_box.IsInBox(point)) {
324
        return FLT_MAX;  // must be inside bbox
325
    }
326

327
    std::vector<unsigned long> indices;
328
    //_pGrid->GetElements(point, indices);
329
    if (indices.empty()) {
330
        std::set<unsigned long> inds;
331
        _pGrid->MeshGrid::SearchNearestFromPoint(point, inds);
332
        indices.insert(indices.begin(), inds.begin(), inds.end());
333
    }
334

335
    float fMinDist = FLT_MAX;
336
    bool positive = true;
337
    for (unsigned long it : indices) {
338
        MeshCore::MeshGeomFacet geomFace = _mesh.GetFacet(it);
339
        if (_bApply) {
340
            geomFace.Transform(_clTrf);
341
        }
342

343
        float fDist = geomFace.DistanceToPoint(point);
344
        if (fabs(fDist) < fabs(fMinDist)) {
345
            fMinDist = fDist;
346
            positive = point.DistanceToPlane(geomFace._aclPoints[0], geomFace.GetNormal()) > 0;
347
        }
348
    }
349

350
    if (!positive) {
351
        fMinDist = -fMinDist;
352
    }
353
    return fMinDist;
354
}
355

356
// ----------------------------------------------------------------
357

358
InspectNominalFastMesh::InspectNominalFastMesh(const Mesh::MeshObject& rMesh, float offset)
359
    : _mesh(rMesh.getKernel())
360
{
361
    const MeshCore::MeshKernel& kernel = rMesh.getKernel();
362

363
    Base::Matrix4D tmp;
364
    _clTrf = rMesh.getTransform();
365
    _bApply = _clTrf != tmp;
366

367
    // Max. limit of grid elements
368
    float fMaxGridElements = 8000000.0f;
369
    Base::BoundBox3f box = kernel.GetBoundBox().Transformed(rMesh.getTransform());
370

371
    // estimate the minimum allowed grid length
372
    float fMinGridLen =
373
        (float)pow((box.LengthX() * box.LengthY() * box.LengthZ() / fMaxGridElements), 0.3333f);
374
    float fGridLen = 5.0f * MeshCore::MeshAlgorithm(kernel).GetAverageEdgeLength();
375

376
    // We want to avoid to get too small grid elements otherwise building up the grid structure
377
    // would take too much time and memory. Having quite a dense grid speeds up more the following
378
    // algorithms extremely. Due to the issue above it's always a compromise between speed and
379
    // memory usage.
380
    fGridLen = std::max<float>(fMinGridLen, fGridLen);
381

382
    // build up grid structure to speed up algorithms
383
    _pGrid = new MeshInspectGrid(kernel, fGridLen, rMesh.getTransform());
384
    _box = box;
385
    _box.Enlarge(offset);
386
    max_level = (unsigned long)(offset / fGridLen);
387
}
388

389
InspectNominalFastMesh::~InspectNominalFastMesh()
390
{
391
    delete this->_pGrid;
392
}
393

394
/**
395
 * This algorithm is not that exact as that from InspectNominalMesh but is by
396
 * factors faster and sufficient for many cases.
397
 */
398
float InspectNominalFastMesh::getDistance(const Base::Vector3f& point) const
399
{
400
    if (!_box.IsInBox(point)) {
401
        return FLT_MAX;  // must be inside bbox
402
    }
403

404
    std::set<unsigned long> indices;
405
#if 0  // a point in a neighbour grid can be nearer
406
    std::vector<unsigned long> elements;
407
    _pGrid->GetElements(point, elements);
408
    indices.insert(elements.begin(), elements.end());
409
#else
410
    unsigned long ulX, ulY, ulZ;
411
    _pGrid->Position(point, ulX, ulY, ulZ);
412
    unsigned long ulLevel = 0;
413
    while (indices.empty() && ulLevel <= max_level) {
414
        _pGrid->GetHull(ulX, ulY, ulZ, ulLevel++, indices);
415
    }
416
    if (indices.empty() || ulLevel == 1) {
417
        _pGrid->GetHull(ulX, ulY, ulZ, ulLevel, indices);
418
    }
419
#endif
420

421
    float fMinDist = FLT_MAX;
422
    bool positive = true;
423
    for (unsigned long it : indices) {
424
        MeshCore::MeshGeomFacet geomFace = _mesh.GetFacet(it);
425
        if (_bApply) {
426
            geomFace.Transform(_clTrf);
427
        }
428

429
        float fDist = geomFace.DistanceToPoint(point);
430
        if (fabs(fDist) < fabs(fMinDist)) {
431
            fMinDist = fDist;
432
            positive = point.DistanceToPlane(geomFace._aclPoints[0], geomFace.GetNormal()) > 0;
433
        }
434
    }
435

436
    if (!positive) {
437
        fMinDist = -fMinDist;
438
    }
439
    return fMinDist;
440
}
441

442
// ----------------------------------------------------------------
443

444
InspectNominalPoints::InspectNominalPoints(const Points::PointKernel& Kernel, float /*offset*/)
445
    : _rKernel(Kernel)
446
{
447
    int uGridPerAxis = 50;  // totally 125.000 grid elements
448
    this->_pGrid = new Points::PointsGrid(Kernel, uGridPerAxis);
449
}
450

451
InspectNominalPoints::~InspectNominalPoints()
452
{
453
    delete this->_pGrid;
454
}
455

456
float InspectNominalPoints::getDistance(const Base::Vector3f& point) const
457
{
458
    // TODO: Make faster
459
    std::set<unsigned long> indices;
460
    unsigned long x, y, z;
461
    Base::Vector3d pointd(point.x, point.y, point.z);
462
    _pGrid->Position(pointd, x, y, z);
463
    _pGrid->GetElements(x, y, z, indices);
464

465
    double fMinDist = DBL_MAX;
466
    for (unsigned long it : indices) {
467
        Base::Vector3d pt = _rKernel.getPoint(it);
468
        double fDist = Base::Distance(pointd, pt);
469
        if (fDist < fMinDist) {
470
            fMinDist = fDist;
471
        }
472
    }
473

474
    return (float)fMinDist;
475
}
476

477
// ----------------------------------------------------------------
478

479
InspectNominalShape::InspectNominalShape(const TopoDS_Shape& shape, float /*radius*/)
480
    : _rShape(shape)
481
{
482
    distss = new BRepExtrema_DistShapeShape();
483
    distss->LoadS1(_rShape);
484

485
    // When having a solid then use its shell because otherwise the distance
486
    // for inner points will always be zero
487
    if (!_rShape.IsNull() && _rShape.ShapeType() == TopAbs_SOLID) {
488
        TopExp_Explorer xp;
489
        xp.Init(_rShape, TopAbs_SHELL);
490
        if (xp.More()) {
491
            distss->LoadS1(xp.Current());
492
            isSolid = true;
493
        }
494
    }
495
    // distss->SetDeflection(radius);
496
}
497

498
InspectNominalShape::~InspectNominalShape()
499
{
500
    delete distss;
501
}
502

503
float InspectNominalShape::getDistance(const Base::Vector3f& point) const
504
{
505
    gp_Pnt pnt3d(point.x, point.y, point.z);
506
    BRepBuilderAPI_MakeVertex mkVert(pnt3d);
507
    distss->LoadS2(mkVert.Vertex());
508

509
    float fMinDist = FLT_MAX;
510
    if (distss->Perform() && distss->NbSolution() > 0) {
511
        fMinDist = (float)distss->Value();
512
        // the shape is a solid, check if the vertex is inside
513
        if (isSolid) {
514
            if (isInsideSolid(pnt3d)) {
515
                fMinDist = -fMinDist;
516
            }
517
        }
518
        else if (fMinDist > 0) {
519
            // check if the distance was computed from a face
520
            if (isBelowFace(pnt3d)) {
521
                fMinDist = -fMinDist;
522
            }
523
        }
524
    }
525
    return fMinDist;
526
}
527

528
bool InspectNominalShape::isInsideSolid(const gp_Pnt& pnt3d) const
529
{
530
    const Standard_Real tol = 0.001;
531
    BRepClass3d_SolidClassifier classifier(_rShape);
532
    classifier.Perform(pnt3d, tol);
533
    return (classifier.State() == TopAbs_IN);
534
}
535

536
bool InspectNominalShape::isBelowFace(const gp_Pnt& pnt3d) const
537
{
538
    // check if the distance was computed from a face
539
    for (Standard_Integer index = 1; index <= distss->NbSolution(); index++) {
540
        if (distss->SupportTypeShape1(index) == BRepExtrema_IsInFace) {
541
            TopoDS_Shape face = distss->SupportOnShape1(index);
542
            Standard_Real u, v;
543
            distss->ParOnFaceS1(index, u, v);
544
            // gp_Pnt pnt = distss->PointOnShape1(index);
545
            BRepGProp_Face props(TopoDS::Face(face));
546
            gp_Vec normal;
547
            gp_Pnt center;
548
            props.Normal(u, v, center, normal);
549
            gp_Vec dir(center, pnt3d);
550
            Standard_Real scalar = normal.Dot(dir);
551
            if (scalar < 0) {
552
                return true;
553
            }
554
            break;
555
        }
556
    }
557

558
    return false;
559
}
560

561
// ----------------------------------------------------------------
562

563
TYPESYSTEM_SOURCE(Inspection::PropertyDistanceList, App::PropertyLists)
564

565
PropertyDistanceList::PropertyDistanceList() = default;
566

567
PropertyDistanceList::~PropertyDistanceList() = default;
568

569
void PropertyDistanceList::setSize(int newSize)
570
{
571
    _lValueList.resize(newSize);
572
}
573

574
int PropertyDistanceList::getSize() const
575
{
576
    return static_cast<int>(_lValueList.size());
577
}
578

579
void PropertyDistanceList::setValue(float lValue)
580
{
581
    aboutToSetValue();
582
    _lValueList.resize(1);
583
    _lValueList[0] = lValue;
584
    hasSetValue();
585
}
586

587
void PropertyDistanceList::setValues(const std::vector<float>& values)
588
{
589
    aboutToSetValue();
590
    _lValueList = values;
591
    hasSetValue();
592
}
593

594
PyObject* PropertyDistanceList::getPyObject()
595
{
596
    PyObject* list = PyList_New(getSize());
597
    for (int i = 0; i < getSize(); i++) {
598
        PyList_SetItem(list, i, PyFloat_FromDouble(_lValueList[i]));
599
    }
600
    return list;
601
}
602

603
void PropertyDistanceList::setPyObject(PyObject* value)
604
{
605
    if (PyList_Check(value)) {
606
        Py_ssize_t nSize = PyList_Size(value);
607
        std::vector<float> values;
608
        values.resize(nSize);
609

610
        for (Py_ssize_t i = 0; i < nSize; ++i) {
611
            PyObject* item = PyList_GetItem(value, i);
612
            if (!PyFloat_Check(item)) {
613
                std::string error = std::string("type in list must be float, not ");
614
                error += item->ob_type->tp_name;
615
                throw Py::TypeError(error);
616
            }
617

618
            values[i] = (float)PyFloat_AsDouble(item);
619
        }
620

621
        setValues(values);
622
    }
623
    else if (PyFloat_Check(value)) {
624
        setValue((float)PyFloat_AsDouble(value));
625
    }
626
    else {
627
        std::string error = std::string("type must be float or list of float, not ");
628
        error += value->ob_type->tp_name;
629
        throw Py::TypeError(error);
630
    }
631
}
632

633
void PropertyDistanceList::Save(Base::Writer& writer) const
634
{
635
    if (writer.isForceXML()) {
636
        writer.Stream() << writer.ind() << "<FloatList count=\"" << getSize() << "\">" << std::endl;
637
        writer.incInd();
638
        for (int i = 0; i < getSize(); i++) {
639
            writer.Stream() << writer.ind() << "<F v=\"" << _lValueList[i] << "\"/>" << std::endl;
640
        }
641
        writer.decInd();
642
        writer.Stream() << writer.ind() << "</FloatList>" << std::endl;
643
    }
644
    else {
645
        writer.Stream() << writer.ind() << "<FloatList file=\"" << writer.addFile(getName(), this)
646
                        << "\"/>" << std::endl;
647
    }
648
}
649

650
void PropertyDistanceList::Restore(Base::XMLReader& reader)
651
{
652
    reader.readElement("FloatList");
653
    std::string file(reader.getAttribute("file"));
654

655
    if (!file.empty()) {
656
        // initiate a file read
657
        reader.addFile(file.c_str(), this);
658
    }
659
}
660

661
void PropertyDistanceList::SaveDocFile(Base::Writer& writer) const
662
{
663
    Base::OutputStream str(writer.Stream());
664
    uint32_t uCt = (uint32_t)getSize();
665
    str << uCt;
666
    for (float it : _lValueList) {
667
        str << it;
668
    }
669
}
670

671
void PropertyDistanceList::RestoreDocFile(Base::Reader& reader)
672
{
673
    Base::InputStream str(reader);
674
    uint32_t uCt = 0;
675
    str >> uCt;
676
    std::vector<float> values(uCt);
677
    for (float& it : values) {
678
        str >> it;
679
    }
680
    setValues(values);
681
}
682

683
App::Property* PropertyDistanceList::Copy() const
684
{
685
    PropertyDistanceList* p = new PropertyDistanceList();
686
    p->_lValueList = _lValueList;
687
    return p;
688
}
689

690
void PropertyDistanceList::Paste(const App::Property& from)
691
{
692
    aboutToSetValue();
693
    _lValueList = dynamic_cast<const PropertyDistanceList&>(from)._lValueList;
694
    hasSetValue();
695
}
696

697
unsigned int PropertyDistanceList::getMemSize() const
698
{
699
    return static_cast<unsigned int>(_lValueList.size() * sizeof(float));
700
}
701

702
// ----------------------------------------------------------------
703

704
namespace Inspection
705
{
706
// helper class to use Qt's concurrent framework
707
struct DistanceInspection
708
{
709

710
    DistanceInspection(float radius,
711
                       InspectActualGeometry* a,
712
                       std::vector<InspectNominalGeometry*> n)
713
        : radius(radius)
714
        , actual(a)
715
        , nominal(n)
716
    {}
717
    float mapped(unsigned long index) const
718
    {
719
        Base::Vector3f pnt = actual->getPoint(index);
720

721
        float fMinDist = FLT_MAX;
722
        for (auto it : nominal) {
723
            float fDist = it->getDistance(pnt);
724
            if (fabs(fDist) < fabs(fMinDist)) {
725
                fMinDist = fDist;
726
            }
727
        }
728

729
        if (fMinDist > this->radius) {
730
            fMinDist = FLT_MAX;
731
        }
732
        else if (-fMinDist > this->radius) {
733
            fMinDist = -FLT_MAX;
734
        }
735

736
        return fMinDist;
737
    }
738

739
    float radius;
740
    InspectActualGeometry* actual;
741
    std::vector<InspectNominalGeometry*> nominal;
742
};
743

744
// Helper internal class for QtConcurrent map operation. Holds sums-of-squares and counts for RMS
745
// calculation
746
class DistanceInspectionRMS
747
{
748
public:
749
    DistanceInspectionRMS() = default;
750
    DistanceInspectionRMS& operator+=(const DistanceInspectionRMS& rhs)
751
    {
752
        this->m_numv += rhs.m_numv;
753
        this->m_sumsq += rhs.m_sumsq;
754
        return *this;
755
    }
756
    double getRMS()
757
    {
758
        if (this->m_numv == 0) {
759
            return 0.0;
760
        }
761
        return sqrt(this->m_sumsq / (double)this->m_numv);
762
    }
763
    int m_numv {0};
764
    double m_sumsq {0.0};
765
};
766
}  // namespace Inspection
767

768
PROPERTY_SOURCE(Inspection::Feature, App::DocumentObject)
769

770
Feature::Feature()
771
{
772
    ADD_PROPERTY(SearchRadius, (0.05));
773
    ADD_PROPERTY(Thickness, (0.0));
774
    ADD_PROPERTY(Actual, (nullptr));
775
    ADD_PROPERTY(Nominals, (nullptr));
776
    ADD_PROPERTY(Distances, (0.0));
777
}
778

779
Feature::~Feature() = default;
780

781
short Feature::mustExecute() const
782
{
783
    if (SearchRadius.isTouched()) {
784
        return 1;
785
    }
786
    if (Thickness.isTouched()) {
787
        return 1;
788
    }
789
    if (Actual.isTouched()) {
790
        return 1;
791
    }
792
    if (Nominals.isTouched()) {
793
        return 1;
794
    }
795
    return 0;
796
}
797

798
App::DocumentObjectExecReturn* Feature::execute()
799
{
800
    bool useMultithreading = true;
801

802
    App::DocumentObject* pcActual = Actual.getValue();
803
    if (!pcActual) {
804
        throw Base::ValueError("No actual geometry to inspect specified");
805
    }
806

807
    InspectActualGeometry* actual = nullptr;
808
    if (pcActual->isDerivedFrom<Mesh::Feature>()) {
809
        Mesh::Feature* mesh = static_cast<Mesh::Feature*>(pcActual);
810
        actual = new InspectActualMesh(mesh->Mesh.getValue());
811
    }
812
    else if (pcActual->isDerivedFrom<Points::Feature>()) {
813
        Points::Feature* pts = static_cast<Points::Feature*>(pcActual);
814
        actual = new InspectActualPoints(pts->Points.getValue());
815
    }
816
    else if (pcActual->isDerivedFrom<Part::Feature>()) {
817
        useMultithreading = false;
818
        Part::Feature* part = static_cast<Part::Feature*>(pcActual);
819
        actual = new InspectActualShape(part->Shape.getShape());
820
    }
821
    else {
822
        throw Base::TypeError("Unknown geometric type");
823
    }
824

825
    // clang-format off
826
    // get a list of nominals
827
    std::vector<InspectNominalGeometry*> inspectNominal;
828
    const std::vector<App::DocumentObject*>& nominals = Nominals.getValues();
829
    for (auto it : nominals) {
830
        InspectNominalGeometry* nominal = nullptr;
831
        if (it->isDerivedFrom<Mesh::Feature>()) {
832
            Mesh::Feature* mesh = static_cast<Mesh::Feature*>(it);
833
            nominal = new InspectNominalMesh(mesh->Mesh.getValue(), this->SearchRadius.getValue());
834
        }
835
        else if (it->isDerivedFrom<Points::Feature>()) {
836
            Points::Feature* pts = static_cast<Points::Feature*>(it);
837
            nominal = new InspectNominalPoints(pts->Points.getValue(), this->SearchRadius.getValue());
838
        }
839
        else if (it->isDerivedFrom<Part::Feature>()) {
840
            useMultithreading = false;
841
            Part::Feature* part = static_cast<Part::Feature*>(it);
842
            nominal = new InspectNominalShape(part->Shape.getValue(), this->SearchRadius.getValue());
843
        }
844

845
        if (nominal) {
846
            inspectNominal.push_back(nominal);
847
        }
848
    }
849
    // clang-format on
850

851
#if 0
852
#if 1  // test with some huge data sets
853
    std::vector<unsigned long> index(actual->countPoints());
854
    std::generate(index.begin(), index.end(), Base::iotaGen<unsigned long>(0));
855
    DistanceInspection check(this->SearchRadius.getValue(), actual, inspectNominal);
856
    QFuture<float> future = QtConcurrent::mapped
857
        (index, std::bind(&DistanceInspection::mapped, &check, sp::_1));
858
    //future.waitForFinished(); // blocks the GUI
859
    Base::FutureWatcherProgress progress("Inspecting...", actual->countPoints());
860
    QFutureWatcher<float> watcher;
861
    QObject::connect(&watcher, &QFutureWatcher<float>::progressValueChanged,
862
                     &progress, &Base::FutureWatcherProgress::progressValueChanged);
863

864
    // keep it responsive during computation
865
    QEventLoop loop;
866
    QObject::connect(&watcher, &QFutureWatcher::finished, &loop, &QEventLoop::quit);
867
    watcher.setFuture(future);
868
    loop.exec();
869

870
    std::vector<float> vals;
871
    vals.insert(vals.end(), future.begin(), future.end());
872
#else
873
    DistanceInspection insp(this->SearchRadius.getValue(), actual, inspectNominal);
874
    unsigned long count = actual->countPoints();
875
    std::stringstream str;
876
    str << "Inspecting " << this->Label.getValue() << "...";
877
    Base::SequencerLauncher seq(str.str().c_str(), count);
878

879
    std::vector<float> vals(count);
880
    for (unsigned long index = 0; index < count; index++) {
881
        float fMinDist = insp.mapped(index);
882
        vals[index] = fMinDist;
883
        seq.next();
884
    }
885
#endif
886

887
    Distances.setValues(vals);
888

889
    float fRMS = 0;
890
    int countRMS = 0;
891
    for (std::vector<float>::iterator it = vals.begin(); it != vals.end(); ++it) {
892
        if (fabs(*it) < FLT_MAX) {
893
            fRMS += (*it) * (*it);
894
            countRMS++;
895
        }
896
    }
897

898
    if (countRMS > 0) {
899
        fRMS = fRMS / countRMS;
900
        fRMS = sqrt(fRMS);
901
    }
902

903
    Base::Console().Message("RMS value for '%s' with search radius [%.4f,%.4f] is: %.4f\n",
904
        this->Label.getValue(), -this->SearchRadius.getValue(), this->SearchRadius.getValue(), fRMS);
905
#else
906
    unsigned long count = actual->countPoints();
907
    std::vector<float> vals(count);
908
    std::function<DistanceInspectionRMS(int)> fMap = [&](unsigned int index) {
909
        DistanceInspectionRMS res;
910
        Base::Vector3f pnt = actual->getPoint(index);
911

912
        float fMinDist = FLT_MAX;
913
        for (auto it : inspectNominal) {
914
            float fDist = it->getDistance(pnt);
915
            if (fabs(fDist) < fabs(fMinDist)) {
916
                fMinDist = fDist;
917
            }
918
        }
919

920
        if (fMinDist > this->SearchRadius.getValue()) {
921
            fMinDist = FLT_MAX;
922
        }
923
        else if (-fMinDist > this->SearchRadius.getValue()) {
924
            fMinDist = -FLT_MAX;
925
        }
926
        else {
927
            res.m_sumsq += fMinDist * fMinDist;
928
            res.m_numv++;
929
        }
930

931
        vals[index] = fMinDist;
932
        return res;
933
    };
934

935
    DistanceInspectionRMS res;
936

937
    if (useMultithreading) {
938
        // Build vector of increasing indices
939
        std::vector<unsigned long> index(count);
940
        std::iota(index.begin(), index.end(), 0);
941
        // Perform map-reduce operation : compute distances and update sum of squares for RMS
942
        // computation
943
        QFuture<DistanceInspectionRMS> future =
944
            QtConcurrent::mappedReduced(index, fMap, &DistanceInspectionRMS::operator+=);
945
        // Setup progress bar
946
        Base::FutureWatcherProgress progress("Inspecting...", actual->countPoints());
947
        QFutureWatcher<DistanceInspectionRMS> watcher;
948
        QObject::connect(&watcher,
949
                         &QFutureWatcher<DistanceInspectionRMS>::progressValueChanged,
950
                         &progress,
951
                         &Base::FutureWatcherProgress::progressValueChanged);
952
        // Keep UI responsive during computation
953
        QEventLoop loop;
954
        QObject::connect(&watcher,
955
                         &QFutureWatcher<DistanceInspectionRMS>::finished,
956
                         &loop,
957
                         &QEventLoop::quit);
958
        watcher.setFuture(future);
959
        loop.exec();
960
        res = future.result();
961
    }
962
    else {
963
        // Single-threaded operation
964
        std::stringstream str;
965
        str << "Inspecting " << this->Label.getValue() << "...";
966
        Base::SequencerLauncher seq(str.str().c_str(), count);
967

968
        for (unsigned int i = 0; i < count; i++) {
969
            res += fMap(i);
970
        }
971
    }
972

973
    Base::Console().Message("RMS value for '%s' with search radius [%.4f,%.4f] is: %.4f\n",
974
                            this->Label.getValue(),
975
                            -this->SearchRadius.getValue(),
976
                            this->SearchRadius.getValue(),
977
                            res.getRMS());
978
    Distances.setValues(vals);
979
#endif
980

981
    delete actual;
982
    for (auto it : inspectNominal) {
983
        delete it;
984
    }
985

986
    return nullptr;
987
}
988

989
// ----------------------------------------------------------------
990

991
PROPERTY_SOURCE(Inspection::Group, App::DocumentObjectGroup)
992

993

994
Group::Group() = default;
995

996
Group::~Group() = default;
997