//
//  hole_filling.h
//  test1
//
//  Created by Bac Alexandra on 12/04/2022.
//

#ifndef hole_filling_h
#define hole_filling_h

#include <iostream>
#include <cmath>
#include <Eigen/Core>
#include <Eigen/Geometry>
#include <Eigen/Dense>
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Cartesian_converter.h>
#include <CGAL/Surface_mesh_default_triangulation_3.h>
#include <CGAL/Complex_2_in_triangulation_3.h>
#include <CGAL/make_surface_mesh.h>
#include <CGAL/Implicit_surface_3.h>
#include <CGAL/IO/facets_in_complex_2_to_triangle_mesh.h>
#include <CGAL/boost/graph/IO/polygon_mesh_io.h>
#include <CGAL/IO/Color.h>
#include <CGAL/Surface_mesh.h>
#include <fstream>
#include "starter_defs.h"
#include "courbures.h"
#include "mesh_functions.hpp"

// default triangulation for Surface_mesher
typedef CGAL::Surface_mesh_default_triangulation_3 Tr;
// c2t3
typedef CGAL::Complex_2_in_triangulation_3<Tr> C2t3;
typedef Tr::Geom_traits GT;
typedef GT::Sphere_3 Sphere_3;
typedef GT::Point_3 Point_3;
typedef GT::FT FT;
typedef FT (*Function)(Point_3);
typedef CGAL::Implicit_surface_3<GT, Function> Surface_3;
typedef CGAL::Surface_mesh<Point_3> Surface_mesh;


using namespace std ;

// ******************************
// Function computing polynomial P_i (degree 2 polynomials in 3 variables)
inline float myp(int i, K::Point_3 X) {
	float x = X[0] ;
	float y = X[1] ;
	float z = X[2] ;
	switch (i) {
	case 0:
		return x*x ;
		break;
	case 1:
		return y*y ;
		break;
	case 2:
		return z*z ;
		break;
	case 3:
		return x*y ;
		break;
	case 4:
		return y*z ;
		break;
	case 5:
		return x*z ;
		break;
	case 6:
		return x ;
		break;
	case 7:
		return y ;
		break;
	case 8:
		return z;
		break;
	case 9:
		return 1;
		break;
	default:
		throw std::runtime_error("Error on indice i : unknown polynomial P.");
	}
}

// ******************************
// Function computing thin spline RBF
inline float myphi(float r) { if (r == 0) return 0 ; else return r*r*log(r) ; }

// ******************************
// Class encoding an implicit function built from the basis of RBF
// f(X) = Sum_i=0^n-1 alpha_i * phi(|| X-c_i ||) + Sum_j=0^9 beta_j * P_j(X)
class Implicit_RBF {
private:
	int _n ;
	int _d ;
	vector<double> _alpha, _beta ;
	vector<K::Point_3> _center ;

public:
	Implicit_RBF (vector<double> alpha, vector<double> beta, vector<K::Point_3> center) ;
    Implicit_RBF (const Implicit_RBF& other) : _n(other._n), _d(other._d), _alpha(other._alpha), _beta(other._beta), _center(other._center) {};
    Implicit_RBF& operator=(const Implicit_RBF& other) {
        _n = other._n;
        _d = other._d;
        _alpha = other._alpha;
        _beta = other._beta;
        _center = other._center;
    }
	double val(K::Point_3) const ;
};

// ******************************
// Class based on MyMesh_Norm encoding mesh with vertex boolean property
class MyMesh_Prop
{
public:
    MyMesh_Norm & _mn ;
    Mesh::Property_map<vertex_descriptor,bool> vprop; // Propriété de marquage du bord
    
    

    MyMesh_Prop(MyMesh_Norm &m) : _mn(m)
    {
        // Allocation de _vprop
        bool created_vprop;
        boost::tie(vprop, created_vprop) = _mn._m.add_property_map<vertex_descriptor,bool>("v:prop",false);
        assert(created_vprop);
    }
    
    ~MyMesh_Prop()
    {
        _mn._m.remove_property_map(vprop) ;
    }
} ;


// ******************************
// Class encoding all the functions required for implicit hole filling
class Hole_Filling
{
protected:
	static bool f_true (Mesh::Vertex_index vh) { return true ;}
	static std::function<bool(Mesh::Vertex_index)> _f_select ;
    
    // Propriété de marquage du bord
    Mesh::Property_map<vertex_descriptor,bool> vprop;
    
    void colorize_prop() ;
    void colorize_verts(const vector<Mesh::Vertex_index> &vlist) ;


    // Computation of boundary and its neighborhood
    Mesh::Halfedge_index find_boundary_edge() ;
    vector<Mesh::Vertex_index> find_boundary(Mesh::Halfedge_index heh) ;
    
    void init_mark_boundary(const vector<Mesh::Vertex_index> & bnd) ;
    vector<Mesh::Vertex_index> next_neighbors(const vector<Mesh::Vertex_index> & bnd) ;
    
public:
	MyMesh_Norm& _mesh ;
	
	Hole_Filling(MyMesh_Norm &mesh);
    inline ~Hole_Filling() {
        std::cout << "End hole filling" << std::endl ;
        _mesh._m.remove_property_map(vprop);
    }

	// Computation of RBF
	pair<pair<Eigen::MatrixXd &,Eigen::VectorXd &>,vector<K::Point_3> &> compute_approx_mat(
			vector<Mesh::Vertex_index> vlist, double scale) ;
	pair<vector<double>&, vector<double>&> solve_approx(const pair<Eigen::MatrixXd &, Eigen::VectorXd &> &p, int n, int d) ;

	// Hole filling
    std::shared_ptr<Implicit_RBF> fill_hole(std::vector<Mesh::Vertex_index>& boundary);

	// IO
	GT::Sphere_3 estimate_BS(const vector<Mesh::Vertex_index> &vlist) ;
    void output_mesh_implicit(std::shared_ptr<Implicit_RBF> implicit, const vector<Mesh::Vertex_index>& boundary, string out_file);
	inline void print_normales()
	{
		Mesh::Vertex_range vr = _mesh._m.vertices() ;
		for (Mesh::Vertex_range::iterator vi = vr.begin(); vi != vr.end(); ++vi)
		{
			K::Vector_3 normale(_mesh.vnormal[*vi]) ;
			cout << "normale " << *vi << " : " << normale.x() << " " << normale.y() << " " << normale.z() << endl ;
		}
	}
};




// Other ....

inline K::Point_3 min (K::Point_3 X, K::Point_3 Y) { double P[3] ; for (int i=0; i<3; i++) P[i] = min(X[i], Y[i]) ; return K::Point_3(P[0], P[1], P[2]) ;}
inline K::Point_3 max (K::Point_3 X, K::Point_3 Y) { double P[3] ; for (int i=0; i<3; i++) P[i] = max(X[i], Y[i]) ; return K::Point_3(P[0], P[1], P[2]) ;}

// Calcul du nom d'export en fonction du nom du fichier ouvert
inline std::string to_out_file_name(string out)
{
	unsigned int ext_pos = out.find_last_of(".") ;
	out = out.substr(0, ext_pos) ;
	out = out + "_patch.obj" ;

	std::cout << out << std::endl ;
	return out ;
}

// ///////

//FT implicit_function [&implicit](Point_3 p) {
//	return implicit.val(p) ;
//}

#endif /* hole_filling_hpp */
