src/cairo-spline.c

   1 /* cairo - a vector graphics library with display and print output
   2  *
   3  * Copyright © 2002 University of Southern California
   4  *
   5  * This library is free software; you can redistribute it and/or
   6  * modify it either under the terms of the GNU Lesser General Public
   7  * License version 2.1 as published by the Free Software Foundation
   8  * (the "LGPL") or, at your option, under the terms of the Mozilla
   9  * Public License Version 1.1 (the "MPL"). If you do not alter this
  10  * notice, a recipient may use your version of this file under either
  11  * the MPL or the LGPL.
  12  *
  13  * You should have received a copy of the LGPL along with this library
  14  * in the file COPYING-LGPL-2.1; if not, write to the Free Software
  15  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
  16  * You should have received a copy of the MPL along with this library
  17  * in the file COPYING-MPL-1.1
  18  *
  19  * The contents of this file are subject to the Mozilla Public License
  20  * Version 1.1 (the "License"); you may not use this file except in
  21  * compliance with the License. You may obtain a copy of the License at
  22  * http://www.mozilla.org/MPL/
  23  *
  24  * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
  25  * OF ANY KIND, either express or implied. See the LGPL or the MPL for
  26  * the specific language governing rights and limitations.
  27  *
  28  * The Original Code is the cairo graphics library.
  29  *
  30  * The Initial Developer of the Original Code is University of Southern
  31  * California.
  32  *
  33  * Contributor(s):
  34  *      Carl D. Worth <cworth@cworth.org>
  35  */
  36
  37 #include "cairoint.h"
  38
  39 cairo_bool_t
  40 _cairo_spline_init (cairo_spline_t *spline,
  41                     cairo_spline_add_point_func_t add_point_func,
  42                     void *closure,
  43                     const cairo_point_t *a, const cairo_point_t *b,
  44                     const cairo_point_t *c, const cairo_point_t *d)
  45 {
  46     spline->add_point_func = add_point_func;
  47     spline->closure = closure;
  48
  49     spline->knots.a = *a;
  50     spline->knots.b = *b;
  51     spline->knots.c = *c;
  52     spline->knots.d = *d;
  53
  54     if (a->x != b->x || a->y != b->y)
  55         _cairo_slope_init (&spline->initial_slope, &spline->knots.a, &spline->knots.b);
  56     else if (a->x != c->x || a->y != c->y)
  57         _cairo_slope_init (&spline->initial_slope, &spline->knots.a, &spline->knots.c);
  58     else if (a->x != d->x || a->y != d->y)
  59         _cairo_slope_init (&spline->initial_slope, &spline->knots.a, &spline->knots.d);
  60     else
  61         return FALSE;
  62
  63     if (c->x != d->x || c->y != d->y)
  64         _cairo_slope_init (&spline->final_slope, &spline->knots.c, &spline->knots.d);
  65     else if (b->x != d->x || b->y != d->y)
  66         _cairo_slope_init (&spline->final_slope, &spline->knots.b, &spline->knots.d);
  67     else
  68         _cairo_slope_init (&spline->final_slope, &spline->knots.a, &spline->knots.d);
  69
  70     return TRUE;
  71 }
  72
  73 static cairo_status_t
  74 _cairo_spline_add_point (cairo_spline_t *spline, cairo_point_t *point)
  75 {
  76     cairo_point_t *prev;
  77
  78     prev = &spline->last_point;
  79     if (prev->x == point->x && prev->y == point->y)
  80         return CAIRO_STATUS_SUCCESS;
  81
  82     spline->last_point = *point;
  83     return spline->add_point_func (spline->closure, point);
  84 }
  85
  86 static void
  87 _lerp_half (const cairo_point_t *a, const cairo_point_t *b, cairo_point_t *result)
  88 {
  89     result->x = a->x + ((b->x - a->x) >> 1);
  90     result->y = a->y + ((b->y - a->y) >> 1);
  91 }
  92
  93 static void
  94 _de_casteljau (cairo_spline_knots_t *s1, cairo_spline_knots_t *s2)
  95 {
  96     cairo_point_t ab, bc, cd;
  97     cairo_point_t abbc, bccd;
  98     cairo_point_t final;
  99
 100     _lerp_half (&s1->a, &s1->b, &ab);
 101     _lerp_half (&s1->b, &s1->c, &bc);
 102     _lerp_half (&s1->c, &s1->d, &cd);
 103     _lerp_half (&ab, &bc, &abbc);
 104     _lerp_half (&bc, &cd, &bccd);
 105     _lerp_half (&abbc, &bccd, &final);
 106
 107     s2->a = final;
 108     s2->b = bccd;
 109     s2->c = cd;
 110     s2->d = s1->d;
 111
 112     s1->b = ab;
 113     s1->c = abbc;
 114     s1->d = final;
 115 }
 116
 117 /* Return an upper bound on the error (squared) that could result from
 118  * approximating a spline as a line segment connecting the two endpoints. */
 119 static double
 120 _cairo_spline_error_squared (const cairo_spline_knots_t *knots)
 121 {
 122     double bdx, bdy, berr;
 123     double cdx, cdy, cerr;
 124
 125     /* Intersection point (px):
 126      *      px = p1 + u(p2 - p1)
 127      *      (p - px) ∙ (p2 - p1) = 0
 128      * Thus:
 129      *      u = ((p - p1) ∙ (p2 - p1)) / ∥p2 - p1∥²;
 130      */
 131     bdx = _cairo_fixed_to_double (knots->b.x - knots->a.x);
 132     bdy = _cairo_fixed_to_double (knots->b.y - knots->a.y);
 133
 134     cdx = _cairo_fixed_to_double (knots->c.x - knots->a.x);
 135     cdy = _cairo_fixed_to_double (knots->c.y - knots->a.y);
 136
 137     if (knots->a.x != knots->d.x || knots->a.y != knots->d.y) {
 138         double dx, dy, u, v;
 139
 140         dx = _cairo_fixed_to_double (knots->d.x - knots->a.x);
 141         dy = _cairo_fixed_to_double (knots->d.y - knots->a.y);
 142          v = dx * dx + dy * dy;
 143
 144         u = bdx * dx + bdy * dy;
 145         if (u <= 0) {
 146             /* bdx -= 0;
 147              * bdy -= 0;
 148              */
 149         } else if (u >= v) {
 150             bdx -= dx;
 151             bdy -= dy;
 152         } else {
 153             bdx -= u/v * dx;
 154             bdy -= u/v * dy;
 155         }
 156
 157         u = cdx * dx + cdy * dy;
 158         if (u <= 0) {
 159             /* cdx -= 0;
 160              * cdy -= 0;
 161              */
 162         } else if (u >= v) {
 163             cdx -= dx;
 164             cdy -= dy;
 165         } else {
 166             cdx -= u/v * dx;
 167             cdy -= u/v * dy;
 168         }
 169     }
 170
 171     berr = bdx * bdx + bdy * bdy;
 172     cerr = cdx * cdx + cdy * cdy;
 173     if (berr > cerr)
 174         return berr;
 175     else
 176         return cerr;
 177 }
 178
 179 static cairo_status_t
 180 _cairo_spline_decompose_into (cairo_spline_knots_t *s1, double tolerance_squared, cairo_spline_t *result)
 181 {
 182     cairo_spline_knots_t s2;
 183     cairo_status_t status;
 184
 185     if (_cairo_spline_error_squared (s1) < tolerance_squared)
 186         return _cairo_spline_add_point (result, &s1->a);
 187
 188     _de_casteljau (s1, &s2);
 189
 190     status = _cairo_spline_decompose_into (s1, tolerance_squared, result);
 191     if (unlikely (status))
 192         return status;
 193
 194     return _cairo_spline_decompose_into (&s2, tolerance_squared, result);
 195 }
 196
 197 cairo_status_t
 198 _cairo_spline_decompose (cairo_spline_t *spline, double tolerance)
 199 {
 200     cairo_spline_knots_t s1;
 201     cairo_status_t status;
 202
 203     s1 = spline->knots;
 204     spline->last_point = s1.a;
 205     status = _cairo_spline_decompose_into (&s1, tolerance * tolerance, spline);
 206     if (unlikely (status))
 207         return status;
 208
 209     return _cairo_spline_add_point (spline, &spline->knots.d);
 210 }
 211
 212 /* Note: this function is only good for computing bounds in device space. */
 213 cairo_status_t
 214 _cairo_spline_bound (cairo_spline_add_point_func_t add_point_func,
 215                      void *closure,
 216                      const cairo_point_t *p0, const cairo_point_t *p1,
 217                      const cairo_point_t *p2, const cairo_point_t *p3)
 218 {
 219     double x0, x1, x2, x3;
 220     double y0, y1, y2, y3;
 221     double a, b, c;
 222     double t[4];
 223     int t_num = 0, i;
 224     cairo_status_t status;
 225
 226     x0 = _cairo_fixed_to_double (p0->x);
 227     y0 = _cairo_fixed_to_double (p0->y);
 228     x1 = _cairo_fixed_to_double (p1->x);
 229     y1 = _cairo_fixed_to_double (p1->y);
 230     x2 = _cairo_fixed_to_double (p2->x);
 231     y2 = _cairo_fixed_to_double (p2->y);
 232     x3 = _cairo_fixed_to_double (p3->x);
 233     y3 = _cairo_fixed_to_double (p3->y);
 234
 235     /* The spline can be written as a polynomial of the four points:
 236      *
 237      *   (1-t)³p0 + 3t(1-t)²p1 + 3t²(1-t)p2 + t³p3
 238      *
 239      * for 0≤t≤1.  Now, the X and Y components of the spline follow the
 240      * same polynomial but with x and y replaced for p.  To find the
 241      * bounds of the spline, we just need to find the X and Y bounds.
 242      * To find the bound, we take the derivative and equal it to zero,
 243      * and solve to find the t's that give the extreme points.
 244      *
 245      * Here is the derivative of the curve, sorted on t:
 246      *
 247      *   3t²(-p0+3p1-3p2+p3) + 2t(3p0-6p1+3p2) -3p0+3p1
 248      *
 249      * Let:
 250      *
 251      *   a = -p0+3p1-3p2+p3
 252      *   b =  p0-2p1+p2
 253      *   c = -p0+p1
 254      *
 255      * Gives:
 256      *
 257      *   a.t² + 2b.t + c = 0
 258      *
 259      * With:
 260      *
 261      *   delta = b*b - a*c
 262      *
 263      * the extreme points are at -c/2b if a is zero, at (-b±√delta)/a if
 264      * delta is positive, and at -b/a if delta is zero.
 265      */
 266
 267 #define ADD(t0) \
 268     { \
 269         double _t0 = (t0); \
 270         if (0 < _t0 && _t0 < 1) \
 271             t[t_num++] = _t0; \
 272     }
 273
 274 #define FIND_EXTREMES(a,b,c) \
 275     { \
 276         if (a == 0) { \
 277             if (b != 0) \
 278                 ADD (-c / (2*b)); \
 279         } else { \
 280             double b2 = b * b; \
 281             double delta = b2 - a * c; \
 282             if (delta > 0) { \
 283                 cairo_bool_t feasible; \
 284                 double _2ab = 2 * a * b; \
 285                 /* We are only interested in solutions t that satisfy 0<t<1 \
 286                  * here.  We do some checks to avoid sqrt if the solutions \
 287                  * are not in that range.  The checks can be derived from: \
 288                  * \
 289                  *   0 < (-b±√delta)/a < 1 \
 290                  */ \
 291                 if (_2ab >= 0) \
 292                     feasible = delta > b2 && delta < a*a + b2 + _2ab; \
 293                 else if (-b / a >= 1) \
 294                     feasible = delta < b2 && delta > a*a + b2 + _2ab; \
 295                 else \
 296                     feasible = delta < b2 || delta < a*a + b2 + _2ab; \
 297                 \
 298                 if (unlikely (feasible)) { \
 299                     double sqrt_delta = sqrt (delta); \
 300                     ADD ((-b - sqrt_delta) / a); \
 301                     ADD ((-b + sqrt_delta) / a); \
 302                 } \
 303             } else if (delta == 0) { \
 304                 ADD (-b / a); \
 305             } \
 306         } \
 307     }
 308
 309     /* Find X extremes */
 310     a = -x0 + 3*x1 - 3*x2 + x3;
 311     b =  x0 - 2*x1 + x2;
 312     c = -x0 + x1;
 313     FIND_EXTREMES (a, b, c);
 314
 315     /* Find Y extremes */
 316     a = -y0 + 3*y1 - 3*y2 + y3;
 317     b =  y0 - 2*y1 + y2;
 318     c = -y0 + y1;
 319     FIND_EXTREMES (a, b, c);
 320
 321     status = add_point_func (closure, p0);
 322     if (unlikely (status))
 323         return status;
 324
 325     for (i = 0; i < t_num; i++) {
 326         cairo_point_t p;
 327         double x, y;
 328         double t_1_0, t_0_1;
 329         double t_2_0, t_0_2;
 330         double t_3_0, t_2_1_3, t_1_2_3, t_0_3;
 331
 332         t_1_0 = t[i];          /*      t  */
 333         t_0_1 = 1 - t_1_0;     /* (1 - t) */
 334
 335         t_2_0 = t_1_0 * t_1_0; /*      t  *      t  */
 336         t_0_2 = t_0_1 * t_0_1; /* (1 - t) * (1 - t) */
 337
 338         t_3_0   = t_2_0 * t_1_0;     /*      t  *      t  *      t      */
 339         t_2_1_3 = t_2_0 * t_0_1 * 3; /*      t  *      t  * (1 - t) * 3 */
 340         t_1_2_3 = t_1_0 * t_0_2 * 3; /*      t  * (1 - t) * (1 - t) * 3 */
 341         t_0_3   = t_0_1 * t_0_2;     /* (1 - t) * (1 - t) * (1 - t)     */
 342
 343         /* Bezier polynomial */
 344         x = x0 * t_0_3
 345           + x1 * t_1_2_3
 346           + x2 * t_2_1_3
 347           + x3 * t_3_0;
 348         y = y0 * t_0_3
 349           + y1 * t_1_2_3
 350           + y2 * t_2_1_3
 351           + y3 * t_3_0;
 352
 353         p.x = _cairo_fixed_from_double (x);
 354         p.y = _cairo_fixed_from_double (y);
 355         status = add_point_func (closure, &p);
 356         if (unlikely (status))
 357             return status;
 358     }
 359
 360     return add_point_func (closure, p3);
 361 }