//
//  mesh_functions.cpp
//  test1
//
//  Created by Bac Alexandra on 17/02/2022.
//

#include <algorithm>
#include <random>
#include "mesh_functions.hpp"

vector<Mesh::Vertex_index> first_vneighbors(Mesh::Vertex_index v, Mesh &m)
{
	vector<Mesh::Vertex_index> res ;
	CGAL::Vertex_around_target_circulator<Mesh> vbegin(m.halfedge(v),m), vdone(vbegin) ;
	do {
		res.push_back(*vbegin) ;
		++vbegin ;
	} while (vbegin != vdone) ;
	return res ;
}

vector<Mesh::Face_index> first_fneighbors(Mesh::Vertex_index v, Mesh &m)
{
	vector<Mesh::Face_index> res ;
	CGAL::Face_around_target_circulator<Mesh> fbegin(m.halfedge(v),m), fdone(fbegin) ;
	do {
		res.push_back(*fbegin) ;
		++fbegin ;
	} while (fbegin != fdone) ;
	return res ;
}

vector<Mesh::Halfedge_index> first_heneighbors(Mesh::Vertex_index v, Mesh &m)
{
	vector<Mesh::Halfedge_index> res ;
	CGAL::Halfedge_around_target_circulator<Mesh> hebegin(m.halfedge(v),m), hedone(hebegin) ;
	do {
		res.push_back(*hebegin) ;
		++hebegin ;
	} while (hebegin != hedone) ;
	return res ;
}

int compute_degrees(Mesh &m, Mesh::Property_map<vertex_descriptor,int> &deg)
{
    int deg_max = -1 ;
    Mesh::Vertex_range r = m.vertices() ;
    for(Mesh::Vertex_range::iterator it = r.begin() ;
        it != r.end() ; ++it)
    {
        vector<Mesh::Vertex_index> res = first_vneighbors(*it, m) ;
        deg[*it] = res.size() ;
        if (deg_max < deg[*it])
        {
            deg_max = res.size() ;
        }
    }
    return deg_max ;
}

void gaussian_curv(Mesh &m, Mesh::Property_map<vertex_descriptor,double> &Kcourb)
{
	Mesh::Vertex_range r = m.vertices() ;
	for(Mesh::Vertex_range::iterator it = r.begin() ;
		it != r.end() ; ++it)
	{
		vector<Mesh::Halfedge_index> neigh = first_heneighbors(*it, m) ;
		K::Point_3 p1 = m.point(*it) ;
		double aire = 0, angle = 0 ;
		for (int i=0; i<neigh.size(); ++i)
		{
			K::Point_3 p2, p3 ;
			Mesh::Halfedge_index he = neigh.at(i) ;
			p2 = m.point(m.source(he)) ;
			Mesh::Halfedge_index he2 = neigh.at((i+1)%neigh.size()) ;
			p3 = m.point(m.source(he2)) ;
			//cout << "*it / id1he / id2he / id1he2 / id2he2 " << *it << m.source(he) << m.target(he) << m.source(he2) << m.target(he2) << endl ;
			K::Vector_3 n(CGAL::cross_product(p2-p1, p3-p1)) ;
			aire += (sqrt(n.squared_length()))/2 ;
			angle += CGAL::approximate_angle(p2, p1, p3) ;
		}
		//		cout << "angular defect :" << 2*M_PI-(angle/180*M_PI) << endl ;
		//		cout << "aire : " << aire << endl ;
		Kcourb[*it] = (2*M_PI-(angle/180*M_PI))/aire ;
		//		cout << Kcourb[*it] << endl ;
	}
}

CGAL::Color rand_color(std::uniform_int_distribution<> &distr, std::mt19937 &gen)
{
    CGAL::Color c ;
    c.red() = distr(gen) ;
    c.green() = distr(gen) ;
    c.blue() = distr(gen) ;
    return c ;
}

CGAL::Color color_ramp(double vmin, double vmax, double mean, double stdev, double v, CGAL::Color col1, CGAL::Color col2)
{
	stdev = stdev/2 ;
	if (vmin < mean-stdev)
		vmin = mean-stdev ;
	if (vmax > mean+stdev)
		vmax = mean+stdev ;
	
	if (v < vmin)
		v = vmin ;
	else if (v > vmax)
		v = vmax ;
	double per ;
	CGAL::Color c ;
	CGAL::Color col_first, col_second ;
	if (v <= mean)
	{
		per = (v-vmin)/(mean-vmin) ;
		col_first = col1 ;
		col_second = CGAL::white() ;
	}
	else
	{
		per = (v-mean)/(vmax-mean) ;
		col_first = CGAL::white() ;
		col_second = col2 ;
	}
	
	c.red() = floor(col_first.red() + per*(col_second.red()-col_first.red())) ;
	c.green() = floor(col_first.green() + per*(col_second.green()-col_first.green())) ;
	c.blue() = floor(col_first.blue() + per*(col_second.blue()-col_first.blue())) ;
	return c ;
}



