#include <utilities/log.h>
#include <boost/foreach.hpp>

AlphaSimplex2D::
AlphaSimplex2D(const Delaunay2D::Vertex& v): alpha_(0), attached_(false)
{
    for (int i = 0; i < 3; ++i)
        if (v.face()->vertex(i) != Vertex_handle() && v.face()->vertex(i)->point() == v.point())
            Parent::add(v.face()->vertex(i));
}

AlphaSimplex2D::
AlphaSimplex2D(const Delaunay2D::Edge& e): attached_(false)
{
    Face_handle f = e.first;
    for (int i = 0; i < 3; ++i)
        if (i != e.second)
            Parent::add(f->vertex(i));
}

// AlphaSimplex2D::
// AlphaSimplex2D(const Delaunay2D::Edge& e, const SimplexSet& simplices, const Delaunay2D& Dt): attached_(false)
// {
//     Face_handle f = e.first;
//     for (int i = 0; i < 3; ++i)
//         if (i != e.second)
//             Parent::add(f->vertex(i));
// 
//     VertexSet::const_iterator v = static_cast<const Parent*>(this)->vertices().begin();
//     const Point& p1 = (*v++)->point();
//     const Point& p2 = (*v)->point();
// 
//     Face_handle o = f->neighbor(e.second);
//     if (o == Face_handle())
//     {
//         alpha_ = CGAL::squared_radius(p1, p2);
//         return;
//     }
//     int oi = o->index(f);
// 
//     attached_ = false;
//     if (!Dt.is_infinite(f->vertex(e.second)) &&
//         CGAL::side_of_bounded_circle(p1, p2,
//                                      f->vertex(e.second)->point()) == CGAL::ON_BOUNDED_SIDE)
//         attached_ = true;
//     else if (!Dt.is_infinite(o->vertex(oi)) &&
//              CGAL::side_of_bounded_circle(p1, p2,
//                                           o->vertex(oi)->point()) == CGAL::ON_BOUNDED_SIDE)
//         attached_ = true;
//     else
//         alpha_ = CGAL::squared_radius(p1, p2);
// 
//     if (attached_)
//     {
//         if (Dt.is_infinite(f))
//             alpha_ = simplices.find(AlphaSimplex2D(*o))->alpha();
//         else if (Dt.is_infinite(o))
//             alpha_ = simplices.find(AlphaSimplex2D(*f))->alpha();
//         else
//             alpha_ = std::min(simplices.find(AlphaSimplex2D(*f))->alpha(),
//                               simplices.find(AlphaSimplex2D(*o))->alpha());
//     }
// }

// VR modified for periodic DT 
AlphaSimplex2D::
AlphaSimplex2D(const Delaunay2D::Edge& e, const SimplexSet& simplices, const Delaunay2D& Dt, const RealValue r): attached_(false)
{
    Face_handle f = e.first;
    for (int i = 0; i < 3; ++i)
        if (i != e.second)
            Parent::add(f->vertex(i));

    Face_handle o = f->neighbor(e.second);
    int oi = o->index(f);
    alpha_ = r; 

    attached_ = false;
// First check if the circle centred on the edge contains its f-opposite vertex
	Point fpt = Dt.point(Dt.periodic_point(e.first,e.second));
	Delaunay2D::Segment fseg = Dt.segment(e.first,e.second);
	if (CGAL::side_of_bounded_circle(fseg[0],fseg[1],fpt) == CGAL::ON_BOUNDED_SIDE)
		attached_ = true;

// Then check if the circle centred on the edge contains its o-opposite vertex. 	
	Point opt = Dt.point(Dt.periodic_point(o,oi)); 
	Delaunay2D::Segment oseg = Dt.segment(o,oi);
	if (CGAL::side_of_bounded_circle(oseg[0],oseg[1],opt) == CGAL::ON_BOUNDED_SIDE)
		attached_ = true;


    if (attached_)
    {
    	//std::cout << "found an attached edge" << std::endl;
    	alpha_ = std::min(simplices.find(AlphaSimplex2D(*f))->alpha(),
                              simplices.find(AlphaSimplex2D(*o))->alpha());
    }
}

// AlphaSimplex2D::
// AlphaSimplex2D(const Delaunay2D::Face& f): attached_(false)
// {
//     for (int i = 0; i < 3; ++i)
//         Parent::add(f.vertex(i));
//     VertexSet::const_iterator v = static_cast<const Parent*>(this)->vertices().begin();
//     Point p1 = (*v++)->point();
//     Point p2 = (*v++)->point();
//     Point p3 = (*v)->point();
//     alpha_ = CGAL::squared_radius(p1, p2, p3);
// }

