bezier-curve.cpp

Go to the documentation of this file.

/* Bezier curve implementation
 *
 * Authors:
 *   MenTaLguY <mental@rydia.net>
 *   Marco Cecchetti <mrcekets at gmail.com>
 *   Krzysztof Kosiński <tweenk.pl@gmail.com>
 *
 * Copyright 2007-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.
 */
 
#include <2geom/bezier-curve.h>
#include <2geom/path-sink.h>
#include <2geom/basic-intersection.h>
#include <2geom/nearest-time.h>
#include <2geom/polynomial.h>
 
namespace Geom
{
 
BezierCurve::BezierCurve(std::vector<Point> const &pts)
    : inner(pts)
{
    if (pts.size() < 2) {
        THROW_RANGEERROR("Bezier curve must have at least 2 control points");
    }
}
 
bool BezierCurve::isDegenerate() const
{
    for (unsigned d = 0; d < 2; ++d) {
        Coord ic = inner[d][0];
        for (unsigned i = 1; i < size(); ++i) {
            if (inner[d][i] != ic) return false;
        }
    }
    return true;
}
 
bool BezierCurve::isLineSegment() const
{
    auto const last_idx = size() - 1;
    if (last_idx == 1) {
        return true;
    }
    auto const start = controlPoint(0);
    auto const end = controlPoint(last_idx);
    for (unsigned i = 1; i < last_idx; ++i) {
        auto const pi = controlPoint(i);
        if (pi != start && pi != end) {
            return false;
        }
    }
    return true;
}
 
void BezierCurve::expandToTransformed(Rect &bbox, Affine const &transform) const
{
    bbox |= bounds_exact(inner * transform);
}
 
Coord BezierCurve::length(Coord tolerance) const
{
    switch (order())
    {
    case 0:
        return 0.0;
    case 1:
        return distance(initialPoint(), finalPoint());
    case 2:
        {
            std::vector<Point> pts = controlPoints();
            return bezier_length(pts[0], pts[1], pts[2], tolerance);
        }
    case 3:
        {
            std::vector<Point> pts = controlPoints();
            return bezier_length(pts[0], pts[1], pts[2], pts[3], tolerance);
        }
    default:
        return bezier_length(controlPoints(), tolerance);
    }
}
 
std::vector<CurveIntersection>
BezierCurve::intersect(Curve const &other, Coord eps) const
{
    std::vector<CurveIntersection> result;
 
    // in case we encounter an order-1 curve created from a vector
    // or a degenerate elliptical arc
    if (isLineSegment()) {
        LineSegment ls(initialPoint(), finalPoint());
        result = ls.intersect(other);
        return result;
    }
 
    // here we are sure that this curve is at least a quadratic Bezier
    BezierCurve const *bez = dynamic_cast<BezierCurve const *>(&other);
    if (bez) {
        std::vector<std::pair<double, double> > xs;
        find_intersections(xs, inner, bez->inner, eps);
        for (auto & i : xs) {
            CurveIntersection x(*this, other, i.first, i.second);
            result.push_back(x);
        }
        return result;
    }
 
    // pass other intersection types to the other curve
    result = other.intersect(*this, eps);
    transpose_in_place(result);
    return result;
}
 
bool BezierCurve::isNear(Curve const &c, Coord precision) const
{
    if (this == &c) return true;
 
    BezierCurve const *other = dynamic_cast<BezierCurve const *>(&c);
    if (!other) return false;
 
    if (!are_near(inner.at0(), other->inner.at0(), precision)) return false;
    if (!are_near(inner.at1(), other->inner.at1(), precision)) return false;
 
    if (size() == other->size()) {
        for (unsigned i = 1; i < order(); ++i) {
            if (!are_near(inner.point(i), other->inner.point(i), precision)) {
                return false;
            }
        }
        return true;
    } else {
        // Must equalize the degrees before comparing
        BezierCurve elevated_this, elevated_other;
        for (size_t dim = 0; dim < 2; dim++) {
            unsigned const our_degree = inner[dim].degree();
            unsigned const other_degree = other->inner[dim].degree();
 
            if (our_degree < other_degree) {
                // Elevate our degree
                elevated_this.inner[dim] = inner[dim].elevate_to_degree(other_degree);
                elevated_other.inner[dim] = other->inner[dim];
            } else if (our_degree > other_degree) {
                // Elevate the other's degree
                elevated_this.inner[dim] = inner[dim];
                elevated_other.inner[dim] = other->inner[dim].elevate_to_degree(our_degree);
            } else {
                // Equal degrees: just copy
                elevated_this.inner[dim] = inner[dim];
                elevated_other.inner[dim] = other->inner[dim];
            }
        }
        assert(elevated_other.size() == elevated_this.size());
        return elevated_this.isNear(elevated_other, precision);
    }
}
 
Curve *BezierCurve::portion(Coord f, Coord t) const
{
    if (f == 0.0 && t == 1.0) {
        return duplicate();
    }
    if (f == 1.0 && t == 0.0) {
        return reverse();
    }
    return new BezierCurve(Geom::portion(inner, f, t));
}
 
bool BezierCurve::_equalTo(Curve const &c) const
{
    if (this == &c) return true;
    auto other = dynamic_cast<BezierCurve const *>(&c);
    if (!other) return false;
 
    if (size() != other->size()) return false;
 
    for (unsigned i = 0; i < size(); ++i) {
        if (controlPoint(i) != other->controlPoint(i)) return false;
    }
 
    return true;
}
 
Coord BezierCurve::nearestTime(Point const &p, Coord from, Coord to) const
{
    return nearest_time(p, inner, from, to);
}
 
void BezierCurve::feed(PathSink &sink, bool moveto_initial) const
{
    if (size() > 4) {
        Curve::feed(sink, moveto_initial);
        return;
    }
 
    Point ip = controlPoint(0);
    if (moveto_initial) {
        sink.moveTo(ip);
    }
    switch (size()) {
    case 2:
        sink.lineTo(controlPoint(1));
        break;
    case 3:
        sink.quadTo(controlPoint(1), controlPoint(2));
        break;
    case 4:
        sink.curveTo(controlPoint(1), controlPoint(2), controlPoint(3));
        break;
    default:
        // TODO: add a path sink method that accepts a vector of control points
        //       and converts to cubic spline by default
        assert(false);
        break;
    }
}
 
BezierCurve *BezierCurve::create(std::vector<Point> const &pts)
{
    switch (pts.size()) {
    case 0:
    case 1:
        THROW_LOGICALERROR("BezierCurve::create: too few points in vector");
        return NULL;
    case 2:
        return new LineSegment(pts[0], pts[1]);
    case 3:
        return new QuadraticBezier(pts[0], pts[1], pts[2]);
    case 4:
        return new CubicBezier(pts[0], pts[1], pts[2], pts[3]);
    default:
        return new BezierCurve(pts);
    }
}
 
std::vector<Coord> BezierCurve::timesWithRadiusOfCurvature(double radius) const
{
    if (this->size() <= 2) {
        return {};
    }
 
    std::vector<Coord> res;
 
    auto const dx = Geom::derivative(inner[X]);
    auto const dy = Geom::derivative(inner[Y]);
    auto const ddx = Geom::derivative(dx);
    auto const ddy = Geom::derivative(dy);
    auto const c0 = (dx * ddy - ddx * dy) * radius;
    auto const c1 = dx * dx + dy * dy;
    auto const p = c0 * c0 - c1 * c1 * c1;
    auto const candidates = p.roots();
    // check which candidates have positive nominator
    // as squaring also give negative (spurious) results
    for (Coord const candidate : candidates) {
        if (c0.valueAt(candidate) > 0) {
            res.push_back(candidate);
        }
    }
 
    return res;
}
 
// optimized specializations for LineSegment
 
template <>
Curve *BezierCurveN<1>::derivative() const {
    double dx = inner[X][1] - inner[X][0], dy = inner[Y][1] - inner[Y][0];
    return new BezierCurveN<1>(Point(dx,dy),Point(dx,dy));
}
 
template<>
Coord BezierCurveN<1>::nearestTime(Point const& p, Coord from, Coord to) const
{
    using std::swap;
 
    if ( from > to ) swap(from, to);
    Point ip = pointAt(from);
    Point fp = pointAt(to);
    Point v = fp - ip;
    Coord l2v = L2sq(v);
    if (l2v == 0) return 0;
    Coord t = dot( p - ip, v ) / l2v;
    if ( t <= 0 )       return from;
    else if ( t >= 1 )  return to;
    else return from + t*(to-from);
}
 
/* Specialized intersection routine for line segments.
 *
 * NOTE: if the segments overlap in part or in full, the function returns the start and end
 * of the overlapping subsegment as intersections. This behavior is more useful than throwing
 * Geom::InfinitelyManySolutions.
 */
template <>
std::vector<CurveIntersection> BezierCurveN<1>::intersect(Curve const &other, Coord eps) const
{
    std::vector<CurveIntersection> result;
 
    // only handle intersections with other LineSegments here
    if (!other.isLineSegment()) {
        // pass all other types to the other curve
        result = other.intersect(*this, eps);
        transpose_in_place(result);
        return result;
    }
 
    Point const u = finalPoint() - initialPoint();
    Point const v = other.initialPoint() - other.finalPoint();
    if (u.isZero() || v.isZero()) {
        return {};
    }
    Coord const uv = u.length() * v.length();
    Coord const u_cross_v = cross(u, v);
    bool const segments_are_parallel = std::abs(u_cross_v) < eps * uv;
 
    if (segments_are_parallel) {
        // We check if the segments lie on the same line.
        Coord const distance_between_lines = std::abs(cross(u.normalized(),
                                                            other.initialPoint() - initialPoint()));
        if (distance_between_lines > eps) {
            // Segments are parallel but aren't part of the same line => no intersections.
            return {};
        }
        // The segments are on the same line, so they may overlap in part or in full.
        // Look for the times on this segment's line at which the line passes through
        // the initial and final points of the other segment.
        Coord const ulen_inverse = 1.0 / u.length();
        auto const time_of_passage = [&](Point const &point_on_line) -> Coord {
            return dot(u.normalized(), point_on_line - initialPoint()) * ulen_inverse;
        };
        // Find the range of times on our segment where we travel through the other segment.
        auto time_in_other = Interval(time_of_passage(other.initialPoint()),
                                      time_of_passage(other.finalPoint()));
        Coord const eps_utime = eps * ulen_inverse;
        if (time_in_other.min() > 1 + eps_utime || time_in_other.max() < -eps_utime) {
            return {};
        }
 
        // Create two intersections, one at each end of the overlap interval.
        Coord last_time = infinity();
        for (Coord t : {time_in_other.min(), time_in_other.max()}) {
            t = std::clamp(t, 0.0, 1.0);
            if (t == last_time) {
                continue;
            }
            last_time = t;
            auto const point = pointAt(t);
            Coord const other_t = std::clamp(dot(v.normalized(), other.initialPoint() - point) / v.length(), 0.0, 1.0);
            auto const other_pt = other.pointAt(other_t);
            if (distance(point, other_pt) > eps) {
                continue;
            }
            result.emplace_back(t, other_t, middle_point(point, other_pt));
        }
        return result;
    } else {
        // Segments are not collinear - there is 0 or 1 intersection.
        // This is the typical case so it's important for it to be fast.
        Point const w = other.initialPoint() - initialPoint();
        Coord candidate_time_this = cross(w, v) / u_cross_v;
 
        auto const is_outside_seg = [eps](Coord t, Coord len) -> bool {
            Coord const arclen_coord = t * len;
            return arclen_coord < -eps || arclen_coord > len + eps;
        };
 
        if (is_outside_seg(candidate_time_this, u.length())) {
            return {};
        }
 
        Coord candidate_time_othr = cross(u, w) / u_cross_v;
        if (is_outside_seg(candidate_time_othr, v.length())) {
            return {};
        }
 
        // Even if there was some fuzz in the interval, we must now restrict to [0, 1] and verify
        // that the intersection found is still within epsilon precision (the fuzz may have been too
        // large, depending on the angle between the intervals).
        candidate_time_this = std::clamp(candidate_time_this, 0.0, 1.0);
        candidate_time_othr = std::clamp(candidate_time_othr, 0.0, 1.0);
 
        auto const pt_this = pointAt(candidate_time_this);
        auto const pt_othr = other.pointAt(candidate_time_othr);
        if (distance(pt_this, pt_othr) > eps) {
            return {};
        }
        return { CurveIntersection(candidate_time_this, candidate_time_othr, middle_point(pt_this, pt_othr)) };
    }
}
 
template <unsigned degree>
static std::vector<CurveIntersection> bezier_line_intersections(BezierCurveN<degree> const &curve, Line const &line)
{
    static_assert(degree == 2 || degree == 3, "bezier_line_intersections<degree>() error: degree must be 2 or 3.");
 
    auto const length = distance(line.initialPoint(), line.finalPoint());
    if (length == 0) {
        return {};
    }
    std::vector<CurveIntersection> result;
 
    // Find the isometry mapping the line to the x-axis, taking the initial point to the origin
    // and the final point to (length, 0). Apply this transformation to the Bézier curve and
    // extract the y-coordinate polynomial.
    auto const transform = line.rotationToZero(Y);
    Bezier const y = (curve.fragment() * transform)[Y];
    std::vector<double> roots;
 
    // Find roots of the polynomial y.
    {
        double const c2 = y[0] + y[2] - 2.0 * y[1];
        double const c1 = y[1] - y[0];
        double const c0 = y[0];
 
        if constexpr (degree == 2) {
            roots = solve_quadratic(c2, 2.0 * c1, c0);
        } else if constexpr (degree == 3) {
            double const c3 = y[3] - y[0] + 3.0 * (y[1] - y[2]);
            roots = solve_cubic(c3, 3.0 * c2, 3.0 * c1 , c0);
        }
    }
 
    // Filter the roots and assemble intersections.
    for (double root : roots) {
        if (root < 0.0 || root > 1.0) {
            continue;
        }
        Coord x = (curve.pointAt(root) * transform)[X];
        if (x < 0.0 || x > length) {
            continue;
        }
        result.emplace_back(curve, line, root, x / length);
    }
    return result;
}
 
/* Specialized intersection routine for quadratic Bézier curves. */
template <>
std::vector<CurveIntersection> BezierCurveN<2>::intersect(Curve const &other, Coord eps) const
{
    if (auto other_bezier = dynamic_cast<BezierCurve const *>(&other)) {
        auto const other_degree = other_bezier->order();
        if (other_degree == 1) {
            // Use the exact method to intersect a quadratic Bézier with a line segment.
            auto line = Line(other_bezier->initialPoint(), other_bezier->finalPoint());
            return bezier_line_intersections<2>(*this, line);
        }
        // TODO: implement exact intersection of two quadratic Béziers using the method of resultants.
    }
    return BezierCurve::intersect(other, eps);
}
 
/* Specialized intersection routine for cubic Bézier curves. */
template <>
std::vector<CurveIntersection> BezierCurveN<3>::intersect(Curve const &other, Coord eps) const
{
    if (auto other_bezier = dynamic_cast<BezierCurve const *>(&other)) {
        if (other_bezier->order() == 1) {
            // Use the exact method to intersect a cubic Bézier with a line segment.
            auto line = Line(other_bezier->initialPoint(), other_bezier->finalPoint());
            return bezier_line_intersections<3>(*this, line);
        }
    }
    return BezierCurve::intersect(other, eps);
}
 
template <>
int BezierCurveN<1>::winding(Point const &p) const
{
    Point ip = inner.at0(), fp = inner.at1();
    if (p[Y] == std::max(ip[Y], fp[Y])) return 0;
 
    Point v = fp - ip;
    assert(v[Y] != 0);
    Coord t = (p[Y] - ip[Y]) / v[Y];
    assert(t >= 0 && t <= 1);
    Coord xcross = lerp(t, ip[X], fp[X]);
    if (xcross > p[X]) {
        return v[Y] > 0 ? 1 : -1;
    }
    return 0;
}
 
template <>
void BezierCurveN<1>::feed(PathSink &sink, bool moveto_initial) const
{
    if (moveto_initial) {
        sink.moveTo(controlPoint(0));
    }
    sink.lineTo(controlPoint(1));
}
 
template <>
void BezierCurveN<2>::feed(PathSink &sink, bool moveto_initial) const
{
    if (moveto_initial) {
        sink.moveTo(controlPoint(0));
    }
    sink.quadTo(controlPoint(1), controlPoint(2));
}
 
template <>
void BezierCurveN<3>::feed(PathSink &sink, bool moveto_initial) const
{
    if (moveto_initial) {
        sink.moveTo(controlPoint(0));
    }
    sink.curveTo(controlPoint(1), controlPoint(2), controlPoint(3));
}
 
static void bezier_expand_to_image(Rect &range, Point const &x0, Point const &x1, Point const &x2)
{
    for (auto i : { X, Y }) {
        bezier_expand_to_image(range[i], x0[i], x1[i], x2[i]);
    }
}
 
static void bezier_expand_to_image(Rect &range, Point const &x0, Point const &x1, Point const &x2, Point const &x3)
{
    for (auto i : { X, Y }) {
        bezier_expand_to_image(range[i], x0[i], x1[i], x2[i], x3[i]);
    }
}
 
template <>
void BezierCurveN<1>::expandToTransformed(Rect &bbox, Affine const &transform) const
{
    bbox.expandTo(finalPoint() * transform);
}
 
template <>
void BezierCurveN<2>::expandToTransformed(Rect &bbox, Affine const &transform) const
{
    bezier_expand_to_image(bbox, controlPoint(0) * transform,
                                 controlPoint(1) * transform,
                                 controlPoint(2) * transform);
}
 
template <>
void BezierCurveN<3>::expandToTransformed(Rect &bbox, Affine const &transform) const
{
    bezier_expand_to_image(bbox, controlPoint(0) * transform,
                                 controlPoint(1) * transform,
                                 controlPoint(2) * transform,
                                 controlPoint(3) * transform);
}
 
static Coord bezier_length_internal(std::vector<Point> &v1, Coord tolerance, int level)
{
    /* The Bezier length algorithm used in 2Geom utilizes a simple fact:
     * the Bezier curve is longer than the distance between its endpoints
     * but shorter than the length of the polyline formed by its control
     * points. When the difference between the two values is smaller than the
     * error tolerance, we can be sure that the true value is no further than
     * 0.5 * tolerance from their arithmetic mean. When it's larger, we recursively
     * subdivide the Bezier curve into two parts and add their lengths.
     *
     * We cap the maximum number of subdivisions at 256, which corresponds to 8 levels.
     */
    Coord lower = distance(v1.front(), v1.back());
    Coord upper = 0.0;
    for (size_t i = 0; i < v1.size() - 1; ++i) {
        upper += distance(v1[i], v1[i+1]);
    }
    if (upper - lower <= 2*tolerance || level >= 8) {
        return (lower + upper) / 2;
    }
 
 
    std::vector<Point> v2 = v1;
 
    /* Compute the right subdivision directly in v1 and the left one in v2.
     * Explanation of the algorithm used:
     * We have to compute the left and right edges of this triangle in which
     * the top row are the control points of the Bezier curve, and each cell
     * is equal to the arithmetic mean of the cells directly above it
     * to the right and left. This corresponds to subdividing the Bezier curve
     * at time value 0.5: the left edge has the control points of the first
     * portion of the Bezier curve and the right edge - the second one.
     * In the example we subdivide a curve with 5 control points (order 4).
     *
     * Start:
     * 0 1 2 3 4
     *  ? ? ? ?
     *   ? ? ?
     *    ? ?
     *     ?
     * # means we have overwritten the value, ? means we don't know
     * the value yet. Numbers mean the value is at i-th position in the vector.
     *
     * After loop with i==1
     * # 1 2 3 4
     *  0 ? ? ? -> write 0 to v2[1]
     *   ? ? ?
     *    ? ?
     *     ?
     *
     * After loop with i==2
     * # # 2 3 4
     *  # 1 ? ?
     *   0 ? ? -> write 0 to v2[2]
     *    ? ?
     *     ?
     *
     * After loop with i==3
     * # # # 3 4
     *  # # 2 ?
     *   # 1 ?
     *    0 ? -> write 0 to v2[3]
     *     ?
     *
     * After loop with i==4, we have the right edge of the triangle in v1,
     * and we write the last value needed for the left edge in v2[4].
     */
 
    for (size_t i = 1; i < v1.size(); ++i) {
        for (size_t j = i; j > 0; --j) {
            v1[j-1] = 0.5 * (v1[j-1] + v1[j]);
        }
        v2[i] = v1[0];
    }
 
    return bezier_length_internal(v1, 0.5 * tolerance, level + 1) +
           bezier_length_internal(v2, 0.5 * tolerance, level + 1);
}
 
Coord bezier_length(std::vector<Point> const &points, Coord tolerance)
{
    if (points.size() < 2) return 0.0;
    std::vector<Point> v1 = points;
    return bezier_length_internal(v1, tolerance, 0);
}
 
static Coord bezier_length_internal(Point a0, Point a1, Point a2, Coord tolerance, int level)
{
    Coord lower = distance(a0, a2);
    Coord upper = distance(a0, a1) + distance(a1, a2);
 
    if (upper - lower <= 2*tolerance || level >= 8) {
        return (lower + upper) / 2;
    }
 
    Point // Casteljau subdivision
        // b0 = a0,
        // c0 = a2,
        b1 = 0.5*(a0 + a1),
        c1 = 0.5*(a1 + a2),
        b2 = 0.5*(b1 + c1); // == c2
    return bezier_length_internal(a0, b1, b2, 0.5 * tolerance, level + 1) +
           bezier_length_internal(b2, c1, a2, 0.5 * tolerance, level + 1);
}
 
Coord bezier_length(Point a0, Point a1, Point a2, Coord tolerance)
{
    return bezier_length_internal(a0, a1, a2, tolerance, 0);
}
 
static Coord bezier_length_internal(Point a0, Point a1, Point a2, Point a3, Coord tolerance, int level)
{
    Coord lower = distance(a0, a3);
    Coord upper = distance(a0, a1) + distance(a1, a2) + distance(a2, a3);
 
    if (upper - lower <= 2*tolerance || level >= 8) {
        return (lower + upper) / 2;
    }
 
    Point // Casteljau subdivision
        // b0 = a0,
        // c0 = a3,
        b1 = 0.5*(a0 + a1),
        t0 = 0.5*(a1 + a2),
        c1 = 0.5*(a2 + a3),
        b2 = 0.5*(b1 + t0),
        c2 = 0.5*(t0 + c1),
        b3 = 0.5*(b2 + c2); // == c3
    return bezier_length_internal(a0, b1, b2, b3, 0.5 * tolerance, level + 1) +
           bezier_length_internal(b3, c2, c1, a3, 0.5 * tolerance, level + 1);
}
 
Coord bezier_length(Point a0, Point a1, Point a2, Point a3, Coord tolerance)
{
    return bezier_length_internal(a0, a1, a2, a3, tolerance, 0);
}
 
} // end namespace Geom
 
/*
  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 :