/* *********************************************************************
*
 *
 *  A Simple Model for Hiding Surfaces
 *  Jay Martin Anderson
 *  Franklin & Marshall College, Lancaster, Pennsylvania
 *  Tymlabs Corporation, Austin, Texas
 *
 *  This program draws a three-dimensional object described by a table
 *  of faces using hidden-surface elimination.
 *  The algorithm is described briefly in M. Berger, "Computer Graphics
 *  with Pascal," (Benjamin/Cummings, 1986), and in somewhat greater
 *  detail by R. J. Rost, reported in Creative Computing, February, 1984
.
 *
 *  Implemented for HP CRTs on HP 3000, using Pascal and AGL (HP's graph
ics
 *     language, April, 1984; revised, April, 1986; March, 1988; revised
 *     to use Tymlabs' C/3000 for the HP 3000, May, 1988.
 *
 *  Implemented May, 1988, for the Apple Macintosh (monochrome), using
 *     Lightspeed C and Quickdraw graphics (Macintosh toolbox).
 *
 * ********************************************************************
*/

 /* ********************* #INCLUDEd files *****************************
*/
 #include "stdio.h"
  /* Use the "stdio" window provided by Lightspeed C, which is
   * 500 pixels wide by 288 pixels tall, rather than using window
   * management from Mac toolbox directly.
   */
 #include "unix.h"
 #include "Quickdraw.h"
 #include "math.h"

 /* ********************* Some Common DEFINEs *************************
*/
 #define MAXFACES 60     /* maximum number of faces */
 #define MAXPTS 100                    /* maximum number of vertices */
 #define NVF 4       /* vertices/face; squares  */
 #define HUGE_VAL 1.0E77    /* a very large number */
 #define TRUE  1
 #define FALSE 0
 #define sysPatListID 0     /* Resource ID, patterns */
 typedef float MATRIX[4][4];            /* 4x4 matrix of real numers */
 struct faceS       /* what we know about a face */
 {
  float z;    /* smallest Z, eye coords, for face */
  int color;    /* "color" of this face */
  int v[NVF];    /* list of vertices for this face */
 };

 /* ********************* Global Variables ****************************
*/
 int v[MAXFACES][NVF];      /* vertex list */
 float x0[MAXPTS], y0[MAXPTS], z0[MAXPTS]; /* original coordinates */
 float xt[MAXPTS], yt[MAXPTS], zt[MAXPTS]; /* transformed coordinates */
 float xs[MAXPTS], ys[MAXPTS];    /* screen coordinates */
 struct faceS faceList[MAXFACES];   /* list of faces */
 MATRIX view;        /* viewing transformation */
 float eyex, eyey, eyez;     /* position of eye */
 float fx, fy, fz;       /* where eye is focussing */
 float ds;         /* scale factor */
 float horizRotn;       /* rotation of horizon */
 float n1, n2, n3;       /* normal to a face */
 int numFaces, numVertices, numVisFaces;

 int readFaceData()
 {
  /* Hardcoded for one cube of size 2, centered at origin;
   * Example 1.
   */
  /* for each vertex, its coordinates */
  x0[0] =  1.0, y0[0] =  1.0, z0[0] =  1.0;
  x0[1] = -1.0, y0[1] =  1.0, z0[1] =  1.0;
  x0[2] = -1.0, y0[2] = -1.0, z0[2] =  1.0;
  x0[3] =  1.0, y0[3] = -1.0, z0[3] =  1.0;
  x0[4] =  1.0, y0[4] =  1.0, z0[4] = -1.0;
  x0[5] = -1.0, y0[5] =  1.0, z0[5] = -1.0;
  x0[6] = -1.0, y0[6] = -1.0, z0[6] = -1.0;
  x0[7] =  1.0, y0[7] = -1.0, z0[7] = -1.0;
  /* for each face, its vertex list in CCW order */
  v[0][0] = 0, v[0][1] = 1, v[0][2] = 2, v[0][3] = 3; /* top */
  v[1][0] = 1, v[1][1] = 5, v[1][2] = 6, v[1][3] = 2; /* rear */
  v[2][0] = 3, v[2][1] = 2, v[2][2] = 6, v[2][3] = 7; /* left */
  v[3][0] = 0, v[3][1] = 3, v[3][2] = 7, v[3][3] = 4; /* front */
  v[4][0] = 1, v[4][1] = 0, v[4][2] = 4, v[4][3] = 5; /* right */
  v[5][0] = 4, v[5][1] = 7, v[5][2] = 6, v[5][3] = 5; /* bottom */

  return(TRUE);
 }

 #define COS1 0.540302
 #define SIN1 0.841471
 int getParams()
 {
  /* Hardcoded for a sequence of views of the one cube which is
   * hardcoded in "readDataFile."  Values are initialized in "initStuff"
.
   */

   float tempx, tempy;

   /* Position of eye */
   tempx = eyex;
   tempy = eyey;

   /* Spiral the viewer, or eye, around the object; at each
    * iteration, rotation by 1 radian, and move up by 1/2 unit.
    * To save time, use defined constants for cos(1) and sin(1).
    */
   eyex = COS1*tempx - SIN1*tempy;
   eyey = SIN1*tempx + COS1*tempy;
   eyez += 0.5;

  /* Focus point; where eye is looking */
   fx = fy = fz = 0.0;

   if (eyez > 10.0) return(FALSE);  /* that's enough */
   else return(TRUE);
 }

 void initMat(M)
 MATRIX M;
 {
  /* Construct a 4x4 identity matrix */
  int i, j;       /* subscript; loop index */

  for (i = 0; i < 4; i++)
   for (j = 0; j < 4; j++)
    M[i][j] = (i != j) ? 0.0 : 1.0;
 }

 void multMat(M3, M1, M2)
 MATRIX M1, M2, M3;
 {
  /* Multiply M1 x M2 and put result in M3; 4x4 square matrices */
  int i, j, k;      /* subscript; loop index */

  for (i = 0; i < 4; i++)
   for (j = 0; j < 4; j++)
   {
    M3[i][j] = 0.0;
    for (k = 0; k < 4; k++) M3[i][j] += M1[i][k]*M2[k][j];
   }
 }

 void calcViewMatrix()
 {
  /* Calculate the viewing transformation matrix */
  int i, j;     /* matrix subscripts; loop index */
  MATRIX t1, t2;    /* temp matrices for transformation */
  double d1, d2;    /* temporary results */

  initMat(view);

  /* translate origin to eye position */
  view[3][0] = -eyex;
  view[3][1] = -eyey;
  view[3][2] = -eyez;
  initMat(t1);

  /* rotate about x-axis by 90 degrees */
  t1[1][1] = 0.0;
  t1[2][2] = 0.0;
  t1[1][2] = -1.0;
  t1[2][1] = 1.0;
  multMat(t2, view, t1);
  for (i = 0; i < 4; i++)
   for (j = 0; j < 4; j++) view[i][j] = t2[i][j];
 initMat(t1);

  /* rotate about y-axis by an angle dependent on focus point */
  fx = eyex - fx;
  fy = eyey - fy;
  fz = eyez - fz;
  d1 = sqrt((double)(fx*fx + fy*fy));
  if (fabs(d1) > 0.0001)
  {
   t1[0][0] = -fy/d1;
   t1[2][2] = -fy/d1;
   t1[0][2] =  fx/d1;
   t1[2][0] = -fx/d1;
   multMat(t2, view, t1);
   for (i = 0; i < 4; i++)
    for (j = 0; j < 4; j++) view[i][j] = t2[i][j];
  }
  initMat(t1);

  /* rotate about x-axis by angle dependent on focus point */
  d2 = sqrt((double)(fx*fx + fy*fy  + fz*fz));
  if (fabs(d2) > 0.0001)
  {
   t1[1][1] = d1/d2;
   t1[2][2] = d1/d2;
   t1[1][2] = fz/d2;
   t1[2][1] = -fz/d2;
   multMat(t2, view, t1);
   for (i = 0; i < 4; i++)
    for (j = 0; j < 4; j++) view[i][j] = t2[i][j];
  }
  initMat(t1);

  /* rotate about z-axis to rotate horizon */
  horizRotn *= PI/180.0; /* convert degrees to radians */
  t1[0][0] = t1[1][1] = (float)cos((double)horizRotn);
  t1[0][1] = (float)sin((double)horizRotn);
  t1[1][0] = -t1[0][1];
  multMat(t2, view, t1);
  for (i = 0; i < 4; i++)
   for (j = 0; j < 4; j++) view[i][j] = t2[i][j];
  initMat(t1);

  /* invert the z-axis */
  t1[2][2] = -1.0;
  /* and scale according to D/S ratio */
  t1[0][0] = ds;
  t1[1][1] = ds;
  multMat(t2, view, t1);
  for (i = 0; i < 4; i++)
   for (j = 0; j < 4; j++) view[i][j] = t2[i][j];
}

