1
/*
2
* Copyright (C) 2017 German Aerospace Center (DLR/SC)
3
*
4
* Created: 2017-05-24 Merlin Pelz <Merlin.Pelz@dlr.de>
5
*
6
* Licensed under the Apache License, Version 2.0 (the "License");
7
* you may not use this file except in compliance with the License.
8
* You may obtain a copy of the License at
9
*
10
*     http://www.apache.org/licenses/LICENSE-2.0
11
*
12
* Unless required by applicable law or agreed to in writing, software
13
* distributed under the License is distributed on an "AS IS" BASIS,
14
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
15
* See the License for the specific language governing permissions and
16
* limitations under the License.
17
*/
18

19
#include "CTiglBSplineAlgorithms.h"
20
#include "CTiglCurveNetworkSorter.h"
21
#include "CTiglCurvesToSurface.h"
22
#include "CTiglError.h"
23
#include "CSharedPtr.h"
24
#include "CTiglBSplineApproxInterp.h"
25
#include "CTiglPointsToBSplineInterpolation.h"
26
#include "to_string.h"
27
#include "tiglcommonfunctions.h"
28
#include "Debugging.h"
29

30
#include <Standard_Version.hxx>
31
#include <Geom2d_BSplineCurve.hxx>
32
#include <Geom_BSplineCurve.hxx>
33
#include <Geom_BSplineSurface.hxx>
34
#include <Geom_TrimmedCurve.hxx>
35
#include <GeomConvert.hxx>
36
#include <Geom2dAPI_Interpolate.hxx>
37
#include <GeomAPI_ExtremaCurveCurve.hxx>
38
#include <TColStd_Array2OfReal.hxx>
39
#include <TColStd_HArray1OfReal.hxx>
40
#include <TColStd_HArray1OfInteger.hxx>
41
#include <TColgp_HArray1OfPnt.hxx>
42
#include <TColgp_HArray1OfPnt2d.hxx>
43
#include <TColgp_Array1OfPnt2d.hxx>
44
#include <TColgp_HArray2OfPnt.hxx>
45
#include <BSplCLib.hxx>
46
#include <GeomAPI_PointsToBSpline.hxx>
47
#include <BRepTools.hxx>
48
#include <BRepBuilderAPI_MakeFace.hxx>
49
#include <BRepBuilderAPI_MakeEdge.hxx>
50
#include <Precision.hxx>
51
#include <Geom2dAPI_ProjectPointOnCurve.hxx>
52
#include <GCPnts_AbscissaPoint.hxx>
53

54
#include <cmath>
55
#include <stdexcept>
56
#include <algorithm>
57
#include <cassert>
58

59
namespace
60
{
61

62
    class helper_function_unique
63
    {
64
    public:
65
        helper_function_unique(double tolerance = 1e-15)
66
            : _tol(tolerance)
67
        {}
68

69
        // helper function for std::unique
70
        bool operator()(double a, double b)
71
        {
72
            return (fabs(a - b) < _tol);
73
        }
74
    private:
75
        double _tol;
76
    };
77
    
78
    enum SurfAdapterDir
79
    {
80
        udir = 0,
81
        vdir = 1
82
    };
83
    
84 1
    class SurfAdapterView
85
    {
86
    public:
87 1
        SurfAdapterView(Handle(Geom_BSplineSurface) surf, SurfAdapterDir dir)
88 1
            : _surf(surf), _dir(dir)
89
        {
90
        }
91

92 1
        void insertKnot(double knot, int mult, double tolerance=1e-15)
93
        {
94 1
            if (_dir == udir) {
95 1
                _surf->InsertUKnot(knot, mult, tolerance, false);
96
            }
97
            else {
98 1
                _surf->InsertVKnot(knot, mult, tolerance, false);
99
            }
100
        }
101

102 1
        double getKnot(int idx) const
103
        {
104 1
            if (_dir == udir) {
105 1
                return _surf->UKnot(idx);
106
            }
107
            else {
108 1
                return _surf->VKnot(idx);
109
            }
110
        }
111

112 1
        int getMult(int idx) const
113
        {
114 1
            if (_dir == udir) {
115 1
                return _surf->UMultiplicity(idx);
116
            }
117
            else {
118 1
                return _surf->VMultiplicity(idx);
119
            }
120
        }
121

122 1
        int getNKnots() const
123
        {
124 1
            if (_dir == udir) {
125 1
                return _surf->NbUKnots();
126
            }
127
            else {
128 1
                return _surf->NbVKnots();
129
            }
130
        }
131

132 1
        int getDegree() const
133
        {
134 1
            if (_dir == udir) {
135 1
                return _surf->UDegree();
136
            }
137
            else {
138 1
                return _surf->VDegree();
139
            }
140
        }
141
        
142 1
        void setDir(SurfAdapterDir dir)
143
        {
144 1
            _dir = dir;
145
        }
146

147 1
        operator const Handle(Geom_BSplineSurface)&() const
148
        {
149 1
            return _surf;
150
        }
151
    private:
152
        Handle(Geom_BSplineSurface) _surf;
153
        SurfAdapterDir _dir;
154
    };
155
    
156 1
    class CurveAdapterView
157
    {
158
    public:
159 1
        CurveAdapterView(Handle(Geom_BSplineCurve) curve)
160 1
            : _curve(curve)
161
        {
162
        }
163

164 1
        void insertKnot(double knot, int mult, double tolerance=1e-15)
165
        {
166 1
            _curve->InsertKnot(knot, mult, tolerance, false);
167
        }
168

169 1
        double getKnot(int idx) const
170
        {
171 1
            return _curve->Knot(idx);
172
        }
173

174 1
        int getMult(int idx) const
175
        {
176 1
            return _curve->Multiplicity(idx);
177
        }
178

179 1
        int getNKnots() const
180
        {
181 1
            return _curve->NbKnots();
182
        }
183

184 1
        int getDegree() const
185
        {
186 1
            return _curve->Degree();
187
        }
188

189 1
        operator const Handle(Geom_BSplineCurve)&() const
190
        {
191 1
            return _curve;
192
        }
193
 
194
    private:
195
        Handle(Geom_BSplineCurve) _curve;
196
    };
197

198
    template <class SplineAdapter>
199 1
    bool haveSameRange(const std::vector<SplineAdapter>& splines_vector, double par_tolerance)
200
    {
201 1
        double begin_param_dir = splines_vector[0].getKnot(1);
202 1
        double end_param_dir = splines_vector[0].getKnot(splines_vector[0].getNKnots());
203 1
        for (unsigned int spline_idx = 1; spline_idx < splines_vector.size(); ++spline_idx) {
204 1
            const SplineAdapter& curSpline = splines_vector[spline_idx];
205 1
            double begin_param_dir_surface = curSpline.getKnot(1);
206 1
            double end_param_dir_surface = curSpline.getKnot(curSpline.getNKnots());
207 1
            if (std::abs(begin_param_dir_surface - begin_param_dir) > par_tolerance || std::abs(end_param_dir_surface - end_param_dir) > par_tolerance) {
208 1
                return false;
209
            }
210
        }
211 1
        return true;
212
    }
213

214
    template <class SplineAdapter>
215 1
    bool haveSameDegree(const std::vector<SplineAdapter>& splines)
216
    {
217 1
        int degree = splines[0].getDegree();
218 1
        for (unsigned int splineIdx = 1; splineIdx < splines.size(); ++splineIdx) {
219 1
            if (splines[splineIdx].getDegree() != degree) {
220 0
                return false;
221
            }
222
        }
223 1
        return true;
224
    }
225

226
    /**
227
     * @brief createCommonKnotsVectorImpl:
228
     *          Creates a common knot vector in u- or v-direction of the given vector of B-splines
229
     *          The common knot vector contains all knots in u- or v-direction of all splines with the highest multiplicity of all splines.
230
     * @param old_splines_vector:
231
     *          the given vector of B-spline splines that could have a different knot vector in u- or v-direction
232
     */
233
    template <class SplineAdapter>
234 1
    void makeGeometryCompatibleImpl(std::vector<SplineAdapter>& splines_vector, double par_tolerance)
235
    {
236
        // all B-spline splines must have the same parameter range in the chosen direction
237 1
        if (!haveSameRange(splines_vector, par_tolerance)) {
238 1
            throw tigl::CTiglError("B-splines don't have the same parameter range at least in one direction (u / v) in method createCommonKnotsVectorImpl!", TIGL_MATH_ERROR);
239
        }
240

241
        // all B-spline splines must have the same degree in the chosen direction
242 1
        if (!haveSameDegree(splines_vector)) {
243 0
            throw tigl::CTiglError("B-splines don't have the same degree at least in one direction (u / v) in method createCommonKnotsVectorImpl!", TIGL_MATH_ERROR);
244
        }
245

246
        // The parametric tolerance must be smaller than half of the minimum knot distance
247 1
        for (typename std::vector<SplineAdapter>::const_iterator splineIt = splines_vector.begin(); splineIt != splines_vector.end(); ++splineIt) {
248 1
            const SplineAdapter& spline = *splineIt;
249 1
            for (int iKnot = 1; iKnot < spline.getNKnots(); ++iKnot) {
250 1
                double knotDist = spline.getKnot(iKnot+1) - spline.getKnot(iKnot);
251 1
                par_tolerance = std::min(par_tolerance, knotDist / 2.);
252
            }
253
        }
254

255
        // insert all knots in first spline
256 1
        SplineAdapter& firstSpline = splines_vector[0];
257 1
        for (typename std::vector<SplineAdapter>::const_iterator splineIt = splines_vector.begin()+1; splineIt != splines_vector.end(); ++splineIt) {
258 1
            const SplineAdapter& spline = *splineIt;
259 1
            for (int knot_idx = 2; knot_idx < spline.getNKnots(); ++knot_idx) {
260 1
                double knot = spline.getKnot(knot_idx);
261 1
                int mult = spline.getMult(knot_idx);
262 1
                firstSpline.insertKnot(knot, mult, par_tolerance);
263
            }
264
        }
265

266

267
        // now insert knots from first into all others
268 1
        for (typename std::vector<SplineAdapter>::iterator splineIt = splines_vector.begin()+1; splineIt != splines_vector.end(); ++splineIt) {
269 1
            SplineAdapter& spline = *splineIt;
270 1
            for (int knot_idx = 2; knot_idx < firstSpline.getNKnots(); ++knot_idx) {
271 1
                double knot = firstSpline.getKnot(knot_idx);
272 1
                int mult = firstSpline.getMult(knot_idx);
273 1
                spline.insertKnot(knot, mult, par_tolerance);
274
            }
275 1
            if (spline.getNKnots() != firstSpline.getNKnots()) {
276 0
                throw tigl::CTiglError("Unexpected error in Algorithm makeGeometryCompatibleImpl.\nPlease contact the developers.");
277
            }
278
        }
279

280

281
    } // makeGeometryCompatibleImpl
282
    
283
    template <class OccMatrix, class OccVector, class OccHandleVector>
284 1
    OccHandleVector array2GetColumn(const OccMatrix& matrix, int colIndex)
285
    {
286 1
        OccHandleVector colVector =  new OccVector(matrix.LowerRow(), matrix.UpperRow());
287

288 1
        for (int rowIdx = matrix.LowerRow(); rowIdx <= matrix.UpperRow(); ++rowIdx) {
289 1
            colVector->SetValue(rowIdx, matrix(rowIdx, colIndex));
290
        }
291

292 1
        return colVector;
293
    }
294
    
295
    template <class OccMatrix, class OccVector, class OccHandleVector>
296 1
    OccHandleVector array2GetRow(const OccMatrix& matrix, int rowIndex)
297
    {
298 1
        OccHandleVector rowVector = new OccVector(matrix.LowerCol(), matrix.UpperCol());
299
        
300 1
        for (int colIdx = matrix.LowerCol(); colIdx <= matrix.UpperCol(); ++colIdx) {
301 1
            rowVector->SetValue(colIdx, matrix(rowIndex, colIdx));
302
        }
303
        
304 1
        return rowVector;
305
    }
306

307 1
    Handle_TColgp_HArray1OfPnt pntArray2GetColumn(const TColgp_Array2OfPnt& matrix, int colIndex)
308
    {
309 1
        return array2GetColumn<TColgp_Array2OfPnt, TColgp_HArray1OfPnt, Handle_TColgp_HArray1OfPnt>(matrix, colIndex);
310
    }
311

312 1
    Handle_TColgp_HArray1OfPnt pntArray2GetRow(const TColgp_Array2OfPnt& matrix, int rowIndex)
313
    {
314 1
        return array2GetRow<TColgp_Array2OfPnt, TColgp_HArray1OfPnt, Handle_TColgp_HArray1OfPnt>(matrix, rowIndex);
315
    }
316
    		
317
} // namespace
318

