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