/*****************************************************************************/
/*                                                                           */
/*           F I S T :  Fast, Industrial-Strength Triangulation              */
/*                                                                           */
/*****************************************************************************/
/*                                                                           */
/* (C)                            Martin Held                                */
/* (C)            Universitaet Salzburg, Salzburg, Austria                   */
/*                                                                           */
/* This code is not in the public domain. All rights reserved! Please make   */
/* sure to read the full copyright statement contained in api_functions.cpp. */
/*                                                                           */
/*****************************************************************************/
 
/*                                                                           */
/* get standard libraries                                                    */
/*                                                                           */
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <assert.h>

/*                                                                           */
/* get my header files                                                       */
/*                                                                           */
#include "fpkernel.h"
#include "martin.h"
#include "defs.h"
#include "list.h"
#include "header.h"
#include "basic.h"
#include "vertex.h"

/*                                                                           */
/* prototypes of functions provided in this file                             */
/*                                                                           */
void ProjectFace(global_struct *all, int loop_min, int loop_max);


#define BLOCK_SIZE 32768

typedef struct { 
   double element[4][4];
} Matrix4;


/* note:  pnt  and  new_pnt  must not be the same memory cell! */
#define MatrixTransformation(matrix, pnt, new_pnt) \
{ \
  (new_pnt).x = (pnt).x * (matrix).element[0][0] \
              + (pnt).y * (matrix).element[1][0] \
              + (pnt).z * (matrix).element[2][0] \
              +           (matrix).element[3][0]; \
  (new_pnt).y = (pnt).x * (matrix).element[0][1] \
              + (pnt).y * (matrix).element[1][1] \
              + (pnt).z * (matrix).element[2][1] \
              +           (matrix).element[3][1]; \
  (new_pnt).z = (pnt).x * (matrix).element[0][2] \
              + (pnt).y * (matrix).element[1][2] \
              + (pnt).z * (matrix).element[2][2] \
              +           (matrix).element[3][2]; \
}

/*                                                                           */
/* this function projects the vertices of the polygons referenced by         */
/* loops[loop_min,..,loop_max-1] to an approximating plane.                  */
/*                                                                           */
void ProjectFace(global_struct *all, int loop_min, int loop_max)
{
   boolean DetermineNormal(global_struct *all, list_ind ind, vertex *normal);
   void ProjectPoints(global_struct *all, int i1, int i2, vertex normal);
   vertex normal, nr;
   int i;
   double d;

   normal.x = normal.y = normal.z = C_0_0;

   /*                                                                        */
   /* determine the normal of the plane onto which the points get projected  */
   /*                                                                        */
   for (i = loop_min;  i < loop_max;  ++i) {
      if (DetermineNormal(all, all->c_list.loops[i], &nr)) {
         if (DotProduct(normal, nr) < C_0_0) {
            InvertVector(nr);
         }
         VectorAdd(normal, nr, normal);
      }
   }

   if (loop_min <= (loop_max - 1)) {
      d = Length_l2(normal);
      if (gt(d, ZERO)) {
         DivScalar(d, normal);
      }
      else {
         if (all->rt_opt.verbose) {
            FIST_Warning("ProjectFace() - zero-length normal vector!\n");
         }
         normal.x = normal.y = C_0_0; 
         normal.z = C_1_0;
      }
   }

   /* COREBACKEND: Simplify the expression tree of the normal vector,        */
   /* otherwise the final point coordinates' expression tree is sometimes    */
   /* too big for  at least CORE.version <= 2.0.8                            */
#ifdef WITH_COREBACKEND
   normal.x = TO_MDOUBLE(normal.x);
   normal.y = TO_MDOUBLE(normal.y);
   normal.z = TO_MDOUBLE(normal.z);
#endif

   /*                                                                        */
   /* project the points onto this plane. the projected points are stored in */
   /* the array `points[0,..,num_pnts]'                                      */
   /*                                                                        */
   ProjectPoints(all, loop_min, loop_max, normal);

   return;
}


