#define CGAL_PMP_REPAIR_POLYGON_SOUP_VERBOSE 1

#include "Scene_polygon_soup_item.h"
#include "Scene_surface_mesh_item.h"

#include <CGAL/Three/Viewer_interface.h>
#include <CGAL/Three/Triangle_container.h>
#include <CGAL/Three/Edge_container.h>
#include <CGAL/Three/Point_container.h>
#include <CGAL/Three/Three.h>

#include "triangulate_primitive.h"

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>

#include <CGAL/array.h>
#include <CGAL/IO/OFF.h>
#include <CGAL/Polygon_mesh_processing/compute_normal.h>
#include <CGAL/Polygon_mesh_processing/polygon_soup_to_polygon_mesh.h>
#include <CGAL/Polygon_mesh_processing/polygon_mesh_to_polygon_soup.h>
#include <CGAL/Polygon_mesh_processing/orient_polygon_soup.h>
#include <CGAL/Polygon_mesh_processing/orientation.h>
#include <CGAL/Polygon_mesh_processing/repair.h>
#include <CGAL/Polygon_mesh_processing/repair_polygon_soup.h>
#include <CGAL/Polygon_mesh_processing/triangulate_faces.h>
#include <CGAL/Polygon_2.h>
#include <CGAL/version.h>

#include <QObject>
#include <QApplication>
#include <QtDebug>

#include <boost/accumulators/accumulators.hpp>
#include <boost/accumulators/statistics/stats.hpp>
#include <boost/accumulators/statistics/mean.hpp>
#include <boost/accumulators/statistics/min.hpp>
#include <boost/accumulators/statistics/max.hpp>
#include <boost/accumulators/statistics/median.hpp>

#include <algorithm>
#include <iostream>
#include <map>
#include <queue>
#include <set>
#include <stack>
#include <streambuf>
#include <vector>

using namespace CGAL::Three;
typedef Viewer_interface Vi;
typedef Triangle_container Tc;
typedef Edge_container Ec;
typedef Point_container Pc;

struct Scene_polygon_soup_item_priv{

  typedef Polygon_soup::Polygons::const_iterator Polygons_iterator;
  typedef EPICK::Point_3 Point_3;

  Scene_polygon_soup_item_priv(Scene_polygon_soup_item* parent)
    : soup(nullptr),
      oriented(false)
  {
    item = parent;
    nb_polys = 0;
    nb_lines = 0;
    nb_nm_edges = 0;
    invalidate_stats();
  }
  ~Scene_polygon_soup_item_priv()
  {
    if(soup)
    {
      delete soup;
      soup = nullptr;
    }
  }
  void compute_normals_and_vertices(void) const;
  void triangulate_polygon(Polygons_iterator, int ) const;
  void invalidate_stats();
  void compute_stats();
  mutable QOpenGLShaderProgram *program;


  enum Face_names {
      Flat_facets=0,
      Smooth_facets,
  };
  enum Edge_names {
    Edges = 0,
    NM_edges
  };

  Polygon_soup* soup;
  bool oriented;
  mutable std::vector<float> positions_poly;
  mutable std::vector<float> positions_lines;
  mutable std::vector<float> f_colors;
  mutable std::vector<float> v_colors;
  mutable std::vector<float> normals;
  mutable std::vector<float> positions_nm_lines;
  mutable std::size_t nb_nm_edges;
  mutable std::size_t nb_polys;
  mutable std::size_t nb_lines;
  bool is_triangle, is_quad, stats_computed;
  double minl, maxl, meanl, midl, mini, maxi, ave;
  std::size_t nb_null_edges, nb_degen_faces;

  Scene_polygon_soup_item* item;

};

typedef Scene_polygon_soup_item_priv Priv;


typedef EPICK Traits;
typedef Polygon_soup::Polygon_3 Facet;

