/*********************************************************************/
/*                                                                   */
/*                     Carnegie Mellon University                    */
/*                         Geometric Modeling                        */
/*                            Spring 1998                            */
/*                                                                   */
/*                        Prof.  Kenji Shimada                       */
/*                                                                   */
/*                       Student : Nidal Jaalouk                     */
/*                                                                   */
/*             *******************************************           */
/*                                                                   */
/*                 Topic :  Regular 3D Objects Packing               */
/*                                                                   */
/*                          File : packing.c                         */
/*                                                                   */
/*********************************************************************/


#include<stdio.h>
#include<stdlib.h>
#include<math.h>

void main()
{
  FILE *file;
  FILE *out;
  FILE *temp;
  FILE *get;
  char f[50], g[50], index[20], point[20];
  float x[10000], y[10000], z[10000], facearea[500];
  int p[500][500];
  int test1=0, test2=0, test3, test4, times;
  int i=1, j, k, num, lind;
  float vxi, vyi, vzi, vx, vy, vz, h, d, add, xp, yp, zp;
  int u ,v, count, all;
  float tx,ty,tz,s1,s2,sign;
  float ar, maxarea=0.0, minarea=1000000.0, edge, maxedge;
  float x0,y0,z0,x1,y1,z1,x2,y2,z2,x3,y3,z3;
  int origin, ind1, ind2, maxface, minface, countfaces, countpoints, ed1, ed2;
  float xmin=1000000.0, xmax=-1000000.0, ymin=1000000.0, ymax=-1000000.0, zmin=1000000.0, zmax=-1000000.0;
  float dim1, dim2, dim3, dimension1, dimension2, dimension3, volume = 1000000.0;
  float deti, vxii, vyii, vzii, detii, vxiii, vyiii, vziii, detiii;
  int er, three, plane, countplane, baseface;
  float area(float, float, float, float, float, float, float, float, float);
  
  
  
  
  printf("\n\nEnter the name of the file to read data\n");
  scanf("%s", &f);
  file = fopen(f, "r");
  if (file == NULL)
    {
      printf("Error in opening the file\n");
      exit(0);
    }
  
  
  printf("\n\nEnter a name of a file to save results\n");
  scanf("%s", &g);
  out = fopen(g, "w");
  if (out == NULL)
    {
      printf("Error in opening the file\n");
      exit(0);
    }
  
  
  printf("\n\nEnter tolerance for box dimensions.\n");
  scanf("%f", &h);
  
  
  
  /******************************************************************/
  /*         Scanning and saving data from the input file           */
  /******************************************************************/
  
  fscanf(file,"%s",&index);
  while(!test1)
    {
      if (index[0] == 'v')
	{
	  fscanf(file,"%f", &x[i]);
	  if (x[i] > xmax)
	    xmax = x[i];
	  if (x[i] < xmin)
	    xmin = x[i];
	  fscanf(file,"%f", &y[i]);
	  if (y[i] > ymax)
	    ymax = y[i];
	  if (y[i] < ymin)
	    ymin = y[i];
	  fscanf(file,"%f", &z[i]);
	  if (z[i] > zmax)
	    zmax = z[i];
	  if (z[i] < zmin)
	    zmin = z[i];
	  fscanf(file,"%s", &index);
	  
	  i++;
	  
	  if (index[0] != 'v')
	    test1 = 1;
	}
    }
  
  
  
  
  
  temp = fopen("temp.dat", "w");
  times = 1;
  while (!test2)
    {
      if (fscanf(file,"%s",&point) == 1)
	{
	  if (point[0] != 'f')
	    fprintf(temp,"%s ",point);
	  else if (times != 1)
	    fprintf(temp,"-1\n");
	}
      else
	test2 = 1;
      
      times = 2;
    }
  fclose(temp);
  
  
  get = fopen("temp.dat", "r");
  test2 = 0;
  i = 1;
  j = 1;
  countfaces = 0;
  countpoints = 0;
  while(!test2)
    {
      if (fscanf(get,"%d", &num) == 1)
	{
	  if (num != -1)
	    {
	      p[i][j] = num + 1;	      
	      countpoints++;
	      j++;
	    }
	  else
	    {
	      countfaces++;
	      p[i][0] = countpoints + 1;
	      p[i][j] = p[i][1];
	      countpoints = 0;
	      i++;
	      j = 1;
	    }
	}
      else
	test2 = 1;
    }
  
  countfaces++;
  p[i][0] = countpoints + 1;
  p[i][j] = p[i][1];


  /******************************************************************/
  /* Finding the areas of all faces and the maximum and minimum     */
  /******************************************************************/

  for (i=1; i<=countfaces; i++)
    {
      ar = 0.0;
      for (j=1; j<=(p[i][0] - 3); j++)
	{
	  origin = p[i][1];
	  x0 = x[origin];
	  y0 = y[origin];
	  z0 = z[origin];
	  ind1 = p[i][j+1];
	  x1 = x[ind1];
	  y1 = y[ind1];
	  z1 = z[ind1];
	  ind2 = p[i][j+2];
	  x2 = x[ind2];
	  y2 = y[ind2];
	  z2 = z[ind2];
	  
	  ar = ar + area(x0, y0, z0, x1, y1, z1, x2, y2, z2);
	}

      facearea[i] = ar;
      
      if (ar > maxarea)
	{
	  maxarea = ar;
	  maxface = i;
	}

      if (ar < minarea)
	{
	  minarea = ar;
	  minface = i;
	}
    } 
  

  /******************************************************************/
  /* Looping through all faces to find the box of the minimum size  */
  /******************************************************************/
  
  for (all=1; all<=countfaces; all++)
    {
      three = 0;
      plane = 0;
      countplane = 0;
      maxedge = 0.0;
      dim1 = 0.0;
      dim2 = 0.0;
      dim3 = 0.0;
      
      /******************************************************************/
      /* Looping through all edges of a face to find the longest edge   */
      /******************************************************************/

      for (j=1; j<=(p[all][0] - 1); j++)
	{
	  ind1 = p[all][j];
	  x1 = x[ind1];
	  y1 = y[ind1];
	  z1 = z[ind1];
	  ind2 = p[all][j+1];
	  x2 = x[ind2];
	  y2 = y[ind2];
	  z2 = z[ind2];
	  
	  edge = sqrt((x1-x2)*(x1-x2) + (y1-y2)*(y1-y2) + (z1-z2)*(z1-z2));
	  
	  if (edge > maxedge)
	    {
	      ed1 = ind1;
	      ed2 = ind2;
	      lind = j+1;
	    }
	}

      ind1 = ed1;
      x1 = x[ind1];
      y1 = y[ind1];
      z1 = z[ind1];
      
      
      ind2 = ed2;
      x2 = x[ind2];
      y2 = y[ind2];
      z2 = z[ind2];
      
      if (ed2 == p[all][(p[all][0])])
  	{
	  ind1 = p[all][2];
        }
      else
	{
	  ind1 = p[all][lind + 1];
        }
      x3 = x[ind1];
      y3 = y[ind1];
      z3 = z[ind1];
      
      
      /**********************************************************************/
      /* Finding the normal vectors of the three initial planes of the box  */
      /**********************************************************************/

      vxi = (y1-y2)*(z3-z2) - (y3-y2)*(z1-z2);
      vyi = (x3-x2)*(z1-z2) - (x1-x2)*(z3-z2);
      vzi = (x1-x2)*(y3-y2) - (x3-x2)*(y1-y2);
      
      deti = sqrt(vxi*vxi + vyi*vyi + vzi*vzi);
      
      
      x3 = vxi;
      y3 = vyi;
      z3 = vzi;
      
      vxii = (y1-y2)*(z3) - (y3)*(z1-z2);
      vyii = (x3)*(z1-z2) - (x1-x2)*(z3);
      vzii = (x1-x2)*(y3) - (x3)*(y1-y2);
      
      detii = sqrt(vxii*vxii + vyii*vyii + vzii*vzii);
      
      x1 = vxi;
      y1 = vyi;
      z1 = vzi;
      
      x3 = vxii;
      y3 = vyii;
      z3 = vzii;
      
      vxiii = (y1)*(z3) - (y3)*(z1);
      vyiii = (x3)*(z1) - (x1)*(z3);
      vziii = (x1)*(y3) - (x3)*(y1);
      
      detiii = sqrt(vxiii*vxiii + vyiii*vyiii + vziii*vziii);
      
     
      /******************************************************************/
      /* Loop through the three planes. For each plane, use the slicer  */
      /* concept to find the boundaries of the object in the direction  */
      /* of the normal vector of the plane.
      /******************************************************************/
 
      while (!three)
	{
	  countplane++;
	  plane++;
	  test4 = 0;
	  
	  if (plane ==1)
	    {
	      vx = vxi/deti;
	      vy = vyi/deti;
	      vz = vzi/deti;
	    }
	  else if (plane == 2)
	    {
	      vx = vxii/detii;
	      vy = vyii/detii;
	      vz = vzii/detii;
	    }
	  else
	    {
	      vx = vxiii/detiii;
	      vy = vyiii/detiii;
	      vz = vziii/detiii;
	    }


	  /******************************************************************/
	  /*     Outer while will be executed twice:                        */
	  /*   1- slice plane begins at an arbitrary point and slides (+h)  */
	  /*   2- slice plane begins at the previous point and slides (-h)  */
	  /*                                                                */
	  /*   -  Inner while to slide the plane.                           */
	  /*   -  Inner loop to find intersections between lines forming    */
	  /*      the faces and the slice plane.                            */
	  /******************************************************************/
	  
	  sign = 1.0;
	  while(!test4)
	    {
	      test3 = 0;
	      
	      xp = x2;
	      yp = y2;
	      zp = z2;
	      
	      i = 0;
	      
	      while(!test3)
		{
		  d = vx*xp + vy*yp + vz*zp + sign*i*h;
		  
		  count = 0;
		  
		  for(k=1; k<=countfaces; k++)
		    {
		      for (er=1; er<=(p[k][0] - 1); er++)
			{
			  u=p[k][er];
			  v=p[k][er + 1];
			  
			  s1 = vx*x[u] + vy*y[u] + vz*z[u];
			  s2 = vx*x[v] + vy*y[v] + vz*z[v] - s1;
			  
			  if (s2 != 0.0)
			    {
			      tx = x[u] + (d-s1)*(x[v]-x[u])/s2;
			      ty = y[u] + (d-s1)*(y[v]-y[u])/s2;
			      tz = z[u] + (d-s1)*(z[v]-z[u])/s2;
			      /******************************************************************/
			      /*  The condition is to determine if the intersection between     */
			      /*  the plane and the line is valid or not:                       */
			      /*    - if the intersection point lies between the two ends of    */
			      /*      the line  =>  valid.                                      */
			      /*    - otherwise =>  not valid.                                  */
			      /******************************************************************/
			      if ( ( sqrt ( (tx-x[u])*(tx-x[u]) + (ty-y[u])*(ty-y[u]) + (tz-z[u])*(tz-z[u])) + sqrt ( (tx-x[v])*(tx-x[v]) + (ty-y[v])*(ty-y[v]) + (tz-z[v])*(tz-z[v])) ) <= 1.0001 * sqrt ( (x[v]-x[u])*(x[v]-x[u]) + (y[v]-y[u])*(y[v]-y[u]) + (z[v]-z[u])*(z[v]-z[u])) )
				{
				  count++;
				}
			    }
			}
		    }
		  
		  if (count == 0)
		    {
		      test3 = 1;
		    }
		  
		  i++;
		}
	      
	      /*******************************************************************/
	      /* Adding the steps of the slicer to detrmine the dimension of the */
	      /* box in each of the three perpendicular directions               */
	      /*******************************************************************/
	      
	      if (plane == 1)
		dim1 = dim1 + (i-1)*h;
	      else if (plane == 2)
		dim2 = dim2 + (i-1)*h;
	      else
		dim3 = dim3 + (i-1)*h;
	      
	      sign = sign - 2.0;
	      if (sign != -1.0)
		test4 = 1;
	    }

	  
	  if (countplane == 3)
	    three = 1;

	}


      if ((dim1*dim2*dim3) < volume)
	{
	  dimension1 = dim1;
	  dimension2 = dim2;
	  dimension3 = dim3;
	  baseface = all;
	  volume = dim1*dim2*dim3;
	}


    }

  fprintf(out,"\n");
  fprintf(out,"***************************************************************************\n");
  fprintf(out,"*                                                                         *\n");
  fprintf(out,"*                         Initial Information                             *\n");
  fprintf(out,"*                                                                         *\n");
  fprintf(out,"***************************************************************************\n");
        
  fprintf(out,"\nThe boundaries of the object are :\n  xmin = %f,    xmax = %f \n  ymin = %f,    ymax = %f \n  zmin = %f,    zmax = %f\n",xmin,xmax,ymin,ymax,zmin,zmax);
  fprintf(out,"\nThe number of faces in the object = %d\n",countfaces);
  fprintf(out,"\nThe area of the largest face of the object = %f    \n The number of the largest face of the object = %d\n",maxarea, maxface);
  fprintf(out,"\nThe area of the smallest face of the object = %f    \n The number of the smallest face of the object = %d\n",minarea, minface);

  fprintf(out,"\n\n\n");
  fprintf(out,"***************************************************************************\n");
  fprintf(out,"*                                                                         *\n");
  fprintf(out,"*                               Results                                   *\n");
  fprintf(out,"*                                                                         *\n");
  fprintf(out,"***************************************************************************\n");
  
  fprintf(out,"\nDimensions of the box :\n  Length = %f \n  Width  = %f \n  Height = %f\n",dimension2,dimension3,dimension1);
  
  fprintf(out,"\nFurther information about the face of the object\n");
  fprintf(out,"that will rest on the box base (%f * %f) :\n",dimension2,dimension3);
  fprintf(out,"\nThe area of the face = %f\n",facearea[baseface]);
  fprintf(out,"\nThe number of the face = %d\n\n\n",baseface);
  
  
}





/********************************************************************/
/* This function calculats and returns the area of a triangle in 3D */
/********************************************************************/ 
  
  float area(float x0, float y0, float z0, float x1, float y1, float z1, float x2, float y2, float z2)
    {
      float xx, yy, zz, aa;
      
      xx = (y1-y0)*(z2-z0) - (z1-z0)*(y2-y0);
      yy = (z1-z0)*(x2-x0) - (x1-x0)*(z2-z0);
      zz = (x1-x0)*(y2-y0) - (y1-y0)*(x2-x0);
      
      aa = 0.5 * sqrt( xx*xx + yy*yy + zz*zz );
      
      return(aa);
    }