added AREXT (.a on Unix, .lib on Windows).
[swftools.git] / pdf2swf / spline.cc
1 /* spline.cc
2    Routine to convert cubic splines into quadratic ones.
3
4    Part of the swftools package.
5
6    Copyright (c) 2001,2002,2003 Matthias Kramm <kramm@quiss.org> 
7  
8    This program is free software; you can redistribute it and/or modify
9    it under the terms of the GNU General Public License as published by
10    the Free Software Foundation; either version 2 of the License, or
11    (at your option) any later version.
12
13    This program is distributed in the hope that it will be useful,
14    but WITHOUT ANY WARRANTY; without even the implied warranty of
15    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
16    GNU General Public License for more details.
17
18    You should have received a copy of the GNU General Public License
19    along with this program; if not, write to the Free Software
20    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA */
21
22 #include <stdlib.h>
23 #include <stdio.h>
24 #include <math.h>
25 #include "spline.h"
26
27 static int solve(double a,double b,double c, double*dd)
28 {
29     double det=b*b-4*a*c;
30     int pos = 0;
31     if(det<0) return 0; // we don't do imaginary. not today.
32     if(det==0) { // unlikely, but we have to deal with it.
33      dd[0]=-b/2*a;
34      return (dd[0]>0 && dd[0]<1);
35     }
36
37     dd[pos]=(-b+sqrt(det))/(2*a);
38     if(dd[pos]>0 && dd[pos]<1)
39      pos++;
40     dd[pos]=(-b-sqrt(det))/(2*a);
41     if(dd[pos]>0 && dd[pos]<1)
42      pos++;
43     return pos;
44 }
45
46 struct plotxy splinepos(struct plotxy p0, struct plotxy p1, struct plotxy p2, struct plotxy p3, double d) {
47     struct plotxy p;
48     p.x = (p0.x * d*d*d + p1.x * 3*(1-d)*d*d + p2.x * 3*(1-d)*(1-d)*d + p3.x * (1-d)*(1-d)*(1-d));
49     p.y = (p0.y * d*d*d + p1.y * 3*(1-d)*d*d + p2.y * 3*(1-d)*(1-d)*d + p3.y * (1-d)*(1-d)*(1-d));
50     return p;
51 }
52
53 inline double plotxy_dist(struct plotxy a, struct plotxy b)
54 {
55     double dx = a.x - b.x;
56     double dy = a.y - b.y;
57     return sqrt(dx*dx+dy*dy);
58 }
59
60
61 int wp(double p0,double p1,double p2,double p3,double*dd)
62 {
63     double div= (6*p0-18*p1+18*p2-6*p3);
64     if(!div) return 0;
65     dd[0] = -(6*p1-12*p2+6*p3)/div;
66     return (dd[0]>0 && dd[0]<1);
67 }
68
69 int approximate(struct plotxy p0, struct plotxy p1, struct plotxy p2, struct plotxy p3, struct qspline*q)
70 {
71     double roots[12];
72     int pos = 0;
73     int s,t;
74     struct plotxy myxy[12];
75     struct plotxy last;
76     // the parameters for the solve function are the 1st deviation of a cubic spline
77     roots[pos] = 0;pos++;
78     pos += solve(3*p0.x-9*p1.x+9*p2.x-3*p3.x, 6*p1.x-12*p2.x+6*p3.x,3*p2.x-3*p3.x, &roots[pos]);
79     pos += solve(3*p0.y-9*p1.y+9*p2.y-3*p3.y, 6*p1.y-12*p2.y+6*p3.y,3*p2.y-3*p3.y, &roots[pos]);
80     pos += wp(p0.x,p1.x,p2.x,p3.x,&roots[pos]);
81     pos += wp(p0.x,p1.x,p2.x,p3.x,&roots[pos]);
82     roots[pos] = 1;pos++;
83   
84     // bubblesort - fast enough for 4-6 parameters
85     for(s=0;s<pos;s++)
86     for(t=s+1;t<pos;t++)
87     if(roots[s]>roots[t])
88     {
89        double tmp=roots[s];
90        roots[s]=roots[t];
91        roots[t]=tmp;
92     }
93     for(t=0;t<pos;t++)
94      myxy[t] = splinepos(p0,p1,p2,p3,roots[t]);
95     
96     s=1;
97     last = myxy[0];
98     for(t=1;t<pos;t++)
99     {
100        double dist=plotxy_dist(myxy[t],last);
101        myxy[s]=myxy[t];
102        roots[s]=roots[t];
103        if(dist>0.01 || t==pos-1) 
104        {
105         s++;
106         last=myxy[t];
107        }
108     }
109     pos = s;
110
111     // try 1:curve through 3 points, using the middle of the cubic spline.
112     for(t=0;t<pos-1;t++) {
113 //     circle(myxy[t].x,myxy[t].y,5);
114      struct plotxy control;
115      struct plotxy midpoint = splinepos(p0,p1,p2,p3,(roots[t]+roots[t+1])/2);
116      control.x = midpoint.x + (midpoint.x-(myxy[t].x+myxy[t+1].x)/2);
117      control.y = midpoint.y + (midpoint.y-(myxy[t].y+myxy[t+1].y)/2);
118      //qspline(myxy[t],control,myxy[t+1]);
119      q[t].start=myxy[t];
120      q[t].control=control;
121      q[t].end=myxy[t+1];
122     }
123
124     /*
125     for(t=0;t<pos-1;t++) {
126      plotxy control;
127      vga.setcolor(0xffffff);
128      circle(myxy[t].x,myxy[t].y,5);
129      if(t==0) {
130        //double lenmain = distance(p3,p0);
131        //double lenq = distance(myxy[0],myxy[1]);
132        //control.x = myxy[0].x + (p2.x-p3.x);// /lenmain*lenq;
133        //control.y = myxy[0].y + (p2.y-p3.y);// /lenmain*lenq;
134        plotxy midpoint = splinepos(p0,p1,p2,p3,(roots[t]+roots[t+1])/2);
135        control.x = midpoint.x + (midpoint.x-(myxy[t].x+myxy[t+1].x)/2);
136        control.y = midpoint.y + (midpoint.y-(myxy[t].y+myxy[t+1].y)/2);
137        qspline(myxy[0], control, myxy[1]);
138      } else {
139        control.x = 2*myxy[t].x - last.x;
140        control.y = 2*myxy[t].y - last.y;
141        qspline(myxy[t], control, myxy[t+1]);
142      }
143      last = control;
144     }*/
145     return pos-1;
146 }
147
148 /* move the control point so that the spline runs through the original
149    control point */
150 void fixcp(qspline*s)
151 {
152     plotxy mid,dir;
153     mid.x = (s->end.x + s->start.x)/2;
154     mid.y = (s->end.y + s->start.y)/2;
155     dir.x = s->control.x - mid.x;
156     dir.y = s->control.y - mid.y;
157     s->control.x = mid.x + 2*dir.x;
158     s->control.y = mid.y + 2*dir.y;
159 }
160
161 int approximate2(struct cspline*s, struct qspline*q, double quality, double start, double end);
162
163 void check(struct cspline*s, struct qspline*q, int num)
164 {
165     int t;
166     plotxy p = s->start;
167     for(t=0;t<num;t++) {
168         plotxy p2 = q[t].start;
169         if(plotxy_dist(p,p2) > 0.005) {
170             printf("--\n");
171             exit(1);
172         }
173         p = q[t].end;
174     }
175     if(plotxy_dist(p, s->end) > 0.005) {
176             printf("--\n");
177             exit(1);
178     }
179 }
180
181 int cspline_approximate(struct cspline*s, struct qspline*q, double quality, approximate_method method)
182 {
183     if(method==0) {
184         return approximate(s->start, s->control1, s->control2, s->end, q);
185     } else  {
186         return approximate2(s, q, quality, 0.0, 1.0);
187     }
188 }
189 inline plotxy cspline_getpoint(cspline*s, double t)
190 {
191     plotxy p;
192     p.x= s->end.x*t*t*t + 3*s->control2.x*t*t*(1-t) 
193             + 3*s->control1.x*t*(1-t)*(1-t) + s->start.x*(1-t)*(1-t)*(1-t);
194     p.y= s->end.y*t*t*t + 3*s->control2.y*t*t*(1-t) 
195             + 3*s->control1.y*t*(1-t)*(1-t) + s->start.y*(1-t)*(1-t)*(1-t);
196     return p;
197 }
198 plotxy cspline_getderivative(cspline*s, double t)
199 {
200     plotxy d;
201     d.x = s->end.x*(3*t*t) + 3*s->control2.x*(2*t-3*t*t) + 
202                 3*s->control1.x*(1-4*t+3*t*t) + s->start.x*(-3+6*t-3*t*t);
203     d.y = s->end.y*(3*t*t) + 3*s->control2.y*(2*t-3*t*t) + 
204                 3*s->control1.y*(1-4*t+3*t*t) + s->start.y*(-3+6*t-3*t*t);
205     return d;
206 }
207 plotxy cspline_getderivative2(cspline*s, double t)
208 {
209     plotxy d;
210     d.x = s->end.x*(6*t) + 3*s->control2.x*(2-6*t) + 
211                 3*s->control1.x*(-4+6*t) + s->start.x*(6-6*t);
212     d.y = s->end.y*(6*t) + 3*s->control2.y*(2-6*t) + 
213                 3*s->control1.y*(-4+6*t) + s->start.y*(6-6*t);
214     return d;
215 }
216 plotxy cspline_getderivative3(cspline*s, double t)
217 {
218     plotxy d;
219     d.x = 6*s->end.x - 18*s->control2.x + 18*s->control1.x - 6*s->start.x;
220     d.y = 6*s->end.y - 18*s->control2.y + 18*s->control1.y - 6*s->start.y;
221     return d;
222 }
223 void cspline_getequalspacedpoints(cspline*s, float**p, int*num, double dist)
224 {
225     plotxy d,next;
226     double t = 0;
227     int end = 0;
228     int pos = 0;
229     float*positions = (float*)malloc(1048576);
230     do
231     {
232         if(t>=1.0) {
233             t = 1.0;
234             end = 1;
235         }
236
237         plotxy d = cspline_getderivative(s, t);
238         plotxy d2 = cspline_getderivative2(s, t);
239
240         double dl = sqrt(d.x*d.x+d.y*d.y);
241         double dl2 = sqrt(d2.x*d2.x+d2.y*d2.y);
242     
243         double rdl = dist/dl;
244
245         if(rdl>1.0-t)
246             rdl = 1.0-t;
247
248         plotxy p = cspline_getpoint(s, t);
249         while(plotxy_dist(cspline_getpoint(s, t+rdl), p) > dist) {
250             /* we were ask to divide the spline into dist long fragments,
251                but for the value we estimated even the geometric distance
252                is bigger than 'dist'. Approximate a better value.
253             */
254             rdl = rdl*0.9;
255         }
256         
257         positions[pos] = t;
258         t+=rdl;
259         pos++;
260     }
261     while(!end);
262     *num = pos;
263     *p = positions;
264 }
265
266 plotxy qspline_getpoint(qspline*s, double t)
267 {
268     plotxy p;
269     p.x= s->end.x*t*t + 2*s->control.x*t*(1-t) + s->start.x*(1-t)*(1-t);
270     p.y= s->end.y*t*t + 2*s->control.y*t*(1-t) + s->start.y*(1-t)*(1-t);
271     return p;
272 }
273 plotxy qspline_getderivative(qspline*s, double t)
274 {
275     plotxy p;
276     p.x= s->end.x*2*t + 2*s->control.x*(1-2*t) + s->start.x*(-2+2*t);
277     p.y= s->end.y*2*t + 2*s->control.y*(1-2*t) + s->start.y*(-2+2*t);
278     return p;
279 }
280 plotxy qspline_getderivative2(qspline*s, double t)
281 {
282     plotxy p;
283     p.x= s->end.x*2 + 2*s->control.x*(-2) + s->start.x*(2);
284     p.y= s->end.y*2 + 2*s->control.y*(-2) + s->start.y*(2);
285     return p;
286 }
287 double qspline_getlength(qspline*s)
288 {
289     double t = 0;
290     int end = 0;
291     double len;
292     plotxy last = qspline_getpoint(s, 0.0);
293     do {
294         if(t>=1.0) {
295             t = 1.0;
296             end = 1;
297         }
298         plotxy d2 = qspline_getderivative2(s, t);
299         double dl2 = sqrt(d2.x*d2.x+d2.y*d2.y);
300         double rdl = 1.0/dl2;
301         if(rdl>0.01)
302             rdl = 0.01;
303         t+=rdl;
304         plotxy here = qspline_getpoint(s, t);
305         len += plotxy_dist(last, here);
306         last = here;
307     }
308     while(!end);
309     return len;
310 }
311 void qsplines_getequalspacedpoints(qspline**s, int num, float**p, int*pnum, double acc)
312 {
313 /*    int t;
314     float r[128];
315     for(t=0;t<num;t++) {
316         qspline_getlength();
317     }*/
318     return;
319 }
320
321 void qsplines_getdrawpoints(qspline*s, int num, float**p, int*pnum, double acc)
322 {
323     plotxy d,next;
324     double t = 0;
325     int end = 0;
326     int pos = 0;
327     float*positions = (float*)malloc(1048576);
328     do
329     {
330         if(t>=1.0) {
331             t = 1.0;
332             end = 1;
333         }
334
335         plotxy d = qspline_getderivative(s, t);
336         double dl = sqrt(d.x*d.x+d.y*d.y);
337         double rdl = acc/dl;
338
339         if(rdl>acc)
340             rdl = acc;
341
342         positions[pos] = t;
343         t+=rdl;
344         pos++;
345     }
346     while(!end);
347     *pnum = pos;
348     *p = positions;
349 }
350
351
352 #define TANGENTS
353
354 int approximate2(struct cspline*s, struct qspline*q, double quality, double start, double end)
355 {
356     int num=0;
357     plotxy qr1,qr2,cr1,cr2;
358     double dist1,dist2;
359     int t;
360     int recurse = 0;
361     int probes = 15;
362     qspline test;
363     test.start = cspline_getpoint(s, start);
364     test.control = cspline_getpoint(s, (start+end)/2);
365     test.end = cspline_getpoint(s, end);
366     fixcp(&test);
367
368 #ifdef TANGENTS
369     if(start< 0.5) {
370         test.control = cspline_getderivative(s, start);
371         test.control.x *= (end-start)/2;
372         test.control.y *= (end-start)/2;
373         test.control.x += test.start.x;
374         test.control.y += test.start.y;
375     } else {
376         test.control = cspline_getderivative(s, end);
377         test.control.x *= -(end-start)/2;
378         test.control.y *= -(end-start)/2;
379         test.control.x += test.end.x;
380         test.control.y += test.end.y;
381     }
382 #endif
383
384     for(t=0;t<probes;t++) {
385         double pos = 0.5/(probes*2)*(t*2+1);
386         qr1 = qspline_getpoint(&test, pos);
387         cr1 = cspline_getpoint(s, start+pos*(end-start));
388         dist1 = plotxy_dist(qr1, cr1);
389         if(dist1>quality) {
390             recurse=1;break;
391         }
392         qr2 = qspline_getpoint(&test, (1-pos));
393         cr2 = cspline_getpoint(s, start+(1-pos)*(end-start));
394         dist2 = plotxy_dist(qr2, cr2);
395         if(dist2>quality) {
396             recurse=1;break;
397         }
398     }
399
400     if(recurse && (end-start)>1.0/120) {
401         /* quality is too bad, split it up recursively */
402         num += approximate2(s, q, quality, start, (start+end)/2);
403         q+=num;
404         num += approximate2(s, q, quality, (start+end)/2, end);
405         return num;
406     } else {
407         *q = test;
408         return 1;
409     }
410 }
411