polygon intersector: improved test routines
[swftools.git] / lib / gfxpoly.c
1 /* gfxpoly.c 
2
3    Various boolean polygon functions.
4
5    Part of the swftools package.
6
7    Copyright (c) 2005 Matthias Kramm <kramm@quiss.org> 
8
9    This program is free software; you can redistribute it and/or modify
10    it under the terms of the GNU General Public License as published by
11    the Free Software Foundation; either version 2 of the License, or
12    (at your option) any later version.
13
14    This program is distributed in the hope that it will be useful,
15    but WITHOUT ANY WARRANTY; without even the implied warranty of
16    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
17    GNU General Public License for more details.
18
19    You should have received a copy of the GNU General Public License
20    along with this program; if not, write to the Free Software
21    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA */
22
23 #include "../config.h"
24 #include "gfxdevice.h"
25 #include "gfxtools.h"
26 #include "gfxpoly.h"
27 #include "mem.h"
28 #ifdef INTERNAL_LIBART
29 #include "art/libart.h"
30 #include "art/art_svp_intersect.h"
31 #include "art/art_svp_ops.h"
32 #else
33 #include <libart_lgpl/libart.h>
34 #include <libart_lgpl/art_svp_intersect.h>
35 #include <libart_lgpl/art_svp_ops.h>
36 #endif
37 #include "log.h"
38 #include <assert.h>
39 #include <memory.h>
40 #include <math.h>
41
42 #define PERTURBATE
43 //#define SHEAR
44 //#define DEBUG
45
46 //----------------------------------------- svp renderer ----------------------------------------
47
48 typedef struct {
49     int xmin;
50     int ymin;
51     int xmax;
52     int ymax;
53     int width;
54     int height;
55 } intbbox_t;
56
57 typedef struct _renderpoint
58 {
59     double x;
60     signed char direction;
61 } renderpoint_t;
62
63 typedef struct _renderline
64 {
65     renderpoint_t*points;
66     int size;
67     int num;
68 } renderline_t;
69
70 typedef struct _renderbuf
71 {
72     intbbox_t bbox;
73     int width;
74     int height;
75     double zoom;
76     renderline_t*lines;
77 } renderbuf_t;
78
79 static inline void add_pixel(renderbuf_t*buf, double x, int y, signed char direction)
80 {
81     renderpoint_t p;
82     p.x = x;
83     p.direction = direction;
84
85     if(x >= buf->bbox.xmax || y >= buf->bbox.ymax || y < buf->bbox.ymin) 
86         return;
87     renderline_t*l = &buf->lines[y-buf->bbox.ymin];
88
89     if(l->num == l->size) {
90         l->size += 32;
91         l->points = (renderpoint_t*)rfx_realloc(l->points, l->size * sizeof(renderpoint_t));
92     }
93     l->points[l->num] = p;
94     l->num++;
95 }
96 #define CUT 0.5
97 #define INT(x) ((int)((x)+16)-16)
98 static void add_line(renderbuf_t*buf, double x1, double y1, double x2, double y2, signed char direction)
99 {
100     x1 *= buf->zoom;
101     y1 *= buf->zoom;
102     x2 *= buf->zoom;
103     y2 *= buf->zoom;
104     double diffx, diffy;
105     double ny1, ny2, stepx;
106     if(y2 < y1) {
107         double x,y;
108         x = x1;x1 = x2;x2=x;
109         y = y1;y1 = y2;y2=y;
110     }
111     diffx = x2 - x1;
112     diffy = y2 - y1;
113     ny1 = INT(y1)+CUT;
114     ny2 = INT(y2)+CUT;
115
116     if(ny1 < y1) {
117         ny1 = INT(y1) + 1.0 + CUT;
118     }
119     if(ny2 >= y2) {
120         ny2 = INT(y2) - 1.0 + CUT;
121     }
122     if(ny1 > ny2)
123         return;
124
125     stepx = diffx/diffy;
126     x1 = x1 + (ny1-y1)*stepx;
127     x2 = x2 + (ny2-y2)*stepx;
128
129     int posy=INT(ny1);
130     int endy=INT(ny2);
131     double posx=0;
132     double startx = x1;
133
134     while(posy<=endy) {
135         double xx = startx + posx;
136         add_pixel(buf, xx, posy, direction);
137         posx+=stepx;
138         posy++;
139     }
140 }
141
142 static int compare_renderpoints(const void * _a, const void * _b)
143 {
144     renderpoint_t*a = (renderpoint_t*)_a;
145     renderpoint_t*b = (renderpoint_t*)_b;
146     if(a->x < b->x) return -1;
147     if(a->x > b->x) return 1;
148     return 0;
149 }
150
151 static void fill_bitwise(unsigned char*line, int x1, int x2)
152 {
153     int p1 = x1>>3;
154     int p2 = x2>>3;
155     int b1 = 0xff >> (x1&7);
156     int b2 = 0xff << (1+7-(x2&7));
157     if(p1==p2) {
158         line[p1] |= b1&b2;
159     } else {
160         line[p1] |= b1;
161         memset(&line[p1+1], 255, p2-(p1+1));
162         line[p2] = b2;
163     }
164 }
165
166 unsigned char* render_svp(ArtSVP*svp, intbbox_t*bbox, double zoom, ArtWindRule rule)
167 {
168     renderbuf_t _buf, *buf=&_buf;
169     buf->width = (bbox->xmax - bbox->xmin);
170     buf->height = (bbox->ymax - bbox->ymin);
171     buf->bbox = *bbox;
172     buf->zoom = zoom;
173     int width8 = (buf->width+7) >> 3;
174     unsigned char* image = (unsigned char*)malloc(width8*buf->height);
175     memset(image, 0, width8*buf->height);
176     
177     buf->lines = (renderline_t*)rfx_alloc(buf->height*sizeof(renderline_t));
178     int y;
179     for(y=0;y<buf->height;y++) {
180         memset(&buf->lines[y], 0, sizeof(renderline_t));
181         buf->lines[y].points = 0;
182         buf->lines[y].num = 0;
183     }
184
185     int t;
186     for(t=0;t<svp->n_segs;t++) {
187         ArtSVPSeg* seg = &svp->segs[t];
188         int s;
189         for(s=0;s<seg->n_points-1;s++) {
190             int dir = seg->dir? 1 : -1;
191             add_line(buf, seg->points[s].x, seg->points[s].y, seg->points[s+1].x, seg->points[s+1].y, dir);
192         }
193     }
194     for(y=0;y<buf->height;y++) {
195         renderpoint_t*points = buf->lines[y].points;
196         unsigned char*line = &image[width8*y];
197         int n;
198         int num = buf->lines[y].num;
199         int wind = 0;
200         qsort(points, num, sizeof(renderpoint_t), compare_renderpoints);
201         int lastx = 0;
202         int fill = 0;
203
204         for(n=0;n<num;n++) {
205             renderpoint_t*p = &points[n];
206             int x = (int)(p->x - bbox->xmin);
207             
208             if(x < lastx)
209                 x = lastx;
210             if(x > buf->width) {
211                 break;
212             }
213             if(fill && x!=lastx) {
214                 fill_bitwise(line, lastx, x);
215             }
216             wind += p->direction;
217             if(rule == ART_WIND_RULE_INTERSECT) {
218                 fill = wind>=2;
219             } else if (rule == ART_WIND_RULE_NONZERO) {
220                 fill = wind!=0;
221             } else if (rule == ART_WIND_RULE_ODDEVEN) {
222                 fill = wind&1;
223             } else { // rule == ART_WIND_RULE_POSITIVE
224                 fill = wind>=1;
225             }
226             lastx = x;
227         }
228         if(fill && lastx!=buf->width)
229             fill_bitwise(line, lastx, buf->width);
230     }
231     
232     for(y=0;y<buf->height;y++) {
233         if(buf->lines[y].points) {
234             free(buf->lines[y].points);
235         }
236         memset(&buf->lines[y], 0, sizeof(renderline_t));
237     }
238     free(buf->lines);buf->lines=0;
239     return image;
240 }
241
242 #define MAX_WIDTH 8192
243 #define MAX_HEIGHT 4096
244
245 intbbox_t get_svp_bbox(ArtSVP*svp, double zoom)
246 {
247     int t;
248     intbbox_t b = {0,0,0,0};
249     if(svp->n_segs && svp->segs[0].n_points) {
250         b.xmin = svp->segs[0].points[0].x;
251         b.ymin = svp->segs[0].points[0].y;
252         b.xmax = svp->segs[0].points[0].x;
253         b.ymax = svp->segs[0].points[0].y;
254     }
255     for(t=0;t<svp->n_segs;t++) {
256         ArtSVPSeg* seg = &svp->segs[t];
257         int s;
258         for(s=0;s<seg->n_points;s++) {
259             double x = seg->points[s].x*zoom;
260             double y = seg->points[s].y*zoom;
261             int x1 = floor(x);
262             int x2 = ceil(x);
263             int y1 = floor(y);
264             int y2 = ceil(y);
265             if(x1 < b.xmin) b.xmin = x1;
266             if(y1 < b.ymin) b.ymin = y1;
267             if(x2 > b.xmax) b.xmax = x2;
268             if(y2 > b.ymax) b.ymax = y2;
269         }
270     }
271     if(b.xmax > (int)(MAX_WIDTH*zoom))
272         b.xmax = (int)(MAX_WIDTH*zoom);
273     if(b.ymax > (int)(MAX_HEIGHT*zoom))
274         b.ymax = (int)(MAX_HEIGHT*zoom);
275     if(b.xmin < -(int)(MAX_WIDTH*zoom))
276         b.xmin = -(int)(MAX_WIDTH*zoom);
277     if(b.ymin < -(int)(MAX_HEIGHT*zoom))
278         b.ymin = -(int)(MAX_HEIGHT*zoom);
279     
280     if(b.xmin > b.xmax) 
281         b.xmin = b.xmax;
282     if(b.ymin > b.ymax) 
283         b.ymin = b.ymax;
284
285     b.width = b.xmax - b.xmin;
286     b.height = b.ymax - b.ymin;
287     return b;
288 }
289
290 #define B11100000 0xe0
291 #define B01110000 0x70
292 #define B00111000 0x38
293 #define B00011100 0x1c
294 #define B00001110 0x0e
295 #define B00000111 0x07
296 #define B10000000 0x80
297 #define B01000000 0x40
298 #define B00100000 0x20
299 #define B00010000 0x10
300 #define B00001000 0x08
301 #define B00000100 0x04
302 #define B00000010 0x02
303 #define B00000001 0x01
304
305 int compare_bitmaps(intbbox_t*bbox, unsigned char*data1, unsigned char*data2)
306 {
307     if(!data1 || !data2) 
308         return 0;
309     int x,y;
310     int height = bbox->height;
311     int width = bbox->width;
312     int width8 = (width+7) >> 3;
313     unsigned char*l1 = &data1[width8];
314     unsigned char*l2 = &data2[width8];
315     for(y=1;y<height-1;y++) {
316         for(x=0;x<width8;x++) {
317             unsigned a = l1[x-width8] & l1[x] & l1[x+width8];
318             unsigned b = l2[x-width8] & l2[x] & l2[x+width8];
319             
320             if((a&B11100000) && !(l2[x]&B01000000)) return 0;
321             if((a&B01110000) && !(l2[x]&B00100000)) return 0;
322             if((a&B00111000) && !(l2[x]&B00010000)) return 0;
323             if((a&B00011100) && !(l2[x]&B00001000)) return 0;
324             if((a&B00001110) && !(l2[x]&B00000100)) return 0;
325             if((a&B00000111) && !(l2[x]&B00000010)) return 0;
326
327             if((b&B11100000) && !(l1[x]&B01000000)) return 0;
328             if((b&B01110000) && !(l1[x]&B00100000)) return 0;
329             if((b&B00111000) && !(l1[x]&B00010000)) return 0;
330             if((b&B00011100) && !(l1[x]&B00001000)) return 0;
331             if((b&B00001110) && !(l1[x]&B00000100)) return 0;
332             if((b&B00000111) && !(l1[x]&B00000010)) return 0;
333         }
334         l1 += width8;
335         l2 += width8;
336     }
337     return 1;
338 }
339
340
341 //-----------------------------------------------------------------------------------------------
342
343 static ArtVpath* gfxline_to_ArtVpath(gfxline_t*line, char fill)
344 {
345     ArtVpath *vec = NULL;
346     int pos=0,len=0;
347     gfxline_t*l2;
348     double x=0,y=0;
349
350     /* factor which determines into how many line fragments a spline is converted */
351     double subfraction = 2.4;//0.3
352
353     l2 = line;
354     while(l2) {
355         if(l2->type == gfx_moveTo) {
356             pos ++;
357         } else if(l2->type == gfx_lineTo) {
358             pos ++;
359         } else if(l2->type == gfx_splineTo) {
360             int parts = (int)(sqrt(fabs(l2->x-2*l2->sx+x) + fabs(l2->y-2*l2->sy+y))*subfraction);
361             if(!parts) parts = 1;
362             pos += parts + 1;
363         }
364         x = l2->x;
365         y = l2->y;
366         l2 = l2->next;
367     }
368     pos++;
369     len = pos;
370
371     vec = art_new (ArtVpath, len+1);
372
373     pos = 0;
374     l2 = line;
375     int lastmove=-1;
376     while(l2) {
377         if(l2->type == gfx_moveTo) {
378             vec[pos].code = ART_MOVETO_OPEN;
379             vec[pos].x = l2->x;
380             vec[pos].y = l2->y;
381             lastmove=pos;
382             pos++; 
383             assert(pos<=len);
384         } else if(l2->type == gfx_lineTo) {
385             vec[pos].code = ART_LINETO;
386             vec[pos].x = l2->x;
387             vec[pos].y = l2->y;
388             pos++; 
389             assert(pos<=len);
390         } else if(l2->type == gfx_splineTo) {
391             int i;
392             int parts = (int)(sqrt(fabs(l2->x-2*l2->sx+x) + fabs(l2->y-2*l2->sy+y))*subfraction);
393             if(!parts) parts = 1;
394             double stepsize = 1.0/parts;
395             for(i=0;i<=parts;i++) {
396                 double t = (double)i*stepsize;
397                 vec[pos].code = ART_LINETO;
398                 vec[pos].x = l2->x*t*t + 2*l2->sx*t*(1-t) + x*(1-t)*(1-t);
399                 vec[pos].y = l2->y*t*t + 2*l2->sy*t*(1-t) + y*(1-t)*(1-t);
400                 pos++;
401                 assert(pos<=len);
402             }
403         }
404         x = l2->x;
405         y = l2->y;
406        
407         /* let closed line segments start w/ MOVETO instead of MOVETO_OPEN */
408         if(lastmove>=0 && l2->type!=gfx_moveTo && (!l2->next || l2->next->type == gfx_moveTo)) {
409             if(vec[lastmove].x == l2->x &&
410                vec[lastmove].y == l2->y) {
411                 assert(vec[lastmove].code == ART_MOVETO_OPEN);
412                 vec[lastmove].code = ART_MOVETO;
413             }
414         }
415
416         l2 = l2->next;
417     }
418     vec[pos++].code = ART_END;
419     assert(pos == len);
420
421     if(!fill) {
422         /* Fix "dotted" lines. Those are lines where singular points are created
423            by a moveto x,y lineto x,y combination. We "fix" these by shifting the
424            point in the lineto a little bit to the right 
425            These should only occur in strokes, not in fills, so do this only
426            when we know we're not filling.
427          */
428         int t;
429         for(t=0;vec[t].code!=ART_END;t++) {
430             if(t>0 && (vec[t-1].code==ART_MOVETO_OPEN || vec[t-1].code==ART_MOVETO) 
431                    && vec[t].code==ART_LINETO && vec[t+1].code!=ART_LINETO &&
432                 vec[t-1].x == vec[t].x && 
433                 vec[t-1].y == vec[t].y) {
434                 vec[t].x += 0.01;
435             }
436         }
437     }
438
439     /* Find adjacent identical points. If an ajdacent pair of identical
440        points is found, the second one is removed.
441        So moveto x,y lineto x,y becomes moveto x,y
442           lineto x,y lineto x,y becomes lineto x,y
443           lineto x,y moveto x,y becomes lineto x,y
444           moveto x,y moveto x,y becomes moveto x,y
445           lineto x,y lineto x2,y2 becomes lineto x2,y2 (if dir(x,y) ~= dir(x2,y2))
446      */
447     pos = 0;
448     int outpos = 0;
449     while(1)
450     {
451         if(vec[pos].code == ART_END) {
452             vec[outpos++] = vec[pos++];
453             break;
454         }
455
456         char samedir = 0, samepoint = 0;
457         if(outpos) {
458             double dx = vec[pos].x-vec[pos-1].x;
459             double dy = vec[pos].y-vec[pos-1].y;
460             /*if(pos<len-1 && vec[pos].code == ART_LINETO && vec[pos+1].code == ART_LINETO) {
461                 double dx2 = vec[pos+1].x-vec[pos].x;
462                 double dy2 = vec[pos+1].y-vec[pos].y;
463                 if(fabs(dx*dy2 - dy*dx2) < 0.0001 && dx*dx2 + dy*dy2 >= 0) {
464                     samedir=1;
465                 }
466             }*/
467             if(fabs(dx) + fabs(dy) < 0.0001) {
468                 samepoint=1;
469             }
470         }
471         if(!samepoint && !samedir) {
472             vec[outpos++] = vec[pos++];
473         } else {
474             pos++; // skip
475         }
476     }
477
478     return vec;
479 }
480
481 static void shear_svp(ArtSVP*svp, double v)
482 {
483     /* do a "shearing" on the polygon. We do this to eliminate all
484        horizontal lines (which libart can't handle properly, even though
485        it tries). */
486
487     int t;
488     for(t=0;t<svp->n_segs;t++) {
489         ArtSVPSeg* seg = &svp->segs[t];
490         int s;
491         for(s=0;s<seg->n_points;s++) {
492             ArtPoint* point = &seg->points[s];
493             point->y += point->x*v;
494         }
495     }
496 }
497
498 static double find_shear_value(ArtSVP*svp)
499 {
500     /* We try random values until we find one
501        where there's no slope below a given value, or if that fails,
502        at least no slope of 0 */
503
504     double v = 0;
505     int tries = 0;
506     while(1) {
507         char fail = 0;
508         int t;
509         for(t=0;t<svp->n_segs;t++) {
510             ArtSVPSeg* seg = &svp->segs[t];
511             int s;
512             for(s=0;s<seg->n_points-1;s++) {
513                 ArtPoint* p1 = &seg->points[s];
514                 ArtPoint* p2 = &seg->points[s+1];
515                 double dx = p2->x - p1->x;
516                 double dy = p2->y - p1->y;
517                 dy += dx*v;
518                 if(dy==0) {
519                     fail = 1;
520                     break;
521                 }
522                 if(tries<100 && dx && fabs(dy / dx) < 1e-5) {
523                     fail = 1;
524                     break;
525                 }
526             }
527             if(fail) 
528                 break;
529         }
530         if(!fail) 
531             break;
532 #ifdef HAVE_LRAND48
533         v = lrand48() / 2000000000.0;;
534 #elif HAVE_RAND
535         v = rand() / 2000000000.0;
536 #else
537 #error "no lrand48()/rand() implementation found"
538 #endif
539         tries++;
540     }
541     return v;
542 }
543
544 void show_path(ArtSVP*path)
545 {
546     int t;
547     printf("Segments: %d\n", path->n_segs);
548     for(t=0;t<path->n_segs;t++) {
549         ArtSVPSeg* seg = &path->segs[t];
550         printf("Segment %d: %d points, %s, BBox: (%f,%f,%f,%f)\n", 
551                 t, seg->n_points, seg->dir==0?"UP  ":"DOWN",
552                 seg->bbox.x0, seg->bbox.y0, seg->bbox.x1, seg->bbox.y1);
553         int p;
554         for(p=0;p<seg->n_points;p++) {
555             ArtPoint* point = &seg->points[p];
556             printf("        (%f,%f)\n", point->x, point->y);
557         }
558     }
559     printf("\n");
560 }
561
562
563 typedef struct svp_segment_part
564 {
565     double y1;
566     double y2;
567     char active;
568 } svp_segment_part_t;
569
570 int compare_double(const void*_y1, const void*_y2)
571 {
572     const double*y1 = _y1;
573     const double*y2 = _y2;
574     if(*y1<*y2) return -1;
575     if(*y1>*y2) return 1;
576     else return 0;
577 }
578
579 int compare_seg_start(const void*_s1, const void*_s2)
580 {
581     svp_segment_part_t*s1 = *(svp_segment_part_t**)_s1;
582     svp_segment_part_t*s2 = *(svp_segment_part_t**)_s2;
583     if(s1->y1 < s2->y1) return -1;
584     if(s1->y1 > s2->y1) return 1;
585     else return 0;
586 }
587
588 int compare_seg_end(const void*_s1, const void*_s2)
589 {
590     svp_segment_part_t*s1 = *(svp_segment_part_t**)_s1;
591     svp_segment_part_t*s2 = *(svp_segment_part_t**)_s2;
592     if(s1->y2 < s2->y2) return -1;
593     if(s1->y2 > s2->y2) return 1;
594     else return 0;
595 }
596
597 void clean_svp(ArtSVP*svp)
598 {
599     int t;
600     int oldpoints = 0;
601     int newpoints = 0;
602     int oldsegs = 0;
603     int newsegs = 0;
604     for(t=0;t<svp->n_segs;t++) {
605         ArtSVPSeg* seg = &svp->segs[t];
606         int p;
607         int pos=0;
608         double lasty = 0;
609         oldpoints += seg->n_points;
610         for(p=0;p<seg->n_points;p++) {
611             ArtPoint* p1 = &seg->points[p];
612             if(!p || lasty!=p1->y) {
613                 seg->points[pos] = seg->points[p];
614                 pos++;
615                 lasty = p1->y;
616             }
617         }
618         seg->n_points = pos;
619         newpoints += seg->n_points;
620     }
621     int pos = 0;
622     oldsegs = svp->n_segs;
623     for(t=0;t<svp->n_segs;t++) {
624         if(svp->segs[t].n_points > 1) {
625             svp->segs[pos++] = svp->segs[t];
626         }
627     }
628     svp->n_segs = pos;
629     newsegs = svp->n_segs;
630     if(newsegs < oldsegs || newpoints < oldpoints) {
631         msg("<verbose> Simplified polygon from %d points to %d points, %d to %d segs",
632                 oldpoints, newpoints, oldsegs, newsegs);
633     }
634 }
635
636 int check_svp(ArtSVP*svp)
637 {
638     /* count number of y coordinates and segs */
639     int t;
640     int num_points = 0;
641     int num_segs = 0;
642     for(t=0;t<svp->n_segs;t++) {
643         if(!svp->segs[t].n_points) {
644             msg("<error> svp contains segment with zero points\n");
645             return 0;
646         }
647         num_points += svp->segs[t].n_points;
648         num_segs += svp->segs[t].n_points - 1;
649     }
650
651     /* create segs and ys */
652     double*y = malloc(sizeof(double)*num_points);
653     svp_segment_part_t*segs = malloc(sizeof(svp_segment_part_t)*num_segs);
654     svp_segment_part_t**seg_start = malloc(sizeof(svp_segment_part_t*)*num_segs);
655     svp_segment_part_t**seg_end = malloc(sizeof(svp_segment_part_t*)*num_segs);
656     
657     int c1=0;
658     num_segs = 0;
659     for(t=0;t<svp->n_segs;t++) {
660         ArtSVPSeg* seg = &svp->segs[t];
661         int p;
662         for(p=0;p<seg->n_points;p++) {
663             y[c1++] = seg->points[p].y;
664             assert(c1 <= num_points);
665         }
666         for(p=0;p<seg->n_points-1;p++) {
667             ArtPoint* p1 = &seg->points[p];
668             ArtPoint* p2 = &seg->points[p+1];
669             
670             if(p1->y >= p2->y) {
671                 if(p1->y > p2->y) {
672                     msg("<error> bad svp, y in seg is non-increasing %.16f -> %.16f\n", p1->y, p2->y);
673                 }
674             } else {
675                 segs[num_segs].y1 = p1->y;
676                 segs[num_segs].y2 = p2->y;
677                 segs[num_segs].active = 0;
678                 seg_start[num_segs] = &segs[num_segs];
679                 seg_end[num_segs] = &segs[num_segs];
680                 num_segs++;
681             }
682         }
683     }
684
685     qsort(y, num_points, sizeof(double), compare_double);
686     qsort(seg_start, num_segs, sizeof(svp_segment_part_t*), compare_seg_start);
687     qsort(seg_end, num_segs, sizeof(svp_segment_part_t*), compare_seg_end);
688
689     double lasty = num_points?y[0]+1:0;
690     int num_active = 0;
691     int bleed = 0;
692     double bleedy1=0,bleedy2 = 0;
693     for(t=0;t<num_points;t++) {
694         if(lasty == y[t])
695             continue;
696         int s;
697         for(s=0;s<num_segs;s++) {
698             if(segs[s].y1==y[t]) {
699                 /* segment is starting */
700                 segs[s].active = 1;
701                 num_active++;
702             } else if(segs[s].y2==y[t]) {
703                 segs[s].active = 0;
704                 num_active--;
705             }
706         }
707         if(num_active&1) {
708             bleed = num_active;
709             bleedy1 = y[t];
710         } else {
711             if(bleed) {
712                 bleedy2 = y[t];
713                 break;
714             }
715         }
716         lasty = y[t];
717     }
718     if(bleed) {
719         msg("<verbose> svp bleeds from y=%.16f to y=%.16f (%d/%d active segments)\n", 
720                 bleedy1, bleedy2,
721                 bleed, num_segs);
722         free(y);free(segs);free(seg_start);free(seg_end);
723         return 0;
724     }
725
726     free(y);
727     free(segs);
728     free(seg_start);
729     free(seg_end);
730
731     return 1;
732 }
733
734 void write_svp_postscript(const char*filename, ArtSVP*svp)
735 {
736     if(!svp)
737         return;
738     FILE*fi = fopen(filename, "wb");
739     int i, j;
740     double xmin=0,ymin=0,xmax=0,ymax=0;
741     fprintf(fi, "%% begin\n");
742     for (i = 0; i < svp->n_segs; i++) {
743         for (j = 0; j < svp->segs[i].n_points; j++) {
744             double x = svp->segs[i].points[j].x;
745             double y = svp->segs[i].points[j].y;
746             if(i==0 && j==0) {
747                 xmin = xmax = x;
748                 ymin = ymax = y;
749             } else {
750                 if(x < xmin) xmin = x;
751                 if(x > xmax) xmax = x;
752                 if(y < ymin) ymin = y;
753                 if(y > ymax) ymax = y;
754             }
755         }
756     }
757     if(xmax == xmin) xmax=xmin+1;
758     if(ymax == ymin) ymax=ymin+1;
759
760     for (i = 0; i < svp->n_segs; i++)
761       {
762         fprintf(fi, "%g setgray\n", svp->segs[i].dir ? 0.7 : 0);
763         for (j = 0; j < svp->segs[i].n_points; j++)
764           {
765             //fprintf(fi, "%g %g %s\n",
766             //        20 + 550*(svp->segs[i].points[j].x-xmin)/(xmax-xmin),
767             //        820 - 800*(svp->segs[i].points[j].y-ymin)/(ymax-ymin),
768             //        j ? "lineto" : "moveto");
769             fprintf(fi, "%.32f %.32f %s\n",
770                     svp->segs[i].points[j].x,
771                     svp->segs[i].points[j].y,
772                     j ? "lineto" : "moveto");
773           }
774         fprintf(fi, "stroke\n");
775       }
776
777     fprintf(fi, "showpage\n");
778     fclose(fi);
779 }
780
781 void write_vpath_postscript(const char*filename, ArtVpath*path)
782 {
783     FILE*fi = fopen(filename, "wb");
784     int i, j;
785     double xmin=0,ymin=0,xmax=0,ymax=0;
786     fprintf(fi, "%% begin\n");
787     ArtVpath*p = path;
788     char first = 1;
789     while(p->code != ART_END) {
790         if(p->code == ART_MOVETO || p->code == ART_MOVETO_OPEN) {
791             if(!first)
792                 fprintf(fi, "stroke\n");
793             first = 0;
794             fprintf(fi, "0.0 setgray\n");
795             fprintf(fi, "%.32f %.32f moveto\n", p->x, p->y);
796         } else {
797             fprintf(fi, "%.32f %.32f lineto\n", p->x, p->y);
798         }
799         p++;
800     }
801     if(!first)
802         fprintf(fi, "stroke\n");
803     fprintf(fi, "showpage\n");
804     fclose(fi);
805 }
806
807 void write_gfxline_postscript(const char*filename, gfxline_t*line)
808 {
809     FILE*fi = fopen(filename, "wb");
810     int i, j;
811     fprintf(fi, "%% begin\n");
812     char first = 1;
813     while(line) {
814         if(line->type == gfx_moveTo) {
815             if(!first)
816                 fprintf(fi, "stroke\n");
817             first = 0;
818             fprintf(fi, "0.0 setgray\n");
819             fprintf(fi, "%.32f %.32f moveto\n", line->x, line->y);
820         } else {
821             fprintf(fi, "%.32f %.32f lineto\n", line->x, line->y);
822         }
823         line = line->next;
824     }
825     if(!first)
826         fprintf(fi, "stroke\n");
827     fprintf(fi, "showpage\n");
828     fclose(fi);
829 }
830
831 static int vpath_len(ArtVpath*svp)
832 {
833     int len = 0;
834     while(svp->code != ART_END) {
835         svp ++;
836         len ++;
837     }
838     return len;
839 }
840
841 int gfxline_len(gfxline_t*line)
842 {
843     gfxline_t*i = line;
844     int len = 0;
845     while(i) {
846         len ++;
847         i = i->next;
848     }
849     return len;
850 }
851
852 static ArtVpath*pvpath= 0;
853 static int cmppos(const void*_p1, const void*_p2)
854 {
855     int*p1 = (int*)_p1;
856     int*p2 = (int*)_p2;
857     ArtVpath*v1 = &pvpath[*p1]; 
858     ArtVpath*v2 = &pvpath[*p2]; 
859     if(v1->y < v2->y)
860         return -1;
861     else if(v1->y > v2->y)
862         return 1;
863     else if(v1->x < v2->x)
864         return -2;
865     else if(v1->x > v2->x)
866         return 2;
867     else 
868         return 0;
869 }
870
871 #define PERTURBATION 2e-3
872 static void my_perturb(ArtVpath*vpath)
873 {
874     int t;
875     int len = vpath_len(vpath);
876     int*pos = (int*)malloc(len*sizeof(int));
877     for(t=0;t<len;t++)
878         pos[t] = t;
879     pvpath = vpath;
880     qsort(pos, len, sizeof(int), cmppos);
881     t=0;
882     while(t<len) {
883         int t2=t+1;
884         while(t2<len && !cmppos(&pos[t],&pos[t2])) {
885             t2++;
886         }
887
888         double dx = (PERTURBATION * rand ()) / RAND_MAX - PERTURBATION * 0.5;
889         double dy = (PERTURBATION * rand ()) / RAND_MAX - PERTURBATION * 0.5;
890         int s;
891         for(s=t;s<t2;s++) {
892             vpath[pos[s]].x += dx;
893             vpath[pos[s]].y += dy;
894         }
895         t = t2;
896     }
897     free(pos);
898     pvpath = 0;
899 }
900
901 static ArtSVP* gfxfillToSVP(gfxline_t*line, int perturb)
902 {
903     ArtVpath* vec = gfxline_to_ArtVpath(line, 1);
904     msg("<verbose> Casting gfxline of %d segments (%d line segments) to a gfxpoly", gfxline_len(line), vpath_len(vec));
905
906     if(perturb) {
907         //ArtVpath* vec2 = art_vpath_perturb(vec);
908         //free(vec);
909         //vec = vec2;
910         my_perturb(vec);
911     }
912     ArtSVP *svp = art_svp_from_vpath(vec);
913     free(vec);
914
915     return svp;
916 }
917
918 //#ifdef SHEAR
919 //    double shear = find_shear_value(svp);
920 //    gfxline_t*line =  gfxpoly_to_gfxline((gfxpoly_t*)svp);
921 //    gfxline_t*l = line;
922 //    while(l) {
923 //        l->y += l->x*shear;
924 //        l->sy += l->sx*shear;
925 //        l= l->next;
926 //    }
927 //    svp = (ArtSVP*)gfxpoly_fillToPoly(line);
928 //    printf("shearing svp by %.16f\n", shear);
929 //#endif
930 // ....
931 //#ifdef SHEAR
932 //    art_svp_free(svp);
933 //    shear_svp(result, -shear);
934 //#endif
935
936
937 ArtSVP* run_intersector(ArtSVP*svp, ArtWindRule rule)
938 {
939     ArtSvpWriter * swr = art_svp_writer_rewind_new(rule);
940
941     double zoom = 1.0;
942     intbbox_t bbox = get_svp_bbox(svp, zoom);
943
944     art_svp_intersector(svp, swr);
945     ArtSVP*result = art_svp_writer_rewind_reap(swr);
946     clean_svp(result);
947     if(!check_svp(result)) {
948         current_svp = result;
949         art_report_error(); // might set art_error_in_intersector
950     } else {
951         msg("<verbose> Comparing polygon renderings of size %dx%d and %dx%d", bbox.width, bbox.height, bbox.width, bbox.height);
952         unsigned char*data1 = render_svp(svp, &bbox, zoom, rule);
953         unsigned char*data2 = render_svp(result, &bbox, zoom, ART_WIND_RULE_ODDEVEN);
954         if(!compare_bitmaps(&bbox, data1, data2)) {
955             msg("<verbose> Bad SVP rewinding result- polygons don't match");
956             current_svp = result;
957             art_report_error(); // might set art_error_in_intersector
958         }
959         free(data1);
960         free(data2);
961     }
962
963     if(art_error_in_intersector) {
964         msg("<verbose> Error in polygon processing");
965         art_svp_free(result);
966         art_error_in_intersector=0;
967         return 0;
968     }
969     return result;
970 }
971
972 gfxline_t* gfxpoly_to_gfxline(gfxpoly_t*poly)
973 {
974     ArtSVP*svp = (ArtSVP*)poly;
975     int size = 0;
976     int t;
977     int pos = 0;
978
979     msg("<verbose> Casting polygon of %d segments back to gfxline", svp->n_segs);
980     
981     for(t=0;t<svp->n_segs;t++) {
982         size += svp->segs[t].n_points;
983     }
984     size = size + 1;
985     gfxline_t* lines = (gfxline_t*)rfx_alloc(sizeof(gfxline_t)*size);
986
987     for(t=0;t<svp->n_segs;t++) {
988         ArtSVPSeg* seg = &svp->segs[t];
989         int p;
990         for(p=0;p<seg->n_points;p++) {
991             lines[pos].type = p==0?gfx_moveTo:gfx_lineTo;
992             ArtPoint* point = &seg->points[p];
993             lines[pos].x = point->x;
994             lines[pos].y = point->y;
995             lines[pos].next = &lines[pos+1];
996             pos++;
997         }
998     }
999     if(pos) {
1000         lines[pos-1].next = 0;
1001         return lines;
1002     } else {
1003         return 0;
1004     }
1005 }
1006
1007 gfxpoly_t* gfxpoly_fillToPoly(gfxline_t*line)
1008 {
1009     /* I'm not sure whether doing perturbation here is actually
1010        a good idea- if that line has been run through the machinery
1011        several times already, it might be safer to leave it alone,
1012        since it should already be in a format libart can handle */
1013 #ifdef PERTURBATE
1014     ArtSVP* svp = gfxfillToSVP(line, 1);
1015 #else
1016     ArtSVP* svp = gfxfillToSVP(line, 0);
1017 #endif
1018
1019 #ifdef DEBUG
1020     char filename[80];
1021     static int counter = 0;
1022     sprintf(filename, "svp%d.ps", counter);
1023     write_svp_postscript(filename, svp);
1024     sprintf(filename, "gfxline%d.ps", counter);
1025     write_gfxline_postscript(filename, line);
1026 #endif
1027
1028     /* we do xor-filling by default, so dir is always 1 
1029        (actually for oddeven rewinding it makes no difference, but
1030         it's "cleaner")
1031      */
1032     int t;
1033     for(t=0; t<svp->n_segs; t++) {
1034         svp->segs[t].dir = 1;
1035     }
1036             
1037     /* for some reason, we need to rewind / self-intersect the polygons that gfxfillToSVP
1038        returns- art probably uses a different fill mode (circular?) for vpaths */
1039     ArtSVP*svp_uncrossed=0;
1040    
1041 #ifdef DEBUG
1042     sprintf(filename, "svp%d_in.ps", counter);
1043     write_svp_postscript(filename, svp);
1044     counter++;
1045 #endif
1046
1047     svp_uncrossed = run_intersector(svp, ART_WIND_RULE_ODDEVEN);
1048
1049     art_svp_free(svp);
1050     svp=svp_uncrossed;
1051
1052     return (gfxpoly_t*)svp;
1053 }
1054
1055 gfxpoly_t* gfxpoly_intersect(gfxpoly_t*poly1, gfxpoly_t*poly2)
1056 {
1057     ArtSvpWriter *swr;
1058
1059     static int counter = 0;
1060     
1061     ArtSVP* svp1 = (ArtSVP*)poly1;
1062     ArtSVP* svp2 = (ArtSVP*)poly2;
1063     msg("<verbose> Intersecting two polygons of %d and %d segments", svp1->n_segs, svp2->n_segs);
1064     
1065 #ifdef DEBUG
1066     char filename[80];
1067     sprintf(filename, "isvp%d_src1.ps", counter);
1068     write_svp_postscript(filename, svp1);
1069     sprintf(filename, "isvp%d_src2.ps", counter);
1070     write_svp_postscript(filename, svp2);
1071 #endif
1072     
1073     ArtSVP* svp3 = art_svp_merge (svp1, svp2);
1074
1075 #ifdef DEBUG
1076     sprintf(filename, "isvp%d_src.ps", counter);
1077     write_svp_postscript(filename, svp3);
1078 #endif
1079
1080     //write_svp_postscript("svp.ps", svp3);
1081     ArtSVP*svp_new = run_intersector(svp3, ART_WIND_RULE_INTERSECT);
1082
1083     art_free (svp3); /* shallow free because svp3 contains shared segments */
1084
1085 #ifdef DEBUG
1086     sprintf(filename, "isvp%d.ps", counter);
1087     write_svp_postscript(filename, svp_new);
1088 #endif
1089
1090     counter++;
1091     
1092     //write_svp_postscript("svp_new.ps", svp_new);
1093         
1094     return (gfxpoly_t*)svp_new;
1095 }
1096
1097 gfxpoly_t* gfxpoly_union(gfxpoly_t*poly1, gfxpoly_t*poly2)
1098 {
1099     ArtSVP* svp1 = (ArtSVP*)poly1;
1100     ArtSVP* svp2 = (ArtSVP*)poly2;
1101     msg("<verbose> Unifying two polygons of %d and %d segments", svp1->n_segs, svp2->n_segs);
1102         
1103     ArtSVP* svp = art_svp_union(svp1, svp2);
1104     return (gfxpoly_t*)svp;
1105 }
1106
1107 gfxpoly_t* gfxpoly_strokeToPoly(gfxline_t*line, gfxcoord_t width, gfx_capType cap_style, gfx_joinType joint_style, double miterLimit)
1108 {
1109     ArtVpath* vec = gfxline_to_ArtVpath(line, 0);
1110     msg("<verbose> Casting gfxline of %d segments to a stroke-polygon", gfxline_len(line));
1111
1112     ArtVpath* vec2 = art_vpath_perturb(vec);
1113     free(vec);
1114     vec = vec2;
1115
1116     ArtSVP *svp = art_svp_vpath_stroke (vec,
1117                         (joint_style==gfx_joinMiter)?ART_PATH_STROKE_JOIN_MITER:
1118                         ((joint_style==gfx_joinRound)?ART_PATH_STROKE_JOIN_ROUND:
1119                          ((joint_style==gfx_joinBevel)?ART_PATH_STROKE_JOIN_BEVEL:ART_PATH_STROKE_JOIN_BEVEL)),
1120                         (cap_style==gfx_capButt)?ART_PATH_STROKE_CAP_BUTT:
1121                         ((cap_style==gfx_capRound)?ART_PATH_STROKE_CAP_ROUND:
1122                          ((cap_style==gfx_capSquare)?ART_PATH_STROKE_CAP_SQUARE:ART_PATH_STROKE_CAP_SQUARE)),
1123                         width, //line_width
1124                         miterLimit, //miter_limit
1125                         0.05 //flatness
1126                         );
1127     free(vec);
1128     return (gfxpoly_t*)svp;
1129 }
1130
1131 gfxline_t* gfxline_circularToEvenOdd(gfxline_t*line)
1132 {
1133     msg("<verbose> Converting circular-filled gfxline of %d segments to even-odd filled gfxline", gfxline_len(line));
1134     ArtSVP* svp = gfxfillToSVP(line, 1);
1135
1136     ArtSVP* svp_rewinded;
1137     
1138     svp_rewinded = run_intersector(svp, ART_WIND_RULE_NONZERO);
1139     if(!svp_rewinded) {
1140         art_svp_free(svp);
1141         return 0;
1142     }
1143
1144     gfxline_t* result = gfxpoly_to_gfxline((gfxpoly_t*)svp_rewinded);
1145     art_svp_free(svp);
1146     art_svp_free(svp_rewinded);
1147     return result;
1148 }
1149
1150 gfxpoly_t* gfxpoly_createbox(double x1, double y1,double x2, double y2)
1151 {
1152     ArtVpath *vec = art_new (ArtVpath, 5+1);
1153     vec[0].code = ART_MOVETO;
1154     vec[0].x = x1;
1155     vec[0].y = y1;
1156     vec[1].code = ART_LINETO;
1157     vec[1].x = x1;
1158     vec[1].y = y2;
1159     vec[2].code = ART_LINETO;
1160     vec[2].x = x2;
1161     vec[2].y = y2;
1162     vec[3].code = ART_LINETO;
1163     vec[3].x = x2;
1164     vec[3].y = y1;
1165     vec[4].code = ART_LINETO;
1166     vec[4].x = x1;
1167     vec[4].y = y1;
1168     vec[5].code = ART_END;
1169     vec[5].x = 0;
1170     vec[5].y = 0;
1171     ArtSVP *svp = art_svp_from_vpath(vec);
1172     free(vec);
1173     return (gfxpoly_t*)svp;
1174 }
1175
1176 void gfxpoly_free(gfxpoly_t*poly)
1177 {
1178     ArtSVP*svp = (ArtSVP*)poly;
1179     art_svp_free(svp);
1180 }