//Nathan Gill

#include <iostream.h>
#include <fstream.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>

#define x     (0)
#define y     (1)
#define z     (2)
#define min_step  (0.125)
//#define maxover_allowed  (0.2)  USER DEFINED
#define max_points  (150)
#define max_faces   (300)
#define max_points_on_face  (10)
#define default_gap   (1.0)


//Function Prototypes
void find_start_and_end(double& starth, double& endh, double v[][3], int count, double d[]); 
int check_intersect(double v[][3], double h, int p1, int p2, double d[]);
void find_intersection(double d[], double h, double v[][3], int p0, int p1,double& X,double& Y,double& Z, ofstream& outfile);
void check_mismatch(double v[][3],double numfaces,int points_per_face[],double d[],double h, double gap,ofstream& outfile,int f[][max_points_on_face],double tempv[][3],double& maxover, double& maxunder,ofstream& outtext,int num_per_slice);


//BEGIN MAIN PROGRAM
int main()
{ //open input file 
  ifstream infile;
  infile.open("infile.dat");

//open output files
  ofstream outfile;

  outfile.open("outfile.wrl");  //VRML OUTPUT FILE
  ofstream outtext;
  outtext.open("out.txt");      //BACKGROUND INFO OUTPUT FILE

  double d[3];    //This is the direction vector;
  double gap;     //This is the distance between slices
  double v[max_points][3];  //STORES ALL POINT FROM INFILE.DAT
  int count;      //Counts the number of points
  char next;

 double maxover_allowed;

  cout <<"Enter the max allowable difference value: "<<endl;
  cin >> maxover_allowed;

//write VRML header
  outfile << "#VRML V1.0 ascii\n\n";   // write header or VRML file
  outfile << "Material { diffuseColor 1.0 0.0 0.0 }\n "
	  << "Coordinate3 {\n"
          << "    point [\n      ";
  

//read in the vertices
  infile >> next;
  count = 0; //count is the number of times vertices in the original file
  while (next == 'v')
    { infile >> v[count][x];
      infile >> v[count][y];
      infile >> v[count][z];
      infile >> next;
      count++;
    }          
//
cout << "The number of vertices is: "<<count<<endl;
//
int numvert;
numvert = count;

//Get direction vector
  cout << "Enter the X Y Z cordinate of the direction vector:\n";
  cin >> d[x] >> d[y] >> d[z];

//Normalize the direction vector
  double den;
  den = sqrt(d[x]*d[x]+d[y]*d[y]+d[z]*d[z]);
  d[x] = d[x]/den;
  d[y] = d[y]/den;
  d[z] = d[z]/den;
//Note that d is now the unit normal vector

//Find the staring and ending h values
  double starth, endh;

  find_start_and_end(starth, endh, v, count,d);
  cout << "starth and endh are " << starth << " " << endh << endl;


//Store face data 
  int f[max_faces][max_points_on_face]; 
                       //used to store vertice indecies for each face.
  int numfaces;        //used to count the number of faces.
  int points_per_face[max_faces]; //STORES THE NUMBER OF POINTS IN EACH FACE

  numfaces = 1;

  char item[15];
  count = 0;
  int point_count = -1;   //COUNTS THE NUMBER OF POINTS PER FACE

  while (! infile.eof())
    { if (strcmp(item, "f") == 0)
	{ points_per_face[numfaces] = point_count;
	  point_count = 0;
	  numfaces ++;
	}
      if (strcmp(item, "f") != 0)  //if item is not an f
	{ f[numfaces][point_count] = atoi(item);
	  point_count ++;            //line 100
	}
      infile >> item;
    }
  points_per_face[numfaces] = point_count;


  cout << " the number of faces is " << numfaces << endl;






//TEST TO SEE FI FILES WERE READ IN PROPERLY
	ofstream test;
	test.open("test.out");

	for (int i=0;i<numvert;i++)
	{ test << "v " << v[i][x] <<" "<<v[i][y]<<" "<<v[i][z]<<endl;
	}

	for (int i=1; i<=numfaces;i++)
	{ test << "f ";
	  for (int j=0;j<points_per_face[i];j++)
		{ test << f[i][j]; }
	  test << endl;
	}
	test.close();

//END OF TEST

//INTERSECT PLANE WITH PART
  double h;
  h = starth;
  count = 0;   //count will now count the number of intersections

  double maxover = 0;  //MAX OVERHANG OF THE SLICE AND THE ACTUAL PART
  double maxunder = 0; //MAX AMOUNT OF ORIGINAL PART NOT INCLUDED IN SLICE

  while (h <= endh)
   { //The plane will be the following equation
     //  (p * n) = d  

     maxover = 0;
     maxunder = 0;   
//   
     gap = default_gap;
//
     cout << "H: " << h << endl;
     outtext << "H: " << h << endl;
//
     double tempv[max_points][3];  //POINTS OF INTERSECTION
     int num_per_slice = 0;        //NUMBER OF INTERSECTIONS PER SLICE

     for (int i=1; i <=numfaces; i++)  //GO THROUGH ALL THE FACES
      { int check;
        double X,Y,Z;   //Thecoordiates of the intersection point

        for (int j=0; j<points_per_face[i]; j++)  //GO THROUGH ALL THE POINTS ON EACH FACE
          { int p0, p1;  //TWO POINTS DEFINING LINE SEGMENTS
	    
	    p0 = f[i][j];
	    p1 = f[i][j+1];
	    if ((j+1) == points_per_face[i])
	      p1 = f[i][0];
	    check = check_intersect(v, h, p0, p1, d);  //ACHECK TO SEE FI THERE IS INTERSECTION BETWEEN LINE SEGMENT AND SLICING PLANE
	      if (check == 1)
		{ count++;
		  find_intersection(d,h,v,p0,p1,X,Y,Z,outfile); //FIND COORDINATES OF INTERSECTOIN


		  if (X != -9999)
		    { //write intersection point
		      if (count > 1)
		        outfile << ",\n      ";
		      tempv[num_per_slice][x] = X;
		      tempv[num_per_slice][y] = Y;
		      tempv[num_per_slice][z] = Z;
		      num_per_slice ++;
		      outfile << X << " " << Y << " " << Z;
		    }
		  else count--;  
		}
	  }
      }
//ENTER CHECK IF ACCEPTABLE HERE
     
     int cont = 0;

//
//     cont = 1;
//     cout << gap <<endl;
//

     while (cont == 0)
     { maxover = 0;

       check_mismatch(v,numfaces,points_per_face,d,h,gap,outfile,f,tempv,maxover,maxunder,outtext,num_per_slice);
       cont = 1;
       if (maxover > maxover_allowed)
         cont = 0;
       //if (maxunder > maxunder_allowed)
         //cont = 0;
       if (cont == 0)
         gap = gap * 0.5;
       if (gap < min_step)
         { cont = 1;
           gap = min_step;
	   }
       }
//
  cout << "   GAP: " << gap << endl;   
  outtext << "   GAP: " << gap <<endl;
//
     h = h + gap;
     
   }
    
  //Write output to close point list
  outfile << " ]\n}\n";

  int numsegments;   //This is the number of line segments to be drawn.
  numsegments = count;
  cout  << "number of segments: " << numsegments << endl;

//Write IndexedLineSet data
  outfile << "IndexedLineSet {\n    coordIndex [\n          ";

  for (int loop=0;loop<numsegments;loop++)
   { if (loop > 0)                        //LINE 200
       outfile << ",\n          ";
     outfile << loop << ",";
     loop++;
     outfile << loop << ",-1";
   }
  outfile << " ]\n}";  


  //close the files
  infile.close();
  outfile.close();
  outtext.close();

}

