FreeCAD

1
// This file is released under the BSD license
2
//
3
// Copyright (c) 2009, Daniel Heeks
4
// All rights reserved.
5
//
6
// Redistribution and use in source and binary forms, with or without modification,
7
// are permitted provided that the following conditions are met:
8
//
9
//    * Redistributions of source code must retain the above copyright notice, this
10
//      list of conditions and the following disclaimer.
11
//    * Redistributions in binary form must reproduce the above copyright notice, this
12
//      list of conditions and the following disclaimer in the documentation and/or
13
//      other materials provided with the distribution.
14
//    * Neither the name of Daniel Heeks nor the names of its contributors may be used
15
//      to endorse or promote products derived from this software without specific prior
16
//      written permission.
17
//
18
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
19
// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
20
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
21
// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
22
// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
23
// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
24
// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
25
// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
26
// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
27
// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
28
// POSSIBILITY OF SUCH DAMAGE.
29

30
#include "PreCompiled.h"
31
#ifndef _PreComp_
32
# include <GC_MakeArcOfCircle.hxx>
33
# include <gp_Ax2.hxx>
34
# include <gp_Circ.hxx>
35
# include <gp_Pln.hxx>
36
#endif
37

38
#include "BSplineCurveBiArcs.h"
39
#include "Geometry.h"
40
#include "Tools.h"
41

42

43
using Part::BSplineCurveBiArcs;
44
using Part::GeomBSplineCurve;
45
using Part::Geometry;
46

47
// Algorithm taken from HeeksCAD
48
namespace Part {
49

50
bool tangentialArc(const gp_Pnt& p0, const gp_Vec& v0, const gp_Pnt& p1, gp_Pnt& c, gp_Dir& axis)
51
{
52
    if (p0.Distance(p1) > Precision::Intersection() &&
53
        v0.Magnitude() > Precision::Intersection()){
54
        gp_Vec v1(p0, p1);
55
        gp_Pnt halfway(p0.XYZ() + v1.XYZ() * 0.5);
56
        gp_Pln pln1(halfway, v1);
57
        gp_Pln pln2(p0, v0);
58
        gp_Lin plane_line;
59
        if (intersect(pln1, pln2, plane_line)) {
60
            gp_Lin l1(halfway, v1);
61
            gp_Pnt p2;
62
            closestPointsOnLines(plane_line, l1, c, p2);
63
            axis = -(plane_line.Direction());
64
            return true;
65
        }
66
    }
67

68
    return false;
69
}
70

71
class TangentialArc
72
{
73
public:
74
    gp_Pnt m_p0; // start point
75
    gp_Vec m_v0; // start direction
76
    gp_Pnt m_p1; // end point
77
    gp_Pnt m_c; // centre point
78
    gp_Dir m_a; // axis
79
    bool m_is_a_line;
80
    TangentialArc(const gp_Pnt& p0, const gp_Vec& v0, const gp_Pnt& p1)
81
        : m_p0(p0), m_v0(v0), m_p1(p1)
82
    {
83
        // calculate a tangential arc that goes through p0 and p1, with a direction of v0 at p0
84
        m_is_a_line = !Part::tangentialArc(m_p0, m_v0, m_p1, m_c, m_a);
85
    }
86
    bool isRadiusEqual(const gp_Pnt &p, double tolerance) const
87
    {
88
        if (m_is_a_line)
89
            return true;
90

91
        double point_radius = gp_Vec(m_c.XYZ() - p.XYZ()).Magnitude();
92
        double diff =  fabs(point_radius - radius());
93
        return diff <= tolerance;
94
    }
95
    double radius() const
96
    {
97
        double r0 = gp_Vec(m_p0.XYZ() - m_c.XYZ()).Magnitude();
98
        double r1 = gp_Vec(m_p1.XYZ() - m_c.XYZ()).Magnitude();
99
        double r = (r0 + r1)/2;
100
        return r;
101
    }
102
    Geometry* makeArc() const
103
    {
104
        if (m_is_a_line) {
105
            GeomLineSegment* line = new GeomLineSegment();
106
            line->setPoints(Base::convertTo<Base::Vector3d>(m_p0),Base::convertTo<Base::Vector3d>(m_p1));
107
            return line;
108
        }
109

110
        gp_Circ c(gp_Ax2(m_c, m_a), radius());
111
        GC_MakeArcOfCircle arc(c, m_p0, m_p1, true);
112
        GeomArcOfCircle* new_object = new GeomArcOfCircle();
113
        new_object->setHandle(arc.Value());
114
        return new_object;
115
    }
116
};
117

118
}
119

120
void BSplineCurveBiArcs::createArcs(double tolerance, std::list<Geometry*>& new_spans,
121
                                    const gp_Pnt& p_start, const gp_Vec& v_start,
122
                                    double t_start, double t_end, gp_Pnt& p_end, gp_Vec& v_end) const
