lib/gfxpoly/poly.c

   1 #include <stdlib.h>
   2 #include <assert.h>
   3 #include <memory.h>
   4 #include <math.h>
   5 #include "../q.h"
   6 #include "poly.h"
   7 #include "active.h"
   8 #include "xrow.h"
   9
  10 #define DEBUG
  11
  12 char point_equals(const void*o1, const void*o2)
  13 {
  14     const point_t*p1 = o1;
  15     const point_t*p2 = o2;
  16     return p1->x == p2->x && p1->y == p2->y;
  17 }
  18 unsigned int point_hash(const void*o)
  19 {
  20     const point_t*p = o;
  21     return p->x^p->y;
  22 }
  23 void* point_dup(const void*o)
  24 {
  25     const point_t*p = o;
  26     point_t*n = malloc(sizeof(point_t));
  27     n->x = p->x;
  28     n->y = p->y;
  29     return n;
  30 }
  31 void point_free(void*o)
  32 {
  33     point_t*p = o;
  34     p->x = 0;
  35     p->y = 0;
  36     free(p);
  37 }
  38 type_t point_type = {
  39     equals: point_equals,
  40     hash: point_hash,
  41     dup: point_dup,
  42     free: point_free,
  43 };
  44
  45 typedef struct _status {
  46     int y;
  47     actlist_t*actlist;
  48     heap_t*queue;
  49 #ifdef DEBUG
  50     dict_t*seen_crossings; //list of crossing we saw so far
  51     dict_t*intersecting_segs; //list of segments intersecting in this scanline
  52     dict_t*segs_with_point; //lists of segments that received a point in this scanline
  53 #endif
  54 } status_t;
  55
  56 int compare_events(const void*_a,const void*_b)
  57 {
  58     event_t* a = (event_t*)_a;
  59     event_t* b = (event_t*)_b;
  60     if(a->p.y < b->p.y) {
  61         return 1;
  62     } else if(a->p.y > b->p.y) {
  63         return -1;
  64     } else if(a->type < b->type) {
  65         return 1;
  66     } else if(a->type > b->type) {
  67         return -1;
  68     } else if(a->p.x < b->p.x) {
  69         return 1;
  70     } else if(a->p.x > b->p.x) {
  71         return -1;
  72     } else
  73         return 0;
  74 }
  75
  76 gfxpoly_t* gfxpoly_new()
  77 {
  78     return 0;
  79 }
  80 void gfxpoly_destroy(gfxpoly_t*poly)
  81 {
  82     segment_t* s = poly;
  83     while(s) {
  84         segment_t*right = s->right;
  85         free(s);
  86         s = right;
  87         if(s==poly)
  88             break;
  89     }
  90 }
  91
  92 void gfxpoly_dump(gfxpoly_t*poly)
  93 {
  94     segment_t* s = (segment_t*)poly;
  95     while(s) {
  96         printf("[%d] (%d,%d) -> (%d,%d) %s\n",
  97                 s->nr, s->a.x, s->a.y, s->b.x, s->b.y,
  98                 s->dir==DIR_UP?"up":"down");
  99         s = s->right;
 100         if(s==poly)
 101             break;
 102     }
 103 }
 104
 105 inline static event_t event_new()
 106 {
 107     event_t e;
 108     memset(&e, 0, sizeof(e));
 109     return e;
 110 }
 111
 112 void event_dump(event_t*e)
 113 {
 114     if(e->type == EVENT_START) {
 115         printf("event: segment [%d] starts at (%d,%d)\n", e->s1->nr, e->p.x, e->p.y);
 116     } else if(e->type == EVENT_END) {
 117         printf("event: segment [%d] ends at (%d,%d)\n", e->s1->nr, e->p.x, e->p.y);
 118     } else if(e->type == EVENT_CROSS) {
 119         printf("event: segment [%d] and [%d] intersect at (%d,%d)\n", e->s1->nr, e->s2->nr, e->p.x, e->p.y);
 120     } else assert(0);
 121 }
 122
 123 static inline max32(int32_t v1, int32_t v2) {return v1>v2?v1:v2;}
 124 static inline min32(int32_t v1, int32_t v2) {return v1<v2?v1:v2;}
 125
 126 void segment_init(segment_t*s, int x1, int y1, int x2, int y2)
 127 {
 128     assert(y1!=y2);
 129     if(y1<y2) {
 130         s->dir = DIR_DOWN;
 131     } else {
 132         int x = x1;x1=x2;x2=x;
 133         int y = y1;y1=y2;y2=y;
 134         s->dir = DIR_UP;
 135     }
 136     s->a.x = x1;
 137     s->a.y = y1;
 138     s->b.x = x2;
 139     s->b.y = y2;
 140     s->k = (double)x1*y2-(double)x2*y1;
 141     s->left = s->right = 0;
 142     s->delta.x = x2-x1;
 143     s->delta.y = y2-y1;
 144     s->pos = s->a;
 145     s->tmp = -1;
 146 #ifdef DEBUG
 147     static int segment_count=0; //for debugging
 148     s->nr = segment_count++;
 149 #endif
 150
 151     assert(LINE_EQ(s->a, s) == 0);
 152     assert(LINE_EQ(s->b, s) == 0);
 153
 154     /* check that all signs are in order:
 155        a        a
 156        |\      /|
 157        | \    / |
 158      minx-b  b--maxx
 159      < 0        > 0
 160     */
 161     point_t p = s->b;
 162     p.x = min32(s->a.x, s->b.x);
 163     assert(LINE_EQ(p, s) <= 0);
 164     p.x = max32(s->a.x, s->b.x);
 165     assert(LINE_EQ(p, s) >= 0);
 166
 167     dict_init2(&s->scheduled_crossings, &ptr_type, 0);
 168 }
 169
 170 segment_t*segment_new(int x1, int y1, int x2, int y2)
 171 {
 172     segment_t*s = malloc(sizeof(segment_t));
 173     segment_init(s, x1, y1, x2, y2);
 174     return s;
 175 }
 176
 177 void segment_insert_before(segment_t**list, segment_t*add)
 178 {
 179     if(!*list) {
 180         *list = add;
 181         add->left = add->right = add;
 182     } else {
 183         add->left = (*list)->left;
 184         add->right = (*list);
 185         add->right->left = add;
 186         add->left->right = add;
 187     }
 188 }
 189
 190 void segments_enqueue(segment_t*list, heap_t*queue)
 191 {
 192     segment_t*l = list;
 193     while(l) {
 194         segment_t*right = l->right;
 195         l->left = l->right = 0;
 196         event_t e = event_new();
 197         e.type = EVENT_START;
 198         e.p = l->a;e.s1 = l;e.s2 = 0;
 199         heap_put(queue, &e);
 200         e.type = EVENT_END;
 201         e.p = l->b;e.s1 = l;e.s2 = 0;
 202         heap_put(queue, &e);
 203         l = right;
 204         if(l==list) break;
 205     }
 206 }
 207
 208 void schedule_crossing(status_t*status, segment_t*s1, segment_t*s2, int miny)
 209 {
 210     /* the code that's required (and the checks you can perform) before
 211        it can be said with 100% certainty that we indeed have a valid crossing
 212        amazes me every time. -mk */
 213     assert(s1!=s2);
 214
 215     if(dict_contains(&s1->scheduled_crossings, s2)) {
 216         // we already know about this one
 217         return;
 218     }
 219
 220     /* we probably could precompute these */
 221     int32_t minx1 = min32(s1->a.x,s1->b.x);
 222     int32_t miny1 = min32(s1->a.y,s1->b.y);
 223     int32_t maxx1 = max32(s1->a.x,s1->b.x);
 224     int32_t maxy1 = max32(s1->a.y,s1->b.y);
 225     int32_t minx2 = min32(s2->a.x,s2->b.x);
 226     int32_t miny2 = min32(s2->a.y,s2->b.y);
 227     int32_t maxx2 = max32(s2->a.x,s2->b.x);
 228     int32_t maxy2 = max32(s2->a.y,s2->b.y);
 229
 230     if(maxx1 <= minx2 || maxx2 <= minx1 ||
 231        maxy1 <= miny2 || maxy2 <= miny1) {
 232         /* bounding boxes don't intersect */
 233         return;
 234     }
 235     double adx = s1->delta.x;
 236     double ady = s1->delta.y;
 237     double bdx = s2->delta.x;
 238     double bdy = s2->delta.y;
 239     double det = adx*bdy - ady*bdx;
 240     if(!det) {
 241         if(s1->k == s2->k) {
 242             // lines are exactly on top of each other (ignored)
 243             fprintf(stderr, "Notice: segments [%d] and [%d] are exactly on top of each other\n", s1->nr, s2->nr);
 244             return;
 245         } else {
 246             /* lines are parallel */
 247             return;
 248         }
 249     }
 250     double asign2 = LINE_EQ(s1->a, s2);
 251     double bsign2 = LINE_EQ(s1->b, s2);
 252     if(asign2<0 && bsign2<0) {
 253         // segment1 is completely to the left of segment2
 254         return;
 255     }
 256     if(asign2>0 && bsign2>0)  {
 257         // segment2 is completely to the left of segment1
 258         return;
 259     }
 260     if(asign2==0) {
 261         // segment1 touches segment2 in a single point (ignored)
 262 #ifdef DEBUG
 263         fprintf(stderr, "Notice: segment [%d]'s start point touches segment [%d]\n", s1->nr, s2->nr);
 264 #endif
 265         return;
 266     }
 267     if(bsign2==0) {
 268         // segment1 touches segment2 in a single point (ignored)
 269 #ifdef DEBUG
 270         fprintf(stderr, "Notice: segment [%d]'s end point touches segment [%d]\n", s1->nr, s2->nr);
 271 #endif
 272         return;
 273     }
 274     double asign1 = LINE_EQ(s2->a, s1);
 275     double bsign1 = LINE_EQ(s2->b, s1);
 276     if(asign1<0 && bsign1<0) {
 277         // segment1 is completely to the left of segment2
 278         return;
 279     }
 280     if(asign1>0 && bsign1>0)  {
 281         // segment2 is completely to the left of segment1
 282         return;
 283     }
 284     if(asign1==0) {
 285         // segment2 touches segment1 in a single point (ignored)
 286 #ifdef DEBUG
 287         fprintf(stderr, "Notice: segment [%d]'s start point touches segment [%d]\n", s2->nr, s1->nr);
 288 #endif
 289         return;
 290     }
 291     if(asign2==0) {
 292         // segment2 touches segment1 in a single point (ignored)
 293 #ifdef DEBUG
 294         fprintf(stderr, "Notice: segment [%d]'s end point touches segment [%d]\n", s2->nr, s1->nr);
 295 #endif
 296         return;
 297     }
 298
 299     double la = (double)s1->a.x*(double)s1->b.y - (double)s1->a.y*(double)s1->b.x;
 300     double lb = (double)s2->a.x*(double)s2->b.y - (double)s2->a.y*(double)s2->b.x;
 301
 302     point_t p;
 303     p.x = (int32_t)ceil((-la*bdx +lb*adx) / det);
 304     p.y = (int32_t)ceil((+lb*ady -la*bdy) / det);
 305     /* we only schedule crossings if they are *after* the current y. That way, we don't
 306        schedule the same crossing twice */
 307     if(p.y >= miny) {
 308 #ifdef DEBUG
 309         point_t pair;
 310         pair.x = s1->nr;
 311         pair.y = s2->nr;
 312         assert(!dict_contains(status->seen_crossings, &pair));
 313         dict_put(status->seen_crossings, &pair, 0);
 314 #endif
 315         event_t e = event_new();
 316         e.type = EVENT_CROSS;
 317         e.p = p;
 318         e.s1 = s1;
 319         e.s2 = s2;
 320 #ifdef DEBUG
 321         printf("schedule crossing between [%d] and [%d] at (%d,%d)\n", s1->nr, s2->nr, e.p.x, e.p.y);
 322 #endif
 323         /* we insert into each other's intersection history because these segments might switch
 324            places */
 325         dict_put(&s1->scheduled_crossings, s2, 0);
 326         dict_put(&s2->scheduled_crossings, s1, 0);
 327
 328         heap_put(status->queue, &e);
 329     }
 330     return;
 331 }
 332
 333 /* routine dealing with the special case of a number of line segments intersecting in exactly one
 334    point */
 335 void exchange_many(status_t*status, event_t*e)
 336 {
 337     int size = 16;
 338     int n = 0;
 339     segment_t**segs = (segment_t**)malloc(sizeof(segment_t*)*size);
 340     segs[n] = e->s1;e->s1->tmp = n;n++;
 341     segs[n] = e->s2;e->s2->tmp = n;n++;
 342 #ifdef DEBUG
 343     dict_put(status->intersecting_segs, e->s1, 0);
 344     dict_put(status->intersecting_segs, e->s2, 0);
 345 #endif
 346     while(1) {
 347         event_t*ee = heap_peek(status->queue);
 348         if(ee->p.x != e->p.x ||
 349            ee->p.y != e->p.y ||
 350            ee->type != e->type) {
 351             break;
 352         }
 353         if(n>=size-1) {
 354             size *= 2;
 355             segs = realloc(segs, sizeof(segment_t*)*size);
 356         }
 357         ee = heap_chopmax(status->queue);
 358         if(ee->s1->tmp<0) {
 359             segs[n] = ee->s1;ee->s1->tmp = n;n++;
 360         }
 361         if(ee->s2->tmp<0) {
 362             segs[n] = ee->s2;ee->s2->tmp = n;n++;
 363         }
 364 #ifdef DEBUG
 365         dict_put(status->intersecting_segs, ee->s1, 0);
 366         dict_put(status->intersecting_segs, ee->s2, 0);
 367 #endif
 368     }
 369     int t;
 370 #ifdef DEBUG
 371     printf("event: multi-intersect: ");
 372     for(t=0;t<n;t++) {
 373         printf("[%d] ", segs[t]->nr);
 374     }
 375     printf("\n");
 376 #endif
 377     segment_t*left = 0;
 378     segment_t*right = 0;
 379     char warn = 0;
 380     for(t=0;t<n;t++) {
 381         if(!segs[t]->left || segs[t]->left->tmp<0) {
 382             assert(!left);
 383             left = segs[t];
 384         }
 385         if(!segs[t]->right || segs[t]->right->tmp<0) {
 386             assert(!right);
 387             right = segs[t];
 388         }
 389         if(segs[t]->tmp<0)
 390             warn = 1;
 391     }
 392     if(warn)
 393         fprintf(stderr, "Warning: multi-cross section contains rogue segments\n");
 394
 395     assert(left);
 396     assert(right);
 397     assert(left!=right);
 398     actlist_invert_fromto(status->actlist, left, right);
 399     for(t=0;t<n;t++) {
 400         segs[t]->tmp = -1;
 401     }
 402     free(segs);
 403     if(right->left)
 404         schedule_crossing(status, right->left, right, status->y);
 405     if(left->right)
 406         schedule_crossing(status, left, left->right, status->y);
 407 }
 408
 409 void exchange_two(status_t*status, event_t*e)
 410 {
 411     //exchange two segments in list
 412     segment_t*s1 = e->s1;
 413     segment_t*s2 = e->s2;
 414 #ifdef DEBUG
 415     dict_put(status->intersecting_segs, s1, 0);
 416     dict_put(status->intersecting_segs, s2, 0);
 417 #endif
 418     segment_t*left = actlist_left(status->actlist, s2);
 419     segment_t*right = actlist_right(status->actlist, s1);
 420     assert(left == s1);
 421     assert(right == s2);
 422     actlist_swap(status->actlist, s1, s2);
 423     assert(actlist_right(status->actlist, s2) == s1);
 424     assert(actlist_left(status->actlist, s1) == s2);
 425     left = actlist_left(status->actlist, s2);
 426     right = actlist_right(status->actlist, s1);
 427     if(left)
 428         schedule_crossing(status, left, s2, status->y);
 429     if(right)
 430         schedule_crossing(status, s1, right, status->y);
 431 }
 432
 433 void insert_point_into_segment(segment_t*s, point_t p)
 434 {
 435     if(s->pos.x == p.x && s->pos.y == p.y) {
 436 #ifdef DEBUG
 437         fprintf(stderr, "Notice: tried to add (%d,%d) to segment [%d] twice\n", p.x, p.y, s->nr);
 438 #endif
 439         return;
 440     }
 441 #ifdef DEBUG
 442     printf("[%d] gets extra point (%d,%d)\n", s->nr, p.x, p.y);
 443 #endif
 444     s->pos = p;
 445 }
 446
 447 typedef struct _box {
 448     point_t left1, left2, right1, right2;
 449 } box_t;
 450 static inline box_t box_new(int x, int y)
 451 {
 452     box_t box;
 453     box.right1.x = box.right2.x = x;
 454     box.left1.x = box.left2.x = x-1;
 455     box.left1.y = box.right1.y = y-1;
 456     box.left2.y = box.right2.y = y;
 457     return box;
 458 }
 459
 460 /* possible optimizations:
 461    1.) keep two different active lists around, one for negative and one for
 462        positive slopes
 463    2.) delay starting events until after this function (tricky, because we *do*
 464        need the start coordinates)
 465 */
 466 /*
 467    SLOPE_POSITIVE:
 468        +     \ +
 469 ------ I      \I
 470       -I\----  I
 471        I \   --I\---
 472        I  \    I \  -------
 473        +   \   +  \
 474 */
 475 void add_points_to_positively_sloped_segments(status_t*status, xrow_t*xrow, int32_t y)
 476 {
 477     int t;
 478     for(t=0;t<xrow->num;t++) {
 479         box_t box = box_new(xrow->x[t], y);
 480         segment_t*seg = actlist_find(status->actlist, box.left2, box.left2);
 481         if(seg)
 482             seg = seg->right;
 483         else
 484             seg = actlist_leftmost(status->actlist);
 485         while(seg) {
 486             if(seg->a.y == y) {
 487                 // this segment just started, ignore it
 488             } else if(seg->delta.x < 0) {
 489                 // ignore segment w/ negative slope
 490             } else {
 491                 double d1 = LINE_EQ(box.right1, seg);
 492                 double d2 = LINE_EQ(box.right2, seg);
 493                 if(d1>=0 || d2>=0) {
 494 #ifdef DEBUG
 495                     dict_put(status->segs_with_point, seg, 0);
 496 #endif
 497                     insert_point_into_segment(seg, box.right2);
 498                 } else {
 499                     break;
 500                 }
 501             }
 502             seg = seg->right;
 503         }
 504     }
 505 }
 506 /* SLOPE_NEGATIVE:
 507    |   +   /|  +  /    /
 508    |   I  / |  I /    /
 509    |   I /  |  I/    /
 510    |   I/   |  I    /
 511    |   I    | /I   /
 512    |  /+    |/ +  /
 513 */
 514 void add_points_to_negatively_sloped_segments(status_t*status, xrow_t*xrow, int32_t y)
 515 {
 516     int t;
 517     for(t=xrow->num-1;t>=0;t--) {
 518         box_t box = box_new(xrow->x[t], y);
 519         segment_t*seg = actlist_find(status->actlist, box.right2, box.right2);
 520         while(seg) {
 521             if(seg->a.y == y) {
 522                 // this segment just started, ignore it
 523             } else if(seg->delta.x >= 0) {
 524                 // ignore segment w/ positive slope
 525             } else {
 526                 double d1 = LINE_EQ(box.left1, seg);
 527                 double d2 = LINE_EQ(box.left2, seg);
 528                 if(d1<0 || d2<0) {
 529 #ifdef DEBUG
 530                     dict_put(status->segs_with_point, seg, 0);
 531 #endif
 532                     insert_point_into_segment(seg, box.right2);
 533                 } else {
 534                     break;
 535                 }
 536             }
 537             seg = seg->left;
 538         }
 539     }
 540 }
 541
 542 void event_apply(status_t*status, event_t*e)
 543 {
 544     switch(e->type) {
 545         case EVENT_END: {
 546             //delete segment from list
 547             event_dump(e);
 548             segment_t*s = e->s1;
 549             insert_point_into_segment(s, s->b);
 550             segment_t*left = actlist_left(status->actlist, s);
 551             segment_t*right = actlist_right(status->actlist, s);
 552             actlist_delete(status->actlist, s);
 553             dict_clear(&s->scheduled_crossings);
 554             if(left && right)
 555                 schedule_crossing(status, left, right, status->y+1);
 556             break;
 557         }
 558         case EVENT_START: {
 559             //insert segment into list
 560             event_dump(e);
 561             segment_t*s = e->s1;
 562             actlist_insert(status->actlist, e->p, s);
 563             segment_t*left = actlist_left(status->actlist, s);
 564             segment_t*right = actlist_right(status->actlist, s);
 565             if(left)
 566                 schedule_crossing(status, left, s, status->y+1);
 567             if(right)
 568                 schedule_crossing(status, s, right, status->y+1);
 569             break;
 570         }
 571         case EVENT_CROSS: {
 572             // exchange two (or more) segments
 573             event_t*ee = heap_peek(status->queue);
 574             if(ee->p.x != e->p.x ||
 575                ee->p.y != e->p.y ||
 576                ee->type != e->type) {
 577                 event_dump(e);
 578                 exchange_two(status, e);
 579             } else {
 580                 exchange_many(status, e);
 581             }
 582         }
 583     }
 584 }
 585
 586 #ifdef DEBUG
 587 void check_status(status_t*status)
 588 {
 589     DICT_ITERATE_KEY(status->intersecting_segs, segment_t*, s) {
 590         if(!dict_contains(status->segs_with_point, s)) {
 591             fprintf(stderr, "Error: segment [%d] (%sslope) intersects in scanline %d, but it didn't receive a point\n",
 592                     s->nr,
 593                     s->delta.x<0?"-":"+",
 594                     status->y);
 595             assert(0);
 596         }
 597     }
 598 }
 599 #endif
 600
 601 void gfxpoly_process(segment_t*poly)
 602 {
 603     heap_t* queue = heap_new(sizeof(event_t), compare_events);
 604     segments_enqueue(poly, queue);
 605     status_t status;
 606     memset(&status, 0, sizeof(status_t));
 607     status.queue = queue;
 608     status.actlist = actlist_new();
 609 #ifdef DEBUG
 610     status.seen_crossings = dict_new2(&point_type);
 611 #endif
 612
 613     xrow_t*xrow = xrow_new();
 614
 615     event_t*e = heap_chopmax(queue);
 616     while(e) {
 617         status.y = e->p.y;
 618 #ifdef DEBUG
 619         status.intersecting_segs = dict_new2(&ptr_type);
 620         status.segs_with_point = dict_new2(&ptr_type);
 621         printf("----------------------------------- %d\n", status.y);
 622         actlist_verify_and_dump(status.actlist, status.y-1);
 623 #endif
 624         xrow_reset(xrow);
 625         do {
 626             xrow_add(xrow, e->p.x);
 627             event_apply(&status, e);
 628             free(e);
 629             e = heap_chopmax(queue);
 630         } while(e && status.y == e->p.y);
 631
 632         xrow_sort(xrow);
 633         add_points_to_positively_sloped_segments(&status, xrow, status.y);
 634         add_points_to_negatively_sloped_segments(&status, xrow, status.y);
 635 #ifdef DEBUG
 636         check_status(&status);
 637         dict_destroy(status.intersecting_segs);
 638         dict_destroy(status.segs_with_point);
 639 #endif
 640     }
 641 #ifdef DEBUG
 642     dict_destroy(status.seen_crossings);
 643 #endif
 644     heap_destroy(queue);
 645     gfxpoly_destroy(poly);
 646     xrow_destroy(xrow);
 647 }