319

320
namespace tigl
321
{
322

323
const double CTiglBSplineAlgorithms::REL_TOL_CLOSED = 1e-8;
324

325 1
bool CTiglBSplineAlgorithms::isUDirClosed(const TColgp_Array2OfPnt& points, double tolerance)
326
{
327 1
    bool uDirClosed = true;
328 1
    int ulo = points.LowerRow();
329 1
    int uhi = points.UpperRow();
330
    // check that first row and last row are the same
331 1
    for (int v_idx = points.LowerCol(); v_idx <= points.UpperCol(); ++v_idx) {
332 1
        gp_Pnt pfirst = points.Value(ulo, v_idx);
333 1
        gp_Pnt pLast = points.Value(uhi, v_idx);
334 1
        uDirClosed = uDirClosed & pfirst.IsEqual(pLast, tolerance);
335
    }
336 1
    return uDirClosed;
337
}
338

339 1
bool CTiglBSplineAlgorithms::isVDirClosed(const TColgp_Array2OfPnt& points, double tolerance)
340
{
341 1
    bool vDirClosed = true;
342 1
    int vlo = points.LowerCol();
343 1
    int vhi = points.UpperCol();
344 1
    for (int u_idx = points.LowerRow(); u_idx <= points.UpperRow(); ++u_idx) {
345 1
        vDirClosed = vDirClosed & points.Value(u_idx, vlo).IsEqual(points.Value(u_idx, vhi), tolerance);
346
    }
347 1
    return vDirClosed;
348
}
349

350 1
std::vector<double> CTiglBSplineAlgorithms::knotsFromCurveParameters(std::vector<double> &params, unsigned int degree, bool closedCurve)
351
{
352 1
    if (params.size() < 2) {
353 0
        throw CTiglError("Parameters must contain two or more elements.");
354
    }
355

356 1
    size_t nCP = params.size();
357 1
    if (closedCurve) {
358
        // For each continuity condition, we have to add one control point
359 1
        nCP += degree - 1;
360
    }
361 1
    size_t nInnerKnots = nCP - degree + 1;
362

363 1
    std::vector<double> innerKnots(nInnerKnots);
364 1
    innerKnots.front() = params.front();
365 1
    innerKnots.back() = params.back();
366

367 1
    std::vector<double> knots;
368

369 1
    if (closedCurve && degree % 2 == 0) {
370 1
        size_t m = params.size() - 2;
371

372
        // build difference vector
373 1
        std::vector<double> dparm(m + 1, 0.);
374 1
        for (size_t iparm = 0; iparm <= m; ++iparm) {
375 1
            dparm[iparm] = (params[iparm+1] - params[iparm]);
376
        }
377

378 1
        innerKnots[1] = innerKnots[0] + 0.5 * (dparm[0] + dparm[m]);
379 1
        for (size_t iparm = 1; iparm < m; ++iparm) {
380 1
            innerKnots[iparm + 1] = innerKnots[iparm] + 0.5 * (dparm[iparm-1] + dparm[iparm]);
381
        }
382

383
        // shift parameters
384 1
        for (size_t iparm = 0; iparm < params.size(); ++iparm) {
385 1
            params[iparm] += dparm[m] / 2.;
386
        }
387

388
    }
389 1
    else if (closedCurve) {
390 1
        assert(innerKnots.size() == params.size());
391 1
        innerKnots = params;
392
    }
393
    else {
394
        // averaging
395 1
        for (size_t j = 1; j < params.size() - degree; ++j) {
396 1
            double sum = 0.;
397
            // average
398 1
            for (size_t i = j; i <= j + degree - 1; ++i) {
399 1
                sum += params[i];
400
            }
401

402 1
            innerKnots[j] = sum / static_cast<double>(degree);
403
        }
404
    }
405

406 1
    if (closedCurve) {
407 1
        double offset = innerKnots[0] - innerKnots[nInnerKnots -1];
408 1
        for (size_t iknot = 0; iknot < degree; ++iknot) {
409 1
            knots.push_back(offset + innerKnots[nInnerKnots - degree -1 + iknot]);
410
        }
411 1
        for (size_t iknot = 0; iknot < nInnerKnots; ++iknot) {
412 1
            knots.push_back(innerKnots[iknot]);
413
        }
414 1
        for (size_t iknot = 0; iknot < degree; ++iknot) {
415 1
            knots.push_back(-offset + innerKnots[iknot + 1]);
416
        }
417
    }
418
    else {
419 1
        for (size_t iknot = 0; iknot < degree; ++iknot) {
420 1
            knots.push_back(innerKnots[0]);
421
        }
422 1
        for (size_t iknot = 0; iknot < nInnerKnots; ++iknot) {
423 1
            knots.push_back(innerKnots[iknot]);
424
        }
425 1
        for (size_t iknot = 0; iknot < degree; ++iknot) {
426 1
            knots.push_back(innerKnots[nInnerKnots - 1]);
427
        }
428
    }
429

430 1
    if (closedCurve && degree <= 1) {
431 1
        size_t nKnots = knots.size();
432 1
        knots[0] = knots[1];
433 1
        knots[nKnots - 1] = knots[nKnots - 2];
434
    }
435

436 1
    return knots;
437
}
438

439 1
double CTiglBSplineAlgorithms::scale(const TColgp_Array2OfPnt& points)
440
{
441 1
    double theScale = 0.;
442 1
    for (int uidx = points.LowerRow(); uidx <= points.UpperRow(); ++uidx) {
443 1
        gp_Pnt pFirst = points.Value(uidx, points.LowerCol());
444 1
        for (int vidx = points.LowerCol() + 1; vidx <= points.UpperCol(); ++vidx) {
445 1
            double dist = pFirst.Distance(points.Value(uidx, vidx));
446 1
            theScale = std::max(theScale, dist);
447
        }
448
    }
449 1
    return theScale;
450
}
451

452 1
double CTiglBSplineAlgorithms::scale(const TColgp_Array1OfPnt& points)
453
{
454 1
    double theScale = 0.;
455

456 1
    for (int i = points.Lower(); i <= points.Upper(); ++i) {
457 1
        for (int j = i + 1; j < points.Upper(); ++j) {
458 1
            double dist = points.Value(i).Distance(points.Value(j));
459 1
            theScale = std::max(theScale, dist);
460
        }
461
    }
462 1
    return theScale;
463
}
464

465 1
std::vector<double> CTiglBSplineAlgorithms::computeParamsBSplineCurve(const Handle(TColgp_HArray1OfPnt)& points, const double alpha)
466
{
467 1
    return computeParamsBSplineCurve(points, 0., 1., alpha);
468
}
469

470 1
std::vector<double> CTiglBSplineAlgorithms::computeParamsBSplineCurve(const Handle(TColgp_HArray1OfPnt)& points, double umin, double umax, const double alpha)
471
{
472 1
    if ( umax <= umin ) {
473 0
        throw tigl::CTiglError("The specified start parameter is larger than the specified end parameter");
474
    }
475
    
476 1
    std::vector<double> parameters(static_cast<size_t>(points->Length()));
477

478 1
    parameters[0] = 0.;
479

480 1
    for (size_t i = 1; i < parameters.size(); ++i) {
481 1
        int iArray = static_cast<int>(i) + points->Lower();
482 1
        double length = pow(points->Value(iArray).SquareDistance(points->Value(iArray - 1)), alpha / 2.);
483 1
        parameters[i] = parameters[i - 1] + length;
484
    }
485

486 1
    double totalLength = parameters.back();
487

488

489 1
    for (size_t i = 0; i < parameters.size(); ++i) {
490 1
        double ratio = 0.;
491 1
        if (totalLength < 1e-10) {
492 1
            ratio = static_cast<double>(i) / static_cast<double>(parameters.size()-1);
493
        }
494
        else {
495 1
            ratio = parameters[i] / totalLength;
496
        }
497 1
        parameters[i] = (umax - umin) * ratio + umin;
498
    }
499

500 1
    return parameters;
501
}
502

503
std::pair<std::vector<double>, std::vector<double> >
504 1
CTiglBSplineAlgorithms::computeParamsBSplineSurf(const TColgp_Array2OfPnt& points, double alpha)
505
{
506
    // first for parameters in u-direction:
507 1
    std::vector<double> paramsU(static_cast<size_t>(points.ColLength()), 0.);
508 1
    for (int vIdx = points.LowerCol(); vIdx <= points.UpperCol(); ++vIdx) {
509 1
        std::vector<double> parameters_u_line = computeParamsBSplineCurve(pntArray2GetColumn(points, vIdx), alpha);
510

511
        // average over columns
512 1
        for (size_t uIdx = 0; uIdx < parameters_u_line.size(); ++uIdx) {
513 1
            paramsU[uIdx] += parameters_u_line[uIdx]/(double)points.RowLength();
514
        }
515
    }
516

517

518
    // now for parameters in v-direction:
519 1
    std::vector<double> paramsV(static_cast<size_t>(points.RowLength()), 0.);
520 1
    for (int uIdx = points.LowerRow(); uIdx <= points.UpperRow(); ++uIdx) {
521 1
        std::vector<double> parameters_v_line = computeParamsBSplineCurve(pntArray2GetRow(points, uIdx), alpha);
522

523
        // average over rows
524 1
        for (size_t vIdx = 0; vIdx < parameters_v_line.size(); ++vIdx) {
525 1
            paramsV[vIdx] += parameters_v_line[vIdx]/(double)points.ColLength();
526
        }
527
    }
528

529
    // put computed parameters for both u- and v-direction in output tuple
530 1
    return std::make_pair(paramsU, paramsV);
531

532
}
533

534

535 1
std::vector<Handle(Geom_BSplineCurve)> CTiglBSplineAlgorithms::createCommonKnotsVectorCurve(const std::vector<Handle(Geom_BSplineCurve)>& splines_vector, double tol)
536
{
537
    // Match parameter range
538 1
    matchParameterRange(splines_vector, tol);
539

540
    // Create a copy that we can modify
541 1
    std::vector<CurveAdapterView> splines_adapter;
542 1
    for (size_t i = 0; i < splines_vector.size(); ++i) {
543 1
        splines_adapter.push_back(Handle(Geom_BSplineCurve)::DownCast(splines_vector[i]->Copy()));
544
    }
545

546 1
    makeGeometryCompatibleImpl(splines_adapter, tol);
547

548 1
    return std::vector<Handle(Geom_BSplineCurve)>(splines_adapter.begin(), splines_adapter.end());
549
}
550

551 1
std::vector<Handle(Geom_BSplineSurface) > CTiglBSplineAlgorithms::createCommonKnotsVectorSurface(const std::vector<Handle(Geom_BSplineSurface) >& old_surfaces_vector, SurfaceDirection dir)
552
{
553
    // Create a copy that we can modify
554 1
    std::vector<SurfAdapterView> adapterSplines;
555 1
    for (size_t i = 0; i < old_surfaces_vector.size(); ++i) {
556 1
        adapterSplines.push_back(SurfAdapterView(Handle(Geom_BSplineSurface)::DownCast(old_surfaces_vector[i]->Copy()), udir));
557
    }
558

559 1
    if (dir == SurfaceDirection::u || dir == SurfaceDirection::both) {
560
        // first in u direction
561 1
        makeGeometryCompatibleImpl(adapterSplines, 1e-14);
562
    }
563

564 1
    if (dir == SurfaceDirection::v || dir == SurfaceDirection::both) {
565
         // now in v direction
566 1
        for (size_t i = 0; i < old_surfaces_vector.size(); ++i) adapterSplines[i].setDir(vdir);
567 1
        makeGeometryCompatibleImpl(adapterSplines, 1e-14);
568
    }
569

570 1
    return std::vector<Handle(Geom_BSplineSurface)>(adapterSplines.begin(), adapterSplines.end());
571
}
572

573 1
void CTiglBSplineAlgorithms::matchParameterRange(std::vector<Handle(Geom_BSplineCurve)> const& bsplines, double tolerance)
574
{
575 1
    Standard_Real umin = bsplines[0]->FirstParameter();
576 1
    Standard_Real umax = bsplines[0]->LastParameter();
577 1
    for (unsigned iP=1; iP<bsplines.size(); ++iP) {
578 1
        Handle(Geom_BSplineCurve) bspl = bsplines[iP];
579 1
        if (fabs(bspl->FirstParameter() - umin) > tolerance ||
580 1
            fabs(bspl->LastParameter() - umax) > tolerance ) {
581 1
            tigl::CTiglBSplineAlgorithms::reparametrizeBSpline(*bspl, umin, umax, tolerance);
582
        }
583
    }
584
}
585

586 1
void CTiglBSplineAlgorithms::matchDegree(const std::vector<Handle(Geom_BSplineCurve) >& bsplines)
587
{
588 1
    int maxDegree = 0;
589 1
    for (std::vector<Handle(Geom_BSplineCurve) >::const_iterator it = bsplines.begin(); it != bsplines.end(); ++it) {
590 1
        int curDegree = (*it)->Degree();
591 1
        if (curDegree > maxDegree) {
592 1
            maxDegree = curDegree;
593
        }
594
    }
595

596 1
    for (std::vector<Handle(Geom_BSplineCurve) >::const_iterator it = bsplines.begin(); it != bsplines.end(); ++it) {
597 1
        int curDegree = (*it)->Degree();
598 1
        if (curDegree < maxDegree) {
599 0
            (*it)->IncreaseDegree(maxDegree);
600
        }
601
    }
602
}
603

604 1
CTiglApproxResult CTiglBSplineAlgorithms::reparametrizeBSplineContinuouslyApprox(const Handle(Geom_BSplineCurve) spline,
605
                                                                                 const std::vector<double>& old_parameters,
606
                                                                                 const std::vector<double>& new_parameters,
607
                                                                                 size_t n_control_pnts)
608
{
609 1
    if (old_parameters.size() != new_parameters.size()) {
610 0
        throw CTiglError("parameter sizes dont match");
611
    }
612

613
    // create a B-spline as a function for reparametrization
614 1
    Handle(TColgp_HArray1OfPnt2d) old_parameters_pnts = new TColgp_HArray1OfPnt2d(1, static_cast<Standard_Integer>(old_parameters.size()));
615 1
    for (size_t parameter_idx = 0; parameter_idx < old_parameters.size(); ++parameter_idx) {
616 1
        int occIdx = static_cast<int>(parameter_idx + 1);
617 1
        old_parameters_pnts->SetValue(occIdx, gp_Pnt2d(old_parameters[parameter_idx], 0));
618
    }
619

620 1
    Geom2dAPI_Interpolate interpolationObject(old_parameters_pnts, OccFArray(new_parameters), false, 1e-15);
621 1
    interpolationObject.Perform();
622

623
    // check that interpolation was successful
624 1
    if (!interpolationObject.IsDone()) {
625 0
        throw CTiglError("Cannot reparametrize", TIGL_MATH_ERROR);
626
    }
627

628 1
    Handle(Geom2d_BSplineCurve) reparametrizing_spline = interpolationObject.Curve();
629

630
    // Create a vector of parameters including the intersection parameters
631 1
    std::vector<double> breaks;
632 1
    for (size_t ipar = 1; ipar < new_parameters.size() - 1; ++ipar) {
633 1
        breaks.push_back(new_parameters[ipar]);
634
    }
635

636 1
    double par_tol = 1e-10;
637

638
#define MODEL_KINKS
639
#ifdef MODEL_KINKS
640
    // remove kinks from breaks
641 1
    std::vector<double> kinks = CTiglBSplineAlgorithms::getKinkParameters(spline);
642
    // convert kink parameters into reparamtetrized parameter using the
643
    // inverse reparametrization function
644 1
    for (size_t ikink = 0; ikink < kinks.size(); ++ikink) {
645 1
        kinks[ikink] = Geom2dAPI_ProjectPointOnCurve(gp_Pnt2d(kinks[ikink], 0.), reparametrizing_spline)
646 1
                           .LowerDistanceParameter();
647
    }
648

649 1
    for (size_t ikink = 0; ikink < kinks.size(); ++ikink) {
650 1
        double kink = kinks[ikink];
651 1
        std::vector<double>::iterator it = std::find_if(breaks.begin(), breaks.end(), IsInsideTolerance(kink, par_tol));
652 1
        if (it != breaks.end()) {
653 1
            breaks.erase(it);
654
        }
655
    }
656
#endif
657

658
    // create equidistance array of parameters, including the breaks
659 1
    std::vector<double> parameters = LinspaceWithBreaks(new_parameters.front(),
660 1
                                                        new_parameters.back(),
661 1
                                                        std::max(static_cast<size_t>(101), n_control_pnts*2),
662 1
                                                        breaks);
663
#ifdef MODEL_KINKS
664
    // insert kinks into parameters array at the correct position
665 1
    for (size_t ikink = 0; ikink < kinks.size(); ++ikink) {
666 1
        double kink = kinks[ikink];
667 1
        parameters.insert( 
668 1
            std::upper_bound( parameters.begin(), parameters.end(), kink),
669 1
            kink);
670
    }
671
#endif
672

673
    // Compute points on spline at the new parameters
674
    // Those will be approximated later on
675 1
    TColgp_Array1OfPnt points(1, static_cast<Standard_Integer>(parameters.size()));
676 1
    for (size_t i = 1; i <= parameters.size(); ++i) {
677 1
        double oldParameter = reparametrizing_spline->Value(parameters[i-1]).X();
678 1
        points(static_cast<Standard_Integer>(i)) = spline->Value(oldParameter);
679
    }
680

681 1
    bool makeContinous = spline->IsClosed() &&
682 1
            spline->DN(spline->FirstParameter(), 1).Angle(spline->DN(spline->LastParameter(), 1)) < 6. / 180. * M_PI;
683

684
    // Create the new spline as a interpolation of the old one
685 1
    CTiglBSplineApproxInterp approximationObj(points, static_cast<int>(n_control_pnts), 3, makeContinous);
686

687 1
    breaks.insert(breaks.begin(), new_parameters.front());
688 1
    breaks.push_back(new_parameters.back());
689
    // Interpolate points at breaking parameters (required for gordon surface)
690 1
    for (size_t ibreak = 0; ibreak < breaks.size(); ++ibreak) {
691 1
        double thebreak = breaks[ibreak];
692
        size_t idx = static_cast<size_t>(
693 1
            std::find_if(parameters.begin(), parameters.end(), IsInsideTolerance(thebreak)) -
694 1
            parameters.begin());
695 1
        approximationObj.InterpolatePoint(idx);
696
    }
697

698
#ifdef MODEL_KINKS
699 1
    for (size_t ikink = 0; ikink < kinks.size(); ++ikink) {
700 1
        double kink = kinks[ikink];
701
        size_t idx = static_cast<size_t>(
702 1
            std::find_if(parameters.begin(), parameters.end(), IsInsideTolerance(kink, par_tol)) -
703 1
            parameters.begin());
704 1
        approximationObj.InterpolatePoint(idx, true);
705
    }
706
#endif
707

708 1
    CTiglApproxResult result = approximationObj.FitCurveOptimal(parameters);
709

710 1
    assert(!result.curve.IsNull());
711

712 1
    return result;
713
}
714

715 1
Handle(Geom_BSplineSurface) CTiglBSplineAlgorithms::flipSurface(const Handle(Geom_BSplineSurface) surface)
716
{
717 1
    Handle(Geom_BSplineSurface) result = Handle(Geom_BSplineSurface)::DownCast(surface->Copy());
718 1
    result->ExchangeUV();
719 1
    return result;
720
}
721

722 1
Handle(Geom_BSplineSurface) CTiglBSplineAlgorithms::pointsToSurface(const TColgp_Array2OfPnt& points,
723
                                                                    const std::vector<double>& uParams,
724
                                                                    const std::vector<double>& vParams,
725
                                                                    bool uContinousIfClosed, bool vContinousIfClosed)
726
{
727

728 1
    double tolerance = REL_TOL_CLOSED * scale(points);
729 1
    bool makeVDirClosed = vContinousIfClosed & isVDirClosed(points, tolerance);
730 1
    bool makeUDirClosed = uContinousIfClosed & isUDirClosed(points, tolerance);
731

732
    // first interpolate all points by B-splines in u-direction
733 1
    std::vector<Handle(Geom_Curve)> uSplines;
734 1
    for (int cpVIdx = points.LowerCol(); cpVIdx <= points.UpperCol(); ++cpVIdx) {
735 1
        Handle_TColgp_HArray1OfPnt points_u = pntArray2GetColumn(points, cpVIdx);
736 1
        CTiglPointsToBSplineInterpolation interpolationObject(points_u, uParams, 3, makeUDirClosed);
737

738 1
        Handle(Geom_Curve) curve = interpolationObject.Curve();
739 1
        uSplines.push_back(curve);
740
    }
741

742
    // now create a skinned surface with these B-splines which represents the interpolating surface
743 1
    CTiglCurvesToSurface skinner(uSplines, vParams, makeVDirClosed );
744 1
    Handle(Geom_BSplineSurface) interpolatingSurf = skinner.Surface();
745

746 1
    return interpolatingSurf;
747
}
748

749

750 1
std::vector<std::pair<double, double> > CTiglBSplineAlgorithms::intersections(const Handle(Geom_BSplineCurve) spline1, const Handle(Geom_BSplineCurve) spline2, double tolerance) {
751
    // light weight simple minimizer
752

753
    // check parametrization of B-splines beforehand
754

755
    // find out the average scale of the two B-splines in order to being able to handle a more approximate curves and find its intersections
756 1
    double splines_scale = (CTiglBSplineAlgorithms::scale(spline1) + CTiglBSplineAlgorithms::scale(spline2)) / 2.;
757

758 1
    std::vector<std::pair<double, double> > intersection_params_vector;
759 1
    GeomAPI_ExtremaCurveCurve intersectionObj(spline1, spline2);
760 1
    for (int intersect_idx = 1; intersect_idx <= intersectionObj.NbExtrema(); ++intersect_idx) {
761 1
        double param1 = 0.;
762 1
        double param2 = 0.;
763 1
        intersectionObj.Parameters(intersect_idx, param1, param2);
764

765
        // filter out real intersections
766 1
        gp_Pnt point1 = spline1->Value(param1);
767 1
        gp_Pnt point2 = spline2->Value(param2);
768

769 1
        if (point1.Distance(point2) < tolerance * splines_scale) {
770 1
            intersection_params_vector.push_back(std::make_pair(param1, param2));
771
        }
772
        else {
773 0
            throw CTiglError("Curves do not intersect each other", TIGL_MATH_ERROR);
774
        }
775

776
        // for closed B-splines:
777 1
        if (intersectionObj.NbExtrema() == 1 && spline1->IsClosed() && std::abs(param1 - spline1->Knot(1)) < 1e-6) {
778
            // GeomAPI_ExtremaCurveCurve doesn't find second intersection point at the end of the closed curve, so add it by hand
779 1
            intersection_params_vector.push_back(std::make_pair(spline1->Knot(spline1->NbKnots()), param2));
780
        }
781

782 1
        if (intersectionObj.NbExtrema() == 1 && spline1->IsClosed() && std::abs(param1 - spline1->Knot(spline1->NbKnots())) < 1e-6) {
783
            // GeomAPI_ExtremaCurveCurve doesn't find second intersection point at the beginning of the closed curve, so add it by hand
784 1
            intersection_params_vector.push_back(std::make_pair(spline1->Knot(1), param2));
785
        }
786

787 1
        if (intersectionObj.NbExtrema() == 1 && spline2->IsClosed() && std::abs(param2 - spline2->Knot(1)) < 1e-6) {
788
            // GeomAPI_ExtremaCurveCurve doesn't find second intersection point at the end of the closed curve, so add it by hand
789 0
            intersection_params_vector.push_back(std::make_pair(param1, spline2->Knot(spline2->NbKnots())));
790
        }
791

792 1
        if (intersectionObj.NbExtrema() == 1 && spline2->IsClosed() && std::abs(param2 - spline2->Knot(spline2->NbKnots())) < 1e-6) {
793
            // GeomAPI_ExtremaCurveCurve doesn't find second intersection point at the beginning of the closed curve, so add it by hand
794 0
            intersection_params_vector.push_back(std::make_pair(param1, spline2->Knot(1)));
795
        }
796
    }
797

798

799 1
    return intersection_params_vector;
800
}
801

802 1
double CTiglBSplineAlgorithms::scale(const std::vector<Handle(Geom_BSplineCurve)>& splines_vector)
803
{
804 1
    double maxScale = 0.;
805 1
    for (std::vector<Handle(Geom_BSplineCurve)>::const_iterator it = splines_vector.begin(); it != splines_vector.end(); ++it) {
806 1
        maxScale = std::max(scale(*it), maxScale);
807
    }
808

809 1
    return maxScale;
810
}
811

812 1
double CTiglBSplineAlgorithms::scale(const Handle(Geom_BSplineCurve)& spline)
813
{
814 1
    double scale = 0.;
815 1
    gp_Pnt first_ctrl_pnt = spline->Pole(1);
816 1
    for (int ctrl_pnt_idx = 2; ctrl_pnt_idx <= spline->NbPoles(); ++ctrl_pnt_idx) {
817
        // compute distance of the first control point to the others and save biggest distance
818 1
        double distance = first_ctrl_pnt.Distance(spline->Pole(ctrl_pnt_idx));
819

820 1
        scale = std::max(scale, distance);
821
    }
822 1
    return scale;
823
}
824

825 1
void CTiglBSplineAlgorithms::reparametrizeBSpline(Geom_BSplineCurve& spline, double umin, double umax, double tol)
826
{
827 1
    if (std::abs(spline.Knot(1) - umin) > tol || std::abs(spline.Knot(spline.NbKnots()) - umax) > tol) {
828 1
        TColStd_Array1OfReal aKnots (1, spline.NbKnots());
829 1
        spline.Knots (aKnots);
830 1
        BSplCLib::Reparametrize (umin, umax, aKnots);
831 1
        spline.SetKnots (aKnots);
832
    }
833
}
834

835 1
void CTiglBSplineAlgorithms::reparametrizeBSpline(Geom_BSplineSurface& spline, double umin, double umax, double vmin, double vmax, double tol)
836
{
837 1
    if (std::abs(spline.UKnot(1) - umin) > tol || std::abs(spline.UKnot(spline.NbUKnots()) - umax) > tol) {
838 1
        TColStd_Array1OfReal aKnots (1, spline.NbUKnots());
839 1
        spline.UKnots (aKnots);
840 1
        BSplCLib::Reparametrize (umin, umax, aKnots);
841 1
        spline.SetUKnots (aKnots);
842
    }
843

844 1
    if (std::abs(spline.VKnot(1) - vmin) > tol || std::abs(spline.VKnot(spline.NbVKnots()) - vmax) > tol) {
845 1
        TColStd_Array1OfReal aKnots (1, spline.NbVKnots());
846 1
        spline.VKnots (aKnots);
847 1
        BSplCLib::Reparametrize (vmin, vmax, aKnots);
848 1
        spline.SetVKnots (aKnots);
849
    }
850
}
851

852 1
CTiglApproxResult CTiglBSplineAlgorithms::reparametrizeBSplineNiceKnots(Handle(Geom_BSplineCurve) spline)
853
{
854 1
    if (spline.IsNull()) {
855 1
        throw CTiglError("Null Pointer curve", TIGL_NULL_POINTER);
856
    }
857

858 1
    double umin = spline->FirstParameter();
859 1
    double umax = spline->LastParameter();
860

861
    // compute the number of desired knot segments
862

863 1
    int nSegmentsOld = spline->NbPoles() - spline->Degree();
864
    // The new number of segments is a power of two
865 1
    int nSegmentsNew = pow(2.0, static_cast<int>(log2(static_cast<float>(nSegmentsOld))));
866

867
    // we are using at least 8 segments to be safe with accuracy
868 1
    nSegmentsNew = std::max(nSegmentsNew, 8);
869

870
    // The reparametrization results in a degree 3 spline
871 1
    int target_degree = 3;
872 1
    return CTiglBSplineAlgorithms::reparametrizeBSplineContinuouslyApprox(spline, {umin, umax}, {umin, umax}, nSegmentsNew + target_degree);
873
}
874

875 1
math_Matrix CTiglBSplineAlgorithms::bsplineBasisMat(int degree, const TColStd_Array1OfReal& knots, const TColStd_Array1OfReal& params, unsigned int derivOrder)
876
{
877 1
    Standard_Integer ncp = knots.Length() - degree - 1;
878 1
    math_Matrix mx(1, params.Length(), 1, ncp);
879 1
    mx.Init(0.);
880 1
    math_Matrix bspl_basis(1, derivOrder + 1, 1, degree + 1);
881 1
    bspl_basis.Init(0.);
882 1
    for (Standard_Integer iparm = 1; iparm <= params.Length(); ++iparm) {
883 1
        Standard_Integer basis_start_index = 0;
884
#if OCC_VERSION_HEX >= VERSION_HEX_CODE(7,1,0)
885
        BSplCLib::EvalBsplineBasis(derivOrder, degree + 1, knots, params.Value(iparm), basis_start_index, bspl_basis);
886
#else
887 1
        BSplCLib::EvalBsplineBasis(1, derivOrder, degree + 1, knots, params.Value(iparm), basis_start_index, bspl_basis);
888
#endif
889 1
        if(derivOrder > 0) {
890 1
            math_Vector help_vector(1, ncp);
891 1
            help_vector.Init(0.);
892 1
            help_vector.Set(basis_start_index, basis_start_index + degree, bspl_basis.Row(derivOrder + 1));
893 1
            mx.SetRow(iparm, help_vector);
894
        }
895
        else {
896 1
            mx.Set(iparm, iparm, basis_start_index, basis_start_index + degree, bspl_basis);
897
        }
898
    }
899 1
    return mx;
900
}
901

902 1
std::vector<double> CTiglBSplineAlgorithms::getKinkParameters(const Handle(Geom_BSplineCurve)& curve)
903
{
904 1
    if (curve.IsNull()) {
905 0
        throw CTiglError("Null Pointer curve", TIGL_NULL_POINTER);
906
    }
907

908 1
    double eps = 1e-8;
909

910 1
    std::vector<double> kinks;
911 1
    for (int knotIndex = 2; knotIndex < curve->NbKnots(); ++knotIndex) {
912 1
        if (curve->Multiplicity(knotIndex) == curve->Degree()) {
913 1
            double knot = curve->Knot(knotIndex);
914
            // check if really a kink
915 1
            double angle = curve->DN(knot + eps, 1).Angle(curve->DN(knot - eps, 1));
916 1
            if (angle > 6./180. * M_PI) {
917 1
                kinks.push_back(knot);
918
            }
919
        }
920
    }
921

922 1
    return kinks;
923
}
924

925 1
CTiglBSplineAlgorithms::SurfaceKinks CTiglBSplineAlgorithms::getKinkParameters(const Handle(Geom_BSplineSurface)& surface)
926
{
927 1
    if (surface.IsNull()) {
928 0
        throw CTiglError("Null Pointer curve", TIGL_NULL_POINTER);
929
    }
930

931 1
    SurfaceKinks kinks;
932

933 1
    for (int knotIndex = 2; knotIndex < surface->NbUKnots(); ++knotIndex) {
934 1
        if (surface->UMultiplicity(knotIndex) == surface->UDegree()) {
935 0
            double knot = surface->UKnot(knotIndex);
936 0
            kinks.u.push_back(knot);
937
        }
938
    }
939

940 1
    for (int knotIndex = 2; knotIndex < surface->NbVKnots(); ++knotIndex) {
941 1
        if (surface->VMultiplicity(knotIndex) == surface->VDegree()) {
942 0
            double knot = surface->VKnot(knotIndex);
943 0
            kinks.v.push_back(knot);
944
        }
945
    }
946

947 1
    return kinks;
948
}
949

950 1
Handle(Geom_BSplineSurface) CTiglBSplineAlgorithms::trimSurface(const Handle(Geom_Surface)& surface, double umin, double umax, double vmin, double vmax)
951
{
952

953 1
    Handle(Geom_BSplineSurface) trimmedSurface = GeomConvert::SurfaceToBSplineSurface(surface);
954

955 1
    trimmedSurface->SetUNotPeriodic();
956 1
    trimmedSurface->SetVNotPeriodic();
957

958
    // workaround for OCCT bug https://tracker.dev.opencascade.org/view.php?id=31402
959
    // We set the trimming parameters to a knot, if they are close to them
960 1
    double parTol = 1e-10;
961
    int i1, i2;
962 1
    trimmedSurface->LocateU(umin, parTol, i1, i2);
963 1
    if (i1 == i2) {
964
        // v is a knot
965 1
        umin = trimmedSurface->UKnot(i1);
966
    }
967 1
    trimmedSurface->LocateU(umax, parTol, i1, i2);
968 1
    if (i1 == i2) {
969
        // v is a knot
970 1
        umax = trimmedSurface->UKnot(i1);
971
    }
972 1
    trimmedSurface->LocateV(vmin, parTol, i1, i2);
973 1
    if (i1 == i2) {
974
        // v is a knot
975 1
        vmin = trimmedSurface->VKnot(i1);
976
    }
977 1
    trimmedSurface->LocateV(vmax, parTol, i1, i2);
978 1
    if (i1 == i2) {
979
        // v is a knot
980 1
        vmax = trimmedSurface->VKnot(i1);
981
    }
982

983
    // Perform the trimming
984 1
    trimmedSurface->Segment(umin, umax, vmin, vmax);
985

986
#ifdef DEBUG
987
    double u1, u2, v1, v2;
988 1
    trimmedSurface->Bounds(u1, u2, v1, v2);
989

990 1
    double tol_check = 1e-12;
991 1
    auto intol = [tol_check](double a, double b) {
992 1
        return fabs(a - b) <= tol_check;
993 1
    };
994

995 1
    assert(intol(u1, umin) && intol(u2, umax) && intol(v1, vmin) && intol(v2, vmax));
996
#endif
997 1
    return trimmedSurface;
998
}
999

1000 1
Handle(Geom_BSplineCurve) CTiglBSplineAlgorithms::trimCurve(const Handle(Geom_BSplineCurve)& curve, double umin, double umax)
1001
{
1002 1
    Handle(Geom_BSplineCurve) copy = Handle(Geom_BSplineCurve)::DownCast(curve->Copy());
1003 1
    copy->Segment(umin, umax);
1004 1
    return copy;
1005
}
1006

1007 1
Handle(Geom_BSplineCurve) CTiglBSplineAlgorithms::concatCurves(std::vector<Handle(Geom_BSplineCurve)> curves,
1008
                                                               bool parByLength, double tolerance)
1009
{
1010

1011 1
    if (curves.size() == 0) {
1012 0
        LOG(ERROR) << "Empty curve vector in CTiglBSplineAlgorithms::concatCurves";
1013 0
        return nullptr;
1014
    }
1015 1
    else if (curves.size() == 1) {
1016 1
        return curves[0];
1017
    }
1018

1019 1
    std::vector<double> lengths;
1020 1
    double totalLen  = 0;
1021 1
    int    maxDegree = 0;
1022

1023
    // get the bsplines of each edge,
1024
    // compute the lengths of each edge,
1025
    // determine maximum degree of the curves
1026 1
    for (auto curve : curves) {
1027

1028 1
        if (parByLength) {
1029
            // find out length of current curve
1030 1
            Standard_Real umin = curve->FirstParameter();
1031 1
            Standard_Real umax = curve->LastParameter();
1032 1
            GeomAdaptor_Curve adaptorCurve(curve, umin, umax);
1033 1
            double len = GCPnts_AbscissaPoint::Length(adaptorCurve, umin, umax);
1034 1
            lengths.push_back(len);
1035 1
            totalLen += len;
1036
        }
1037

1038
        // find out maximum degree
1039 1
        if (curve->Degree() > maxDegree) {
1040 1
            maxDegree = curve->Degree();
1041
        }
1042
    }
1043

1044

1045
    // check connectivities
1046 1
    for (unsigned int icurve = 1; icurve < curves.size(); ++icurve) {
1047 1
        Handle(Geom_BSplineCurve) c1 = curves[icurve-1];
1048 1
        Handle(Geom_BSplineCurve) c2 = curves[icurve];
1049

1050 1
        gp_Pnt p1 = c1->Pole(c1->NbPoles());
1051 1
        gp_Pnt p2 = c2->Pole(1);
1052

1053 1
        if (p1.Distance(p2) > tolerance) {
1054
            // error
1055 0
            LOG(ERROR) << "Curves not connected within tolerance in concatCurves";
1056 0
            return nullptr;
1057
        }
1058
    }
1059

1060
    // elevate degree of all curves to maxDegree
1061 1
    for (unsigned int icurve = 0; icurve < curves.size(); ++icurve) {
1062 1
        Handle(Geom_BSplineCurve) curve = curves[icurve];
1063 1
        curve->IncreaseDegree(maxDegree);
1064
    }
1065

1066
#ifdef DEBUG
1067
    // check that each curve is at maxDegree
1068 1
    for (unsigned int icurve = 0; icurve < curves.size(); ++icurve) {
1069 1
        Handle(Geom_BSplineCurve) curve = curves[icurve];
1070 1
        assert(curve->Degree() == maxDegree);
1071
    }
1072
#endif
1073

1074
    // shift knots of curves
1075 1
    double startPar = 0;
1076 1
    for (unsigned int icurve = 0; icurve < curves.size(); ++icurve) {
1077 1
        Handle(Geom_BSplineCurve) curve = curves[icurve];
1078 1
        double stopPar = startPar;
1079 1
        if (parByLength) {
1080 1
            stopPar += lengths[icurve]/totalLen;
1081
        }
1082
        else {
1083 1
            stopPar += curve->LastParameter() - curve->FirstParameter();
1084
        }
1085 1
        CTiglBSplineAlgorithms::reparametrizeBSpline(*curve, startPar, stopPar);
1086 1
        curves[icurve] = curve;
1087

1088 1
        startPar = stopPar;
1089
    }
1090

1091
    // count number of knots and control points for the final b-spline
1092 1
    int nbknots = 1;
1093 1
    int nbcp    = 1;
1094 1
    for (unsigned int icurve = 0; icurve < curves.size(); ++icurve) {
1095 1
        Handle(Geom_BSplineCurve) curve = curves[icurve];
1096 1
        nbknots += curve->NbKnots() - 1;
1097 1
        nbcp    += curve->NbPoles() - 1;
1098
    }
1099

1100
    // allocate arrays
1101 1
    TColgp_Array1OfPnt      cpoints(1, nbcp);
1102 1
    TColStd_Array1OfReal    weights(1, nbcp);
1103 1
    TColStd_Array1OfReal    knots(1, nbknots);
1104 1
    TColStd_Array1OfInteger mults(1, nbknots);
1105

1106
    // concatenate everything
1107 1
    int iknotT = 1, imultT = 1, icpT = 1, iweightT = 1;
1108 1
    for (unsigned int icurve = 0; icurve < curves.size(); ++icurve) {
1109 1
        Handle(Geom_BSplineCurve) curve = curves[icurve];
1110

1111
        // special handling of the first knot, control point
1112 1
        knots.SetValue(iknotT++, curve->Knot(1));
1113 1
        if (icurve == 0) {
1114
            // we just copy the data of the very first point/knot
1115 1
            mults.SetValue(imultT++, curve->Multiplicity(1));
1116 1
            cpoints.SetValue(icpT++, curve->Pole(1));
1117 1
            weights.SetValue(iweightT++, curve->Weight(1));
1118
        }
1119
        else {
1120
            // set multiplicity to maxDegree to allow c0 concatenation
1121 1
            mults.SetValue(imultT++, maxDegree);
1122

1123
            // compute midpoint between endpoint of previous
1124
            // curve and startpoint of current curve
1125 1
            Handle(Geom_BSplineCurve) lastCurve = curves[icurve-1];
1126 1
            gp_Pnt endPoint   = lastCurve->Pole(lastCurve->NbPoles());
1127 1
            gp_Pnt startPoint = curve->Pole(1);
1128 1
            gp_Pnt midPoint = (endPoint.XYZ() + startPoint.XYZ())/2.;
1129 1
            cpoints.SetValue(icpT++, midPoint);
1130

1131
            // we use the average weight of previous curve and current curve
1132
            // This is probably wrong and could change the shape of the curve
1133
            // Instead, one could scale all weights of the curve to match the weight of
1134
            // the previous curve
1135 1
            weights.SetValue(iweightT++, (lastCurve->Weight(lastCurve->NbPoles()) + curve->Weight(1))/2.);
1136
        }
1137

1138
        // just copy control points, weights, knots and multiplicites
1139 1
        for (int iknot = 2; iknot < curve->NbKnots(); ++iknot) {
1140 1
            knots.SetValue(iknotT++, curve->Knot(iknot));
1141 1
            mults.SetValue(imultT++, curve->Multiplicity(iknot));
1142
        }
1143 1
        for (int icp = 2; icp < curve->NbPoles(); ++icp) {
1144 1
            cpoints.SetValue(icpT++, curve->Pole(icp));
1145 1
            weights.SetValue(iweightT++, curve->Weight(icp));
1146
        }
1147

1148
    }
1149

1150
    // special handling of the last point and knot
1151 1
    Handle(Geom_BSplineCurve) lastCurve = curves[curves.size()-1];
1152 1
    knots.SetValue(iknotT, lastCurve->Knot(lastCurve->NbKnots()));
1153 1
    mults.SetValue(imultT,  lastCurve->Multiplicity(lastCurve->NbKnots()));
1154 1
    cpoints.SetValue(icpT, lastCurve->Pole(lastCurve->NbPoles()));
1155 1
    weights.SetValue(iweightT, lastCurve->Weight(lastCurve->NbPoles()));
1156

1157
#ifdef DEBUG
1158
    // check that we have the correct number of knots, control points etc...
1159 1
    int nkn = 0;
1160 1
    for (int ik = knots.Lower(); ik <= knots.Upper(); ++ik) {
1161 1
        nkn += mults.Value(ik);
1162
    }
1163

1164
    // check validity of bspline
1165 1
    assert (cpoints.Length() + maxDegree + 1 == nkn);
1166
#endif
1167

1168
    // build the resulting B-Spline
1169 1
    Handle(Geom_BSplineCurve) result = new Geom_BSplineCurve(cpoints, weights, knots, mults, maxDegree, false);
1170 1
    return result;
1171
}
1172

1173 1
Handle(Geom_BSplineSurface) CTiglBSplineAlgorithms::concatSurfacesUDir(Handle(Geom_BSplineSurface) bspl1,
1174
                                                                       Handle(Geom_BSplineSurface) bspl2,
1175
                                                                       double space_tol)
1176
{
1177

1178 1
    DEBUG_SCOPE(debug);
1179 1
    debug.addShape(BRepBuilderAPI_MakeFace(bspl1, 1e-6).Face(), "surf1");
1180 1
    debug.addShape(BRepBuilderAPI_MakeFace(bspl2, 1e-6).Face(), "surf2");
1181

1182
    // the surfaces must not be periodic in u direction
1183 1
    assert (bspl1->IsUPeriodic() == false);
1184 1
    assert (bspl2->IsUPeriodic() == false);
1185

1186
    // we dont have nurbs implemented right now
1187 1
    assert (bspl1->IsURational() == false);
1188 1
    assert (bspl2->IsURational() == false);
1189 1
    assert (bspl1->IsVRational() == false);
1190 1
    assert (bspl2->IsVRational() == false);
1191

1192 1
    double param_tol = 1e-14;
1193

1194
    // check that surfaces are following parametrically
1195
    double umin1, umax1, vmin1, vmax1;
1196 1
    bspl1->Bounds(umin1, umax1, vmin1, vmax1);
1197
    double umin2, umax2, vmin2, vmax2;
1198 1
    bspl2->Bounds(umin2, umax2, vmin2, vmax2);
1199

1200 1
    if (std::abs(umax1 - umin2) > param_tol) {
1201 1
        std::stringstream str;
1202 1
        str << "Surfaces do not follow in u-parametric direction in CTiglBSplineAlgorithms::concatSurfacesUDir. ";
1203 1
        str << "Surface 1 ends at " << umax1 << ", Surface 2 begins at " << umin2 << "!";
1204 1
        throw CTiglError(str.str(), TIGL_MATH_ERROR);
1205
    }
1206

1207 1
    bspl1->SetVNotPeriodic();
1208 1
    bspl2->SetVNotPeriodic();
1209

1210 1
    auto u_degree = std::max(bspl1->UDegree(), bspl2->UDegree());
1211 1
    auto v_degree = std::max(bspl1->VDegree(), bspl2->VDegree());
1212

1213 1
    bspl1->IncreaseDegree(u_degree, v_degree);
1214 1
    bspl2->IncreaseDegree(u_degree, v_degree);
1215

1216
    // check that corner points match
1217 1
    auto leftCornerDist = bspl1->Value(umax1, vmin1).Distance(bspl2->Value(umin2, vmin2));
1218 1
    auto rightCornerDist = bspl1->Value(umax1, vmax1).Distance(bspl2->Value(umin2, vmax2)) ;
1219 1
    if (leftCornerDist > space_tol || rightCornerDist > space_tol) {
1220 1
        throw CTiglError("Surfaces don't match within tolerances in CTiglBSplineAlgorithms::concatSurfacesUDir.", TIGL_MATH_ERROR);
1221
    }
1222

1223 1
    auto spl_vec = CTiglBSplineAlgorithms::createCommonKnotsVectorSurface({bspl1, bspl2}, SurfaceDirection::v);
1224 1
    bspl1 = spl_vec[0];
1225 1
    bspl2 = spl_vec[1];
1226

1227 1
    assert(bspl1->NbVPoles() == bspl2->NbVPoles());
1228

1229
    // concat control points
1230 1
    TColgp_Array2OfPnt cpsNew(1, bspl1->NbUPoles() + bspl2->NbUPoles() - 1, 1, bspl1->NbVPoles());
1231 1
    for (int iv = 0; iv < bspl1->NbVPoles(); ++iv) {
1232 1
        for (int iu = 0; iu < bspl1->NbUPoles() - 1; ++iu) {
1233 1
            cpsNew.SetValue(iu+1, iv+1, bspl1->Pole(iu+1, iv+1));
1234
        }
1235 1
        auto offset = bspl1->NbUPoles();
1236 1
        for (int iu = 0; iu < bspl2->NbUPoles(); ++iu) {
1237 1
            cpsNew.SetValue(iu + offset, iv+1, bspl2->Pole(iu+1, iv+1));
1238
        }
1239
    }
1240

1241
    // concat knots and mults
1242 1
    TColStd_Array1OfReal uknots_new(1, bspl1->NbUKnots() + bspl2->NbUKnots() - 1);
1243 1
    TColStd_Array1OfInteger umults_new(1, bspl1->NbUKnots() + bspl2->NbUKnots() - 1);
1244

1245 1
    for (int i = 0; i < bspl1->NbUKnots()-1; ++i) {
1246 1
        uknots_new.SetValue(i+1, bspl1->UKnot(i+1));
1247 1
        umults_new.SetValue(i+1, bspl1->UMultiplicity(i+1));
1248
    }
1249

1250 1
    int middleIdx = bspl1->NbUKnots();
1251 1
    uknots_new.SetValue(middleIdx, 0.5 * (bspl1->UKnot(middleIdx) + bspl2->UKnot(1)));
1252 1
    umults_new.SetValue(middleIdx, bspl1->UMultiplicity(middleIdx) - 1);
1253

1254 1
    for (int i = 1; i < bspl2->NbUKnots(); ++i) {
1255 1
        uknots_new.SetValue(i + middleIdx, bspl2->UKnot(i+1));
1256 1
        umults_new.SetValue(i + middleIdx, bspl2->UMultiplicity(i+1));
1257
    }
1258

1259
    // we simply take the v-knots of the first spline
1260 1
    TColStd_Array1OfReal vknots_new(1, bspl1->NbVKnots());
1261 1
    TColStd_Array1OfInteger vmults_new(1, bspl1->NbVKnots());
1262 1
    bspl1->VKnots(vknots_new);
1263 1
    bspl1->VMultiplicities(vmults_new);
1264

1265
#ifdef DEBUG
1266 1
    int sum_umults = 0;
1267 1
    for (int i = umults_new.Lower(); i <= umults_new.Upper(); ++i) {
1268 1
        sum_umults += umults_new.Value(i);
1269
    }
1270

1271 1
    int sum_vmults = 0;
1272 1
    for (int i = vmults_new.Lower(); i <= vmults_new.Upper(); ++i) {
1273 1
        sum_vmults += vmults_new.Value(i);
1274
    }
1275

1276 1
    assert(cpsNew.ColLength() + u_degree + 1 == sum_umults);
1277 1
    assert(cpsNew.RowLength() + v_degree + 1 == sum_vmults);
1278
#endif
1279

1280
    return new Geom_BSplineSurface(cpsNew, uknots_new, vknots_new,
1281
                                   umults_new, vmults_new,
1282 1
                                   u_degree, v_degree, false, false);
1283
}
1284

1285 1
Handle(Geom_BSplineSurface) CTiglBSplineAlgorithms::makeKnotsUniform(Handle(Geom_BSplineSurface) surf, unsigned int nseg_u, unsigned int nseg_v)
1286
{
1287
    double u1, u2, v1, v2;
1288 1
    surf->Bounds(u1, u2, v1, v2);
1289

1290
    // we create a degree 3 approximation in all directions
1291 1
    auto ncp_u = std::max(4, static_cast<int>(nseg_u + 1));
1292 1
    auto ncp_v = std::max(4, static_cast<int>(nseg_v + 1));
1293

1294 1
    auto u = LinspaceWithBreaks(u1, u2, ncp_u, {});
1295 1
    auto v = LinspaceWithBreaks(v1, v2, ncp_v, {});
1296

1297 1
    TColgp_Array2OfPnt points(1, ncp_u, 1, ncp_v);
1298

1299 1
    for (size_t i = 0; i < u.size(); ++i) {
1300 1
        for (size_t j = 0; j < v.size(); ++j) {
1301 1
            auto pnt = surf->Value(u[i], v[j]);
1302 1
            points.SetValue(static_cast<Standard_Integer>(i+1), static_cast<Standard_Integer>(j+1), pnt);
1303
        }
1304
    }
1305

1306 1
    return CTiglBSplineAlgorithms::pointsToSurface(points, u, v, true, true);
1307
}
1308

1309
} // namespace tigl

Read our documentation on viewing source code .

Loading