src/cairo-bentley-ottmann.c

   1 /*
   2  * Copyright © 2004 Carl Worth
   3  * Copyright © 2006 Red Hat, Inc.
   4  * Copyright © 2008 Chris Wilson
   5  *
   6  * This library is free software; you can redistribute it and/or
   7  * modify it either under the terms of the GNU Lesser General Public
   8  * License version 2.1 as published by the Free Software Foundation
   9  * (the "LGPL") or, at your option, under the terms of the Mozilla
  10  * Public License Version 1.1 (the "MPL"). If you do not alter this
  11  * notice, a recipient may use your version of this file under either
  12  * the MPL or the LGPL.
  13  *
  14  * You should have received a copy of the LGPL along with this library
  15  * in the file COPYING-LGPL-2.1; if not, write to the Free Software
  16  * Foundation, Inc., 51 Franklin Street, Suite 500, Boston, MA 02110-1335, USA
  17  * You should have received a copy of the MPL along with this library
  18  * in the file COPYING-MPL-1.1
  19  *
  20  * The contents of this file are subject to the Mozilla Public License
  21  * Version 1.1 (the "License"); you may not use this file except in
  22  * compliance with the License. You may obtain a copy of the License at
  23  * http://www.mozilla.org/MPL/
  24  *
  25  * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
  26  * OF ANY KIND, either express or implied. See the LGPL or the MPL for
  27  * the specific language governing rights and limitations.
  28  *
  29  * The Original Code is the cairo graphics library.
  30  *
  31  * The Initial Developer of the Original Code is Carl Worth
  32  *
  33  * Contributor(s):
  34  *      Carl D. Worth <cworth@cworth.org>
  35  *      Chris Wilson <chris@chris-wilson.co.uk>
  36  */
  37
  38 /* Provide definitions for standalone compilation */
  39 #include "cairoint.h"
  40
  41 #include "cairo-error-private.h"
  42 #include "cairo-freelist-private.h"
  43 #include "cairo-combsort-inline.h"
  44 #include "cairo-traps-private.h"
  45
  46 #define DEBUG_PRINT_STATE 0
  47 #define DEBUG_EVENTS 0
  48 #define DEBUG_TRAPS 0
  49
  50 typedef cairo_point_t cairo_bo_point32_t;
  51
  52 typedef struct _cairo_bo_intersect_ordinate {
  53     int32_t ordinate;
  54     enum { EXACT, INEXACT } exactness;
  55 } cairo_bo_intersect_ordinate_t;
  56
  57 typedef struct _cairo_bo_intersect_point {
  58     cairo_bo_intersect_ordinate_t x;
  59     cairo_bo_intersect_ordinate_t y;
  60 } cairo_bo_intersect_point_t;
  61
  62 typedef struct _cairo_bo_edge cairo_bo_edge_t;
  63 typedef struct _cairo_bo_trap cairo_bo_trap_t;
  64
  65 /* A deferred trapezoid of an edge */
  66 struct _cairo_bo_trap {
  67     cairo_bo_edge_t *right;
  68     int32_t top;
  69 };
  70
  71 struct _cairo_bo_edge {
  72     cairo_edge_t edge;
  73     cairo_bo_edge_t *prev;
  74     cairo_bo_edge_t *next;
  75     cairo_bo_edge_t *colinear;
  76     cairo_bo_trap_t deferred_trap;
  77 };
  78
  79 /* the parent is always given by index/2 */
  80 #define PQ_PARENT_INDEX(i) ((i) >> 1)
  81 #define PQ_FIRST_ENTRY 1
  82
  83 /* left and right children are index * 2 and (index * 2) +1 respectively */
  84 #define PQ_LEFT_CHILD_INDEX(i) ((i) << 1)
  85
  86 typedef enum {
  87     CAIRO_BO_EVENT_TYPE_STOP,
  88     CAIRO_BO_EVENT_TYPE_INTERSECTION,
  89     CAIRO_BO_EVENT_TYPE_START
  90 } cairo_bo_event_type_t;
  91
  92 typedef struct _cairo_bo_event {
  93     cairo_bo_event_type_t type;
  94     cairo_point_t point;
  95 } cairo_bo_event_t;
  96
  97 typedef struct _cairo_bo_start_event {
  98     cairo_bo_event_type_t type;
  99     cairo_point_t point;
 100     cairo_bo_edge_t edge;
 101 } cairo_bo_start_event_t;
 102
 103 typedef struct _cairo_bo_queue_event {
 104     cairo_bo_event_type_t type;
 105     cairo_point_t point;
 106     cairo_bo_edge_t *e1;
 107     cairo_bo_edge_t *e2;
 108 } cairo_bo_queue_event_t;
 109
 110 typedef struct _pqueue {
 111     int size, max_size;
 112
 113     cairo_bo_event_t **elements;
 114     cairo_bo_event_t *elements_embedded[1024];
 115 } pqueue_t;
 116
 117 typedef struct _cairo_bo_event_queue {
 118     cairo_freepool_t pool;
 119     pqueue_t pqueue;
 120     cairo_bo_event_t **start_events;
 121 } cairo_bo_event_queue_t;
 122
 123 typedef struct _cairo_bo_sweep_line {
 124     cairo_bo_edge_t *head;
 125     cairo_bo_edge_t *stopped;
 126     int32_t current_y;
 127     cairo_bo_edge_t *current_edge;
 128 } cairo_bo_sweep_line_t;
 129
 130 #if DEBUG_TRAPS
 131 static void
 132 dump_traps (cairo_traps_t *traps, const char *filename)
 133 {
 134     FILE *file;
 135     cairo_box_t extents;
 136     int n;
 137
 138     if (getenv ("CAIRO_DEBUG_TRAPS") == NULL)
 139         return;
 140
 141 #if 0
 142     if (traps->has_limits) {
 143         printf ("%s: limits=(%d, %d, %d, %d)\n",
 144                 filename,
 145                 traps->limits.p1.x, traps->limits.p1.y,
 146                 traps->limits.p2.x, traps->limits.p2.y);
 147     }
 148 #endif
 149     _cairo_traps_extents (traps, &extents);
 150     printf ("%s: extents=(%d, %d, %d, %d)\n",
 151             filename,
 152             extents.p1.x, extents.p1.y,
 153             extents.p2.x, extents.p2.y);
 154
 155     file = fopen (filename, "a");
 156     if (file != NULL) {
 157         for (n = 0; n < traps->num_traps; n++) {
 158             fprintf (file, "%d %d L:(%d, %d), (%d, %d) R:(%d, %d), (%d, %d)\n",
 159                      traps->traps[n].top,
 160                      traps->traps[n].bottom,
 161                      traps->traps[n].left.p1.x,
 162                      traps->traps[n].left.p1.y,
 163                      traps->traps[n].left.p2.x,
 164                      traps->traps[n].left.p2.y,
 165                      traps->traps[n].right.p1.x,
 166                      traps->traps[n].right.p1.y,
 167                      traps->traps[n].right.p2.x,
 168                      traps->traps[n].right.p2.y);
 169         }
 170         fprintf (file, "\n");
 171         fclose (file);
 172     }
 173 }
 174
 175 static void
 176 dump_edges (cairo_bo_start_event_t *events,
 177             int num_edges,
 178             const char *filename)
 179 {
 180     FILE *file;
 181     int n;
 182
 183     if (getenv ("CAIRO_DEBUG_TRAPS") == NULL)
 184         return;
 185
 186     file = fopen (filename, "a");
 187     if (file != NULL) {
 188         for (n = 0; n < num_edges; n++) {
 189             fprintf (file, "(%d, %d), (%d, %d) %d %d %d\n",
 190                      events[n].edge.edge.line.p1.x,
 191                      events[n].edge.edge.line.p1.y,
 192                      events[n].edge.edge.line.p2.x,
 193                      events[n].edge.edge.line.p2.y,
 194                      events[n].edge.edge.top,
 195                      events[n].edge.edge.bottom,
 196                      events[n].edge.edge.dir);
 197         }
 198         fprintf (file, "\n");
 199         fclose (file);
 200     }
 201 }
 202 #endif
 203
 204 static cairo_fixed_t
 205 _line_compute_intersection_x_for_y (const cairo_line_t *line,
 206                                     cairo_fixed_t y)
 207 {
 208     cairo_fixed_t x, dy;
 209
 210     if (y == line->p1.y)
 211         return line->p1.x;
 212     if (y == line->p2.y)
 213         return line->p2.x;
 214
 215     x = line->p1.x;
 216     dy = line->p2.y - line->p1.y;
 217     if (dy != 0) {
 218         x += _cairo_fixed_mul_div_floor (y - line->p1.y,
 219                                          line->p2.x - line->p1.x,
 220                                          dy);
 221     }
 222
 223     return x;
 224 }
 225
 226 static inline int
 227 _cairo_bo_point32_compare (cairo_bo_point32_t const *a,
 228                            cairo_bo_point32_t const *b)
 229 {
 230     int cmp;
 231
 232     cmp = a->y - b->y;
 233     if (cmp)
 234         return cmp;
 235
 236     return a->x - b->x;
 237 }
 238
 239 /* Compare the slope of a to the slope of b, returning 1, 0, -1 if the
 240  * slope a is respectively greater than, equal to, or less than the
 241  * slope of b.
 242  *
 243  * For each edge, consider the direction vector formed from:
 244  *
 245  *      top -> bottom
 246  *
 247  * which is:
 248  *
 249  *      (dx, dy) = (line.p2.x - line.p1.x, line.p2.y - line.p1.y)
 250  *
 251  * We then define the slope of each edge as dx/dy, (which is the
 252  * inverse of the slope typically used in math instruction). We never
 253  * compute a slope directly as the value approaches infinity, but we
 254  * can derive a slope comparison without division as follows, (where
 255  * the ? represents our compare operator).
 256  *
 257  * 1.      slope(a) ? slope(b)
 258  * 2.       adx/ady ? bdx/bdy
 259  * 3.   (adx * bdy) ? (bdx * ady)
 260  *
 261  * Note that from step 2 to step 3 there is no change needed in the
 262  * sign of the result since both ady and bdy are guaranteed to be
 263  * greater than or equal to 0.
 264  *
 265  * When using this slope comparison to sort edges, some care is needed
 266  * when interpreting the results. Since the slope compare operates on
 267  * distance vectors from top to bottom it gives a correct left to
 268  * right sort for edges that have a common top point, (such as two
 269  * edges with start events at the same location). On the other hand,
 270  * the sense of the result will be exactly reversed for two edges that
 271  * have a common stop point.
 272  */
 273 static inline int
 274 _slope_compare (const cairo_bo_edge_t *a,
 275                 const cairo_bo_edge_t *b)
 276 {
 277     /* XXX: We're assuming here that dx and dy will still fit in 32
 278      * bits. That's not true in general as there could be overflow. We
 279      * should prevent that before the tessellation algorithm
 280      * begins.
 281      */
 282     int32_t adx = a->edge.line.p2.x - a->edge.line.p1.x;
 283     int32_t bdx = b->edge.line.p2.x - b->edge.line.p1.x;
 284
 285     /* Since the dy's are all positive by construction we can fast
 286      * path several common cases.
 287      */
 288
 289     /* First check for vertical lines. */
 290     if (adx == 0)
 291         return -bdx;
 292     if (bdx == 0)
 293         return adx;
 294
 295     /* Then where the two edges point in different directions wrt x. */
 296     if ((adx ^ bdx) < 0)
 297         return adx;
 298
 299     /* Finally we actually need to do the general comparison. */
 300     {
 301         int32_t ady = a->edge.line.p2.y - a->edge.line.p1.y;
 302         int32_t bdy = b->edge.line.p2.y - b->edge.line.p1.y;
 303         cairo_int64_t adx_bdy = _cairo_int32x32_64_mul (adx, bdy);
 304         cairo_int64_t bdx_ady = _cairo_int32x32_64_mul (bdx, ady);
 305
 306         return _cairo_int64_cmp (adx_bdy, bdx_ady);
 307     }
 308 }
 309
 310 /*
 311  * We need to compare the x-coordinates of a pair of lines for a particular y,
 312  * without loss of precision.
 313  *
 314  * The x-coordinate along an edge for a given y is:
 315  *   X = A_x + (Y - A_y) * A_dx / A_dy
 316  *
 317  * So the inequality we wish to test is:
 318  *   A_x + (Y - A_y) * A_dx / A_dy ∘ B_x + (Y - B_y) * B_dx / B_dy,
 319  * where ∘ is our inequality operator.
 320  *
 321  * By construction, we know that A_dy and B_dy (and (Y - A_y), (Y - B_y)) are
 322  * all positive, so we can rearrange it thus without causing a sign change:
 323  *   A_dy * B_dy * (A_x - B_x) ∘ (Y - B_y) * B_dx * A_dy
 324  *                                 - (Y - A_y) * A_dx * B_dy
 325  *
 326  * Given the assumption that all the deltas fit within 32 bits, we can compute
 327  * this comparison directly using 128 bit arithmetic. For certain, but common,
 328  * input we can reduce this down to a single 32 bit compare by inspecting the
 329  * deltas.
 330  *
 331  * (And put the burden of the work on developing fast 128 bit ops, which are
 332  * required throughout the tessellator.)
 333  *
 334  * See the similar discussion for _slope_compare().
 335  */
 336 static int
 337 edges_compare_x_for_y_general (const cairo_bo_edge_t *a,
 338                                const cairo_bo_edge_t *b,
 339                                int32_t y)
 340 {
 341     /* XXX: We're assuming here that dx and dy will still fit in 32
 342      * bits. That's not true in general as there could be overflow. We
 343      * should prevent that before the tessellation algorithm
 344      * begins.
 345      */
 346     int32_t dx;
 347     int32_t adx, ady;
 348     int32_t bdx, bdy;
 349     enum {
 350        HAVE_NONE    = 0x0,
 351        HAVE_DX      = 0x1,
 352        HAVE_ADX     = 0x2,
 353        HAVE_DX_ADX  = HAVE_DX | HAVE_ADX,
 354        HAVE_BDX     = 0x4,
 355        HAVE_DX_BDX  = HAVE_DX | HAVE_BDX,
 356        HAVE_ADX_BDX = HAVE_ADX | HAVE_BDX,
 357        HAVE_ALL     = HAVE_DX | HAVE_ADX | HAVE_BDX
 358     } have_dx_adx_bdx = HAVE_ALL;
 359
 360     /* don't bother solving for abscissa if the edges' bounding boxes
 361      * can be used to order them. */
 362     {
 363            int32_t amin, amax;
 364            int32_t bmin, bmax;
 365            if (a->edge.line.p1.x < a->edge.line.p2.x) {
 366                    amin = a->edge.line.p1.x;
 367                    amax = a->edge.line.p2.x;
 368            } else {
 369                    amin = a->edge.line.p2.x;
 370                    amax = a->edge.line.p1.x;
 371            }
 372            if (b->edge.line.p1.x < b->edge.line.p2.x) {
 373                    bmin = b->edge.line.p1.x;
 374                    bmax = b->edge.line.p2.x;
 375            } else {
 376                    bmin = b->edge.line.p2.x;
 377                    bmax = b->edge.line.p1.x;
 378            }
 379            if (amax < bmin) return -1;
 380            if (amin > bmax) return +1;
 381     }
 382
 383     ady = a->edge.line.p2.y - a->edge.line.p1.y;
 384     adx = a->edge.line.p2.x - a->edge.line.p1.x;
 385     if (adx == 0)
 386         have_dx_adx_bdx &= ~HAVE_ADX;
 387
 388     bdy = b->edge.line.p2.y - b->edge.line.p1.y;
 389     bdx = b->edge.line.p2.x - b->edge.line.p1.x;
 390     if (bdx == 0)
 391         have_dx_adx_bdx &= ~HAVE_BDX;
 392
 393     dx = a->edge.line.p1.x - b->edge.line.p1.x;
 394     if (dx == 0)
 395         have_dx_adx_bdx &= ~HAVE_DX;
 396
 397 #define L _cairo_int64x32_128_mul (_cairo_int32x32_64_mul (ady, bdy), dx)
 398 #define A _cairo_int64x32_128_mul (_cairo_int32x32_64_mul (adx, bdy), y - a->edge.line.p1.y)
 399 #define B _cairo_int64x32_128_mul (_cairo_int32x32_64_mul (bdx, ady), y - b->edge.line.p1.y)
 400     switch (have_dx_adx_bdx) {
 401     default:
 402     case HAVE_NONE:
 403         return 0;
 404     case HAVE_DX:
 405         /* A_dy * B_dy * (A_x - B_x) ∘ 0 */
 406         return dx; /* ady * bdy is positive definite */
 407     case HAVE_ADX:
 408         /* 0 ∘  - (Y - A_y) * A_dx * B_dy */
 409         return adx; /* bdy * (y - a->top.y) is positive definite */
 410     case HAVE_BDX:
 411         /* 0 ∘ (Y - B_y) * B_dx * A_dy */
 412         return -bdx; /* ady * (y - b->top.y) is positive definite */
 413     case HAVE_ADX_BDX:
 414         /*  0 ∘ (Y - B_y) * B_dx * A_dy - (Y - A_y) * A_dx * B_dy */
 415         if ((adx ^ bdx) < 0) {
 416             return adx;
 417         } else if (a->edge.line.p1.y == b->edge.line.p1.y) { /* common origin */
 418             cairo_int64_t adx_bdy, bdx_ady;
 419
 420             /* ∴ A_dx * B_dy ∘ B_dx * A_dy */
 421
 422             adx_bdy = _cairo_int32x32_64_mul (adx, bdy);
 423             bdx_ady = _cairo_int32x32_64_mul (bdx, ady);
 424
 425             return _cairo_int64_cmp (adx_bdy, bdx_ady);
 426         } else
 427             return _cairo_int128_cmp (A, B);
 428     case HAVE_DX_ADX:
 429         /* A_dy * (A_x - B_x) ∘ - (Y - A_y) * A_dx */
 430         if ((-adx ^ dx) < 0) {
 431             return dx;
 432         } else {
 433             cairo_int64_t ady_dx, dy_adx;
 434
 435             ady_dx = _cairo_int32x32_64_mul (ady, dx);
 436             dy_adx = _cairo_int32x32_64_mul (a->edge.line.p1.y - y, adx);
 437
 438             return _cairo_int64_cmp (ady_dx, dy_adx);
 439         }
 440     case HAVE_DX_BDX:
 441         /* B_dy * (A_x - B_x) ∘ (Y - B_y) * B_dx */
 442         if ((bdx ^ dx) < 0) {
 443             return dx;
 444         } else {
 445             cairo_int64_t bdy_dx, dy_bdx;
 446
 447             bdy_dx = _cairo_int32x32_64_mul (bdy, dx);
 448             dy_bdx = _cairo_int32x32_64_mul (y - b->edge.line.p1.y, bdx);
 449
 450             return _cairo_int64_cmp (bdy_dx, dy_bdx);
 451         }
 452     case HAVE_ALL:
 453         /* XXX try comparing (a->edge.line.p2.x - b->edge.line.p2.x) et al */
 454         return _cairo_int128_cmp (L, _cairo_int128_sub (B, A));
 455     }
 456 #undef B
 457 #undef A
 458 #undef L
 459 }
 460
 461 /*
 462  * We need to compare the x-coordinate of a line for a particular y wrt to a
 463  * given x, without loss of precision.
 464  *
 465  * The x-coordinate along an edge for a given y is:
 466  *   X = A_x + (Y - A_y) * A_dx / A_dy
 467  *
 468  * So the inequality we wish to test is:
 469  *   A_x + (Y - A_y) * A_dx / A_dy ∘ X
 470  * where ∘ is our inequality operator.
 471  *
 472  * By construction, we know that A_dy (and (Y - A_y)) are
 473  * all positive, so we can rearrange it thus without causing a sign change:
 474  *   (Y - A_y) * A_dx ∘ (X - A_x) * A_dy
 475  *
 476  * Given the assumption that all the deltas fit within 32 bits, we can compute
 477  * this comparison directly using 64 bit arithmetic.
 478  *
 479  * See the similar discussion for _slope_compare() and
 480  * edges_compare_x_for_y_general().
 481  */
 482 static int
 483 edge_compare_for_y_against_x (const cairo_bo_edge_t *a,
 484                               int32_t y,
 485                               int32_t x)
 486 {
 487     int32_t adx, ady;
 488     int32_t dx, dy;
 489     cairo_int64_t L, R;
 490
 491     if (x < a->edge.line.p1.x && x < a->edge.line.p2.x)
 492         return 1;
 493     if (x > a->edge.line.p1.x && x > a->edge.line.p2.x)
 494         return -1;
 495
 496     adx = a->edge.line.p2.x - a->edge.line.p1.x;
 497     dx = x - a->edge.line.p1.x;
 498
 499     if (adx == 0)
 500         return -dx;
 501     if (dx == 0 || (adx ^ dx) < 0)
 502         return adx;
 503
 504     dy = y - a->edge.line.p1.y;
 505     ady = a->edge.line.p2.y - a->edge.line.p1.y;
 506
 507     L = _cairo_int32x32_64_mul (dy, adx);
 508     R = _cairo_int32x32_64_mul (dx, ady);
 509
 510     return _cairo_int64_cmp (L, R);
 511 }
 512
 513 static int
 514 edges_compare_x_for_y (const cairo_bo_edge_t *a,
 515                        const cairo_bo_edge_t *b,
 516                        int32_t y)
 517 {
 518     /* If the sweep-line is currently on an end-point of a line,
 519      * then we know its precise x value (and considering that we often need to
 520      * compare events at end-points, this happens frequently enough to warrant
 521      * special casing).
 522      */
 523     enum {
 524        HAVE_NEITHER = 0x0,
 525        HAVE_AX      = 0x1,
 526        HAVE_BX      = 0x2,
 527        HAVE_BOTH    = HAVE_AX | HAVE_BX
 528     } have_ax_bx = HAVE_BOTH;
 529     int32_t ax, bx;
 530
 531     if (y == a->edge.line.p1.y)
 532         ax = a->edge.line.p1.x;
 533     else if (y == a->edge.line.p2.y)
 534         ax = a->edge.line.p2.x;
 535     else
 536         have_ax_bx &= ~HAVE_AX;
 537
 538     if (y == b->edge.line.p1.y)
 539         bx = b->edge.line.p1.x;
 540     else if (y == b->edge.line.p2.y)
 541         bx = b->edge.line.p2.x;
 542     else
 543         have_ax_bx &= ~HAVE_BX;
 544
 545     switch (have_ax_bx) {
 546     default:
 547     case HAVE_NEITHER:
 548         return edges_compare_x_for_y_general (a, b, y);
 549     case HAVE_AX:
 550         return -edge_compare_for_y_against_x (b, y, ax);
 551     case HAVE_BX:
 552         return edge_compare_for_y_against_x (a, y, bx);
 553     case HAVE_BOTH:
 554         return ax - bx;
 555     }
 556 }
 557
 558 static inline int
 559 _line_equal (const cairo_line_t *a, const cairo_line_t *b)
 560 {
 561     return a->p1.x == b->p1.x && a->p1.y == b->p1.y &&
 562            a->p2.x == b->p2.x && a->p2.y == b->p2.y;
 563 }
 564
 565 static inline int
 566 _cairo_bo_sweep_line_compare_edges (const cairo_bo_sweep_line_t *sweep_line,
 567                                     const cairo_bo_edge_t       *a,
 568                                     const cairo_bo_edge_t       *b)
 569 {
 570     int cmp;
 571
 572     /* compare the edges if not identical */
 573     if (! _line_equal (&a->edge.line, &b->edge.line)) {
 574         if (MAX (a->edge.line.p1.x, a->edge.line.p2.x) <
 575             MIN (b->edge.line.p1.x, b->edge.line.p2.x))
 576             return -1;
 577         else if (MIN (a->edge.line.p1.x, a->edge.line.p2.x) >
 578                  MAX (b->edge.line.p1.x, b->edge.line.p2.x))
 579             return 1;
 580
 581         cmp = edges_compare_x_for_y (a, b, sweep_line->current_y);
 582         if (cmp)
 583             return cmp;
 584
 585         /* The two edges intersect exactly at y, so fall back on slope
 586          * comparison. We know that this compare_edges function will be
 587          * called only when starting a new edge, (not when stopping an
 588          * edge), so we don't have to worry about conditionally inverting
 589          * the sense of _slope_compare. */
 590         cmp = _slope_compare (a, b);
 591         if (cmp)
 592             return cmp;
 593     }
 594
 595     /* We've got two collinear edges now. */
 596     return b->edge.bottom - a->edge.bottom;
 597 }
 598
 599 static inline cairo_int64_t
 600 det32_64 (int32_t a, int32_t b,
 601           int32_t c, int32_t d)
 602 {
 603     /* det = a * d - b * c */
 604     return _cairo_int64_sub (_cairo_int32x32_64_mul (a, d),
 605                              _cairo_int32x32_64_mul (b, c));
 606 }
 607
 608 static inline cairo_int128_t
 609 det64x32_128 (cairo_int64_t a, int32_t       b,
 610               cairo_int64_t c, int32_t       d)
 611 {
 612     /* det = a * d - b * c */
 613     return _cairo_int128_sub (_cairo_int64x32_128_mul (a, d),
 614                               _cairo_int64x32_128_mul (c, b));
 615 }
 616
 617 /* Compute the intersection of two lines as defined by two edges. The
 618  * result is provided as a coordinate pair of 128-bit integers.
 619  *
 620  * Returns %CAIRO_BO_STATUS_INTERSECTION if there is an intersection or
 621  * %CAIRO_BO_STATUS_PARALLEL if the two lines are exactly parallel.
 622  */
 623 static cairo_bool_t
 624 intersect_lines (cairo_bo_edge_t                *a,
 625                  cairo_bo_edge_t                *b,
 626                  cairo_bo_intersect_point_t     *intersection)
 627 {
 628     cairo_int64_t a_det, b_det;
 629
 630     /* XXX: We're assuming here that dx and dy will still fit in 32
 631      * bits. That's not true in general as there could be overflow. We
 632      * should prevent that before the tessellation algorithm begins.
 633      * What we're doing to mitigate this is to perform clamping in
 634      * cairo_bo_tessellate_polygon().
 635      */
 636     int32_t dx1 = a->edge.line.p1.x - a->edge.line.p2.x;
 637     int32_t dy1 = a->edge.line.p1.y - a->edge.line.p2.y;
 638
 639     int32_t dx2 = b->edge.line.p1.x - b->edge.line.p2.x;
 640     int32_t dy2 = b->edge.line.p1.y - b->edge.line.p2.y;
 641
 642     cairo_int64_t den_det;
 643     cairo_int64_t R;
 644     cairo_quorem64_t qr;
 645
 646     den_det = det32_64 (dx1, dy1, dx2, dy2);
 647
 648      /* Q: Can we determine that the lines do not intersect (within range)
 649       * much more cheaply than computing the intersection point i.e. by
 650       * avoiding the division?
 651       *
 652       *   X = ax + t * adx = bx + s * bdx;
 653       *   Y = ay + t * ady = by + s * bdy;
 654       *   ∴ t * (ady*bdx - bdy*adx) = bdx * (by - ay) + bdy * (ax - bx)
 655       *   => t * L = R
 656       *
 657       * Therefore we can reject any intersection (under the criteria for
 658       * valid intersection events) if:
 659       *   L^R < 0 => t < 0, or
 660       *   L<R => t > 1
 661       *
 662       * (where top/bottom must at least extend to the line endpoints).
 663       *
 664       * A similar substitution can be performed for s, yielding:
 665       *   s * (ady*bdx - bdy*adx) = ady * (ax - bx) - adx * (ay - by)
 666       */
 667     R = det32_64 (dx2, dy2,
 668                   b->edge.line.p1.x - a->edge.line.p1.x,
 669                   b->edge.line.p1.y - a->edge.line.p1.y);
 670     if (_cairo_int64_negative (den_det)) {
 671         if (_cairo_int64_ge (den_det, R))
 672             return FALSE;
 673     } else {
 674         if (_cairo_int64_le (den_det, R))
 675             return FALSE;
 676     }
 677
 678     R = det32_64 (dy1, dx1,
 679                   a->edge.line.p1.y - b->edge.line.p1.y,
 680                   a->edge.line.p1.x - b->edge.line.p1.x);
 681     if (_cairo_int64_negative (den_det)) {
 682         if (_cairo_int64_ge (den_det, R))
 683             return FALSE;
 684     } else {
 685         if (_cairo_int64_le (den_det, R))
 686             return FALSE;
 687     }
 688
 689     /* We now know that the two lines should intersect within range. */
 690
 691     a_det = det32_64 (a->edge.line.p1.x, a->edge.line.p1.y,
 692                       a->edge.line.p2.x, a->edge.line.p2.y);
 693     b_det = det32_64 (b->edge.line.p1.x, b->edge.line.p1.y,
 694                       b->edge.line.p2.x, b->edge.line.p2.y);
 695
 696     /* x = det (a_det, dx1, b_det, dx2) / den_det */
 697     qr = _cairo_int_96by64_32x64_divrem (det64x32_128 (a_det, dx1,
 698                                                        b_det, dx2),
 699                                          den_det);
 700     if (_cairo_int64_eq (qr.rem, den_det))
 701         return FALSE;
 702 #if 0
 703     intersection->x.exactness = _cairo_int64_is_zero (qr.rem) ? EXACT : INEXACT;
 704 #else
 705     intersection->x.exactness = EXACT;
 706     if (! _cairo_int64_is_zero (qr.rem)) {
 707         if (_cairo_int64_negative (den_det) ^ _cairo_int64_negative (qr.rem))
 708             qr.rem = _cairo_int64_negate (qr.rem);
 709         qr.rem = _cairo_int64_mul (qr.rem, _cairo_int32_to_int64 (2));
 710         if (_cairo_int64_ge (qr.rem, den_det)) {
 711             qr.quo = _cairo_int64_add (qr.quo,
 712                                        _cairo_int32_to_int64 (_cairo_int64_negative (qr.quo) ? -1 : 1));
 713         } else
 714             intersection->x.exactness = INEXACT;
 715     }
 716 #endif
 717     intersection->x.ordinate = _cairo_int64_to_int32 (qr.quo);
 718
 719     /* y = det (a_det, dy1, b_det, dy2) / den_det */
 720     qr = _cairo_int_96by64_32x64_divrem (det64x32_128 (a_det, dy1,
 721                                                        b_det, dy2),
 722                                          den_det);
 723     if (_cairo_int64_eq (qr.rem, den_det))
 724         return FALSE;
 725 #if 0
 726     intersection->y.exactness = _cairo_int64_is_zero (qr.rem) ? EXACT : INEXACT;
 727 #else
 728     intersection->y.exactness = EXACT;
 729     if (! _cairo_int64_is_zero (qr.rem)) {
 730         if (_cairo_int64_negative (den_det) ^ _cairo_int64_negative (qr.rem))
 731             qr.rem = _cairo_int64_negate (qr.rem);
 732         qr.rem = _cairo_int64_mul (qr.rem, _cairo_int32_to_int64 (2));
 733         if (_cairo_int64_ge (qr.rem, den_det)) {
 734             qr.quo = _cairo_int64_add (qr.quo,
 735                                        _cairo_int32_to_int64 (_cairo_int64_negative (qr.quo) ? -1 : 1));
 736         } else
 737             intersection->y.exactness = INEXACT;
 738     }
 739 #endif
 740     intersection->y.ordinate = _cairo_int64_to_int32 (qr.quo);
 741
 742     return TRUE;
 743 }
 744
 745 static int
 746 _cairo_bo_intersect_ordinate_32_compare (cairo_bo_intersect_ordinate_t  a,
 747                                          int32_t                        b)
 748 {
 749     /* First compare the quotient */
 750     if (a.ordinate > b)
 751         return +1;
 752     if (a.ordinate < b)
 753         return -1;
 754     /* With quotient identical, if remainder is 0 then compare equal */
 755     /* Otherwise, the non-zero remainder makes a > b */
 756     return INEXACT == a.exactness;
 757 }
 758
 759 /* Does the given edge contain the given point. The point must already
 760  * be known to be contained within the line determined by the edge,
 761  * (most likely the point results from an intersection of this edge
 762  * with another).
 763  *
 764  * If we had exact arithmetic, then this function would simply be a
 765  * matter of examining whether the y value of the point lies within
 766  * the range of y values of the edge. But since intersection points
 767  * are not exact due to being rounded to the nearest integer within
 768  * the available precision, we must also examine the x value of the
 769  * point.
 770  *
 771  * The definition of "contains" here is that the given intersection
 772  * point will be seen by the sweep line after the start event for the
 773  * given edge and before the stop event for the edge. See the comments
 774  * in the implementation for more details.
 775  */
 776 static cairo_bool_t
 777 _cairo_bo_edge_contains_intersect_point (cairo_bo_edge_t                *edge,
 778                                          cairo_bo_intersect_point_t     *point)
 779 {
 780     int cmp_top, cmp_bottom;
 781
 782     /* XXX: When running the actual algorithm, we don't actually need to
 783      * compare against edge->top at all here, since any intersection above
 784      * top is eliminated early via a slope comparison. We're leaving these
 785      * here for now only for the sake of the quadratic-time intersection
 786      * finder which needs them.
 787      */
 788
 789     cmp_top = _cairo_bo_intersect_ordinate_32_compare (point->y,
 790                                                        edge->edge.top);
 791     cmp_bottom = _cairo_bo_intersect_ordinate_32_compare (point->y,
 792                                                           edge->edge.bottom);
 793
 794     if (cmp_top < 0 || cmp_bottom > 0)
 795     {
 796         return FALSE;
 797     }
 798
 799     if (cmp_top > 0 && cmp_bottom < 0)
 800     {
 801         return TRUE;
 802     }
 803
 804     /* At this stage, the point lies on the same y value as either
 805      * edge->top or edge->bottom, so we have to examine the x value in
 806      * order to properly determine containment. */
 807
 808     /* If the y value of the point is the same as the y value of the
 809      * top of the edge, then the x value of the point must be greater
 810      * to be considered as inside the edge. Similarly, if the y value
 811      * of the point is the same as the y value of the bottom of the
 812      * edge, then the x value of the point must be less to be
 813      * considered as inside. */
 814
 815     if (cmp_top == 0) {
 816         cairo_fixed_t top_x;
 817
 818         top_x = _line_compute_intersection_x_for_y (&edge->edge.line,
 819                                                     edge->edge.top);
 820         return _cairo_bo_intersect_ordinate_32_compare (point->x, top_x) > 0;
 821     } else { /* cmp_bottom == 0 */
 822         cairo_fixed_t bot_x;
 823
 824         bot_x = _line_compute_intersection_x_for_y (&edge->edge.line,
 825                                                     edge->edge.bottom);
 826         return _cairo_bo_intersect_ordinate_32_compare (point->x, bot_x) < 0;
 827     }
 828 }
 829
 830 /* Compute the intersection of two edges. The result is provided as a
 831  * coordinate pair of 128-bit integers.
 832  *
 833  * Returns %CAIRO_BO_STATUS_INTERSECTION if there is an intersection
 834  * that is within both edges, %CAIRO_BO_STATUS_NO_INTERSECTION if the
 835  * intersection of the lines defined by the edges occurs outside of
 836  * one or both edges, and %CAIRO_BO_STATUS_PARALLEL if the two edges
 837  * are exactly parallel.
 838  *
 839  * Note that when determining if a candidate intersection is "inside"
 840  * an edge, we consider both the infinitesimal shortening and the
 841  * infinitesimal tilt rules described by John Hobby. Specifically, if
 842  * the intersection is exactly the same as an edge point, it is
 843  * effectively outside (no intersection is returned). Also, if the
 844  * intersection point has the same
 845  */
 846 static cairo_bool_t
 847 _cairo_bo_edge_intersect (cairo_bo_edge_t       *a,
 848                           cairo_bo_edge_t       *b,
 849                           cairo_bo_point32_t    *intersection)
 850 {
 851     cairo_bo_intersect_point_t quorem;
 852
 853     if (! intersect_lines (a, b, &quorem))
 854         return FALSE;
 855
 856     if (! _cairo_bo_edge_contains_intersect_point (a, &quorem))
 857         return FALSE;
 858
 859     if (! _cairo_bo_edge_contains_intersect_point (b, &quorem))
 860         return FALSE;
 861
 862     /* Now that we've correctly compared the intersection point and
 863      * determined that it lies within the edge, then we know that we
 864      * no longer need any more bits of storage for the intersection
 865      * than we do for our edge coordinates. We also no longer need the
 866      * remainder from the division. */
 867     intersection->x = quorem.x.ordinate;
 868     intersection->y = quorem.y.ordinate;
 869
 870     return TRUE;
 871 }
 872
 873 static inline int
 874 cairo_bo_event_compare (const cairo_bo_event_t *a,
 875                         const cairo_bo_event_t *b)
 876 {
 877     int cmp;
 878
 879     cmp = _cairo_bo_point32_compare (&a->point, &b->point);
 880     if (cmp)
 881         return cmp;
 882
 883     cmp = a->type - b->type;
 884     if (cmp)
 885         return cmp;
 886
 887     return a - b;
 888 }
 889
 890 static inline void
 891 _pqueue_init (pqueue_t *pq)
 892 {
 893     pq->max_size = ARRAY_LENGTH (pq->elements_embedded);
 894     pq->size = 0;
 895
 896     pq->elements = pq->elements_embedded;
 897 }
 898
 899 static inline void
 900 _pqueue_fini (pqueue_t *pq)
 901 {
 902     if (pq->elements != pq->elements_embedded)
 903         free (pq->elements);
 904 }
 905
 906 static cairo_status_t
 907 _pqueue_grow (pqueue_t *pq)
 908 {
 909     cairo_bo_event_t **new_elements;
 910     pq->max_size *= 2;
 911
 912     if (pq->elements == pq->elements_embedded) {
 913         new_elements = _cairo_malloc_ab (pq->max_size,
 914                                          sizeof (cairo_bo_event_t *));
 915         if (unlikely (new_elements == NULL))
 916             return _cairo_error (CAIRO_STATUS_NO_MEMORY);
 917
 918         memcpy (new_elements, pq->elements_embedded,
 919                 sizeof (pq->elements_embedded));
 920     } else {
 921         new_elements = _cairo_realloc_ab (pq->elements,
 922                                           pq->max_size,
 923                                           sizeof (cairo_bo_event_t *));
 924         if (unlikely (new_elements == NULL))
 925             return _cairo_error (CAIRO_STATUS_NO_MEMORY);
 926     }
 927
 928     pq->elements = new_elements;
 929     return CAIRO_STATUS_SUCCESS;
 930 }
 931
 932 static inline cairo_status_t
 933 _pqueue_push (pqueue_t *pq, cairo_bo_event_t *event)
 934 {
 935     cairo_bo_event_t **elements;
 936     int i, parent;
 937
 938     if (unlikely (pq->size + 1 == pq->max_size)) {
 939         cairo_status_t status;
 940
 941         status = _pqueue_grow (pq);
 942         if (unlikely (status))
 943             return status;
 944     }
 945
 946     elements = pq->elements;
 947
 948     for (i = ++pq->size;
 949          i != PQ_FIRST_ENTRY &&
 950          cairo_bo_event_compare (event,
 951                                  elements[parent = PQ_PARENT_INDEX (i)]) < 0;
 952          i = parent)
 953     {
 954         elements[i] = elements[parent];
 955     }
 956
 957     elements[i] = event;
 958
 959     return CAIRO_STATUS_SUCCESS;
 960 }
 961
 962 static inline void
 963 _pqueue_pop (pqueue_t *pq)
 964 {
 965     cairo_bo_event_t **elements = pq->elements;
 966     cairo_bo_event_t *tail;
 967     int child, i;
 968
 969     tail = elements[pq->size--];
 970     if (pq->size == 0) {
 971         elements[PQ_FIRST_ENTRY] = NULL;
 972         return;
 973     }
 974
 975     for (i = PQ_FIRST_ENTRY;
 976          (child = PQ_LEFT_CHILD_INDEX (i)) <= pq->size;
 977          i = child)
 978     {
 979         if (child != pq->size &&
 980             cairo_bo_event_compare (elements[child+1],
 981                                     elements[child]) < 0)
 982         {
 983             child++;
 984         }
 985
 986         if (cairo_bo_event_compare (elements[child], tail) >= 0)
 987             break;
 988
 989         elements[i] = elements[child];
 990     }
 991     elements[i] = tail;
 992 }
 993
 994 static inline cairo_status_t
 995 _cairo_bo_event_queue_insert (cairo_bo_event_queue_t    *queue,
 996                               cairo_bo_event_type_t      type,
 997                               cairo_bo_edge_t           *e1,
 998                               cairo_bo_edge_t           *e2,
 999                               const cairo_point_t        *point)
