summaryrefslogtreecommitdiffstats
path: root/src/2geom/conicsec.h
blob: b84edc699894edd98d13f30cb2cddb6ec2b0bebe (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
/** @file
 * @brief Conic Section
 *//*
 * Authors:
 *      Nathan Hurst <njh@njhurst.com>
 *
 * Copyright 2009  authors
 *
 * This library is free software; you can redistribute it and/or
 * modify it either under the terms of the GNU Lesser General Public
 * License version 2.1 as published by the Free Software Foundation
 * (the "LGPL") or, at your option, under the terms of the Mozilla
 * Public License Version 1.1 (the "MPL"). If you do not alter this
 * notice, a recipient may use your version of this file under either
 * the MPL or the LGPL.
 *
 * You should have received a copy of the LGPL along with this library
 * in the file COPYING-LGPL-2.1; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
 * You should have received a copy of the MPL along with this library
 * in the file COPYING-MPL-1.1
 *
 * The contents of this file are subject to the Mozilla Public License
 * Version 1.1 (the "License"); you may not use this file except in
 * compliance with the License. You may obtain a copy of the License at
 * http://www.mozilla.org/MPL/
 *
 * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
 * OF ANY KIND, either express or implied. See the LGPL or the MPL for
 * the specific language governing rights and limitations.
 */


#ifndef LIB2GEOM_SEEN_CONICSEC_H
#define LIB2GEOM_SEEN_CONICSEC_H

#include <2geom/exception.h>
#include <2geom/angle.h>
#include <2geom/rect.h>
#include <2geom/affine.h>
#include <2geom/point.h>
#include <2geom/line.h>
#include <2geom/bezier-curve.h>
#include <2geom/numeric/linear_system.h>
#include <2geom/numeric/symmetric-matrix-fs.h>
#include <2geom/numeric/symmetric-matrix-fs-operation.h>

#include <boost/optional/optional.hpp>

#include <string>
#include <vector>
#include <ostream>




namespace Geom
{

class RatQuad{
    /**
     * A curve of the form B02*A + B12*B*w + B22*C/(B02 + B12*w + B22)
     * These curves can exactly represent a piece conic section less than a certain angle (find out)
     *
     **/

public:
    Point P[3];
    double w;
    RatQuad() {}
    RatQuad(Point a, Point b, Point c, double w) : w(w) {
        P[0] = a;
        P[1] = b;
        P[2] = c;
    }
    double lambda() const;

    static RatQuad fromPointsTangents(Point P0, Point dP0,
                                      Point P,
                                      Point P2, Point dP2);
    static RatQuad circularArc(Point P0, Point P1, Point P2);

    CubicBezier toCubic() const;
    CubicBezier toCubic(double lam) const;

    Point pointAt(double t) const;
    Point at0() const {return P[0];}
    Point at1() const {return P[2];}

    void split(RatQuad &a, RatQuad &b) const;

    D2<SBasis> hermite() const;
    std::vector<SBasis> homogeneous() const;
};




class xAx{
public:
    double c[6];

    enum kind_t
    {
        PARABOLA,
        CIRCLE,
        REAL_ELLIPSE,
        IMAGINARY_ELLIPSE,
        RECTANGULAR_HYPERBOLA,
        HYPERBOLA,
        DOUBLE_LINE,
        TWO_REAL_PARALLEL_LINES,
        TWO_IMAGINARY_PARALLEL_LINES,
        TWO_REAL_CROSSING_LINES,
        TWO_IMAGINARY_CROSSING_LINES, // crossing at a real point
        SINGLE_POINT = TWO_IMAGINARY_CROSSING_LINES,
        UNKNOWN
    };


    xAx() {}

    /*
     *  Define the conic section by its algebraic equation coefficients
     *
     *  c0, .., c5: equation coefficients
     */
    xAx (double c0, double c1, double c2, double c3, double c4, double c5)
    {
        set (c0, c1, c2, c3, c4, c5);
    }

    /*
     *  Define a conic section by its related symmetric matrix
     */
    xAx (const NL::ConstSymmetricMatrixView<3> & C)
    {
        set(C);
    }

    /*
     *  Define a conic section by computing the one that fits better with
     *  N points.
     *
     *  points: points to fit
     *
     *  precondition: there must be at least 5 non-overlapping points
     */
    xAx (std::vector<Point> const& points)
    {
        set (points);
    }

    /*
     *  Define a section conic by providing the coordinates of one of its
     *  vertex,the major axis inclination angle and the coordinates of its foci
     *  with respect to the unidimensional system defined by the major axis with
     *  origin set at the provided vertex.
     *
     *  _vertex :   section conic vertex V
     *  _angle :    section conic major axis angle
     *  _dist1:     +/-distance btw V and nearest focus
     *  _dist2:     +/-distance btw V and farest focus
     *
     *  prerequisite: _dist1 <= _dist2
     */
    xAx (const Point& _vertex, double _angle, double _dist1, double _dist2)
    {
        set (_vertex, _angle, _dist1, _dist2);
    }

    /*
     *  Define a conic section by providing one of its vertex and its foci.
     *
     *  _vertex: section conic vertex
     *  _focus1: section conic focus
     *  _focus2: section conic focus
     */
    xAx (const Point& _vertex, const Point& _focus1, const Point& _focus2)
    {
        set(_vertex, _focus1, _focus2);
    }

    /*
     *  Define a conic section by passing a focus, the related directrix,
     *  and the eccentricity (e)
     *  (e < 1 -> ellipse; e = 1 -> parabola; e > 1 -> hyperbola)
     *
     *  _focus:         a focus of the conic section
     *  _directrix:     the directrix related to the given focus
     *  _eccentricity:  the eccentricity parameter of the conic section
     */
    xAx (const Point & _focus, const Line & _directrix, double _eccentricity)
    {
        set (_focus, _directrix, _eccentricity);
    }

    /*
     *  Made up a degenerate conic section as a pair of lines
     *
     *  l1, l2: lines that made up the conic section
     */
    xAx (const Line& l1, const Line& l2)
    {
        set (l1, l2);
    }

    /*
     *  Define the conic section by its algebraic equation coefficients
     *  c0, ..., c5: equation coefficients
     */
    void set (double c0, double c1, double c2, double c3, double c4, double c5)
    {
        c[0] = c0;   c[1] = c1;   c[2] = c2; // xx, xy, yy
        c[3] = c3;   c[4] = c4; // x, y
        c[5] = c5; // 1
    }

    /*
     *  Define a conic section by its related symmetric matrix
     */
    void set (const NL::ConstSymmetricMatrixView<3> & C)
    {
        set(C(0,0), 2*C(1,0), C(1,1), 2*C(2,0), 2*C(2,1), C(2,2));
    }

    void set (std::vector<Point> const& points);

    void set (const Point& _vertex, double _angle, double _dist1, double _dist2);

    void set (const Point& _vertex, const Point& _focus1, const Point& _focus2);

    void set (const Point & _focus, const Line & _directrix, double _eccentricity);

    void set (const Line& l1, const Line& l2);


    static xAx fromPoint(Point p);
    static xAx fromDistPoint(Point p, double d);
    static xAx fromLine(Point n, double d);
    static xAx fromLine(Line l);
    static xAx fromPoints(std::vector<Point> const &pts);


    template<typename T>
    T evaluate_at(T x, T y) const {
        return c[0]*x*x + c[1]*x*y + c[2]*y*y + c[3]*x + c[4]*y + c[5];
    }

    double valueAt(Point P) const;

    std::vector<double> implicit_form_coefficients() const {
        return std::vector<double>(c, c+6);
    }

    template<typename T>
    T evaluate_at(T x, T y, T w) const {
        return c[0]*x*x + c[1]*x*y + c[2]*y*y + c[3]*x*w + c[4]*y*w + c[5]*w*w;
    }

    xAx scale(double sx, double sy) const;

    Point gradient(Point p) const;

    xAx operator-(xAx const &b) const;
    xAx operator+(xAx const &b) const;
    xAx operator+(double const &b) const;
    xAx operator*(double const &b) const;

    std::vector<Point> crossings(Rect r) const;
    boost::optional<RatQuad> toCurve(Rect const & bnd) const;
    std::vector<double> roots(Point d, Point o) const;

    std::vector<double> roots(Line const &l) const;

    static Interval quad_ex(double a, double b, double c, Interval ivl);

    Geom::Affine hessian() const;

    boost::optional<Point> bottom() const;

    Interval extrema(Rect r) const;


    /*
     *  Return the symmetric matrix related to the conic section.
     *  Modifying the matrix does not modify the conic section
     */
    NL::SymmetricMatrix<3> get_matrix() const
    {
        NL::SymmetricMatrix<3> C(c);
        C(1,0) *= 0.5; C(2,0) *= 0.5; C(2,1) *= 0.5;
        return C;
    }

    /*
     *  Return the i-th coefficient of the conic section algebraic equation
     *  Modifying the returned value does not modify the conic section coefficient
     */
    double coeff (size_t i) const
    {
        return c[i];
    }

    /*
     *  Return the i-th coefficient of the conic section algebraic equation
     *  Modifying the returned value modifies the conic section coefficient
     */
    double& coeff (size_t i)
    {
        return c[i];
    }

    kind_t kind () const;

    std::string categorise() const;

    /*
     *  Return true if the equation:
     *  c0*x^2 + c1*xy + c2*y^2 + c3*x + c4*y +c5 == 0
     *  really defines a conic, false otherwise
     */
    bool is_quadratic() const
    {
        return (coeff(0) != 0 || coeff(1) != 0 || coeff(2) != 0);
    }

    /*
     *  Return true if the conic is degenerate, i.e. if the related matrix
     *  determinant is null, false otherwise
     */
    bool isDegenerate() const
    {
        return (det_sgn (get_matrix()) == 0);
    }

    /*
     *  Compute the centre of symmetry of the conic section when it exists,
     *  else it return an uninitialized boost::optional<Point> instance.
     */
    boost::optional<Point> centre() const
    {
        typedef boost::optional<Point> opt_point_t;

        double d = coeff(1) * coeff(1) - 4 * coeff(0) * coeff(2);
        if (are_near (d, 0))  return opt_point_t();
        NL::Matrix Q(2, 2);
        Q(0,0) = coeff(0);
        Q(1,1) = coeff(2);
        Q(0,1) = Q(1,0) = coeff(1) * 0.5;
        NL::Vector T(2);
        T[0] = - coeff(3) * 0.5;
        T[1] = - coeff(4) * 0.5;

        NL::LinearSystem ls (Q, T);
        NL::Vector sol = ls.SV_solve();
        Point C;
        C[0] = sol[0];
        C[1] = sol[1];

        return opt_point_t(C);
    }

    double axis_angle() const;

    void roots (std::vector<double>& sol, Coord v, Dim2 d) const;

    xAx translate (const Point & _offset) const;

    xAx rotate (double angle) const;

    /*
     *  Rotate the conic section by the given angle wrt the provided point.
     *
     *  _rot_centre:   the rotation centre
     *  _angle:        the rotation angle
     */
    xAx rotate (const Point & _rot_centre, double _angle) const
    {
        xAx result
            = translate (-_rot_centre).rotate (_angle).translate (_rot_centre);
        return result;
    }

    /*
     *  Compute the tangent line of the conic section at the provided point
     *
     *  _point: the conic section point the tangent line pass through
     */
    Line tangent (const Point & _point) const
    {
        NL::Vector pp(3);
        pp[0] = _point[0]; pp[1] = _point[1]; pp[2] = 1;
        NL::SymmetricMatrix<3> C = get_matrix();
        NL::Vector line = C * pp;
        return Line(line[0], line[1], line[2]);
    }

    /*
     *  For a non degenerate conic compute the dual conic.
     *  TODO: investigate degenerate case
     */
    xAx dual () const
    {
        //assert (! isDegenerate());
        NL::SymmetricMatrix<3> C = get_matrix();
        NL::SymmetricMatrix<3> D = adj(C);
        xAx dc(D);
        return dc;
    }

    bool decompose (Line& l1, Line& l2) const;

    /*
     *  Generate a RatQuad object from a conic arc.
     *
     *  p0: the initial point of the arc
     *  p1: the inner point of the arc
     *  p2: the final point of the arc
     */
    RatQuad toRatQuad (const Point & p0,
                       const Point & p1,
                       const Point & p2) const
    {
        Point dp0 = gradient (p0);
        Point dp2 = gradient (p2);
        return
            RatQuad::fromPointsTangents (p0, rot90 (dp0), p1, p2, rot90 (dp2));
    }

    /*
     *  Return the angle related to the normal gradient computed at the passed
     *  point.
     *
     *  _point: the point at which computes the angle
     *
     *  prerequisite: the passed point must lie on the conic
     */
    double angle_at (const Point & _point) const
    {
        double angle = atan2 (gradient (_point));
        if (angle < 0)  angle += (2*M_PI);
        return angle;
    }

    /*
     *  Return true if the given point is contained in the conic arc determined
     *  by the passed points.
     *
     *  _point:     the point to be tested
     *  _initial:   the initial point of the arc
     *  _inner:     an inner point of the arc
     *  _final:     the final point of the arc
     *
     *  prerequisite: the passed points must lie on the conic, the inner point
     *                has to be strictly contained in the arc, except when the
     *                initial and final points are equal: in such a case if the
     *                inner point is also equal to them, then they define an arc
     *                made up by a single point.
     *
     */
    bool arc_contains (const Point & _point, const Point & _initial,
                       const Point & _inner, const Point & _final) const
    {
        AngleInterval ai(angle_at(_initial), angle_at(_inner), angle_at(_final));
        return ai.contains(angle_at(_point));
    }

    Rect arc_bound (const Point & P1, const Point & Q, const Point & P2) const;

    std::vector<Point> allNearestTimes (const Point &P) const;

    /*
     *  Return the point on the conic section nearest to the passed point "P".
     *
     *  P: the point to compute the nearest one
     */
    Point nearestTime (const Point &P) const
    {
        std::vector<Point> points = allNearestTimes (P);
        if ( !points.empty() )
        {
            return points.front();
        }
        // else
        THROW_LOGICALERROR ("nearestTime: no nearest point found");
        return Point();
    }

};

std::vector<Point> intersect(const xAx & C1, const xAx & C2);

bool clip (std::vector<RatQuad> & rq, const xAx & cs, const Rect & R);

inline std::ostream &operator<< (std::ostream &out_file, const xAx &x) {
    for(int i = 0; i < 6; i++) {
        out_file << x.c[i] << ", ";
    }
    return out_file;
}

};


#endif // LIB2GEOM_SEEN_CONICSEC_H

/*
  Local Variables:
  mode:c++
  c-file-style:"stroustrup"
  c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
  indent-tabs-mode:nil
  fill-column:99
  End:
*/
// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99 :