// determine perpendicular direction from model node
// determine contour model nodes

// =====================================================================
// PRELIMINARIES

#include <math.h>

#define BACK    0
#define FRONT   1
#define CONTOUR 2
#define UNKNOWN 3

typedef struct vertex_struct {
  float  pos[3];              // object space points
  float  mImage[2];            // 'model' image
  float  iImage[2];           // 'image' image
  float  der[3][2];              // first derivative
  float  der2[3][2];             // second derivative
  float  perpendicular[2];
  int    frontback;
  int    symmetric;
} vertex_td;

typedef struct face_struct {
  float area;
  int   vindex[3];
  int   frontback;
} face_td;

// =====================================================================
// OBJECT
// vertex list and face list
// each face is triangle with 3 indices into vertex list

int  numFaces = 12;
int  numVerts = 8;
int xres=20,yres=20;       // x and y res of image

static vertex_td   varray[8];
static face_td     farray[12];
float  mTransform[3];          // model translation transformation
float  iTransform[3];          // image translation transformation

static int image[200][200];        // zero or one for now

int  debug0 = 0;
int  debug1 = 1;


void init_image(void);
void generate_images(void);
void vertex_get_perpendicular(int,float *,float *);
float face_get_area(int);

// =====================================================================
// =====================================================================
// IMAGE ROUTINES

int main()
{
  init_image();
  generate_images();
}


// =================================================================
// INIT IMAGE
void init_image()
{
  int  i,j;


  iTransform[0] = 0.0;
  iTransform[1] = 0.0;
  iTransform[2] = 4.0;

  mTransform[0] = 1.0;
  mTransform[1] = 0.0;
  mTransform[2] = 5.0;

  varray[0].pos[0] =  1.0;
  varray[0].pos[1] =  1.0;
  varray[0].pos[2] =  1.0;
  varray[0].symmetric = 2;

  varray[1].pos[0] =  1.0;
  varray[1].pos[1] =  1.0;
  varray[1].pos[2] = -1.0;
  varray[1].symmetric = 1;

  varray[2].pos[0] = -1.0;
  varray[2].pos[1] =  1.0;
  varray[2].pos[2] = -1.0;
  varray[2].symmetric = 0;

  varray[3].pos[0] = -1.0;
  varray[3].pos[1] =  1.0;
  varray[3].pos[2] =  1.0;
  varray[3].symmetric = 3;

  varray[4].pos[0] =  1.0;
  varray[4].pos[1] = -1.0;
  varray[4].pos[2] =  1.0;
  varray[4].symmetric = 6;

  varray[5].pos[0] =  1.0;
  varray[5].pos[1] = -1.0;
  varray[5].pos[2] = -1.0;
  varray[5].symmetric = 7;

  varray[6].pos[0] = -1.0;
  varray[6].pos[1] = -1.0;
  varray[6].pos[2] = -1.0;
  varray[6].symmetric = 4;

  varray[7].pos[0] = -1.0;
  varray[7].pos[1] = -1.0;
  varray[7].pos[2] =  1.0;
  varray[7].symmetric = 5;

  // top face
  farray[0].vindex[0] = 0;
  farray[0].vindex[1] = 1;
  farray[0].vindex[2] = 3;

  farray[1].vindex[0] = 1;
  farray[1].vindex[1] = 2;
  farray[1].vindex[2] = 3;

  // x = 1 face
  farray[2].vindex[0] = 0;
  farray[2].vindex[1] = 4;
  farray[2].vindex[2] = 5;

  farray[3].vindex[0] = 5;
  farray[3].vindex[1] = 1;
  farray[3].vindex[2] = 0;

  // z = -1 face
  farray[4].vindex[0] = 1;
  farray[4].vindex[1] = 5;
  farray[4].vindex[2] = 6;

  farray[5].vindex[0] = 6;
  farray[5].vindex[1] = 2;
  farray[5].vindex[2] = 1;

  // x = -1 face
  farray[6].vindex[0] = 2;
  farray[6].vindex[1] = 6;
  farray[6].vindex[2] = 3;

  farray[7].vindex[0] = 3;
  farray[7].vindex[1] = 6;
  farray[7].vindex[2] = 7;

  // z = 1 face
  farray[8].vindex[0] = 0;
  farray[8].vindex[1] = 3;
  farray[8].vindex[2] = 7;

  farray[9].vindex[0] = 0;
  farray[9].vindex[1] = 7;
  farray[9].vindex[2] = 4;
  
  // y = -1 face
  farray[10].vindex[0] = 5;
  farray[10].vindex[1] = 4;
  farray[10].vindex[2] = 7;

  farray[11].vindex[0] = 6;
  farray[11].vindex[1] = 5;
  farray[11].vindex[2] = 7;

  // initialize image array to 0's
  for (i=0; i<xres; i++) 
    for (j=0; j<yres; j++)
      image[i][j] = 0;

}


