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