/*******************************************************************************
+
+  LEDA 3.5
+
+  d3_hull.c
+
+  This file is part of the LEDA research version (LEDA-R) that can be 
+  used free of charge in academic research and teaching. Any commercial
+  use of this software requires a license which is distributed by the
+  LEDA Software GmbH, Postfach 151101, 66041 Saarbruecken, FRG
+  (fax +49 681 31104).
+
+  Copyright (c) 1991-1997  by  Max-Planck-Institut fuer Informatik
+  Im Stadtwald, 66123 Saarbruecken, Germany     
+  All rights reserved.
+ 
*******************************************************************************/
#include <LEDA/d3_hull.h>
#include <LEDA/slist.h>
#include <math.h>

static int vis_count = 0;

bool CHECK_HULL(const GRAPH<POINT,int>& H)
{  
   VECTOR c_vec(3);  // zero vector

   int n = 0;

   node v;
   forall_nodes(v,H) 
     if (H.degree(v) > 0) 
     { c_vec += H[v].to_vector();
       n++;
      }

   c_vec = c_vec/n;

   edge e;
   forall_edges(e,H)
    { edge e1 = H.cyclic_adj_succ(e);
      edge e2 = H.cyclic_adj_pred(e);
      POINT a = H[source(e)];
      POINT b = H[target(e)];
      POINT c = H[target(e1)];
      POINT d = H[target(e2)];

      if (orientation(a,b,c,c_vec) < 0)
      { error_handler(0,"orientation error");
        return false;
       }

      if (orientation(a,b,c,d) < 0)
      { error_handler(0,"local non-convexity");
        return false;
       }

     }

   return true;
}
     


void del_visible_faces(GRAPH<POINT,int>& H, edge& border_e, POINT p)
{ 
  slist<edge> S;
  slist<edge> del_edges;

  vis_count++;

  // initialize S with a visible face cycle

  edge x = border_e;
  do { S.append(x);
       H[x] = vis_count;
       x = H.face_cycle_succ(x);
  } while (x != border_e);


  while (!S.empty())
  { edge e = S.pop();
    edge r = H.reversal(e);
    if (H[r] == vis_count) continue;
    edge r1 = H.face_cycle_succ(r);
    POINT A = H[source(r)];
    POINT B = H[source(r1)];
    POINT C = H[target(r1)];
    if (orientation(A,B,C,p) > 0)
       border_e = e;
    else
      { edge x = r;
        do { H[x] = vis_count;
             edge y = H.reversal(x);
             if (H[y] != vis_count)
                S.append(x);
             else 
              { del_edges.append(x);
                del_edges.append(y);
               }
             x = H.face_cycle_succ(x);
         } while (x != r);
       }
   }

   // remove visible edges

  edge e;
  forall(e,del_edges)
  { node x = source(e);
    node y = target(e);
    H.del_edge(e);
    if (H.degree(x) == 0) H.del_node(x);
    if (H.degree(y) == 0) H.del_node(y);
  }

}


void join_coplanar_faces(GRAPH<POINT,int>& H)
{
  list<edge> L;
  edge_array<bool> considered(H,false);

  edge e;
  forall_edges(e,H)
  { if (considered[e]) continue;
    edge r = H.reversal(e);
    edge e1 = H.face_cycle_succ(e);
    edge r1 = H.face_cycle_succ(r);

    POINT A = H[source(e)];
    POINT B = H[source(e1)];
    POINT C = H[target(e1)];
    POINT D = H[target(r1)];

    if (orientation(A,B,C,D) == 0) 
    { L.append(e);
      L.append(r);
      considered[e] = considered[r] = true;
     }
   }

   forall(e,L) H.del_edge(e);
}

void d2_add_point(GRAPH<POINT,int>& H, edge& border_e, const POINT& D, 
                                                       const POINT& p)
{
  node v = H.new_node(p);

  // construct upper tangent
  edge up = border_e;
  do  up = H.face_cycle_pred(up); 
  while (orientation(p,H[source(up)],H[target(up)],D) >= 0);

  up = H.face_cycle_succ(up);
  edge x = H.new_edge(up,v);
  edge y = H.new_edge(v,source(up));
  H.set_reversal(x,y);

  // construct lower tangent
  edge down = border_e;
  while (orientation(p,H[source(down)],H[target(down)],D) >= 0)
       down = H.face_cycle_succ(down);

  x = H.new_edge(down,v);
  y = H.new_edge(v,source(down));
  H.set_reversal(x,y);

  // remove visible edges
  while (target(up) != v)
  { edge e = H.face_cycle_succ(up);
    node u = source(up);
    H.del_edge(H.reversal(up));
    H.del_edge(up);
    if (H.degree(u) == 0) H.del_node(u);
    up = e;
   }

  border_e = H.last_adj_edge(v);
}




