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

161
private:
162
    Py::Object approxSurface(const Py::Tuple& args, const Py::Dict& kwds)
163
    {
164
        PyObject *o;
165
        PyObject *uvdirs = nullptr;
166
        // spline parameters
167
        int uDegree = 3;
168
        int vDegree = 3;
169
        int uPoles = 6;
170
        int vPoles = 6;
171
        // smoothing
172
        PyObject* smooth = Py_True;
173
        double weight = 0.1;
174
        double grad = 1.0;  //0.5
175
        double bend = 0.0; //0.2
176
        double curv = 0.0; //0.3
177
        // other parameters
178
        int iteration = 5;
179
        PyObject* correction = Py_True;
180
        double factor = 1.0;
181

182
        static const std::array<const char *, 15> kwds_approx{"Points", "UDegree", "VDegree", "NbUPoles", "NbVPoles",
183
                                                              "Smooth", "Weight", "Grad", "Bend", "Curv", "Iterations",
184
                                                              "Correction", "PatchFactor", "UVDirs", nullptr};
185
        if (!Base::Wrapped_ParseTupleAndKeywords(args.ptr(), kwds.ptr(), "O|iiiiO!ddddiO!dO!", kwds_approx,
186
                                                 &o, &uDegree, &vDegree, &uPoles, &vPoles,
187
                                                 &PyBool_Type, &smooth, &weight, &grad, &bend, &curv,
188
                                                 &iteration, &PyBool_Type, &correction, &factor,
189
                                                 &PyTuple_Type, &uvdirs)) {
190
            throw Py::Exception();
191
        }
192

193
        int uOrder = uDegree + 1;
194
        int vOrder = vDegree + 1;
195

196
        // error checking
197
        if (grad < 0.0 || grad > 1.0) {
198
            throw Py::ValueError("Value of Grad out of range [0,1]");
199
        }
200
        if (bend < 0.0 || bend > 1.0) {
201
            throw Py::ValueError("Value of Bend out of range [0,1]");
202
        }
203
        if (curv < 0.0 || curv > 1.0) {
204
            throw Py::ValueError("Value of Curv out of range [0,1]");
205
        }
206
        if (uDegree < 1 || uOrder > uPoles) {
207
            throw Py::ValueError("Value of uDegree out of range [1,NbUPoles-1]");
208
        }
209
        if (vDegree < 1 || vOrder > vPoles) {
210
            throw Py::ValueError("Value of vDegree out of range [1,NbVPoles-1]");
211
        }
212

213
        double sum = (grad + bend + curv);
214
        if (sum > 0)
215
            weight = weight / sum;
216

217
        try {
218
            std::vector<Base::Vector3f> pts;
219
            if (PyObject_TypeCheck(o, &(Points::PointsPy::Type))) {
220
                Points::PointsPy* pPoints = static_cast<Points::PointsPy*>(o);
221
                Points::PointKernel* points = pPoints->getPointKernelPtr();
222
                pts = points->getBasicPoints();
223
            }
224
            else if (PyObject_TypeCheck(o, &(Mesh::MeshPy::Type))) {
225
                const Mesh::MeshObject* mesh = static_cast<Mesh::MeshPy*>(o)->getMeshObjectPtr();
226
                const MeshCore::MeshPointArray& points = mesh->getKernel().GetPoints();
227
                pts.insert(pts.begin(), points.begin(), points.end());
228
            }
229
            else {
230
                Py::Sequence l(o);
231
                pts.reserve(l.size());
232
                for (Py::Sequence::iterator it = l.begin(); it != l.end(); ++it) {
233
                    Py::Tuple t(*it);
234
                    pts.emplace_back(
235
                        Py::Float(t.getItem(0)),
236
                        Py::Float(t.getItem(1)),
237
                        Py::Float(t.getItem(2))
238
                    );
239
                }
240
            }
241

242
            TColgp_Array1OfPnt clPoints(0, pts.size()-1);
243
            if (clPoints.Length() < uPoles * vPoles) {
244
                throw Py::ValueError("Too less data points for the specified number of poles");
245
            }
246

247
            int index=0;
248
            for (const auto & pt : pts) {
249
                clPoints(index++) = gp_Pnt(pt.x, pt.y, pt.z);
250
            }
251

252
            Reen::BSplineParameterCorrection pc(uOrder,vOrder,uPoles,vPoles);
253
            Handle(Geom_BSplineSurface) hSurf;
254

255
            if (uvdirs) {
256
                Py::Tuple t(uvdirs);
257
                Base::Vector3d u = Py::Vector(t.getItem(0)).toVector();
258
                Base::Vector3d v = Py::Vector(t.getItem(1)).toVector();
259
                pc.SetUV(u, v);
260
            }
261
            pc.EnableSmoothing(Base::asBoolean(smooth), weight, grad, bend, curv);
262
            hSurf = pc.CreateSurface(clPoints, iteration, Base::asBoolean(correction), factor);
263
            if (!hSurf.IsNull()) {
264
                return Py::asObject(new Part::BSplineSurfacePy(new Part::GeomBSplineSurface(hSurf)));
265
            }
266

267
            throw Py::RuntimeError("Computation of B-spline surface failed");
268
        }
269
        catch (const Py::Exception&) {
270
            // re-throw
271
            throw;
272
        }
273
        catch (Standard_Failure &e) {
274
            std::string str;
275
            Standard_CString msg = e.GetMessageString();
276
            str += typeid(e).name();
277
            str += " ";
278
            if (msg) {str += msg;}
279
            else     {str += "No OCCT Exception Message";}
280
            throw Py::RuntimeError(str);
281
        }
282
        catch (const Base::Exception &e) {
283
            throw Py::RuntimeError(e.what());
284
        }
285
        catch (...) {
286
            throw Py::RuntimeError("Unknown C++ exception");
287
        }
288
    }