/*                                                                           */
/* this function computes the average of all normals defined by triples of   */
/* successive vertices of the polygon. we'll see whether this is a good      */
/* heuristic for finding a suitable plane normal...                          */
/*                                                                           */
boolean DetermineNormal(global_struct *all, list_ind ind, vertex *normal)
{
	listdef   *list   = &all->c_list;
	vertexdef *vert = &all->c_vertex;

#ifdef SUM_OF_NORMALS
   vertex nr, pq, pr;
   list_ind ind0, ind1, ind2;
   int i0, i1, i2;
   double d;

   ind1 = ind;
   i1   = GetIndex(list, ind1);
   ind0 = GetPrevNode(list, ind1);
   i0   = GetIndex(list, ind0);
   ind2 = GetNextNode(list, ind1);
   i2   = GetIndex(list, ind2);
   VectorSub(vert->vertices[i0], vert->vertices[i1], pq);
   VectorSub(vert->vertices[i2], vert->vertices[i1], pr);
   VectorProduct(pq, pr, nr);
   d = Length_l2(nr);
   if (gt(d, ZERO)) {
      DivScalar(d, nr);
      *normal = nr;
   }
   else {
      normal->x = normal->y = normal->z = C_0_0;
   }
   pq   = pr;
   ind1 = ind2;
   ind2 = GetNextNode(list, ind1);
   i2   = GetIndex(list, ind2);
   while (ind1 != ind) {
      VectorSub(vert->vertices[i2], vert->vertices[i1], pr);
      VectorProduct(pq, pr, nr);
      d = Length_l2(nr);
      if (gt(d, ZERO))  {
         DivScalar(d, nr);
         if (DotProduct(*normal, nr) < C_0_0) {
            InvertVector(nr);
         }
         VectorAdd(*normal, nr , *normal);
      }
      pq   = pr;
      ind1 = ind2;
      ind2 = GetNextNode(list, ind1);
      ind2 = GetIndex(list, ind2);
   }   

   d = Length_l2(*normal);
   if (gt(d, ZERO)) {
      DivScalar(d, *normal);
      return  true;
   }
   else {
      if (all->rt_opt.verbose) {
         FIST_Warning("DetermineNormal() - zero-length normal vector!\n");
      }
      normal->x = normal->y = C_0_0; 
      normal->z = C_1_0;
      return  false;
   }
#else
   vertex nr, pq, pr;
   vertex x_normal, y_normal, z_normal;
   list_ind ind0, ind1, ind2;
   int i0, i1, i2;
   double d;

   x_normal.x = x_normal.y = x_normal.z = C_0_0;
   y_normal.x = y_normal.y = y_normal.z = C_0_0;
   z_normal.x = z_normal.y = z_normal.z = C_0_0;

   ind1 = ind;
   i1   = GetIndex(list, ind1);
   ind0 = GetPrevNode(list, ind1);
   i0   = GetIndex(list, ind0);
   VectorSub(vert->vertices[i0], vert->vertices[i1], pq);
   ind2 = GetNextNode(list, ind1);
   i2   = GetIndex(list, ind2);
   do {
      VectorSub(vert->vertices[i2], vert->vertices[i1], pr);
      VectorProduct(pq, pr, nr);
      d = Length_l2(nr);
      if (gt(d, ZERO))  {
         DivScalar(d, nr);
         if (Abs(nr.x) <= Abs(nr.y)) {
            if (Abs(nr.y) <= Abs(nr.z)) {   
               if (nr.z >= C_0_0)  {VectorAdd(z_normal, nr, z_normal);}
               else                {VectorSub(z_normal, nr, z_normal);}
            }
            else {
               if (nr.y >= C_0_0)  {VectorAdd(y_normal, nr, y_normal);}
               else                {VectorSub(y_normal, nr, y_normal);}
            }
         }
         else {
            if (Abs(nr.x) <= Abs(nr.z)) {   
               if (nr.z >= C_0_0)  {VectorAdd(z_normal, nr, z_normal);}
               else                {VectorSub(z_normal, nr, z_normal);}
            }
            else {
               if (nr.x >= C_0_0)  {VectorAdd(x_normal, nr, x_normal);}
               else                {VectorSub(x_normal, nr, x_normal);}
            }
         }
      }
      pq   = pr;
      ind1 = ind2;
      ind2 = GetNextNode(list, ind1);
      i2   = GetIndex(list, ind2);
   } while (ind1 != ind);

   if (DotProduct(z_normal, x_normal) < C_0_0) {
      InvertVector(x_normal);
   }
   VectorAdd(z_normal, x_normal, z_normal);
   if (DotProduct(z_normal, y_normal) < C_0_0) {
      InvertVector(y_normal);
   }
   VectorAdd(z_normal, y_normal, *normal);
   d = Length_l2(*normal);
   if (gt(d, ZERO)) {
      DivScalar(d, *normal);
      return  true;
   }
   else {
      if (all->rt_opt.verbose) {
         FIST_Warning("DetermineNormal() - zero-length normal vector!\n");
      }
      normal->x = normal->y = C_0_0; 
      normal->z = C_1_0;
      return  false;
   }
#endif
}


