#include <stdio.h>
#include <math.h>
#include "uigui.h"

char *buf;
typedef struct { short x, y; } XPoint;

int nRealPoints, maxRealPoints;
XPoint *realPoints;

void MydrawCurve(Window w, XPoint *p, int n);

void populaterealpoints()
{
   int x, y;
   for(y=0; y< 600; y++)
      for(x=0; x < 600; x++)
	 if (buf[y * 600 + x] == 0) {
	    realPoints[nRealPoints].x = x;
	    realPoints[nRealPoints].y = y;
	    nRealPoints++;
	    if (nRealPoints >= maxRealPoints) {
	       maxRealPoints *= 2;
	       realPoints = realloc(realPoints, maxRealPoints * sizeof(XPoint));
	       }
	    }
}

void real_floodfill_4(short *buf, int top, int bottom);

/* for interior-defined regions. */
void floodfill_4(int x,int y)
{
   short q[72000];
   int top=0, bottom=2;
   q[0]=x; q[1]=y;
   real_floodfill_4(q, top, bottom);
   }

void real_floodfill_4(short *q, int top, int bottom)
   {
      int x, y;
      while (top != bottom) {
   x = q[top++]; y = q[top++];
   if (top >= 72000) top = 0;
   if (x < 0 || y < 0 || x > 599 || y > 599) continue;
   if (buf[y * 600 + x] == 0) {
      buf[y * 600 + x] = 1;
/*      DrawPoint(x,y); */
      q[bottom++] = x; q[bottom++] = y-1;
      if (bottom >= 72000) bottom = 0;
      q[bottom++] = x; q[bottom++] = y+1;
      if (bottom >= 72000) bottom = 0;
      q[bottom++] = x-1; q[bottom++] = y;
      if (bottom >= 72000) bottom = 0;
      q[bottom++] = x+1; q[bottom++] = y;
      if (bottom >= 72000) bottom = 0;
      }
      }
   }


int main()
{
   short p[] = { 100, 200, 100, 100, 200, 140, 300, 100, 300, 200, 200, 160, 100, 200, 100, 100, 200, 140 };
   buf = calloc(600,600);
   WOpen("drj");
   WAttrib("size=600,600");
   MydrawCurve(stdscr, p, 9);
   Event();
}

/*
 * genCurve - draw a smooth curve through a set of points.
 *  Algorithm from Barry, Phillip J., and Goldman, Ronald N. (1988).
 *  A Recursive Evaluation Algorithm for a class of Catmull-Rom Splines.
 *  Computer Graphics 22(4), 199-204.
 */
void MygenCurve(w, p, n, helper)
Window w;
XPoint *p;
int n;
void (*helper)	(Window, XPoint [], int);
   {
   int    i, j, steps;
   float  ax, ay, bx, by, stepsize, stepsize2, stepsize3;
   float  x, dx, d2x, d3x, y, dy, d2y, d3y;
   XPoint *thepoints = NULL;
   long npoints = 0;

   for (i = 3; i < n; i++) {
      /*
       * build the coefficients ax, ay, bx and by, using:
       *                             _              _   _    _
       *   i                 i    1 | -1   3  -3   1 | | Pi-3 |
       *  Q (t) = T * M   * G   = - |  2  -5   4  -1 | | Pi-2 |
       *               CR    Bs   2 | -1   0   1   0 | | Pi-1 |
       *                            |_ 0   2   0   0_| |_Pi  _|
       */

      ax = p[i].x - 3 * p[i-1].x + 3 * p[i-2].x - p[i-3].x;
      ay = p[i].y - 3 * p[i-1].y + 3 * p[i-2].y - p[i-3].y;
      bx = 2 * p[i-3].x - 5 * p[i-2].x + 4 * p[i-1].x - p[i].x;
      by = 2 * p[i-3].y - 5 * p[i-2].y + 4 * p[i-1].y - p[i].y;

      /*
       * calculate the forward differences for the function using
       * intervals of size 0.1
       */
#ifndef abs
#define abs(x) ((x)<0?-(x):(x))
#endif
#ifndef max
#define max(x,y) ((x>y)?x:y)
#endif

      { double xsquare, ysquare;
      xsquare = abs(p[i-1].x - p[i-2].x);
      xsquare *= xsquare;
      ysquare = abs(p[i-1].y - p[i-2].y);
      ysquare *= ysquare;
      steps = (int) sqrt(xsquare + ysquare) * 2;
	}

      if (steps+4 > npoints) {
         if (thepoints != NULL) free(thepoints);
	 thepoints = (XPoint *)calloc((steps+4) , sizeof(XPoint));
	 npoints = steps+4;
         }

      stepsize = 1.0/steps;
      stepsize2 = stepsize * stepsize;
      stepsize3 = stepsize * stepsize2;

      x = thepoints[0].x = p[i-2].x;
      y = thepoints[0].y = p[i-2].y;
      dx = (stepsize3*0.5)*ax + (stepsize2*0.5)*bx + (stepsize*0.5)*(p[i-1].x-p[i-3].x);
      dy = (stepsize3*0.5)*ay + (stepsize2*0.5)*by + (stepsize*0.5)*(p[i-1].y-p[i-3].y);
      d2x = (stepsize3*3) * ax + stepsize2 * bx;
      d2y = (stepsize3*3) * ay + stepsize2 * by;
      d3x = (stepsize3*3) * ax;
      d3y = (stepsize3*3) * ay;

      /* calculate the points for drawing the curve */

      int jj = 1;
      for (j = 0; j < steps; j++) {
	 x = x + dx;
	 y = y + dy;
	 dx = dx + d2x;
	 dy = dy + d2y;
	 d2x = d2x + d3x;
	 d2y = d2y + d3y;
	 if ((thepoints[jj-1].x != (int)x) ||(thepoints[jj-1].y != (int)y)) {
         thepoints[jj].x = (int)x;
         thepoints[jj++].y = (int)y;
	 }
         }
      helper(w, thepoints, jj);
      }
   if (thepoints != NULL) {
      free(thepoints);
      thepoints = NULL;
      }
   }

static void curveHelper(Window w, XPoint *thepoints, int n)
   {
      int i;
      for(i=0; i<n; i++) {
	 buf[thepoints[i].y * 600 + thepoints[i].x] = 1;
	 }
   /*
    * Could use drawpoints(w, thepoints, n)
    *  but that ignores the linewidth and linestyle attributes...
    * Might make linestyle work a little better by "compressing" straight
    *  sections produced by genCurve into single drawline points.
    */
   drawpts(w, thepoints, n);
   }

/*
 * draw a smooth curve through the array of points
 */
void MydrawCurve(Window w, XPoint *p, int n)
   {
   nRealPoints = 0;
   maxRealPoints = 72000;
   realPoints = calloc( 72000, sizeof(XPoint));
   printf("gencurve...\n"); fflush(stdout);
   MygenCurve(w, p, n, curveHelper);
   printf("floodfill reverse...\n"); fflush(stdout);
   floodfill_4(0,0);
   printf("floodfill done...\n"); fflush(stdout);
   populaterealpoints();
   printf("realpoints done...\n"); fflush(stdout);
   drawpts(w, realPoints, nRealPoints);
   }