289
#if defined(HAVE_PCL_SURFACE)
290
    /*
291
import ReverseEngineering as Reen
292
import Points
293
import Mesh
294
import random
295

296
r=random.Random()
297

298
p=Points.Points()
299
pts=[]
300
for i in range(21):
301
  for j in range(21):
302
    pts.append(App.Vector(i,j,r.gauss(5,0.05)))
303

304
p.addPoints(pts)
305
m=Reen.triangulate(Points=p,SearchRadius=2.2)
306
Mesh.show(m)
307
    */
308
    Py::Object triangulate(const Py::Tuple& args, const Py::Dict& kwds)
309
    {
310
        PyObject *pts;
311
        double searchRadius;
312
        PyObject *vec = 0;
313
        int ksearch=5;
314
        double mu=2.5;
315

316
        static const std::array<const char*,6> kwds_greedy {"Points", "SearchRadius", "Mu", "KSearch",
317
                                      "Normals", NULL};
318
        if (!Base::Wrapped_ParseTupleAndKeywords(args.ptr(), kwds.ptr(), "O!d|diO", kwds_greedy,
319
                                        &(Points::PointsPy::Type), &pts,
320
                                        &searchRadius, &mu, &ksearch, &vec))
321
            throw Py::Exception();
322

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

325
        Mesh::MeshObject* mesh = new Mesh::MeshObject();
326
        SurfaceTriangulation tria(*points, *mesh);
327
        tria.setMu(mu);
328
        tria.setSearchRadius(searchRadius);
329
        if (vec) {
330
            Py::Sequence list(vec);
331
            std::vector<Base::Vector3f> normals;
332
            normals.reserve(list.size());
333
            for (Py::Sequence::iterator it = list.begin(); it != list.end(); ++it) {
334
                Base::Vector3d v = Py::Vector(*it).toVector();
335
                normals.push_back(Base::convertTo<Base::Vector3f>(v));
336
            }
337
            tria.perform(normals);
338
        }
339
        else {
340
            tria.perform(ksearch);
341
        }
342

343
        return Py::asObject(new Mesh::MeshPy(mesh));
344
    }