K::Point_3 barycenter (const vector<Mesh::Vertex_index> &verts, Mesh &m)
{
	// Calcul du barycentre
	K::Vector_3 bv(0,0,0) ;
	K::Point_3 Zero(0,0,0), p ;
	for (int i=0; i<verts.size(); ++i)
	{
		p = m.point((verts.at(i))) ;
		bv += (p-K::Point_3(0,0,0)) ;
	}
	bv /= verts.size() ;
	return Zero+bv ;
}

vector<face_descriptor> bary_sub(Mesh::Face_index f, Mesh &m)
{
	// Récupération des sommets de la face
	vector<vertex_descriptor> verts ;
	vector<face_descriptor> res ;
	CGAL::Vertex_around_face_circulator<Mesh> vit(m.halfedge(f), m), vend(vit) ;
	do
	{
		verts.push_back(*vit) ;
		++vit ;
	} while (vit != vend) ;
	cout << "nombre de sommets sur la face : " << verts.size() << endl ;
	K::Point_3 bary(barycenter(verts, m)) ;
	cout << "Barycentre : " << bary  << endl ;
	
	// Suppression de la face
	CGAL::Euler::remove_face(m.halfedge(f), m) ;
	
	// Ajout du barycentre
	Mesh::Vertex_index bi = m.add_vertex(bary) ;
	// Ajout des faces
	for (int i=0; i<verts.size(); ++i)
	{
		cout << bi << " " << verts.at(i) << " " << verts.at((i+1)%verts.size()) << endl ;
		Mesh::Face_index ftmp = m.add_face(verts.at(i), verts.at((i+1)%verts.size()), bi) ;
		assert(ftmp != Mesh::null_face()) ;
		res.push_back(ftmp) ;
		cout << ftmp << " " ;
	}
	cout << endl ;
	cout << "nombre de faces ajoutées : " << res.size() << endl ;
	
	return res ;
}

vector<Mesh::Halfedge_index> subdiv_curved(Mesh &m, const Mesh::Property_map<vertex_descriptor,double> &Kcourb, double kmean, double kstdev)
{
	vector<Mesh::Vertex_index> sorted_vert ;
	Mesh::Vertex_range r = m.vertices() ;
	
	for (Mesh::Vertex_range::iterator it = r.begin(); it != r.end(); ++it)
	{
		sorted_vert.push_back(*it) ;
	}
	// Tri des sommets par ordre décroissant de |Kcourb|
	std::sort(sorted_vert.begin(), sorted_vert.end(), // Lambda expression begins
			  [&Kcourb](vertex_descriptor &v1, vertex_descriptor &v2) {
		 return abs(Kcourb[v1]) > abs(Kcourb[v2]) ;
	 } // end of lambda expression) ;
			  );
	
	// Ajout d'un booléen sur les faces pour marquage de subdivision
	Mesh::Property_map<Mesh::halfedge_index,bool> fsub;
	bool created_fsub;
	boost::tie(fsub, created_fsub) = m.add_property_map<Mesh::halfedge_index,bool>("f:flag",false);
	assert(created_fsub);
	
	int i =0 ;
	double thresh = kmean+kstdev/2 ;
	while (abs(Kcourb[sorted_vert.at(i)]) > thresh) {
		// Sommet à subdiviser -> marquage des face voisines
		cout << "sommet : " << sorted_vert.at(i) << " / " << abs(Kcourb[sorted_vert.at(i)]) << endl ;
		CGAL::Face_around_target_circulator<Mesh> fit(m.halfedge(sorted_vert.at(i)), m), fend(fit) ;
		do {
			fsub[m.halfedge(*fit)] = true ;
			++fit ;
		} while (fit != fend) ;
		++i ;
	}
	
	vector<Mesh::Halfedge_index> res ;
//	// Subdivision de toutes les faces marquées
	Mesh::Halfedge_range fr = m.halfedges() ;
	for(Mesh::Halfedge_range::iterator it = fr.begin(); it != fr.end(); ++it)
	{
		if(fsub[*it])
		{
			vector<face_descriptor> tmp = bary_sub(m.face(*it), m) ;
			for (int j=0; j<tmp.size(); ++j)
				res.push_back(m.halfedge(tmp.at(j))) ;
		}
	}
	
	// Ajout de toutes les faces marquées au vecteur de sortie
//	Mesh::Halfedge_range fr = m.halfedges() ;
//	for(Mesh::Halfedge_range::iterator it = fr.begin(); it != fr.end(); ++it)
//	{
//		if(fsub[*it])
//			res.push_back(*it) ;
//	}
	m.remove_property_map(fsub) ;
	return res ;
}

