src/gallium/state_trackers/vega/arc.c

   1 /**************************************************************************
   2  *
   3  * Copyright 2009 VMware, Inc.  All Rights Reserved.
   4  *
   5  * Permission is hereby granted, free of charge, to any person obtaining a
   6  * copy of this software and associated documentation files (the
   7  * "Software"), to deal in the Software without restriction, including
   8  * without limitation the rights to use, copy, modify, merge, publish,
   9  * distribute, sub license, and/or sell copies of the Software, and to
  10  * permit persons to whom the Software is furnished to do so, subject to
  11  * the following conditions:
  12  *
  13  * The above copyright notice and this permission notice (including the
  14  * next paragraph) shall be included in all copies or substantial portions
  15  * of the Software.
  16  *
  17  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
  18  * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
  19  * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NON-INFRINGEMENT.
  20  * IN NO EVENT SHALL VMWARE AND/OR ITS SUPPLIERS BE LIABLE FOR
  21  * ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
  22  * TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
  23  * SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
  24  *
  25  **************************************************************************/
  26
  27 #include "arc.h"
  28
  29 #include "matrix.h"
  30 #include "bezier.h"
  31 #include "polygon.h"
  32 #include "stroker.h"
  33 #include "path.h"
  34
  35 #include "util/u_debug.h"
  36
  37 #include <math.h>
  38
  39 #ifndef M_PI
  40 #define M_PI 3.14159265358979323846
  41 #endif
  42
  43 #define DEBUG_ARCS 0
  44
  45 static const VGfloat two_pi = M_PI * 2;
  46
  47
  48 static const double coeffs3Low[2][4][4] = {
  49    {
  50       {  3.85268,   -21.229,      -0.330434,    0.0127842  },
  51       { -1.61486,     0.706564,    0.225945,    0.263682   },
  52       { -0.910164,    0.388383,    0.00551445,  0.00671814 },
  53       { -0.630184,    0.192402,    0.0098871,   0.0102527  }
  54    },
  55    {
  56       { -0.162211,    9.94329,     0.13723,     0.0124084  },
  57       { -0.253135,    0.00187735,  0.0230286,   0.01264    },
  58       { -0.0695069,  -0.0437594,   0.0120636,   0.0163087  },
  59       { -0.0328856,  -0.00926032, -0.00173573,  0.00527385 }
  60    }
  61 };
  62
  63 /* coefficients for error estimation
  64    while using cubic Bézier curves for approximation
  65    1/4 <= b/a <= 1 */
  66 static const double coeffs3High[2][4][4] = {
  67    {
  68       {  0.0899116, -19.2349,     -4.11711,     0.183362   },
  69       {  0.138148,   -1.45804,     1.32044,     1.38474    },
  70       {  0.230903,   -0.450262,    0.219963,    0.414038   },
  71       {  0.0590565,  -0.101062,    0.0430592,   0.0204699  }
  72    },
  73    {
  74       {  0.0164649,   9.89394,     0.0919496,   0.00760802 },
  75       {  0.0191603,  -0.0322058,   0.0134667,  -0.0825018  },
  76       {  0.0156192,  -0.017535,    0.00326508, -0.228157   },
  77       { -0.0236752,   0.0405821,  -0.0173086,   0.176187   }
  78    }
  79 };
  80
  81 /* safety factor to convert the "best" error approximation
  82    into a "max bound" error */
  83 static const double safety3[] = {
  84    0.001, 4.98, 0.207, 0.0067
  85 };
  86
  87 /* The code below is from the OpenVG 1.1 Spec
  88  * Section 18.4 */
  89
  90 /* Given: Points (x0, y0) and (x1, y1)
  91  * Return: TRUE if a solution exists, FALSE otherwise
  92  *         Circle centers are written to (cx0, cy0) and (cx1, cy1)
  93  */
  94 static VGboolean
  95 find_unit_circles(double x0, double y0, double x1, double y1,
  96                   double *cx0, double *cy0,
  97                   double *cx1, double *cy1)
  98 {
  99    /* Compute differences and averages */
 100    double dx = x0 - x1;
 101    double dy = y0 - y1;
 102    double xm = (x0 + x1)/2;
 103    double ym = (y0 + y1)/2;
 104    double dsq, disc, s, sdx, sdy;
 105
 106    /* Solve for intersecting unit circles */
 107    dsq = dx*dx + dy*dy;
 108    if (dsq == 0.0) return VG_FALSE; /* Points are coincident */
 109    disc = 1.0/dsq - 1.0/4.0;
 110
 111    /* the precision we care about here is around float so if we're
 112     * around the float defined zero then make it official to avoid
 113     * precision problems later on */
 114    if (floatIsZero(disc))
 115       disc = 0.0;
 116
 117    if (disc < 0.0) return VG_FALSE; /* Points are too far apart */
 118    s = sqrt(disc);
 119    sdx = s*dx;
 120    sdy = s*dy;
 121    *cx0 = xm + sdy;
 122    *cy0 = ym - sdx;
 123    *cx1 = xm - sdy;
 124    *cy1 = ym + sdx;
 125    return VG_TRUE;
 126 }
 127
 128
 129 /* Given:  Ellipse parameters rh, rv, rot (in degrees),
 130  *         endpoints (x0, y0) and (x1, y1)
 131  * Return: TRUE if a solution exists, FALSE otherwise
 132  *         Ellipse centers are written to (cx0, cy0) and (cx1, cy1)
 133  */
 134 static VGboolean
 135 find_ellipses(double rh, double rv, double rot,
 136               double x0, double y0, double x1, double y1,
 137               double *cx0, double *cy0, double *cx1, double *cy1)
 138 {
 139    double COS, SIN, x0p, y0p, x1p, y1p, pcx0, pcy0, pcx1, pcy1;
 140    /* Convert rotation angle from degrees to radians */
 141    rot *= M_PI/180.0;
 142    /* Pre-compute rotation matrix entries */
 143    COS = cos(rot); SIN = sin(rot);
 144    /* Transform (x0, y0) and (x1, y1) into unit space */
 145    /* using (inverse) rotate, followed by (inverse) scale   */
 146    x0p = (x0*COS + y0*SIN)/rh;
 147    y0p = (-x0*SIN + y0*COS)/rv;
 148    x1p = (x1*COS + y1*SIN)/rh;
 149    y1p = (-x1*SIN + y1*COS)/rv;
 150    if (!find_unit_circles(x0p, y0p, x1p, y1p,
 151                           &pcx0, &pcy0, &pcx1, &pcy1)) {
 152       return VG_FALSE;
 153    }
 154    /* Transform back to original coordinate space */
 155    /* using (forward) scale followed by (forward) rotate */
 156    pcx0 *= rh; pcy0 *= rv;
 157    pcx1 *= rh; pcy1 *= rv;
 158    *cx0 = pcx0*COS - pcy0*SIN;
 159    *cy0 = pcx0*SIN + pcy0*COS;
 160    *cx1 = pcx1*COS - pcy1*SIN;
 161    *cy1 = pcx1*SIN + pcy1*COS;
 162    return VG_TRUE;
 163 }
 164
 165 static INLINE VGboolean
 166 try_to_fix_radii(struct arc *arc)
 167 {
 168    double COS, SIN, rot, x0p, y0p, x1p, y1p;
 169    double dx, dy, dsq, scale;
 170
 171    /* Convert rotation angle from degrees to radians */
 172    rot = DEGREES_TO_RADIANS(arc->theta);
 173
 174    /* Pre-compute rotation matrix entries */
 175    COS = cos(rot); SIN = sin(rot);
 176
 177    /* Transform (x0, y0) and (x1, y1) into unit space */
 178    /* using (inverse) rotate, followed by (inverse) scale   */
 179    x0p = (arc->x1*COS + arc->y1*SIN)/arc->a;
 180    y0p = (-arc->x1*SIN + arc->y1*COS)/arc->b;
 181    x1p = (arc->x2*COS + arc->y2*SIN)/arc->a;
 182    y1p = (-arc->x2*SIN + arc->y2*COS)/arc->b;
 183    /* Compute differences and averages */
 184    dx = x0p - x1p;
 185    dy = y0p - y1p;
 186
 187    dsq = dx*dx + dy*dy;
 188 #if 0
 189    if (dsq <= 0.001) {
 190       debug_printf("AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAaaaaa\n");
 191    }
 192 #endif
 193    scale = 1/(2/sqrt(dsq));
 194    arc->a *= scale;
 195    arc->b *= scale;
 196    return VG_TRUE;
 197 }
 198
 199 static INLINE double vector_normalize(double *v)
 200 {
 201    double sq = v[0] * v[0] + v[1] * v[1];
 202    return sqrt(sq);
 203 }
 204 static INLINE double vector_orientation(double *v)
 205 {
 206    double norm = vector_normalize(v);
 207    double cosa = v[0] / norm;
 208    double sina = v[1] / norm;
 209    return (sina>=0 ? acos(cosa) : 2*M_PI - acos(cosa));
 210 }
 211 static INLINE double vector_dot(double *v0,
 212                                 double *v1)
 213 {
 214    return v0[0] * v1[0] + v0[1] * v1[1];
 215 }
 216
 217 static INLINE double vector_angles(double *v0,
 218                                    double *v1)
 219 {
 220    double dot = vector_dot(v0, v1);
 221    double norm0 = vector_normalize(v0);
 222    double norm1 = vector_normalize(v1);
 223
 224    return acos(dot / (norm0 * norm1));
 225 }
 226
 227 static VGboolean find_angles(struct arc *arc)
 228 {
 229    double vec0[2], vec1[2];
 230    double lambda1, lambda2;
 231    double angle;
 232    struct matrix matrix;
 233
 234    if (floatIsZero(arc->a) || floatIsZero(arc->b)) {
 235       return VG_FALSE;
 236    }
 237    /* map the points to an identity circle */
 238    matrix_load_identity(&matrix);
 239    matrix_scale(&matrix, 1.f, arc->a/arc->b);
 240    matrix_rotate(&matrix, -arc->theta);
 241    matrix_map_point(&matrix,
 242                     arc->x1, arc->y1,
 243                     &arc->x1, &arc->y1);
 244    matrix_map_point(&matrix,
 245                     arc->x2, arc->y2,
 246                     &arc->x2, &arc->y2);
 247    matrix_map_point(&matrix,
 248                     arc->cx, arc->cy,
 249                     &arc->cx, &arc->cy);
 250
 251 #if DEBUG_ARCS
 252    debug_printf("Matrix 3 [%f, %f, %f| %f, %f, %f| %f, %f, %f]\n",
 253                 matrix.m[0], matrix.m[1], matrix.m[2],
 254                 matrix.m[3], matrix.m[4], matrix.m[5],
 255                 matrix.m[6], matrix.m[7], matrix.m[8]);
 256    debug_printf("Endpoints [%f, %f], [%f, %f]\n",
 257                 arc->x1, arc->y1, arc->x2, arc->y2);
 258 #endif
 259
 260    vec0[0] = arc->x1 - arc->cx;
 261    vec0[1] = arc->y1 - arc->cy;
 262    vec1[0] = arc->x2 - arc->cx;
 263    vec1[1] = arc->y2 - arc->cy;
 264
 265 #if DEBUG_ARCS
 266    debug_printf("Vec is [%f, %f], [%f, %f], [%f, %f]\n",
 267                 vec0[0], vec0[1], vec1[0], vec1[1], arc->cx, arc->cy);
 268 #endif
 269
 270    lambda1 = vector_orientation(vec0);
 271
 272    if (isnan(lambda1))
 273       lambda1 = 0.f;
 274
 275    if (arc->type == VG_SCWARC_TO ||
 276        arc->type == VG_SCCWARC_TO)
 277       angle = vector_angles(vec0, vec1);
 278    else if (arc->type == VG_LCWARC_TO ||
 279             arc->type == VG_LCCWARC_TO) {
 280       angle = 2*M_PI - vector_angles(vec0, vec1);
 281    } else
 282       abort();
 283
 284    if (isnan(angle))
 285       angle = M_PI;
 286
 287
 288    if (arc->type == VG_SCWARC_TO ||
 289        arc->type == VG_LCWARC_TO)
 290       lambda2 = lambda1 - angle;
 291    else
 292       lambda2 = lambda1 + angle;
 293
 294 #if DEBUG_ARCS
 295    debug_printf("Angle is %f and (%f, %f)\n", angle, lambda1, lambda2);
 296 #endif
 297
 298 #if 0
 299    arc->eta1 = atan2(sin(lambda1) / arc->b,
 300                      cos(lambda1) / arc->a);
 301    arc->eta2 = atan2(sin(lambda2) / arc->b,
 302                      cos(lambda2) / arc->a);
 303
 304    /* make sure we have eta1 <= eta2 <= eta1 + 2 PI */
 305    arc->eta2 -= two_pi * floor((arc->eta2 - arc->eta1) / two_pi);
 306
 307    /* the preceding correction fails if we have exactly et2 - eta1 = 2 PI
 308       it reduces the interval to zero length */
 309    if ((lambda2 - lambda1 > M_PI) && (arc->eta2 - arc->eta1 < M_PI)) {
 310       arc->eta2 += 2 * M_PI;
 311    }
 312 #else
 313    arc->eta1 = lambda1;
 314    arc->eta2 = lambda2;
 315 #endif
 316
 317    return VG_TRUE;
 318 }
 319
 320 #if DEBUG_ARCS
 321 static void check_endpoints(struct arc *arc)
 322 {
 323    double x1, y1, x2, y2;
 324
 325    double a_cos_eta1 = arc->a * cos(arc->eta1);
 326    double b_sin_eta1 = arc->b * sin(arc->eta1);
 327    x1 = arc->cx + a_cos_eta1 * arc->cos_theta -
 328         b_sin_eta1 * arc->sin_theta;
 329    y1 = arc->cy + a_cos_eta1 * arc->sin_theta +
 330         b_sin_eta1 * arc->cos_theta;
 331
 332    double a_cos_eta2 = arc->a * cos(arc->eta2);
 333    double b_sin_eta2 = arc->b * sin(arc->eta2);
 334    x2 = arc->cx + a_cos_eta2 * arc->cos_theta -
 335         b_sin_eta2 * arc->sin_theta;
 336    y2 = arc->cy + a_cos_eta2 * arc->sin_theta +
 337         b_sin_eta2 * arc->cos_theta;
 338
 339    debug_printf("Computed (%f, %f), (%f, %f)\n",
 340                 x1, y1, x2, y2);
 341    debug_printf("Real     (%f, %f), (%f, %f)\n",
 342                 arc->x1, arc->y1,
 343                 arc->x2, arc->y2);
 344 }
 345 #endif
 346
 347 void arc_init(struct arc *arc,
 348               VGPathSegment type,
 349               VGfloat x1, VGfloat y1,
 350               VGfloat x2, VGfloat y2,
 351               VGfloat rh, VGfloat rv,
 352               VGfloat rot)
 353 {
 354    assert(type == VG_SCCWARC_TO ||
 355           type == VG_SCWARC_TO ||
 356           type == VG_LCCWARC_TO ||
 357           type == VG_LCWARC_TO);
 358    arc->type = type;
 359    arc->x1  = x1;
 360    arc->y1  = y1;
 361    arc->x2  = x2;
 362    arc->y2  = y2;
 363    arc->a   = rh;
 364    arc->b   = rv;
 365    arc->theta = rot;
 366    arc->cos_theta = cos(arc->theta);
 367    arc->sin_theta = sin(arc->theta);
 368    {
 369       double cx0, cy0, cx1, cy1;
 370       double cx, cy;
 371       arc->is_valid =  find_ellipses(rh, rv, rot, x1, y1, x2, y2,
 372                                      &cx0, &cy0, &cx1, &cy1);
 373
 374       if (!arc->is_valid && try_to_fix_radii(arc)) {
 375          rh = arc->a;
 376          rv = arc->b;
 377          arc->is_valid =
 378             find_ellipses(rh, rv, rot, x1, y1, x2, y2,
 379                           &cx0, &cy0, &cx1, &cy1);
 380       }
 381
 382       if (type == VG_SCWARC_TO ||
 383           type == VG_LCCWARC_TO) {
 384          cx = cx1;
 385          cy = cy1;
 386       } else {
 387          cx = cx0;
 388          cy = cy0;
 389       }
 390 #if DEBUG_ARCS
 391       debug_printf("Centers are : (%f, %f) , (%f, %f). Real (%f, %f)\n",
 392                    cx0, cy0, cx1, cy1, cx, cy);
 393 #endif
 394       arc->cx = cx;
 395       arc->cy = cy;
 396       if (arc->is_valid) {
 397          arc->is_valid = find_angles(arc);
 398 #if DEBUG_ARCS
 399          check_endpoints(arc);
 400 #endif
 401          /* remap a few points. find_angles requires
 402           * rot in angles, the rest of the code
 403           * will need them in radians. and find_angles
 404           * modifies the center to match an identity
 405           * circle so lets reset it */
 406          arc->theta = DEGREES_TO_RADIANS(rot);
 407          arc->cos_theta = cos(arc->theta);
 408          arc->sin_theta = sin(arc->theta);
 409          arc->cx = cx;
 410          arc->cy = cy;
 411       }
 412    }
 413 }
 414
 415 static INLINE double rational_function(double x, const double *c)
 416 {
 417    return (x * (x * c[0] + c[1]) + c[2]) / (x + c[3]);
 418 }
 419
 420 static double estimate_error(struct arc *arc,
 421                              double etaA, double etaB)
 422 {
 423    double eta  = 0.5 * (etaA + etaB);
 424
 425    double x    = arc->b / arc->a;
 426    double dEta = etaB - etaA;
 427    double cos2 = cos(2 * eta);
 428    double cos4 = cos(4 * eta);
 429    double cos6 = cos(6 * eta);
 430    double c0, c1;
 431
 432    /* select the right coeficients set according to degree and b/a */
 433    const double (*coeffs)[4][4];
 434    const double *safety;
 435    coeffs = (x < 0.25) ? coeffs3Low : coeffs3High;
 436    safety = safety3;
 437
 438    c0 = rational_function(x, coeffs[0][0])
 439         + cos2 * rational_function(x, coeffs[0][1])
 440         + cos4 * rational_function(x, coeffs[0][2])
 441         + cos6 * rational_function(x, coeffs[0][3]);
 442
 443    c1 = rational_function(x, coeffs[1][0])
 444         + cos2 * rational_function(x, coeffs[1][1])
 445         + cos4 * rational_function(x, coeffs[1][2])
 446         + cos6 * rational_function(x, coeffs[1][3]);
 447
 448    return rational_function(x, safety) * arc->a * exp(c0 + c1 * dEta);
 449 }
 450
 451 struct arc_cb {
 452    void (*move)(struct arc_cb *cb, VGfloat x, VGfloat y);
 453    void (*point)(struct arc_cb *cb, VGfloat x, VGfloat y);
 454    void (*bezier)(struct arc_cb *cb, struct bezier *bezier);
 455
 456    void *user_data;
 457 };
 458
 459 static void cb_null_move(struct arc_cb *cb, VGfloat x, VGfloat y)
 460 {
 461 }
 462
 463 static void polygon_point(struct arc_cb *cb, VGfloat x, VGfloat y)
 464 {
 465    struct polygon *poly = (struct polygon*)cb->user_data;
 466    polygon_vertex_append(poly, x, y);
 467 }
 468
 469 static void polygon_bezier(struct arc_cb *cb, struct bezier *bezier)
 470 {
 471    struct polygon *poly = (struct polygon*)cb->user_data;
 472    bezier_add_to_polygon(bezier, poly);
 473 }
 474
 475 static void stroke_point(struct arc_cb *cb, VGfloat x, VGfloat y)
 476 {
 477    struct stroker *stroker = (struct stroker*)cb->user_data;
 478    stroker_line_to(stroker, x, y);
 479 }
 480
 481 static void stroke_curve(struct arc_cb *cb, struct bezier *bezier)
 482 {
 483    struct stroker *stroker = (struct stroker*)cb->user_data;
 484    stroker_curve_to(stroker,
 485                     bezier->x2, bezier->y2,
 486                     bezier->x3, bezier->y3,
 487                     bezier->x4, bezier->y4);
 488 }
 489
 490 static void stroke_emit_point(struct arc_cb *cb, VGfloat x, VGfloat y)
 491 {
 492    struct stroker *stroker = (struct stroker*)cb->user_data;
 493    stroker_emit_line_to(stroker, x, y);
 494 }
 495
 496 static void stroke_emit_curve(struct arc_cb *cb, struct bezier *bezier)
 497 {
 498    struct stroker *stroker = (struct stroker*)cb->user_data;
 499    stroker_emit_curve_to(stroker,
 500                          bezier->x2, bezier->y2,
 501                          bezier->x3, bezier->y3,
 502                          bezier->x4, bezier->y4);
 503 }
 504
 505 static void arc_path_move(struct arc_cb *cb, VGfloat x, VGfloat y)
 506 {
 507    struct path *path = (struct path*)cb->user_data;
 508    path_move_to(path, x, y);
 509 }
 510
 511 static void arc_path_point(struct arc_cb *cb, VGfloat x, VGfloat y)
 512 {
 513    struct path *path = (struct path*)cb->user_data;
 514    path_line_to(path, x, y);
 515 }
 516
 517 static void arc_path_bezier(struct arc_cb *cb, struct bezier *bezier)
 518 {
 519    struct path *path = (struct path*)cb->user_data;
 520    path_cubic_to(path,
 521                  bezier->x2, bezier->y2,
 522                  bezier->x3, bezier->y3,
 523                  bezier->x4, bezier->y4);
 524 }
 525
 526 static INLINE int num_beziers_needed(struct arc *arc)
 527 {
 528    double threshold = 0.05;
 529    VGboolean found = VG_FALSE;
 530    int n = 1;
 531    double min_eta, max_eta;
 532
 533    min_eta = MIN2(arc->eta1, arc->eta2);
 534    max_eta = MAX2(arc->eta1, arc->eta2);
 535
 536    while ((! found) && (n < 1024)) {
 537       double d_eta = (max_eta - min_eta) / n;
 538       if (d_eta <= 0.5 * M_PI) {
 539          double eta_b = min_eta;
 540          int i;
 541          found = VG_TRUE;
 542          for (i = 0; found && (i < n); ++i) {
 543             double etaA = eta_b;
 544             eta_b += d_eta;
 545             found = (estimate_error(arc, etaA, eta_b) <= threshold);
 546          }
 547       }
 548       n = n << 1;
 549    }
 550
 551    return n;
 552 }
 553
 554 static void arc_to_beziers(struct arc *arc,
 555                            struct arc_cb cb,
 556                            struct matrix *matrix)
 557 {
 558    int i;
 559    int n = 1;
 560    double d_eta, eta_b, cos_eta_b,
 561       sin_eta_b, a_cos_eta_b, b_sin_eta_b, a_sin_eta_b,
 562       b_cos_eta_b, x_b, y_b, x_b_dot, y_b_dot, lx, ly;
 563    double t, alpha;
 564
 565    { /* always move to the start of the arc */
 566       VGfloat x = arc->x1;
 567       VGfloat y = arc->y1;
 568       matrix_map_point(matrix, x, y, &x, &y);
 569       cb.move(&cb, x, y);
 570    }
 571
 572    if (!arc->is_valid) {
 573       VGfloat x = arc->x2;
 574       VGfloat y = arc->y2;
 575       matrix_map_point(matrix, x, y, &x, &y);
 576       cb.point(&cb, x, y);
 577       return;
 578    }
 579
 580    /* find the number of Bézier curves needed */
 581    n = num_beziers_needed(arc);
 582
 583    d_eta = (arc->eta2 - arc->eta1) / n;
 584    eta_b = arc->eta1;
 585
 586    cos_eta_b  = cos(eta_b);
 587    sin_eta_b  = sin(eta_b);
 588    a_cos_eta_b = arc->a * cos_eta_b;
 589    b_sin_eta_b = arc->b * sin_eta_b;
 590    a_sin_eta_b = arc->a * sin_eta_b;
 591    b_cos_eta_b = arc->b * cos_eta_b;
 592    x_b       = arc->cx + a_cos_eta_b * arc->cos_theta -
 593                b_sin_eta_b * arc->sin_theta;
 594    y_b       = arc->cy + a_cos_eta_b * arc->sin_theta +
 595                b_sin_eta_b * arc->cos_theta;
 596    x_b_dot    = -a_sin_eta_b * arc->cos_theta -
 597                 b_cos_eta_b * arc->sin_theta;
 598    y_b_dot    = -a_sin_eta_b * arc->sin_theta +
 599                 b_cos_eta_b * arc->cos_theta;
 600
 601    {
 602       VGfloat x = x_b, y = y_b;
 603       matrix_map_point(matrix, x, y, &x, &y);
 604       cb.point(&cb, x, y);
 605    }
 606    lx = x_b;
 607    ly = y_b;
 608
 609    t     = tan(0.5 * d_eta);
 610    alpha = sin(d_eta) * (sqrt(4 + 3 * t * t) - 1) / 3;
 611
 612    for (i = 0; i < n; ++i) {
 613       struct bezier bezier;
 614       double xA    = x_b;
 615       double yA    = y_b;
 616       double xADot = x_b_dot;
 617       double yADot = y_b_dot;
 618
 619       eta_b    += d_eta;
 620       cos_eta_b  = cos(eta_b);
 621       sin_eta_b  = sin(eta_b);
 622       a_cos_eta_b = arc->a * cos_eta_b;
 623       b_sin_eta_b = arc->b * sin_eta_b;
 624       a_sin_eta_b = arc->a * sin_eta_b;
 625       b_cos_eta_b = arc->b * cos_eta_b;
 626       x_b       = arc->cx + a_cos_eta_b * arc->cos_theta -
 627                   b_sin_eta_b * arc->sin_theta;
 628       y_b       = arc->cy + a_cos_eta_b * arc->sin_theta +
 629                   b_sin_eta_b * arc->cos_theta;
 630       x_b_dot    = -a_sin_eta_b * arc->cos_theta -
 631                    b_cos_eta_b * arc->sin_theta;
 632       y_b_dot    = -a_sin_eta_b * arc->sin_theta +
 633                    b_cos_eta_b * arc->cos_theta;
 634
 635       bezier_init(&bezier,
 636                   lx, ly,
 637                   (float) (xA + alpha * xADot), (float) (yA + alpha * yADot),
 638                   (float) (x_b - alpha * x_b_dot), (float) (y_b - alpha * y_b_dot),
 639                   (float) x_b,                   (float) y_b);
 640 #if 0
 641       debug_printf("%d) Bezier (%f, %f), (%f, %f), (%f, %f), (%f, %f)\n",
 642                    i,
 643                    bezier.x1, bezier.y1,
 644                    bezier.x2, bezier.y2,
 645                    bezier.x3, bezier.y3,
 646                    bezier.x4, bezier.y4);
 647 #endif
 648       bezier_transform(&bezier, matrix);
 649       cb.bezier(&cb, &bezier);
 650       lx = x_b;
 651       ly = y_b;
 652    }
 653 }
 654
 655
 656 void arc_add_to_polygon(struct arc *arc,
 657                         struct polygon *poly,
 658                         struct matrix *matrix)
 659 {
 660    struct arc_cb cb;
 661
 662    cb.move = cb_null_move;
 663    cb.point = polygon_point;
 664    cb.bezier = polygon_bezier;
 665    cb.user_data = poly;
 666
 667    arc_to_beziers(arc, cb, matrix);
 668 }
 669
 670 void arc_stroke_cb(struct arc *arc,
 671                    struct stroker *stroke,
 672                    struct matrix *matrix)
 673 {
 674    struct arc_cb cb;
 675
 676    cb.move = cb_null_move;
 677    cb.point = stroke_point;
 678    cb.bezier = stroke_curve;
 679    cb.user_data = stroke;
 680
 681    arc_to_beziers(arc, cb, matrix);
 682 }
 683
 684 void arc_stroker_emit(struct arc *arc,
 685                       struct stroker *stroker,
 686                       struct matrix *matrix)
 687 {
 688    struct arc_cb cb;
 689
 690    cb.move = cb_null_move;
 691    cb.point = stroke_emit_point;
 692    cb.bezier = stroke_emit_curve;
 693    cb.user_data = stroker;
 694
 695    arc_to_beziers(arc, cb, matrix);
 696 }
 697
 698 void arc_to_path(struct arc *arc,
 699                  struct path *path,
 700                  struct matrix *matrix)
 701 {
 702    struct arc_cb cb;
 703
 704    cb.move = arc_path_move;
 705    cb.point = arc_path_point;
 706    cb.bezier = arc_path_bezier;
 707    cb.user_data = path;
 708
 709    arc_to_beziers(arc, cb, matrix);
 710 }