+ return fcrossraysxx(&curve[1]);
+}
+
+/* the front-end getting the arguments from gentries:
+ * rays are defined as beginning of curve1 and end of curve 2
+ */
+
+int
+fcrossraysge(
+ GENTRY *ge1,
+ GENTRY *ge2,
+ double *max1,
+ double *max2,
+ double crossdot[2][2]
+)
+{
+ ray[0].x1 = ge1->prev->fx3;
+ ray[0].y1 = ge1->prev->fy3;
+ ray[0].x2 = ge1->fpoints[X][ge1->ftg];
+ ray[0].y2 = ge1->fpoints[Y][ge1->ftg];
+ ray[0].maxp = max1;
+
+ ray[1].x1 = ge2->fx3;
+ ray[1].y1 = ge2->fy3;
+ if(ge2->rtg < 0) {
+ ray[1].x2 = ge2->prev->fx3;
+ ray[1].y2 = ge2->prev->fy3;
+ } else {
+ ray[1].x2 = ge2->fpoints[X][ge2->rtg];
+ ray[1].y2 = ge2->fpoints[Y][ge2->rtg];
+ }
+ ray[1].maxp = max2;
+
+ return fcrossraysxx(crossdot);
+}
+
+/* debugging printout functions */
+
+#if defined(DEBUG_DOTSEG) || defined(DEBUG_DOTCURVE) || defined(DEBUG_APPROXCV)
+
+/* for debugging */
+static
+printdot(
+ double dot[2]
+)
+{
+ fprintf(stderr, "(%g,%g)", dot[0], dot[1]);
+}
+
+static
+printseg(
+ double seg[2][2]
+)
+{
+ putc('[', stderr);
+ printdot(seg[0]);
+ putc(' ', stderr);
+ printdot(seg[1]);
+ putc(']', stderr);
+}
+
+#endif /* DEBUG_* */
+
+/*
+ * Find squared distance from a dot to a line segment
+ */
+
+double
+fdotsegdist2(
+ double seg[2][2 /*X,Y*/],
+ double dot[2 /*X,Y*/]
+)
+{
+#define x1 seg[0][X]
+#define y1 seg[0][Y]
+#define x2 seg[1][X]
+#define y2 seg[1][Y]
+#define xdot dot[X]
+#define ydot dot[Y]
+
+ double dx, dy; /* segment dimensions */
+ double kline, bline; /* segment line formula is y=k*x+b */
+ double kperp, bperp; /* perpendicular from the dot to the line */
+ double xcross, ycross; /* where the perpendicular crosses the segment */
+
+/* handle the situation where the nearest point of the segment is its end */
+#define HANDLE_LIMITS(less12, lesscr1, lesscr2) \
+ if( less12 ) { \
+ if( lesscr1 ) { \
+ xcross = x1; \
+ ycross = y1; \
+ } else if( !(lesscr2) ) { \
+ xcross = x2; \
+ ycross = y2; \
+ } \
+ } else { \
+ if( !(lesscr1) ) { \
+ xcross = x1; \
+ ycross = y1; \
+ } else if( lesscr2 ) { \
+ xcross = x2; \
+ ycross = y2; \
+ } \
+ } \
+ /* end of macro */
+
+
+ dx = x2 - x1;
+ dy = y2 - y1;
+
+ if(fabs(dx) < FEPS) {
+ /* special case - vertical line */
+#ifdef DEBUG_DOTSEG
+ printf("vertical line!\n");
+#endif
+ xcross = x1;
+ ycross = ydot;
+ HANDLE_LIMITS( y1 < y2, ycross < y1, ycross < y2);
+ } else if(fabs(dy) < FEPS) {
+ /* special case - horizontal line */
+#ifdef DEBUG_DOTSEG
+ printf("horizontal line!\n");
+#endif
+ xcross = xdot;
+ ycross = y1;
+ HANDLE_LIMITS( x1 < x2, xcross < x1, xcross < x2)
+ } else {
+ kline = dy/dx;
+ bline = y1 - x1*kline;
+ kperp = -1./kline;
+ bperp = ydot - xdot*kperp;
+
+ xcross = (bline-bperp) / (kperp-kline);
+ ycross = kline*xcross + bline;
+
+ HANDLE_LIMITS( x1 < x2, xcross < x1, xcross < x2)
+ }
+#ifdef DEBUG_DOTSEG
+ printf("crossover at (%g,%g)\n", xcross, ycross);
+#endif
+
+ dx = xdot-xcross;
+ dy = ydot-ycross;
+ return dx*dx+dy*dy;
+#undef x1
+#undef y1
+#undef x2
+#undef y2
+#undef xdot
+#undef ydot
+#undef HANDLE_LIMITS
+}
+
+/* find the weighted quadratic average for the distance of a set
+ * of dots from the curve; also fills out the individual distances
+ * for each dot; if maxp!=NULL then returns the maximal squared
+ * distance in there
+ */
+
+double
+fdotcurvdist2(
+ double curve[4][2 /*X,Y*/ ],
+ struct dot_dist *dots,
+ int ndots, /* number of entries in dots */
+ double *maxp
+)
+{
+ /* a curve is approximated by this many straight segments */
+#define NAPSECT 16
+ /* a curve is divided into this many sections with equal weight each */
+#define NWSECT 4
+ /* table of coefficients for finding the dots on the curve */
+ /* tt[0] is left unused */
+ static double tt[NAPSECT][4];
+ static int havett = 0; /* flag: tt is initialized */
+ /* dots on the curve */
+ double cvd[NAPSECT+1][2 /*X,Y*/];
+ /* sums by sections */
+ double sum[NWSECT];
+ /* counts by sections */
+ double count[NWSECT];
+ int d, i, j;
+ int id1, id2;
+ double dist1, dist2, dist3, dx, dy, x, y;
+ double max = 0.;
+
+ if(!havett) {
+ double t, nt, t2, nt2, step;
+
+ havett++;
+ step = 1. / NAPSECT;
+ t = 0;
+ for(i=1; i<NAPSECT; i++) {
+ t += step;
+ nt = 1 - t;
+ t2 = t*t;
+ nt2 = nt*nt;
+ tt[i][0] = nt2*nt; /* (1-t)^3 */
+ tt[i][1] = 3*nt2*t; /* 3*(1-t)^2*t */
+ tt[i][2] = 3*nt*t2; /* 3*(1-t)*t^2 */
+ tt[i][3] = t2*t; /* t^3 */
+ }
+ }
+
+ for(i=0; i<NWSECT; i++) {
+ sum[i] = 0.;
+ count[i] = 0;
+ }
+
+ /* split the curve into segments */
+ for(d=0; d<2; d++) { /* X and Y */
+ cvd[0][d] = curve[0][d]; /* endpoints */
+ cvd[NAPSECT][d] = curve[3][d];
+ for(i=1; i<NAPSECT; i++) {
+ cvd[i][d] = curve[0][d] * tt[i][0]
+ + curve[1][d] * tt[i][1]
+ + curve[2][d] * tt[i][2]
+ + curve[3][d] * tt[i][3];
+ }
+ }
+
+ for(d=0; d<ndots; d++) {
+#ifdef DEBUG_DOTCURVE
+ printf("dot %d ", d); printdot(dots[d].p); printf(":\n");
+
+ /* for debugging */
+ for(i=0; i< NAPSECT; i++) {
+ dist1 = fdotsegdist2(&cvd[i], dots[d].p);
+ printf(" seg %d ",i); printseg(&cvd[i]); printf(" dist=%g\n", sqrt(dist1));
+ }
+#endif
+
+ x = dots[d].p[X];
+ y = dots[d].p[Y];
+
+ /* find the nearest dot on the curve
+ * there may be up to 2 local minimums - so we start from the
+ * ends of curve and go to the center
+ */
+
+ id1 = 0;
+ dx = x - cvd[0][X];
+ dy = y - cvd[0][Y];
+ dist1 = dx*dx + dy*dy;
+#ifdef DEBUG_DOTCURVE
+ printf(" dot 0 "); printdot(cvd[id1]); printf(" dist=%g\n", sqrt(dist1));
+#endif
+ for(i = 1; i<=NAPSECT; i++) {
+ dx = x - cvd[i][X];
+ dy = y - cvd[i][Y];
+ dist3 = dx*dx + dy*dy;
+#ifdef DEBUG_DOTCURVE
+ printf(" dot %d ",i); printdot(cvd[i]); printf(" dist=%g\n", sqrt(dist3));
+#endif
+ if(dist3 < dist1) {
+ dist1 = dist3;
+ id1 = i;
+ } else
+ break;
+ }
+
+ if(id1 < NAPSECT-1) {
+ id2 = NAPSECT;
+ dx = x - cvd[NAPSECT][X];
+ dy = y - cvd[NAPSECT][Y];
+ dist2 = dx*dx + dy*dy;
+#ifdef DEBUG_DOTCURVE
+ printf(" +dot %d ", id2); printdot(cvd[id2]); printf(" dist=%g\n", sqrt(dist2));
+#endif
+ for(i = NAPSECT-1; i>id1+1; i--) {
+ dx = x - cvd[i][X];
+ dy = y - cvd[i][Y];
+ dist3 = dx*dx + dy*dy;
+#ifdef DEBUG_DOTCURVE
+ printf(" dot %d ",i); printdot(cvd[i]); printf(" dist=%g\n", sqrt(dist3));
+#endif
+ if(dist3 < dist2) {
+ dist2 = dist3;
+ id2 = i;
+ } else
+ break;
+ }
+
+ /* now find which of the local minimums is smaller */
+ if(dist2 < dist1) {
+ id1 = id2;
+ }
+ }
+
+ /* the nearest segment must include the nearest dot */
+ if(id1==0) {
+ dots[d].seg = 0;
+ dots[d].dist2 = fdotsegdist2(&cvd[0], dots[d].p);
+ } else if(id1==NAPSECT) {
+ dots[d].seg = NAPSECT-1;
+ dots[d].dist2 = fdotsegdist2(&cvd[NAPSECT-1], dots[d].p);
+ } else {
+ dist1 = fdotsegdist2(&cvd[id1], dots[d].p);
+ dist2 = fdotsegdist2(&cvd[id1-1], dots[d].p);
+ if(dist2 < dist1) {
+ dots[d].seg = id1-1;
+ dots[d].dist2 = dist2;
+ } else {
+ dots[d].seg = id1;
+ dots[d].dist2 = dist1;
+ }
+ }
+
+ i = dots[d].seg % NWSECT;
+ sum[i] += dots[d].dist2;
+ if(dots[d].dist2 > max)
+ max = dots[d].dist2;
+ count[i]++;
+#ifdef DEBUG_DOTCURVE
+ printf(" best seg %d sect %d dist=%g\n", dots[d].seg, i, sqrt(dots[d].dist2));
+#endif
+ }
+
+ /* calculate the weighted average */
+ id1=0;
+ dist1=0.;
+ for(i=0; i<NWSECT; i++) {
+ if(count[i]==0)
+ continue;
+ id1++;
+ dist1 += sum[i]/count[i];
+ }
+ if(maxp)
+ *maxp = max;
+ if(id1==0) /* no dots, strange */
+ return 0.;
+ else
+ return dist1/id1; /* to get the average distance apply sqrt() */
+}
+
+/*
+ * Approximate a curve matching the giving set of points and with
+ * middle reference points going along the given segments (and no farther
+ * than these segments).
+ */
+
+void
+fapproxcurve(
+ double cv[4][2 /*X,Y*/ ], /* points 0-3 are passed in, points 1,2 - out */
+ struct dot_dist *dots, /* the dots to approximate - distances returned
+ * there may be invalid */
+ int ndots
+)
+{
+ /* b and c are the middle control points */
+#define B 0
+#define C 1
+ /* maximal number of sections on each axis - used for the first step */
+#define MAXSECT 2
+ /* number of sections used for the other steps */
+#define NORMSECT 2
+ /* when the steps become less than this many points, it's time to stop */
+#define STEPEPS 1.
+ double from[2 /*B,C*/], to[2 /*B,C*/];
+ double middf[2 /*B,C*/][2 /*X,Y*/], df;
+ double coef[2 /*B,C*/][MAXSECT];
+ double res[MAXSECT][MAXSECT], thisres, bestres, goodres;
+ int ncoef[2 /*B,C*/], best[2 /*B,C*/], good[2 /*B,C*/];
+ int i, j, k, keepsym;
+ char bc[]="BC";
+ char xy[]="XY";
+
+#ifdef DEBUG_APPROXCV
+ fprintf(stderr, "Curve points:");
+ for(i=0; i<4; i++) {
+ fprintf(stderr, " ");
+ printdot(cv[i]);
+ }
+ fprintf(stderr, "\nDots:");
+ for(i=0; i<ndots; i++) {
+ fprintf(stderr, " ");
+ printdot(dots[i].p);
+ }
+ fprintf(stderr, "\n");
+#endif
+
+ /* load the endpoints and calculate differences */
+ for(i=0; i<2; i++) {
+ /* i is X, Y */
+ middf[B][i] = cv[1][i]-cv[0][i];
+ middf[C][i] = cv[2][i]-cv[3][i];
+
+ /* i is B, C */
+ from[i] = 0.;
+ to[i] = 1.;
+ ncoef[i] = MAXSECT;
+ }
+
+ while(ncoef[B] != 1 || ncoef[C] != 1) {
+ /* prepare the values of coefficients */
+ for(i=0; i<2; i++) { /*B,C*/
+#ifdef DEBUG_APPROXCV
+ fprintf(stderr, "Coefficients by %c(%g,%g):", bc[i], from[i], to[i]);
+#endif
+ df = (to[i]-from[i]) / (ncoef[i]*2);
+ for(j=0; j<ncoef[i]; j++) {
+ coef[i][j] = from[i] + df*(2*j+1);
+#ifdef DEBUG_APPROXCV
+ fprintf(stderr, " %g", coef[i][j]);
+#endif
+ }
+#ifdef DEBUG_APPROXCV
+ fprintf(stderr, "\n");
+#endif
+ }
+ bestres = FBIGVAL;
+ /* i iterates by ncoef[B], j iterates by ncoef[C] */
+ for(i=0; i<ncoef[B]; i++) {
+ for(j=0; j<ncoef[C]; j++) {
+ for(k=0; k<2; k++) { /*X, Y*/
+ cv[1][k] = cv[0][k] + middf[B][k]*coef[B][i];
+ cv[2][k] = cv[3][k] + middf[C][k]*coef[C][j];
+ }
+ res[i][j] = thisres = fdotcurvdist2(cv, dots, ndots, NULL);
+ if(thisres < bestres) {
+ goodres = bestres;
+ good[B] = best[B];
+ good[C] = best[C];
+ bestres = thisres;
+ best[B] = i;
+ best[C] = j;
+ } else if(thisres < goodres) {
+ goodres = thisres;
+ good[B] = i;
+ good[C] = j;
+ }
+#ifdef DEBUG_APPROXCV
+ fprintf(stderr, " at (%g,%g) dist=%g %s\n", coef[B][i], coef[C][j], sqrt(thisres),
+ (best[B]==i && best[C]==j)? "(BEST)":"");
+#endif
+ }
+ }
+#ifdef DEBUG_APPROXCV
+ fprintf(stderr, " best: at (%g, %g) dist=%g\n",
+ coef[B][best[B]], coef[C][best[C]], sqrt(bestres));
+ fprintf(stderr, " B:%d,%d C:%d,%d -- 2nd best: at (%g, %g) dist=%g\n",
+ best[B], good[B], best[C], good[C], coef[B][good[B]], coef[C][good[C]], sqrt(goodres));
+#endif
+
+ if(bestres < (0.1*0.1)) { /* consider it close enough */
+ /* calculate the coordinates to return */
+ for(k=0; k<2; k++) { /*X, Y*/
+ cv[1][k] = cv[0][k] + middf[B][k]*coef[B][best[B]];
+ cv[2][k] = cv[3][k] + middf[C][k]*coef[C][best[C]];
+ }
+#ifdef DEBUG_APPROXCV
+ fprintf(stderr, "quick approximated middle points "); printdot(cv[1]);
+ fprintf(stderr, " "); printdot(cv[2]); fprintf(stderr, "\n");
+#endif
+ return;
+ }
+ keepsym = 0;
+ if(best[B] != best[C] && best[B]-best[C] == good[C]-good[B]) {
+ keepsym = 1;
+#ifdef DEBUG_APPROXCV
+ fprintf(stderr, "keeping symmetry!\n");
+#endif
+ }
+ for(i=0; i<2; i++) { /*B,C*/
+ if(ncoef[i]==1)
+ continue;
+ if(keepsym) {
+ /* try to keep the symmetry */
+ if(best[i] < good[i]) {
+ from[i] = coef[i][best[i]];
+ to[i] = coef[i][good[i]];
+ } else {
+ from[i] = coef[i][good[i]];
+ to[i] = coef[i][best[i]];
+ }
+ } else {
+ df = (to[i]-from[i]) / ncoef[i];
+ from[i] += df*best[i];
+ to[i] = from[i] + df;
+ }
+ if( fabs(df*middf[i][0]) < STEPEPS && fabs(df*middf[i][1]) < STEPEPS) {
+ /* this side has converged */
+ from[i] = to[i] = (from[i]+to[i]) / 2.;
+ ncoef[i] = 1;
+ } else
+ ncoef[i] = NORMSECT;
+ }
+
+ }
+ /* calculate the coordinates to return */
+ for(k=0; k<2; k++) { /*X, Y*/
+ cv[1][k] = cv[0][k] + middf[B][k]*from[B];
+ cv[2][k] = cv[3][k] + middf[C][k]*from[C];
+ }
+#ifdef DEBUG_APPROXCV
+ fprintf(stderr, "approximated middle points "); printdot(cv[1]);
+ fprintf(stderr, " "); printdot(cv[2]); fprintf(stderr, "\n");
+#endif
+#undef B
+#undef C
+#undef MAXSECT
+#undef NORMSECT
+#undef STEPEPS
+}
+
+/*
+ * Find the squared value of the sinus of the angle between the
+ * end of ge1 and the beginning of ge2
+ * The curve must be normalized.
+ */
+
+static double
+fjointsin2(
+ GENTRY *ge1,
+ GENTRY *ge2
+)
+{
+ double d[3][2 /*X,Y*/];
+ double scale1, scale2, len1, len2;
+ int axis;
+
+ if(ge1->rtg < 0) {
+ d[1][X] = ge1->fx3 - ge1->prev->fx3;
+ d[1][Y] = ge1->fy3 - ge1->prev->fy3;
+ } else {
+ d[1][X] = ge1->fx3 - ge1->fpoints[X][ge1->rtg];
+ d[1][Y] = ge1->fy3 - ge1->fpoints[Y][ge1->rtg];
+ }
+ d[2][X] = ge2->fpoints[X][ge2->ftg] - ge2->prev->fx3;
+ d[2][Y] = ge2->fpoints[Y][ge2->ftg] - ge2->prev->fy3;
+
+ len1 = sqrt( d[1][X]*d[1][X] + d[1][Y]*d[1][Y] );
+ len2 = sqrt( d[2][X]*d[2][X] + d[2][Y]*d[2][Y] );
+ /* scale the 2nd segment to the length of 1
+ * and to make sure that the 1st segment is longer scale it to
+ * the length of 2 and extend to the same distance backwards
+ */
+ scale1 = 2./len1;
+ scale2 = 1./len2;
+
+ for(axis=0; axis <2; axis++) {
+ d[0][axis] = -( d[1][axis] *= scale1 );
+ d[2][axis] *= scale2;
+ }
+ return fdotsegdist2(d, d[2]);
+}
+
+#if 0
+/* find the area covered by the curve
+ * (limited by the projections to the X axis)
+ */
+
+static double
+fcvarea(
+ GENTRY *ge
+)
+{
+ double Ly, My, Ny, Py, Qx, Rx, Sx;
+ double area;
+
+ /* y = Ly*t^3 + My*t^2 + Ny*t + Py */
+ Ly = -ge->prev->fy3 + 3*(ge->fy1 - ge->fy2) + ge->fy3;
+ My = 3*ge->prev->fy3 - 6*ge->fy1 + 3*ge->fy2;
+ Ny = 3*(-ge->prev->fy3 + ge->fy1);
+ Py = ge->prev->fy3;
+
+ /* dx/dt = Qx*t^2 + Rx*t + Sx */
+ Qx = 3*(-ge->prev->fx3 + 3*(ge->fx1 - ge->fx2) + ge->fx3);
+ Rx = 6*(ge->prev->fx3 - 2*ge->fx1 + ge->fx2);
+ Sx = 3*(-ge->prev->fx3 + ge->fx1);