345
    Py::Object poissonReconstruction(const Py::Tuple& args, const Py::Dict& kwds)
346
    {
347
        PyObject *pts;
348
        PyObject *vec = 0;
349
        int ksearch=5;
350
        int octreeDepth=-1;
351
        int solverDivide=-1;
352
        double samplesPerNode=-1.0;
353

354
        static const std::array<const char*,7> kwds_poisson {"Points", "KSearch", "OctreeDepth", "SolverDivide",
355
                                      "SamplesPerNode", "Normals", NULL};
356
        if (!Base::Wrapped_ParseTupleAndKeywords(args.ptr(), kwds.ptr(), "O!|iiidO", kwds_poisson,
357
                                        &(Points::PointsPy::Type), &pts,
358
                                        &ksearch, &octreeDepth, &solverDivide, &samplesPerNode, &vec))
359
            throw Py::Exception();
360

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

363
        Mesh::MeshObject* mesh = new Mesh::MeshObject();
364
        Reen::PoissonReconstruction poisson(*points, *mesh);
365
        poisson.setDepth(octreeDepth);
366
        poisson.setSolverDivide(solverDivide);
367
        poisson.setSamplesPerNode(samplesPerNode);
368
        if (vec) {
369
            Py::Sequence list(vec);
370
            std::vector<Base::Vector3f> normals;
371
            normals.reserve(list.size());
372
            for (Py::Sequence::iterator it = list.begin(); it != list.end(); ++it) {
373
                Base::Vector3d v = Py::Vector(*it).toVector();
374
                normals.push_back(Base::convertTo<Base::Vector3f>(v));
375
            }
376
            poisson.perform(normals);
377
        }
378
        else {
379
            poisson.perform(ksearch);
380
        }
381

382
        return Py::asObject(new Mesh::MeshPy(mesh));
383
    }
384
   /*
385
import ReverseEngineering as Reen
386
import Points
387
import Mesh
388
import random
389
import math
390
r=random.Random()
391
p=Points.Points()
392
pts=[]
393
for i in range(21):
394
  for j in range(21):
395
    pts.append(App.Vector(i,j,r.random()))
396
p.addPoints(pts)
397
m=Reen.viewTriangulation(p,21,21)
398
Mesh.show(m)
399
def boxmueller():
400
  r1,r2=random.random(),random.random()
401
  return math.sqrt(-2*math.log(r1))*math.cos(2*math.pi*r2)
402
p=Points.Points()
403
pts=[]
404
for i in range(21):
405
  for j in range(21):
406
    pts.append(App.Vector(i,j,r.gauss(5,0.05)))
407
p.addPoints(pts)
408
m=Reen.viewTriangulation(p,21,21)
409
Mesh.show(m)
410
    */
411
    Py::Object viewTriangulation(const Py::Tuple& args, const Py::Dict& kwds)
412
    {
413
        PyObject *pts;
414
        int width;
415
        int height;
416

417
        static const std::array<const char*,4> kwds_view {"Points", "Width", "Height", NULL};
418
        if (!Base::Wrapped_ParseTupleAndKeywords(args.ptr(), kwds.ptr(), "O!|ii", kwds_view,
419
                                        &(Points::PointsPy::Type), &pts,
420
                                        &width, &height))
421
            throw Py::Exception();
422

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

425
        try {
426
            Mesh::MeshObject* mesh = new Mesh::MeshObject();
427
            ImageTriangulation view(width, height, *points, *mesh);
428
            view.perform();
429

430
            return Py::asObject(new Mesh::MeshPy(mesh));
431
        }
432
        catch (const Base::Exception& e) {
433
            throw Py::RuntimeError(e.what());
434
        }
435
    }
