FreeCAD

1
/***************************************************************************
2
 *   Copyright (c) 2008 Jürgen Riegel <juergen.riegel@web.de>              *
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 <Geom_BSplineSurface.hxx>
26
#include <TColgp_Array1OfPnt.hxx>
27
#endif
28

29
#include <Base/Console.h>
30
#include <Base/Converter.h>
31
#include <Base/GeometryPyCXX.h>
32
#include <Base/Interpreter.h>
33
#include <Base/PyWrapParseTupleAndKeywords.h>
34
#include <Mod/Mesh/App/MeshPy.h>
35
#include <Mod/Part/App/BSplineSurfacePy.h>
36
#include <Mod/Points/App/PointsPy.h>
37
#if defined(HAVE_PCL_FILTERS)
38
#include <pcl/filters/passthrough.h>
39
#include <pcl/filters/voxel_grid.h>
40
#include <pcl/point_types.h>
41
#endif
42

43
#include "ApproxSurface.h"
44
#include "BSplineFitting.h"
45
#include "RegionGrowing.h"
46
#include "SampleConsensus.h"
47
#include "Segmentation.h"
48
#include "SurfaceTriangulation.h"
49

50
// clang-format off
51
/*
52
Dependency of pcl components:
53
common: none
54
features: common, kdtree, octree, search, (range_image)
55
filters: common, kdtree, octree, sample_consenus, search
56
geometry: common
57
io: common, octree
58
kdtree: common
59
keypoints: common, features, filters, kdtree, octree, search, (range_image)
60
octree: common
61
recognition: common, features, search
62
registration: common, features, kdtree, sample_consensus
63
sample_consensus: common
64
search: common, kdtree, octree
65
segmentation: common, kdtree, octree, sample_consensus, search
66
surface: common, kdtree, octree, search
67
*/
68

69
using namespace Reen;
70