void
Scene_polygon_soup_item_priv::triangulate_polygon(Polygons_iterator pit, int polygon_id) const
{
  const CGAL::qglviewer::Vec off = static_cast<CGAL::Three::Viewer_interface*>(CGAL::QGLViewer::QGLViewerPool().first())->offset();
  EPICK::Vector_3 offset(off.x,off.y,off.z);

    //Computes the normal of the facet
    Traits::Vector_3 normal = CGAL::NULL_VECTOR;

    //Newell's method
    for (std::size_t i = 0; i < pit->size() ; ++ i){
      const Point_3& pa = soup->points[pit->at(i)];
      const Point_3& pb = soup->points[pit->at((i+1)%pit->size())];
      double x = normal.x() + (pa.y()-pb.y())*(pa.z()+pb.z());
      double y = normal.y() + (pa.z()-pb.z())*(pa.x()+pb.x());
      double z = normal.z() + (pa.x()-pb.x())*(pa.y()+pb.y());
      normal = Traits::Vector_3(x,y,z);
    }
    if (normal == CGAL::NULL_VECTOR) // No normal could be computed, return
      return;

    typedef std::pair<EPICK::Point_3, std::size_t> PointAndId;

    std::size_t it = 0;
    std::size_t it_end =pit->size();
    std::vector<PointAndId> pointIds;
    do {
      PointAndId pointId;

      pointId.first = soup->points[pit->at(it)]+offset;
      pointId.second = pit->at(it);
      pointIds.push_back(pointId);
    } while( ++it != it_end );

    //iterates on the internal faces to add the vertices to the positions
    //and the normals to the appropriate vectors
    auto f = [&](auto& ffit, auto& v2v) {
      if (ffit.info().is_external)
        return;

      positions_poly.push_back(ffit.vertex(0)->point().x());
      positions_poly.push_back(ffit.vertex(0)->point().y());
      positions_poly.push_back(ffit.vertex(0)->point().z());

      positions_poly.push_back(ffit.vertex(1)->point().x());
      positions_poly.push_back(ffit.vertex(1)->point().y());
      positions_poly.push_back(ffit.vertex(1)->point().z());

      positions_poly.push_back(ffit.vertex(2)->point().x());
      positions_poly.push_back(ffit.vertex(2)->point().y());
      positions_poly.push_back(ffit.vertex(2)->point().z());

      CGAL::IO::Color color;
      if (!soup->fcolors.empty())
        color = soup->fcolors[polygon_id];
      for (int i = 0; i < 3; i++)
      {
        normals.push_back(normal.x());
        normals.push_back(normal.y());
        normals.push_back(normal.z());
        if (!soup->fcolors.empty())
        {
          f_colors.push_back(static_cast<float>(color.red()) / 255);
          f_colors.push_back(static_cast<float>(color.green()) / 255);
          f_colors.push_back(static_cast<float>(color.blue()) / 255);
        }
        if (!soup->vcolors.empty())
        {
          CGAL::IO::Color vcolor = soup->vcolors[v2v[ffit.vertex(i)]];
          v_colors.push_back(static_cast<float>(vcolor.red()) / 255);
          v_colors.push_back(static_cast<float>(vcolor.green()) / 255);
          v_colors.push_back(static_cast<float>(vcolor.blue()) / 255);
        }
      }
      };

    try {
      FacetTriangulator<SMesh, EPICK, std::size_t> triangulation(pointIds, normal);
      triangulation.per_face(f);
    }
    catch (...) {
      FacetTriangulator<SMesh, EPICK, std::size_t, CGAL::Exact_intersections_tag> triangulation(pointIds, normal);
      triangulation.per_face(f);
    }
}

void
Scene_polygon_soup_item_priv::compute_normals_and_vertices() const{

    //get the vertices and normals
    const CGAL::qglviewer::Vec offset = static_cast<CGAL::Three::Viewer_interface*>(CGAL::QGLViewer::QGLViewerPool().first())->offset();

    typedef Polygon_soup::Polygons::size_type size_type;
    positions_poly.resize(0);
    positions_lines.resize(0);
    normals.resize(0);
    positions_nm_lines.resize(0);
    f_colors.resize(0);
    v_colors.resize(0);
    soup->fill_edges();
    int nb = 0;
    for(Polygons_iterator it = soup->polygons.begin();
        it != soup->polygons.end(); ++it)
    {
        if(it->size()!=3)
        {
            triangulate_polygon(it, nb);
        }
        else{
            const Point_3& pa = soup->points[it->at(0)];
            const Point_3& pb = soup->points[it->at(1)];
            const Point_3& pc = soup->points[it->at(2)];

            EPICK::Vector_3 n = CGAL::cross_product(pb-pa, pc -pa);
            n = n / std::sqrt(n * n);

            normals.push_back(n.x());
            normals.push_back(n.y());
            normals.push_back(n.z());

            normals.push_back(n.x());
            normals.push_back(n.y());
            normals.push_back(n.z());


            normals.push_back(n.x());
            normals.push_back(n.y());
            normals.push_back(n.z());


            for(size_type i = 0; i < it->size(); ++i)
            {
                const Point_3& p = soup->points[it->at(i)];
                positions_poly.push_back(p.x()+offset.x);
                positions_poly.push_back(p.y()+offset.y);
                positions_poly.push_back(p.z()+offset.z);
                if(!soup->fcolors.empty())
                {
                  const CGAL::IO::Color color = soup->fcolors[nb];
                    f_colors.push_back(static_cast<float>(color.red())/255);
                    f_colors.push_back(static_cast<float>(color.green())/255);
                    f_colors.push_back(static_cast<float>(color.blue())/255);
                }

                if(!soup->vcolors.empty())
                {
                  const CGAL::IO::Color color = soup->vcolors[it->at(i)];
                  v_colors.push_back(static_cast<float>(color.red())/255);
                  v_colors.push_back(static_cast<float>(color.green())/255);
                  v_colors.push_back(static_cast<float>(color.blue())/255);
                }
            }
        }
        nb++;

        //Lines
        for(size_type i = 0; i < it->size(); ++i)
        {

            const Point_3& pa = soup->points[it->at(i)];
            const Point_3& pb = soup->points[it->at((i+1)%it->size())];
            positions_lines.push_back(pa.x()+offset.x);
            positions_lines.push_back(pa.y()+offset.y);
            positions_lines.push_back(pa.z()+offset.z);

            positions_lines.push_back(pb.x()+offset.x);
            positions_lines.push_back(pb.y()+offset.y);
            positions_lines.push_back(pb.z()+offset.z);
        }
    }

    //Non manifold edges
    for(const Polygon_soup::Edge& edge :
                    soup->non_manifold_edges)
    {

        const Point_3& a = soup->points[edge[0]];
        const Point_3& b = soup->points[edge[1]];
        positions_nm_lines.push_back(a.x()+offset.x);
        positions_nm_lines.push_back(a.y()+offset.y);
        positions_nm_lines.push_back(a.z()+offset.z);

        positions_nm_lines.push_back(b.x()+offset.x);
        positions_nm_lines.push_back(b.y()+offset.y);
        positions_nm_lines.push_back(b.z()+offset.z);
    }

}


