Bug fixes for the polygon rendering code where naming the same pixel twice
[imager.git] / polygon.c
1 #include "image.h"
2 #include "draw.h"
3 #include "log.h"
4
5
6 #define IMTRUNC(x) ((int)((x)*16))
7
8 #define coarse(x) ((x)/16)
9 #define fine(x)   ((x)%16)
10
11 #define POLY_DEB(x)
12
13
14
15 typedef int pcord;
16
17 typedef struct {
18   int n;
19   pcord x,y;
20 } p_point;
21
22 typedef struct {
23   int n;
24   pcord x1,y1;
25   pcord x2,y2;
26   pcord miny,maxy;
27   pcord minx,maxx;
28   int updown; /* -1 means down, 0 vertical, 1 up */
29 } p_line;
30
31 typedef struct {
32   int n;
33   double x;
34 } p_slice;
35
36 typedef struct {
37   int start;
38   int stop;
39 } ss_pair;
40
41 typedef struct {
42   int *line;            /* temporary buffer for scanline */
43   int linelen;          /* length of scanline */
44   ss_pair *ss_list;     /* list of start stop linepairs */
45   int ssnext;           /* index of the next pair to use */
46   int sslen;            /* maximum number of start stop pairs */
47 } ss_scanline;
48
49
50
51
52
53
54
55
56 static
57 int
58 p_compy(const p_point *p1, const p_point *p2) {
59   if (p1->y > p2->y) return 1;
60   if (p1->y < p2->y) return -1;
61   return 0;
62 }
63
64 static
65 int
66 p_compx(const p_slice *p1, const p_slice *p2) {
67   if (p1->x > p2->x) return 1;
68   if (p1->x < p2->x) return -1;
69   return 0;
70 }
71
72 /* Change this to int? and round right goddamn it! */
73
74 static
75 double
76 p_eval_aty(p_line *l, pcord y) {
77   int t;
78   t=l->y2-l->y1;
79   if (t) return ( (y-l->y1)*l->x2 + (l->y2-y)*l->x1 )/t;
80   return (l->x1+l->x2)/2.0;
81 }
82
83 static
84 double
85 p_eval_atx(p_line *l, pcord x) {
86   int t;
87   t = l->x2-l->x1;
88   if (t) return ( (x-l->x1)*l->y2 + (l->x2-x)*l->y1 )/t;
89   return (l->y1+l->y2)/2.0;
90 }
91
92 static
93 p_line *
94 line_set_new(double *x, double *y, int l) {
95   int i;
96   p_line *lset = mymalloc(sizeof(p_line) * l);
97
98   for(i=0; i<l; i++) {
99     lset[i].n=i;
100     lset[i].x1 = IMTRUNC(x[i]);
101     lset[i].y1 = IMTRUNC(y[i]);
102     lset[i].x2 = IMTRUNC(x[(i+1)%l]);
103     lset[i].y2 = IMTRUNC(y[(i+1)%l]);
104     lset[i].miny=min(lset[i].y1,lset[i].y2);
105     lset[i].maxy=max(lset[i].y1,lset[i].y2);
106     lset[i].minx=min(lset[i].x1,lset[i].x2);
107     lset[i].maxx=max(lset[i].x1,lset[i].x2);
108   }
109   return lset;
110 }
111
112 static
113 p_point *
114 point_set_new(double *x, double *y, int l) {
115   int i;
116   p_point *pset = mymalloc(sizeof(p_point) * l);
117   
118   for(i=0; i<l; i++) {
119     pset[i].n=i;
120     pset[i].x=IMTRUNC(x[i]);
121     pset[i].y=IMTRUNC(y[i]);
122   }
123   return pset;
124 }
125
126 static
127 void
128 p_line_dump(p_line *l) {
129   printf("%d (%d,%d)->(%d,%d) [%d-%d,%d-%d]\n", l->n, l->x1, l->y1, l->x2, l->y2, 
130          l->minx, l->maxx, l->miny, l->maxy);
131 }
132
133
134 static
135 void
136 ss_scanline_reset(ss_scanline *ss) {
137   ss->ssnext = 0;
138   memset(ss->line, 0, sizeof(int) * ss->linelen);
139 }
140
141 static
142 void
143 ss_scanline_init(ss_scanline *ss, int linelen, int linepairs) {
144   ss->line    = mymalloc( sizeof(int) * linelen );
145   ss->linelen = linelen;
146   ss->ss_list = mymalloc( sizeof(ss_pair) * linepairs );
147   ss->sslen   = linepairs;
148   ss_scanline_reset(ss);
149 }
150
151 static
152 void
153 ss_scanline_exorcise(ss_scanline *ss) {
154   myfree(ss->line);
155   myfree(ss->ss_list);
156 }
157   
158                      
159
160
161 /* returns the number of matches */
162
163 static
164 int
165 lines_in_interval(p_line *lset, int l, p_slice *tllist, pcord minc, pcord maxc) {
166   int k;
167   int count = 0;
168   for(k=0; k<l; k++) {
169     if (lset[k].maxy > minc && lset[k].miny < maxc) {
170       if (lset[k].miny == lset[k].maxy) {
171         POLY_DEB( printf(" HORIZONTAL - skipped\n") );
172       } else {
173         tllist[count].x=p_eval_aty(&lset[k],(minc+maxc)/2.0 );
174         tllist[count].n=k;
175         count++;
176       }
177     }
178   }
179   return count;
180 }
181
182 static
183 int
184 lines_in_interval_old(p_line *lset, int l, p_slice *tllist, pcord cc) {
185   int k;
186   int count = 0;
187   for(k=0; k<l; k++) {
188     if (cc >= lset[k].miny && cc <=  lset[k].maxy) {
189       if (lset[k].miny == lset[k].maxy) {
190         POLY_DEB( printf(" HORIZONTAL - skipped\n") );
191       }
192       else {
193         tllist[count].x=p_eval_aty(&lset[k],cc);
194         tllist[count].n=k;
195         count++;
196       }
197     }
198   }
199   return count;
200 }
201
202 /* marks the up variable for all lines in a slice */
203
204 static
205 void
206 mark_updown_slices(p_line *lset, p_slice *tllist, int count) {
207   p_line *l, *r;
208   int k;
209   for(k=0; k<count; k+=2) {
210     l = lset + tllist[k].n;
211
212     if (l->y1 == l->y2) {
213       mm_log((1, "mark_updown_slices: horizontal line being marked: internal error!\n"));
214       exit(3);
215     }
216
217     l->updown = (l->x1 == l->x2) ?
218       0 :
219       (l->x1 > l->x2)
220       ? 
221       (l->y1 > l->y2) ? -1 : 1
222       : 
223       (l->y1 > l->y2) ? 1 : -1;
224
225     POLY_DEB( printf("marking left line %d as %s(%d)\n", l->n,
226                      l->updown ?  l->updown == 1 ? "up" : "down" : "vert", l->updown, l->updown)
227               );
228
229     if (k+1 >= count) {
230       mm_log((1, "Invalid polygon spec, odd number of line crossings.\n"));
231       return;
232     }
233
234     r = lset + tllist[k+1].n;
235     if (r->y1 == r->y2) {
236       mm_log((1, "mark_updown_slices: horizontal line being marked: internal error!\n"));
237       exit(3);
238     }
239
240     r->updown = (r->x1 == r->x2) ?
241       0 :
242       (r->x1 > r->x2)
243       ? 
244       (r->y1 > r->y2) ? -1 : 1
245       : 
246       (r->y1 > r->y2) ? 1 : -1;
247     
248     POLY_DEB( printf("marking right line %d as %s(%d)\n", r->n,
249                      r->updown ?  r->updown == 1 ? "up" : "down" : "vert", r->updown, r->updown)
250               );
251   }
252 }
253
254
255
256 static
257 unsigned char
258 saturate(int in) {
259   if (in>255) { return 255; }
260   else if (in>0) return in;
261   return 0;
262 }
263
264
265 /* This function must be modified later to do proper blending */
266
267 void
268 scanline_flush(i_img *im, ss_scanline *ss, int y, i_color *val) {
269   int x, ch, tv;
270   i_color t;
271   for(x=0; x<im->xsize; x++) {
272     tv = saturate(ss->line[x]);
273     i_gpix(im, x, y, &t);
274     for(ch=0; ch<im->channels; ch++) 
275       t.channel[ch] = tv/255.0 * val->channel[ch] + (1.0-tv/255.0) * t.channel[ch];
276     i_ppix(im, x, y, &t);
277   }
278 }
279
280
281
282 static
283 int
284 trap_square(pcord xlen, pcord ylen, double xl, double yl) {
285   POLY_DEB( printf("trap_square: %d %d %.2f %.2f\n", xlen, ylen, xl, yl) );
286   return xlen*ylen-(xl*yl)/2.0;
287 }
288
289
290 /* 
291    pixel_coverage calculates the 'left side' pixel coverage of a pixel that is
292    within the min/max ranges.  The shape always corresponds to a square with some
293    sort of a triangle cut from it (which can also yield a triangle).
294 */
295
296
297 static
298 int 
299 pixel_coverage(p_line *line, pcord minx, pcord maxx, pcord  miny, pcord maxy) {
300   double lycross, rycross;
301   int l, r;
302
303   double xs, ys;
304   
305   if (!line->updown) {
306     l = r = 0;
307   } else {
308     lycross = p_eval_atx(line, minx);
309     rycross = p_eval_atx(line, maxx);
310     l = lycross <= maxy && lycross >= miny; /* true if it enters through left side */
311     r = rycross <= maxy && rycross >= miny; /* true if it enters through left side */
312   }
313   POLY_DEB(
314            printf("%4s(%+d): ", line->updown ?  line->updown == 1 ? "up" : "down" : "vert", line->updown);
315            printf("(%2d,%2d) [%3d-%3d, %3d-%3d] lycross=%.2f rycross=%.2f", coarse(minx), coarse(miny), minx, maxx, miny, maxy, lycross, rycross);
316            printf("    l=%d r=%d\n", l, r)
317            );
318   
319   if (l && r) 
320     return line->updown == 1 ? 
321       (double)(maxx-minx) * (2.0*maxy-lycross-rycross)/2.0  /* up case */
322       :
323       (double)(maxx-minx) * (lycross+rycross-2*miny)/2.0;  /* down case */
324   
325   if (!l && !r) return (maxy-miny)*(maxx*2-p_eval_aty(line, miny)-p_eval_aty(line, maxy))/2.0;
326
327   if (l && !r)
328     return line->updown == 1 ?
329       trap_square(maxx-minx, maxy-miny, p_eval_aty(line, miny)-minx, p_eval_atx(line, minx)-miny) : 
330       trap_square(maxx-minx, maxy-miny, p_eval_aty(line, maxy)-minx, maxy-p_eval_atx(line, minx));
331   
332
333   if (!l && r) {
334     int r = line->updown == 1 ?
335       (maxx-p_eval_aty(line, maxy))*(maxy-p_eval_atx(line, maxx))/2.0 : 
336       (maxx-p_eval_aty(line, miny))*(p_eval_atx(line, maxx)-miny)/2.0;
337     return r;
338   }
339 }
340
341
342
343
344
345 /* 
346    handle the scanline slice in three steps 
347    
348    1.  Where only the left edge is inside a pixel
349    2a. Where both left and right edge are inside a pixel
350    2b. Where neither left or right edge are inside a pixel
351    3.  Where only the right edge is inside a pixel
352 */
353
354 static
355 void
356 render_slice_scanline(ss_scanline *ss, int y, p_line *l, p_line *r) {
357   
358   pcord miny, maxy;     /* y bounds in fine coordinates */
359   pcord lminx, lmaxx;   /* left line min/max within y bounds in fine coords */
360   pcord rminx, rmaxx;   /* right line min/max within y bounds in fine coords */
361   int cpix;             /* x-coordinate of current pixel */
362   int thin;             /* boolean for thin/thick segment */
363   int startpix;         /* temporary variable for "start of this interval" */
364   int stoppix;          /* temporary variable for "end of this interval" */
365   int step2end;         /* temporary variable to mark where step2 ends */
366
367   /* Find the y bounds of scanline_slice */
368
369   maxy = min( l->maxy, r->maxy );
370   miny = max( l->miny, r->miny );
371
372   maxy = min( maxy, (y+1)*16 );
373   miny = max( miny,  y*16 );
374
375   lminx = min( p_eval_aty(l, maxy), p_eval_aty(l, miny) );
376   lmaxx = max( p_eval_aty(l, maxy), p_eval_aty(l, miny) );
377
378   rminx = min( p_eval_aty(r, maxy), p_eval_aty(r, miny) );
379   rmaxx = max( p_eval_aty(r, maxy), p_eval_aty(r, miny) );
380
381   thin = coarse(lmaxx) >= coarse(rminx);
382
383   startpix = max( coarse(lminx), 0 );
384   stoppix  = min( coarse(rmaxx-1), ss->linelen-1 );
385   
386   for(cpix=startpix; cpix<=stoppix; cpix++) {
387     int lt = coarse(lmaxx-1) >= cpix;
388     int rt = coarse(rminx) <= cpix;
389     
390     int A, B, C;
391     
392     POLY_DEB( printf("(%d,%d) lt=%d rt=%d\n", cpix, y, lt, rt) );
393
394     A = lt ? pixel_coverage(l, cpix*16, cpix*16+16, miny, maxy) : 0;
395     B = lt ? 0 : 16*(maxy-miny);
396     C = rt ? pixel_coverage(r, cpix*16, cpix*16+16, miny, maxy) : 0;
397
398     POLY_DEB( printf("A=%d B=%d C=%d\n", A, B, C) );
399
400     ss->line[cpix] += A+B-C;
401
402   }
403   
404 }
405
406
407
408 static
409 void
410 render_slice_scanline_old(ss_scanline *ss, int y, p_line *l, p_line *r) {
411   
412   pcord miny, maxy;     /* y bounds in fine coordinates */
413   pcord lminx, lmaxx;   /* left line min/max within y bounds in fine coords */
414   pcord rminx, rmaxx;   /* right line min/max within y bounds in fine coords */
415   int cpix;             /* x-coordinate of current pixel */
416   int thin;             /* boolean for thin/thick segment */
417   int startpix;         /* temporary variable for "start of this interval" */
418   int stoppix;          /* temporary variable for "end of this interval" */
419   int step2end;         /* temporary variable to mark where step2 ends */
420
421   /* Find the y bounds of scanline_slice */
422
423   maxy = min( l->maxy, r->maxy );
424   miny = max( l->miny, r->miny );
425
426   maxy = min( maxy, (y+1)*16 );
427   miny = max( miny,  y*16 );
428
429   lminx = min( p_eval_aty(l, maxy), p_eval_aty(l, miny) );
430   lmaxx = max( p_eval_aty(l, maxy), p_eval_aty(l, miny) );
431
432   rminx = min( p_eval_aty(r, maxy), p_eval_aty(r, miny) );
433   rmaxx = max( p_eval_aty(r, maxy), p_eval_aty(r, miny) );
434
435   thin = coarse(lmaxx) >= coarse(rminx);
436
437
438   /* First step */
439   startpix = coarse(lminx);                             /* includes tricky starting pixel */
440   stoppix  = min(coarse(lmaxx), coarse(rminx) );        /* last pixel is tricky */
441   
442   /* handle start pixel */
443
444   cpix = startpix;
445   if (cpix < stoppix) {
446     ss->line[cpix] += pixel_coverage(l, cpix*16, cpix*16+16, miny, maxy);
447     printf("%2d: step1 - start pixel\n", cpix);
448   }
449   
450   for(cpix=startpix+1; cpix<stoppix; cpix++) {
451     printf("%2d: step1 pixel\n", cpix);
452     ss->line[cpix] += l->updown == 1 ? 
453       8.0 * (2*maxy-p_eval_atx(l, 16*cpix)-p_eval_atx(l, 16*cpix+16))   /* up case */
454       :
455       8.0 * (p_eval_atx(l, 16*cpix)+p_eval_atx(l, 16*cpix+16)-2*miny);  /* down case */
456   }
457   
458   
459   /* handle stop pixel */
460
461   if (thin) { /* step 2a */
462     startpix = coarse(rminx);
463     stoppix = coarse(lmaxx+15); /* one more than needed */
464     
465     for(cpix=startpix; cpix<stoppix; cpix++) {
466       printf("%2d: step2a pixel\n", cpix);
467       ss->line[cpix] += 
468         pixel_coverage(l, cpix*16, cpix*16+16, miny, maxy)
469         +(cpix*16+16-min(cpix*16+16, l->maxx))*(maxy-miny)
470         -pixel_coverage(r, cpix*16, cpix*16+16, miny, maxy);
471     }
472   } else { /* step 2b */
473     stoppix = coarse(rminx);
474     for(/* cpix already correct */; cpix<stoppix; cpix++) {
475       printf("%2d: step2b pixel\n", cpix);
476       ss->line[cpix] += 16.0*(maxy-miny);
477     }
478   }
479   
480   /* step 3 */
481
482   cpix = max(coarse(rminx), coarse(lmaxx+15));
483   stoppix = coarse(rmaxx-15);
484   
485   printf("step3 from %d to %d\n", cpix, stoppix);
486
487   for(; cpix<stoppix; cpix++) {
488     printf("%2d: step3 pixel\n", cpix);
489     ss->line[cpix] += 0+ 
490       (l->updown == 1 ?
491        8.0 * (2*maxy-p_eval_atx(r, 16*cpix)-p_eval_atx(r, 16*cpix+16))  /* up case */
492        :
493        8.0 * (p_eval_atx(r, 16*cpix)+p_eval_atx(r, 16*cpix+16)-2*miny));  /* down case */
494   }
495   
496   ss->line[cpix] += (16.0)*(maxy-miny) - pixel_coverage(r, cpix*16, cpix*16+16, miny, maxy);
497 }
498
499
500
501
502
503
504 /* Antialiasing polygon algorithm 
505    specs:
506    1. only nice polygons - no crossovers
507    2. 1/16 pixel resolution
508    3. full antialiasing ( complete spectrum of blends )
509    4. uses hardly any memory
510    5. no subsampling phase
511    
512
513    Algorithm outline:
514    1. Split into vertical intervals.
515    2. handle each interval 
516
517    For each interval we must: 
518    1. find which lines are in it
519    2. order the lines from in increasing x order.
520       since we are assuming no crossovers it is sufficent
521       to check a single point on each line.
522 */
523
524 /*
525   Definitions:
526   
527   1. Interval:  A vertical segment in which no lines cross nor end.
528   2. Scanline:  A physical line, contains 16 subpixels in the horizontal direction
529   3. Slice:     A start stop line pair.
530   
531  */
532
533
534 void
535 i_poly_aa(i_img *im, int l, double *x, double *y, i_color *val) {
536   int i ,k;                     /* Index variables */
537   int clc;                      /* Lines inside current interval */
538   pcord miny ,maxy;             /* Min and max values of the current slice in the subcord system */
539   pcord tempy;
540   int cscl;                     /* Current scanline */
541
542   ss_scanline templine;         /* scanline accumulator */
543   p_point *pset;                /* List of points in polygon */
544   p_line  *lset;                /* List of lines in polygon */
545   p_slice *tllist;              /* List of slices */
546
547   mm_log((1, "i_poly_aa(im %p, l %d, x %p, y %p, val %p)\n", im, l, x, y, val));
548
549   for(i=0; i<l; i++) {
550     mm_log((2, "(%.2f, %.2f)\n", x[i], y[i]));
551   }
552
553
554   POLY_DEB(
555            fflush(stdout);
556            setbuf(stdout, NULL);
557            );
558
559   tllist   = mymalloc(sizeof(p_slice)*l);
560   
561   ss_scanline_init(&templine, im->xsize, l);
562
563   pset     = point_set_new(x, y, l);
564   lset     = line_set_new(x, y, l);
565
566
567   qsort(pset, l, sizeof(p_point), (int(*)(const void *,const void *))p_compy);
568   
569   POLY_DEB(
570            for(i=0;i<l;i++) {
571              printf("%d [ %d ] (%d , %d) -> (%d , %d) yspan ( %d , %d )\n",
572                     i, lset[i].n, lset[i].x1, lset[i].y1, lset[i].x2, lset[i].y2, lset[i].miny, lset[i].maxy);
573            }
574            printf("MAIN LOOP\n\n");
575            );
576   
577
578   /* loop on intervals */
579   for(i=0; i<l-1; i++) {
580     int startscan = max( coarse(pset[i].y), 0);
581     int stopscan = min( coarse(pset[i+1].y+15), im->ysize);
582     pcord cc = (pset[i].y + pset[i+1].y)/2;
583
584     if (pset[i].y == pset[i+1].y) {
585       POLY_DEB( printf("current slice thickness = 0 => skipping\n") );
586       continue;
587     }
588
589     POLY_DEB(
590              printf("current slice is %d: %d to %d ( cpoint %d ) scanlines %d to %d\n", 
591                     i, pset[i].y, pset[i+1].y, cc, startscan, stopscan)
592              );
593     
594     
595     clc = lines_in_interval(lset, l, tllist, pset[i].y, pset[i+1].y);
596     qsort(tllist, clc, sizeof(p_slice), (int(*)(const void *,const void *))p_compx);
597
598     mark_updown_slices(lset, tllist, clc);
599
600     POLY_DEB( printf("Interval contains %d lines\n", clc) );
601
602     for(k=0; k<clc; k++) {
603       int lno = tllist[k].n;
604       p_line *ln = lset+lno;
605       POLY_DEB(
606                printf("%d:  line #%2d: (%2d, %2d)->(%2d, %2d) (%2d/%2d, %2d/%2d) -> (%2d/%2d, %2d/%2d) alignment=%s\n",
607                       k, lno, ln->x1, ln->y1, ln->x2, ln->y2, 
608                       coarse(ln->x1), fine(ln->x1), 
609                       coarse(ln->y1), fine(ln->y1), 
610                       coarse(ln->x2), fine(ln->x2), 
611                       coarse(ln->y2), fine(ln->y2),
612                       ln->updown == 0 ? "vert" : ln->updown == 1 ? "up" : "down")
613                );
614     }
615     for(cscl=startscan; cscl<stopscan; cscl++) {
616       tempy = min(cscl*16+16, pset[i+1].y);
617       POLY_DEB( printf("evaluating scan line %d \n", cscl) );
618       for(k=0; k<clc-1; k+=2) {
619         POLY_DEB( printf("evaluating slice %d\n", k) );
620         render_slice_scanline(&templine, cscl, lset+tllist[k].n, lset+tllist[k+1].n);
621       }
622       if (16*coarse(tempy) == tempy) {
623         POLY_DEB( printf("flushing scan line %d\n", cscl) );
624         scanline_flush(im, &templine, cscl, val);
625         ss_scanline_reset(&templine);
626       }
627       /*
628         else {
629         scanline_flush(im, &templine, cscl, val);
630         ss_scanline_reset(&templine);
631         return 0;
632         }
633       */
634     }
635   } /* Intervals */
636   if (16*coarse(tempy) != tempy) 
637     scanline_flush(im, &templine, cscl-1, val);
638
639   ss_scanline_exorcise(&templine);
640   myfree(pset);
641   myfree(lset);
642   myfree(tllist);
643   
644 } /* Function */