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., 51 Franklin Street, Suite 500, Boston, MA 02110-1335, 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 #include "cairo-box-inline.h"
  40 #include "cairo-slope-private.h"
  41
  42 cairo_bool_t
  43 _cairo_spline_intersects (const cairo_point_t *a,
  44                           const cairo_point_t *b,
  45                           const cairo_point_t *c,
  46                           const cairo_point_t *d,
  47                           const cairo_box_t *box)
  48 {
  49     cairo_box_t bounds;
  50
  51     if (_cairo_box_contains_point (box, a) ||
  52         _cairo_box_contains_point (box, b) ||
  53         _cairo_box_contains_point (box, c) ||
  54         _cairo_box_contains_point (box, d))
  55     {
  56         return TRUE;
  57     }
  58
  59     bounds.p2 = bounds.p1 = *a;
  60     _cairo_box_add_point (&bounds, b);
  61     _cairo_box_add_point (&bounds, c);
  62     _cairo_box_add_point (&bounds, d);
  63
  64     if (bounds.p2.x <= box->p1.x || bounds.p1.x >= box->p2.x ||
  65         bounds.p2.y <= box->p1.y || bounds.p1.y >= box->p2.y)
  66     {
  67         return FALSE;
  68     }
  69
  70 #if 0 /* worth refining? */
  71     bounds.p2 = bounds.p1 = *a;
  72     _cairo_box_add_curve_to (&bounds, b, c, d);
  73     if (bounds.p2.x <= box->p1.x || bounds.p1.x >= box->p2.x ||
  74         bounds.p2.y <= box->p1.y || bounds.p1.y >= box->p2.y)
  75     {
  76         return FALSE;
  77     }
  78 #endif
  79
  80     return TRUE;
  81 }
  82
  83 cairo_bool_t
  84 _cairo_spline_init (cairo_spline_t *spline,
  85                     cairo_spline_add_point_func_t add_point_func,
  86                     void *closure,
  87                     const cairo_point_t *a, const cairo_point_t *b,
  88                     const cairo_point_t *c, const cairo_point_t *d)
  89 {
  90     /* If both tangents are zero, this is just a straight line */
  91     if (a->x == b->x && a->y == b->y && c->x == d->x && c->y == d->y)
  92         return FALSE;
  93
  94     spline->add_point_func = add_point_func;
  95     spline->closure = closure;
  96
  97     spline->knots.a = *a;
  98     spline->knots.b = *b;
  99     spline->knots.c = *c;
 100     spline->knots.d = *d;
 101
 102     if (a->x != b->x || a->y != b->y)
 103         _cairo_slope_init (&spline->initial_slope, &spline->knots.a, &spline->knots.b);
 104     else if (a->x != c->x || a->y != c->y)
 105         _cairo_slope_init (&spline->initial_slope, &spline->knots.a, &spline->knots.c);
 106     else if (a->x != d->x || a->y != d->y)
 107         _cairo_slope_init (&spline->initial_slope, &spline->knots.a, &spline->knots.d);
 108     else
 109         return FALSE;
 110
 111     if (c->x != d->x || c->y != d->y)
 112         _cairo_slope_init (&spline->final_slope, &spline->knots.c, &spline->knots.d);
 113     else if (b->x != d->x || b->y != d->y)
 114         _cairo_slope_init (&spline->final_slope, &spline->knots.b, &spline->knots.d);
 115     else
 116         return FALSE; /* just treat this as a straight-line from a -> d */
 117
 118     /* XXX if the initial, final and vector are all equal, this is just a line */
 119
 120     return TRUE;
 121 }
 122
 123 static cairo_status_t
 124 _cairo_spline_add_point (cairo_spline_t *spline,
 125                          const cairo_point_t *point,
 126                          const cairo_point_t *knot)
 127 {
 128     cairo_point_t *prev;
 129     cairo_slope_t slope;
 130
 131     prev = &spline->last_point;
 132     if (prev->x == point->x && prev->y == point->y)
 133         return CAIRO_STATUS_SUCCESS;
 134
 135     _cairo_slope_init (&slope, point, knot);
 136
 137     spline->last_point = *point;
 138     return spline->add_point_func (spline->closure, point, &slope);
 139 }
 140
 141 static void
 142 _lerp_half (const cairo_point_t *a, const cairo_point_t *b, cairo_point_t *result)
 143 {
 144     result->x = a->x + ((b->x - a->x) >> 1);
 145     result->y = a->y + ((b->y - a->y) >> 1);
 146 }
 147
 148 static void
 149 _de_casteljau (cairo_spline_knots_t *s1, cairo_spline_knots_t *s2)
 150 {
 151     cairo_point_t ab, bc, cd;
 152     cairo_point_t abbc, bccd;
 153     cairo_point_t final;
 154
 155     _lerp_half (&s1->a, &s1->b, &ab);
 156     _lerp_half (&s1->b, &s1->c, &bc);
 157     _lerp_half (&s1->c, &s1->d, &cd);
 158     _lerp_half (&ab, &bc, &abbc);
 159     _lerp_half (&bc, &cd, &bccd);
 160     _lerp_half (&abbc, &bccd, &final);
 161
 162     s2->a = final;
 163     s2->b = bccd;
 164     s2->c = cd;
 165     s2->d = s1->d;
 166
 167     s1->b = ab;
 168     s1->c = abbc;
 169     s1->d = final;
 170 }
 171
 172 /* Return an upper bound on the error (squared) that could result from
 173  * approximating a spline as a line segment connecting the two endpoints. */
 174 static double
 175 _cairo_spline_error_squared (const cairo_spline_knots_t *knots)
 176 {
 177     double bdx, bdy, berr;
 178     double cdx, cdy, cerr;
 179
 180     /* We are going to compute the distance (squared) between each of the the b
 181      * and c control points and the segment a-b. The maximum of these two
 182      * distances will be our approximation error. */
 183
 184     bdx = _cairo_fixed_to_double (knots->b.x - knots->a.x);
 185     bdy = _cairo_fixed_to_double (knots->b.y - knots->a.y);
 186
 187     cdx = _cairo_fixed_to_double (knots->c.x - knots->a.x);
 188     cdy = _cairo_fixed_to_double (knots->c.y - knots->a.y);
 189
 190     if (knots->a.x != knots->d.x || knots->a.y != knots->d.y) {
 191         /* Intersection point (px):
 192          *     px = p1 + u(p2 - p1)
 193          *     (p - px) ∙ (p2 - p1) = 0
 194          * Thus:
 195          *     u = ((p - p1) ∙ (p2 - p1)) / ∥p2 - p1∥²;
 196          */
 197
 198         double dx, dy, u, v;
 199
 200         dx = _cairo_fixed_to_double (knots->d.x - knots->a.x);
 201         dy = _cairo_fixed_to_double (knots->d.y - knots->a.y);
 202          v = dx * dx + dy * dy;
 203
 204         u = bdx * dx + bdy * dy;
 205         if (u <= 0) {
 206             /* bdx -= 0;
 207              * bdy -= 0;
 208              */
 209         } else if (u >= v) {
 210             bdx -= dx;
 211             bdy -= dy;
 212         } else {
 213             bdx -= u/v * dx;
 214             bdy -= u/v * dy;
 215         }
 216
 217         u = cdx * dx + cdy * dy;
 218         if (u <= 0) {
 219             /* cdx -= 0;
 220              * cdy -= 0;
 221              */
 222         } else if (u >= v) {
 223             cdx -= dx;
 224             cdy -= dy;
 225         } else {
 226             cdx -= u/v * dx;
 227             cdy -= u/v * dy;
 228         }
 229     }
 230
 231     berr = bdx * bdx + bdy * bdy;
 232     cerr = cdx * cdx + cdy * cdy;
 233     if (berr > cerr)
 234         return berr;
 235     else
 236         return cerr;
 237 }
 238
 239 static cairo_status_t
 240 _cairo_spline_decompose_into (cairo_spline_knots_t *s1,
 241                               double tolerance_squared,
 242                               cairo_spline_t *result)
 243 {
 244     cairo_spline_knots_t s2;
 245     cairo_status_t status;
 246
 247     if (_cairo_spline_error_squared (s1) < tolerance_squared)
 248         return _cairo_spline_add_point (result, &s1->a, &s1->b);
 249
 250     _de_casteljau (s1, &s2);
 251
 252     status = _cairo_spline_decompose_into (s1, tolerance_squared, result);
 253     if (unlikely (status))
 254         return status;
 255
 256     return _cairo_spline_decompose_into (&s2, tolerance_squared, result);
 257 }
 258
 259 cairo_status_t
 260 _cairo_spline_decompose (cairo_spline_t *spline, double tolerance)
 261 {
 262     cairo_spline_knots_t s1;
 263     cairo_status_t status;
 264
 265     /* this is the entry point for spline decompose, we adjust the
 266        final_slope if b, c, d are very close
 267      */
 268     if (_cairo_spline_error_squared (&spline->knots) < tolerance * tolerance)
 269         spline->final_slope = spline->initial_slope;
 270
 271     s1 = spline->knots;
 272     spline->last_point = s1.a;
 273     status = _cairo_spline_decompose_into (&s1, tolerance * tolerance, spline);
 274     if (unlikely (status))
 275         return status;
 276
 277     return spline->add_point_func (spline->closure,
 278                                    &spline->knots.d, &spline->final_slope);
 279 }
 280
 281 /* Note: this function is only good for computing bounds in device space. */
 282 cairo_status_t
 283 _cairo_spline_bound (cairo_spline_add_point_func_t add_point_func,
 284                      void *closure,
 285                      const cairo_point_t *p0, const cairo_point_t *p1,
 286                      const cairo_point_t *p2, const cairo_point_t *p3)
 287 {
 288     double x0, x1, x2, x3;
 289     double y0, y1, y2, y3;
 290     double a, b, c;
 291     double t[4];
 292     int t_num = 0, i;
 293     cairo_status_t status;
 294
 295     x0 = _cairo_fixed_to_double (p0->x);
 296     y0 = _cairo_fixed_to_double (p0->y);
 297     x1 = _cairo_fixed_to_double (p1->x);
 298     y1 = _cairo_fixed_to_double (p1->y);
 299     x2 = _cairo_fixed_to_double (p2->x);
 300     y2 = _cairo_fixed_to_double (p2->y);
 301     x3 = _cairo_fixed_to_double (p3->x);
 302     y3 = _cairo_fixed_to_double (p3->y);
 303
 304     /* The spline can be written as a polynomial of the four points:
 305      *
 306      *   (1-t)³p0 + 3t(1-t)²p1 + 3t²(1-t)p2 + t³p3
 307      *
 308      * for 0≤t≤1.  Now, the X and Y components of the spline follow the
 309      * same polynomial but with x and y replaced for p.  To find the
 310      * bounds of the spline, we just need to find the X and Y bounds.
 311      * To find the bound, we take the derivative and equal it to zero,
 312      * and solve to find the t's that give the extreme points.
 313      *
 314      * Here is the derivative of the curve, sorted on t:
 315      *
 316      *   3t²(-p0+3p1-3p2+p3) + 2t(3p0-6p1+3p2) -3p0+3p1
 317      *
 318      * Let:
 319      *
 320      *   a = -p0+3p1-3p2+p3
 321      *   b =  p0-2p1+p2
 322      *   c = -p0+p1
 323      *
 324      * Gives:
 325      *
 326      *   a.t² + 2b.t + c = 0
 327      *
 328      * With:
 329      *
 330      *   delta = b*b - a*c
 331      *
 332      * the extreme points are at -c/2b if a is zero, at (-b±√delta)/a if
 333      * delta is positive, and at -b/a if delta is zero.
 334      */
 335
 336 #define ADD(t0) \
 337     { \
 338         double _t0 = (t0); \
 339         if (0 < _t0 && _t0 < 1) \
 340             t[t_num++] = _t0; \
 341     }
 342
 343 #define FIND_EXTREMES(a,b,c) \
 344     { \
 345         if (a == 0) { \
 346             if (b != 0) \
 347                 ADD (-c / (2*b)); \
 348         } else { \
 349             double b2 = b * b; \
 350             double delta = b2 - a * c; \
 351             if (delta > 0) { \
 352                 cairo_bool_t feasible; \
 353                 double _2ab = 2 * a * b; \
 354                 /* We are only interested in solutions t that satisfy 0<t<1 \
 355                  * here.  We do some checks to avoid sqrt if the solutions \
 356                  * are not in that range.  The checks can be derived from: \
 357                  * \
 358                  *   0 < (-b±√delta)/a < 1 \
 359                  */ \
 360                 if (_2ab >= 0) \
 361                     feasible = delta > b2 && delta < a*a + b2 + _2ab; \
 362                 else if (-b / a >= 1) \
 363                     feasible = delta < b2 && delta > a*a + b2 + _2ab; \
 364                 else \
 365                     feasible = delta < b2 || delta < a*a + b2 + _2ab; \
 366                 \
 367                 if (unlikely (feasible)) { \
 368                     double sqrt_delta = sqrt (delta); \
 369                     ADD ((-b - sqrt_delta) / a); \
 370                     ADD ((-b + sqrt_delta) / a); \
 371                 } \
 372             } else if (delta == 0) { \
 373                 ADD (-b / a); \
 374             } \
 375         } \
 376     }
 377
 378     /* Find X extremes */
 379     a = -x0 + 3*x1 - 3*x2 + x3;
 380     b =  x0 - 2*x1 + x2;
 381     c = -x0 + x1;
 382     FIND_EXTREMES (a, b, c);
 383
 384     /* Find Y extremes */
 385     a = -y0 + 3*y1 - 3*y2 + y3;
 386     b =  y0 - 2*y1 + y2;
 387     c = -y0 + y1;
 388     FIND_EXTREMES (a, b, c);
 389
 390     status = add_point_func (closure, p0, NULL);
 391     if (unlikely (status))
 392         return status;
 393
 394     for (i = 0; i < t_num; i++) {
 395         cairo_point_t p;
 396         double x, y;
 397         double t_1_0, t_0_1;
 398         double t_2_0, t_0_2;
 399         double t_3_0, t_2_1_3, t_1_2_3, t_0_3;
 400
 401         t_1_0 = t[i];          /*      t  */
 402         t_0_1 = 1 - t_1_0;     /* (1 - t) */
 403
 404         t_2_0 = t_1_0 * t_1_0; /*      t  *      t  */
 405         t_0_2 = t_0_1 * t_0_1; /* (1 - t) * (1 - t) */
 406
 407         t_3_0   = t_2_0 * t_1_0;     /*      t  *      t  *      t      */
 408         t_2_1_3 = t_2_0 * t_0_1 * 3; /*      t  *      t  * (1 - t) * 3 */
 409         t_1_2_3 = t_1_0 * t_0_2 * 3; /*      t  * (1 - t) * (1 - t) * 3 */
 410         t_0_3   = t_0_1 * t_0_2;     /* (1 - t) * (1 - t) * (1 - t)     */
 411
 412         /* Bezier polynomial */
 413         x = x0 * t_0_3
 414           + x1 * t_1_2_3
 415           + x2 * t_2_1_3
 416           + x3 * t_3_0;
 417         y = y0 * t_0_3
 418           + y1 * t_1_2_3
 419           + y2 * t_2_1_3
 420           + y3 * t_3_0;
 421
 422         p.x = _cairo_fixed_from_double (x);
 423         p.y = _cairo_fixed_from_double (y);
 424         status = add_point_func (closure, &p, NULL);
 425         if (unlikely (status))
 426             return status;
 427     }
 428
 429     return add_point_func (closure, p3, NULL);
 430 }