optimized spline approximation
[swftools.git] / lib / gfxtools.c
1 /* gfxtools.c 
2
3    Various utility functions for dealing with gfxdevices.
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 <stdio.h>
24 #include <memory.h>
25 #include <math.h>
26 #include <assert.h>
27 #include "gfxtools.h"
28
29 typedef struct _linedraw_internal
30 {
31     gfxline_t*start;
32     gfxline_t*next;
33 } linedraw_internal_t;
34
35 static void linedraw_moveTo(gfxdrawer_t*d, gfxcoord_t x, gfxcoord_t y)
36 {
37     linedraw_internal_t*i = (linedraw_internal_t*)d->internal;
38     gfxline_t*l = rfx_alloc(sizeof(gfxline_t));
39     l->type = gfx_moveTo;
40     if((int)((d->x * 5120) == (int)(x * 5120)) &&
41        (int)((d->y * 5120) == (int)(y * 5120))) {
42         /* never mind- we're already there */
43         return;
44
45     }
46     d->x = l->x = x;
47     d->y = l->y = y;
48     l->next = 0;
49     if(i->next)
50         i->next->next = l;
51     i->next = l;
52     if(!i->start)
53         i->start = l;
54 }
55 static void linedraw_lineTo(gfxdrawer_t*d, gfxcoord_t x, gfxcoord_t y)
56 {
57     linedraw_internal_t*i = (linedraw_internal_t*)d->internal;
58     gfxline_t*l = rfx_alloc(sizeof(gfxline_t));
59     l->type = gfx_lineTo;
60     d->x = l->x = x;
61     d->y = l->y = y;
62     l->next = 0;
63     if(i->next)
64         i->next->next = l;
65     i->next = l;
66     if(!i->start)
67         i->start = l;
68 }
69 static void linedraw_splineTo(gfxdrawer_t*d, gfxcoord_t sx, gfxcoord_t sy, gfxcoord_t x, gfxcoord_t y)
70 {
71     linedraw_internal_t*i = (linedraw_internal_t*)d->internal;
72     gfxline_t*l = rfx_alloc(sizeof(gfxline_t));
73     l->type = gfx_splineTo;
74     d->x = l->x = x; 
75     d->y = l->y = y;
76     l->sx = sx; 
77     l->sy = sy;
78     l->next = 0;
79     if(i->next)
80         i->next->next = l;
81     i->next = l;
82     if(!i->start)
83         i->start = l;
84 }
85 static void* linedraw_result(gfxdrawer_t*d)
86 {
87     linedraw_internal_t*i = (linedraw_internal_t*)d->internal;
88     void*result = (void*)i->start;
89     rfx_free(i);
90     memset(d, 0, sizeof(gfxdrawer_t));
91     return result;
92 }
93
94 void gfxdrawer_target_gfxline(gfxdrawer_t*d)
95 {
96     linedraw_internal_t*i = (linedraw_internal_t*)rfx_calloc(sizeof(linedraw_internal_t));
97     d->internal = i;
98     d->moveTo = linedraw_moveTo;
99     d->lineTo = linedraw_lineTo;
100     d->splineTo = linedraw_splineTo;
101     d->result = linedraw_result;
102 }
103
104 typedef struct _qspline_abc
105 {
106     double ax,bx,cx;
107     double ay,by,cy;
108 } qspline_abc_t;
109
110 typedef struct qspline_t
111 {
112     gfxpoint_t start;
113     gfxpoint_t control;
114     gfxpoint_t end;
115 } qspline_t;
116
117 typedef struct cspline_t
118 {
119     gfxpoint_t start;
120     gfxpoint_t control1;
121     gfxpoint_t control2;
122     gfxpoint_t end;
123 } cspline_t;
124
125 static void mkspline(qspline_abc_t*s, double x, double y, gfxline_t*l)
126 {
127     /* 
128        Form 1: x = t*t*l->x + 2*t*(1-t)*l->sx + (1-t)*(1-t)*x;
129        Form 2: x = a*t*t + b*t + c
130     */
131     s->cx = x; s->bx = 2*l->sx - 2*x; s->ax = l->x - 2*l->sx + x;
132     s->cy = y; s->by = 2*l->sy - 2*y; s->ay = l->y - 2*l->sy + y;
133 }
134
135 static void spline_get_controlpoint(qspline_abc_t*q, double t1, double t2, double*dx, double*dy)
136 {
137     double dt = t2-t1;
138     double nax = q->ax*dt*dt;
139     double nay = q->ay*dt*dt;
140     double nbx = 2*q->ax*dt*t1 + q->bx*dt;
141     double nby = 2*q->ay*dt*t1 + q->by*dt;
142     double ncx = q->ax*t1*t1 + q->bx*t1 + q->cx;
143     double ncy = q->ay*t1*t1 + q->by*t1 + q->cy;
144     *dx = ncx + nbx/2;
145     *dy = ncy + nby/2;
146 }
147
148 static double get_spline_len(qspline_abc_t*s)
149 {
150     int parts = (int)(sqrt(fabs(s->ax) + fabs(s->ay))*3);
151     int i;
152     double len = 0;
153     double r;
154     double r2;
155     if(parts < 3) parts = 3;
156     r = 1.0/parts;
157     r2 = 1.0/(parts*parts);
158     for(i=0;i<parts;i++)
159     {
160         double dx = s->ax*(2*i+1)*r2 + s->bx*r;
161         double dy = s->ay*(2*i+1)*r2 + s->by*r;
162         len += sqrt(dx*dx+dy*dy);
163     }
164     /*printf("Spline from %f,%f to %f,%f has len %f (%f)\n", s->cx, s->cy, 
165             s->cx + s->bx + s->ax,
166             s->cy + s->by + s->ay, len,
167             sqrt((s->bx + s->ax)*(s->bx + s->ax) + (s->by + s->ay)*(s->by + s->ay))
168             );
169     assert(len+0.5 >= sqrt((s->bx + s->ax)*(s->bx + s->ax) + (s->by + s->ay)*(s->by + s->ay)));
170      */
171     return len;
172 }
173
174 void gfxtool_draw_dashed_line(gfxdrawer_t*d, gfxline_t*line, float*r, float phase)
175 {
176     double x=0,y=0;
177     double linepos,nextpos;
178     char on;
179     int apos;
180
181     if(line && line->type != gfx_moveTo) {
182         fprintf(stderr, "gfxtool: outline doesn't start with a moveTo");
183         return;
184     }
185     if(!r || r[0]<0 || phase<0) {
186         fprintf(stderr, "gfxtool: invalid dashes");
187         return;
188     }
189
190     for(;line;line=line->next) {
191         if(line->type == gfx_moveTo) {
192             d->moveTo(d, line->x, line->y);
193             on = 1; nextpos = r[0]; apos = 0; linepos = 0;
194             x = line->x; y = line->y;
195             while(linepos < phase) {
196                 //printf("[+] linepos: %f, phase: %f, on:%d, apos:%d nextpos:%f\n", linepos, phase, on, apos, nextpos);
197                 linepos += r[apos];
198                 if(linepos < phase) {
199                     on ^= 1;
200                     if(r[++apos]<0)
201                         apos = 0;
202                     nextpos += r[apos];
203                 }
204             }
205             linepos = phase;
206             //printf("[k] linepos: %f, phase: %f, on:%d, apos:%d nextpos:%f \n", linepos, phase, on, apos, nextpos);
207         } else if(line->type == gfx_lineTo) {
208             double dx = line->x - x;
209             double dy = line->y - y;
210             double len = sqrt(dx*dx+dy*dy);
211             double vx;
212             double vy;
213             double lineend = linepos+len;
214             if(len==0)
215                 continue;
216             vx = dx/len;
217             vy = dy/len;
218             assert(nextpos>=linepos);
219             //printf("(line) on:%d apos: %d nextpos: %f, line pos: %f, line end: %f\n", on, apos, nextpos, linepos, linepos+len);
220             while(nextpos<lineend) {
221                 double nx = x + vx*(nextpos-linepos);
222                 double ny = y + vy*(nextpos-linepos);
223                 if(on) {d->lineTo(d, nx,ny);/*printf("lineTo %f\n", nextpos);*/}
224                 else   {d->moveTo(d, nx,ny);/*printf("moveTo %f\n", nextpos);*/}
225                 on^=1;
226                 if(r[++apos]<0)
227                     apos = 0;
228                 nextpos+=r[apos];
229             }
230             linepos = lineend;
231             if(on) {
232                 //printf("lineTo %f\n", 1.0);
233                 d->lineTo(d, line->x,line->y);
234             }
235             x = line->x; y = line->y;
236         } else if(line->type == gfx_splineTo) {
237             qspline_abc_t q;
238             double len, lineend,lastt;
239             mkspline(&q, x, y, line);
240
241             len = get_spline_len(&q);
242             //printf("%f %f -> %f %f, len: %f\n", x, y, line->x, line->y, len);
243             if(len==0)
244                 continue;
245             lineend = linepos+len;
246             lastt = 0;
247             if(nextpos<linepos)
248                 printf("%f !< %f\n", nextpos, linepos);
249             assert(nextpos>=linepos);
250             //printf("(spline) on:%d apos: %d nextpos: %f, line pos: %f, line end: %f\n", on, apos, nextpos, linepos, linepos+len);
251             while(nextpos<lineend) {
252                 double t = (nextpos-linepos)/len;
253                 //printf("%f (%f-%f) apos=%d r[apos]=%f\n", t, nextpos, linepos, apos, r[apos]);
254                 double nx = q.ax*t*t+q.bx*t+q.cx;
255                 double ny = q.ay*t*t+q.by*t+q.cy;
256                 if(on) {
257                     double sx,sy;
258                     spline_get_controlpoint(&q, lastt, t, &sx, &sy);
259                     d->splineTo(d, sx, sy, nx,ny);
260                     //printf("splineTo %f\n", nextpos);
261                 } else  {
262                     d->moveTo(d, nx,ny);
263                     //printf("moveTo %f\n", nextpos);
264                 }
265                 lastt =  t;
266                 on^=1;
267                 if(r[++apos]<0)
268                     apos = 0;
269                 nextpos+=r[apos];
270             }
271             linepos = lineend;
272             if(on) {
273                 double sx,sy;
274                 spline_get_controlpoint(&q, lastt, 1, &sx, &sy);
275                 d->splineTo(d, sx, sy, line->x,line->y);
276                 //printf("splineTo %f\n", 1.0);
277             }
278             x = line->x; y = line->y;
279         }
280     }
281 }
282
283 gfxline_t * gfxline_clone(gfxline_t*line)
284 {
285     gfxline_t*dest = 0;
286     gfxline_t*pos = 0;
287     while(line) {
288         gfxline_t*n = rfx_calloc(sizeof(gfxline_t));
289         *n = *line;
290         n->next = 0;
291         if(!pos) {
292             dest = pos = n;
293         } else {
294             pos->next = n;
295             pos = n;
296         }
297         line = line->next;
298     }
299     return dest;
300 }
301
302 gfxline_t* gfxtool_dash_line(gfxline_t*line, float*dashes, float phase)
303 {
304     gfxdrawer_t d;
305     gfxline_t*result;
306     gfxdrawer_target_gfxline(&d);
307     gfxtool_draw_dashed_line(&d, line, dashes, phase);
308     result= (gfxline_t*)d.result(&d);
309     return result;
310 }
311
312 void gfxline_show(gfxline_t*l, FILE*fi)
313 {
314     while(l) {
315         if(l->type == moveTo) {
316             fprintf(fi, "moveTo %.2f,%.2f\n", l->x, l->y);
317         }
318         if(l->type == lineTo) {
319             fprintf(fi, "lineTo %.2f,%.2f\n", l->x, l->y);
320         }
321         if(l->type == splineTo) {
322             fprintf(fi, "splineTo %.2f,%.2f %.2f,%.2f\n", l->sx, l->sy, l->x, l->y);
323         }
324         l = l->next;
325     }
326 }
327
328 void gfxline_free(gfxline_t*l)
329 {
330     gfxline_t*next;
331     while(l) {
332         next = l->next;
333         l->next = 0;
334         rfx_free(l);
335         l = next;
336     }
337 }
338
339 static inline gfxpoint_t cspline_getpoint(const struct cspline_t*s, double t)
340 {
341     gfxpoint_t p;
342     double tt = t*t;
343     double ttt = tt*t;
344     double mt = (1-t);
345     double mtmt = mt*(1-t);
346     double mtmtmt = mtmt*(1-t);
347     p.x= s->end.x*ttt + 3*s->control2.x*tt*mt
348             + 3*s->control1.x*t*mtmt + s->start.x*mtmtmt;
349     p.y= s->end.y*ttt + 3*s->control2.y*tt*mt
350             + 3*s->control1.y*t*mtmt + s->start.y*mtmtmt;
351     return p;
352 }
353 static gfxpoint_t qspline_getpoint(const qspline_t*s, double t)
354 {
355     gfxpoint_t p;
356     p.x= s->end.x*t*t + 2*s->control.x*t*(1-t) + s->start.x*(1-t)*(1-t);
357     p.y= s->end.y*t*t + 2*s->control.y*t*(1-t) + s->start.y*(1-t)*(1-t);
358     return p;
359 }
360
361 static int approximate3(const cspline_t*s, qspline_t*q, int size, double quality2)
362 {
363     unsigned int gran = 0;
364     unsigned int istep = 0x80000000;
365     unsigned int istart = 0;
366     int num = 0;
367     int level = 0;
368     
369     while(istart<0x80000000)
370     {
371         unsigned int iend = istart + istep;
372         double start = istart/(double)0x80000000;
373         double end = iend/(double)0x80000000;
374         qspline_t test;
375         double pos,qpos;
376         char left = 0,recurse=0;
377         int t;
378         int probes = 15;
379         double dx,dy;
380
381         /* create simple approximation: a qspline_t which run's through the
382            qspline_t point at 0.5 */
383         test.start = cspline_getpoint(s, start);
384         test.control = cspline_getpoint(s, (start+end)/2);
385         test.end = cspline_getpoint(s, end);
386         /* fix the control point:
387            move it so that the new spline does runs through it */
388         test.control.x = -(test.end.x + test.start.x)/2 + 2*(test.control.x);
389         test.control.y = -(test.end.y + test.start.y)/2 + 2*(test.control.y);
390
391         /* depending on where we are in the spline, we either try to match
392            the left or right tangent */
393         if(start<0.5) 
394             left=1;
395         /* get derivative */
396         pos = left?start:end;
397         qpos = pos*pos;
398         test.control.x = s->end.x*(3*qpos) + 3*s->control2.x*(2*pos-3*qpos) + 
399                     3*s->control1.x*(1-4*pos+3*qpos) + s->start.x*(-3+6*pos-3*qpos);
400         test.control.y = s->end.y*(3*qpos) + 3*s->control2.y*(2*pos-3*qpos) + 
401                     3*s->control1.y*(1-4*pos+3*qpos) + s->start.y*(-3+6*pos-3*qpos);
402         if(left) {
403             test.control.x *= (end-start)/2;
404             test.control.y *= (end-start)/2;
405             test.control.x += test.start.x;
406             test.control.y += test.start.y;
407         } else {
408             test.control.x *= -(end-start)/2;
409             test.control.y *= -(end-start)/2;
410             test.control.x += test.end.x;
411             test.control.y += test.end.y;
412         }
413
414 //#define PROBES
415 #ifdef PROBES
416         /* measure the spline's accurancy, by taking a number of probes */
417         for(t=0;t<probes;t++) {
418             gfxpoint_t qr1,qr2,cr1,cr2;
419             double pos = 0.5/(probes*2)*(t*2+1);
420             double dx,dy;
421             double dist1,dist2;
422             qr1 = qspline_getpoint(&test, pos);
423             cr1 = cspline_getpoint(s, start+pos*(end-start));
424
425             dx = qr1.x - cr1.x;
426             dy = qr1.y - cr1.y;
427             dist1 = dx*dx+dy*dy;
428
429             if(dist1>quality2) {
430                 recurse=1;break;
431             }
432             qr2 = qspline_getpoint(&test, (1-pos));
433             cr2 = cspline_getpoint(s, start+(1-pos)*(end-start));
434
435             dx = qr2.x - cr2.x;
436             dy = qr2.y - cr2.y;
437             dist2 = dx*dx+dy*dy;
438
439             if(dist2>quality2) {
440                 recurse=1;break;
441             }
442         }
443 #else // quadratic error: *much* faster!
444
445         /* convert control point representation to 
446            d*x^3 + c*x^2 + b*x + a */
447         dx= s->end.x  - s->control2.x*3 + s->control1.x*3 - s->start.x;
448         dy= s->end.y  - s->control2.y*3 + s->control1.y*3 - s->start.y;
449         
450         /* we need to do this for the subspline between [start,end], not [0,1] 
451            as a transformation of t->a*t+b does nothing to highest coefficient
452            of the spline except multiply it with a^3, we just need to modify
453            d here. */
454         {double m = end-start;
455          dx*=m*m*m;
456          dy*=m*m*m;
457         }
458         
459         /* use the integral over (f(x)-g(x))^2 between 0 and 1
460            to measure the approximation quality. 
461            (it boils down to const*d^2) */
462         recurse = (dx*dx + dy*dy > quality2);
463 #endif
464
465         if(recurse && istep>1 && size-level > num) {
466             istep >>= 1;
467             level++;
468         } else {
469             *q++ = test;
470             num++;
471             istart += istep;
472             while(!(istart & istep)) {
473                 level--;
474                 istep <<= 1;
475             }
476         }
477     }
478     return num;
479 }
480
481 void gfxdraw_conicTo(gfxdrawer_t*draw, double cx, double cy, double tox, double toy)
482 {
483     double c1x = (draw->x + 2 * cx) / 3;
484     double c1y = (draw->y + 2 * cy) / 3;
485     double c2x = (2 * cx + tox) / 3;
486     double c2y = (2 * cy + toy) / 3;
487     gfxdraw_cubicTo(draw, c1x, c1y, c2x, c2y, tox, toy);
488 }
489
490
491 void gfxdraw_cubicTo(gfxdrawer_t*draw, double c1x, double c1y, double c2x, double c2y, double x, double y)
492 {
493     qspline_t q[128];
494     cspline_t c;
495     double maxerror = 0.01;
496     int t,num;
497
498     c.start.x = draw->x;
499     c.start.y = draw->y;
500     c.control1.x = c1x;
501     c.control1.y = c1y;
502     c.control2.x = c2x;
503     c.control2.y = c2y;
504     c.end.x = x;
505     c.end.y = y;
506     
507     num = approximate3(&c, q, 128, maxerror);
508
509     for(t=0;t<num;t++) {
510         gfxpoint_t mid;
511         gfxpoint_t to;
512         mid.x = q[t].control.x;
513         mid.y = q[t].control.y;
514         to.x = q[t].end.x;
515         to.y = q[t].end.y;
516         draw->splineTo(draw, mid.x, mid.y, to.x, to.y);
517     }
518 }
519
520 gfxbbox_t gfxbbox_expand_to_point(gfxbbox_t box, gfxcoord_t x, gfxcoord_t y)
521 {
522     if(box.xmin==0 && box.ymin==0 && box.xmax==0 && box.ymax==0) {
523         box.xmin = x;
524         box.ymin = y;
525         box.xmax = x;
526         box.ymax = y;
527         if(x==0 && y==0) box.xmax = 0.0000001;
528         return box;
529     }
530     if(x < box.xmin)
531         box.xmin = x;
532     if(x > box.xmax)
533         box.xmax = x;
534     if(y < box.ymin)
535         box.ymin = y;
536     if(y > box.ymax)
537         box.ymax = y;
538     return box;
539 }
540
541 gfxbbox_t gfxline_getbbox(gfxline_t*line)
542 {
543     gfxcoord_t x=0,y=0;
544     gfxbbox_t bbox = {0,0,0,0};
545     char last = 0;
546     while(line) {
547         if(line->type == gfx_moveTo) {
548             last = 1;
549         } else if(line->type == gfx_lineTo) {
550             if(last) bbox = gfxbbox_expand_to_point(bbox, x, y);
551             bbox = gfxbbox_expand_to_point(bbox, line->x, line->y);
552             last = 0;
553         } else if(line->type == gfx_splineTo) {
554             if(last) bbox = gfxbbox_expand_to_point(bbox, x, y);
555             bbox = gfxbbox_expand_to_point(bbox, line->sx, line->sy);
556             bbox = gfxbbox_expand_to_point(bbox, line->x, line->y);
557             last = 0;
558         }
559         x = line->x;
560         y = line->y;
561         line = line->next;
562     }
563     return bbox;
564 }
565
566 void gfxline_dump(gfxline_t*line, FILE*fi, char*prefix)
567 {
568     while(line) {
569         if(line->type == gfx_moveTo) {
570             fprintf(fi, "%smoveTo %.2f %.2f\n", prefix, line->x, line->y);
571         } else if(line->type == gfx_lineTo) {
572             fprintf(fi, "%slineTo %.2f %.2f\n", prefix, line->x, line->y);
573         } else if(line->type == gfx_splineTo) {
574             fprintf(fi, "%ssplineTo (%.2f %.2f) %.2f %.2f\n", prefix, line->sx, line->sy, line->x, line->y);
575         }
576         line = line->next;
577     }
578 }
579
580 gfxline_t* gfxline_append(gfxline_t*line1, gfxline_t*line2)
581 {
582     gfxline_t*l = line1;;
583     if(!l)
584         return line2;
585     while(l->next) {
586         l = l->next;
587     }
588     l->next = line2;
589     return line1;
590 }
591
592 void gfxline_transform(gfxline_t*line, gfxmatrix_t*matrix)
593 {
594     while(line) {
595         double x = matrix->m00*line->x + matrix->m10*line->y + matrix->tx;
596         double y = matrix->m01*line->x + matrix->m11*line->y + matrix->ty;
597         line->x = x;
598         line->y = y;
599         if(line->type == gfx_splineTo) {
600             double sx = matrix->m00*line->sx + matrix->m10*line->sy + matrix->tx;
601             double sy = matrix->m01*line->sx + matrix->m11*line->sy + matrix->ty;
602             line->sx = sx;
603             line->sy = sy;
604         }
605         line = line->next;
606     }
607 }
608
609 void gfxmatrix_dump(gfxmatrix_t*m, FILE*fi, char*prefix)
610 {
611     fprintf(fi, "%f %f | %f\n", m->m00, m->m10, m->tx);
612     fprintf(fi, "%f %f | %f\n", m->m01, m->m11, m->ty);
613 }