1000 {
1001     cairo_bo_queue_event_t *event;
1002
1003     event = _cairo_freepool_alloc (&queue->pool);
1004     if (unlikely (event == NULL))
1005         return _cairo_error (CAIRO_STATUS_NO_MEMORY);
1006
1007     event->type = type;
1008     event->e1 = e1;
1009     event->e2 = e2;
1010     event->point = *point;
1011
1012     return _pqueue_push (&queue->pqueue, (cairo_bo_event_t *) event);
1013 }
1014
1015 static void
1016 _cairo_bo_event_queue_delete (cairo_bo_event_queue_t *queue,
1017                               cairo_bo_event_t       *event)
1018 {
1019     _cairo_freepool_free (&queue->pool, event);
1020 }
1021
1022 static cairo_bo_event_t *
1023 _cairo_bo_event_dequeue (cairo_bo_event_queue_t *event_queue)
1024 {
1025     cairo_bo_event_t *event, *cmp;
1026
1027     event = event_queue->pqueue.elements[PQ_FIRST_ENTRY];
1028     cmp = *event_queue->start_events;
1029     if (event == NULL ||
1030         (cmp != NULL && cairo_bo_event_compare (cmp, event) < 0))
1031     {
1032         event = cmp;
1033         event_queue->start_events++;
1034     }
1035     else
1036     {
1037         _pqueue_pop (&event_queue->pqueue);
1038     }
1039
1040     return event;
1041 }
1042
1043 CAIRO_COMBSORT_DECLARE (_cairo_bo_event_queue_sort,
1044                         cairo_bo_event_t *,
1045                         cairo_bo_event_compare)
1046
1047 static void
1048 _cairo_bo_event_queue_init (cairo_bo_event_queue_t       *event_queue,
1049                             cairo_bo_event_t            **start_events,
1050                             int                           num_events)
1051 {
1052     event_queue->start_events = start_events;
1053
1054     _cairo_freepool_init (&event_queue->pool,
1055                           sizeof (cairo_bo_queue_event_t));
1056     _pqueue_init (&event_queue->pqueue);
1057     event_queue->pqueue.elements[PQ_FIRST_ENTRY] = NULL;
1058 }
1059
1060 static cairo_status_t
1061 _cairo_bo_event_queue_insert_stop (cairo_bo_event_queue_t       *event_queue,
1062                                    cairo_bo_edge_t              *edge)
1063 {
1064     cairo_bo_point32_t point;
1065
1066     point.y = edge->edge.bottom;
1067     point.x = _line_compute_intersection_x_for_y (&edge->edge.line,
1068                                                   point.y);
1069     return _cairo_bo_event_queue_insert (event_queue,
1070                                          CAIRO_BO_EVENT_TYPE_STOP,
1071                                          edge, NULL,
1072                                          &point);
1073 }
1074
1075 static void
1076 _cairo_bo_event_queue_fini (cairo_bo_event_queue_t *event_queue)
1077 {
1078     _pqueue_fini (&event_queue->pqueue);
1079     _cairo_freepool_fini (&event_queue->pool);
1080 }
1081
1082 static inline cairo_status_t
1083 _cairo_bo_event_queue_insert_if_intersect_below_current_y (cairo_bo_event_queue_t       *event_queue,
1084                                                            cairo_bo_edge_t      *left,
1085                                                            cairo_bo_edge_t *right)
1086 {
1087     cairo_bo_point32_t intersection;
1088
1089     if (MAX (left->edge.line.p1.x, left->edge.line.p2.x) <=
1090         MIN (right->edge.line.p1.x, right->edge.line.p2.x))
1091         return CAIRO_STATUS_SUCCESS;
1092
1093     if (_line_equal (&left->edge.line, &right->edge.line))
1094         return CAIRO_STATUS_SUCCESS;
1095
1096     /* The names "left" and "right" here are correct descriptions of
1097      * the order of the two edges within the active edge list. So if a
1098      * slope comparison also puts left less than right, then we know
1099      * that the intersection of these two segments has already
1100      * occurred before the current sweep line position. */
1101     if (_slope_compare (left, right) <= 0)
1102         return CAIRO_STATUS_SUCCESS;
1103
1104     if (! _cairo_bo_edge_intersect (left, right, &intersection))
1105         return CAIRO_STATUS_SUCCESS;
1106
1107     return _cairo_bo_event_queue_insert (event_queue,
1108                                          CAIRO_BO_EVENT_TYPE_INTERSECTION,
1109                                          left, right,
1110                                          &intersection);
1111 }
1112
1113 static void
1114 _cairo_bo_sweep_line_init (cairo_bo_sweep_line_t *sweep_line)
1115 {
1116     sweep_line->head = NULL;
1117     sweep_line->stopped = NULL;
1118     sweep_line->current_y = INT32_MIN;
1119     sweep_line->current_edge = NULL;
1120 }
1121
1122 static void
1123 _cairo_bo_sweep_line_insert (cairo_bo_sweep_line_t      *sweep_line,
1124                              cairo_bo_edge_t            *edge)
1125 {
1126     if (sweep_line->current_edge != NULL) {
1127         cairo_bo_edge_t *prev, *next;
1128         int cmp;
1129
1130         cmp = _cairo_bo_sweep_line_compare_edges (sweep_line,
1131                                                   sweep_line->current_edge,
1132                                                   edge);
1133         if (cmp < 0) {
1134             prev = sweep_line->current_edge;
1135             next = prev->next;
1136             while (next != NULL &&
1137                    _cairo_bo_sweep_line_compare_edges (sweep_line,
1138                                                        next, edge) < 0)
1139             {
1140                 prev = next, next = prev->next;
1141             }
1142
1143             prev->next = edge;
1144             edge->prev = prev;
1145             edge->next = next;
1146             if (next != NULL)
1147                 next->prev = edge;
1148         } else if (cmp > 0) {
1149             next = sweep_line->current_edge;
1150             prev = next->prev;
1151             while (prev != NULL &&
1152                    _cairo_bo_sweep_line_compare_edges (sweep_line,
1153                                                        prev, edge) > 0)
1154             {
1155                 next = prev, prev = next->prev;
1156             }
1157
1158             next->prev = edge;
1159             edge->next = next;
1160             edge->prev = prev;
1161             if (prev != NULL)
1162                 prev->next = edge;
1163             else
1164                 sweep_line->head = edge;
1165         } else {
1166             prev = sweep_line->current_edge;
1167             edge->prev = prev;
1168             edge->next = prev->next;
1169             if (prev->next != NULL)
1170                 prev->next->prev = edge;
1171             prev->next = edge;
1172         }
1173     } else {
1174         sweep_line->head = edge;
1175         edge->next = NULL;
1176     }
1177
1178     sweep_line->current_edge = edge;
1179 }
1180
1181 static void
1182 _cairo_bo_sweep_line_delete (cairo_bo_sweep_line_t      *sweep_line,
1183                              cairo_bo_edge_t    *edge)
1184 {
1185     if (edge->prev != NULL)
1186         edge->prev->next = edge->next;
1187     else
1188         sweep_line->head = edge->next;
1189
1190     if (edge->next != NULL)
1191         edge->next->prev = edge->prev;
1192
1193     if (sweep_line->current_edge == edge)
1194         sweep_line->current_edge = edge->prev ? edge->prev : edge->next;
1195 }
1196
1197 static void
1198 _cairo_bo_sweep_line_swap (cairo_bo_sweep_line_t        *sweep_line,
1199                            cairo_bo_edge_t              *left,
1200                            cairo_bo_edge_t              *right)
1201 {
1202     if (left->prev != NULL)
1203         left->prev->next = right;
1204     else
1205         sweep_line->head = right;
1206
1207     if (right->next != NULL)
1208         right->next->prev = left;
1209
1210     left->next = right->next;
1211     right->next = left;
1212
1213     right->prev = left->prev;
1214     left->prev = right;
1215 }
1216
1217 #if DEBUG_PRINT_STATE
1218 static void
1219 _cairo_bo_edge_print (cairo_bo_edge_t *edge)
1220 {
1221     printf ("(0x%x, 0x%x)-(0x%x, 0x%x)",
1222             edge->edge.line.p1.x, edge->edge.line.p1.y,
1223             edge->edge.line.p2.x, edge->edge.line.p2.y);
1224 }
1225
1226 static void
1227 _cairo_bo_event_print (cairo_bo_event_t *event)
1228 {
1229     switch (event->type) {
1230     case CAIRO_BO_EVENT_TYPE_START:
1231         printf ("Start: ");
1232         break;
1233     case CAIRO_BO_EVENT_TYPE_STOP:
1234         printf ("Stop: ");
1235         break;
1236     case CAIRO_BO_EVENT_TYPE_INTERSECTION:
1237         printf ("Intersection: ");
1238         break;
1239     }
1240     printf ("(%d, %d)\t", event->point.x, event->point.y);
1241     _cairo_bo_edge_print (event->e1);
1242     if (event->type == CAIRO_BO_EVENT_TYPE_INTERSECTION) {
1243         printf (" X ");
1244         _cairo_bo_edge_print (event->e2);
1245     }
1246     printf ("\n");
1247 }
1248
1249 static void
1250 _cairo_bo_event_queue_print (cairo_bo_event_queue_t *event_queue)
1251 {
1252     /* XXX: fixme to print the start/stop array too. */
1253     printf ("Event queue:\n");
1254 }
1255
1256 static void
1257 _cairo_bo_sweep_line_print (cairo_bo_sweep_line_t *sweep_line)
1258 {
1259     cairo_bool_t first = TRUE;
1260     cairo_bo_edge_t *edge;
1261
1262     printf ("Sweep line from edge list: ");
1263     first = TRUE;
1264     for (edge = sweep_line->head;
1265          edge;
1266          edge = edge->next)
1267     {
1268         if (!first)
1269             printf (", ");
1270         _cairo_bo_edge_print (edge);
1271         first = FALSE;
1272     }
1273     printf ("\n");
1274 }
1275
1276 static void
1277 print_state (const char                 *msg,
1278              cairo_bo_event_t           *event,
1279              cairo_bo_event_queue_t     *event_queue,
1280              cairo_bo_sweep_line_t      *sweep_line)
1281 {
1282     printf ("%s ", msg);
1283     _cairo_bo_event_print (event);
1284     _cairo_bo_event_queue_print (event_queue);
1285     _cairo_bo_sweep_line_print (sweep_line);
1286     printf ("\n");
1287 }
1288 #endif
1289
1290 #if DEBUG_EVENTS
1291 static void CAIRO_PRINTF_FORMAT (1, 2)
1292 event_log (const char *fmt, ...)
1293 {
1294     FILE *file;
1295
1296     if (getenv ("CAIRO_DEBUG_EVENTS") == NULL)
1297         return;
1298
1299     file = fopen ("bo-events.txt", "a");
1300     if (file != NULL) {
1301         va_list ap;
1302
1303         va_start (ap, fmt);
1304         vfprintf (file, fmt, ap);
1305         va_end (ap);
1306
1307         fclose (file);
1308     }
1309 }
1310 #endif
1311
1312 #define HAS_COLINEAR(a, b) ((cairo_bo_edge_t *)(((uintptr_t)(a))&~1) == (b))
1313 #define IS_COLINEAR(e) (((uintptr_t)(e))&1)
1314 #define MARK_COLINEAR(e, v) ((cairo_bo_edge_t *)(((uintptr_t)(e))|(v)))
1315
1316 static inline cairo_bool_t
1317 edges_colinear (cairo_bo_edge_t *a, const cairo_bo_edge_t *b)
1318 {
1319     unsigned p;
1320
1321     if (HAS_COLINEAR(a->colinear, b))
1322         return IS_COLINEAR(a->colinear);
1323
1324     if (HAS_COLINEAR(b->colinear, a)) {
1325         p = IS_COLINEAR(b->colinear);
1326         a->colinear = MARK_COLINEAR(b, p);
1327         return p;
1328     }
1329
1330     p = 0;
1331     p |= (a->edge.line.p1.x == b->edge.line.p1.x) << 0;
1332     p |= (a->edge.line.p1.y == b->edge.line.p1.y) << 1;
1333     p |= (a->edge.line.p2.x == b->edge.line.p2.x) << 3;
1334     p |= (a->edge.line.p2.y == b->edge.line.p2.y) << 4;
1335     if (p == ((1 << 0) | (1 << 1) | (1 << 3) | (1 << 4))) {
1336         a->colinear = MARK_COLINEAR(b, 1);
1337         return TRUE;
1338     }
1339
1340     if (_slope_compare (a, b)) {
1341         a->colinear = MARK_COLINEAR(b, 0);
1342         return FALSE;
1343     }
1344
1345     /* The choice of y is not truly arbitrary since we must guarantee that it
1346      * is greater than the start of either line.
1347      */
1348     if (p != 0) {
1349         /* colinear if either end-point are coincident */
1350         p = (((p >> 1) & p) & 5) != 0;
1351     } else if (a->edge.line.p1.y < b->edge.line.p1.y) {
1352         p = edge_compare_for_y_against_x (b,
1353                                           a->edge.line.p1.y,
1354                                           a->edge.line.p1.x) == 0;
1355     } else {
1356         p = edge_compare_for_y_against_x (a,
1357                                           b->edge.line.p1.y,
1358                                           b->edge.line.p1.x) == 0;
1359     }
1360
1361     a->colinear = MARK_COLINEAR(b, p);
1362     return p;
1363 }
1364
1365 /* Adds the trapezoid, if any, of the left edge to the #cairo_traps_t */
1366 static void
1367 _cairo_bo_edge_end_trap (cairo_bo_edge_t        *left,
1368                          int32_t                 bot,
1369                          cairo_traps_t          *traps)
1370 {
1371     cairo_bo_trap_t *trap = &left->deferred_trap;
1372
1373     /* Only emit (trivial) non-degenerate trapezoids with positive height. */
1374     if (likely (trap->top < bot)) {
1375         _cairo_traps_add_trap (traps,
1376                                trap->top, bot,
1377                                &left->edge.line, &trap->right->edge.line);
1378
1379 #if DEBUG_PRINT_STATE
1380         printf ("Deferred trap: left=(%x, %x)-(%x,%x) "
1381                 "right=(%x,%x)-(%x,%x) top=%x, bot=%x\n",
1382                 left->edge.line.p1.x, left->edge.line.p1.y,
1383                 left->edge.line.p2.x, left->edge.line.p2.y,
1384                 trap->right->edge.line.p1.x, trap->right->edge.line.p1.y,
1385                 trap->right->edge.line.p2.x, trap->right->edge.line.p2.y,
1386                 trap->top, bot);
1387 #endif
1388 #if DEBUG_EVENTS
1389         event_log ("end trap: %lu %lu %d %d\n",
1390                    (long) left,
1391                    (long) trap->right,
1392                    trap->top,
1393                    bot);
1394 #endif
1395     }
1396
1397     trap->right = NULL;
1398 }
1399
1400
1401 /* Start a new trapezoid at the given top y coordinate, whose edges
1402  * are `edge' and `edge->next'. If `edge' already has a trapezoid,
1403  * then either add it to the traps in `traps', if the trapezoid's
1404  * right edge differs from `edge->next', or do nothing if the new
1405  * trapezoid would be a continuation of the existing one. */
1406 static inline void
1407 _cairo_bo_edge_start_or_continue_trap (cairo_bo_edge_t  *left,
1408                                        cairo_bo_edge_t  *right,
1409                                        int               top,
1410                                        cairo_traps_t    *traps)
1411 {
1412     if (left->deferred_trap.right == right)
1413         return;
1414
1415     assert (right);
1416     if (left->deferred_trap.right != NULL) {
1417         if (edges_colinear (left->deferred_trap.right, right))
1418         {
1419             /* continuation on right, so just swap edges */
1420             left->deferred_trap.right = right;
1421             return;
1422         }
1423
1424         _cairo_bo_edge_end_trap (left, top, traps);
1425     }
1426
1427     if (! edges_colinear (left, right)) {
1428         left->deferred_trap.top = top;
1429         left->deferred_trap.right = right;
1430
1431 #if DEBUG_EVENTS
1432         event_log ("begin trap: %lu %lu %d\n",
1433                    (long) left,
1434                    (long) right,
1435                    top);
1436 #endif
1437     }
1438 }
1439
1440 static inline void
1441 _active_edges_to_traps (cairo_bo_edge_t *pos,
1442                         int32_t          top,
1443                         unsigned         mask,
1444                         cairo_traps_t        *traps)
1445 {
1446     cairo_bo_edge_t *left;
1447     int in_out;
1448
1449
1450 #if DEBUG_PRINT_STATE
1451     printf ("Processing active edges for %x\n", top);
1452 #endif
1453
1454     in_out = 0;
1455     left = pos;
1456     while (pos != NULL) {
1457         if (pos != left && pos->deferred_trap.right) {
1458             /* XXX It shouldn't be possible to here with 2 deferred traps
1459              * on colinear edges... See bug-bo-rictoz.
1460              */
1461             if (left->deferred_trap.right == NULL &&
1462                 edges_colinear (left, pos))
1463             {
1464                 /* continuation on left */
1465                 left->deferred_trap = pos->deferred_trap;
1466                 pos->deferred_trap.right = NULL;
1467             }
1468             else
1469             {
1470                 _cairo_bo_edge_end_trap (pos, top, traps);
1471             }
1472         }
1473
1474         in_out += pos->edge.dir;
1475         if ((in_out & mask) == 0) {
1476             /* skip co-linear edges */
1477             if (pos->next == NULL || ! edges_colinear (pos, pos->next)) {
1478                 _cairo_bo_edge_start_or_continue_trap (left, pos, top, traps);
1479                 left = pos->next;
1480             }
1481         }
1482
1483         pos = pos->next;
1484     }
1485 }
1486
1487 /* Execute a single pass of the Bentley-Ottmann algorithm on edges,
1488  * generating trapezoids according to the fill_rule and appending them
1489  * to traps. */
1490 static cairo_status_t
1491 _cairo_bentley_ottmann_tessellate_bo_edges (cairo_bo_event_t   **start_events,
1492                                             int                  num_events,
1493                                             unsigned             fill_rule,
1494                                             cairo_traps_t       *traps,
1495                                             int                 *num_intersections)
1496 {
1497     cairo_status_t status;
1498     int intersection_count = 0;
1499     cairo_bo_event_queue_t event_queue;
1500     cairo_bo_sweep_line_t sweep_line;
1501     cairo_bo_event_t *event;
1502     cairo_bo_edge_t *left, *right;
1503     cairo_bo_edge_t *e1, *e2;
1504
1505     /* convert the fill_rule into a winding mask */
1506     if (fill_rule == CAIRO_FILL_RULE_WINDING)
1507         fill_rule = (unsigned) -1;
1508     else
1509         fill_rule = 1;
1510
1511 #if DEBUG_EVENTS
1512     {
1513         int i;
1514
1515         for (i = 0; i < num_events; i++) {
1516             cairo_bo_start_event_t *event =
1517                 ((cairo_bo_start_event_t **) start_events)[i];
1518             event_log ("edge: %lu (%d, %d) (%d, %d) (%d, %d) %d\n",
1519                        (long) &events[i].edge,
1520                        event->edge.edge.line.p1.x,
1521                        event->edge.edge.line.p1.y,
1522                        event->edge.edge.line.p2.x,
1523                        event->edge.edge.line.p2.y,
1524                        event->edge.top,
1525                        event->edge.bottom,
1526                        event->edge.edge.dir);
1527         }
1528     }
1529 #endif
1530
1531     _cairo_bo_event_queue_init (&event_queue, start_events, num_events);
1532     _cairo_bo_sweep_line_init (&sweep_line);
1533
1534     while ((event = _cairo_bo_event_dequeue (&event_queue))) {
1535         if (event->point.y != sweep_line.current_y) {
1536             for (e1 = sweep_line.stopped; e1; e1 = e1->next) {
1537                 if (e1->deferred_trap.right != NULL) {
1538                     _cairo_bo_edge_end_trap (e1,
1539                                              e1->edge.bottom,
1540                                              traps);
1541                 }
1542             }
1543             sweep_line.stopped = NULL;
1544
1545             _active_edges_to_traps (sweep_line.head,
1546                                     sweep_line.current_y,
1547                                     fill_rule, traps);
1548
1549             sweep_line.current_y = event->point.y;
1550         }
1551
1552 #if DEBUG_EVENTS
1553         event_log ("event: %d (%ld, %ld) %lu, %lu\n",
1554                    event->type,
1555                    (long) event->point.x,
1556                    (long) event->point.y,
1557                    (long) event->e1,
1558                    (long) event->e2);
1559 #endif
1560
1561         switch (event->type) {
1562         case CAIRO_BO_EVENT_TYPE_START:
1563             e1 = &((cairo_bo_start_event_t *) event)->edge;
1564
1565             _cairo_bo_sweep_line_insert (&sweep_line, e1);
1566
1567             status = _cairo_bo_event_queue_insert_stop (&event_queue, e1);
1568             if (unlikely (status))
1569                 goto unwind;
1570
1571             /* check to see if this is a continuation of a stopped edge */
1572             /* XXX change to an infinitesimal lengthening rule */
1573             for (left = sweep_line.stopped; left; left = left->next) {
1574                 if (e1->edge.top <= left->edge.bottom &&
1575                     edges_colinear (e1, left))
1576                 {
1577                     e1->deferred_trap = left->deferred_trap;
1578                     if (left->prev != NULL)
1579                         left->prev = left->next;
1580                     else
1581                         sweep_line.stopped = left->next;
1582                     if (left->next != NULL)
1583                         left->next->prev = left->prev;
1584                     break;
1585                 }
1586             }
1587
1588             left = e1->prev;
1589             right = e1->next;
1590
1591             if (left != NULL) {
1592                 status = _cairo_bo_event_queue_insert_if_intersect_below_current_y (&event_queue, left, e1);
1593                 if (unlikely (status))
1594                     goto unwind;
1595             }
1596
1597             if (right != NULL) {
1598                 status = _cairo_bo_event_queue_insert_if_intersect_below_current_y (&event_queue, e1, right);
1599                 if (unlikely (status))
1600                     goto unwind;
1601             }
1602
1603             break;
1604
1605         case CAIRO_BO_EVENT_TYPE_STOP:
1606             e1 = ((cairo_bo_queue_event_t *) event)->e1;
1607             _cairo_bo_event_queue_delete (&event_queue, event);
1608
1609             left = e1->prev;
1610             right = e1->next;
1611
1612             _cairo_bo_sweep_line_delete (&sweep_line, e1);
1613
1614             /* first, check to see if we have a continuation via a fresh edge */
1615             if (e1->deferred_trap.right != NULL) {
1616                 e1->next = sweep_line.stopped;
1617                 if (sweep_line.stopped != NULL)
1618                     sweep_line.stopped->prev = e1;
1619                 sweep_line.stopped = e1;
1620                 e1->prev = NULL;
1621             }
1622
1623             if (left != NULL && right != NULL) {
1624                 status = _cairo_bo_event_queue_insert_if_intersect_below_current_y (&event_queue, left, right);
1625                 if (unlikely (status))
1626                     goto unwind;
1627             }
1628
1629             break;
1630
1631         case CAIRO_BO_EVENT_TYPE_INTERSECTION:
1632             e1 = ((cairo_bo_queue_event_t *) event)->e1;
1633             e2 = ((cairo_bo_queue_event_t *) event)->e2;
1634             _cairo_bo_event_queue_delete (&event_queue, event);
1635
1636             /* skip this intersection if its edges are not adjacent */
1637             if (e2 != e1->next)
1638                 break;
1639
1640             intersection_count++;
1641
1642             left = e1->prev;
1643             right = e2->next;
1644
1645             _cairo_bo_sweep_line_swap (&sweep_line, e1, e2);
1646
1647             /* after the swap e2 is left of e1 */
1648
1649             if (left != NULL) {
1650                 status = _cairo_bo_event_queue_insert_if_intersect_below_current_y (&event_queue, left, e2);
1651                 if (unlikely (status))
1652                     goto unwind;
1653             }
1654
1655             if (right != NULL) {
1656                 status = _cairo_bo_event_queue_insert_if_intersect_below_current_y (&event_queue, e1, right);
1657                 if (unlikely (status))
1658                     goto unwind;
1659             }
1660
1661             break;
1662         }
1663     }
1664
1665     *num_intersections = intersection_count;
1666     for (e1 = sweep_line.stopped; e1; e1 = e1->next) {
1667         if (e1->deferred_trap.right != NULL) {
1668             _cairo_bo_edge_end_trap (e1, e1->edge.bottom, traps);
1669         }
1670     }
1671     status = traps->status;
1672  unwind:
1673     _cairo_bo_event_queue_fini (&event_queue);
1674
1675 #if DEBUG_EVENTS
1676     event_log ("\n");
1677 #endif
1678
1679     return status;
1680 }
1681
1682 cairo_status_t
1683 _cairo_bentley_ottmann_tessellate_polygon (cairo_traps_t         *traps,
1684                                            const cairo_polygon_t *polygon,
1685                                            cairo_fill_rule_t      fill_rule)
1686 {
1687     int intersections;
1688     cairo_bo_start_event_t stack_events[CAIRO_STACK_ARRAY_LENGTH (cairo_bo_start_event_t)];
1689     cairo_bo_start_event_t *events;
1690     cairo_bo_event_t *stack_event_ptrs[ARRAY_LENGTH (stack_events) + 1];
1691     cairo_bo_event_t **event_ptrs;
1692     cairo_bo_start_event_t *stack_event_y[64];
1693     cairo_bo_start_event_t **event_y = NULL;
1694     int i, num_events, y, ymin = 0, ymax = 0;
1695     cairo_status_t status;
1696
1697     num_events = polygon->num_edges;
1698     if (unlikely (0 == num_events))
1699         return CAIRO_STATUS_SUCCESS;
1700
1701     if (polygon->num_limits) {
1702         ymin = _cairo_fixed_integer_floor (polygon->limit.p1.y);
1703         ymax = _cairo_fixed_integer_ceil (polygon->limit.p2.y) - ymin;
1704
1705         if (ymax > 64) {
1706             event_y = _cairo_malloc_ab(sizeof (cairo_bo_event_t*), ymax);
1707             if (unlikely (event_y == NULL))
1708                 return _cairo_error (CAIRO_STATUS_NO_MEMORY);
1709         } else {
1710             event_y = stack_event_y;
1711         }
1712         memset (event_y, 0, ymax * sizeof(cairo_bo_event_t *));
1713     }
1714
1715     events = stack_events;
1716     event_ptrs = stack_event_ptrs;
1717     if (num_events > ARRAY_LENGTH (stack_events)) {
1718         events = _cairo_malloc_ab_plus_c (num_events,
1719                                           sizeof (cairo_bo_start_event_t) +
1720                                           sizeof (cairo_bo_event_t *),
1721                                           sizeof (cairo_bo_event_t *));
1722         if (unlikely (events == NULL)) {
1723             if (event_y != stack_event_y)
1724                 free (event_y);
1725             return _cairo_error (CAIRO_STATUS_NO_MEMORY);
1726         }
1727
1728         event_ptrs = (cairo_bo_event_t **) (events + num_events);
1729     }
1730
1731     for (i = 0; i < num_events; i++) {
1732         events[i].type = CAIRO_BO_EVENT_TYPE_START;
1733         events[i].point.y = polygon->edges[i].top;
1734         events[i].point.x =
1735             _line_compute_intersection_x_for_y (&polygon->edges[i].line,
1736                                                 events[i].point.y);
1737
1738         events[i].edge.edge = polygon->edges[i];
1739         events[i].edge.deferred_trap.right = NULL;
1740         events[i].edge.prev = NULL;
1741         events[i].edge.next = NULL;
1742         events[i].edge.colinear = NULL;
1743
1744         if (event_y) {
1745             y = _cairo_fixed_integer_floor (events[i].point.y) - ymin;
1746             events[i].edge.next = (cairo_bo_edge_t *) event_y[y];
1747             event_y[y] = (cairo_bo_start_event_t *) &events[i];
1748         } else
1749             event_ptrs[i] = (cairo_bo_event_t *) &events[i];
1750     }
1751
1752     if (event_y) {
1753         for (y = i = 0; y < ymax && i < num_events; y++) {
1754             cairo_bo_start_event_t *e;
1755             int j = i;
1756             for (e = event_y[y]; e; e = (cairo_bo_start_event_t *)e->edge.next)
1757                 event_ptrs[i++] = (cairo_bo_event_t *) e;
1758             if (i > j + 1)
1759                 _cairo_bo_event_queue_sort (event_ptrs+j, i-j);
1760         }
1761         if (event_y != stack_event_y)
1762             free (event_y);
1763     } else
1764         _cairo_bo_event_queue_sort (event_ptrs, i);
1765     event_ptrs[i] = NULL;
1766
1767 #if DEBUG_TRAPS
1768     dump_edges (events, num_events, "bo-polygon-edges.txt");
1769 #endif
1770
1771     /* XXX: This would be the convenient place to throw in multiple
1772      * passes of the Bentley-Ottmann algorithm. It would merely
1773      * require storing the results of each pass into a temporary
1774      * cairo_traps_t. */
1775     status = _cairo_bentley_ottmann_tessellate_bo_edges (event_ptrs, num_events,
1776                                                          fill_rule, traps,
1777                                                          &intersections);
1778 #if DEBUG_TRAPS
1779     dump_traps (traps, "bo-polygon-out.txt");
1780 #endif
1781
1782     if (events != stack_events)
1783         free (events);
1784
1785     return status;
1786 }
1787
1788 cairo_status_t
1789 _cairo_bentley_ottmann_tessellate_traps (cairo_traps_t *traps,
1790                                          cairo_fill_rule_t fill_rule)
1791 {
1792     cairo_status_t status;
1793     cairo_polygon_t polygon;
1794     int i;
1795
1796     if (unlikely (0 == traps->num_traps))
1797         return CAIRO_STATUS_SUCCESS;
1798
1799 #if DEBUG_TRAPS
1800     dump_traps (traps, "bo-traps-in.txt");
1801 #endif
1802
1803     _cairo_polygon_init (&polygon, traps->limits, traps->num_limits);
1804
1805     for (i = 0; i < traps->num_traps; i++) {
1806         status = _cairo_polygon_add_line (&polygon,
1807                                           &traps->traps[i].left,
1808                                           traps->traps[i].top,
1809                                           traps->traps[i].bottom,
1810                                           1);
1811         if (unlikely (status))
1812             goto CLEANUP;
1813
1814         status = _cairo_polygon_add_line (&polygon,
1815                                           &traps->traps[i].right,
1816                                           traps->traps[i].top,
1817                                           traps->traps[i].bottom,
1818                                           -1);
1819         if (unlikely (status))
1820             goto CLEANUP;
1821     }
1822
1823     _cairo_traps_clear (traps);
1824     status = _cairo_bentley_ottmann_tessellate_polygon (traps,
1825                                                         &polygon,
1826                                                         fill_rule);
1827
1828 #if DEBUG_TRAPS
1829     dump_traps (traps, "bo-traps-out.txt");
1830 #endif
1831
1832   CLEANUP:
1833     _cairo_polygon_fini (&polygon);
1834
1835     return status;
1836 }
1837
1838 #if 0
1839 static cairo_bool_t
1840 edges_have_an_intersection_quadratic (cairo_bo_edge_t   *edges,
1841                                       int                num_edges)
1842
1843 {
1844     int i, j;
1845     cairo_bo_edge_t *a, *b;
1846     cairo_bo_point32_t intersection;
1847
1848     /* We must not be given any upside-down edges. */
1849     for (i = 0; i < num_edges; i++) {
1850         assert (_cairo_bo_point32_compare (&edges[i].top, &edges[i].bottom) < 0);
1851         edges[i].line.p1.x <<= CAIRO_BO_GUARD_BITS;
1852         edges[i].line.p1.y <<= CAIRO_BO_GUARD_BITS;
1853         edges[i].line.p2.x <<= CAIRO_BO_GUARD_BITS;
1854         edges[i].line.p2.y <<= CAIRO_BO_GUARD_BITS;
1855     }
1856
1857     for (i = 0; i < num_edges; i++) {
1858         for (j = 0; j < num_edges; j++) {
1859             if (i == j)
1860                 continue;
1861
1862             a = &edges[i];
1863             b = &edges[j];
1864
1865             if (! _cairo_bo_edge_intersect (a, b, &intersection))
1866                 continue;
1867
1868             printf ("Found intersection (%d,%d) between (%d,%d)-(%d,%d) and (%d,%d)-(%d,%d)\n",
1869                     intersection.x,
1870                     intersection.y,
1871                     a->line.p1.x, a->line.p1.y,
1872                     a->line.p2.x, a->line.p2.y,
1873                     b->line.p1.x, b->line.p1.y,
1874                     b->line.p2.x, b->line.p2.y);
1875
1876             return TRUE;
1877         }
1878     }
1879     return FALSE;
1880 }
1881
1882 #define TEST_MAX_EDGES 10
1883
1884 typedef struct test {
1885     const char *name;
1886     const char *description;
1887     int num_edges;
1888     cairo_bo_edge_t edges[TEST_MAX_EDGES];
1889 } test_t;
1890
1891 static test_t
1892 tests[] = {
1893     {
1894         "3 near misses",
1895         "3 edges all intersecting very close to each other",
1896         3,
1897         {
1898             { { 4, 2}, {0, 0}, { 9, 9}, NULL, NULL },
1899             { { 7, 2}, {0, 0}, { 2, 3}, NULL, NULL },
1900             { { 5, 2}, {0, 0}, { 1, 7}, NULL, NULL }
1901         }
1902     },
1903     {
1904         "inconsistent data",
1905         "Derived from random testing---was leading to skip list and edge list disagreeing.",
1906         2,
1907         {
1908             { { 2, 3}, {0, 0}, { 8, 9}, NULL, NULL },
1909             { { 2, 3}, {0, 0}, { 6, 7}, NULL, NULL }
1910         }
1911     },
1912     {
1913         "failed sort",
1914         "A test derived from random testing that leads to an inconsistent sort --- looks like we just can't attempt to validate the sweep line with edge_compare?",
1915         3,
1916         {
1917             { { 6, 2}, {0, 0}, { 6, 5}, NULL, NULL },
1918             { { 3, 5}, {0, 0}, { 5, 6}, NULL, NULL },
1919             { { 9, 2}, {0, 0}, { 5, 6}, NULL, NULL },
1920         }
1921     },
1922     {
1923         "minimal-intersection",
1924         "Intersection of a two from among the smallest possible edges.",
1925         2,
1926         {
1927             { { 0, 0}, {0, 0}, { 1, 1}, NULL, NULL },
1928             { { 1, 0}, {0, 0}, { 0, 1}, NULL, NULL }
1929         }
1930     },
1931     {
1932         "simple",
1933         "A simple intersection of two edges at an integer (2,2).",
1934         2,
1935         {
1936             { { 1, 1}, {0, 0}, { 3, 3}, NULL, NULL },
1937             { { 2, 1}, {0, 0}, { 2, 3}, NULL, NULL }
1938         }
1939     },
1940     {
1941         "bend-to-horizontal",
1942         "With intersection truncation one edge bends to horizontal",
1943         2,
1944         {
1945             { { 9, 1}, {0, 0}, {3, 7}, NULL, NULL },
1946             { { 3, 5}, {0, 0}, {9, 9}, NULL, NULL }
1947         }
1948     }
1949 };
1950
1951 /*
1952     {
1953         "endpoint",
1954         "An intersection that occurs at the endpoint of a segment.",
1955         {
1956             { { 4, 6}, { 5, 6}, NULL, { { NULL }} },
1957             { { 4, 5}, { 5, 7}, NULL, { { NULL }} },
1958             { { 0, 0}, { 0, 0}, NULL, { { NULL }} },
1959         }
1960     }
1961     {
1962         name = "overlapping",
1963         desc = "Parallel segments that share an endpoint, with different slopes.",
1964         edges = {
1965             { top = { x = 2, y = 0}, bottom = { x = 1, y = 1}},
1966             { top = { x = 2, y = 0}, bottom = { x = 0, y = 2}},
1967             { top = { x = 0, y = 3}, bottom = { x = 1, y = 3}},
1968             { top = { x = 0, y = 3}, bottom = { x = 2, y = 3}},
1969             { top = { x = 0, y = 4}, bottom = { x = 0, y = 6}},
1970             { top = { x = 0, y = 5}, bottom = { x = 0, y = 6}}
1971         }
1972     },
1973     {
1974         name = "hobby_stage_3",
1975         desc = "A particularly tricky part of the 3rd stage of the 'hobby' test below.",
1976         edges = {
1977             { top = { x = -1, y = -2}, bottom = { x =  4, y = 2}},
1978             { top = { x =  5, y =  3}, bottom = { x =  9, y = 5}},
1979             { top = { x =  5, y =  3}, bottom = { x =  6, y = 3}},
1980         }
1981     },
1982     {
1983         name = "hobby",
1984         desc = "Example from John Hobby's paper. Requires 3 passes of the iterative algorithm.",
1985         edges = {
1986             { top = { x =   0, y =   0}, bottom = { x =   9, y =   5}},
1987             { top = { x =   0, y =   0}, bottom = { x =  13, y =   6}},
1988             { top = { x =  -1, y =  -2}, bottom = { x =   9, y =   5}}
1989         }
1990     },
1991     {
1992         name = "slope",
1993         desc = "Edges with same start/stop points but different slopes",
1994         edges = {
1995             { top = { x = 4, y = 1}, bottom = { x = 6, y = 3}},
1996             { top = { x = 4, y = 1}, bottom = { x = 2, y = 3}},
1997             { top = { x = 2, y = 4}, bottom = { x = 4, y = 6}},
1998             { top = { x = 6, y = 4}, bottom = { x = 4, y = 6}}
1999         }
2000     },
2001     {
2002         name = "horizontal",
2003         desc = "Test of a horizontal edge",
2004         edges = {
2005             { top = { x = 1, y = 1}, bottom = { x = 6, y = 6}},
2006             { top = { x = 2, y = 3}, bottom = { x = 5, y = 3}}
2007         }
2008     },
2009     {
2010         name = "vertical",
2011         desc = "Test of a vertical edge",
2012         edges = {
2013             { top = { x = 5, y = 1}, bottom = { x = 5, y = 7}},
2014             { top = { x = 2, y = 4}, bottom = { x = 8, y = 5}}
2015         }
2016     },
2017     {
2018         name = "congruent",
2019         desc = "Two overlapping edges with the same slope",
2020         edges = {
2021             { top = { x = 5, y = 1}, bottom = { x = 5, y = 7}},
2022             { top = { x = 5, y = 2}, bottom = { x = 5, y = 6}},
2023             { top = { x = 2, y = 4}, bottom = { x = 8, y = 5}}
2024         }
2025     },
2026     {
2027         name = "multi",
2028         desc = "Several segments with a common intersection point",
2029         edges = {
2030             { top = { x = 1, y = 2}, bottom = { x = 5, y = 4} },
2031             { top = { x = 1, y = 1}, bottom = { x = 5, y = 5} },
2032             { top = { x = 2, y = 1}, bottom = { x = 4, y = 5} },
2033             { top = { x = 4, y = 1}, bottom = { x = 2, y = 5} },
2034             { top = { x = 5, y = 1}, bottom = { x = 1, y = 5} },
2035             { top = { x = 5, y = 2}, bottom = { x = 1, y = 4} }
2036         }
2037     }
2038 };
2039 */
2040
2041 static int
2042 run_test (const char            *test_name,
2043           cairo_bo_edge_t       *test_edges,
2044           int                    num_edges)
2045 {
2046     int i, intersections, passes;
2047     cairo_bo_edge_t *edges;
2048     cairo_array_t intersected_edges;
2049
2050     printf ("Testing: %s\n", test_name);
2051
2052     _cairo_array_init (&intersected_edges, sizeof (cairo_bo_edge_t));
2053
2054     intersections = _cairo_bentley_ottmann_intersect_edges (test_edges, num_edges, &intersected_edges);
2055     if (intersections)
2056         printf ("Pass 1 found %d intersections:\n", intersections);
2057
2058
2059     /* XXX: Multi-pass Bentley-Ottmmann. Preferable would be to add a
2060      * pass of Hobby's tolerance-square algorithm instead. */
2061     passes = 1;
2062     while (intersections) {
2063         int num_edges = _cairo_array_num_elements (&intersected_edges);
2064         passes++;
2065         edges = _cairo_malloc_ab (num_edges, sizeof (cairo_bo_edge_t));
2066         assert (edges != NULL);
2067         memcpy (edges, _cairo_array_index (&intersected_edges, 0), num_edges * sizeof (cairo_bo_edge_t));
2068         _cairo_array_fini (&intersected_edges);
2069         _cairo_array_init (&intersected_edges, sizeof (cairo_bo_edge_t));
2070         intersections = _cairo_bentley_ottmann_intersect_edges (edges, num_edges, &intersected_edges);
2071         free (edges);
2072
2073         if (intersections){
2074             printf ("Pass %d found %d remaining intersections:\n", passes, intersections);
2075         } else {
2076             if (passes > 3)
2077                 for (i = 0; i < passes; i++)
2078                     printf ("*");
2079             printf ("No remainining intersections found after pass %d\n", passes);
2080         }
2081     }
2082
2083     if (edges_have_an_intersection_quadratic (_cairo_array_index (&intersected_edges, 0),
2084                                               _cairo_array_num_elements (&intersected_edges)))
2085         printf ("*** FAIL ***\n");
2086     else
2087         printf ("PASS\n");
2088
2089     _cairo_array_fini (&intersected_edges);
2090
2091     return 0;
2092 }
2093
2094 #define MAX_RANDOM 300
2095
2096 int
2097 main (void)
2098 {
2099     char random_name[] = "random-XX";
2100     cairo_bo_edge_t random_edges[MAX_RANDOM], *edge;
2101     unsigned int i, num_random;
2102     test_t *test;
2103
2104     for (i = 0; i < ARRAY_LENGTH (tests); i++) {
2105         test = &tests[i];
2106         run_test (test->name, test->edges, test->num_edges);
2107     }
2108
2109     for (num_random = 0; num_random < MAX_RANDOM; num_random++) {
2110         srand (0);
2111         for (i = 0; i < num_random; i++) {
2112             do {
2113                 edge = &random_edges[i];
2114                 edge->line.p1.x = (int32_t) (10.0 * (rand() / (RAND_MAX + 1.0)));
2115                 edge->line.p1.y = (int32_t) (10.0 * (rand() / (RAND_MAX + 1.0)));
2116                 edge->line.p2.x = (int32_t) (10.0 * (rand() / (RAND_MAX + 1.0)));
2117                 edge->line.p2.y = (int32_t) (10.0 * (rand() / (RAND_MAX + 1.0)));
2118                 if (edge->line.p1.y > edge->line.p2.y) {
2119                     int32_t tmp = edge->line.p1.y;
2120                     edge->line.p1.y = edge->line.p2.y;
2121                     edge->line.p2.y = tmp;
2122                 }
2123             } while (edge->line.p1.y == edge->line.p2.y);
2124         }
2125
2126         sprintf (random_name, "random-%02d", num_random);
2127
2128         run_test (random_name, random_edges, num_random);
2129     }
2130
2131     return 0;
2132 }
2133 #endif