436
    Py::Object gridProjection(const Py::Tuple& args, const Py::Dict& kwds)
437
    {
438
        PyObject *pts;
439
        PyObject *vec = 0;
440
        int ksearch=5;
441

442
        static const std::array<const char*,4> kwds_greedy {"Points", "KSearch", "Normals", NULL};
443
        if (!Base::Wrapped_ParseTupleAndKeywords(args.ptr(), kwds.ptr(), "O!|iO", kwds_greedy,
444
                                        &(Points::PointsPy::Type), &pts,
445
                                        &ksearch, &vec))
446
            throw Py::Exception();
447

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

450
        Mesh::MeshObject* mesh = new Mesh::MeshObject();
451
        GridReconstruction tria(*points, *mesh);
452
        if (vec) {
453
            Py::Sequence list(vec);
454
            std::vector<Base::Vector3f> normals;
455
            normals.reserve(list.size());
456
            for (Py::Sequence::iterator it = list.begin(); it != list.end(); ++it) {
457
                Base::Vector3d v = Py::Vector(*it).toVector();
458
                normals.push_back(Base::convertTo<Base::Vector3f>(v));
459
            }
460
            tria.perform(normals);
461
        }
462
        else {
463
            tria.perform(ksearch);
464
        }
465

466
        return Py::asObject(new Mesh::MeshPy(mesh));
467
    }
468
    Py::Object marchingCubesRBF(const Py::Tuple& args, const Py::Dict& kwds)
469
    {
470
        PyObject *pts;
471
        PyObject *vec = 0;
472
        int ksearch=5;
473

474
        static const std::array<const char*,4> kwds_greedy {"Points", "KSearch", "Normals", NULL};
475
        if (!Base::Wrapped_ParseTupleAndKeywords(args.ptr(), kwds.ptr(), "O!|iO", kwds_greedy,
476
                                        &(Points::PointsPy::Type), &pts,
477
                                        &ksearch, &vec))
478
            throw Py::Exception();
479

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

482
        Mesh::MeshObject* mesh = new Mesh::MeshObject();
483
        MarchingCubesRBF tria(*points, *mesh);
484
        if (vec) {
485
            Py::Sequence list(vec);
486
            std::vector<Base::Vector3f> normals;
487
            normals.reserve(list.size());
488
            for (Py::Sequence::iterator it = list.begin(); it != list.end(); ++it) {
489
                Base::Vector3d v = Py::Vector(*it).toVector();
490
                normals.push_back(Base::convertTo<Base::Vector3f>(v));
491
            }
492
            tria.perform(normals);
493
        }
494
        else {
495
            tria.perform(ksearch);
496
        }
497

498
        return Py::asObject(new Mesh::MeshPy(mesh));
499
    }
500
/*
501
import ReverseEngineering as Reen
502
import Points
503
import Mesh
504
import random
505
r=random.Random()
506
p=Points.Points()
507
pts=[]
508
for i in range(21):
509
  for j in range(21):
510
    pts.append(App.Vector(i,j,r.gauss(5,0.05)))
511
p.addPoints(pts)
512
m=Reen.marchingCubesHoppe(Points=p)
513
Mesh.show(m)
514
    */
515
    Py::Object marchingCubesHoppe(const Py::Tuple& args, const Py::Dict& kwds)
516
    {
517
        PyObject *pts;
518
        PyObject *vec = 0;
519
        int ksearch=5;
520

521
        static const std::array<const char*,4> kwds_greedy {"Points", "KSearch", "Normals", NULL};
522
        if (!Base::Wrapped_ParseTupleAndKeywords(args.ptr(), kwds.ptr(), "O!|iO", kwds_greedy,
523
                                        &(Points::PointsPy::Type), &pts,
524
                                        &ksearch, &vec))
525
            throw Py::Exception();
526

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

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

545
        return Py::asObject(new Mesh::MeshPy(mesh));
546
    }