//END MAIN PROGRAM
//******************************************************
//Functions

void find_start_and_end(double& starth, double& endh, double v[][3], int count, double d[])
{  double temp, max, min;
   temp = 0.0;
   max = 0.0;
   min = 10000.0;
   
   for (int i=0; i<=count; i++)
     { temp = v[i][x]*d[x] + v[i][y]*d[y] + v[i][z]*d[z];
       if (temp > max) 
         max = temp;
       if (temp < min) 
         min = temp;
     }
   starth = min - min_step; //The min step is to make sure the intersection 
   endh = max + min_step;   //isn't a single point
}

//******************************************************************

int check_intersect(double v[][3], double h, int p1, int p2, double d[])
{  //This procedure checks for intersections between the slicing
   //plane and a line segment.  
   //Note: it does not deal with the case where either of the
   //points defining the line segment lie on the plane.

   double a,b;
   a = v[p1][x]*d[x]+v[p1][y]*d[y]+v[p1][z]*d[z];
   b = v[p2][x]*d[x]+v[p2][y]*d[y]+v[p2][z]*d[z];

   if (a < h)
    { if (b < h)
       return 0;
      else
       return 1;
    }
   else
    { if (b > h)
       return 0;
      else
       return 1;
    }
}
   
//****************************************************************