123
{
124
    this->myCurve->D1(t_end, p_end, v_end);
125

126
    gp_Pnt p1, p2, p3;
127

128
    Type can_do_spline_whole = calculateBiArcPoints(t_start, p_start, v_start, t_end, p_end, v_end, p1, p2, p3);
129

130
    Geometry* arc_object1 = nullptr;
131
    Geometry* arc_object2 = nullptr;
132

133
    if (can_do_spline_whole == Type::SingleArc) {
134
        Part::TangentialArc arc1(p_start, v_start, p2);
135
        Part::TangentialArc arc2(p2, gp_Vec(p3.XYZ() - p2.XYZ()), p_end);
136

137
        gp_Pnt p_middle1, p_middle2;
138
        this->myCurve->D0(t_start + ((t_end - t_start) * 0.25), p_middle1);
139
        this->myCurve->D0(t_start + ((t_end - t_start) * 0.75), p_middle2);
140

141
        if (!arc1.isRadiusEqual(p_middle1, tolerance) ||
142
            !arc2.isRadiusEqual(p_middle2, tolerance)) {
143
            can_do_spline_whole = Type::SplitCurve;
144
        }
145
        else {
146
            arc_object1 = arc1.makeArc();
147
            arc_object2 = arc2.makeArc();
148
        }
149
    }
150

151
    if (can_do_spline_whole == Type::SingleArc) {
152
        new_spans.push_back(arc_object1);
153
        new_spans.push_back(arc_object2);
154
    }
155
    else if (can_do_spline_whole == Type::SplitCurve) {
156
        double t_middle = t_start + ((t_end - t_start) * 0.5);
157
        gp_Pnt p_middle;
158
        gp_Vec v_middle;
159
        createArcs(tolerance, new_spans, p_start, v_start, t_start, t_middle, p_middle, v_middle);// recursive
160
        gp_Pnt new_p_end;
161
        gp_Vec new_v_end;
162
        createArcs(tolerance, new_spans, p_middle, v_middle, t_middle, t_end, new_p_end, new_v_end);
163
    }
164
    else {
165
        // calculate_biarc_points failed, just add a line
166
        Part::GeomLineSegment* line = new Part::GeomLineSegment();
167
        line->setPoints(Base::convertTo<Base::Vector3d>(p_start),Base::convertTo<Base::Vector3d>(p_end));
168
        new_spans.push_back(line);
169
    }
170
}
171

172
BSplineCurveBiArcs::Type
173
BSplineCurveBiArcs::calculateBiArcPoints(double t_start, const gp_Pnt& p0, gp_Vec v_start,
174
                                         double t_end, const gp_Pnt& p4, gp_Vec v_end,
175
                                         gp_Pnt& p1, gp_Pnt& p2, gp_Pnt& p3) const
176
{
177
    if (v_start.Magnitude() < Precision::Intersection())
178
        v_start = gp_Vec(p0, p1);
179
    if (v_end.Magnitude() < Precision::Intersection())
180
        v_end = gp_Vec(p3, p4);
181

182
    v_start.Normalize();
183
    v_end.Normalize();
184

185
    gp_Vec v = p0.XYZ() - p4.XYZ();
186

187
    double a = 2*(v_start*v_end-1);
188
    double c = v*v;
189
    double b = (v*2)*(v_start+v_end);
190
    if (fabs(a) < Precision::Intersection()) {
191
        // Check the tangent of a value between t_start and t_end
192
        double t_mid = 0.9 * t_start + 0.1 * t_end;
193
        if (fabs(t_mid) > 0.1) {
194
            gp_Pnt p_mid;
195
            gp_Vec v_mid;
196
            this->myCurve->D1(t_mid, p_mid, v_mid);
197
            v_mid.Normalize();
198
            double a = 2*(v_start*v_mid-1);
199
            if (fabs(a) >= Precision::Intersection()) {
200
                return Type::SplitCurve;
201
            }
202
        }
203
        return Type::SingleLine;
204
    }
205

206

207
    double d = b*b-4*a*c;
208
    if (d < 0.0)
209
        return Type::SingleLine;
210

211
    double sd = sqrt(d);
212
    double e1 = (-b - sd) / (2.0 * a);
213
    double e2 = (-b + sd) / (2.0 * a);
214
    if (e1 > 0 && e2 > 0)
215
        return Type::SingleLine;
216

217
    double e = e1;
218
    if (e2 > e)
219
        e = e2;
220
    if (e < 0)
221
        return Type::SingleLine;
222

223
    p1 = p0.XYZ() + v_start.XYZ() * e;
224
    p3 = p4.XYZ() - v_end.XYZ() * e;
225
    p2 = p1.XYZ() * 0.5 + p3.XYZ() * 0.5;
226

227

228
    return Type::SingleArc;
229
}
230

231
BSplineCurveBiArcs::BSplineCurveBiArcs(const Handle(Geom_Curve)& c)
232
    : myCurve(c)
233
{
234

235
}
236

237
std::list<Geometry*> BSplineCurveBiArcs::toBiArcs(double tolerance) const
238
{
239
    gp_Pnt p_start;
240
    gp_Vec v_start;
241
    gp_Pnt p_end;
242
    gp_Vec v_end;
243
    this->myCurve->D0(this->myCurve->FirstParameter(), p_start);
244
    this->myCurve->D0(this->myCurve->LastParameter(), p_end);
245

246
    std::list<Geometry*> list;
247

248
    // the spline is closed
249
    if (p_start.Distance(p_end) < Precision::Intersection()) {
250
        this->myCurve->D1(this->myCurve->FirstParameter(), p_start, v_start);
251
        createArcs(tolerance, list, p_start, v_start, this->myCurve->FirstParameter(),
252
                   this->myCurve->LastParameter()/2, p_end, v_end);
253
        this->myCurve->D1(this->myCurve->LastParameter()/2, p_start, v_start);
254
        createArcs(tolerance, list, p_start, v_start, this->myCurve->LastParameter()/2,
255
                   this->myCurve->LastParameter(), p_end, v_end);
256
    }
257
    else {
258
        this->myCurve->D1(this->myCurve->FirstParameter(), p_start, v_start);
259
        createArcs(tolerance, list, p_start, v_start, this->myCurve->FirstParameter(),
260
                   this->myCurve->LastParameter(), p_end, v_end);
261
    }
262

263
    return list;
264
}
265