547
#endif
548
#if defined(HAVE_PCL_OPENNURBS)
549
    Py::Object fitBSpline(const Py::Tuple& args, const Py::Dict& kwds)
550
    {
551
        PyObject *pts;
552
        int degree = 2;
553
        int refinement = 4;
554
        int iterations = 10;
555
        double interiorSmoothness = 0.2;
556
        double interiorWeight = 1.0;
557
        double boundarySmoothness = 0.2;
558
        double boundaryWeight = 0.0;
559

560
        static const std::array<const char*,9> kwds_approx {"Points", "Degree", "Refinement", "Iterations",
561
                                      "InteriorSmoothness", "InteriorWeight", "BoundarySmoothness", "BoundaryWeight", NULL};
562
        if (!Base::Wrapped_ParseTupleAndKeywords(args.ptr(), kwds.ptr(), "O!|iiidddd", kwds_approx,
563
                                        &(Points::PointsPy::Type), &pts,
564
                                        &degree, &refinement, &iterations,
565
                                        &interiorSmoothness, &interiorWeight,
566
                                        &boundarySmoothness, &boundaryWeight))
567
            throw Py::Exception();
568

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

571
        BSplineFitting fit(points->getBasicPoints());
572
        fit.setOrder(degree+1);
573
        fit.setRefinement(refinement);
574
        fit.setIterations(iterations);
575
        fit.setInteriorSmoothness(interiorSmoothness);
576
        fit.setInteriorWeight(interiorWeight);
577
        fit.setBoundarySmoothness(boundarySmoothness);
578
        fit.setBoundaryWeight(boundaryWeight);
579
        Handle(Geom_BSplineSurface) hSurf = fit.perform();
580

581
        if (!hSurf.IsNull()) {
582
            return Py::asObject(new Part::BSplineSurfacePy(new Part::GeomBSplineSurface(hSurf)));
583
        }
584

585
        throw Py::RuntimeError("Computation of B-spline surface failed");
586
    }
587
#endif
588
#if defined(HAVE_PCL_FILTERS)
589
    Py::Object filterVoxelGrid(const Py::Tuple& args, const Py::Dict& kwds)
590
    {
591
        PyObject *pts;
592
        double voxDimX = 0;
593
        double voxDimY = 0;
594
        double voxDimZ = 0;
595

596
        static const std::array<const char*,5>  kwds_voxel {"Points", "DimX", "DimY", "DimZ", NULL};
597
        if (!Base::Wrapped_ParseTupleAndKeywords(args.ptr(), kwds.ptr(), "O!d|dd", kwds_voxel,
598
                                        &(Points::PointsPy::Type), &pts,
599
                                        &voxDimX, &voxDimY, &voxDimZ))
600
            throw Py::Exception();
601

602
        if (voxDimY == 0)
603
            voxDimY = voxDimX;
604

605
        if (voxDimZ == 0)
606
            voxDimZ = voxDimX;
607

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

610
        pcl::PointCloud<pcl::PointXYZ>::Ptr cloud (new pcl::PointCloud<pcl::PointXYZ>);
611
        cloud->reserve(points->size());
612
        for (Points::PointKernel::const_iterator it = points->begin(); it != points->end(); ++it) {
613
            cloud->push_back(pcl::PointXYZ(it->x, it->y, it->z));
614
        }
615

616
        // Create the filtering object
617
        pcl::PointCloud<pcl::PointXYZ>::Ptr cloud_downSmpl (new pcl::PointCloud<pcl::PointXYZ>);
618
        pcl::VoxelGrid<pcl::PointXYZ> voxG;
619
        voxG.setInputCloud (cloud);
620
        voxG.setLeafSize (voxDimX, voxDimY, voxDimZ);
621
        voxG.filter (*cloud_downSmpl);
622

623
        Points::PointKernel* points_sample = new Points::PointKernel();
624
        points_sample->reserve(cloud_downSmpl->size());
625
        for (pcl::PointCloud<pcl::PointXYZ>::const_iterator it = cloud_downSmpl->begin();it!=cloud_downSmpl->end();++it) {
626
            points_sample->push_back(Base::Vector3d(it->x,it->y,it->z));
627
        }
628

629
        return Py::asObject(new Points::PointsPy(points_sample));
630
    }
