/*******************************************************************************
+
+  LEDA 3.5
+
+  polygon.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/plane_alg.h>
#include <LEDA/graph.h>
#include <math.h>
#include <assert.h>

//------------------------------------------------------------------------------
// POLYGON
//
// by S. Naeher (1995,1996)
//------------------------------------------------------------------------------


//extern void SWEEP_SEGMENTS(const list<SEGMENT>&,GRAPH<POINT,SEGMENT>&);
//extern void MULMULEY_SEGMENTS(const list<SEGMENT>&,GRAPH<POINT,SEGMENT>&);



POLYGON_REP::POLYGON_REP(const list<SEGMENT>& sl) : seg_list(sl)
{ SEGMENT s;
  forall(s,sl) pt_list.append(s.start());
 }


polygon POLYGON::to_polygon() const
{ list<point> L;
  SEGMENT s;
  forall(s,ptr()->seg_list) L.append(s.source().to_point());
  return polygon(L);
}


ostream& operator<<(ostream& out, const POLYGON& p) 
{ p.vertices().print(out);
  return out << endl;
 } 

istream& operator>>(istream& in,  POLYGON& p) 
{ 
#if defined(POLYGON)
  list<POINT> L; 
  L.read(in,'\n'); 
  p = POLYGON(L); 
#endif
  return in;
}


NUMBER POLYGON::compute_area(const list<SEGMENT>& seg_list) const
{

  if (seg_list.length() < 3) return 0;

  list_item it = seg_list.get_item(1);
  POINT     p  = seg_list[it].source();

  it = seg_list.succ(it);

  NUMBER A  = 0;

  while (it)
  { SEGMENT s = seg_list[it];
    A += ::area(p,s.source(),s.target());
    it = seg_list.succ(it);
   }

  return A;
}


static void check_simplicity(const list<SEGMENT>& seg_list)
{ GRAPH<POINT,SEGMENT> G;
  SWEEP_SEGMENTS(seg_list,G);
  //MULMULEY_SEGMENTS(seg_list,G);
  node v;
  forall_nodes(v,G)
    if (G.degree(v) != 2)
      error_handler(1,"POLYGON: POLYGON must be simple");
}


POLYGON::POLYGON(const list<SEGMENT>& sl)
{ PTR = new POLYGON_REP(sl); }
 

POLYGON::POLYGON(const list<POINT>& pl)
{ 
  list<SEGMENT> seglist;

  for(list_item it = pl.first(); it; it = pl.succ(it))
     seglist.append(SEGMENT(pl[it],pl[pl.cyclic_succ(it)]));

  if (compute_area(seglist) < 0)
  { // reverse edge list
    seglist.clear();
    for(list_item it = pl.last(); it; it = pl.pred(it))
      seglist.append(SEGMENT(pl[it],pl[pl.cyclic_pred(it)]));
   }
  check_simplicity(seglist);

  PTR = new POLYGON_REP(seglist);
}


POLYGON POLYGON::translate(NUMBER dx, NUMBER dy) const
{ list<SEGMENT> L;
  SEGMENT s;
  forall(s,ptr()->seg_list) L.append(s.translate(dx,dy));
  return POLYGON(L);
}


POLYGON POLYGON::translate(const VECTOR& v) const
{ list<SEGMENT> L;
  SEGMENT s;
  forall(s,ptr()->seg_list) L.append(s.translate(v));
  return POLYGON(L);
}


POLYGON POLYGON::rotate90(const POINT& p) const
{ list<SEGMENT> L;
  SEGMENT s;
  forall(s,ptr()->seg_list) L.append(s.rotate90(p));
  return POLYGON(L);
}


POLYGON POLYGON::reflect(const POINT& p, const POINT& q) const
{ list<SEGMENT> L;
  SEGMENT s;
  forall(s,ptr()->seg_list) L.append(s.reflect(p,q));
  return POLYGON(L);
}


POLYGON POLYGON::reflect(const POINT& p) const
{ list<SEGMENT> L;
  SEGMENT s;
  forall(s,ptr()->seg_list) L.append(s.reflect(p));
  return POLYGON(L);
}



list<POINT> POLYGON::intersection(const SEGMENT& s) const
{ list<POINT> result;
  SEGMENT t;
  POINT p;

  forall(t,ptr()->seg_list) 
    if (s.intersection(t,p))
     if (result.empty()) result.append(p);
     else if (p != result.tail() ) result.append(p);

  return result;
}


list<POINT> POLYGON::intersection(const LINE& l) const
{ list<POINT> result;
  SEGMENT t;
  POINT p;

  forall(t,ptr()->seg_list) 
    if (l.intersection(t,p))
     if (result.empty()) result.append(p);
     else if (p != result.tail() ) result.append(p);

  return result;
}


// intersection or union with POLYGON

static bool right_oriented(SEGMENT s, POINT p)
{
  int orient = orientation(s,p);
  if (orient < 0) return true;
  if (orient > 0) return false;

  VECTOR v1 = s.target() - s.source();
  VECTOR v2 = p - s.source();

  return (v1*v2 > 0);
}


static bool left_oriented(SEGMENT s, POINT p)
{
  int orient = orientation(s,p);
  if (orient > 0) return true;
  if (orient < 0) return false;

  VECTOR v1 = s.target() - s.source();
  VECTOR v2 = p - s.source();

  return (v1*v2 > 0);
}



static bool test_edge(GRAPH<POINT,SEGMENT>& G, edge i1, int mode)
{ 
  node v  = G.target(i1);
  edge i2 = G.cyclic_in_succ(i1);

  node u1 = G.source(i1);
  node u2 = G.source(i2);

  if (i1 == i2) return false;

  // parallel edges
  if (u1 == u2) return i1 == G.first_in_edge(v);

  edge o1 = G.first_adj_edge(v);
  edge o2 = G.last_adj_edge(v);

  // anti-parallel edges
  if (target(o1) == u1 || target(o2) == u1) return false;


  POINT p1 = G[o1].target();
  POINT p2 = G[o2].target();

  SEGMENT si1 = G[i1];
  SEGMENT si2 = G[i2];

  if (mode == 0) // intersection
  { if (orientation(si1,si2.source()) >= 0)
      return orientation(si1,p1) > 0 && orientation(si2,p1) < 0 && 
             orientation(si1,p2) > 0 && orientation(si2,p2) < 0;
   else
      return (orientation(si1,p1) > 0 || orientation(si2,p1) < 0) && 
             (orientation(si1,p2) > 0 || orientation(si2,p2) < 0);
  }
  else // union
  { if (orientation(si1,si2.source()) <= 0)
        return right_oriented(si1,p1) && left_oriented(si2,p1) ||
               right_oriented(si1,p2) && left_oriented(si2,p2);
   else
        return right_oriented(si1,p1) || left_oriented(si2,p1) ||
               right_oriented(si1,p2) || left_oriented(si2,p2);
   }
}



static edge next_edge(GRAPH<POINT,SEGMENT>& G, edge i1, int dir)
{ 
  // if dir = +1 turn left
  // if dir = -1 turn right

  node v = G.target(i1);

  edge o1 = G.first_adj_edge(v);
  edge o2 = G.last_adj_edge(v);

  if (o1 == o2) return o1;

  node w1 = G.target(o1);
  node w2 = G.target(o2);

  if (w1 == w2) return (o1 == G.first_in_edge(w1)) ? o1 : o2;

  SEGMENT si1 = G[i1];
  SEGMENT so1 = G[o1];
  SEGMENT so2 = G[o2];

  int orient1 = orientation(si1,so1.target());
  int orient2 = orientation(si1,so2.target());

  if (orient1 == orient2)
     return (orientation(so1,so2.target()) == dir) ? o2 : o1;
  else
     if (orient1 < orient2)
        return (dir == -1) ? o1 : o2;
     else
        return (dir == -1) ? o2 : o1;
}



LIST_POLYGON_ POLYGON::intersection(const POLYGON& P) const
{
  list<POLYGON> result;
  list<SEGMENT> seg_list;
  GRAPH<POINT,SEGMENT> G;

  SEGMENT s;
  forall(s,ptr()->seg_list) seg_list.append(s);
  forall(s,P.ptr()->seg_list) seg_list.append(s);

  SWEEP_SEGMENTS(seg_list,G);
  //MULMULEY_SEGMENTS(seg_list,G);

  bool borders_intersect = false;

  node v;
  forall_nodes(v,G) 
  { int deg = G.degree(v);
    assert(deg == 2 || deg == 4);
    if (deg == 4) borders_intersect = true;
   }
    

  if ( ! borders_intersect )
  { // no intersections between edges of (*this) and P
    // check for inclusion

    SEGMENT s1 = ptr()->seg_list.head();
    SEGMENT s2 = P.ptr()->seg_list.head();

    if ( P.inside(s1.start()) ) result.append(*this);
    if ( inside(s2.start())   ) result.append(P);                   
   }
  else
  { edge_array<bool> marked(G,false);
    edge e;
    forall_edges(e,G) 
    { node start = source(e);
      if ( !marked[e] && test_edge(G,e,0) )
      { list<SEGMENT> pol;
        do { node v = source(e);
             node w = target(e);
             marked[e] = true;
             pol.append(SEGMENT(G[v],G[w]));
             e = next_edge(G,e,+1);
            } while (source(e) != start);
        result.append(POLYGON(pol));
       }
     }
   }

 return result;

}


LIST_POLYGON_ POLYGON::unite(const POLYGON& P) const
{
  list<POLYGON> result;
  list<SEGMENT> seg_list;
  GRAPH<POINT,SEGMENT> G;

  SEGMENT s;
  forall(s,ptr()->seg_list) seg_list.append(s);
  forall(s,P.ptr()->seg_list) seg_list.append(s);

  SWEEP_SEGMENTS(seg_list,G);
  //MULMULEY_SEGMENTS(seg_list,G);

  bool borders_intersect = false;

  node v;
  forall_nodes(v,G) 
  { int deg = G.degree(v);
    assert(deg == 2 || deg == 4);
    if (deg == 4) borders_intersect = true;
   }

  if (!borders_intersect)
  { 
    // check for inclusion

    SEGMENT s1 = ptr()->seg_list.head();
    SEGMENT s2 = P.ptr()->seg_list.head();

    if ( ! P.inside(s1.start())) result.append(*this);
    if ( ! inside(s2.start())) result.append(P);                   

    return result;
   }

  edge_array<bool> marked(G,false);

  edge e;
  forall_edges(e,G) 
  { node start = source(e);
    if ( !marked[e] && test_edge(G,e,1) )
    { list<SEGMENT> pol;
      do { node v = source(e);
           node w = target(e);
           marked[e] = true;
           pol.append(SEGMENT(G[v],G[w]));
           e = next_edge(G,e,-1);
          } while (source(e) != start);
      result.append(POLYGON(pol));
     }
   }

  if (result.length() > 1)
  { // find outer POLYGON and move it to front of result


    list_item min_it = result.first();
    SEGMENT        s = (result[min_it].ptr()->seg_list).head();
    NUMBER     min_x = s.xcoord1();

    list_item it;
    forall_items(it,result)
      forall(s,result[it].ptr()->seg_list) 
        if (s.xcoord1() < min_x)
        { min_x = s.xcoord1(); 
          min_it = it;
         }
  
    result.move_to_front(min_it);
  }

 return result;
}



bool POLYGON::inside(const POINT& p) const
{
  list<SEGMENT>& seglist = ptr()->seg_list;

  int count = 0;

  NUMBER px = p.xcoord();

  list_item it0 = seglist.first();
  list_item it1 = seglist.first();

  NUMBER x0 = seglist[it0].xcoord2();
  NUMBER x1 = seglist[it1].xcoord2();

  list_item it;

  forall_items(it,seglist)
  { SEGMENT s = seglist[it];
    if (s.xcoord2() < x0) 
    { it0 = it;
      x0 = s.xcoord2();
     }
    if (s.xcoord2() > x1) 
    { it1 = it1;
      x1 = s.xcoord2();
     }
   }

  if (px <= x0 || x1 <= px) return false;

  while (seglist[it0].xcoord2() <= px) it0 = seglist.cyclic_succ(it0);

  it = it0;
  do
  { while (seglist[it].xcoord2() >= px) it = seglist.cyclic_succ(it);
    if (orientation(seglist[it],p) > 0) count++;
    while (seglist[it].xcoord2() <= px) it = seglist.cyclic_succ(it);
    if (orientation(seglist[it],p) < 0) count++;
  } while (it != it0);

  return count % 2;

}


#if defined(USE_FLOAT_KERNEL)

polygon polygon::translate_by_angle(double alpha, double d) const
{ list<segment> L;
  segment s;
  forall(s,ptr()->seg_list) L.append(s.translate_by_angle(alpha,d));
  return polygon(L);
}

polygon polygon::rotate(const point& p, double alpha) const
{ list<segment> L;
  segment s;
  forall(s,ptr()->seg_list) L.append(s.rotate(p,alpha));
  return polygon(L);
}

polygon polygon::rotate(double alpha) const
{ return rotate(point(0,0),alpha); }

#endif


#if defined(USE_RATIONAL_KERNEL)

rat_polygon rat_polygon::translate(integer dx, integer dy, integer dw) const
{ list<rat_segment> L;
  rat_segment s;
  forall(s,ptr()->seg_list) L.append(s.translate(dx,dy,dw));
  return rat_polygon(L);
}

#endif