Scene_polygon_soup_item::Scene_polygon_soup_item()
{
  d = new Scene_polygon_soup_item_priv(this);
  for(int i = 1; i>=0; --i)
  {
    setTriangleContainer(i,
                         new Tc(Vi::PROGRAM_WITH_LIGHT, false));
    setEdgeContainer(i,
                     new Ec(Vi::PROGRAM_NO_SELECTION, false));
  }
  setPointContainer(0,
                    new Pc(Vi::PROGRAM_NO_SELECTION, false));
}

Scene_polygon_soup_item::~Scene_polygon_soup_item()
{
  delete d;
}

Scene_polygon_soup_item*
Scene_polygon_soup_item::clone() const {
  Scene_polygon_soup_item* new_soup = new Scene_polygon_soup_item();
  new_soup->d->soup = d->soup->clone();
  new_soup->d->oriented = d->oriented;
  return new_soup;
}


bool
Scene_polygon_soup_item::load(std::istream& in)
{
  if (!d->soup)
    d->soup = new Polygon_soup();
  else
    d->soup->clear();

  bool result = CGAL::IO::read_OFF(in, d->soup->points, d->soup->polygons,
                                      CGAL::parameters::vertex_color_output_iterator(std::back_inserter(d->soup->vcolors))
                                                        .face_color_output_iterator(std::back_inserter(d->soup->fcolors)));

  invalidateOpenGLBuffers();
  return result;
}

void Scene_polygon_soup_item::init_polygon_soup(std::size_t nb_pts, std::size_t nb_polygons){
  if(!d->soup)
    d->soup = new Polygon_soup;
  d->soup->clear();
  d->soup->points.reserve(nb_pts);
  d->soup->polygons.reserve(nb_polygons);
  d->soup->vcolors.resize(0);
  d->soup->fcolors.resize(0);
  d->oriented = false;
}

template<class PolygonMesh>
void polygon_mesh_to_soup(PolygonMesh& mesh, Polygon_soup& soup)
{
  soup.clear();
  CGAL::Polygon_mesh_processing::polygon_mesh_to_polygon_soup(mesh, soup.points, soup.polygons);
  soup.fill_edges();
}

void Scene_polygon_soup_item::load(Scene_surface_mesh_item* sm_item) {
  if(!sm_item) return;
  if(!sm_item->face_graph()) return;

  if(!d->soup)
    d->soup = new Polygon_soup;
  polygon_mesh_to_soup(*sm_item->face_graph(), *d->soup);
  invalidateOpenGLBuffers();
}
void
Scene_polygon_soup_item::setDisplayNonManifoldEdges(const bool b)
{

  d->soup->display_non_manifold_edges = b;
}

bool
Scene_polygon_soup_item::displayNonManifoldEdges() const {

  return d->soup->display_non_manifold_edges;
}

void Scene_polygon_soup_item::shuffle_orientations()
{
  for(Polygon_soup::size_type i = 0, end = d->soup->polygons.size();
      i < end; ++i)
  {
    if(std::rand() % 2 == 0) d->soup->inverse_orientation(i);
  }
  invalidateOpenGLBuffers();
}

