/***************************************************************************
 *   Copyright (C) 2010 by Martin Weusten                                  *
 *   martin dot weusten (ät" rwth minus aachen.de                          *
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 *   This program is distributed in the hope that it will be useful,       *
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
 *   GNU General Public License for more details.                          *
 *                                                                         *
 *   You should have received a copy of the GNU General Public License     *
 *   along with this program; if not, write to the                         *
 *   Free Software Foundation, Inc.,                                       *
 *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
 ***************************************************************************/

void BuildingPolyhedron::triangulate(CGAL::Polyhedron polyhedron,
                                     std::vector<CML::Vector3f>& triangleList)
{
   CGAL::Polyhedron::Facet_iterator itrFace = polyhedron.facets_begin();
   CGAL::Polyhedron::Facet_iterator endFace = polyhedron.facets_end();
   for (; itrFace != endFace; ++itrFace) {
      CGAL::Polyhedron::Facet face = (*itrFace);

      // extract polygon vertices
      Ring<CML::Vector3f> polygon(face.facet_degree());
      CGAL::Polyhedron::Facet::Halfedge_around_facet_const_circulator circHalfEdge, endHalfEdge;
      circHalfEdge = endHalfEdge = face.facet_begin();
      for (; circHalfEdge != endHalfEdge; ++circHalfEdge) {
         polygon.add(toVector3f(circHalfEdge->vertex()->point()));
      }
   
      // triangulate polygon
      triangulate(polygon, triangleList);
   }
}

void BuildingPolyhedron::triangulate(CDL::Ring<CML::Vector3f> polygon,
                                     std::vector<CML::Vector3f>& triangleList)
{
   int size = polygon.size();
   Q_ASSERT_X(size >= 3, "BuildingPolyhedron", "A polygon with less than 3 vertices?");

   // compute face normal and projection
   CML::Vector3f normal = CML::faceNormal(polygon[0], polygon[1], polygon[2]);
   CML::Vector3f axis0(1,0,0);
   if (abs(axis0 * normal) > 0.8) { axis0 = CML::Vector3f(0,1,0); }
   CML::Vector3f axis1 = axis0 % normal;
   axis0 = axis1 % normal;
   CML::Matrix33f projection(axis0, axis1, normal);

   // project vertices into 2D space and find minimum x value
   float minValue = CML::FLOAT_MAX;
   int minIndex = -1;
   Ring<CML::Vector2f> polygon2D(size);
   for (int i=0; i<size; ++i) {
      CML::Vector3f vertex = polygon[i] * projection;
      polygon2D.add(CML::Vector2f(vertex.x, vertex.y));
      if (vertex.x < minValue) {
         minValue = vertex.x;
         minIndex = i;
      }
   }

   // get sign of convex vertices (minimum vertex must be convex)
   int convexSign = CML::sgn(adjacentEdgeCross(polygon2D, minIndex));

   // generate triangles using curner cutting algorithm
   for (int i=0; size>3; ++i) {
      // skip concave corners
      if (adjacentEdgeCross(polygon2D, i) * convexSign < 0) { continue; }

      // extract the triangle
      const CML::Vector2f& A = polygon2D[i-1];
      const CML::Vector2f& B = polygon2D[i  ];
      const CML::Vector2f& C = polygon2D[i+1];

      // check if other vertex is inside the triangle
      bool vertexInside = false;
      for (int j=i+2; j<i+size-1 && !vertexInside; ++j) {
         // skip convex corners
         if (adjacentEdgeCross(polygon2D, j) * convexSign > 0) { continue; }
         // check if vertex j is inside triangle ABC
         vertexInside = CML::pointInTriangle2D(A,B,C,polygon2D[j]);
      }
      if (vertexInside) { continue; }

      // ok, let's cut the corner
      triangleList.push_back(polygon[i-1]);
      triangleList.push_back(polygon[i  ]);
      triangleList.push_back(polygon[i+1]);
      polygon.remove(i);
      polygon2D.remove(i);
      --size;
      --i;
   }
   // add remaining triangle
   triangleList.push_back(polygon[0]);
   triangleList.push_back(polygon[1]);
   triangleList.push_back(polygon[2]);
}

inline float BuildingPolyhedron::adjacentEdgeCross(const Ring<CML::Vector2f>& vertices2D, const int corner) {
   CML::Vector2f edge1 = vertices2D[corner-1] - vertices2D[corner];
   CML::Vector2f edge2 = vertices2D[corner+1] - vertices2D[corner];
   return edge1 % edge2;
}


namespace CDL {

   /** A dynamic array with circular access methods.
    * Circular means that index access is always modulo the current size.
    * E.g index \c -1 does map to the last element.
    * @author Martin Weusten */
   template<typename T>
   class Ring {
   public:
      /** Create new Ring with given initial capacity. */
      Ring(const int capacity = 5);
      Ring(const Ring<T>& r);
      ~Ring();

      Ring<T>& operator=(const Ring<T>& r);

      /** Add element at the end. Operation takes constant runtime if capacity is sufficient */
      void add(const T& x);

      /** Remove a specific element.
        * Note that this method works circular. Operation takes linear runtime. */
      void remove(int i);