631
#endif
632
#if defined(HAVE_PCL_FILTERS)
633
    Py::Object normalEstimation(const Py::Tuple& args, const Py::Dict& kwds)
634
    {
635
        PyObject *pts;
636
        int ksearch=0;
637
        double searchRadius=0;
638

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

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

647
        std::vector<Base::Vector3d> normals;
648
        NormalEstimation estimate(*points);
649
        estimate.setKSearch(ksearch);
650
        estimate.setSearchRadius(searchRadius);
651
        estimate.perform(normals);
652

653
        Py::List list;
654
        for (std::vector<Base::Vector3d>::iterator it = normals.begin(); it != normals.end(); ++it) {
655
            list.append(Py::Vector(*it));
656
        }
657

658
        return list;
659
    }
660
#endif
661
#if defined(HAVE_PCL_SEGMENTATION)
662
    Py::Object regionGrowingSegmentation(const Py::Tuple& args, const Py::Dict& kwds)
663
    {
664
        PyObject *pts;
665
        PyObject *vec = 0;
666
        int ksearch=5;
667

668
        static const std::array<const char*,4> kwds_segment {"Points", "KSearch", "Normals", NULL};
669
        if (!Base::Wrapped_ParseTupleAndKeywords(args.ptr(), kwds.ptr(), "O!|iO", kwds_segment,
670
                                        &(Points::PointsPy::Type), &pts,
671
                                        &ksearch, &vec))
672
            throw Py::Exception();
673

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

676
        std::list<std::vector<int> > clusters;
677
        RegionGrowing segm(*points, clusters);
678
        if (vec) {
679
            Py::Sequence list(vec);
680
            std::vector<Base::Vector3f> normals;
681
            normals.reserve(list.size());
682
            for (Py::Sequence::iterator it = list.begin(); it != list.end(); ++it) {
683
                Base::Vector3d v = Py::Vector(*it).toVector();
684
                normals.push_back(Base::convertTo<Base::Vector3f>(v));
685
            }
686
            segm.perform(normals);
687
        }
688
        else {
689
            segm.perform(ksearch);
690
        }
691

692
        Py::List lists;
693
        for (std::list<std::vector<int> >::iterator it = clusters.begin(); it != clusters.end(); ++it) {
694
            Py::Tuple tuple(it->size());
695
            for (std::size_t i = 0; i < it->size(); i++) {
696
                tuple.setItem(i, Py::Long((*it)[i]));
697
            }
698
            lists.append(tuple);
699
        }
700

701
        return lists;
702
    }
703
    Py::Object featureSegmentation(const Py::Tuple& args, const Py::Dict& kwds)
704
    {
705
        PyObject *pts;
706
        int ksearch=5;
707

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

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

715
        std::list<std::vector<int> > clusters;
716
        Segmentation segm(*points, clusters);
717
        segm.perform(ksearch);
718

719
        Py::List lists;
720
        for (std::list<std::vector<int> >::iterator it = clusters.begin(); it != clusters.end(); ++it) {
721
            Py::Tuple tuple(it->size());
722
            for (std::size_t i = 0; i < it->size(); i++) {
723
                tuple.setItem(i, Py::Long((*it)[i]));
724
            }
725
            lists.append(tuple);
726
        }
727

728
        return lists;
729
    }