void Scene_polygon_soup_item::inside_out()
{
  for(Polygon_soup::size_type i = 0, end = d->soup->polygons.size();
      i < end; ++i)
  {
    d->soup->inverse_orientation(i);
  }
  invalidateOpenGLBuffers();
}

void Scene_polygon_soup_item::repair(bool erase_dup, bool req_same_orientation)
{
  QApplication::setOverrideCursor(Qt::BusyCursor);
  CGAL::Polygon_mesh_processing::repair_polygon_soup(
        d->soup->points,
        d->soup->polygons,
        CGAL::parameters::erase_all_duplicates(erase_dup)
                         .require_same_orientation(req_same_orientation));
  QApplication::restoreOverrideCursor();
  invalidateOpenGLBuffers();
}

bool Scene_polygon_soup_item::triangulate()
{
  QApplication::setOverrideCursor(Qt::BusyCursor);

  bool success = true;

  CGAL::Polygon_mesh_processing::triangulate_polygons(d->soup->points, d->soup->polygons);

  QApplication::restoreOverrideCursor();
  invalidateOpenGLBuffers();

  return success;
}

bool
Scene_polygon_soup_item::orient(std::vector<std::size_t>& non_manifold_vertices)
{
  struct Visitor : public CGAL::Polygon_mesh_processing::Default_orientation_visitor
  {
    const Polygon_soup::Polygons& polygons;
    std::set<std::size_t>& nm_vertices;
    std::set< std::pair<std::size_t, std::size_t> > nm_edges;
    std::vector<CGAL::IO::Color>& vcolors;

    Visitor(const Polygon_soup::Polygons& polygons,
            std::set<std::size_t>& nm_vertices,
            std::vector<CGAL::IO::Color>& vcolors)
      : polygons(polygons)
      , nm_vertices(nm_vertices)
      , vcolors(vcolors)
    {}

    void non_manifold_edge(std::size_t v1, std::size_t v2, std::size_t)
    {
      nm_edges.insert(CGAL::make_sorted_pair(v1, v2));
    }

    void link_connected_polygons(std::size_t v, const std::vector<std::size_t>& polygons_in_current_cycle)
    {
      // check if the current components of polygon incident to the link of the vertex contains
      // a non-manifold edge. If no, then it is a "pure" non-manifold vertex in this component
      if (!nm_edges.empty())
        for(std::size_t pid : polygons_in_current_cycle)
        {
          for (std::size_t i=0; i<polygons[pid].size(); ++i)
          {
            if (polygons[pid][i]==v)
            {
              std::size_t vn1 = i==0?(polygons[pid].size()-1):(i-1);
              std::size_t vn2 = i+1;
              if (vn2==polygons[pid].size()) vn2=0;

              if (nm_edges.count(CGAL::make_sorted_pair(polygons[pid][i], polygons[pid][vn1]))==1) return;
              if (nm_edges.count(CGAL::make_sorted_pair(polygons[pid][i], polygons[pid][vn2]))==1) return;
              break;
            }
            else
              CGAL_assertion(i!=polygons[pid].size()-1); // the vertex is incident to the polygon
          }
        }
      nm_vertices.insert(v);
    }

    void non_manifold_vertex(std::size_t, std::size_t nb_link_ccs)
    {
      if (!vcolors.empty())
        vcolors.resize(vcolors.size()+nb_link_ccs-1);
    }

    void vertex_id_in_polygon_replaced(std::size_t, std::size_t input_id, std::size_t new_id)
    {
      if (!vcolors.empty())
        vcolors[new_id]=vcolors[input_id];
    }
  };
  if(isEmpty() || d->oriented)
    return true; // nothing to do
  QApplication::setOverrideCursor(Qt::WaitCursor);
  d->oriented=true;

  //first skip degenerate polygons
  Polygon_soup::Polygons valid_polygons;
  valid_polygons.reserve(d->soup->polygons.size());
  for(Polygon_soup::Polygon_3& polygon : d->soup->polygons)
  {
    std::set<std::size_t> vids;
    bool to_remove=false;
    for(std::size_t id : polygon)
    {
      if (!vids.insert(id).second){
        to_remove=true;
        break;
      }
    }
    if (!to_remove) valid_polygons.push_back(polygon);
  }
  QApplication::restoreOverrideCursor();
  if (valid_polygons.size()!=d->soup->polygons.size())
    d->soup->polygons.swap(valid_polygons);

  bool res;
  std::set<std::size_t> nm_v_set;
  Visitor visitor(valid_polygons, nm_v_set, d->soup->vcolors);
  QApplication::setOverrideCursor(Qt::WaitCursor);
  res =  CGAL::Polygon_mesh_processing::
    orient_polygon_soup(d->soup->points, d->soup->polygons, CGAL::parameters::visitor(visitor));
  non_manifold_vertices.assign(nm_v_set.begin(), nm_v_set.end());
  QApplication::restoreOverrideCursor();
  return res;
}