void initStuff()
{
 /* Whatever needs to be initialized; may be device-dependent;
  * may be implementation-dependent.  Includes fixed parameters for
  * this example.
  */

  numFaces = 6;
  numVertices = 8;

 /* Write at upper left of console window */
 gotoxy(0,0);
 puts("Example 1.  Cube.");

 /* Initial position of eye */
 eyez = -10.0;
 eyex = 7.50;
 eyey = 13.0;

   /* Horizon rotation */
   horizRotn = 0.0;
   /* Scaling factor */
   ds = 1.5;
}

void finishStuff()
{
 /* Whatever needs to be closed, etc.  May be device-dependent;
  * may be implementation-dependent.
  *
  * Nothing to do in this example.
  */
}

void transform(px, py, pz)
float *px, *py, *pz;
{
 /* Transforms a point, pointed to by px, py, pz, by view transform */

 float temp1, temp2, temp3;
 temp1 = *px, temp2 = *py, temp3 = *pz;
 *px = view[0][0]*temp1+view[1][0]*temp2+view[2][0]*temp3+view[3][0];
 *py = view[0][1]*temp1+view[1][1]*temp2+view[2][1]*temp3+view[3][1];
 *pz = view[0][2]*temp1+view[1][2]*temp2+view[2][2]*temp3+view[3][2];
}

