summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorAndrea Canciani <ranma42@gmail.com>2010-06-28 16:47:52 +0200
committerAndrea Canciani <ranma42@gmail.com>2010-10-23 10:52:47 +0200
commit0e02921b8b0f674d8570fd4b929909b25d48fa02 (patch)
treedf72cc1221d6f2f9e2b542746f89eed1035a0cd9
parente3dd9f3856201c6b38d747971a7d5dfd18a7f635 (diff)
-rw-r--r--src/cairo-spline-offset-jm.c488
-rw-r--r--src/cairo-spline-offset.c992
-rw-r--r--src/cairo-spline-offset.h20
-rw-r--r--src/cairo-spline-private.h292
-rw-r--r--src/cairo-stroke-to-path.c186
-rw-r--r--src/cairo-vector-private.h134
6 files changed, 1015 insertions, 1097 deletions
diff --git a/src/cairo-spline-offset-jm.c b/src/cairo-spline-offset-jm.c
new file mode 100644
index 000000000..01283d95e
--- /dev/null
+++ b/src/cairo-spline-offset-jm.c
@@ -0,0 +1,488 @@
+/* -*- Mode: c; tab-width: 8; c-basic-offset: 4; indent-tabs-mode: t; -*- */
+/* cairo - a vector graphics library with display and print output
+ *
+ * Copyright © 2009 Mozilla Foundation
+ * Copyright © 2010 Jeff Muizelaar
+ *
+ * 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.
+ *
+ * The Original Code is the cairo graphics library.
+ *
+ * The Initial Developer of the Original Code is Jeff Muizelaar
+ *
+ */
+
+#define _GNU_SOURCE
+
+#include "cairoint.h"
+
+#include "cairo-spline-offset.h"
+
+/* Finds the coefficient with the maximum error. This should be
+ * the coefficient closest to area of maximum error.
+ * This returns an upper bound of the error. */
+static double
+error_upper_bound (const double bezier[7], int degree, double stop_value, double *t)
+{
+ double worst_t = .5;
+ double max_error = 0;
+ int i;
+
+ for (i=0; i<7; i++)
+ if (fabs(bezier[i]) > error) {
+ error = fabs(bezier[i]);
+ worst_t = i/6.;
+ }
+
+ *t = worst_t;
+
+ /* we shouldn't ever split at the edges because the approximation
+ * should be perfect there. XXX: we might be able to take
+ * advantage of this when computing the error polynomial. */
+#if 0
+ if (*t == 0 || *t == 1) {
+ for (i=0; i<=degree; i++) {
+ error = fabs(bezier[i]);
+ printf("error: %f\n", error);
+ }
+ assert(0);
+ }
+#endif
+ /* *t = 0.5; */
+ return error >= stop_value;
+}
+
+static cairo_bool_t
+error_within_bounds_elber (knots_t bz1, knots_t bz_offset_approx,
+ double desired_offset,
+ double max_error,
+ double *t)
+{
+ double bzprod[3+3 +1] = {};
+
+ /* Elber and Cohnen propose the following method of computing the
+ error between an curve and it's offset:
+ e(t) = || C(t) - Coffset_approx(t) ||^2 - offset^2
+
+ It is important to note that this function doesn't represent a unique offset curve.
+ For example, it's possible for the offset curve to cross the original spline
+ provided it stays the appropriate distance away:
+
+ xxx
+ x
+ ****
+ **x
+ * x
+ * x
+ * x
+ **
+ ****
+ x
+ xxx
+
+ "Planar Curve Offset Base on Circle Approximation" suggests further problems
+ with this error measure. They give the following example: C(t) + (r, 0) will
+ have an error mesaure of 0, but is not even close to the actual offset curve.
+
+ This paper suggests using the angle between the difference vector and the normal
+ at point t, but doesn't go into much detail.
+ */
+
+ /* Now we compute this error function.
+ * First, subtract the approx */
+ bz1.a.x -= bz_offset_approx.a.x;
+ bz1.a.y -= bz_offset_approx.a.y;
+ bz1.b.x -= bz_offset_approx.b.x;
+ bz1.b.y -= bz_offset_approx.b.y;
+ bz1.c.x -= bz_offset_approx.c.x;
+ bz1.c.y -= bz_offset_approx.c.y;
+ bz1.d.x -= bz_offset_approx.d.x;
+ bz1.d.y -= bz_offset_approx.d.y;
+/*
+ Next we compute the dot product of bz1 with itself This involves squaring
+ the x and y polynomials and adding the result together.
+
+ However, we want to multiply the coefficients so that we can treat the
+ result as a higher degree bezier curve.
+
+ The following method is given in the Symbolic Computation section of
+ "Comparing Offset Curve Approximation Methods" by Elber, Lee and Kim
+
+ k l ⎛ k⎞ i k-i⎛ l⎞ j l-j
+ B (t) B (t) = ⎜ ⎟ t (1-t) ⎜ ⎟ t (1-t)
+ i j ⎝ i⎠ ⎝ j⎠
+ ⎛ k⎞ ⎛ l⎞ i+j k+l-i-j
+ = ⎜ ⎟ ⎜ ⎟ t (1-t)
+ ⎝ i⎠ ⎝ j⎠
+
+ ⎛ k⎞ ⎛ l⎞
+ = ⎝ i⎠ ⎝ j⎠ k+l
+ ──────── B (t)
+ ⎛ k + l⎞ i+j
+ ⎝ i + j⎠
+
+ m n
+C1(t) C2(t) = ∑ P_i B_i,m (t) ∑ P_j B_j,n (t)
+ i=0 j=0
+
+ m n
+ = ∑ ∑ P_i P_j B_i,m (t) B_j,n (t)
+ i=0 j=0
+
+ m+n m + n
+ = ∑ Q_k B_k (t)
+ k=0
+
+Where Q_k is:
+ ⎛ m⎞ ⎛ n⎞
+ ⎝ i⎠ ⎝ j⎠
+ Q_k = P_i * P_j ─────────
+ ⎛ m+n⎞
+ ⎝ i+j⎠
+
+And for the case of degree 3 beziers
+
+ ⎛ 3⎞ ⎛ 3⎞
+ ⎝ i⎠ ⎝ j⎠
+ Q_k = P_i * P_j ─────────
+ ⎛ 6 ⎞
+ ⎝ i+j⎠
+
+ This expands to:
+
+ i j
+ 0 0 0 = (1 * 1) / 1 = 1
+
+ 0 1 1 = (1 * 3) / 6 = .5
+ 1 0 1 = (3 * 1) / 6 = .5
+
+ 0 2 2 = (1 * 3) / 15 = .2
+ 1 1 2 = (3 * 3) / 15 = .6
+ 2 0 2 = (3 * 1) / 15 = .2
+
+ 0 3 3 = (1 * 1) / 20 = .05
+ 1 2 3 = (3 * 3) / 20 = .45
+ 2 1 3 = (3 * 3) / 20 = .45
+ 3 0 3 = (1 * 1) / 20 = .05
+
+ 1 3 4 = (3 * 1) / 15 = .2
+ 2 2 4 = (3 * 3) / 15 = .6
+ 3 1 4 = (1 * 3) / 15 = .2
+
+ 2 3 5 = (3 * 1) / 6 = .5
+ 3 2 5 = (1 * 3) / 6 = .5
+
+ 3 3 6 = (1 * 1) / 1 = 1
+
+Which simplifies to the following:
+ */
+
+
+ bzprod[0] += bz1.a.x * bz1.a.x;
+ bzprod[0] += bz1.a.y * bz1.a.y;
+
+ bzprod[1] += bz1.a.x * bz1.b.x;
+ bzprod[1] += bz1.a.y * bz1.b.y;
+
+ bzprod[2] += bz1.a.x * bz1.c.x * 0.4;
+ bzprod[2] += bz1.b.x * bz1.b.x * 0.6;
+ bzprod[2] += bz1.a.y * bz1.c.y * 0.4;
+ bzprod[2] += bz1.b.y * bz1.b.y * 0.6;
+
+ bzprod[3] += bz1.d.x * bz1.a.x * 0.1;
+ bzprod[3] += bz1.b.x * bz1.c.x * 0.9;
+ bzprod[3] += bz1.d.y * bz1.a.y * 0.1;
+ bzprod[3] += bz1.b.y * bz1.c.y * 0.9;
+
+ bzprod[4] += bz1.b.x * bz1.d.x * 0.4;
+ bzprod[4] += bz1.c.x * bz1.c.x * 0.6;
+ bzprod[4] += bz1.d.y * bz1.b.y * 0.4;
+ bzprod[4] += bz1.c.y * bz1.c.y * 0.6;
+
+ bzprod[5] += bz1.c.x * bz1.d.x;
+ bzprod[5] += bz1.d.y * bz1.c.y;
+
+ bzprod[6] += bz1.d.x * bz1.d.x;
+ bzprod[6] += bz1.d.y * bz1.d.y;
+
+
+ /* the result is a bezier polynomial that represents the distance squared between
+ * the two input polynomials */
+
+ /* we don't need to subtract the desired offset because our metric doesn't depend on being centered around 0 */
+#if 1
+ bzprod[0] -= (desired_offset*desired_offset);
+ bzprod[1] -= (desired_offset*desired_offset);
+ bzprod[2] -= (desired_offset*desired_offset);
+ bzprod[3] -= (desired_offset*desired_offset);
+ bzprod[4] -= (desired_offset*desired_offset);
+ bzprod[5] -= (desired_offset*desired_offset);
+ bzprod[6] -= (desired_offset*desired_offset);
+#endif
+ *t = 0.5;
+ return !is_curvy_t(bzprod, sizeof(bzprod)/sizeof(bzprod[0])-1, max_error, t);
+}
+
+void
+print_knot (knots_t c) {
+ /*printf("(%.13a %.13a) (%.13a %.13a) (%.13a %.13a) (%.13a %.13a)\n", c.a.x, c.a.y, c.b.x, c.b.y, c.c.x, c.c.y, c.d.x, c.d.y);*/
+ printf("(%.5f %.5f) (%.5f %.5f) (%.5f %.5f) (%.5f %.5f)\n", c.a.x, c.a.y, c.b.x, c.b.y, c.c.x, c.c.y, c.d.x, c.d.y);
+}
+
+static knots_t
+convolve_with_circle (knots_t spline, double dist)
+{
+ vector_t in_normal = normalize(perp(begin_tangent(spline)));
+ vector_t out_normal = normalize(perp(end_tangent(spline)));
+ /* find an arc the goes from the input normal to the output_normal */
+ knots_t arc_spline = arc_segment_cart(in_normal, out_normal, 1, in_normal, out_normal);
+
+ /* approximate convolving 'spline' with the arc spline (XXX: is convolve the right word?) */
+ spline.a.x += dist * arc_spline.a.x;
+ spline.a.y += dist * arc_spline.a.y;
+ spline.b.x += dist * arc_spline.b.x;
+ spline.b.y += dist * arc_spline.b.y;
+ spline.c.x += dist * arc_spline.c.x;
+ spline.c.y += dist * arc_spline.c.y;
+ spline.d.x += dist * arc_spline.d.x;
+ spline.d.y += dist * arc_spline.d.y;
+
+ return spline;
+}
+
+/* This approach is inspired by "Approximation of circular arcs and offset curves by Bezier curves of high degree"
+ * by YJ Ahn, Y soo Kim, Y Shin.
+ *
+ * This is a similar idea to:
+ * "Offset approximation based on reparameterizaing the path of a moving point along the base circle", ZHAO, WANG
+ * and a bunch of papers by Lee et al.
+ * "Planar curve offset based on circle approximation"
+ * "New approximation methods of planar offset and convolution curves"
+ * "Polynomial/rational approximation of Minkowski sum boundary curves"
+ *
+ * However these other approaches produce curves of higher degree than the input curve making them unsuitable
+ * for use here.
+ * */
+static cairo_status_t
+approximate_offset_curve_with_shin (const cairo_spline_knots_t *self, double dist, cairo_status_t (*curve_fn)(void *, const cairo_spline_knots_t *), void *closure)
+{
+ /* since we're going to approximating the offset curve with arcs
+ we want to ensure that we don't have any inflection points in our spline */
+ double t_inflect[2] = _cairo_spline_inflection_points (k);
+
+ /* split the curve at the inflection points and convolve each with an approximation of
+ * a circle */
+ cairo_spline_knots_t remaining = *self;
+ double adjusted_inflect = t_inflect[1];
+
+ /* check that the first inflection point is in range and not equal to t2
+ * (we avoid reparametrizing the curve and just do the second split if
+ * t1==t2) */
+ if (t_inflect[0] > 0 && t_inflect[0] < 1 && t_inflect[0] < t_inflect[1]) {
+ /* split at the inflection point */
+ cairo_spline_knots_t first;
+ _cairo_spline_de_casteljau (&remaining, &first, &remaining, t_inflect[0]);
+
+ /* approximate */
+ if (!_cairo_spline_isdegenerate (&first)) {
+ convolve_with_circle (&first, &first, dist);
+ status = curve_fn(closure, &first);
+ if (unlikely (status))
+ return status;
+ }
+
+ /* reparamaterize the second inflection point over the 'remaining' spline:
+ * (domain between inflection points)/(domain of remaining) */
+ t_inflect[1] = (t_inflect[1] - t_inflect[0]) / (1 - t_inflect[0]);
+ }
+
+ if (t_inflect[1] > 0 && t_inflect[1] < 1) {
+ /* split at the inflection point */
+ cairo_spline_knots_t second;
+ _cairo_spline_de_casteljau (&remaining, &second, &remaining, t_inflect[1]);
+
+ /* approximate */
+ if (!_cairo_spline_isdegenerate (&second)) {
+ convolve_with_circle (&second, &second, dist);
+ status = curve_fn(closure, &second);
+ if (unlikely (status))
+ return status;
+ }
+ }
+
+ /* approximate the last part of the spline */
+ if (!_cairo_spline_isdegenerate (&remaining)) {
+ convolve_with_circle (&remaining, &remaining, dist);
+ status = curve_fn(closure, &remaining);
+ }
+ return status;
+}
+
+/* Computes the offset polygon for the control polygon of spline. We can
+ * use this is an approximation of the offset spline. This approximation is
+ * give by Tiller W, Hanson E G. in "Offsets of two-dimensional profile" */
+static cairo_bool_t
+knot_offset (cairo_spline_knots_t *dest, const cairo_spline_knots_t *src, double width)
+{
+ /* this code is inspired by the libart mitreing code */
+
+ /* difference vectors between each point in the knot */
+ cairo_vector_t diff[3];
+ /* offset and average offset vectors */
+ cairo_vector_t off[3], mid[2];
+
+ diff[0] = _cairo_vector_sub (self->b, self->a);
+ diff[1] = _cairo_vector_sub (self->c, self->b);
+ diff[2] = _cairo_vector_sub (self->d, self->c);
+
+ /* compute normalized normals */
+ /* deal with any degeneracies in the spline by treating
+ * degenerate segments as having the same normal
+ * as their neighbour. If all of the segments are degenerate
+ * than we will fail the is_parallel test below */
+ for (i=0; i<3; i++)
+ if (!_cairo_vector_iszero (diff[i]))
+ off[i] = _cairo_vector_perpendicular (_cairo_vector_normalize (diff[i]));
+
+ if (_cairo_vector_iszero (diff[1]))
+ off[1] = _cairo_vector_iszero(diff[2]) ? off[0] : off[2]
+
+ if (_cairo_vector_iszero (diff[2]))
+ off[2] = off[1];
+
+ if (_cairo_vector_iszero (diff[0]))
+ off[0] = off[1];
+
+ /* mid-point vector sum */
+ mid[0] = _cairo_vector_add (off[0], off[1]);
+ mid[1] = _cairo_vector_add (off[1], off[2]);
+
+ /* TODO: check nongeneracy */
+ /* the choice of this EPSILON is arbitrary and could use further study */
+#define PARALLEL_EPSILON 1e-6
+ if (fabs (_cairo_vector_dot(mid[0], mid[0])) < PARALLEL_EPSILON || fabs (_cairo_vector_dot(mid[1], mid[1])) < PARALLEL_EPSILON)
+ return FALSE;
+
+ mid[0] = _cairo_vector_normalize (mid[0]);
+ mid[1] = _cairo_vector_normalize (mid[1]);
+
+ /* write out the result */
+ dest->a = _cairo_vector_add (src->a, _cairo_vector_scale (off[0], width));
+ dest->b = _cairo_vector_add (src->b, _cairo_vector_scale (mid[0], 2 * width));
+ dest->c = _cairo_vector_add (src->c, _cairo_vector_scale (mid[1], 2 * width));
+ dest->d = _cairo_vector_add (src->d, _cairo_vector_scale (off[2], width));
+
+ return TRUE;
+}
+
+/* plan:
+ * if we approximate the original bezier curve with line segments
+ * and then produce an offset to those lines. We will end up
+ * with sharp angles at places of high curvature.
+ *
+ * It is at these places of high curvature that we want to approximate
+ * the offset with an arc instead of continuing to subdivide.
+ *
+ * How do we find the areas of high curvature?
+ * - Ideas:
+ * As we subdivide the original bezier
+ * we can get an idea of how flat it is.
+ * If the original bezier is flat but the offset approximation is not within our bounds
+ * we should stop recursing and approximate the segment with an arc.
+ *
+ * The degree of required flatness of the original curve should be related to the
+ * the offset length. The greater the offset distance the less flat the original
+ * bezier needs to be.
+ *
+ * pros:
+ * - this gives us an easier condition to reason about the termination of the
+ * algorithm because we terminate the recursion when the input bezier has
+ * been subdivided to 'flat' instead of when the offset is good.
+ * This is valuable because it's possible to create an input curve that does not have a
+ * good offset approximation down to the resolution of a double.
+ */
+cairo_status_t
+curve_offset (const cairo_spline_knots_t *k, double dist, double tolerance, cairo_status_t (*curve_fn)(void *, const cairo_spline_knots_t *), void *closure)
+{
+ cairo_knots_t offset;
+ cairo_bool_t recurse = knot_offset(&offset, k, dist);
+ double break_point = 0.5;
+ /* TODO: we could probably choose a better place to split than halfway
+ * for times when we know we're going to recurse */
+
+ if (!recurse) {
+ assert (_cairo_spline_isfinite (offset));
+ /* error is too large. subdivide on max t and try again. */
+ recurse = error_upper_bound (k, &offset, dist, &break_point) > tolerance;
+ }
+
+ if (!recurse) {
+ return curve_fn (closure, &offset);
+ } else if (_cairo_spline_curvyness (k) < tolerance * 0.125) {
+ /* We need to do something to handle regions of high curvature
+ *
+ * skia uses the first and second derivitives at the points of 'max curvature' to check for
+ * pinchynes.
+ * qt tests whether the points are close and reverse direction.
+ * float l = (orig->x1 - orig->x2)*(orig->x1 - orig->x2) + (orig->y3 - orig->y4)*(orig->y3 - orig->y4);
+ * float dot = (orig->x1 - orig->x2)*(orig->x3 - orig->x4) + (orig->y1 - orig->y2)*(orig->y3 - orig->y4);
+ *
+ * we just check if the knot is flat and has large error.
+ */
+ /* we don't want to recurse anymore because we're pretty flat,
+ * instead approximate the offset curve with an arc.
+ *
+ * The previously generated offset curve can become really bad if the the intersections
+ * are far away from the original curve (self) */
+ return approximate_offset_curve_with_shin (self, dist, curve_fn, closure);
+ } else {
+ cairo_bool_t is_continuous = _cairo_spline_de_casteljau (k, &left, &right, break_point);
+
+ status = curve_offset(&left, dist, tolerance, curve_fn, closure);
+ if (unlikely (status))
+ return status;
+
+ if (!is_continuous) {
+ /* add circle:
+ we can use the exit the normal from the left side
+ as the begining normal of our cusp */
+ cairo_vector_t mid = _cairo_vector_normalize (_cairo_spline_end_tangent (&left));
+ cairo_vector_t start = _cairo_vector_perpendicular (mid);
+ cairo_vector_t end = _cairo_vector_opposite (start);
+ cairo_spline_knots_t arc;
+
+ _cairo_spline_arc (&arc, left.d, dist, &start, &mid);
+ status = curve_fn (&arc, closure);
+ if (unlikely (status))
+ return status;
+
+ _cairo_spline_arc (&arc, left.d, dist, &mid, &end);
+ status = curve_fn (&arc, closure);
+ if (unlikely (status))
+ return status;
+ }
+
+ return curve_offset (&right, dist, tolerance, curve_fn, closure);
+ }
+}
diff --git a/src/cairo-spline-offset.c b/src/cairo-spline-offset.c
index b15495ac3..747ed42e2 100644
--- a/src/cairo-spline-offset.c
+++ b/src/cairo-spline-offset.c
@@ -33,941 +33,135 @@
*
*/
-#define _GNU_SOURCE
-
-#include "cairoint.h"
-
-#include <math.h>
-#include <stdio.h>
-#include <assert.h>
-#include <stdlib.h>
-
-#if _XOPEN_SOURCE >= 600 || defined (_ISOC99_SOURCE)
-#define ISFINITE(x) isfinite (x)
-#else
-#define ISFINITE(x) ((x) * (x) >= 0.) /* check for NaNs */
-#endif
-
-#include "cairoint.h"
+#include "cairo-spline-private.h"
#include "cairo-spline-offset.h"
+/* assumes that the input spline is not degenerate and that its normal sweeps less than pi.
+ src and dest MUST be different splines */
static void
-_lerp (const point_t *a, const point_t *b, point_t *result, double t)
-{
- result->x = (1-t)*a->x + t*b->x;
- result->y = (1-t)*a->y + t*b->y;
-}
-
-/* returns false if there could be a discontinuity
- * between the split curves */
-static cairo_bool_t
-_de_casteljau_t (knots_t *k1, knots_t *k2, double t)
-{
- point_t ab, bc, cd;
- point_t abbc, bccd;
- point_t final;
-
- _lerp (&k1->a, &k1->b, &ab, t);
- _lerp (&k1->b, &k1->c, &bc, t);
- _lerp (&k1->c, &k1->d, &cd, t);
- _lerp (&ab, &bc, &abbc, t);
- _lerp (&bc, &cd, &bccd, t);
- _lerp (&abbc, &bccd, &final, t);
-
- k2->a = final;
- k2->b = bccd;
- k2->c = cd;
- k2->d = k1->d;
-
- k1->b = ab;
- k1->c = abbc;
- k1->d = final;
-
- /* we can test whether the split curve will have a cusp
- * inbetween the two parts by checking the length of the
- * line abbc-bccd. If this line has a non-zero length
- * then the two curves will share a common normal
- * at the points where they meet. Therefore the
- * length must be zero for the normals to be different.
- * However this condition is necessary but not sufficient.
- * XXX: does this matter and how should we deal with it if it does? */
- return abbc.x != final.x || abbc.y != final.y;
-}
-
-typedef point_t vector_t;
-
-static vector_t
-begin_tangent (knots_t c)
-{
- vector_t t;
- if (c.a.x != c.b.x || c.a.y != c.b.y) {
- t.x = c.b.x - c.a.x;
- t.y = c.b.y - c.a.y;
- return t;
- }
- if (c.a.x != c.c.x || c.a.y != c.c.y) {
- t.x = c.c.x - c.a.x;
- t.y = c.c.y - c.a.y;
- return t;
- }
- t.x = c.d.x - c.a.x;
- t.y = c.d.y - c.a.y;
- return t;
-}
-
-static vector_t
-end_tangent (knots_t c)
-{
- vector_t t;
- if (c.d.x != c.c.x || c.d.y != c.c.y) {
- t.x = c.d.x - c.c.x;
- t.y = c.d.y - c.c.y;
- return t;
- }
- if (c.d.x != c.b.x || c.d.y != c.b.y) {
- t.x = c.d.x - c.b.x;
- t.y = c.d.y - c.b.y;
- return t;
- }
- t.x = c.d.x - c.a.x;
- t.y = c.d.y - c.a.y;
- return t;
-}
-
-static void
-ensure_finite (knots_t k)
-{
-
- if (!(
- isfinite(k.a.x) &&
- isfinite(k.a.y) &&
- isfinite(k.b.x) &&
- isfinite(k.b.y) &&
- isfinite(k.c.x) &&
- isfinite(k.c.y) &&
- isfinite(k.d.x) &&
- isfinite(k.d.y))) {
- print_knot(k);
- assert(0);
- }
-}
-
-static int
-degenerate_knots (knots_t k)
-{
-
- return k.a.x == k.b.x &&
- k.a.y == k.b.y &&
- k.a.x == k.c.x &&
- k.a.y == k.c.y &&
- k.a.x == k.d.x &&
- k.a.y == k.d.y;
-}
-
-#if 0
-static void ensure_smoothness (knots_t a, knots_t b)
-{
- double dist_len = hypot(a.d.x - b.a.x, a.d.y - b.a.y);
- print_knot(a);
- print_knot(b);
- printf("dist: %f\n", dist_len);
- vector_t in_norm = end_tangent(a);
- vector_t out_norm = begin_tangent(b);
- double in_a = atan2 (in_norm.y, in_norm.x);
- double out_a = atan2 (out_norm.y, out_norm.x);
- if (in_a < out_a)
- in_a += 2*M_PI;
- double diff = in_a - out_a;
- if (diff > M_PI)
- diff -= M_PI*2;
- if (!(fabs(diff) < 0.001) && !(fabs(diff - M_PI) < 0.001) && !(fabs(diff + M_PI) < 0.001)) {
- /* The angle difference between normals can be either 0 or +-M_PI for a change in direction.
- * The change in direction can occur in a region of high curvature with a large offset direction.
- * We move in reverse if the offset distance exceeds the radius of curvature. Note: that you need to
- * use signed curvature for this too work out properly. Because we only encounter this phenomenom in
- * the convex areas of the curve.
- *
- * XXX: explain why... */
- printf("angle diff: %f %f %f\n", diff, in_a, out_a);
- print_knot(a);
- print_knot(b);
- assert(0);
- }
- if (fabs(dist_len) > 0.01) {
- printf("distlen: %f\n", dist_len);
- print_knot(a);
- print_knot(b);
- assert(0);
- }
-}
-
-static void
-ensure_smoothness2 (knots_t a, knots_t b, cairo_bool_t ccw) {
- vector_t in_norm = end_tangent(a);
- vector_t out_norm = begin_tangent(b);
- if (!ccw) {
- in_norm = end_tangent(a);
- out_norm = begin_tangent(b);
- }
- double in_a = atan2 (in_norm.y, in_norm.x);
- double out_a = atan2 (out_norm.y, out_norm.x);
- if (in_a < out_a)
- in_a += 2*M_PI;
- double diff = in_a - out_a;
- if (diff > M_PI)
- diff -= M_PI*2;
- if (fabs(diff) > 0.01) {
- printf("gangle diff: %f\n", diff);
- print_knot(a);
- print_knot(b);
- }
-}
-#endif
-
-static inline vector_t
-flip (vector_t a)
-{
- vector_t ar;
- ar.x = -a.x;
- ar.y = -a.y;
- return ar;
-}
-
-static vector_t
-perp (vector_t v)
-{
- vector_t p = {-v.y, v.x};
- return p;
-}
-
-static inline double
-dot (vector_t a, vector_t b)
+_cairo_spline_offset_circle (cairo_spline_knots_double_t *dest, const cairo_spline_knots_double_t *src, double dist)
{
- return a.x * b.x + a.y * b.y;
-}
+ cairo_vector_t in_normal = _cairo_vector_perpendicular (_cairo_spline_begin_tangent (src));
+ cairo_vector_t out_normal = _cairo_vector_perpendicular (_cairo_spline_begin_tangent (src));
-static vector_t
-normalize (vector_t v)
-{
- double len = hypot(v.x, v.y);
- v.x /= len;
- v.y /= len;
- return v;
-}
+ assert (src != dest);
-/* given a normal rotate the vector 90 degrees to the right clockwise
- * This function has a period of 4. e.g. swap(swap(swap(swap(x) == x */
-static vector_t
-swap (vector_t a)
-{
- vector_t ar;
- /* one of these needs to be negative. We choose a.x so that we rotate to the right instead of negating */
- ar.x = a.y;
- ar.y = -a.x;
- return ar;
-}
+ /* find an arc the goes from the input normal to the output_normal */
+ _cairo_spline_arc (dest, &_cairo_vector_zero, dist, &in_normal, &out_normal);
-static vector_t
-unperp (vector_t a)
-{
- return swap (a);
+ /* approximate convolving 'spline' with the arc spline */
+ _cairo_spline_scale (dest, dest, dist);
+ _cairo_spline_add (dest, dest, src);
}
-typedef struct {
- double t1;
- double t2;
-} inflection_points_t;
-
-/* Find the inflection points for a knots_t.
- * The lower inflection point is returned in t1, the higher one in t2.
- *
- * This method for computing inflections points is from:
- * "Fast, precise flattening of cubic Bezier path and offset curves", Hain et. al */
-static inflection_points_t
-inflection_points (knots_t s)
+static cairo_status_t
+_cairo_spline_do_arc (const cairo_vector_t *center, double radius, const cairo_vector_t *start, const cairo_vector_t *end, double tolerance, cairo_status_t (*curve_fn)(void *, const cairo_spline_knots_double_t *), void *closure)
{
- inflection_points_t ret;
-
- point_t a = {x: -s.a.x + 3*s.b.x - 3*s.c.x + s.d.x, -s.a.y + 3*s.b.y - 3*s.c.y + s.d.y};
- point_t b = {x: 3*s.a.x - 6*s.b.x + 3*s.c.x, 3*s.a.y - 6*s.b.y + 3*s.c.y};
- point_t c = {x:-3*s.a.x + 3*s.b.x, -3*s.a.y + 3*s.b.y};
- /* we don't need 'd'
- point_t d = {x: s.a.x, s.a.y};
- */
-
- /* we're solving for the roots of this equation:
- dx/dt * d2y/dt2 - d2x/dt2 * dy/dt
- == 6*(a.y*b.x - a.x*b.y)*t^2 +6*(a.y*c.x-a.x*c.y)*t + 2*(b.y*c.x - b.x*c.y)
- */
-
- if ((a.y*b.x - a.x*b.y) == 0.0) {
- ret.t1 = ret.t2 = 1.0;
- return ret;
- }
+ cairo_spline_knots_double_t arc;
+ cairo_vector_t tmp, tmp2;
+ double maxsqradius;
+ int subdivision;
+ cairo_status_t status;
- double t_cusp = (-1./2)*((a.y*c.x - a.x*c.y)/(a.y*b.x - a.x*b.y));
- double t_inner = (t_cusp*t_cusp - ((b.y*c.x - b.x*c.y)/(a.y*b.x - a.x*b.y))/3);
+ maxsqradius = radius + tolerance;
+ maxsqradius *= maxsqradius;
- if (t_inner > 0) {
- double t_adj = sqrt (t_inner);
- ret.t1 = t_cusp - t_adj;
- ret.t2 = t_cusp + t_adj;
- } else {
- ret.t1 = ret.t2 = t_cusp;
- }
- return ret;
-}
+ for (subdivision = 1; _cairo_spline_arc (&arc, center, radius, start, &tmp) > maxsqradius; subdivision <<= 1)
+ tmp = _cairo_vector_bisect (*start, tmp);
-/* this function could use some tests.
- I believe it only works for angles from 0 to pi. It will always use the smaller arc between the two vectors? */
-static knots_t
-arc_segment_cart (point_t start, point_t end, double radius, vector_t a, vector_t b)
-{
- knots_t arc;
- double r_sin_A = radius * a.y;
- double r_cos_A = radius * a.x;
- double r_sin_B = radius * b.y;
- double r_cos_B = radius * b.x;
+ switch (subdivision) {
+ default:
+ for (tmp2 = *start; subdivision > 2; subdivision -= 2) {
+ status = curve_fn (closure, &arc);
+ if (unlikely (status))
+ return status;
- point_t mid, mid2;
- double l, h;
+ tmp2 = _cairo_vector_add (tmp, _cairo_vector_opposite (_cairo_vector_sub (tmp2, tmp)));
+ _cairo_spline_arc (&arc, center, radius, &tmp, &tmp2);
+ status = curve_fn (closure, &arc);
+ if (unlikely (status))
+ return status;
- mid.x = a.x + b.x;
- mid.y = a.y + b.y;
- if (dot (a, b) < 0) {
- /* otherwise, we can flip a, add it
- * and then use the perpendicular of the result */
- a = flip (a);
- mid.x = a.x + b.x;
- mid.y = a.y + b.y;
- if (dot (a, perp(b)) <= 0) {
- mid = unperp (mid);
- } else {
- mid = perp(mid);
+ tmp = _cairo_vector_add (tmp2, _cairo_vector_opposite (_cairo_vector_sub (tmp, tmp2)));
+ _cairo_spline_arc (&arc, center, radius, &tmp2, &tmp);
}
+ /* fall through */
+ case 2:
+ status = curve_fn (closure, &arc);
+ if (unlikely (status))
+ return status;
+ _cairo_spline_arc (&arc, center, radius, &tmp, end);
+ /* fall through */
+ case 1:
+ return curve_fn (closure, &arc);
}
- /* normalize */
- l = sqrt(mid.x*mid.x + mid.y*mid.y);
- mid.x /= l;
- mid.y /= l;
-
- mid2.x = a.x + mid.x;
- mid2.y = a.y + mid.y;
-
- h = (4./3.)*dot(perp(a), mid2)/dot(a, mid2);
-
- arc.a.x = start.x;
- arc.a.y = start.y;
-
- arc.b.x = start.x - h * r_sin_A;
- arc.b.y = start.y + h * r_cos_A;
-
- arc.c.x = end.x + h * r_sin_B;
- arc.c.y = end.y - h * r_cos_B;
-
- arc.d.x = end.x;
- arc.d.y = end.y;
-
- return arc;
-}
-
-static knots_t
-quarter_circle (double xc, double yc, double radius, vector_t normal)
-{
- double r_sin_A, r_cos_A;
- double r_sin_B, r_cos_B;
- double h;
- knots_t arc;
-
- r_sin_A = radius * normal.y;
- r_cos_A = radius * normal.x;
- r_sin_B = radius * -normal.x; /* sin (angle_B); */
- r_cos_B = radius * normal.y; /* cos (angle_B); */
-
- h = -0.55228474983079345; /* 4.0/3.0 * tan ((-M_PI/2) / 4.0); */
-
- arc.a.x = xc + r_cos_A;
- arc.a.y = yc + r_sin_A;
-
- arc.b.x = xc + r_cos_A - h * r_sin_A;
- arc.b.y = yc + r_sin_A + h * r_cos_A;
-
- arc.c.x = xc + r_cos_B + h * r_sin_B;
- arc.c.y = yc + r_sin_B - h * r_cos_B;
-
- arc.d.x = xc + r_cos_B;
- arc.d.y = yc + r_sin_B;
-
- return arc;
-}
-
-static void semi_circle (double xc, double yc, double radius, vector_t out_normal, void (*curve_fn)(void *, knots_t), void *closure)
-{
- double len = hypot (out_normal.y, out_normal.x);
- vector_t mid_normal;
-
- out_normal.x /= len;
- out_normal.y /= len;
-
- mid_normal.x = out_normal.y; /* we need to rotate by -M_PI/2 */
- mid_normal.y = -out_normal.x;
-
- curve_fn(closure, quarter_circle (xc, yc, radius, out_normal));
- curve_fn(closure, quarter_circle (xc, yc, radius, mid_normal));
}
static cairo_bool_t
-is_knot_flat (knots_t c)
+_cairo_spline_offset_needs_recursion (const cairo_spline_knots_double_t *k, double dist, double tolerance, double *t)
{
- /*
- A well-known flatness test is both cheaper and more reliable
- than the ones you have tried. The essential observation is that
- when the curve is a uniform speed straight line from end to end,
- the control points are evenly spaced from beginning to end.
- Therefore, our measure of how far we deviate from that ideal
- uses distance of the middle controls, not from the line itself,
- but from their ideal *arrangement*. Point 2 should be halfway
- between points 1 and 3; point 3 should be halfway between points
- 2 and 4.
- This, too, can be improved. Yes, we can eliminate the square roots
- in the distance tests, retaining the Euclidean metric; but the
- taxicab metric is faster, and also safe. The length of displacement
- (x,y) in this metric is |x|+|y|.
- */
-
- /* XXX: this isn't necessarily the best test because it
- tests for flatness instead of smallness.
- however, flat things will be offseted easily so we sort of implicitly do the right thing
- */
- /* this value should be much less than the tolerance because we only want to deal with situations where
- * the curvature is extreme. */
- double stop_value = 0.05;
- double error = 0;
- error += fabs((c.c.x - c.b.x) - (c.b.x - c.a.x));
- error += fabs((c.c.y - c.b.y) - (c.b.y - c.a.y));
- error += fabs((c.c.x - c.b.x) - (c.d.x - c.c.x));
- error += fabs((c.c.y - c.b.y) - (c.d.y - c.c.y));
- return error <= stop_value;
-}
+ double sweep[2];
+
+ _cairo_spline_normal_sweep (k, sweep);
-/* Finds the coefficient with the maximum error. This should be
- * the coefficient closest to area of maximum error.
- * This returns an upper bound of the error. */
-static cairo_bool_t
-is_curvy_t (double *bezier, int degree, double stop_value, double *t)
-{
- double highest_error = .5;
- double error = 0;
- int i=0;
- assert(degree == 6);
- for (i=0; i<=degree; i++) {
- if (fabs(bezier[i]) > error) {
- error = fabs(bezier[i]);
- highest_error = i/(double)degree;
- }
- }
- *t = highest_error;
+ *t = 0.5;
- /* we shouldn't ever split at the edges because the approximation
- * should be perfect there. XXX: we might be able to take
- * advantage of this when computing the error polynomial. */
-#if 0
- if (*t == 0 || *t == 1) {
- for (i=0; i<=degree; i++) {
- error = fabs(bezier[i]);
- printf("error: %f\n", error);
- }
- assert(0);
- }
-#endif
- /* *t = 0.5; */
- return error >= stop_value;
+ return fabs (sweep[0]) > 1.0 || fabs (sweep[1]) > 1.0; /* TODO: additional checks for tolerance */
}
-static cairo_bool_t
-error_within_bounds_elber (knots_t bz1, knots_t bz_offset_approx,
- double desired_offset,
- double max_error,
- double *t)
+/* k is modified by the function */
+static cairo_status_t
+_cairo_spline_offset_simple (cairo_spline_knots_double_t *k, double dist, double tolerance, cairo_status_t (*curve_fn)(void *, const cairo_spline_knots_double_t *), void *closure)
{
- double bzprod[3+3 +1] = {};
-
- /* Elber and Cohnen propose the following method of computing the
- error between an curve and it's offset:
- e(t) = || C(t) - Coffset_approx(t) ||^2 - offset^2
-
- It is important to note that this function doesn't represent a unique offset curve.
- For example, it's possible for the offset curve to cross the original spline
- provided it stays the appropriate distance away:
-
- xxx
- x
- ****
- **x
- * x
- * x
- * x
- **
- ****
- x
- xxx
+ cairo_spline_knots_double_t tmp;
+ cairo_status_t status;
+ double t;
- "Planar Curve Offset Base on Circle Approximation" suggests further problems
- with this error measure. They give the following example: C(t) + (r, 0) will
- have an error mesaure of 0, but is not even close to the actual offset curve.
+ while (1) {
+ if (_cairo_spline_isdegenerate (k))
+ return CAIRO_STATUS_SUCCESS;
- This paper suggests using the angle between the difference vector and the normal
- at point t, but doesn't go into much detail.
- */
-
- /* Now we compute this error function.
- * First, subtract the approx */
- bz1.a.x -= bz_offset_approx.a.x;
- bz1.a.y -= bz_offset_approx.a.y;
- bz1.b.x -= bz_offset_approx.b.x;
- bz1.b.y -= bz_offset_approx.b.y;
- bz1.c.x -= bz_offset_approx.c.x;
- bz1.c.y -= bz_offset_approx.c.y;
- bz1.d.x -= bz_offset_approx.d.x;
- bz1.d.y -= bz_offset_approx.d.y;
-/*
- Next we compute the dot product of bz1 with itself This involves squaring
- the x and y polynomials and adding the result together.
-
- However, we want to multiply the coefficients so that we can treat the
- result as a higher degree bezier curve.
-
- The following method is given in the Symbolic Computation section of
- "Comparing Offset Curve Approximation Methods" by Elber, Lee and Kim
-
- k l ⎛ k⎞ i k-i⎛ l⎞ j l-j
- B (t) B (t) = ⎜ ⎟ t (1-t) ⎜ ⎟ t (1-t)
- i j ⎝ i⎠ ⎝ j⎠
- ⎛ k⎞ ⎛ l⎞ i+j k+l-i-j
- = ⎜ ⎟ ⎜ ⎟ t (1-t)
- ⎝ i⎠ ⎝ j⎠
-
- ⎛ k⎞ ⎛ l⎞
- = ⎝ i⎠ ⎝ j⎠ k+l
- ──────── B (t)
- ⎛ k + l⎞ i+j
- ⎝ i + j⎠
-
- m n
-C1(t) C2(t) = ∑ P_i B_i,m (t) ∑ P_j B_j,n (t)
- i=0 j=0
-
- m n
- = ∑ ∑ P_i P_j B_i,m (t) B_j,n (t)
- i=0 j=0
-
- m+n m + n
- = ∑ Q_k B_k (t)
- k=0
-
-Where Q_k is:
- ⎛ m⎞ ⎛ n⎞
- ⎝ i⎠ ⎝ j⎠
- Q_k = P_i * P_j ─────────
- ⎛ m+n⎞
- ⎝ i+j⎠
-
-And for the case of degree 3 beziers
-
- ⎛ 3⎞ ⎛ 3⎞
- ⎝ i⎠ ⎝ j⎠
- Q_k = P_i * P_j ─────────
- ⎛ 6 ⎞
- ⎝ i+j⎠
-
- This expands to:
-
- i j
- 0 0 0 = (1 * 1) / 1 = 1
-
- 0 1 1 = (1 * 3) / 6 = .5
- 1 0 1 = (3 * 1) / 6 = .5
-
- 0 2 2 = (1 * 3) / 15 = .2
- 1 1 2 = (3 * 3) / 15 = .6
- 2 0 2 = (3 * 1) / 15 = .2
-
- 0 3 3 = (1 * 1) / 20 = .05
- 1 2 3 = (3 * 3) / 20 = .45
- 2 1 3 = (3 * 3) / 20 = .45
- 3 0 3 = (1 * 1) / 20 = .05
-
- 1 3 4 = (3 * 1) / 15 = .2
- 2 2 4 = (3 * 3) / 15 = .6
- 3 1 4 = (1 * 3) / 15 = .2
-
- 2 3 5 = (3 * 1) / 6 = .5
- 3 2 5 = (1 * 3) / 6 = .5
-
- 3 3 6 = (1 * 1) / 1 = 1
-
-Which simplifies to the following:
- */
-
-
- bzprod[0] += bz1.a.x * bz1.a.x;
- bzprod[0] += bz1.a.y * bz1.a.y;
-
- bzprod[1] += bz1.a.x * bz1.b.x;
- bzprod[1] += bz1.a.y * bz1.b.y;
-
- bzprod[2] += bz1.a.x * bz1.c.x * 0.4;
- bzprod[2] += bz1.b.x * bz1.b.x * 0.6;
- bzprod[2] += bz1.a.y * bz1.c.y * 0.4;
- bzprod[2] += bz1.b.y * bz1.b.y * 0.6;
-
- bzprod[3] += bz1.d.x * bz1.a.x * 0.1;
- bzprod[3] += bz1.b.x * bz1.c.x * 0.9;
- bzprod[3] += bz1.d.y * bz1.a.y * 0.1;
- bzprod[3] += bz1.b.y * bz1.c.y * 0.9;
-
- bzprod[4] += bz1.b.x * bz1.d.x * 0.4;
- bzprod[4] += bz1.c.x * bz1.c.x * 0.6;
- bzprod[4] += bz1.d.y * bz1.b.y * 0.4;
- bzprod[4] += bz1.c.y * bz1.c.y * 0.6;
-
- bzprod[5] += bz1.c.x * bz1.d.x;
- bzprod[5] += bz1.d.y * bz1.c.y;
-
- bzprod[6] += bz1.d.x * bz1.d.x;
- bzprod[6] += bz1.d.y * bz1.d.y;
-
-
- /* the result is a bezier polynomial that represents the distance squared between
- * the two input polynomials */
-
- /* we don't need to subtract the desired offset because our metric doesn't depend on being centered around 0 */
-#if 1
- bzprod[0] -= (desired_offset*desired_offset);
- bzprod[1] -= (desired_offset*desired_offset);
- bzprod[2] -= (desired_offset*desired_offset);
- bzprod[3] -= (desired_offset*desired_offset);
- bzprod[4] -= (desired_offset*desired_offset);
- bzprod[5] -= (desired_offset*desired_offset);
- bzprod[6] -= (desired_offset*desired_offset);
-#endif
- *t = 0.5;
- return !is_curvy_t(bzprod, sizeof(bzprod)/sizeof(bzprod[0])-1, max_error, t);
-}
-
-void
-print_knot (knots_t c) {
- /*printf("(%.13a %.13a) (%.13a %.13a) (%.13a %.13a) (%.13a %.13a)\n", c.a.x, c.a.y, c.b.x, c.b.y, c.c.x, c.c.y, c.d.x, c.d.y);*/
- printf("(%.5f %.5f) (%.5f %.5f) (%.5f %.5f) (%.5f %.5f)\n", c.a.x, c.a.y, c.b.x, c.b.y, c.c.x, c.c.y, c.d.x, c.d.y);
-}
-
-static knots_t
-convolve_with_circle (knots_t spline, double dist)
-{
- vector_t in_normal = normalize(perp(begin_tangent(spline)));
- vector_t out_normal = normalize(perp(end_tangent(spline)));
- /* find an arc the goes from the input normal to the output_normal */
- knots_t arc_spline = arc_segment_cart(in_normal, out_normal, 1, in_normal, out_normal);
-
- /* approximate convolving 'spline' with the arc spline (XXX: is convolve the right word?) */
- spline.a.x += dist * arc_spline.a.x;
- spline.a.y += dist * arc_spline.a.y;
- spline.b.x += dist * arc_spline.b.x;
- spline.b.y += dist * arc_spline.b.y;
- spline.c.x += dist * arc_spline.c.x;
- spline.c.y += dist * arc_spline.c.y;
- spline.d.x += dist * arc_spline.d.x;
- spline.d.y += dist * arc_spline.d.y;
-
- return spline;
-}
-
-/* This approach is inspired by "Approximation of circular arcs and offset curves by Bezier curves of high degree"
- * by YJ Ahn, Y soo Kim, Y Shin.
- *
- * This is a similar idea to:
- * "Offset approximation based on reparameterizaing the path of a moving point along the base circle", ZHAO, WANG
- * and a bunch of papers by Lee et al.
- * "Planar curve offset based on circle approximation"
- * "New approximation methods of planar offset and convolution curves"
- * "Polynomial/rational approximation of Minkowski sum boundary curves"
- *
- * However these other approaches produce curves of higher degree than the input curve making them unsuitable
- * for use here.
- * */
-static void
-approximate_offset_curve_with_shin (knots_t self, double dist, void (*curve_fn)(void *, knots_t), void *closure)
-{
- /* since we're going to approximating the offset curve with arcs
- we want to ensure that we don't have any inflection points in our spline */
- inflection_points_t t_inflect = inflection_points(self);
-
- /* split the curve at the inflection points and convolve each with an approximation of
- * a circle */
- knots_t remaining = self;
- double adjusted_inflect = t_inflect.t2;
-
- /* check that the first inflection point is in range */
- if (t_inflect.t1 > 0 && t_inflect.t1 < 1) {
- /* split at the inflection point */
- _de_casteljau_t (&self, &remaining, t_inflect.t1);
-
- /* approximate */
- if (!degenerate_knots (self))
- curve_fn(closure, convolve_with_circle(self, dist));
-
- /* reparamaterize the second inflection point over the 'remaining' spline:
- * (domain between inflection points)/(domain of remaining) */
- adjusted_inflect = (t_inflect.t2-t_inflect.t1)/(1 - t_inflect.t1);
- }
-
- /* check that the second inflection point is in range and not equal to t1.
- * we don't use the reparameterized value so that test remains simple */
- if (t_inflect.t2 > t_inflect.t1 && t_inflect.t2 > 0 && t_inflect.t2 < 1) {
- /* split into one segment */
- knots_t self = remaining;
- _de_casteljau_t (&self, &remaining, adjusted_inflect);
- if (!degenerate_knots (self))
- curve_fn(closure, convolve_with_circle(self, dist));
- }
-
- /* deal with what's left of the spline */
- if (!degenerate_knots (remaining))
- curve_fn(closure, convolve_with_circle(remaining, dist));
-}
-
-/* Computes the offset polygon for the control polygon of spline. We can
- * use this is an approximation of the offset spline. This approximation is
- * give by Tiller W, Hanson E G. in "Offsets of two-dimensional profile" */
-static knots_t
-knot_offset (knots_t self, double width, cairo_bool_t *is_parallel)
-{
- /* this code is inspired by the libart mitreing code */
-
- /* difference vectors between each point in the knot */
- double dx0, dy0, dx1, dy1, dx2, dy2;
- /* difference vector lengths */
- double ld0, ld1, ld2;
- /* temporary vector for dealing with degeneracies */
- double last_x, last_y;
-
- double scale;
-
- /* normalized normal vectors for each difference vector */
- double dlx0, dly0, dlx1, dly1, dlx2, dly2;
-
- /* bisected/midpoint vectors */
- double dm1x, dm1y, dm2x, dm2y;
-
- double dmr2_1, dmr2_2;
-
- knots_t result;
-
- dx0 = self.b.x - self.a.x;
- dy0 = self.b.y - self.a.y;
- ld0 = sqrt(dx0*dx0 + dy0*dy0);
-
- dx1 = self.c.x - self.b.x;
- dy1 = self.c.y - self.b.y;
- ld1 = sqrt(dx1*dx1 + dy1*dy1);
-
- dx2 = self.d.x - self.c.x;
- dy2 = self.d.y - self.c.y;
- ld2 = sqrt(dx2*dx2 + dy2*dy2);
-
- /* compute normalized normals */
+ if (_cairo_spline_curvyness (k) < tolerance) {
+ /* TODO: respect tolerance */
+ _cairo_spline_offset_circle (&tmp, k, dist);
+ return curve_fn (closure, &tmp);
+ }
- /* deal with any degeneracies in the spline by treating
- * degenerate segments as having the same normal
- * as their neighbour. If all of the segments are degenerate
- * than we will fail the is_parallel test below */
- last_x = 0;
- last_y = 0;
+ if (_cairo_spline_offset_needs_recursion (k, dist, tolerance, &t)) {
+ _cairo_spline_de_casteljau (k, &tmp, k, t);
- /* first pass */
- if (ld0) {
- scale = 1. / ld0;
- dlx0 = -dy0 * scale;
- dly0 = dx0 * scale;
- last_x = dlx0;
- last_y = dly0;
- }
- if (ld1) {
- scale = 1. / ld1;
- dlx1 = -dy1 * scale;
- dly1 = dx1 * scale;
- last_x = dlx1;
- last_y = dly1;
- }
- if (ld2) {
- scale = 1. / ld2;
- dlx2 = -dy2 * scale;
- dly2 = dx2 * scale;
- last_x = dlx2;
- last_y = dly2;
- }
+ status = _cairo_spline_offset_simple (&tmp, dist, tolerance, curve_fn, closure);
+ if (unlikely (status))
+ return status;
- /* second pass */
- if (!ld2) {
- dlx2 = last_x;
- dly2 = last_y;
- } else {
- scale = 1. / ld2;
- dlx2 = -dy2 * scale;
- dly2 = dx2 * scale;
- last_x = dlx2;
- last_y = dly2;
- }
- if (!ld1) {
- dlx1 = last_x;
- dly1 = last_y;
- } else {
- scale = 1. / ld1;
- dlx1 = -dy1 * scale;
- dly1 = dx1 * scale;
- last_x = dlx1;
- last_y = dly1;
- }
- if (!ld0) {
- dlx0 = last_x;
- dly0 = last_y;
- } else {
- scale = 1. / ld0;
- dlx0 = -dy0 * scale;
- dly0 = dx0 * scale;
- last_x = dlx0;
- last_y = dly0;
+ /* tail-recursion*/
+ } else {
+ _cairo_spline_offset_circle (&tmp, k, dist);
+ return curve_fn (closure, &tmp);
+ }
}
-
- /* mid-point vector sum */
- dm1x = (dlx0 + dlx1);
- dm1y = (dly0 + dly1);
-
- dm2x = (dlx1 + dlx2);
- dm2y = (dly1 + dly2);
-
- /* length squared of the mid-point vector sum */
- dmr2_1 = (dm1x * dm1x + dm1y * dm1y);
- dmr2_2 = (dm2x * dm2x + dm2y * dm2y);
-
- /* the choice of this EPSILON is arbitrary and could use further study */
-#define PARALLEL_EPSILON 1e-6
- if (fabs(dmr2_1) < PARALLEL_EPSILON || fabs(dmr2_2) < PARALLEL_EPSILON)
- *is_parallel = TRUE;
- else
- *is_parallel = FALSE;
-
- scale = width * 2 / dmr2_1;
- dm1x *= scale;
- dm1y *= scale;
-
- scale = width * 2 / dmr2_2;
- dm2x *= scale;
- dm2y *= scale;
-
- /* write out the result */
- result.a.x = self.a.x + dlx0 * width;
- result.a.y = self.a.y + dly0 * width;
-
- result.b.x = self.b.x + dm1x;
- result.b.y = self.b.y + dm1y;
-
- result.c.x = self.c.x + dm2x;
- result.c.y = self.c.y + dm2y;
-
- result.d.x = self.d.x + dlx2 * width;
- result.d.y = self.d.y + dly2 * width;
-
- return result;
}
-/* plan:
- * if we approximate the original bezier curve with line segments
- * and then produce an offset to those lines. We will end up
- * with sharp angles at places of high curvature.
- *
- * It is at these places of high curvature that we want to approximate
- * the offset with an arc instead of continuing to subdivide.
- *
- * How do we find the areas of high curvature?
- * - Ideas:
- * As we subdivide the original bezier
- * we can get an idea of how flat it is.
- * If the original bezier is flat but the offset approximation is not within our bounds
- * we should stop recursing and approximate the segment with an arc.
- *
- * The degree of required flatness of the original curve should be related to the
- * the offset length. The greater the offset distance the less flat the original
- * bezier needs to be.
- *
- * pros:
- * - this gives us an easier condition to reason about the termination of the
- * algorithm because we terminate the recursion when the input bezier has
- * been subdivided to 'flat' instead of when the offset is good.
- * This is valuable because it's possible to create an input curve that does not have a
- * good offset approximation down to the resolution of a double.
- */
-void
-curve_offset (knots_t self, double dist, double tolerance, void (*curve_fn)(void *, knots_t), void *closure)
+cairo_status_t
+_cairo_spline_offset (const cairo_spline_knots_double_t *k, double dist, double tolerance, cairo_status_t (*curve_fn)(void *, const cairo_spline_knots_double_t *), void *closure)
{
- cairo_bool_t recurse;
- double break_point;
- double max_error = tolerance;
- knots_t offset = knot_offset(self, dist, &recurse);
+ cairo_spline_knots_double_t split[3];
+ cairo_bool_t cusp[2];
+ int i, parts = _cairo_spline_split_inflection (split, cusp, k);
+ cairo_status_t status;
- if (!recurse) {
- /* we need to make sure we have a finite offset before checking the error */
- ensure_finite(offset);
- recurse = !error_within_bounds_elber(self, offset, dist, max_error, &break_point);
- } else {
- /* TODO: we could probably choose a better place to split than halfway
- * for times when we know we're going to recurse */
- break_point = .5;
+ for (i = 0; i < parts - 1; i++) {
+ status = _cairo_spline_offset_simple (split+i, dist, tolerance, curve_fn, closure);
+ if (unlikely (status))
+ return status;
+ if (cusp[i]) {
+ cairo_vector_t begin = _cairo_vector_perpendicular (_cairo_spline_end_tangent (split+i));
+ cairo_vector_t end = _cairo_vector_perpendicular (_cairo_spline_begin_tangent (split+i+1));
+ status = _cairo_spline_do_arc (&split[i].d, dist, &begin, &end, tolerance, curve_fn, closure);
+ if (unlikely (status))
+ return status;
+ }
}
- if (recurse) {
- /* error is too large. subdivide on max t and try again. */
- knots_t left = self;
- knots_t right;
-
- cairo_bool_t is_continuous;
-
- /* We need to do something to handle regions of high curvature
- *
- * skia uses the first and second derivitives at the points of 'max curvature' to check for
- * pinchynes.
- * qt tests whether the points are close and reverse direction.
- * float l = (orig->x1 - orig->x2)*(orig->x1 - orig->x2) + (orig->y3 - orig->y4)*(orig->y3 - orig->y4);
- * float dot = (orig->x1 - orig->x2)*(orig->x3 - orig->x4) + (orig->y1 - orig->y2)*(orig->y3 - orig->y4);
- *
- * we just check if the knot is flat and has large error.
- */
- if (is_knot_flat(self)) {
- /* we don't want to recurse anymore because we're pretty flat,
- * instead approximate the offset curve with an arc.
- *
- * The previously generated offset curve can become really bad if the the intersections
- * are far away from the original curve (self) */
-
- approximate_offset_curve_with_shin (self, dist, curve_fn, closure);
- return;
- }
-
- is_continuous = _de_casteljau_t(&left, &right, break_point);
-
- curve_offset(left, dist, tolerance, curve_fn, closure);
-
- if (!is_continuous) {
- /* add circle:
- we can use the exit the normal from the left side
- as the begining normal of our cusp */
- vector_t normal = perp (end_tangent (left));
- semi_circle (left.d.x, left.d.y, dist, normal, curve_fn, closure);
- }
-
- curve_offset(right, dist, tolerance, curve_fn, closure);
- } else {
- curve_fn(closure, offset);
- }
+ return _cairo_spline_offset_simple (split+parts-1, dist, tolerance, curve_fn, closure);;
}
-
diff --git a/src/cairo-spline-offset.h b/src/cairo-spline-offset.h
index 9faa2cb49..37f712942 100644
--- a/src/cairo-spline-offset.h
+++ b/src/cairo-spline-offset.h
@@ -1,18 +1,8 @@
-typedef cairo_point_double_t point_t;
-typedef struct _knots {
- point_t a,b,c,d;
-} knots_t;
+#ifndef CAIRO_SPLINE_OFFSET_H
+#define CAIRO_SPLINE_OFFSET_H
-void print_knot(knots_t k);
+cairo_status_t
+_cairo_spline_offset (const cairo_spline_knots_double_t *k, double offset, double tolerance, cairo_status_t (*curve_fn)(void *, const cairo_spline_knots_double_t *), void *closure);
-struct curve_list {
- knots_t curve;
- double dist;
- struct curve_list *next;
- struct curve_list *prev;
-};
+#endif /*CAIRO_DIRECTFB_H*/
-typedef struct curve_list curve_list_t;
-
-void
-curve_offset(knots_t self, double offset, double tolerance, void (*curve_to_fn)(void *path, knots_t k), void *closure);
diff --git a/src/cairo-spline-private.h b/src/cairo-spline-private.h
new file mode 100644
index 000000000..d368c22ae
--- /dev/null
+++ b/src/cairo-spline-private.h
@@ -0,0 +1,292 @@
+/* -*- Mode: c; tab-width: 8; c-basic-offset: 4; indent-tabs-mode: t; -*- */
+#ifndef CAIRO_SPLINE_PRIVATE_H
+#define CAIRO_SPLINE_PRIVATE_H
+
+#include "cairo-vector-private.h"
+
+typedef struct _cairo_spline_knots_double {
+ cairo_point_double_t a, b, c, d;
+} cairo_spline_knots_double_t;
+
+static cairo_bool_t
+_cairo_spline_isdegenerate (const cairo_spline_knots_double_t *k)
+{
+ return
+ _cairo_vector_equals (k->a, k->b) &&
+ _cairo_vector_equals (k->a, k->c) &&
+ _cairo_vector_equals (k->a, k->d);
+}
+
+static cairo_vector_t
+_cairo_spline_begin_tangent (const cairo_spline_knots_double_t *c)
+{
+ cairo_vector_t r;
+ assert (!_cairo_spline_isdegenerate (c));
+
+ if (!_cairo_vector_equals (c->a, c->b))
+ r = _cairo_vector_sub (c->b, c->a);
+ else if (!_cairo_vector_equals (c->a, c->c))
+ r = _cairo_vector_sub (c->c, c->a);
+ else
+ r = _cairo_vector_sub (c->d, c->a);
+
+ return _cairo_vector_normalize (r);
+}
+
+static cairo_vector_t
+_cairo_spline_end_tangent (const cairo_spline_knots_double_t *c)
+{
+ cairo_vector_t r;
+ assert (!_cairo_spline_isdegenerate (c));
+
+ if (!_cairo_vector_equals (c->d, c->c))
+ r = _cairo_vector_sub (c->d, c->c);
+ else if (!_cairo_vector_equals (c->d, c->b))
+ r = _cairo_vector_sub (c->d, c->b);
+ else
+ r = _cairo_vector_sub (c->d, c->a);
+ return _cairo_vector_normalize (r);
+}
+
+static cairo_bool_t
+_cairo_spline_isfinite (const cairo_spline_knots_double_t *k)
+{
+ return
+ _cairo_vector_isfinite (k->a) &&
+ _cairo_vector_isfinite (k->b) &&
+ _cairo_vector_isfinite (k->c) &&
+ _cairo_vector_isfinite (k->d);
+}
+
+static double
+_cairo_spline_curvyness (const cairo_spline_knots_double_t *k)
+{
+ /*
+ A well-known flatness test is both cheaper and more reliable
+ than the ones you have tried. The essential observation is that
+ when the curve is a uniform speed straight line from end to end,
+ the control points are evenly spaced from beginning to end.
+ Therefore, our measure of how far we deviate from that ideal
+ uses distance of the middle controls, not from the line itself,
+ but from their ideal *arrangement*. Point 2 should be halfway
+ between points 1 and 3; point 3 should be halfway between points
+ 2 and 4.
+ This, too, can be improved. Yes, we can eliminate the square roots
+ in the distance tests, retaining the Euclidean metric; but the
+ taxicab metric is faster, and also safe. The length of displacement
+ (x,y) in this metric is |x|+|y|.
+ */
+ cairo_vector_t ba = _cairo_vector_sub (k->b, k->a);
+ cairo_vector_t cb = _cairo_vector_sub (k->c, k->b);
+ cairo_vector_t dc = _cairo_vector_sub (k->d, k->c);
+ cairo_vector_t e1 = _cairo_vector_sub (cb, ba);
+ cairo_vector_t e2 = _cairo_vector_sub (cb, dc);
+ return fabs (e1.x) + fabs (e1.y) + fabs (e2.x) + fabs (e2.y);
+}
+
+
+static cairo_bool_t
+_cairo_spline_de_casteljau (const cairo_spline_knots_double_t *k, cairo_spline_knots_double_t *k1, cairo_spline_knots_double_t *k2, double t)
+{
+ cairo_vector_t tmp = _cairo_vector_lerp (k->b, k->c, t); /* bc */
+
+ k1->a = k->a;
+ k2->d = k->d;
+
+ k1->b = _cairo_vector_lerp (k->a, k->b, t); /* ab */
+ k2->c = _cairo_vector_lerp (k->c, k->d, t); /* cd */
+
+ k1->c = _cairo_vector_lerp (k1->b, tmp, t); /* abbc */
+ k2->b = _cairo_vector_lerp (tmp, k2->c, t); /* bccd */
+
+ k1->d = _cairo_vector_lerp (k1->c, k2->b, t);
+ k2->a = k1->d;
+
+ /* we can test whether the split curve will have a cusp
+ * inbetween the two parts by checking abbc and bccd.
+ * A cusp can be present iff they are equal (unless the
+ * spline is degenerate) */
+ return _cairo_vector_equals (k1->c, k2->b);
+}
+
+/* Find the inflection points for a knots_t.
+ * The lower inflection point is returned in t1, the higher one in t2.
+ *
+ * This method for computing inflections points is from:
+ * "Fast, precise flattening of cubic Bezier path and offset curves", Hain et. al */
+static int
+_cairo_spline_inflection_points (const cairo_spline_knots_double_t *k, double *t_inflection)
+{
+ double den;
+ int r = 0;
+
+ cairo_vector_t a = { -k->a.x + 3*k->b.x - 3*k->c.x + k->d.x, -k->a.y + 3*k->b.y - 3*k->c.y + k->d.y};
+ cairo_vector_t b = { 3*k->a.x - 6*k->b.x + 3*k->c.x, 3*k->a.y - 6*k->b.y + 3*k->c.y};
+ cairo_vector_t c = {-3*k->a.x + 3*k->b.x, -3*k->a.y + 3*k->b.y};
+ /* we don't need 'd'
+ cairo_vector_t d = { k->a.x, k->a.y};
+ */
+
+ /* we're solving for the roots of this equation:
+ dx/dt * d2y/dt2 - d2x/dt2 * dy/dt
+ == 6*(a.y*b.x - a.x*b.y)*t^2 +6*(a.y*c.x-a.x*c.y)*t + 2*(b.y*c.x - b.x*c.y)
+ */
+ if ((den = _cairo_vector_cross(a, b)) != 0.0) {
+ double invden = 1 / den;
+ double t_cusp = -0.5 * _cairo_vector_cross (a, c) * invden;
+ double t_inner = t_cusp*t_cusp - _cairo_vector_cross (b, c) * invden / 3;
+
+ if (t_inner > 0.0) {
+ double t_adj = sqrt (t_inner);
+ t_inflection[r] = t_cusp - t_adj;
+ if (t_inflection[r] > 0.0 && t_inflection[r] < 1.0)
+ r++;
+ t_inflection[r] = t_cusp + t_adj;
+ } else {
+ t_inflection[r] = t_cusp;
+ }
+ } else if ((den = _cairo_vector_cross(a, c)) != 0.0) {
+ t_inflection[r] = - _cairo_vector_cross (b, c) / (den * 3);
+ } else
+ return 0;
+
+ if (t_inflection[r] > 0.0 && t_inflection[r] < 1.0)
+ r++;
+
+ return r;
+}
+
+static void
+_cairo_spline_translate (cairo_spline_knots_double_t *dest, const cairo_spline_knots_double_t *src, const cairo_vector_t *t)
+{
+ dest->a = _cairo_vector_add (src->a, *t);
+ dest->b = _cairo_vector_add (src->b, *t);
+ dest->c = _cairo_vector_add (src->c, *t);
+ dest->d = _cairo_vector_add (src->d, *t);
+}
+
+static void
+_cairo_spline_scale (cairo_spline_knots_double_t *dest, const cairo_spline_knots_double_t *src, double scale)
+{
+ dest->a = _cairo_vector_scale (src->a, scale);
+ dest->b = _cairo_vector_scale (src->b, scale);
+ dest->c = _cairo_vector_scale (src->c, scale);
+ dest->d = _cairo_vector_scale (src->d, scale);
+}
+
+static void
+_cairo_spline_add (cairo_spline_knots_double_t *dest, const cairo_spline_knots_double_t *a, const cairo_spline_knots_double_t *b)
+{
+ dest->a = _cairo_vector_add (a->a, b->a);
+ dest->b = _cairo_vector_add (a->b, b->b);
+ dest->c = _cairo_vector_add (a->c, b->c);
+ dest->d = _cairo_vector_add (a->d, b->d);
+}
+
+/* Compute a spline approximation of the arc
+ centered at xc, yc from the angle a to the angle b
+
+ The angle between a and b should not be more than a
+ quarter circle (pi/2)
+
+ The approximation is similar to an approximation given in:
+ "Approximation of a cubic bezier curve by circular arcs and vice versa"
+ by Alekas Riškus. However that approximation becomes unstable when the
+ angle of the arc approaches 0.
+
+ This approximation is inspired by a discusion with Boris Zbarsky
+ and essentially just computes:
+
+ h = 4.0/3.0 * tan ((angle_B - angle_A) / 4.0);
+
+ without converting to polar coordinates.
+
+ A different way to do this is covered in "Approximation of a cubic bezier
+ curve by circular arcs and vice versa" by Alekas Riškus. However, the method
+ presented there doesn't handle arcs with angles close to 0 because it
+ divides by the perp dot product of the two angle vectors.
+
+ NB: start and end must be normalized and cannot be elements of the output knots structure
+*/
+static double
+_cairo_spline_arc (cairo_spline_knots_double_t *k, const cairo_vector_t *center, double radius, const cairo_vector_t *start, const cairo_vector_t *end)
+{
+ cairo_vector_t half = _cairo_vector_bisect (*start, *end);
+ double coshalf = _cairo_vector_dot (*start, half);
+ /* no need to use bisect since the angle between nstart and half is <= pi/2 < pi
+ * and the vector doesn't need to be normalized */
+ cairo_vector_t quart = _cairo_vector_add (*start, half);
+ /* h = |start| |quart| sin(angle) / (|start| |quart| cos (angle) = tan (angle) */
+ double h = (4./3.) * _cairo_vector_cross (*start, quart) / _cairo_vector_dot (*start, quart);
+
+ k->a = _cairo_vector_scale (*start, radius);
+ k->d = _cairo_vector_scale (*end, radius);
+ k->b = _cairo_vector_add (k->a, _cairo_vector_scale(_cairo_vector_perpendicular (k->a), +h));
+ k->c = _cairo_vector_add (k->d, _cairo_vector_scale(_cairo_vector_perpendicular (k->d), -h));
+ _cairo_spline_translate (k, k, center);
+
+ return radius * (7.-coshalf) * (2.+coshalf) * (2.+coshalf) / (27. * (1.+coshalf));
+}
+
+static void
+_cairo_spline_normal_sweep (const cairo_spline_knots_double_t *k, double *sweep)
+{
+ cairo_vector_t hodograph[3];
+ double dot;
+
+ hodograph[0] = _cairo_vector_normalize (_cairo_vector_sub (k->b, k->a));
+ hodograph[1] = _cairo_vector_normalize (_cairo_vector_sub (k->c, k->b));
+ hodograph[2] = _cairo_vector_normalize (_cairo_vector_sub (k->d, k->c));
+
+ /* the angle covered by each segment of the hodograph is < pi, thus its
+ * tangent uniquely determines it */
+ /* tan(a/2) = sin(a)/(1+cos(a)) */
+ dot = _cairo_vector_dot (hodograph[0], hodograph[1]);
+ if (dot != -1.0)
+ sweep[0] = _cairo_vector_cross (hodograph[0], hodograph[1]) / (1.0 + dot);
+ else
+ sweep[0] = HUGE_VAL;
+
+ dot = _cairo_vector_dot (hodograph[2], hodograph[1]);
+ if (dot != -1.0)
+ sweep[1] = _cairo_vector_cross (hodograph[2], hodograph[1]) / (1 + dot);
+ else
+ sweep[1] = HUGE_VAL;
+
+}
+
+/* split must be an array sufficient for the split curves (at most 3 splines, depending on
+ the inflection points */
+static int
+_cairo_spline_split_inflection (cairo_spline_knots_double_t *split, cairo_bool_t *cusp, const cairo_spline_knots_double_t *k)
+{
+ double t_inflect[2];
+ int r = 0, inflections = _cairo_spline_inflection_points (k, t_inflect);
+ cairo_spline_knots_double_t remaining = *k;
+
+ switch (inflections) {
+ case 2:
+ cusp[r] = _cairo_spline_de_casteljau (&remaining, split+r, &remaining, t_inflect[r]);
+ r++;
+
+ /* reparamaterize the second inflection point over the 'remaining' spline:
+ * (domain between inflection points)/(domain of remaining) */
+ t_inflect[1] = (t_inflect[1] - t_inflect[0]) / (1.0 - t_inflect[0]);
+
+ /* fall through */
+ case 1:
+ cusp[r] = _cairo_spline_de_casteljau (&remaining, split+r, &remaining, t_inflect[r]);
+ r++;
+
+ /* fall through */
+ case 0:
+ split[r++] = remaining;
+ return r;
+
+ default:
+ ASSERT_NOT_REACHED;
+ return 0;
+ }
+}
+
+#endif /* CAIRO_SPLINE_PRIVATE_H */
diff --git a/src/cairo-stroke-to-path.c b/src/cairo-stroke-to-path.c
index 8c0da821a..5b78b014b 100644
--- a/src/cairo-stroke-to-path.c
+++ b/src/cairo-stroke-to-path.c
@@ -90,151 +90,6 @@ static void close_path (path_output_t *path)
_cairo_path_fixed_new_sub_path (path->path);
}
-typedef point_t vector_t;
-
-static inline double
-dot (vector_t a, vector_t b)
-{
- return a.x * b.x + a.y * b.y;
-}
-
-static cairo_bool_t
-point_eq (point_t a, point_t b)
-{
- return a.x == b.x && a.y == b.y;
-}
-
-static inline vector_t
-perp (vector_t v)
-{
- vector_t p = {-v.y, v.x};
- return p;
-}
-
-static inline vector_t
-flip (vector_t a)
-{
- vector_t ar;
- ar.x = -a.x;
- ar.y = -a.y;
- return ar;
-}
-
-/* given a normal rotate the vector 90 degrees to the right clockwise
- * This function has a period of 4. e.g. swap(swap(swap(swap(x) == x */
-static vector_t
-swap (vector_t a)
-{
- vector_t ar;
- /* one of these needs to be negative. We choose a.x so that we rotate to the right instead of negating */
- ar.x = a.y;
- ar.y = -a.x;
- return ar;
-}
-
-static vector_t
-unperp (vector_t a)
-{
- return swap (a);
-}
-
-/* Compute a spline approximation of the arc
- centered at xc, yc from the angle a to the angle b
-
- The angle between a and b should not be more than a
- quarter circle (pi/2)
-
- The approximation is similar to an approximation given in:
- "Approximation of a cubic bezier curve by circular arcs and vice versa"
- by Alekas Riškus. However that approximation becomes unstable when the
- angle of the arc approaches 0.
-
- This approximation is inspired by a discusion with Boris Zbarsky
- and essentially just computes:
-
- h = 4.0/3.0 * tan ((angle_B - angle_A) / 4.0);
-
- without converting to polar coordinates.
-
- A different way to do this is covered in "Approximation of a cubic bezier
- curve by circular arcs and vice versa" by Alekas Riškus. However, the method
- presented there doesn't handle arcs with angles close to 0 because it
- divides by the perp dot product of the two angle vectors.
- */
-
-static void
-arc_segment (path_output_t *path,
- double xc,
- double yc,
- double radius,
- vector_t a,
- vector_t b)
-{
- double r_sin_A, r_cos_A;
- double r_sin_B, r_cos_B;
- double h;
- vector_t mid, mid2;
- double l;
-
- r_sin_A = radius * a.y;
- r_cos_A = radius * a.x;
- r_sin_B = radius * b.y;
- r_cos_B = radius * b.x;
-
- /* bisect the angle between 'a' and 'b' with 'mid' */
- mid.x = a.x + b.x;
- mid.y = a.y + b.y;
- l = sqrt(mid.x*mid.x + mid.y*mid.y);
- mid.x /= l;
- mid.y /= l;
-
- /* bisect the angle between 'a' and 'mid' with 'mid2' this is parallel to a
- * line with angle (B - A)/4 */
- mid2.x = a.x + mid.x;
- mid2.y = a.y + mid.y;
-
- h = (4./3.)*dot(perp(a), mid2)/dot(a, mid2);
- //h = (4./3.)*(-a.y*mid.x + a.x*mid.y)/(a.x*mid2.x + a.y*mid2.y);
-
- curve_to (path,
- xc + r_cos_A - h * r_sin_A,
- yc + r_sin_A + h * r_cos_A,
- xc + r_cos_B + h * r_sin_B,
- yc + r_sin_B - h * r_cos_B,
- xc + r_cos_B,
- yc + r_sin_B);
-}
-
-
-static vector_t
-bisect (vector_t a, vector_t b)
-{
- vector_t mid;
- double len;
-
- if (dot (a, b) >= 0) {
- /* if the angle between a and b is accute, then we can
- * just add the vectors and normalize */
- mid.x = a.x + b.x;
- mid.y = a.y + b.y;
- } else {
- /* otherwise, we can flip a, add it
- * and then use the perpendicular of the result */
- a = flip (a);
- mid.x = a.x + b.x;
- mid.y = a.y + b.y;
- mid = perp (mid);
- }
-
- /* normalize */
- /* because we assume that 'a' and 'b' are normalized, we can use
- * sqrt instead of hypot because the range of mid is limited */
- len = sqrt (mid.x*mid.x + mid.y*mid.y);
- mid.x /= len;
- mid.y /= len;
- return mid;
-}
-
static void
arc (path_output_t *path, double xc, double yc, double radius, vector_t a, vector_t b)
{
@@ -351,50 +206,13 @@ compute_normal (point_t p0, point_t p1)
-/* computes the outgoing normal for curve given by a,b,c,d
- * returning the original normal if the curve is point */
-static vector_t
-curve_normal (vector_t original_normal, point_t a, point_t b, point_t c, point_t d)
-{
- vector_t n = original_normal;
- if (point_eq (c, d)) {
- if (point_eq (b, c)) {
- if (!point_eq (a, b)) {
- n = compute_normal (a, b);
- }
- } else {
- n = compute_normal (b, c);
- }
- } else {
- n = compute_normal (c, d);
- }
- return n;
-}
-static void ensure_finite (knots_t k)
-{
-
- if (!(
- isfinite(k.a.x) &&
- isfinite(k.a.y) &&
- isfinite(k.b.x) &&
- isfinite(k.b.y) &&
- isfinite(k.c.x) &&
- isfinite(k.c.y) &&
- isfinite(k.d.x) &&
- isfinite(k.d.y))) {
- print_knot(k);
- assert(0);
- }
-}
-
-
static void offset_curve_to(void *closure, knots_t c)
{
path_output_t *path = closure;
ensure_finite(c);
//XXX: it would be nice if we didn't have to check this for every point
if (!path->has_current_point)
- line_to(path, c.a.x, c.a.y);
+ move_to(path, c.a.x, c.a.y);
curve_to(path, c.b.x, c.b.y, c.c.x, c.c.y, c.d.x, c.d.y);
}
@@ -429,6 +247,7 @@ line_intersection (point_t A, vector_t a_perp, point_t B, vector_t b_perp)
vector_t a = unperp(a_perp);
vector_t c = {B.x - A.x, B.y - A.y};
denom = dot(b_perp, a);
+ /* CHECK ME!!! */
if (denom == 0.0) {
assert(0 && "TROUBLE");
}
@@ -1247,6 +1066,7 @@ dash_subpath (dash_point_t start_point, path_output_t *path_out, path_ittr *inde
closed_p1 = c2p (ctm_i, i.points[0]);
closed_p2 = end_point.point;
closed_n1 = n1;
+ /* CHECK ME!!! */
closed_n2 = compute_normal (closed_p1, closed_p2);
closed_n3 = compute_normal (closed_p2, c2p(ctm_i, ittr_next (next_i).points[0]));
join_segment_line (path_out, style, closed_p1, closed_n1, closed_n2);
diff --git a/src/cairo-vector-private.h b/src/cairo-vector-private.h
new file mode 100644
index 000000000..13a1eeee7
--- /dev/null
+++ b/src/cairo-vector-private.h
@@ -0,0 +1,134 @@
+/* -*- Mode: c; tab-width: 8; c-basic-offset: 4; indent-tabs-mode: t; -*- */
+#ifndef CAIRO_VECTOR_PRIVATE_H
+#define CAIRO_VECTOR_PRIVATE_H
+
+#include <math.h>
+
+#if _XOPEN_SOURCE >= 600 || defined (_ISOC99_SOURCE)
+#define ISFINITE(x) isfinite (x)
+#else
+#define ISFINITE(x) ((x) * (x) >= 0.) /* check for NaNs */
+#endif
+
+#include "cairoint.h"
+
+typedef cairo_point_double_t cairo_vector_t;
+
+static const cairo_vector_t _cairo_vector_zero = { 0, 0 };
+
+static inline cairo_bool_t
+_cairo_vector_equals (cairo_vector_t a, cairo_vector_t b)
+{
+ return a.x == b.x && a.y == b.y; /* TODO: fixed point precision? */
+}
+
+static inline cairo_bool_t
+_cairo_vector_iszero (cairo_vector_t a)
+{
+ return _cairo_vector_equals (a, _cairo_vector_zero);
+}
+
+static inline cairo_bool_t
+_cairo_vector_isfinite (cairo_vector_t a)
+{
+ return ISFINITE (a.x) && ISFINITE (a.y);
+}
+
+static inline double
+_cairo_vector_cross (cairo_vector_t a, cairo_vector_t b)
+{
+ return a.x * b.y - a.y * b.x;
+}
+
+static inline double
+_cairo_vector_dot (cairo_vector_t a, cairo_vector_t b)
+{
+ return a.x * b.x + a.y * b.y;
+}
+
+static inline double
+_cairo_vector_length (cairo_vector_t v)
+{
+ return hypot (v.x, v.y);
+}
+
+static inline cairo_vector_t
+_cairo_vector_opposite (cairo_vector_t p)
+{
+ cairo_vector_t r = { -p.x, -p.y };
+ return r;
+}
+
+static inline cairo_vector_t
+_cairo_vector_perpendicular (cairo_vector_t p)
+{
+ cairo_vector_t r = { -p.y, p.x};
+ return r;
+}
+
+static inline cairo_vector_t
+_cairo_vector_add (cairo_vector_t a, cairo_vector_t b)
+{
+ cairo_vector_t r = { a.x + b.x, a.y + b.y };
+ return r;
+}
+
+static inline cairo_vector_t
+_cairo_vector_sub (cairo_vector_t a, cairo_vector_t b)
+{
+ cairo_vector_t r = { a.x - b.x, a.y - b.y };
+ return r;
+}
+
+static inline cairo_vector_t
+_cairo_vector_scale (cairo_vector_t v, double scale)
+{
+ cairo_vector_t r = { v.x * scale, v.y * scale };
+ return r;
+}
+
+static inline cairo_vector_t
+_cairo_vector_normalize (cairo_vector_t v)
+{
+ assert (!_cairo_vector_iszero (v));
+ return _cairo_vector_scale (v, 1 / _cairo_vector_length (v));
+}
+
+static inline cairo_vector_t
+_cairo_vector_lerp (cairo_vector_t start, cairo_vector_t end, double t)
+{
+ return _cairo_vector_add (_cairo_vector_scale (start, 1-t), _cairo_vector_scale (end, t));
+}
+
+/* a and b must be normalized */
+static cairo_vector_t
+_cairo_vector_bisect_ccw (cairo_vector_t a, cairo_vector_t b)
+{
+ cairo_vector_t mid = _cairo_vector_add (a, b);
+
+ if (_cairo_vector_iszero (mid))
+ return _cairo_vector_perpendicular (a);
+
+ mid = _cairo_vector_normalize (mid);
+ /* If the angle between a and b is a reflex angle (> pi),
+ * then the direction of the bisector is opposite to the
+ * one we just computed. */
+ if (_cairo_vector_cross (a, b) < 0)
+ mid = _cairo_vector_opposite (mid);
+
+ return mid;
+}
+
+/* a and b must be normalized */
+static cairo_vector_t
+_cairo_vector_bisect (cairo_vector_t a, cairo_vector_t b)
+{
+ cairo_vector_t mid = _cairo_vector_add (a, b);
+
+ if (_cairo_vector_iszero (mid))
+ return _cairo_vector_perpendicular (a);
+ else
+ return _cairo_vector_normalize (mid);
+}
+
+#endif /* CAIRO_VECTOR_PRIVATE_H */