bool
Scene_polygon_soup_item::save(std::ostream& out) const
{

  typedef Polygon_soup::size_type size_type;
  CGAL::File_writer_OFF writer;
  writer.write_header(out,
                      d->soup->points.size(),
                      0,
                      d->soup->polygons.size());
  for(size_type i = 0, end = d->soup->points.size();
      i < end; ++i)
  {
    const Point_3& p = d->soup->points[i];
    writer.write_vertex( p.x(), p.y(), p.z() );
  }
  writer.write_facet_header();
  for(size_type i = 0, end = d->soup->polygons.size();
      i < end; ++i)
  {
    const Polygon_soup::Polygon_3& polygon = d->soup->polygons[i];
    const size_type size = polygon.size();
    writer.write_facet_begin(size);
    for(size_type j = 0; j < size; ++j) {
      writer.write_facet_vertex_index(polygon[j]);
    }
    writer.write_facet_end();
  }
  writer.write_footer();

  return (bool) out;
}

bool
Scene_polygon_soup_item::exportAsSurfaceMesh(SMesh *out_surface_mesh)
{
  std::vector<std::size_t> dum;
  orient(dum);
  CGAL::Polygon_mesh_processing::polygon_soup_to_polygon_mesh< CGAL::Surface_mesh<Point_3> >(
    d->soup->points, d->soup->polygons, *out_surface_mesh);
  std::size_t rv = CGAL::Polygon_mesh_processing::remove_isolated_vertices(*out_surface_mesh);
  if(rv > 0){
    std::cerr << "Ignore isolated vertices: " << rv << std::endl;
    out_surface_mesh->collect_garbage();
  }
  if(out_surface_mesh->vertices().size() > 0) {
    return true;
  }
  return false;
}
QString
Scene_polygon_soup_item::toolTip() const
{

  if(!d->soup)
    return QString();

  QString str = QObject::tr("<p><b>%1</b> (mode: %5, color: %6)<br />"
                     "<i>Polygon soup</i></p>"
                     "<p>Number of vertices: %2<br />"
                     "Number of polygons: %3</p>")
    .arg(this->name())
    .arg(d->soup->points.size())
    .arg(d->soup->polygons.size())
    .arg(this->renderingModeName())
    .arg(this->color().name());
    str += QString("<br />Number of isolated vertices: %1<br />").arg(getNbIsolatedvertices());
    return str;
}

void
Scene_polygon_soup_item::draw(CGAL::Three::Viewer_interface* viewer) const {
    if(d->soup == nullptr) return;
    if(!isInit(viewer))
      initGL(viewer);
    if ( getBuffersFilled() &&
         ! getBuffersInit(viewer))
    {
      initializeBuffers(viewer);
      setBuffersInit(viewer, true);
    }
    if(!getBuffersFilled())
    {
      computeElements();
      initializeBuffers(viewer);
    }

    if(renderingMode() == Flat || renderingMode() == FlatPlusEdges)
    {

      if(d->soup->fcolors.empty())
        getTriangleContainer(Priv::Flat_facets)->setColor(this->color());
      getTriangleContainer(Priv::Flat_facets)->setAlpha(alpha());
      getTriangleContainer(Priv::Flat_facets)->draw(viewer, d->soup->fcolors.empty());
    }
    else if(renderingMode() == Gouraud)
    {
      if(d->soup->vcolors.empty())
        getTriangleContainer(Priv::Smooth_facets)->setColor(this->color());
      getTriangleContainer(Priv::Smooth_facets)->setAlpha(alpha());
      getTriangleContainer(Priv::Smooth_facets)->draw(viewer, d->soup->vcolors.empty());
    }
  }

void
Scene_polygon_soup_item::drawPoints(CGAL::Three::Viewer_interface* viewer) const {

    if(d->soup == nullptr) return;
    if(!isInit(viewer))
      initGL(viewer);
    if ( getBuffersFilled() &&
         ! getBuffersInit(viewer))
    {
      initializeBuffers(viewer);
      setBuffersInit(viewer, true);
    }
    if(!getBuffersFilled())
    {
      computeElements();
      initializeBuffers(viewer);
    }
    getPointContainer(0)->setColor(this->color());
    getPointContainer(0)->draw(viewer, true);
}

void
Scene_polygon_soup_item::drawEdges(CGAL::Three::Viewer_interface* viewer) const {
    if(d->soup == nullptr) return;
    if(!isInit(viewer))
      initGL(viewer);
    if ( getBuffersFilled() &&
         ! getBuffersInit(viewer))
    {
      initializeBuffers(viewer);
      setBuffersInit(viewer, true);
    }
    if(!getBuffersFilled())
    {
      computeElements();
      initializeBuffers(viewer);
    }
    getEdgeContainer(Priv::Edges)->setColor(QColor(Qt::black));
    getEdgeContainer(Priv::Edges)->draw(viewer, true);
    if(displayNonManifoldEdges())
    {
      getEdgeContainer(Priv::NM_edges)->setColor(QColor(Qt::red));
      //draw the edges
      getEdgeContainer(Priv::NM_edges)->draw(viewer, true);
    }
}

