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