Adding upstream version 1.3.upstream/1.3upstream

Signed-off-by: Daniel Baumann <daniel.baumann@progress-linux.org>

1 files changed, 537 insertions, 0 deletions

diff --git a/include/2geom/conicsec.h b/include/2geom/conicsec.h
new file mode 100644
index 0000000..bfd5f36
--- /dev/null
+++ b/include/2geom/conicsec.h
@@ -0,0 +1,537 @@
+/** @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 <optional>
+
+#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;
+    std::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;
+
+    std::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 std::optional<Point> instance.
+     */
+    std::optional<Point> centre() const
+    {
+        typedef std::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;
+
+    /**
+     * @brief Division-free decomposition of a degenerate conic section, without
+     *        degeneration test.
+     *
+     * When the conic is degenerate, it consists of 0, 1 or 2 lines in the xy-plane.
+     * This function returns these lines. But it does not check whether the conic
+     * is really degenerate, so calling it on a non-degenerate conic produces a
+     * meaningless result.
+     *
+     * If the number of lines is less than two, the trailing lines in the returned
+     * array will be degenerate. Use Line::isDegenerate() to test for that.
+     *
+     * This version of the decomposition is division-free, which improves numerical
+     * stability compared to decompose().
+     *
+     * @param epsilon The numerical threshold for floating point comparison of discriminants.
+     */
+    std::array<Line, 2> decompose_df(Coord epsilon = EPSILON) 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(double i : x.c) {
+        out_file << 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 :
+