bool
Scene_polygon_soup_item::isEmpty() const {

  return (d->soup == nullptr || d->soup->points.empty());
}
void
Scene_polygon_soup_item::invalidateOpenGLBuffers()
{
    compute_bbox();
    for(int i=0; i<2; ++i)
    {
      getTriangleContainer(i)->reset_vbos(ALL);
      getEdgeContainer(i)->reset_vbos(ALL);
    }
    getPointContainer(0)->reset_vbos(ALL);
    setBuffersFilled(false);
    d->invalidate_stats();
}

void Scene_polygon_soup_item::compute_bbox() const {

  if (isEmpty())
    return;
  const Point_3& p = *(d->soup->points.begin());
  CGAL::Bbox_3 bbox(p.x(), p.y(), p.z(), p.x(), p.y(), p.z());
  for(Polygon_soup::Points::const_iterator it = d->soup->points.begin();
      it != d->soup->points.end();
      ++it) {
    bbox = bbox + it->bbox();
  }
  setBbox(Bbox(bbox.xmin(),bbox.ymin(),bbox.zmin(),
               bbox.xmax(),bbox.ymax(),bbox.zmax()));
}

void
Scene_polygon_soup_item::new_vertex(const double& x,
                                    const double& y,
                                    const double& z)
{

    d->soup->points.push_back(Point_3(x, y, z));
}

void
Scene_polygon_soup_item::new_triangle(const std::size_t i,
                                      const std::size_t j,
                                      const std::size_t k)
{

  Polygon_soup::Polygon_3 new_polygon(3);
  new_polygon[0] = i;
  new_polygon[1] = j;
  new_polygon[2] = k;
  d->soup->polygons.push_back(new_polygon);
}

template <class Point, typename Polygon>
void Scene_polygon_soup_item::load(const std::vector<Point>& points, const std::vector<Polygon>& polygons)
{
    if(!d->soup)
        d->soup = new Polygon_soup;
    d->soup->clear();

    /// add points
    d->soup->points.reserve(points.size());
    for(const Point& p : points)
            d->soup->points.push_back( Point_3(p[0], p[1], p[2]) );

    /// add polygons
    std::size_t nb_polygons=polygons.size();
    d->soup->polygons.resize(nb_polygons);
    for(std::size_t i=0; i<nb_polygons; ++i)
        d->soup->polygons[i].assign(polygons[i].begin(), polygons[i].end());

    /// fill non-manifold edges container
    //soup->fill_edges();
    d->oriented = false;
    invalidateOpenGLBuffers();
}

template <class Point, class Polygon>
void Scene_polygon_soup_item::load(const std::vector<Point>& points, const std::vector<Polygon>& polygons,
                                   const std::vector<CGAL::IO::Color>& fcolors,
                                   const std::vector<CGAL::IO::Color>& vcolors)
{
    load (points, polygons);

    d->soup->fcolors.reserve (fcolors.size());
    std::copy (fcolors.begin(), fcolors.end(), std::back_inserter (d->soup->fcolors));

    d->soup->vcolors.reserve (vcolors.size());
    std::copy (vcolors.begin(), vcolors.end(), std::back_inserter (d->soup->vcolors));
}
// Force the instantiation of the template function for the types used in the STL_io_plugin. This is needed
// because the d-pointer forbid the definition in the .h for this function.
template SCENE_POLYGON_SOUP_ITEM_EXPORT void Scene_polygon_soup_item::load<EPICK::Point_3, std::vector<int> >
(const std::vector<EPICK::Point_3>& points, const std::vector<std::vector<int> >& polygons);
template SCENE_POLYGON_SOUP_ITEM_EXPORT void Scene_polygon_soup_item::load<CGAL::Epick::Point_3, std::vector<std::size_t> >
(const std::vector<CGAL::Epick::Point_3>& points, const std::vector<std::vector<std::size_t> >& polygons);
template SCENE_POLYGON_SOUP_ITEM_EXPORT void Scene_polygon_soup_item::load<CGAL::Epick::Point_3, std::vector<std::size_t> >
(const std::vector<CGAL::Epick::Point_3>& points, const std::vector<std::vector<std::size_t> >& polygons,
 const std::vector<CGAL::IO::Color>& fcolors,
 const std::vector<CGAL::IO::Color>& vcolors);

// Local Variables:
// c-basic-offset: 4
// End:

const Scene_polygon_soup_item::Points& Scene_polygon_soup_item::points() const { return d->soup->points; }
const Scene_polygon_soup_item::Polygons& Scene_polygon_soup_item::polygons() const { return d->soup->polygons; }
bool Scene_polygon_soup_item::isDataColored() { return d->soup->fcolors.size()>0 || d->soup->vcolors.size()>0;}
std::vector<CGAL::IO::Color> Scene_polygon_soup_item::getVColors() const{return d->soup->vcolors;}
std::vector<CGAL::IO::Color> Scene_polygon_soup_item::getFColors() const{return d->soup->fcolors;}

void Scene_polygon_soup_item::itemAboutToBeDestroyed(Scene_item *item)
{
  Scene_item::itemAboutToBeDestroyed(item);
  if(d && item == this)
  {
    if(d->soup)
    {
      delete d->soup;
      d->soup=nullptr;
    }
  }
}

const Polygon_soup::Edges&
Scene_polygon_soup_item::non_manifold_edges() const
{
  return d->soup->non_manifold_edges;
}

void Scene_polygon_soup_item::initializeBuffers(Viewer_interface *v) const
{
  getTriangleContainer(Priv::Flat_facets)->initializeBuffers(v);
  getTriangleContainer(Priv::Flat_facets)->setFlatDataSize(d->nb_polys);

  getTriangleContainer(Priv::Smooth_facets)->initializeBuffers(v);
  getTriangleContainer(Priv::Smooth_facets)->setFlatDataSize(d->nb_polys);

  getEdgeContainer(Priv::Edges)->initializeBuffers(v);
  getEdgeContainer(Priv::Edges)->setFlatDataSize(d->nb_lines);

  getEdgeContainer(Priv::NM_edges)->initializeBuffers(v);
  getEdgeContainer(Priv::NM_edges)->setFlatDataSize(d->nb_nm_edges);

  getPointContainer(0)->initializeBuffers(v);
  getPointContainer(0)->setFlatDataSize(d->nb_lines);

  d->normals.resize(0);
  d->positions_poly.resize(0);
  d->normals.shrink_to_fit();
  d->positions_poly.shrink_to_fit();
  d->v_colors.resize(0);
  d->v_colors.shrink_to_fit();
  d->positions_lines.resize(0);
  d->positions_lines.shrink_to_fit();
  d->positions_nm_lines.resize(0);
  d->positions_nm_lines.shrink_to_fit();
}

void Scene_polygon_soup_item::computeElements() const
{
  QApplication::setOverrideCursor(Qt::WaitCursor);
  d->compute_normals_and_vertices();

  getTriangleContainer(Priv::Flat_facets)->allocate(
        Tc::Flat_vertices,
        d->positions_poly.data(),
        static_cast<int>(d->positions_poly.size()*sizeof(float)));

  getTriangleContainer(Priv::Flat_facets)->allocate(
        Tc::Flat_normals,
        d->normals.data(),
        static_cast<int>(d->normals.size()*sizeof(float)));
  if(!d->f_colors.empty())
  {
    getTriangleContainer(Priv::Flat_facets)->allocate(
          Tc::FColors,
          d->f_colors.data(),
          static_cast<int>(d->f_colors.size()*sizeof(float)));
  }
  getTriangleContainer(Priv::Smooth_facets)->allocate(
        Tc::Flat_vertices,
        d->positions_poly.data(),
        static_cast<int>(d->positions_poly.size()*sizeof(float)));

  getTriangleContainer(Priv::Smooth_facets)->allocate(
        Tc::Flat_normals,
        d->normals.data(),
        static_cast<int>(d->normals.size()*sizeof(float)));

  if(!d->v_colors.empty())
  {
    getTriangleContainer(Priv::Smooth_facets)->allocate(
          Tc::VColors,
          d->v_colors.data(),
          static_cast<int>(d->v_colors.size()*sizeof(float)));
  }

  d->nb_polys = d->positions_poly.size();

  getEdgeContainer(Priv::Edges)->allocate(
        Ec::Vertices,
        d->positions_lines.data(),
        static_cast<int>(d->positions_lines.size()*sizeof(float)));

  getPointContainer(0)->allocate(
        Pc::Vertices,
        d->positions_lines.data(),
        static_cast<int>(d->positions_lines.size()*sizeof(float)));


  getEdgeContainer(Priv::NM_edges)->allocate(
        Ec::Vertices,
        d->positions_nm_lines.data(),
        static_cast<int>(d->positions_nm_lines.size()*sizeof(float)));

  d->nb_nm_edges = d->positions_nm_lines.size();
  d->nb_lines = d->positions_lines.size();

  setBuffersFilled(true);
  QApplication::restoreOverrideCursor();
}

