FreeCAD
991 строка · 31.0 Кб
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
60using namespace Inspection;
61namespace sp = std::placeholders;
62
63InspectActualMesh::InspectActualMesh(const Mesh::MeshObject& rMesh)
64: _mesh(rMesh.getKernel())
65{
66Base::Matrix4D tmp;
67_clTrf = rMesh.getTransform();
68_bApply = _clTrf != tmp;
69}
70
71InspectActualMesh::~InspectActualMesh() = default;
72
73unsigned long InspectActualMesh::countPoints() const
74{
75return _mesh.CountPoints();
76}
77
78Base::Vector3f InspectActualMesh::getPoint(unsigned long index) const
79{
80Base::Vector3f point = _mesh.GetPoint(index);
81if (_bApply) {
82_clTrf.multVec(point, point);
83}
84return point;
85}
86
87// ----------------------------------------------------------------
88
89InspectActualPoints::InspectActualPoints(const Points::PointKernel& rPoints)
90: _rKernel(rPoints)
91{}
92
93unsigned long InspectActualPoints::countPoints() const
94{
95return _rKernel.size();
96}
97
98Base::Vector3f InspectActualPoints::getPoint(unsigned long index) const
99{
100Base::Vector3d pnt = _rKernel.getPoint(index);
101return Base::Vector3f(float(pnt.x), float(pnt.y), float(pnt.z));
102}
103
104// ----------------------------------------------------------------
105
106InspectActualShape::InspectActualShape(const Part::TopoShape& shape)
107: _rShape(shape)
108{
109Standard_Real deflection = _rShape.getAccuracy();
110fetchPoints(deflection);
111}
112
113void InspectActualShape::fetchPoints(double deflection)
114{
115// get points from faces or sub-sampled edges
116TopTools_IndexedMapOfShape mapOfShapes;
117TopExp::MapShapes(_rShape.getShape(), TopAbs_FACE, mapOfShapes);
118if (!mapOfShapes.IsEmpty()) {
119std::vector<Data::ComplexGeoData::Facet> faces;
120_rShape.getFaces(points, faces, deflection);
121}
122else {
123TopExp::MapShapes(_rShape.getShape(), TopAbs_EDGE, mapOfShapes);
124if (!mapOfShapes.IsEmpty()) {
125std::vector<Data::ComplexGeoData::Line> lines;
126_rShape.getLines(points, lines, deflection);
127}
128else {
129std::vector<Base::Vector3d> normals;
130_rShape.getPoints(points, normals, deflection);
131}
132}
133}
134
135unsigned long InspectActualShape::countPoints() const
136{
137return points.size();
138}
139
140Base::Vector3f InspectActualShape::getPoint(unsigned long index) const
141{
142return Base::toVector<float>(points[index]);
143}
144
145// ----------------------------------------------------------------
146
147namespace Inspection
148{
149class MeshInspectGrid: public MeshCore::MeshGrid
150{
151public:
152MeshInspectGrid(const MeshCore::MeshKernel& mesh, float fGridLen, const Base::Matrix4D& mat)
153: MeshCore::MeshGrid(mesh)
154, _transform(mat)
155{
156Base::BoundBox3f clBBMesh = _pclMesh->GetBoundBox().Transformed(mat);
157Rebuild(std::max<unsigned long>((unsigned long)(clBBMesh.LengthX() / fGridLen), 1),
158std::max<unsigned long>((unsigned long)(clBBMesh.LengthY() / fGridLen), 1),
159std::max<unsigned long>((unsigned long)(clBBMesh.LengthZ() / fGridLen), 1));
160}
161
162void Validate(const MeshCore::MeshKernel& kernel) override
163{
164// do nothing
165boost::ignore_unused(kernel);
166}
167
168void Validate()
169{
170// do nothing
171}
172
173bool Verify() const override
174{
175// do nothing
176return true;
177}
178
179protected:
180void CalculateGridLength(int /*iCtGridPerAxis*/) override
181{
182// do nothing
183}
184
185unsigned long HasElements() const override
186{
187return _pclMesh->CountFacets();
188}
189
190void Pos(const Base::Vector3f& rclPoint,
191unsigned long& rulX,
192unsigned long& rulY,
193unsigned long& rulZ) const
194{
195rulX = (unsigned long)((rclPoint.x - _fMinX) / _fGridLenX);
196rulY = (unsigned long)((rclPoint.y - _fMinY) / _fGridLenY);
197rulZ = (unsigned long)((rclPoint.z - _fMinZ) / _fGridLenZ);
198
199assert((rulX < _ulCtGridsX) && (rulY < _ulCtGridsY) && (rulZ < _ulCtGridsZ));
200}
201
202void AddFacet(const MeshCore::MeshGeomFacet& rclFacet, unsigned long ulFacetIndex)
203{
204unsigned long ulX1;
205unsigned long ulY1;
206unsigned long ulZ1;
207unsigned long ulX2;
208unsigned long ulY2;
209unsigned long ulZ2;
210
211Base::BoundBox3f clBB;
212clBB.Add(rclFacet._aclPoints[0]);
213clBB.Add(rclFacet._aclPoints[1]);
214clBB.Add(rclFacet._aclPoints[2]);
215
216Pos(Base::Vector3f(clBB.MinX, clBB.MinY, clBB.MinZ), ulX1, ulY1, ulZ1);
217Pos(Base::Vector3f(clBB.MaxX, clBB.MaxY, clBB.MaxZ), ulX2, ulY2, ulZ2);
218
219
220if ((ulX1 < ulX2) || (ulY1 < ulY2) || (ulZ1 < ulZ2)) {
221for (unsigned long ulX = ulX1; ulX <= ulX2; ulX++) {
222for (unsigned long ulY = ulY1; ulY <= ulY2; ulY++) {
223for (unsigned long ulZ = ulZ1; ulZ <= ulZ2; ulZ++) {
224if (rclFacet.IntersectBoundingBox(GetBoundBox(ulX, ulY, ulZ))) {
225_aulGrid[ulX][ulY][ulZ].insert(ulFacetIndex);
226}
227}
228}
229}
230}
231else {
232_aulGrid[ulX1][ulY1][ulZ1].insert(ulFacetIndex);
233}
234}
235
236void InitGrid() override
237{
238unsigned long i, j;
239
240Base::BoundBox3f clBBMesh = _pclMesh->GetBoundBox().Transformed(_transform);
241
242float fLengthX = clBBMesh.LengthX();
243float fLengthY = clBBMesh.LengthY();
244float fLengthZ = clBBMesh.LengthZ();
245
246_fGridLenX = (1.0f + fLengthX) / float(_ulCtGridsX);
247_fMinX = clBBMesh.MinX - 0.5f;
248
249_fGridLenY = (1.0f + fLengthY) / float(_ulCtGridsY);
250_fMinY = clBBMesh.MinY - 0.5f;
251
252_fGridLenZ = (1.0f + fLengthZ) / float(_ulCtGridsZ);
253_fMinZ = clBBMesh.MinZ - 0.5f;
254
255_aulGrid.clear();
256_aulGrid.resize(_ulCtGridsX);
257for (i = 0; i < _ulCtGridsX; i++) {
258_aulGrid[i].resize(_ulCtGridsY);
259for (j = 0; j < _ulCtGridsY; j++) {
260_aulGrid[i][j].resize(_ulCtGridsZ);
261}
262}
263}
264
265void RebuildGrid() override
266{
267_ulCtElements = _pclMesh->CountFacets();
268InitGrid();
269
270unsigned long i = 0;
271MeshCore::MeshFacetIterator clFIter(*_pclMesh);
272clFIter.Transform(_transform);
273for (clFIter.Init(); clFIter.More(); clFIter.Next()) {
274AddFacet(*clFIter, i++);
275}
276}
277
278private:
279Base::Matrix4D _transform;
280};
281} // namespace Inspection
282
283InspectNominalMesh::InspectNominalMesh(const Mesh::MeshObject& rMesh, float offset)
284: _mesh(rMesh.getKernel())
285{
286Base::Matrix4D tmp;
287_clTrf = rMesh.getTransform();
288_bApply = _clTrf != tmp;
289
290// Max. limit of grid elements
291float fMaxGridElements = 8000000.0f;
292Base::BoundBox3f box = _mesh.GetBoundBox().Transformed(rMesh.getTransform());
293
294// estimate the minimum allowed grid length
295float fMinGridLen =
296(float)pow((box.LengthX() * box.LengthY() * box.LengthZ() / fMaxGridElements), 0.3333f);
297float fGridLen = 5.0f * MeshCore::MeshAlgorithm(_mesh).GetAverageEdgeLength();
298
299// We want to avoid to get too small grid elements otherwise building up the grid structure
300// would take too much time and memory. Having quite a dense grid speeds up more the following
301// algorithms extremely. Due to the issue above it's always a compromise between speed and
302// memory usage.
303fGridLen = std::max<float>(fMinGridLen, fGridLen);
304
305// build up grid structure to speed up algorithms
306_pGrid = new MeshInspectGrid(_mesh, fGridLen, rMesh.getTransform());
307_box = box;
308_box.Enlarge(offset);
309}
310
311InspectNominalMesh::~InspectNominalMesh()
312{
313delete this->_pGrid;
314}
315
316float InspectNominalMesh::getDistance(const Base::Vector3f& point) const
317{
318if (!_box.IsInBox(point)) {
319return FLT_MAX; // must be inside bbox
320}
321
322std::vector<unsigned long> indices;
323//_pGrid->GetElements(point, indices);
324if (indices.empty()) {
325std::set<unsigned long> inds;
326_pGrid->MeshGrid::SearchNearestFromPoint(point, inds);
327indices.insert(indices.begin(), inds.begin(), inds.end());
328}
329
330float fMinDist = FLT_MAX;
331bool positive = true;
332for (unsigned long it : indices) {
333MeshCore::MeshGeomFacet geomFace = _mesh.GetFacet(it);
334if (_bApply) {
335geomFace.Transform(_clTrf);
336}
337
338float fDist = geomFace.DistanceToPoint(point);
339if (fabs(fDist) < fabs(fMinDist)) {
340fMinDist = fDist;
341positive = point.DistanceToPlane(geomFace._aclPoints[0], geomFace.GetNormal()) > 0;
342}
343}
344
345if (!positive) {
346fMinDist = -fMinDist;
347}
348return fMinDist;
349}
350
351// ----------------------------------------------------------------
352
353InspectNominalFastMesh::InspectNominalFastMesh(const Mesh::MeshObject& rMesh, float offset)
354: _mesh(rMesh.getKernel())
355{
356const MeshCore::MeshKernel& kernel = rMesh.getKernel();
357
358Base::Matrix4D tmp;
359_clTrf = rMesh.getTransform();
360_bApply = _clTrf != tmp;
361
362// Max. limit of grid elements
363float fMaxGridElements = 8000000.0f;
364Base::BoundBox3f box = kernel.GetBoundBox().Transformed(rMesh.getTransform());
365
366// estimate the minimum allowed grid length
367float fMinGridLen =
368(float)pow((box.LengthX() * box.LengthY() * box.LengthZ() / fMaxGridElements), 0.3333f);
369float fGridLen = 5.0f * MeshCore::MeshAlgorithm(kernel).GetAverageEdgeLength();
370
371// We want to avoid to get too small grid elements otherwise building up the grid structure
372// would take too much time and memory. Having quite a dense grid speeds up more the following
373// algorithms extremely. Due to the issue above it's always a compromise between speed and
374// memory usage.
375fGridLen = std::max<float>(fMinGridLen, fGridLen);
376
377// build up grid structure to speed up algorithms
378_pGrid = new MeshInspectGrid(kernel, fGridLen, rMesh.getTransform());
379_box = box;
380_box.Enlarge(offset);
381max_level = (unsigned long)(offset / fGridLen);
382}
383
384InspectNominalFastMesh::~InspectNominalFastMesh()
385{
386delete this->_pGrid;
387}
388
389/**
390* This algorithm is not that exact as that from InspectNominalMesh but is by
391* factors faster and sufficient for many cases.
392*/
393float InspectNominalFastMesh::getDistance(const Base::Vector3f& point) const
394{
395if (!_box.IsInBox(point)) {
396return FLT_MAX; // must be inside bbox
397}
398
399std::set<unsigned long> indices;
400#if 0 // a point in a neighbour grid can be nearer
401std::vector<unsigned long> elements;
402_pGrid->GetElements(point, elements);
403indices.insert(elements.begin(), elements.end());
404#else
405unsigned long ulX, ulY, ulZ;
406_pGrid->Position(point, ulX, ulY, ulZ);
407unsigned long ulLevel = 0;
408while (indices.empty() && ulLevel <= max_level) {
409_pGrid->GetHull(ulX, ulY, ulZ, ulLevel++, indices);
410}
411if (indices.empty() || ulLevel == 1) {
412_pGrid->GetHull(ulX, ulY, ulZ, ulLevel, indices);
413}
414#endif
415
416float fMinDist = FLT_MAX;
417bool positive = true;
418for (unsigned long it : indices) {
419MeshCore::MeshGeomFacet geomFace = _mesh.GetFacet(it);
420if (_bApply) {
421geomFace.Transform(_clTrf);
422}
423
424float fDist = geomFace.DistanceToPoint(point);
425if (fabs(fDist) < fabs(fMinDist)) {
426fMinDist = fDist;
427positive = point.DistanceToPlane(geomFace._aclPoints[0], geomFace.GetNormal()) > 0;
428}
429}
430
431if (!positive) {
432fMinDist = -fMinDist;
433}
434return fMinDist;
435}
436
437// ----------------------------------------------------------------
438
439InspectNominalPoints::InspectNominalPoints(const Points::PointKernel& Kernel, float /*offset*/)
440: _rKernel(Kernel)
441{
442int uGridPerAxis = 50; // totally 125.000 grid elements
443this->_pGrid = new Points::PointsGrid(Kernel, uGridPerAxis);
444}
445
446InspectNominalPoints::~InspectNominalPoints()
447{
448delete this->_pGrid;
449}
450
451float InspectNominalPoints::getDistance(const Base::Vector3f& point) const
452{
453// TODO: Make faster
454std::set<unsigned long> indices;
455unsigned long x, y, z;
456Base::Vector3d pointd(point.x, point.y, point.z);
457_pGrid->Position(pointd, x, y, z);
458_pGrid->GetElements(x, y, z, indices);
459
460double fMinDist = DBL_MAX;
461for (unsigned long it : indices) {
462Base::Vector3d pt = _rKernel.getPoint(it);
463double fDist = Base::Distance(pointd, pt);
464if (fDist < fMinDist) {
465fMinDist = fDist;
466}
467}
468
469return (float)fMinDist;
470}
471
472// ----------------------------------------------------------------
473
474InspectNominalShape::InspectNominalShape(const TopoDS_Shape& shape, float /*radius*/)
475: _rShape(shape)
476{
477distss = new BRepExtrema_DistShapeShape();
478distss->LoadS1(_rShape);
479
480// When having a solid then use its shell because otherwise the distance
481// for inner points will always be zero
482if (!_rShape.IsNull() && _rShape.ShapeType() == TopAbs_SOLID) {
483TopExp_Explorer xp;
484xp.Init(_rShape, TopAbs_SHELL);
485if (xp.More()) {
486distss->LoadS1(xp.Current());
487isSolid = true;
488}
489}
490// distss->SetDeflection(radius);
491}
492
493InspectNominalShape::~InspectNominalShape()
494{
495delete distss;
496}
497
498float InspectNominalShape::getDistance(const Base::Vector3f& point) const
499{
500gp_Pnt pnt3d(point.x, point.y, point.z);
501BRepBuilderAPI_MakeVertex mkVert(pnt3d);
502distss->LoadS2(mkVert.Vertex());
503
504float fMinDist = FLT_MAX;
505if (distss->Perform() && distss->NbSolution() > 0) {
506fMinDist = (float)distss->Value();
507// the shape is a solid, check if the vertex is inside
508if (isSolid) {
509if (isInsideSolid(pnt3d)) {
510fMinDist = -fMinDist;
511}
512}
513else if (fMinDist > 0) {
514// check if the distance was computed from a face
515if (isBelowFace(pnt3d)) {
516fMinDist = -fMinDist;
517}
518}
519}
520return fMinDist;
521}
522
523bool InspectNominalShape::isInsideSolid(const gp_Pnt& pnt3d) const
524{
525const Standard_Real tol = 0.001;
526BRepClass3d_SolidClassifier classifier(_rShape);
527classifier.Perform(pnt3d, tol);
528return (classifier.State() == TopAbs_IN);
529}
530
531bool InspectNominalShape::isBelowFace(const gp_Pnt& pnt3d) const
532{
533// check if the distance was computed from a face
534for (Standard_Integer index = 1; index <= distss->NbSolution(); index++) {
535if (distss->SupportTypeShape1(index) == BRepExtrema_IsInFace) {
536TopoDS_Shape face = distss->SupportOnShape1(index);
537Standard_Real u, v;
538distss->ParOnFaceS1(index, u, v);
539// gp_Pnt pnt = distss->PointOnShape1(index);
540BRepGProp_Face props(TopoDS::Face(face));
541gp_Vec normal;
542gp_Pnt center;
543props.Normal(u, v, center, normal);
544gp_Vec dir(center, pnt3d);
545Standard_Real scalar = normal.Dot(dir);
546if (scalar < 0) {
547return true;
548}
549break;
550}
551}
552
553return false;
554}
555
556// ----------------------------------------------------------------
557
558TYPESYSTEM_SOURCE(Inspection::PropertyDistanceList, App::PropertyLists)
559
560PropertyDistanceList::PropertyDistanceList() = default;
561
562PropertyDistanceList::~PropertyDistanceList() = default;
563
564void PropertyDistanceList::setSize(int newSize)
565{
566_lValueList.resize(newSize);
567}
568
569int PropertyDistanceList::getSize() const
570{
571return static_cast<int>(_lValueList.size());
572}
573
574void PropertyDistanceList::setValue(float lValue)
575{
576aboutToSetValue();
577_lValueList.resize(1);
578_lValueList[0] = lValue;
579hasSetValue();
580}
581
582void PropertyDistanceList::setValues(const std::vector<float>& values)
583{
584aboutToSetValue();
585_lValueList = values;
586hasSetValue();
587}
588
589PyObject* PropertyDistanceList::getPyObject()
590{
591PyObject* list = PyList_New(getSize());
592for (int i = 0; i < getSize(); i++) {
593PyList_SetItem(list, i, PyFloat_FromDouble(_lValueList[i]));
594}
595return list;
596}
597
598void PropertyDistanceList::setPyObject(PyObject* value)
599{
600if (PyList_Check(value)) {
601Py_ssize_t nSize = PyList_Size(value);
602std::vector<float> values;
603values.resize(nSize);
604
605for (Py_ssize_t i = 0; i < nSize; ++i) {
606PyObject* item = PyList_GetItem(value, i);
607if (!PyFloat_Check(item)) {
608std::string error = std::string("type in list must be float, not ");
609error += item->ob_type->tp_name;
610throw Py::TypeError(error);
611}
612
613values[i] = (float)PyFloat_AsDouble(item);
614}
615
616setValues(values);
617}
618else if (PyFloat_Check(value)) {
619setValue((float)PyFloat_AsDouble(value));
620}
621else {
622std::string error = std::string("type must be float or list of float, not ");
623error += value->ob_type->tp_name;
624throw Py::TypeError(error);
625}
626}
627
628void PropertyDistanceList::Save(Base::Writer& writer) const
629{
630if (writer.isForceXML()) {
631writer.Stream() << writer.ind() << "<FloatList count=\"" << getSize() << "\">" << std::endl;
632writer.incInd();
633for (int i = 0; i < getSize(); i++) {
634writer.Stream() << writer.ind() << "<F v=\"" << _lValueList[i] << "\"/>" << std::endl;
635}
636writer.decInd();
637writer.Stream() << writer.ind() << "</FloatList>" << std::endl;
638}
639else {
640writer.Stream() << writer.ind() << "<FloatList file=\"" << writer.addFile(getName(), this)
641<< "\"/>" << std::endl;
642}
643}
644
645void PropertyDistanceList::Restore(Base::XMLReader& reader)
646{
647reader.readElement("FloatList");
648std::string file(reader.getAttribute("file"));
649
650if (!file.empty()) {
651// initiate a file read
652reader.addFile(file.c_str(), this);
653}
654}
655
656void PropertyDistanceList::SaveDocFile(Base::Writer& writer) const
657{
658Base::OutputStream str(writer.Stream());
659uint32_t uCt = (uint32_t)getSize();
660str << uCt;
661for (float it : _lValueList) {
662str << it;
663}
664}
665
666void PropertyDistanceList::RestoreDocFile(Base::Reader& reader)
667{
668Base::InputStream str(reader);
669uint32_t uCt = 0;
670str >> uCt;
671std::vector<float> values(uCt);
672for (float& it : values) {
673str >> it;
674}
675setValues(values);
676}
677
678App::Property* PropertyDistanceList::Copy() const
679{
680PropertyDistanceList* p = new PropertyDistanceList();
681p->_lValueList = _lValueList;
682return p;
683}
684
685void PropertyDistanceList::Paste(const App::Property& from)
686{
687aboutToSetValue();
688_lValueList = dynamic_cast<const PropertyDistanceList&>(from)._lValueList;
689hasSetValue();
690}
691
692unsigned int PropertyDistanceList::getMemSize() const
693{
694return static_cast<unsigned int>(_lValueList.size() * sizeof(float));
695}
696
697// ----------------------------------------------------------------
698
699namespace Inspection
700{
701// helper class to use Qt's concurrent framework
702struct DistanceInspection
703{
704
705DistanceInspection(float radius,
706InspectActualGeometry* a,
707std::vector<InspectNominalGeometry*> n)
708: radius(radius)
709, actual(a)
710, nominal(n)
711{}
712float mapped(unsigned long index) const
713{
714Base::Vector3f pnt = actual->getPoint(index);
715
716float fMinDist = FLT_MAX;
717for (auto it : nominal) {
718float fDist = it->getDistance(pnt);
719if (fabs(fDist) < fabs(fMinDist)) {
720fMinDist = fDist;
721}
722}
723
724if (fMinDist > this->radius) {
725fMinDist = FLT_MAX;
726}
727else if (-fMinDist > this->radius) {
728fMinDist = -FLT_MAX;
729}
730
731return fMinDist;
732}
733
734float radius;
735InspectActualGeometry* actual;
736std::vector<InspectNominalGeometry*> nominal;
737};
738
739// Helper internal class for QtConcurrent map operation. Holds sums-of-squares and counts for RMS
740// calculation
741class DistanceInspectionRMS
742{
743public:
744DistanceInspectionRMS() = default;
745DistanceInspectionRMS& operator+=(const DistanceInspectionRMS& rhs)
746{
747this->m_numv += rhs.m_numv;
748this->m_sumsq += rhs.m_sumsq;
749return *this;
750}
751double getRMS()
752{
753if (this->m_numv == 0) {
754return 0.0;
755}
756return sqrt(this->m_sumsq / (double)this->m_numv);
757}
758int m_numv {0};
759double m_sumsq {0.0};
760};
761} // namespace Inspection
762
763PROPERTY_SOURCE(Inspection::Feature, App::DocumentObject)
764
765Feature::Feature()
766{
767ADD_PROPERTY(SearchRadius, (0.05));
768ADD_PROPERTY(Thickness, (0.0));
769ADD_PROPERTY(Actual, (nullptr));
770ADD_PROPERTY(Nominals, (nullptr));
771ADD_PROPERTY(Distances, (0.0));
772}
773
774Feature::~Feature() = default;
775
776short Feature::mustExecute() const
777{
778if (SearchRadius.isTouched()) {
779return 1;
780}
781if (Thickness.isTouched()) {
782return 1;
783}
784if (Actual.isTouched()) {
785return 1;
786}
787if (Nominals.isTouched()) {
788return 1;
789}
790return 0;
791}
792
793App::DocumentObjectExecReturn* Feature::execute()
794{
795bool useMultithreading = true;
796
797App::DocumentObject* pcActual = Actual.getValue();
798if (!pcActual) {
799throw Base::ValueError("No actual geometry to inspect specified");
800}
801
802InspectActualGeometry* actual = nullptr;
803if (pcActual->isDerivedFrom<Mesh::Feature>()) {
804Mesh::Feature* mesh = static_cast<Mesh::Feature*>(pcActual);
805actual = new InspectActualMesh(mesh->Mesh.getValue());
806}
807else if (pcActual->isDerivedFrom<Points::Feature>()) {
808Points::Feature* pts = static_cast<Points::Feature*>(pcActual);
809actual = new InspectActualPoints(pts->Points.getValue());
810}
811else if (pcActual->isDerivedFrom<Part::Feature>()) {
812useMultithreading = false;
813Part::Feature* part = static_cast<Part::Feature*>(pcActual);
814actual = new InspectActualShape(part->Shape.getShape());
815}
816else {
817throw Base::TypeError("Unknown geometric type");
818}
819
820// clang-format off
821// get a list of nominals
822std::vector<InspectNominalGeometry*> inspectNominal;
823const std::vector<App::DocumentObject*>& nominals = Nominals.getValues();
824for (auto it : nominals) {
825InspectNominalGeometry* nominal = nullptr;
826if (it->isDerivedFrom<Mesh::Feature>()) {
827Mesh::Feature* mesh = static_cast<Mesh::Feature*>(it);
828nominal = new InspectNominalMesh(mesh->Mesh.getValue(), this->SearchRadius.getValue());
829}
830else if (it->isDerivedFrom<Points::Feature>()) {
831Points::Feature* pts = static_cast<Points::Feature*>(it);
832nominal = new InspectNominalPoints(pts->Points.getValue(), this->SearchRadius.getValue());
833}
834else if (it->isDerivedFrom<Part::Feature>()) {
835useMultithreading = false;
836Part::Feature* part = static_cast<Part::Feature*>(it);
837nominal = new InspectNominalShape(part->Shape.getValue(), this->SearchRadius.getValue());
838}
839
840if (nominal) {
841inspectNominal.push_back(nominal);
842}
843}
844// clang-format on
845
846#if 0
847#if 1 // test with some huge data sets
848std::vector<unsigned long> index(actual->countPoints());
849std::generate(index.begin(), index.end(), Base::iotaGen<unsigned long>(0));
850DistanceInspection check(this->SearchRadius.getValue(), actual, inspectNominal);
851QFuture<float> future = QtConcurrent::mapped
852(index, std::bind(&DistanceInspection::mapped, &check, sp::_1));
853//future.waitForFinished(); // blocks the GUI
854Base::FutureWatcherProgress progress("Inspecting...", actual->countPoints());
855QFutureWatcher<float> watcher;
856QObject::connect(&watcher, &QFutureWatcher<float>::progressValueChanged,
857&progress, &Base::FutureWatcherProgress::progressValueChanged);
858
859// keep it responsive during computation
860QEventLoop loop;
861QObject::connect(&watcher, &QFutureWatcher::finished, &loop, &QEventLoop::quit);
862watcher.setFuture(future);
863loop.exec();
864
865std::vector<float> vals;
866vals.insert(vals.end(), future.begin(), future.end());
867#else
868DistanceInspection insp(this->SearchRadius.getValue(), actual, inspectNominal);
869unsigned long count = actual->countPoints();
870std::stringstream str;
871str << "Inspecting " << this->Label.getValue() << "...";
872Base::SequencerLauncher seq(str.str().c_str(), count);
873
874std::vector<float> vals(count);
875for (unsigned long index = 0; index < count; index++) {
876float fMinDist = insp.mapped(index);
877vals[index] = fMinDist;
878seq.next();
879}
880#endif
881
882Distances.setValues(vals);
883
884float fRMS = 0;
885int countRMS = 0;
886for (std::vector<float>::iterator it = vals.begin(); it != vals.end(); ++it) {
887if (fabs(*it) < FLT_MAX) {
888fRMS += (*it) * (*it);
889countRMS++;
890}
891}
892
893if (countRMS > 0) {
894fRMS = fRMS / countRMS;
895fRMS = sqrt(fRMS);
896}
897
898Base::Console().Message("RMS value for '%s' with search radius [%.4f,%.4f] is: %.4f\n",
899this->Label.getValue(), -this->SearchRadius.getValue(), this->SearchRadius.getValue(), fRMS);
900#else
901unsigned long count = actual->countPoints();
902std::vector<float> vals(count);
903std::function<DistanceInspectionRMS(int)> fMap = [&](unsigned int index) {
904DistanceInspectionRMS res;
905Base::Vector3f pnt = actual->getPoint(index);
906
907float fMinDist = FLT_MAX;
908for (auto it : inspectNominal) {
909float fDist = it->getDistance(pnt);
910if (fabs(fDist) < fabs(fMinDist)) {
911fMinDist = fDist;
912}
913}
914
915if (fMinDist > this->SearchRadius.getValue()) {
916fMinDist = FLT_MAX;
917}
918else if (-fMinDist > this->SearchRadius.getValue()) {
919fMinDist = -FLT_MAX;
920}
921else {
922res.m_sumsq += fMinDist * fMinDist;
923res.m_numv++;
924}
925
926vals[index] = fMinDist;
927return res;
928};
929
930DistanceInspectionRMS res;
931
932if (useMultithreading) {
933// Build vector of increasing indices
934std::vector<unsigned long> index(count);
935std::iota(index.begin(), index.end(), 0);
936// Perform map-reduce operation : compute distances and update sum of squares for RMS
937// computation
938QFuture<DistanceInspectionRMS> future =
939QtConcurrent::mappedReduced(index, fMap, &DistanceInspectionRMS::operator+=);
940// Setup progress bar
941Base::FutureWatcherProgress progress("Inspecting...", actual->countPoints());
942QFutureWatcher<DistanceInspectionRMS> watcher;
943QObject::connect(&watcher,
944&QFutureWatcher<DistanceInspectionRMS>::progressValueChanged,
945&progress,
946&Base::FutureWatcherProgress::progressValueChanged);
947// Keep UI responsive during computation
948QEventLoop loop;
949QObject::connect(&watcher,
950&QFutureWatcher<DistanceInspectionRMS>::finished,
951&loop,
952&QEventLoop::quit);
953watcher.setFuture(future);
954loop.exec();
955res = future.result();
956}
957else {
958// Single-threaded operation
959std::stringstream str;
960str << "Inspecting " << this->Label.getValue() << "...";
961Base::SequencerLauncher seq(str.str().c_str(), count);
962
963for (unsigned int i = 0; i < count; i++) {
964res += fMap(i);
965}
966}
967
968Base::Console().Message("RMS value for '%s' with search radius [%.4f,%.4f] is: %.4f\n",
969this->Label.getValue(),
970-this->SearchRadius.getValue(),
971this->SearchRadius.getValue(),
972res.getRMS());
973Distances.setValues(vals);
974#endif
975
976delete actual;
977for (auto it : inspectNominal) {
978delete it;
979}
980
981return nullptr;
982}
983
984// ----------------------------------------------------------------
985
986PROPERTY_SOURCE(Inspection::Group, App::DocumentObjectGroup)
987
988
989Group::Group() = default;
990
991Group::~Group() = default;
992