730
#endif
731
#if defined(HAVE_PCL_SAMPLE_CONSENSUS)
732
/*
733
import ReverseEngineering as reen
734
import Points
735
import Part
736
p = App.ActiveDocument.Points.Points
737
data = p.Points
738
n = reen.normalEstimation(p, 10)
739
model = reen.sampleConsensus(SacModel="Plane", Points=p)
740
indices = model["Model"]
741
param = model["Parameters"]
742
plane = Part.Plane()
743
plane.Axis = param[0:3]
744
plane.Position = -plane.Axis * param[3]
745
np = Points.Points()
746
np.addPoints([data[i] for i in indices])
747
Points.show(np)
748
# sort in descending order
749
indices = list(indices)
750
indices.sort(reverse=True)
751
# remove points of segment
752
for i in indices:
753
    del data[i]
754
    del n[i]
755
p = Points.Points()
756
p.addPoints(data)
757
model = reen.sampleConsensus(SacModel="Cylinder", Points=p, Normals=n)
758
indices = model["Model"]
759
np = Points.Points()
760
np.addPoints([data[i] for i in indices])
761
Points.show(np)
762
    */
763
    Py::Object sampleConsensus(const Py::Tuple& args, const Py::Dict& kwds)
764
    {
765
        PyObject *pts;
766
        PyObject *vec = nullptr;
767
        const char* sacModelType = nullptr;
768

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

774
        Points::PointKernel* points = static_cast<Points::PointsPy*>(pts)->getPointKernelPtr();
775
        std::vector<Base::Vector3d> normals;
776
        if (vec) {
777
            Py::Sequence list(vec);
778
            normals.reserve(list.size());
779
            for (Py::Sequence::iterator it = list.begin(); it != list.end(); ++it) {
780
                Base::Vector3d v = Py::Vector(*it).toVector();
781
                normals.push_back(v);
782
            }
783
        }
784

785
        SampleConsensus::SacModel sacModel = SampleConsensus::SACMODEL_PLANE;
786
        if (sacModelType) {
787
            if (strcmp(sacModelType, "Cylinder") == 0)
788
                sacModel = SampleConsensus::SACMODEL_CYLINDER;
789
            else if (strcmp(sacModelType, "Sphere") == 0)
790
                sacModel = SampleConsensus::SACMODEL_SPHERE;
791
            else if (strcmp(sacModelType, "Cone") == 0)
792
                sacModel = SampleConsensus::SACMODEL_CONE;
793
        }
794

795
        std::vector<float> parameters;
796
        SampleConsensus sample(sacModel, *points, normals);
797
        std::vector<int> model;
798
        double probability = sample.perform(parameters, model);
799

800
        Py::Dict dict;
801
        Py::Tuple tuple(parameters.size());
802
        for (std::size_t i = 0; i < parameters.size(); i++)
803
            tuple.setItem(i, Py::Float(parameters[i]));
804
        Py::Tuple data(model.size());
805
        for (std::size_t i = 0; i < model.size(); i++)
806
            data.setItem(i, Py::Long(model[i]));
807
        dict.setItem(Py::String("Probability"), Py::Float(probability));
808
        dict.setItem(Py::String("Parameters"), tuple);
809
        dict.setItem(Py::String("Model"), data);
810

811
        return dict;
812
    }
813
#endif
814
};
815

816
PyObject* initModule()
817
{
818
    return Base::Interpreter().addModule(new Module);
819
}
820

821
} // namespace Reen
822

823

824
/* Python entry */
825
PyMOD_INIT_FUNC(ReverseEngineering)
826
{
827
    // load dependent module
828
    try {
829
        Base::Interpreter().loadModule("Part");
830
        Base::Interpreter().loadModule("Mesh");
831
    }
832
    catch(const Base::Exception& e) {
833
        PyErr_SetString(PyExc_ImportError, e.what());
834
        PyMOD_Return(nullptr);
835
    }
836

837
    PyObject* mod = Reen::initModule();
838
    Base::Console().Log("Loading ReverseEngineering module... done\n");
839
    PyMOD_Return(mod);
840
}
841
// clang-format on
842