+/* move the control point so that the spline runs through the original
+ control point */
+void fixcp(qspline*s)
+{
+ plotxy mid,dir;
+ mid.x = (s->end.x + s->start.x)/2;
+ mid.y = (s->end.y + s->start.y)/2;
+ dir.x = s->control.x - mid.x;
+ dir.y = s->control.y - mid.y;
+ s->control.x = mid.x + 2*dir.x;
+ s->control.y = mid.y + 2*dir.y;
+}
+
+int approximate2(struct cspline*s, struct qspline*q, double quality, double start, double end);
+
+void check(struct cspline*s, struct qspline*q, int num)
+{
+ int t;
+ plotxy p = s->start;
+ for(t=0;t<num;t++) {
+ plotxy p2 = q[t].start;
+ if(plotxy_dist(p,p2) > 0.005) {
+ printf("--\n");
+ exit(1);
+ }
+ p = q[t].end;
+ }
+ if(plotxy_dist(p, s->end) > 0.005) {
+ printf("--\n");
+ exit(1);
+ }
+}
+
+int cspline_approximate(struct cspline*s, struct qspline*q, double quality, approximate_method method)
+{
+ if(method==0) {
+ return approximate(s->start, s->control1, s->control2, s->end, q);
+ } else {
+ return approximate2(s, q, quality, 0.0, 1.0);
+ }
+}
+inline plotxy cspline_getpoint(cspline*s, double t)
+{
+ plotxy p;
+ p.x= s->end.x*t*t*t + 3*s->control2.x*t*t*(1-t)
+ + 3*s->control1.x*t*(1-t)*(1-t) + s->start.x*(1-t)*(1-t)*(1-t);
+ p.y= s->end.y*t*t*t + 3*s->control2.y*t*t*(1-t)
+ + 3*s->control1.y*t*(1-t)*(1-t) + s->start.y*(1-t)*(1-t)*(1-t);
+ return p;
+}
+plotxy cspline_getderivative(cspline*s, double t)
+{
+ plotxy d;
+ d.x = s->end.x*(3*t*t) + 3*s->control2.x*(2*t-3*t*t) +
+ 3*s->control1.x*(1-4*t+3*t*t) + s->start.x*(-3+6*t-3*t*t);
+ d.y = s->end.y*(3*t*t) + 3*s->control2.y*(2*t-3*t*t) +
+ 3*s->control1.y*(1-4*t+3*t*t) + s->start.y*(-3+6*t-3*t*t);
+ return d;
+}
+plotxy cspline_getderivative2(cspline*s, double t)
+{
+ plotxy d;
+ d.x = s->end.x*(6*t) + 3*s->control2.x*(2-6*t) +
+ 3*s->control1.x*(-4+6*t) + s->start.x*(6-6*t);
+ d.y = s->end.y*(6*t) + 3*s->control2.y*(2-6*t) +
+ 3*s->control1.y*(-4+6*t) + s->start.y*(6-6*t);
+ return d;
+}
+plotxy cspline_getderivative3(cspline*s, double t)
+{
+ plotxy d;
+ d.x = 6*s->end.x - 18*s->control2.x + 18*s->control1.x - 6*s->start.x;
+ d.y = 6*s->end.y - 18*s->control2.y + 18*s->control1.y - 6*s->start.y;
+ return d;
+}
+void cspline_getequalspacedpoints(cspline*s, float**p, int*num, double dist)
+{
+ plotxy d,next;
+ double t = 0;
+ int end = 0;
+ int pos = 0;
+ float*positions = (float*)malloc(1048576);
+ do
+ {
+ if(t>=1.0) {
+ t = 1.0;
+ end = 1;
+ }
+
+ plotxy d = cspline_getderivative(s, t);
+ plotxy d2 = cspline_getderivative2(s, t);
+
+ double dl = sqrt(d.x*d.x+d.y*d.y);
+ double dl2 = sqrt(d2.x*d2.x+d2.y*d2.y);
+
+ double rdl = dist/dl;
+
+ if(rdl>1.0-t)
+ rdl = 1.0-t;
+
+ plotxy p = cspline_getpoint(s, t);
+ while(plotxy_dist(cspline_getpoint(s, t+rdl), p) > dist) {
+ /* we were ask to divide the spline into dist long fragments,
+ but for the value we estimated even the geometric distance
+ is bigger than 'dist'. Approximate a better value.
+ */
+ rdl = rdl*0.9;
+ }
+
+ positions[pos] = t;
+ t+=rdl;
+ pos++;
+ }
+ while(!end);
+ *num = pos;
+ *p = positions;
+}
+
+plotxy qspline_getpoint(qspline*s, double t)
+{
+ plotxy p;
+ p.x= s->end.x*t*t + 2*s->control.x*t*(1-t) + s->start.x*(1-t)*(1-t);
+ p.y= s->end.y*t*t + 2*s->control.y*t*(1-t) + s->start.y*(1-t)*(1-t);
+ return p;
+}
+plotxy qspline_getderivative(qspline*s, double t)
+{
+ plotxy p;
+ p.x= s->end.x*2*t + 2*s->control.x*(1-2*t) + s->start.x*(-2+2*t);
+ p.y= s->end.y*2*t + 2*s->control.y*(1-2*t) + s->start.y*(-2+2*t);
+ return p;
+}
+plotxy qspline_getderivative2(qspline*s, double t)
+{
+ plotxy p;
+ p.x= s->end.x*2 + 2*s->control.x*(-2) + s->start.x*(2);
+ p.y= s->end.y*2 + 2*s->control.y*(-2) + s->start.y*(2);
+ return p;
+}
+double qspline_getlength(qspline*s)
+{
+ double t = 0;
+ int end = 0;
+ double len;
+ plotxy last = qspline_getpoint(s, 0.0);
+ do {
+ if(t>=1.0) {
+ t = 1.0;
+ end = 1;
+ }
+ plotxy d2 = qspline_getderivative2(s, t);
+ double dl2 = sqrt(d2.x*d2.x+d2.y*d2.y);
+ double rdl = 1.0/dl2;
+ if(rdl>0.01)
+ rdl = 0.01;
+ t+=rdl;
+ plotxy here = qspline_getpoint(s, t);
+ len += plotxy_dist(last, here);
+ last = here;
+ }
+ while(!end);
+ return len;
+}
+void qsplines_getequalspacedpoints(qspline**s, int num, float**p, int*pnum, double acc)
+{
+/* int t;
+ float r[128];
+ for(t=0;t<num;t++) {
+ qspline_getlength();
+ }*/
+ return;
+}
+
+void qsplines_getdrawpoints(qspline*s, int num, float**p, int*pnum, double acc)
+{
+ plotxy d,next;
+ double t = 0;
+ int end = 0;
+ int pos = 0;
+ float*positions = (float*)malloc(1048576);
+ do
+ {
+ if(t>=1.0) {
+ t = 1.0;
+ end = 1;
+ }
+
+ plotxy d = qspline_getderivative(s, t);
+ double dl = sqrt(d.x*d.x+d.y*d.y);
+ double rdl = acc/dl;
+
+ if(rdl>acc)
+ rdl = acc;
+
+ positions[pos] = t;
+ t+=rdl;
+ pos++;
+ }
+ while(!end);
+ *pnum = pos;
+ *p = positions;
+}
+
+
+#define TANGENTS
+
+int approximate2(struct cspline*s, struct qspline*q, double quality, double start, double end)
+{
+ int num=0;
+ plotxy qr1,qr2,cr1,cr2;
+ double dist1,dist2;
+ int t;
+ int recurse = 0;
+ int probes = 15;
+ qspline test;
+ test.start = cspline_getpoint(s, start);
+ test.control = cspline_getpoint(s, (start+end)/2);
+ test.end = cspline_getpoint(s, end);
+ fixcp(&test);
+
+#ifdef TANGENTS
+ if(start< 0.5) {
+ test.control = cspline_getderivative(s, start);
+ test.control.x *= (end-start)/2;
+ test.control.y *= (end-start)/2;
+ test.control.x += test.start.x;
+ test.control.y += test.start.y;
+ } else {
+ test.control = cspline_getderivative(s, end);
+ test.control.x *= -(end-start)/2;
+ test.control.y *= -(end-start)/2;
+ test.control.x += test.end.x;
+ test.control.y += test.end.y;
+ }
+#endif
+
+ for(t=0;t<probes;t++) {
+ double pos = 0.5/(probes*2)*(t*2+1);
+ qr1 = qspline_getpoint(&test, pos);
+ cr1 = cspline_getpoint(s, start+pos*(end-start));
+ dist1 = plotxy_dist(qr1, cr1);
+ if(dist1>quality) {
+ recurse=1;break;
+ }
+ qr2 = qspline_getpoint(&test, (1-pos));
+ cr2 = cspline_getpoint(s, start+(1-pos)*(end-start));
+ dist2 = plotxy_dist(qr2, cr2);
+ if(dist2>quality) {
+ recurse=1;break;
+ }
+ }
+
+ if(recurse && (end-start)>1.0/120) {
+ /* quality is too bad, split it up recursively */
+ num += approximate2(s, q, quality, start, (start+end)/2);
+ q+=num;
+ num += approximate2(s, q, quality, (start+end)/2, end);
+ return num;
+ } else {
+ *q = test;
+ return 1;
+ }
+}
+