void compute_hull(list<POINT>& L, GRAPH<POINT,int>& H)
{
  POINT A = L.pop();
  node a = H.new_node(A);

  if (L.empty()) return;


  POINT B = L.pop();
  while ( !L.empty() && collinear(A,B,L.head()) ) B = L.pop();
  node b = H.new_node(B);

  if (L.empty()) // all points are collinear
  { edge x = H.new_edge(a,b,0);
    edge y = H.new_edge(b,a,0);
    H.set_reversal(x,y);
    return;
   }


  // construct a triangle

  POINT C = L.pop();
  node c = H.new_node(C);

  POINT D = point_on_positive_side(A,B,C);

  H.new_edge(a,b,0);
  H.new_edge(a,c,0);
  H.new_edge(b,a,0);
  H.new_edge(b,c,0);

  if (orientation(A,B,C,D) > 0)
   { H.new_edge(c,a,0);
     H.new_edge(c,b,0);
    }
  else
   { H.new_edge(c,b,0);
     H.new_edge(c,a,0);
    }

  H.make_map();

  if (L.empty()) return;

  edge border_e = H.last_adj_edge(c);

  int  dim = 2;

  while (!L.empty())
  {
    POINT p = L.pop();

    if (dim == 2) 
    { int orient = orientation(A,B,C,p);
      if (orient == 0)
      { d2_add_point(H,border_e,D,p);
        continue;
       }
      if (orient < 0) border_e = H.reversal(border_e);
      dim = 3;
    }


    // 3-dimensional case
    // remove all faces visible from p


    del_visible_faces(H,border_e,p);

    // and re-triangulate

    node v = H.new_node(p);

    POINT q;

    if (!L.empty()) 
       q = L.head(); 
    else 
       q = p;

    bool vflag = orientation(p,H[source(border_e)],H[target(border_e)],q) < 0;
    edge be = H.face_cycle_succ(border_e);


    edge e = be;
    do { bool construct_edge;
         if (orientation(p,H[source(e)],H[target(e)],q) < 0)
         { construct_edge = !vflag;
           vflag = true;
           border_e = e;
          }
         else
         { construct_edge = true;
           vflag = false;
          }

         if (construct_edge)
         { edge x = H.new_edge(e,v,0,after);
           edge y = H.new_edge(v,source(e),0);
           H.set_reversal(x,y);
          }
         e = H.face_cycle_succ(e);
     } while (source(e) != source(be));

    if (H.outdeg(v) == 0) H.del_node(v);
  }

  //if (dim == 3) join_coplanar_faces(H);
}


static void random_sample(int n, const list<POINT>& L, list<POINT>& L1)
{
  int N        = L.length();
  list_item* A = new list_item[N];
  int i        = 0;

  list_item it;
  forall_items(it,L) A[i++] = it;

  L1.clear();

  while (n--)
    L1.append(L.inf(A[rand_int(0,N-1)]));

  delete[] A;
}
    


int simplify(list<POINT>& L, int n)
{

float t = used_time();

  GRAPH<POINT,int> G0;
  list<POINT> L0;

  random_sample(n,L,L0);
  L0.sort();
  L0.unique();
  compute_hull(L0,G0);

  if (G0.number_of_nodes() >= 3*n/4) return 0;
  if (G0.number_of_nodes() <  6)     return 0;

  int M = G0.number_of_edges()/3;

  double* AX = new double[M];
  double* AY = new double[M];
  double* AZ = new double[M];
  double* AW = new double[M];
  double* NX = new double[M];
  double* NY = new double[M];
  double* NZ = new double[M];

  int m = 0;

  edge_array<bool>  considered(G0,false);
  edge e0;
  forall_edges(e0,G0)
  { if (considered[e0]) continue;

    edge x = e0;
    do { considered[x] = true;
         x = G0.face_cycle_succ(x);
    } while (x != e0);

    edge e1 = G0.face_cycle_succ(e0);
    edge e2 = G0.face_cycle_succ(e1);
    POINT a = G0[source(e0)];
    POINT b = G0[source(e1)];
    POINT c = G0[source(e2)];
    AX[m] = a.XD();
    AY[m] = a.YD();
    AZ[m] = a.ZD();
    AW[m] = a.WD();
    double X1 = b.XD()*a.WD() - a.XD()*b.WD();
    double Y1 = b.YD()*a.WD() - a.YD()*b.WD();
    double Z1 = b.ZD()*a.WD() - a.ZD()*b.WD();
    double X2 = c.XD()*a.WD() - a.XD()*c.WD();
    double Y2 = c.YD()*a.WD() - a.YD()*c.WD();
    double Z2 = c.ZD()*a.WD() - a.ZD()*c.WD();
    NX[m] = Z1*Y2 - Y1*Z2;
    NY[m] = X1*Z2 - Z1*X2;
    NZ[m] = Y1*X2 - X1*Y2;
    m++;
  }

  int count = 0;

  list_item it = L.first();
  while (it)
  { 
    list_item next = L.succ(it);
    POINT p = L[it];
    double px = p.XD();
    double py = p.YD();
    double pz = p.ZD();
    double pw = p.WD();
    int i = 0;
    while (i < m)
    { double aw = AW[i];
      if ( NX[i]*(px*aw-AX[i]*pw) +
           NY[i]*(py*aw-AY[i]*pw) +
           NZ[i]*(pz*aw-AZ[i]*pw) <= 0 ) break;
      i++;
     }
    if (i == m) 
    { count++;
      L.del_item(it); 
     }
    it = next;
   }

  delete[] AX;
  delete[] AY;
  delete[] AZ;
  delete[] AW;
  delete[] NX;
  delete[] NY;
  delete[] NZ;

//cout << string("simplify:  |L| = %d   %5.2f",L.length(),used_time(t));

  return count;
} 



void D3_HULL(const list<POINT>& L0, GRAPH<POINT,int>& H, int n)
{ H.clear();
  if (L0.empty()) return;
  list<POINT> L = L0;

  if (n == -1)
  { double l = L0.length();
    n = int(pow(l,0.6));
   }

  if (n > 6) simplify(L,n);

  L.sort();
  L.unique();
  compute_hull(L,H);
}