71
namespace Reen {
72
class Module : public Py::ExtensionModule<Module>
73
{
74
public:
75
    Module() : Py::ExtensionModule<Module>("ReverseEngineering")
76
    {
77
        add_keyword_method("approxCurve", &Module::approxCurve, "Approximate curve");
78
        add_keyword_method("approxSurface",&Module::approxSurface,
79
            "approxSurface(Points, UDegree=3, VDegree=3, NbUPoles=6, NbVPoles=6,\n"
80
            "Smooth=True, Weight=0.1, Grad=1.0, Bend=0.0, Curv=0.0\n"
81
            "Iterations=5, Correction=True, PatchFactor=1.0, UVDirs=((ux, uy, uz), (vx, vy, vz)))\n\n"
82
            "Points: the input data (e.g. a point cloud or mesh)\n"
83
            "UDegree: the degree in u parametric direction\n"
84
            "VDegree: the degree in v parametric direction\n"
85
            "NbUPoles: the number of control points in u parametric direction\n"
86
            "NbVPoles: the number of control points in v parametric direction\n"
87
            "Smooth: use energy terms to create a smooth surface\n"
88
            "Weight: weight of the energy terms altogether\n"
89
            "Grad: weight of the gradient term\n"
90
            "Bend: weight of the bending energy term\n"
91
            "Curv: weight of the deviation of curvature term\n"
92
            "Iterations: number of iterations\n"
93
            "Correction: perform a parameter correction of each iteration step\n"
94
            "PatchFactor: create an extended surface\n"
95
            "UVDirs: set the u,v parameter directions as tuple of two vectors\n"
96
            "        If not set then they will be determined by computing a best-fit plane\n"
97
        );
98
#if defined(HAVE_PCL_SURFACE)
99
        add_keyword_method("triangulate",&Module::triangulate,
100
            "triangulate(PointKernel,searchRadius[,mu=2.5])."
101
        );
102
        add_keyword_method("poissonReconstruction",&Module::poissonReconstruction,
103
            "poissonReconstruction(PointKernel)."
104
        );
105
        add_keyword_method("viewTriangulation",&Module::viewTriangulation,
106
            "viewTriangulation(PointKernel, width, height)."
107
        );
108
        add_keyword_method("gridProjection",&Module::gridProjection,
109
            "gridProjection(PointKernel)."
110
        );
111
        add_keyword_method("marchingCubesRBF",&Module::marchingCubesRBF,
112
            "marchingCubesRBF(PointKernel)."
113
        );
114
        add_keyword_method("marchingCubesHoppe",&Module::marchingCubesHoppe,
115
            "marchingCubesHoppe(PointKernel)."
116
        );
117
#endif
118
#if defined(HAVE_PCL_OPENNURBS)
119
        add_keyword_method("fitBSpline",&Module::fitBSpline,
120
            "fitBSpline(PointKernel)."
121
        );
122
#endif
123
#if defined(HAVE_PCL_FILTERS)
124
        add_keyword_method("filterVoxelGrid",&Module::filterVoxelGrid,
125
            "filterVoxelGrid(dim)."
126
        );
127
        add_keyword_method("normalEstimation",&Module::normalEstimation,
128
            "normalEstimation(Points,[KSearch=0, SearchRadius=0]) -> Normals\n"
129
            "KSearch is an int and used to search the k-nearest neighbours in\n"
130
            "the k-d tree. Alternatively, SearchRadius (a float) can be used\n"
131
            "as spatial distance to determine the neighbours of a point\n"
132
            "Example:\n"
133
            "\n"
134
            "import ReverseEngineering as Reen\n"
135
            "pts=App.ActiveDocument.ActiveObject.Points\n"
136
            "nor=Reen.normalEstimation(pts,KSearch=5)\n"
137
            "\n"
138
            "f=App.ActiveDocument.addObject('Points::FeaturePython','Normals')\n"
139
            "f.addProperty('Points::PropertyNormalList','Normal')\n"
140
            "f.Points=pts\n"
141
            "f.Normal=nor\n"
142
            "f.ViewObject.Proxy=0\n"
143
            "f.ViewObject.DisplayMode=1\n"
144
        );
145
#endif
146
#if defined(HAVE_PCL_SEGMENTATION)
147
        add_keyword_method("regionGrowingSegmentation",&Module::regionGrowingSegmentation,
148
            "regionGrowingSegmentation()."
149
        );
150
        add_keyword_method("featureSegmentation",&Module::featureSegmentation,
151
            "featureSegmentation()."
152
        );
153
#endif
154
#if defined(HAVE_PCL_SAMPLE_CONSENSUS)
155
        add_keyword_method("sampleConsensus",&Module::sampleConsensus,
156
            "sampleConsensus()."
157
        );
158
#endif
159
        initialize("This module is the ReverseEngineering module."); // register with Python
160
    }
161

162
private:
163
    static std::vector<Base::Vector3d> getPoints(PyObject* pts, bool closed)
164
    {
165
        std::vector<Base::Vector3d> data;
166
        if (PyObject_TypeCheck(pts, &(Points::PointsPy::Type))) {
167
            std::vector<Base::Vector3d> normal;
168
            auto pypts = static_cast<Points::PointsPy*>(pts);
169
            Points::PointKernel* points = pypts->getPointKernelPtr();
170
            points->getPoints(data, normal, 0.0);
171
        }
172
        else {
173
            Py::Sequence l(pts);
174
            data.reserve(l.size());
175
            for (Py::Sequence::iterator it = l.begin(); it != l.end(); ++it) {
176
                Py::Tuple t(*it);
177
                data.emplace_back(
178
                    Py::Float(t.getItem(0)),
179
                    Py::Float(t.getItem(1)),
180
                    Py::Float(t.getItem(2))
181
                );
182
            }
183
        }
184

185
        if (closed) {
186
            if (!data.empty()) {
187
                data.push_back(data.front());
188
            }
189
        }
190

191
        return data;
192
    }
193

194
    static PyObject* approx1(const Py::Tuple& args, const Py::Dict& kwds)
195
    {
196
        PyObject* pts {};
197
        PyObject* closed = Py_False;
198
        int minDegree = 3;  // NOLINT
199
        int maxDegree = 8;  // NOLINT
200
        int cont = int(GeomAbs_C2);
201
        double tol3d = 1.0e-3;  // NOLINT
202

203
        static const std::array<const char *, 7> kwds_approx{"Points",
204
                                                             "Closed",
205
                                                             "MinDegree",
206
                                                             "MaxDegree",
207
                                                             "Continuity",
208
                                                             "Tolerance",
209
                                                             nullptr};
210
        if (!Base::Wrapped_ParseTupleAndKeywords(args.ptr(), kwds.ptr(), "O|O!iiid", kwds_approx,
211
                                                 &pts, &PyBool_Type, &closed, &minDegree,
212
                                                 &maxDegree, &cont, &tol3d)) {
213
            return nullptr;
214
        }
215

216
        std::vector<Base::Vector3d> data = getPoints(pts, Base::asBoolean(closed));
217

218
        Part::GeomBSplineCurve curve;
219
        curve.approximate(data, minDegree, maxDegree, GeomAbs_Shape(cont), tol3d);
220
        return curve.getPyObject();
221
    }
222

223
    static PyObject* approx2(const Py::Tuple& args, const Py::Dict& kwds)
224
    {
225
        PyObject* pts {};
226
        char* parType {};
227
        PyObject* closed = Py_False;
228
        int minDegree = 3;  // NOLINT
229
        int maxDegree = 8;  // NOLINT
230
        int cont = int(GeomAbs_C2);
231
        double tol3d = 1.0e-3;  // NOLINT
232

233
        static const std::array<const char *, 8> kwds_approx{"Points",
234
                                                             "ParametrizationType",
235
                                                             "Closed",
236
                                                             "MinDegree",
237
                                                             "MaxDegree",
238
                                                             "Continuity",
239
                                                             "Tolerance",
240
                                                             nullptr};
241
        if (!Base::Wrapped_ParseTupleAndKeywords(args.ptr(), kwds.ptr(), "Os|O!iiid", kwds_approx,
242
                                                 &pts, &parType, &PyBool_Type, &closed, &minDegree,
243
                                                 &maxDegree, &cont, &tol3d)) {
244
            return nullptr;
245
        }
246

247
        std::vector<Base::Vector3d> data = getPoints(pts, Base::asBoolean(closed));
248

249
        Approx_ParametrizationType pt {Approx_ChordLength};
250
        std::string pstr = parType;
251
        if (pstr == "Uniform") {
252
            pt = Approx_IsoParametric;
253
        }
254
        else if (pstr == "Centripetal") {
255
            pt = Approx_Centripetal;
256
        }
257

258
        Part::GeomBSplineCurve curve;
259
        curve.approximate(data, pt, minDegree, maxDegree, GeomAbs_Shape(cont), tol3d);
260
        return curve.getPyObject();
261
    }
262

263
    static PyObject* approx3(const Py::Tuple& args, const Py::Dict& kwds)
264
    {
265
        PyObject* pts {};
266
        double weight1 {};
267
        double weight2 {};
268
        double weight3 {};
269
        PyObject* closed = Py_False;
270
        int maxDegree = 8;  // NOLINT
271
        int cont = int(GeomAbs_C2);
272
        double tol3d = 1.0e-3;  // NOLINT
273

274
        static const std::array<const char *, 9> kwds_approx{"Points",
275
                                                             "Weight1",
276
                                                             "Weight2",
277
                                                             "Weight3",
278
                                                             "Closed",
279
                                                             "MaxDegree",
280
                                                             "Continuity",
281
                                                             "Tolerance",
282
                                                             nullptr};
283
        if (!Base::Wrapped_ParseTupleAndKeywords(args.ptr(), kwds.ptr(), "Oddd|O!iid", kwds_approx,
284
                                                 &pts, &weight1, &weight2, &weight3,
285
                                                 &PyBool_Type, &closed,
286
                                                 &maxDegree, &cont, &tol3d)) {
287
            return nullptr;
288
        }
289

290
        std::vector<Base::Vector3d> data = getPoints(pts, Base::asBoolean(closed));
291

292
        Part::GeomBSplineCurve curve;
293
        curve.approximate(data, weight1, weight2, weight3, maxDegree, GeomAbs_Shape(cont), tol3d);
294
        return curve.getPyObject();
295
    }
296

297
    Py::Object approxCurve(const Py::Tuple& args, const Py::Dict& kwds)
298
    {
299
        try {
300
            using approxFunc = std::function<PyObject*(const Py::Tuple& args, const Py::Dict& kwds)>;
301

302
            std::vector<approxFunc> funcs;
303
            funcs.emplace_back(approx3);
304
            funcs.emplace_back(approx2);
305
            funcs.emplace_back(approx1);
306

307
            for (const auto& func : funcs) {
308
                if (PyObject* py = func(args, kwds)) {
309
                    return Py::asObject(py);
310
                }
311

312
                PyErr_Clear();
313
            }
314

315
            throw Py::ValueError("Wrong arguments ReverseEngineering.approxCurve()");
316
        }
317
        catch (const Base::Exception& e) {
318
            std::string msg = e.what();
319
            if (msg.empty()) {
320
                msg = "ReverseEngineering.approxCurve() failed";
321
            }
322
            throw Py::RuntimeError(msg);
323
        }
324
    }
325

326
    Py::Object approxSurface(const Py::Tuple& args, const Py::Dict& kwds)
327
    {
328
        PyObject *o;
329
        PyObject *uvdirs = nullptr;
330
        // spline parameters
331
        int uDegree = 3;
332
        int vDegree = 3;
333
        int uPoles = 6;
334
        int vPoles = 6;
335
        // smoothing
336
        PyObject* smooth = Py_True;
337
        double weight = 0.1;
338
        double grad = 1.0;  //0.5
339
        double bend = 0.0; //0.2
340
        double curv = 0.0; //0.3
341
        // other parameters
342
        int iteration = 5;
343
        PyObject* correction = Py_True;
344
        double factor = 1.0;
345

346
        static const std::array<const char *, 15> kwds_approx{"Points", "UDegree", "VDegree", "NbUPoles", "NbVPoles",
347
                                                              "Smooth", "Weight", "Grad", "Bend", "Curv", "Iterations",
348
                                                              "Correction", "PatchFactor", "UVDirs", nullptr};
349
        if (!Base::Wrapped_ParseTupleAndKeywords(args.ptr(), kwds.ptr(), "O|iiiiO!ddddiO!dO!", kwds_approx,
350
                                                 &o, &uDegree, &vDegree, &uPoles, &vPoles,
351
                                                 &PyBool_Type, &smooth, &weight, &grad, &bend, &curv,
352
                                                 &iteration, &PyBool_Type, &correction, &factor,
353
                                                 &PyTuple_Type, &uvdirs)) {
354
            throw Py::Exception();
355
        }
356

357
        int uOrder = uDegree + 1;
358
        int vOrder = vDegree + 1;
359

360
        // error checking
361
        if (grad < 0.0 || grad > 1.0) {
362
            throw Py::ValueError("Value of Grad out of range [0,1]");
363
        }
364
        if (bend < 0.0 || bend > 1.0) {
365
            throw Py::ValueError("Value of Bend out of range [0,1]");
366
        }
367
        if (curv < 0.0 || curv > 1.0) {
368
            throw Py::ValueError("Value of Curv out of range [0,1]");
369
        }
370
        if (uDegree < 1 || uOrder > uPoles) {
371
            throw Py::ValueError("Value of uDegree out of range [1,NbUPoles-1]");
372
        }
373
        if (vDegree < 1 || vOrder > vPoles) {
374
            throw Py::ValueError("Value of vDegree out of range [1,NbVPoles-1]");
375
        }
376

377
        double sum = (grad + bend + curv);
378
        if (sum > 0)
379
            weight = weight / sum;
380

381
        try {
382
            std::vector<Base::Vector3f> pts;
383
            if (PyObject_TypeCheck(o, &(Points::PointsPy::Type))) {
384
                Points::PointsPy* pPoints = static_cast<Points::PointsPy*>(o);
385
                Points::PointKernel* points = pPoints->getPointKernelPtr();
386
                pts = points->getBasicPoints();
387
            }
388
            else if (PyObject_TypeCheck(o, &(Mesh::MeshPy::Type))) {
389
                const Mesh::MeshObject* mesh = static_cast<Mesh::MeshPy*>(o)->getMeshObjectPtr();
390
                const MeshCore::MeshPointArray& points = mesh->getKernel().GetPoints();
391
                pts.insert(pts.begin(), points.begin(), points.end());
392
            }
393
            else {
394
                Py::Sequence l(o);
395
                pts.reserve(l.size());
396
                for (Py::Sequence::iterator it = l.begin(); it != l.end(); ++it) {
397
                    Py::Tuple t(*it);
398
                    pts.emplace_back(
399
                        Py::Float(t.getItem(0)),
400
                        Py::Float(t.getItem(1)),
401
                        Py::Float(t.getItem(2))
402
                    );
403
                }
404
            }
405

406
            TColgp_Array1OfPnt clPoints(0, pts.size()-1);
407
            if (clPoints.Length() < uPoles * vPoles) {
408
                throw Py::ValueError("Too less data points for the specified number of poles");
409
            }
410

411
            int index=0;
412
            for (const auto & pt : pts) {
413
                clPoints(index++) = gp_Pnt(pt.x, pt.y, pt.z);
414
            }
415

416
            Reen::BSplineParameterCorrection pc(uOrder,vOrder,uPoles,vPoles);
417
            Handle(Geom_BSplineSurface) hSurf;
418

419
            if (uvdirs) {
420
                Py::Tuple t(uvdirs);
421
                Base::Vector3d u = Py::Vector(t.getItem(0)).toVector();
422
                Base::Vector3d v = Py::Vector(t.getItem(1)).toVector();
423
                pc.SetUV(u, v);
424
            }
425
            pc.EnableSmoothing(Base::asBoolean(smooth), weight, grad, bend, curv);
426
            hSurf = pc.CreateSurface(clPoints, iteration, Base::asBoolean(correction), factor);
427
            if (!hSurf.IsNull()) {
428
                return Py::asObject(new Part::BSplineSurfacePy(new Part::GeomBSplineSurface(hSurf)));
429
            }
430

431
            throw Py::RuntimeError("Computation of B-spline surface failed");
432
        }
433
        catch (const Py::Exception&) {
434
            // re-throw
435
            throw;
436
        }
437
        catch (Standard_Failure &e) {
438
            std::string str;
439
            Standard_CString msg = e.GetMessageString();
440
            str += typeid(e).name();
441
            str += " ";
442
            if (msg) {str += msg;}
443
            else     {str += "No OCCT Exception Message";}
444
            throw Py::RuntimeError(str);
445
        }
446
        catch (const Base::Exception &e) {
447
            throw Py::RuntimeError(e.what());
448
        }
449
        catch (...) {
450
            throw Py::RuntimeError("Unknown C++ exception");
451
        }
452
    }
453
#if defined(HAVE_PCL_SURFACE)
454
    /*
455
import ReverseEngineering as Reen
456
import Points
457
import Mesh
458
import random
459

460
r=random.Random()
461

462
p=Points.Points()
463
pts=[]
464
for i in range(21):
465
  for j in range(21):
466
    pts.append(App.Vector(i,j,r.gauss(5,0.05)))
467

468
p.addPoints(pts)
469
m=Reen.triangulate(Points=p,SearchRadius=2.2)
470
Mesh.show(m)
471
    */
472
    Py::Object triangulate(const Py::Tuple& args, const Py::Dict& kwds)
473
    {
474
        PyObject *pts;
475
        double searchRadius;
476
        PyObject *vec = 0;
477
        int ksearch=5;
478
        double mu=2.5;
479

480
        static const std::array<const char*,6> kwds_greedy {"Points", "SearchRadius", "Mu", "KSearch",
481
                                      "Normals", NULL};
482
        if (!Base::Wrapped_ParseTupleAndKeywords(args.ptr(), kwds.ptr(), "O!d|diO", kwds_greedy,
483
                                        &(Points::PointsPy::Type), &pts,
484
                                        &searchRadius, &mu, &ksearch, &vec))
485
            throw Py::Exception();
486

487
        Points::PointKernel* points = static_cast<Points::PointsPy*>(pts)->getPointKernelPtr();
488

489
        Mesh::MeshObject* mesh = new Mesh::MeshObject();
490
        SurfaceTriangulation tria(*points, *mesh);
491
        tria.setMu(mu);
492
        tria.setSearchRadius(searchRadius);
493
        if (vec) {
494
            Py::Sequence list(vec);
495
            std::vector<Base::Vector3f> normals;
496
            normals.reserve(list.size());
497
            for (Py::Sequence::iterator it = list.begin(); it != list.end(); ++it) {
498
                Base::Vector3d v = Py::Vector(*it).toVector();
499
                normals.push_back(Base::convertTo<Base::Vector3f>(v));
500
            }
501
            tria.perform(normals);
502
        }
503
        else {
504
            tria.perform(ksearch);
505
        }
506

507
        return Py::asObject(new Mesh::MeshPy(mesh));
508
    }
509
    Py::Object poissonReconstruction(const Py::Tuple& args, const Py::Dict& kwds)
510
    {
511
        PyObject *pts;
512
        PyObject *vec = 0;
513
        int ksearch=5;
514
        int octreeDepth=-1;
515
        int solverDivide=-1;
516
        double samplesPerNode=-1.0;
517

518
        static const std::array<const char*,7> kwds_poisson {"Points", "KSearch", "OctreeDepth", "SolverDivide",
519
                                      "SamplesPerNode", "Normals", NULL};
520
        if (!Base::Wrapped_ParseTupleAndKeywords(args.ptr(), kwds.ptr(), "O!|iiidO", kwds_poisson,
521
                                        &(Points::PointsPy::Type), &pts,
522
                                        &ksearch, &octreeDepth, &solverDivide, &samplesPerNode, &vec))
523
            throw Py::Exception();
524

525
        Points::PointKernel* points = static_cast<Points::PointsPy*>(pts)->getPointKernelPtr();
526

527
        Mesh::MeshObject* mesh = new Mesh::MeshObject();
528
        Reen::PoissonReconstruction poisson(*points, *mesh);
529
        poisson.setDepth(octreeDepth);
530
        poisson.setSolverDivide(solverDivide);
531
        poisson.setSamplesPerNode(samplesPerNode);
532
        if (vec) {
533
            Py::Sequence list(vec);
534
            std::vector<Base::Vector3f> normals;
535
            normals.reserve(list.size());
536
            for (Py::Sequence::iterator it = list.begin(); it != list.end(); ++it) {
537
                Base::Vector3d v = Py::Vector(*it).toVector();
538
                normals.push_back(Base::convertTo<Base::Vector3f>(v));
539
            }
540
            poisson.perform(normals);
541
        }
542
        else {
543
            poisson.perform(ksearch);
544
        }
545

546
        return Py::asObject(new Mesh::MeshPy(mesh));
547
    }
548
   /*
549
import ReverseEngineering as Reen
550
import Points
551
import Mesh
552
import random
553
import math
554
r=random.Random()
555
p=Points.Points()
556
pts=[]
557
for i in range(21):
558
  for j in range(21):
559
    pts.append(App.Vector(i,j,r.random()))
560
p.addPoints(pts)
561
m=Reen.viewTriangulation(p,21,21)
562
Mesh.show(m)
563
def boxmueller():
564
  r1,r2=random.random(),random.random()
565
  return math.sqrt(-2*math.log(r1))*math.cos(2*math.pi*r2)
566
p=Points.Points()
567
pts=[]
568
for i in range(21):
569
  for j in range(21):
570
    pts.append(App.Vector(i,j,r.gauss(5,0.05)))
571
p.addPoints(pts)
572
m=Reen.viewTriangulation(p,21,21)
573
Mesh.show(m)
574
    */
575
    Py::Object viewTriangulation(const Py::Tuple& args, const Py::Dict& kwds)
576
    {
577
        PyObject *pts;
578
        int width;
579
        int height;
580

581
        static const std::array<const char*,4> kwds_view {"Points", "Width", "Height", NULL};
582
        if (!Base::Wrapped_ParseTupleAndKeywords(args.ptr(), kwds.ptr(), "O!|ii", kwds_view,
583
                                        &(Points::PointsPy::Type), &pts,
584
                                        &width, &height))
585
            throw Py::Exception();
586

587
        Points::PointKernel* points = static_cast<Points::PointsPy*>(pts)->getPointKernelPtr();
588

589
        try {
590
            Mesh::MeshObject* mesh = new Mesh::MeshObject();
591
            ImageTriangulation view(width, height, *points, *mesh);
592
            view.perform();
593

594
            return Py::asObject(new Mesh::MeshPy(mesh));
595
        }
596
        catch (const Base::Exception& e) {
597
            throw Py::RuntimeError(e.what());
598
        }
599
    }
600
    Py::Object gridProjection(const Py::Tuple& args, const Py::Dict& kwds)
601
    {
602
        PyObject *pts;
603
        PyObject *vec = 0;
604
        int ksearch=5;
605

606
        static const std::array<const char*,4> kwds_greedy {"Points", "KSearch", "Normals", NULL};
607
        if (!Base::Wrapped_ParseTupleAndKeywords(args.ptr(), kwds.ptr(), "O!|iO", kwds_greedy,
608
                                        &(Points::PointsPy::Type), &pts,
609
                                        &ksearch, &vec))
610
            throw Py::Exception();
611

612
        Points::PointKernel* points = static_cast<Points::PointsPy*>(pts)->getPointKernelPtr();
613

614
        Mesh::MeshObject* mesh = new Mesh::MeshObject();
615
        GridReconstruction tria(*points, *mesh);
616
        if (vec) {
617
            Py::Sequence list(vec);
618
            std::vector<Base::Vector3f> normals;
619
            normals.reserve(list.size());
620
            for (Py::Sequence::iterator it = list.begin(); it != list.end(); ++it) {
621
                Base::Vector3d v = Py::Vector(*it).toVector();
622
                normals.push_back(Base::convertTo<Base::Vector3f>(v));
623
            }
624
            tria.perform(normals);
625
        }
626
        else {
627
            tria.perform(ksearch);
628
        }
629

630
        return Py::asObject(new Mesh::MeshPy(mesh));
631
    }
632
    Py::Object marchingCubesRBF(const Py::Tuple& args, const Py::Dict& kwds)
633
    {
634
        PyObject *pts;
635
        PyObject *vec = 0;
636
        int ksearch=5;
637

638
        static const std::array<const char*,4> kwds_greedy {"Points", "KSearch", "Normals", NULL};
639
        if (!Base::Wrapped_ParseTupleAndKeywords(args.ptr(), kwds.ptr(), "O!|iO", kwds_greedy,
640
                                        &(Points::PointsPy::Type), &pts,
641
                                        &ksearch, &vec))
642
            throw Py::Exception();
643

644
        Points::PointKernel* points = static_cast<Points::PointsPy*>(pts)->getPointKernelPtr();
645

646
        Mesh::MeshObject* mesh = new Mesh::MeshObject();
647
        MarchingCubesRBF tria(*points, *mesh);
648
        if (vec) {
649
            Py::Sequence list(vec);
650
            std::vector<Base::Vector3f> normals;
651
            normals.reserve(list.size());
652
            for (Py::Sequence::iterator it = list.begin(); it != list.end(); ++it) {
653
                Base::Vector3d v = Py::Vector(*it).toVector();
654
                normals.push_back(Base::convertTo<Base::Vector3f>(v));
655
            }
656
            tria.perform(normals);
657
        }
658
        else {
659
            tria.perform(ksearch);
660
        }
661

662
        return Py::asObject(new Mesh::MeshPy(mesh));
663
    }
664
/*
665
import ReverseEngineering as Reen
666
import Points
667
import Mesh
668
import random
669
r=random.Random()
670
p=Points.Points()
671
pts=[]
672
for i in range(21):
673
  for j in range(21):
674
    pts.append(App.Vector(i,j,r.gauss(5,0.05)))
675
p.addPoints(pts)
676
m=Reen.marchingCubesHoppe(Points=p)
677
Mesh.show(m)
678
    */
679
    Py::Object marchingCubesHoppe(const Py::Tuple& args, const Py::Dict& kwds)
680
    {
681
        PyObject *pts;
682
        PyObject *vec = 0;
683
        int ksearch=5;
684

685
        static const std::array<const char*,4> kwds_greedy {"Points", "KSearch", "Normals", NULL};
686
        if (!Base::Wrapped_ParseTupleAndKeywords(args.ptr(), kwds.ptr(), "O!|iO", kwds_greedy,
687
                                        &(Points::PointsPy::Type), &pts,
688
                                        &ksearch, &vec))
689
            throw Py::Exception();
690

691
        Points::PointKernel* points = static_cast<Points::PointsPy*>(pts)->getPointKernelPtr();
692

693
        Mesh::MeshObject* mesh = new Mesh::MeshObject();
694
        MarchingCubesHoppe tria(*points, *mesh);
695
        if (vec) {
696
            Py::Sequence list(vec);
697
            std::vector<Base::Vector3f> normals;
698
            normals.reserve(list.size());
699
            for (Py::Sequence::iterator it = list.begin(); it != list.end(); ++it) {
700
                Base::Vector3d v = Py::Vector(*it).toVector();
701
                normals.push_back(Base::convertTo<Base::Vector3f>(v));
702
            }
703
            tria.perform(normals);
704
        }
705
        else {
706
            tria.perform(ksearch);
707
        }
708

709
        return Py::asObject(new Mesh::MeshPy(mesh));
710
    }
711
#endif
712
#if defined(HAVE_PCL_OPENNURBS)
713
    Py::Object fitBSpline(const Py::Tuple& args, const Py::Dict& kwds)
714
    {
715
        PyObject *pts;
716
        int degree = 2;
717
        int refinement = 4;
718
        int iterations = 10;
719
        double interiorSmoothness = 0.2;
720
        double interiorWeight = 1.0;
721
        double boundarySmoothness = 0.2;
722
        double boundaryWeight = 0.0;
723

724
        static const std::array<const char*,9> kwds_approx {"Points", "Degree", "Refinement", "Iterations",
725
                                      "InteriorSmoothness", "InteriorWeight", "BoundarySmoothness", "BoundaryWeight", NULL};
726
        if (!Base::Wrapped_ParseTupleAndKeywords(args.ptr(), kwds.ptr(), "O!|iiidddd", kwds_approx,
727
                                        &(Points::PointsPy::Type), &pts,
728
                                        &degree, &refinement, &iterations,
729
                                        &interiorSmoothness, &interiorWeight,
730
                                        &boundarySmoothness, &boundaryWeight))
731
            throw Py::Exception();
732

733
        Points::PointKernel* points = static_cast<Points::PointsPy*>(pts)->getPointKernelPtr();
734

735
        BSplineFitting fit(points->getBasicPoints());
736
        fit.setOrder(degree+1);
737
        fit.setRefinement(refinement);
738
        fit.setIterations(iterations);
739
        fit.setInteriorSmoothness(interiorSmoothness);
740
        fit.setInteriorWeight(interiorWeight);
741
        fit.setBoundarySmoothness(boundarySmoothness);
742
        fit.setBoundaryWeight(boundaryWeight);
743
        Handle(Geom_BSplineSurface) hSurf = fit.perform();
744

745
        if (!hSurf.IsNull()) {
746
            return Py::asObject(new Part::BSplineSurfacePy(new Part::GeomBSplineSurface(hSurf)));
747
        }
748

749
        throw Py::RuntimeError("Computation of B-spline surface failed");
750
    }
751
#endif
752
#if defined(HAVE_PCL_FILTERS)
753
    Py::Object filterVoxelGrid(const Py::Tuple& args, const Py::Dict& kwds)
754
    {
755
        PyObject *pts;
756
        double voxDimX = 0;
757
        double voxDimY = 0;
758
        double voxDimZ = 0;
759

760
        static const std::array<const char*,5>  kwds_voxel {"Points", "DimX", "DimY", "DimZ", NULL};
761
        if (!Base::Wrapped_ParseTupleAndKeywords(args.ptr(), kwds.ptr(), "O!d|dd", kwds_voxel,
762
                                        &(Points::PointsPy::Type), &pts,
763
                                        &voxDimX, &voxDimY, &voxDimZ))
764
            throw Py::Exception();
765

766
        if (voxDimY == 0)
767
            voxDimY = voxDimX;
768

769
        if (voxDimZ == 0)
770
            voxDimZ = voxDimX;
771

772
        Points::PointKernel* points = static_cast<Points::PointsPy*>(pts)->getPointKernelPtr();
773

774
        pcl::PointCloud<pcl::PointXYZ>::Ptr cloud (new pcl::PointCloud<pcl::PointXYZ>);
775
        cloud->reserve(points->size());
776
        for (Points::PointKernel::const_iterator it = points->begin(); it != points->end(); ++it) {
777
            cloud->push_back(pcl::PointXYZ(it->x, it->y, it->z));
778
        }
779

780
        // Create the filtering object
781
        pcl::PointCloud<pcl::PointXYZ>::Ptr cloud_downSmpl (new pcl::PointCloud<pcl::PointXYZ>);
782
        pcl::VoxelGrid<pcl::PointXYZ> voxG;
783
        voxG.setInputCloud (cloud);
784
        voxG.setLeafSize (voxDimX, voxDimY, voxDimZ);
785
        voxG.filter (*cloud_downSmpl);
786

787
        Points::PointKernel* points_sample = new Points::PointKernel();
788
        points_sample->reserve(cloud_downSmpl->size());
789
        for (pcl::PointCloud<pcl::PointXYZ>::const_iterator it = cloud_downSmpl->begin();it!=cloud_downSmpl->end();++it) {
790
            points_sample->push_back(Base::Vector3d(it->x,it->y,it->z));
791
        }
792

793
        return Py::asObject(new Points::PointsPy(points_sample));
794
    }
795
#endif
796
#if defined(HAVE_PCL_FILTERS)
797
    Py::Object normalEstimation(const Py::Tuple& args, const Py::Dict& kwds)
798
    {
799
        PyObject *pts;
800
        int ksearch=0;
801
        double searchRadius=0;
802

803
        static const std::array<const char*,4> kwds_normals {"Points", "KSearch", "SearchRadius", NULL};
804
        if (!Base::Wrapped_ParseTupleAndKeywords(args.ptr(), kwds.ptr(), "O!|id", kwds_normals,
805
                                        &(Points::PointsPy::Type), &pts,
806
                                        &ksearch, &searchRadius))
807
            throw Py::Exception();
808

809
        Points::PointKernel* points = static_cast<Points::PointsPy*>(pts)->getPointKernelPtr();
810

811
        std::vector<Base::Vector3d> normals;
812
        NormalEstimation estimate(*points);
813
        estimate.setKSearch(ksearch);
814
        estimate.setSearchRadius(searchRadius);
815
        estimate.perform(normals);
816

817
        Py::List list;
818
        for (std::vector<Base::Vector3d>::iterator it = normals.begin(); it != normals.end(); ++it) {
819
            list.append(Py::Vector(*it));
820
        }
821

822
        return list;
823
    }
824
#endif
825
#if defined(HAVE_PCL_SEGMENTATION)
826
    Py::Object regionGrowingSegmentation(const Py::Tuple& args, const Py::Dict& kwds)
827
    {
828
        PyObject *pts;
829
        PyObject *vec = 0;
830
        int ksearch=5;
831

832
        static const std::array<const char*,4> kwds_segment {"Points", "KSearch", "Normals", NULL};
833
        if (!Base::Wrapped_ParseTupleAndKeywords(args.ptr(), kwds.ptr(), "O!|iO", kwds_segment,
834
                                        &(Points::PointsPy::Type), &pts,
835
                                        &ksearch, &vec))
836
            throw Py::Exception();
837

838
        Points::PointKernel* points = static_cast<Points::PointsPy*>(pts)->getPointKernelPtr();
839

840
        std::list<std::vector<int> > clusters;
841
        RegionGrowing segm(*points, clusters);
842
        if (vec) {
843
            Py::Sequence list(vec);
844
            std::vector<Base::Vector3f> normals;
845
            normals.reserve(list.size());
846
            for (Py::Sequence::iterator it = list.begin(); it != list.end(); ++it) {
847
                Base::Vector3d v = Py::Vector(*it).toVector();
848
                normals.push_back(Base::convertTo<Base::Vector3f>(v));
849
            }
850
            segm.perform(normals);
851
        }
852
        else {
853
            segm.perform(ksearch);
854
        }
855

856
        Py::List lists;
857
        for (std::list<std::vector<int> >::iterator it = clusters.begin(); it != clusters.end(); ++it) {
858
            Py::Tuple tuple(it->size());
859
            for (std::size_t i = 0; i < it->size(); i++) {
860
                tuple.setItem(i, Py::Long((*it)[i]));
861
            }
862
            lists.append(tuple);
863
        }
864

865
        return lists;
866
    }
867
    Py::Object featureSegmentation(const Py::Tuple& args, const Py::Dict& kwds)
868
    {
869
        PyObject *pts;
870
        int ksearch=5;
871

872
        static const std::array<const char*,3> kwds_segment {"Points", "KSearch", NULL};
873
        if (!Base::Wrapped_ParseTupleAndKeywords(args.ptr(), kwds.ptr(), "O!|i", kwds_segment,
874
                                        &(Points::PointsPy::Type), &pts, &ksearch))
875
            throw Py::Exception();
876

877
        Points::PointKernel* points = static_cast<Points::PointsPy*>(pts)->getPointKernelPtr();
878

879
        std::list<std::vector<int> > clusters;
880
        Segmentation segm(*points, clusters);
881
        segm.perform(ksearch);
882

883
        Py::List lists;
884
        for (std::list<std::vector<int> >::iterator it = clusters.begin(); it != clusters.end(); ++it) {
885
            Py::Tuple tuple(it->size());
886
            for (std::size_t i = 0; i < it->size(); i++) {
887
                tuple.setItem(i, Py::Long((*it)[i]));
888
            }
889
            lists.append(tuple);
890
        }
891

892
        return lists;
893
    }
894
#endif
895
#if defined(HAVE_PCL_SAMPLE_CONSENSUS)
896
/*
897
import ReverseEngineering as reen
898
import Points
899
import Part
900
p = App.ActiveDocument.Points.Points
901
data = p.Points
902
n = reen.normalEstimation(p, 10)
903
model = reen.sampleConsensus(SacModel="Plane", Points=p)
904
indices = model["Model"]
905
param = model["Parameters"]
906
plane = Part.Plane()
907
plane.Axis = param[0:3]
908
plane.Position = -plane.Axis * param[3]
909
np = Points.Points()
910
np.addPoints([data[i] for i in indices])
911
Points.show(np)
912
# sort in descending order
913
indices = list(indices)
914
indices.sort(reverse=True)
915
# remove points of segment
916
for i in indices:
917
    del data[i]
918
    del n[i]
919
p = Points.Points()
920
p.addPoints(data)
921
model = reen.sampleConsensus(SacModel="Cylinder", Points=p, Normals=n)
922
indices = model["Model"]
923
np = Points.Points()
924
np.addPoints([data[i] for i in indices])
925
Points.show(np)
926
    */
927
    Py::Object sampleConsensus(const Py::Tuple& args, const Py::Dict& kwds)
928
    {
929
        PyObject *pts;
930
        PyObject *vec = nullptr;
931
        const char* sacModelType = nullptr;
932

933
        static const std::array<const char*,4> kwds_sample {"SacModel", "Points", "Normals", NULL};
934
        if (!Base::Wrapped_ParseTupleAndKeywords(args.ptr(), kwds.ptr(), "sO!|O", kwds_sample,
935
                                        &sacModelType, &(Points::PointsPy::Type), &pts, &vec))
936
            throw Py::Exception();
937

938
        Points::PointKernel* points = static_cast<Points::PointsPy*>(pts)->getPointKernelPtr();
939
        std::vector<Base::Vector3d> normals;
940
        if (vec) {
941
            Py::Sequence list(vec);
942
            normals.reserve(list.size());
943
            for (Py::Sequence::iterator it = list.begin(); it != list.end(); ++it) {
944
                Base::Vector3d v = Py::Vector(*it).toVector();
945
                normals.push_back(v);
946
            }
947
        }
948

949
        SampleConsensus::SacModel sacModel = SampleConsensus::SACMODEL_PLANE;
950
        if (sacModelType) {
951
            if (strcmp(sacModelType, "Cylinder") == 0)
952
                sacModel = SampleConsensus::SACMODEL_CYLINDER;
953
            else if (strcmp(sacModelType, "Sphere") == 0)
954
                sacModel = SampleConsensus::SACMODEL_SPHERE;
955
            else if (strcmp(sacModelType, "Cone") == 0)
956
                sacModel = SampleConsensus::SACMODEL_CONE;
957
        }
958

959
        std::vector<float> parameters;
960
        SampleConsensus sample(sacModel, *points, normals);
961
        std::vector<int> model;
962
        double probability = sample.perform(parameters, model);
963

964
        Py::Dict dict;
965
        Py::Tuple tuple(parameters.size());
966
        for (std::size_t i = 0; i < parameters.size(); i++)
967
            tuple.setItem(i, Py::Float(parameters[i]));
968
        Py::Tuple data(model.size());
969
        for (std::size_t i = 0; i < model.size(); i++)
970
            data.setItem(i, Py::Long(model[i]));
971
        dict.setItem(Py::String("Probability"), Py::Float(probability));
972
        dict.setItem(Py::String("Parameters"), tuple);
973
        dict.setItem(Py::String("Model"), data);
974

975
        return dict;
976
    }
977
#endif
978
};
979

980
PyObject* initModule()
981
{
982
    return Base::Interpreter().addModule(new Module);
983
}
984

985
} // namespace Reen
986

987

988
/* Python entry */
989
PyMOD_INIT_FUNC(ReverseEngineering)
990
{
991
    // load dependent module
992
    try {
993
        Base::Interpreter().loadModule("Part");
994
        Base::Interpreter().loadModule("Mesh");
995
    }
996
    catch(const Base::Exception& e) {
997
        PyErr_SetString(PyExc_ImportError, e.what());
998
        PyMOD_Return(nullptr);
999
    }
1000

1001
    PyObject* mod = Reen::initModule();
1002
    Base::Console().Log("Loading ReverseEngineering module... done\n");
1003
    PyMOD_Return(mod);
1004
}
1005
// clang-format on
1006