      /** Set size to zero. */
      inline void clear() { m_size = 0; }

      inline bool empty() const { return m_size == 0; }
      inline int size() const { return m_size; }
      inline int capacity() const { return m_capacity; }
      
      /** Resize the internal array. Elements at the end are dropped if the
       * capacity is smaller than the current size. Operation takes linear runtime. */
      void capacity(const int n);

      /** Access element at index. Note that this operator works circular. */
      inline T& operator[](const int i) {
         return m_data[mod(i, m_size)];
      }

      /** Access element at index. Note that this operator works circular. */
      inline const T& operator[](const int i) const {
         return m_data[mod(i, m_size)];
      }
   private:
      T* m_data;
      int m_size;
      int m_capacity;
      float m_growFactor;

      /** Simple modulo which does also work for negative a. */
      static inline int mod(const int a, const int n) {
         int r = a - n * (a/n);
         return r < 0 ? r+n : r;
      }
   };

   // ...

}

//********************************************************************
Function TTesselator.SplitPolygon(APolygon: TIndexList): TBool32;
//----------------------------------------------------
Procedure CreateContour(AInternPoly: TIndexList;
                     ANewCount,AOldIndex: TTessIndex);
Var NewIndex,OldIndex: TTessIndex;
Begin
   OldIndex:= AOldIndex;
   For NewIndex:= 0 To ANewCount-1 Do Begin
      AInternPoly[NewIndex]:= APolygon[OldIndex];
      OldIndex:= Successor(OldIndex,Length(APolygon));
   End;
End;
//----------------------------------------------------
Var
   Polygon0,Polygon1: TIndexList;
   Count0,Count1,Index0,Index1,Actual,I: TTessIndex;
   Found: TBool32;
Begin
   Result:= FALSE;
    Found:= GetDiagonal(APolygon,Index0,Index1);

   If Found
   Then Begin

      // Create first polygon ----------------------------------
      If Index1 > Index0
      Then Count0:= Index1-Index0+1
      Else Count0:= Length(APolygon)-(Index0-Index1)+1;
      SetLength(Polygon0,Count0);
      CreateContour(Polygon0,Count0,Index0);

      If Length(Polygon0) > 3
      Then SplitPolygon(Polygon0) // RECURSE with first polygon
      Else Begin  // End: write triangle to output
         Move(Polygon0[0],fTriangles[fActSlot],3*SizeOf(TTessIndex));
         Inc(fActSlot,3);
      End;
      Finalize(Polygon0);

      // Create second polygon ---------------------------------
      Count1:= Length(APolygon)-Count0+2;
      SetLength(Polygon1,Count1);
      CreateContour(Polygon1,Count1,Index1);

      If Length(Polygon1) > 3
      Then SplitPolygon(Polygon1) // RECURSE with second polygon
      Else Begin  // End: write triangle to output
         Move(Polygon1[0],fTriangles[fActSlot],3*SizeOf(TTessIndex));
         Inc(fActSlot,3);
      End;
      Finalize(Polygon1);
     
      Result:= TRUE;
   End
   Else Raise ETesselError.Create('Could not split polygon');
End;
//********************************************************************
Function TTesselator.GetDiagonal(APolygon: TIndexList;
                   Var AIndex0,AIndex1: TTessIndex): TBool32;
Var
   RotationType: TRotationType;
   Prior,Actual,Next,Count,Trouble,I: TTessIndex;
   FoundTrouble: TBool32;
   ActP,PriorP,NextP,SomeP: TTessVertex2D;
Begin
   // Find first ear
   Count:= Length(APolygon);
   RotationType:= rotNone;

   I:= 0;
   While RotationType <> rotCCW Do
   Begin
      Actual:= I;
      Prior:= Predecessor(Actual,Count);
      Next:= Successor(Actual,Count);

      PriorP:= fContour[APolygon[Prior]];
      ActP  := fContour[APolygon[Actual]];
      NextP := fContour[APolygon[Next]];

      RotationType:= GetRotationType
         (fContour[APolygon[Prior]],
          fContour[APolygon[Actual]],
          fContour[APolygon[Next]]);

      Inc(I); //============>
   End;

   // If this is an ear, then we have to look for trouble :)
   If RotationType = rotCCW
   Then Begin
      I:= Successor(Next,Count);
      FoundTrouble:= FALSE;

      While (I <> Prior) And (Not FoundTrouble) Do
      Begin
         SomeP:= fContour[APolygon[I]];
         If (I <> Actual) And (I <> Prior) And (I <> Next) And
            (PointInTriangle(fContour[APolygon[I]],
                fContour[APolygon[Prior]],
                fContour[APolygon[Actual]],
                fContour[APolygon[Next]]))
            Then Begin
               Trouble:= I;
               FoundTrouble:= TRUE;
            End
            Else I:= Successor(I,Count); //=======>
      End;
   End;

   If RotationType = rotCCW  // This is an ear!
   Then Begin
      If Not FoundTrouble
         Then Begin AIndex0:= Prior; AIndex1:= Next; End
         Else Begin AIndex0:= Actual; AIndex1:= Trouble; End;
      Result:= TRUE;
   End
   Else Result:= FALSE;
End;
//********************************************************************