diff options
author | Andrea Canciani <ranma42@gmail.com> | 2010-06-28 16:47:52 +0200 |
---|---|---|
committer | Andrea Canciani <ranma42@gmail.com> | 2010-10-23 10:52:47 +0200 |
commit | 0e02921b8b0f674d8570fd4b929909b25d48fa02 (patch) | |
tree | df72cc1221d6f2f9e2b542746f89eed1035a0cd9 | |
parent | e3dd9f3856201c6b38d747971a7d5dfd18a7f635 (diff) |
-rw-r--r-- | src/cairo-spline-offset-jm.c | 488 | ||||
-rw-r--r-- | src/cairo-spline-offset.c | 992 | ||||
-rw-r--r-- | src/cairo-spline-offset.h | 20 | ||||
-rw-r--r-- | src/cairo-spline-private.h | 292 | ||||
-rw-r--r-- | src/cairo-stroke-to-path.c | 186 | ||||
-rw-r--r-- | src/cairo-vector-private.h | 134 |
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 */ |