// =================================================================
// GENERATE IMAGES
void generate_images()
{
  int     i,j;
  float   x,y,z;
  int     i1,i2,i3;
  float   va[2],vb[2];
  float   nz;
  float   xa,ya,za,xb,yb,zb,xc,yc,zc;
  float   x1,y1,x2,y2,x3,y3;
  float   t,xl,xr,xL,xR,dxl,dxr,dy;
  int     xi,yi;
  int     vcwi,vccwi,cwi,ccwi;


  // -------------------------------------------------
  // generate the 'real' image

  if (debug0) printf("iTransform: %f %f %f\n",
		    iTransform[0],iTransform[1],iTransform[2]);

  // generate the image image space vertices
  for (i=0; i<numVerts; i++) {

    if (debug0) printf("%d:",i);

    // transform object to image space
    x = varray[i].pos[0] + iTransform[0];
    y = varray[i].pos[1] + iTransform[1];
    z = varray[i].pos[2] + iTransform[2];

    if (debug0) printf(" %f %f %f:",x,y,z);

    // do perspective
    x = x/z; y = y/z;

    if (debug0) printf(" %f %f:",x,y);

    // map to image space
    varray[i].iImage[0] = xres*(x+1)/2;
    varray[i].iImage[1] = yres*(y+1)/2;


    if (debug0) printf("\n");

  // -------------------------------------------------
  // generate the 'model' image vertices

    // transform object to model-image space
    x = varray[i].pos[0] + mTransform[0];
    y = varray[i].pos[1] + mTransform[1];
    z = varray[i].pos[2] + mTransform[2];

    // do perspective
    x = x/z; y = y/z;

    // map to image space
    varray[i].mImage[0] = xres*(x+1)/2;
    varray[i].mImage[1] = yres*(y+1)/2;

    // compute image space derivative for the model vertices
    varray[i].der[0][0] = 1/z;
    varray[i].der[0][1] = 0.0;
    varray[i].der[1][0] = 0.0;
    varray[i].der[1][1] = 1/z;
    varray[i].der[2][0] = mTransform[0]/(z+1) - mTransform[0]/z;
    varray[i].der[2][1] = mTransform[1]/(z+1) - mTransform[1]/z;
    varray[i].der2[0][0] = 0.0;
    varray[i].der2[0][1] = 0.0;    
    varray[i].der2[1][0] = 0.0;
    varray[i].der2[1][1] = 0.0;    
    varray[i].der2[2][0] = 0.0;
    varray[i].der2[2][1] = 0.0;    
  }

  // --------------------------------------------------
  if (debug0) {
    printf("image vertices\n");
    // dump out vertices
    for (i=0; i<numVerts; i++) {
      printf("%d: %f %f %f\n",i,
               varray[i].pos[0],varray[i].pos[1],varray[i].pos[2]);
      printf("   %f %f\n",varray[i].iImage[0],varray[i].iImage[1]);
    }
  }


  // --------------------------------------------------
  // scan convert the image-image vertices to fill image buffer

  for (i=0; i<numFaces; i++) {

    if (debug0) printf("face %d\n",i);

    i1 = farray[i].vindex[0];
    i2 = farray[i].vindex[1];
    i3 = farray[i].vindex[2];
    x1 = varray[i1].iImage[0];
    y1 = varray[i1].iImage[1];
    x2 = varray[i2].iImage[0];
    y2 = varray[i2].iImage[1];
    x3 = varray[i3].iImage[0];
    y3 = varray[i3].iImage[1];

    // calculate the z coordinate of the normal to see if front face
    va[0] = x1-x2; va[1] = y1-y2; 
    vb[0] = x3-x2; vb[1] = y3-y2;
    nz = vb[0]*va[1] - vb[1]*va[0];

    // determine y-ordering of 3 vertices, call them pa, pb, pc
    if (nz < 0.0) {

      if (debug0) {
	printf(" front face %d\n",i);
        printf(" :(%f,%f) : (%f,%f) : (%f,%f) \n",
              x1,y1,x2,y2,x3,y3);
      }

      if (y1>y2) {    //  1-2
        if (y3>y1) {  //  3-1-2
          xa = x3; ya = y3; xb = x1; yb = y1; xc = x2; yc = y2;
        }
        else {        // 1-2-3 or 1-3-2
          if (y2>y3) { // 1-2-3
            xa = x1; ya = y1; xb = x2; yb = y2; xc = x3; yc = y3;
          }
          else {       // 1-3-2
            xa = x1; ya = y1; xb = x3; yb = y3; xc = x2; yc = y2;
          }
        }
      }
      else {          // 2-1
        if (y3>y2) {  //  3-2-1
          xa = x3; ya = y3; xb = x2; yb = y2; xc = x1; yc = y1;
        }
        else {        // 2-3-1 or 2-1-3
          if (y1>y3) { // 2-1-3
            xa = x2; ya = y2; xb = x1; yb = y1; xc = x3; yc = y3;
          }
          else {       // 2-3-1
            xa = x2; ya = y2; xb = x3; yb = y3; xc = x1; yc = y1;
          }
        }
      }

      // split triangle at yb
      // calculate the left and right x values
      t = xa + (xc-xa)*(ya-yb)/(ya-yc);
      if (t>xb) { xl = xb; xr = t; }
      else { xl = t; xr = xb; }

      if (debug0) {
        printf(" %f %f\n %f:%f %f \n %f %f\n",
              xa,ya,xl,xr,yb,xc,yc);
      }

      // fill in the upper triangle
      dy = ya-yb;
      dxl = (xl-xa)/dy; dxr = (xr-xa)/dy;
      xL = xa; xR = xa;
      yi = (int)(ya);
      xL = xa + dxl*(ya-yi);
      xR = xa + dxr*(ya-yi);
      while (yi > yb) {
        for (xi=(int)(xL+1); xi<=(int)(xR); xi++) {
          image[xi][yi] = 1;
        }
        xL += dxl;
        xR += dxr;
        yi--;
      }

      // fill in the lower triangle
      dy = yb-yc;
      dxl = (xl-xc)/dy; dxr = (xr-xc)/dy;
      xL = xc; xR = xc;
      yi = (int)(yc+1.0);
      xL = xc + dxl*(yi-yc);
      xR = xc + dxr*(yi-yc);
      while (yi < yb) {
        for (xi=(int)(xL+1); xi<=(int)(xR); xi++) {
          image[xi][yi] = 1;
        }
        xL += dxl;
        xR += dxr;
        yi++;
      }

    }
  }

  // --------------------------------------------------
  // process model-image faces

  for (i=0; i<numFaces; i++) {

    if (debug0) printf("face %d\n",i);

    i1 = farray[i].vindex[0];
    i2 = farray[i].vindex[1];
    i3 = farray[i].vindex[2];
    x1 = varray[i1].mImage[0];
    y1 = varray[i1].mImage[1];
    x2 = varray[i2].mImage[0];
    y2 = varray[i2].mImage[1];
    x3 = varray[i3].mImage[0];
    y3 = varray[i3].mImage[1];

    // calculate the z coordinate of the normal to see if front face
    va[0] = x1-x2; va[1] = y1-y2; 
    vb[0] = x3-x2; vb[1] = y3-y2;
    nz = vb[0]*va[1] - vb[1]*va[0];

    // determine y-ordering of 3 vertices, call them pa, pb, pc
    if (nz < 0.0) {
      if (debug0) {
	printf(" front face %d\n",i);
        printf(" :(%f,%f) : (%f,%f) : (%f,%f) \n",
              x1,y1,x2,y2,x3,y3);
      }
      farray[i].frontback = FRONT;
      
      farray[i].area = (x1-x2)*(y1+y2)/2 + 
                       (x2-x3)*(y2+y3)/2 + 
                       (x3-x1)*(y3+y1)/2;
    }
    else { // face is backface
      farray[i].frontback = BACK;
    }
  }

  // ----------------------------------------
  // determine front/back status of vertex

  for (i=0; i<numVerts; i++) {
    varray[i].frontback = UNKNOWN;
    for (j=0; j<numFaces; j++) {
      if ( (i==farray[j].vindex[0]) || 
           (i==farray[j].vindex[1]) || 	
           (i==farray[j].vindex[2])    ) {	
        if (varray[i].frontback == UNKNOWN) {	
          if (farray[j].frontback == FRONT) varray[i].frontback = FRONT;
          else varray[i].frontback = BACK;
        }
        else {
          if (varray[i].frontback == FRONT) {
            if (farray[j].frontback == BACK) varray[i].frontback = CONTOUR;
          }
          else {
            if (farray[j].frontback == FRONT) varray[i].frontback = CONTOUR;
          }
        }		
      }
    }
  }

  // ----------------------------------------
  // determine contour segment around contour vertex
  float  v1[2],v2[2];
  float  len;
  for (i=0; i<numVerts; i++) {
    // if this vertex is a CONTOUR vertex
    if (varray[i].frontback == CONTOUR) {
      // loop through faces to find the faces it's part of
      for (j=0; j<numFaces; j++) {
        if ( (i==farray[j].vindex[0]) ||
             (i==farray[j].vindex[1]) ||
             (i==farray[j].vindex[2])   ) {
          // get the other two vertices in CW and CCW order from this vertex
          if (i==farray[j].vindex[0]) {
            vcwi = farray[j].vindex[1];
            vccwi = farray[j].vindex[2];
          }
          else if (i==farray[j].vindex[1]) {
            vcwi = farray[j].vindex[2];
            vccwi = farray[j].vindex[0];
          }
          else {
            vcwi = farray[j].vindex[0];
            vccwi = farray[j].vindex[1];
	  }
          // if either one is a CONTOUR vertex, record which direction it's in
          if (varray[vcwi].frontback == CONTOUR) cwi = vcwi;
          if (varray[vccwi].frontback == CONTOUR) ccwi = vccwi;
        }
      }
      // normalize the edges first
      // only works with convex shapes right now
      v1[0] = varray[ccwi].mImage[0] - varray[i].mImage[0];
      v1[1] = varray[ccwi].mImage[1] - varray[i].mImage[1];
      len = sqrt(v1[0]*v1[0]+v1[1]*v1[1]);
      v1[0] = v1[0]/len; v1[1] = v1[1]/len;
      v2[0] = varray[cwi].mImage[0] - varray[i].mImage[0];
      v2[1] = varray[cwi].mImage[1] - varray[i].mImage[1];
      len = sqrt(v2[0]*v2[0]+v2[1]*v2[1]);
      v2[0] = v2[0]/len; v2[1] = v2[1]/len;
      varray[i].perpendicular[0] = v1[0] + v2[0];
      varray[i].perpendicular[1] = v1[1] + v2[1];
    }
  }

  // --------------------------------------------------
  // DUMP STUFF OUT

  if (debug1) {
    // dump out image
    for (j=yres-1; j>=0; j--) {
      printf("\n");
      for (i=0; i<xres; i++) {
        printf("%d",image[i][j]);
      }
    }

    printf("\n");
    printf("\n");

    // dump out vertices
    for (i=0; i<numVerts; i++) {
      printf("%d: %f %f %f\n", i,
               varray[i].pos[0],varray[i].pos[1],varray[i].pos[2]);
      printf("   %f %f\n",varray[i].mImage[0],varray[i].mImage[1]);
    }
  }
}


// =================================================================
// VERTEX: GET PERPENDICULAR
void vertex_get_perpendicular(int vi, float *x, float *y)
{
  *x = varray[vi].perpendicular[0];  
  *y = varray[vi].perpendicular[1];  
}

// =================================================================
// FACE: GET AREA
float face_get_area(int fi)
{
  return (farray[fi].area);
}