CGAL::Three::Scene_item::Header_data Scene_polygon_soup_item::header() const
{
  CGAL::Three::Scene_item::Header_data data;

  //categories
  data.categories.append(std::pair<QString,int>(QString("Properties"),2));
  data.categories.append(std::pair<QString,int>(QString("Vertices"),1));
  data.categories.append(std::pair<QString,int>(QString("Polygons"),2));
  data.categories.append(std::pair<QString,int>(QString("Edges"),6));
  data.categories.append(std::pair<QString,int>(QString("Angles"),3));

  //titles
  data.titles.append(QString("Pure Triangle"));
  data.titles.append(QString("Pure Quad"));

  data.titles.append(QString("#Points"));

  data.titles.append(QString("#Polygons"));
  data.titles.append(QString("#Degenerate Polygons"));

  data.titles.append(QString("#Edges"));
  data.titles.append(QString("#Degenerate Edges"));
  data.titles.append(QString("Minimum Length"));
  data.titles.append(QString("Maximum Length"));
  data.titles.append(QString("Median Length"));
  data.titles.append(QString("Mean Length"));

  data.titles.append(QString("Minimum"));
  data.titles.append(QString("Maximum"));
  data.titles.append(QString("Average"));

  return data;
}

QString Scene_polygon_soup_item::computeStats(int type)
{
  if(!d->stats_computed)
    d->compute_stats();

  switch(type)
  {
  case NB_VERTICES:
    return QString::number(d->soup->points.size());
  case NB_FACETS:
    return QString::number(d->soup->polygons.size());
  case NB_EDGES:
    return QString::number(d->nb_lines/6);

  case NB_DEGENERATE_FACES:
  {
    if(d->is_triangle)
    {
      return QString::number(d->nb_degen_faces);
    }
    else
      return QString("n/a");
  }

  case MIN_LENGTH:
    return QString::number(d->minl);
  case MAX_LENGTH:
    return QString::number(d->maxl);
  case MED_LENGTH:
    return QString::number(d->midl);
  case MEAN_LENGTH:
    return QString::number(d->meanl);
  case NB_DEGENERATE_EDGES:
    return QString::number(d->nb_null_edges);

  case MIN_ANGLE:
    return QString::number(d->mini);
  case MAX_ANGLE:
    return QString::number(d->maxi);
  case MEAN_ANGLE:
    return QString::number(d->ave);

  case IS_PURE_TRIANGLE:
    if(d->is_triangle)
      return QString("yes");
    else
      return QString("no");
  case IS_PURE_QUAD:
    if (d->is_quad)
      return QString("yes");
    else
      return QString("no");
  }
  return QString();
}

void
Scene_polygon_soup_item_priv::
invalidate_stats()
{
  is_triangle = true;
  is_quad = true;
  minl=0;
  maxl=0;
  meanl=0;
  midl=0;
  mini=0;
  maxi=0;
  ave=0;
  nb_null_edges=0;
  nb_degen_faces=0;
  stats_computed = false;
}


void
Scene_polygon_soup_item_priv::compute_stats()
{
  using namespace boost::accumulators;
  accumulator_set< double,
    features< tag::min, tag::max, tag::mean , tag::median> > edges_acc;
  accumulator_set< double,
    features< tag::min, tag::max, tag::mean > > angles_acc;
  double rad_to_deg = 180. / CGAL_PI;


  for(auto poly : soup->polygons)
  {
    if(poly.size() != 3)
      is_triangle = false;
    if(poly.size() != 4)
      is_quad = false;
    for(std::size_t i = 0; i< poly.size(); ++i)
    {
      Polygon_soup::Point_3 a(soup->points[poly[i]]),
          b(soup->points[poly[(i+1)%poly.size()]]),
          c(soup->points[poly[(i+2)%poly.size()]]);
      if (a == b)
        ++nb_null_edges;
      edges_acc(CGAL::sqrt(CGAL::squared_distance(a, b)));
      typename Traits::Vector_3 ba(b, a);
      typename Traits::Vector_3 bc(b, c);
      double cos_angle = (ba * bc)
        / std::sqrt(ba.squared_length() * bc.squared_length());
      if(cos_angle == CGAL_PI || cos_angle == 0)
        ++nb_degen_faces;
      angles_acc(std::acos(cos_angle) * rad_to_deg);
    }
  }

  minl = extract_result< tag::min >(edges_acc);
  maxl = extract_result< tag::max >(edges_acc);
  meanl = extract_result< tag::mean >(edges_acc);
  midl =  extract_result< tag::median >(edges_acc);
  mini = extract_result< tag::min >(angles_acc);
  maxi = extract_result< tag::max >(angles_acc);
  ave = extract_result< tag::mean >(angles_acc);

  stats_computed = true;
}