/*                                                                           */
/* this function maps the vertices of the polygons referenced by             */
/* loops[i1,..,i2-1] to the plane  n3.x * x + n3.y * y + n3.z * z = 0.       */
/* every mapped vertex  (x,y,z)  is then expressed in terms of  (x',y',z'),  */
/* where  z'=0.  this is achieved by transforming the original vertices into */
/* a coordinate system whose z-axis coincides with  n3,  and whose two other */
/* coordinate axes  n1  and  n2  are orthonormal on  n3. note that n3 is     */
/* supposed to be of unit length!                                            */
/*                                                                           */
void ProjectPoints(global_struct *all, int i1, int i2, vertex n3)
{
   listdef   *list = &all->c_list;
   vertexdef *vert = &all->c_vertex;
   datadef   *data = &all->c_data;

   Matrix4 matrix;
   vertex vtx, n1, n2;
   double d;
   list_ind ind, ind1;
   int i, j1;

   /*                                                                        */
   /* choose  n1  and  n2  appropriately                                     */
   /*                                                                        */
   if ((Abs(n3.x) > C_0_1)  || (Abs(n3.y) > C_0_1)) { 
      n1.x = -n3.y; 
      n1.y =  n3.x; 
      n1.z =  C_0_0; 
   } 
   else { 
      n1.x =  n3.z; 
      n1.z = -n3.x; 
      n1.y =  C_0_0; 
   } 
   d = Length_l2(n1); 
   DivScalar(d, n1); 
   VectorProduct(n1, n3, n2); 
   d = Length_l2(n2); 
   DivScalar(d, n2); 
   
   /*                                                                        */
   /* initialize the transformation matrix                                   */
   /*                                                                        */
   matrix.element[0][0] = n1.x; 
   matrix.element[1][0] = n1.y; 
   matrix.element[2][0] = n1.z; 
   matrix.element[3][0] = C_0_0;     /* translation of the coordinate system */
   matrix.element[0][1] = n2.x; 
   matrix.element[1][1] = n2.y; 
   matrix.element[2][1] = n2.z; 
   matrix.element[3][1] = C_0_0;     /* translation of the coordinate system */
   matrix.element[0][2] = n3.x; 
   matrix.element[1][2] = n3.y; 
   matrix.element[2][2] = n3.z; 
   matrix.element[3][2] = C_0_0;     /* translation of the coordinate system */
   matrix.element[0][3] = C_0_0; 
   matrix.element[1][3] = C_0_0; 
   matrix.element[2][3] = C_0_0; 
   matrix.element[3][3] = C_1_0;      

   /*                                                                        */
   /* transform the vertices and store the transformed vertices in the array */
   /* `points'                                                               */
   /*                                                                        */
  
   // MODIF : alloco solo lo spazio necessario per l'array dei punti ( verrą eventualmente esteso da StorePnt)
   // InitPnts(data, BLOCK_SIZE);
   InitPnts( data, data->max_num_pnts) ;
   for (i = i1;  i < i2;  ++i) {
      ind  = list->loops[i];
      ind1 = ind;
      j1   = GetIndex(list, ind1);
      MatrixTransformation(matrix, vert->vertices[j1], vtx);
#ifdef EXT_APPL_SITES
      j1 = StorePnt(data, vtx.x, vtx.y, GetExtApplVtx(j1));
#else
      j1 = StorePnt(data, vtx.x, vtx.y);
#endif
      UpdateIndex(list, ind1, j1);
      ind1 = GetNextNode(list, ind1);
      j1   = GetIndex(list, ind1);
      while (ind1 != ind) {
         MatrixTransformation(matrix, vert->vertices[j1], vtx);
#ifdef EXT_APPL_SITES
         j1 = StorePnt(data, vtx.x, vtx.y, GetExtApplVtx(j1));
#else
         j1 = StorePnt(data, vtx.x, vtx.y);
#endif
         UpdateIndex(list, ind1, j1);
         ind1 = GetNextNode(list, ind1);
         j1   = GetIndex(list, ind1);
      }   
   }

   return;
}