vector<Mesh::Halfedge_index> edges_from_faces (vector<Mesh::Halfedge_index> &hef, Mesh &m)
{
	vector<Mesh::Halfedge_index> res ;
	// Ajout d'un booléen sur les he pour marquage de flip
	Mesh::Property_map<Mesh::halfedge_index,bool> flip;
	bool created_flip;
	boost::tie(flip, created_flip) = m.add_property_map<Mesh::halfedge_index,bool>("he:flag",false);
	assert(created_flip);
	
	// Marquage des 1/2 aretes à flipper
	for (int i=0; i<hef.size(); ++i)
	{
		// pour chaque 1/2 arete de face, marquage des 1/2 arêtes de son bord
		CGAL::Halfedge_around_target_circulator<Mesh> heit(hef.at(i), m), heend(heit) ;
		do {
			flip[*heit] = true ;
			++heit ;
		} while(heit != heend) ;
	}
	// Parcours de 1/2 arête et insertion dans res des 1/2 arêtes marquées
	Mesh::Halfedge_range her(m.halfedges()) ;
	for(Mesh::Halfedge_range::iterator it = her.begin(); it != her.end(); ++it)
	{
		res.push_back(*it) ;
	}
	m.remove_property_map(flip) ;
	return res ;
}

vector<Mesh::Vertex_index> flip_edges(vector<Mesh::Halfedge_index> &hef, Mesh &m)
{
	// Récupération des 1/2 arêtes à flipper
	vector<Mesh::Halfedge_index> res = edges_from_faces (hef, m) ;
	cout << "nombre de he a flipper " << res.size() << endl ;
	
	// Ajout d'un booléen sur les sommets pour marquage de flip
	Mesh::Property_map<Mesh::Vertex_index,bool> flip;
	bool created_flip;
	boost::tie(flip, created_flip) = m.add_property_map<Mesh::Vertex_index,bool>("v:flag",false);
	assert(created_flip);
	
	// Flip des arêtes pour
	for (int i=0; i<res.size(); ++i)
	{
		Mesh::Halfedge_index he = res.at(i);
		if (test_flip(he, m))
		{
			flip[m.source(he)] = true ;
			flip[m.target(he)] = true ;
			
			CGAL::Euler::flip_edge(res.at(i), m) ;
			
			flip[m.source(he)] = true ;
			flip[m.target(he)] = true ;
		}
	}
	// Liste des sommets touchés
	Mesh::Vertex_range r = m.vertices() ;
	vector<Mesh::Vertex_index> vlist ;
	for(Mesh::Vertex_range::iterator it = r.begin(); it != r.end(); ++it)
	{
		if (flip[*it])
			vlist.push_back(*it) ;
	}
	m.remove_property_map(flip) ;
	return vlist ;
}

bool test_flip(Mesh::Halfedge_index he, Mesh &m)
{
	/*      P
	  /     |      \
	 /      |       \
 Pa /alpha  |he  beta\ Pb
	\       |        /
	 \      |       /
	        PP
	 Teste si alpha + beta > pi -> he doit être échangée
	 */
	Mesh::Halfedge_index hea, heb ;
	heb = m.next(he) ;
	hea = m.next(m.opposite(he)) ;
	K::Point_3 Pa,P, PP, Pb ;
	P = m.point(m.target(he)) ;
	PP = m.point(m.source(he)) ;
	Pa = m.point(m.target(hea)) ;
	Pb = m.point(m.target(heb)) ;
	double alpha, beta ;
	alpha = CGAL::approximate_angle(P, Pa, PP);
	beta = CGAL::approximate_angle(PP, Pb, P) ;
	return ((alpha+beta) > 180. ) ;
}

K::Vector_3 laplac_discr(Mesh::Vertex_index v, Mesh &m)
{
	CGAL::Vertex_around_target_circulator<Mesh> vit(m.halfedge(v), m), vend(vit) ;
	K::Point_3 zero(0,0,0), P ;
	K::Vector_3 Psum(0,0,0) ;
	double cpt = 0. ;
	do
	{
		Psum += (m.point(*vit)-zero) ;
		++vit ;
		++cpt ;
	} while (vit != vend) ;
	Psum = Psum / cpt ;
	P = zero+Psum ;
	return P-m.point(v) ;
}

void laplac_smooth(vector<Mesh::Vertex_index> l, double lambda, Mesh &m)
{
	K::Vector_3 lapl ;
	cout << "lissage laplacien de : " << l.size() << " sommets" << endl ;
	for(int i=0; i<l.size(); ++i)
	{
		lapl = laplac_discr(l.at(i), m) ;
		m.point(l.at(i)) += lambda * lapl ;
	}
}