void transformAll()
{
 /* Transforms all points */
 int i;
 for (i = 0; i < numVertices; i++) transform(&xt[i],&yt[i],&zt[i]);
}

void getCrossProduct(nv1, nv2, nv3)
int nv1,nv2,nv3;
{
 /* Computes cross-product of vectors representing two sides of a face;
  * the two sides are the vectors from vertex nv2-nv1 and nv3-nv1.  The
  * cross-product is the outward normal to this face; its three
  * components are stored in the global variables n1, n2, n3.
  */
 float v1, v2, v3, w1, w2, w3;

 v1 = xt[nv2] - xt[nv1];
 v2 = yt[nv2] - yt[nv1];
 v3 = zt[nv2] - zt[nv1];
 w1 = xt[nv3] - xt[nv1];
 w2 = yt[nv3] - yt[nv1];
 w3 = zt[nv3] - zt[nv1];
 n1 = v2*w3 - v3*w2;
 n2 = v3*w1 - v1*w3;
 n3 = v1*w2 - v2*w1;
}

void getDotProduct(pdot, nv1)
float *pdot;      /* the dot product */
int nv1;       /* which vertex */
{
 /* Computes the dot product of the outbound normal from a face
  * with vector from eye to face.
  */
 *pdot = n1*xt[nv1] + n2*yt[nv1] + n3*zt[nv1];
}

void eyeToScreen(x, y, z, px, py)
float x, y, z, *px, *py;
{
 /* Transforms a point from x, y, z eye coordinates to
  * x, y screen coordinates.  This includes the projection
  * and the use of the physical screen size in pixels
  */
 float xC, yC;     /* center of screen */
 float xR, yR;     /* width of screen */

 /* Hardcoded for Lightspeed "stdio window" */
 xC = 250, yC = 144;
 xR = 500, yR = 288;
 *px = xR*(x/z) + xC;
 *py = yR - (yR*(y/z) + yC);
}

int zcompare(pFace1, pFace2)
struct faceS *pFace1, *pFace2;
{
 /* Comparison function used by C library function "qsort" to do
  * sorting.  This function compares the minimum z coordinate of
  * faces, so that faces are sorted in the order of distance from
  * eye.  This is a DESCENDING sort!!
  */
 if (pFace1->z < pFace2->z) return(1);
 else if (pFace1->z > pFace2->z) return(-1);
 else return(0);
}

void sortFaces()
{
 /* Sorts the faceList in descending order of Z for a Z-buffer
  * display; uses quicksort as implemented in library.
  * Requires the preceding function which compares z values for
  * faces.
  */

 qsort(faceList, numVisFaces, sizeof(struct faceS), zcompare);
 }