void find_intersection(double d[], double h, double v[][3], int p0, int p1,double& X,double& Y,double& Z, ofstream& outfile)
{ //THIS FINDS THE INTERSECTION BETWEEN A PLANE AND A LINE SEGMENT
  /* PLANE:  (p*d)=h
     LINE:   p(t) = p0 + tv;   v = p1 - p0
  */

  double x0,y0,z0;  //FIRST POINT OF HTE LINE SEGMENT
  double vx,vy,vz;  //DIRECTION VECTOR OF LINE SEGMENT
  double t;         //PARAMETER VALUE OF LINE SEGMENT
  
  x0 = v[p0][x];
  y0 = v[p0][y];
  z0 = v[p0][z];
  
  vx = v[p1][x] - v[p0][x];
  vy = v[p1][y] - v[p0][y];
  vz = v[p1][z] - v[p0][z];
  
  double temp1, temp2;
  temp1 = d[x]*x0 + d[y]*y0 + d[z]*z0;
  temp1 = h - temp1;
  temp2 = d[x]*vx + d[y]*vy + d[z]*vz;
  t = temp1/temp2;   //THIS IS THE PARAMETER VALUE


  if ((vx==0)&&(vy==0)&&(vz==0))
    X = -9999;
  else
    { X = x0 + t*vx;
      Y = y0 + t*vy;   
      Z = z0 + t*vz; 
    }
}

//****************************************************************





void check_mismatch(double v[][3],double numfaces,int points_per_face[],double d[],double h,double gap,ofstream& outfile,int f[][max_points_on_face],double tempv[][3],double& maxover,double& maxunder,ofstream& outtext,int num_per_slice)
{ //THIS FINDS THE DIFFERENCES BETWEEN THE ACTUAL PART AND THE SLICES
  //TO DETERMINE IS THE CURRENT GAP VALUE IS ACCEPTABLE.

  double r[3], point[3];  //THE POINT TO BE CHECKED AGAINST SLICES
  double X,Y,Z; 
  double dist = 999999; //THE DISTANCE FROM THE POINT TO THE LINE SEGMENT
  maxover = 0;

//
	outtext << "CHECKING NEXT PLANE"<<endl;
//

  for (int i=1; i <=numfaces; i++)
      { for (int j=0; j<points_per_face[i]; j++)
          { int p0, p1;  //TWO POINTS DEFINING LINE SEGMENTS
	    
	    p0 = f[i][j];
	    p1 = f[i][j+1];
	    if ((j+1) == points_per_face[i])
	      p1 = f[i][0];
	    int check = 0;
	    check = check_intersect(v, h+gap, p0, p1, d);
	    if (check == 1)
		{ find_intersection(d,h+gap,v,p0,p1,X,Y,Z,outfile);
	      point[x] = X;
		point[y] = Y;
		point[z] = Z;	      
		double n[3];  //This is p1 - p0;
//
//	outtext <<"p0: "<<
		n[x] = tempv[p1][x] - tempv[p0][x];
		n[y] = tempv[p1][y] - tempv[p0][z];
		n[z] = tempv[p1][z] - tempv[p0][z];

		r[x] = point[x];// - (gap)*d[x];
		r[y] = point[y];// - (gap)*d[y];
		r[z] = point[z];// - (gap)*d[z];

		double p[3];  //THIS IS AN INTERMEDIATE VARIABLE
		double length;  //	perpendicular DISTANCE FOMR POINT TO LINE
		double temp_dist;
		double t;    //PARAMETER VALUE OF THE LINE
		double term; //USED IN CALCULATIONS

		term = r[x]*n[x]+r[y]*n[y]+r[z]*n[z];
		t = term;
		term = tempv[p0][x]*n[x]+tempv[p0][y]*n[y]+tempv[p0][z]*n[z];
		t = t - term;
		term = n[x]*n[x]+n[y]*n[y]+n[z]*n[z];
		t = t / term;

		p[x] = tempv[p0][x] + t*n[x];
		p[y] = tempv[p0][y] + t*n[y];
		p[z] = tempv[p0][z] + t*n[z];

		length = sqrt((p[x]-r[x])*(p[x]-r[x])+(p[y]-r[y])*(p[y]-r[y])+(p[z]-r[z])*(p[z]-r[z]));
		temp_dist = length;

		
		if ((t > 1) || (t < 0))
		  {	temp_dist = sqrt(length*length + t*(sqrt(n[x]*n[x]+n[y]*n[y]+n[z]*n[z])));
		  }
            if (temp_dist < dist)
		  dist = temp_dist;


//SOME OUTPUTS
		outtext << " length: "<<length<<endl<<" t: "<<t<<endl<<" temp_dist: "<<temp_dist<<endl;
		outtext << " dist##: "<<dist<<endl;
//
            }


            j++;
	   }
	   if (dist > 9999)
		dist = 0;
	   if (dist > maxover)
            maxover = dist;
	//OUTPUTS
		outtext << " dist: "<<dist<<endl<<" maxover: "<<maxover<<endl;
//

  }
}