// VR: new initialisation for periodic DT faces. 

AlphaSimplex2D::
AlphaSimplex2D(const Delaunay2D::Face& f)
{
    for (int i = 0; i < 3; ++i)
        Parent::add(f.vertex(i));
}

AlphaSimplex2D::
AlphaSimplex2D(const Delaunay2D::Face& f, const RealValue r): attached_(false)
{
    for (int i = 0; i < 3; ++i)
        Parent::add(f.vertex(i));

    alpha_ = r;
}


bool
AlphaSimplex2D::AlphaOrder::
operator()(const AlphaSimplex2D& first, const AlphaSimplex2D& second) const
{
    if (first.alpha() == second.alpha())
        return (first.dimension() < second.dimension());
    else
        return (first.alpha() < second.alpha());
}

std::ostream&
AlphaSimplex2D::
operator<<(std::ostream& out) const
{
    for (VertexSet::const_iterator cur = Parent::vertices().begin();
                                   cur != Parent::vertices().end(); ++cur)
        out << **cur << ", ";
    out << "value = " << value();

    return out;
}

// void fill_simplex_set(const Delaunay2D& Dt, AlphaSimplex2D::SimplexSet& simplices)
// {
//     for(Face_iterator cur = Dt.finite_faces_begin(); cur != Dt.finite_faces_end(); ++cur)
//         simplices.insert(AlphaSimplex2D(*cur));
//     rInfo("Faces inserted");
//     for(Edge_iterator cur = Dt.finite_edges_begin(); cur != Dt.finite_edges_end(); ++cur)
//         simplices.insert(AlphaSimplex2D(*cur, simplices, Dt));
//     rInfo("Edges inserted");
//     for(Vertex_iterator cur = Dt.finite_vertices_begin(); cur != Dt.finite_vertices_end(); ++cur)
//         simplices.insert(AlphaSimplex2D(*cur));
//     rInfo("Vertices inserted");
// }

//VR: modified function for periodic DT

void fill_simplex_set(const Delaunay2D& Dt, AlphaSimplex2D::SimplexSet& simplices)
{
    // Compute all simplices with their alpha values and attachment information
	int counter = 0;
    for(Face_iterator cur = Dt.faces_begin(); cur != Dt.faces_end(); ++cur){
        //Delaunay2D::Periodic_triangle ptri = Dt.periodic_triangle(cur);
        //Delaunay2D::Triangle tri = Dt.triangle(ptri);
        Delaunay2D::Triangle tri = Dt.triangle(cur);
        RealValue sqrad = CGAL::squared_radius(tri[0], tri[1], tri[2]);
        simplices.insert(AlphaSimplex2D(*cur, sqrad));
        counter++;
    }
    //std::cout << "Number of triangles inserted in filtration " << counter << std::endl;  
    rInfo("Facets inserted");
    counter = 0;
    for(Edge_iterator cur = Dt.edges_begin(); cur != Dt.edges_end(); ++cur){
        //Delaunay2D::Periodic_segment pseg = Dt.periodic_segment(*cur);
        //Delaunay2D::Segment seg = Dt.segment(pseg);
        Delaunay2D::Segment seg = Dt.segment(cur);
        RealValue sqrad = CGAL::squared_radius(seg[0],seg[1]);
        simplices.insert(AlphaSimplex2D(*cur, simplices, Dt, sqrad));
        counter++;
    }
    //std::cout << "Number of edges inserted in filtration " << counter << std::endl;  
    rInfo("Edges inserted");      
    for(Vertex_iterator cur = Dt.vertices_begin(); cur != Dt.vertices_end(); ++cur)
        simplices.insert(AlphaSimplex2D(*cur));
    rInfo("Vertices inserted");
    
}



template<class Filtration>
void fill_complex(const Delaunay2D& Dt, Filtration& filtration)
{
    // Compute all simplices with their alpha values and attachment information
    // TODO: this can be optimized; the new Filtration can act as a SimplexSet
    AlphaSimplex2D::SimplexSet simplices;
    fill_simplex_set(Dt, simplices);
    BOOST_FOREACH(const AlphaSimplex2D& s, simplices)
        filtration.push_back(s);
}