void drawFace(f)
int f;
{
 /* This procedure draws a face.  It is device- and implementation-
  * dependent.  In this implementation, it uses the Lightspeed C
  * "stdio" or console window, and Macintosh Quickdraw procedures.
  *
  * Each face is a polygon, and is developed with the sequence
  * OpenPoly...MoveTo or LineTo...ClosePoly.  The "color" of a
  * face is interpreted as a QuickDraw pen pattern, and is looked
  * up in the system pattern list.  The system pattern list is
  * presumed to be available.  It consists of 38 patterns.  The
  * first six of these patterns are used for the 6 faces of a
  * cube.  Each face is then "framed," that is, outlined in solid
  * black.
  */
 PolyHandle aFace;
 Pattern thePattern;

 /* Develop the 4-sided polygon which is a face */
 aFace = OpenPoly();
 MoveTo((int)xs[faceList[f].v[0]],(int)ys[faceList[f].v[0]]);
 LineTo((int)xs[faceList[f].v[1]],(int)ys[faceList[f].v[1]]);
 LineTo((int)xs[faceList[f].v[2]],(int)ys[faceList[f].v[2]]);
 LineTo((int)xs[faceList[f].v[3]],(int)ys[faceList[f].v[3]]);
 LineTo((int)xs[faceList[f].v[0]],(int)ys[faceList[f].v[0]]);
 ClosePoly();

 /* "Color" the face */
 GetIndPattern(thePattern, sysPatListID, faceList[f].color);
 PenPat(thePattern);
 PaintPoly(aFace);

 /* Outline or frame the face */
 PenPat(black);
 FramePoly(aFace);

 /* That's all for this face; dispose of the polygon */
 KillPoly(aFace);
}

void screenAll()
{
 /* Projects from eye coordinates to screen coordinates */
 int i;

 for (i = 0; i < numVertices; i++)
 {
  eyeToScreen(xt[i], yt[i], zt[i], &xs[i], &ys[i]);
 }
}

void doFace(f, ff)
int f, ff;
{
 /* Adds a face to the visible face list and computes its z */
 float zmin;
 int i;

 zmin = HUGE_VAL;
 for (i = 0; i < NVF; i++)
 {
  faceList[ff].v[i] = v[f][i];
  if (zt[v[f][i]] < zmin) zmin = zt[v[f][i]];
 }
 faceList[ff].z = zmin;
 faceList[ff].color = f + 1;  /* first six Mac patterns */
}

void copyData()
{
 /* Copies the original data so we can transform it without
  * destroying the original
  */
 int i;

 for (i = 0; i < numVertices; i++)
  xt[i] = x0[i], yt[i] = y0[i], zt[i] = z0[i];
}

void drawPicture()
{
 /* This prcoedure draws the scene, with hidden surfaces removed */
 int i;
 float dot;

 transformAll();     /* Transformation to eye coords */
 screenAll();     /* And to screen coords */
 numVisFaces = 0;
 for (i = 0; i < numFaces; i++)
 {
  /* ******  HERE IS THE HIDDEN SURFACE REMOVAL ALGORITHM!!!  ******
   * Compute the cross-product of any two sides of a face (we use
   * the sides 1-0 and 2-0).  This gets an outbound normal from the
   * face. Take the dot product of the outbound normal with a vector
   * from the eye to the face (we use vector from eye to vertex 0).
   * If the dot product is non-negative, the face is visible.
   * ************************************************************ */
  getCrossProduct(v[i][0], v[i][1], v[i][2]);
  getDotProduct(&dot, i);
  if (dot >= 0.0)
  {
   doFace(i, numVisFaces);
   numVisFaces++;
  }
 }

 /* Arrange the faces in decreasing order of distance-from-eye,
  * so that the last face to be drawn is the one nearest the viewer.
  * This extends the basic algorithm to allow scenes to be made up
  * of more than one solid figure, in some cases.  It is called a
  * z-buffer algorithm, because the faces to be drawn are buffered
  * in order of their "z" value.
  */
 sortFaces();

 /* Erase screen before drawing */
 eraseplot();

 /* ... and draw all visible faces */
 for (i = 0; i < numVisFaces; i++) drawFace(i);
}

 main()
 {
  int ok;
  initStuff();     /* Any and all initializations */

  ok = TRUE;
 do
 {
   ok = readFaceData();  /* Get faces and make safe copy */
   copyData();
   if (ok)
   {
    ok = getParams();  /* Get parameters for this picture */
    if (ok)
    {
     calcViewMatrix(); /* Calculate viewing transformation */
     drawPicture();  /* Draw this picture */
    }
   }
  } while (ok);

  finishStuff();     /* Any cleanup and closing */
 }