#define CGAL_USE_BASIC_VIEWER

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>

#include <CGAL/Constrained_Delaunay_triangulation_2.h>
#include <CGAL/Triangulation_face_base_with_info_2.h>
#include <CGAL/Triangulation_vertex_base_with_info_2.h>

#include <CGAL/Delaunay_mesher_2.h>
#include <CGAL/Delaunay_mesh_face_base_2.h>
#include <CGAL/Delaunay_mesh_vertex_base_2.h>
#include <CGAL/Delaunay_mesh_size_criteria_2.h>
#include <CGAL/lloyd_optimize_mesh_2.h>
#include <CGAL/Polygon_2.h>

#include <iostream>
#include <cstdint>
#include <fstream>
#include <iostream>

#include "mex.hpp"
#include "mexAdapter.hpp"
#include "MatlabDataArray.hpp"
#include <omp.h>

struct FaceInfo2
{
	FaceInfo2() {}
	int nesting_level;
	bool in_domain() {
		return nesting_level % 2 == 1;
	}
};

typedef CGAL::Exact_predicates_inexact_constructions_kernel	K;

typedef K::Point_2											Point;
//typedef CGAL::Triangulation_vertex_base_2<K>                      Vb1;
typedef CGAL::Triangulation_vertex_base_with_info_2<unsigned int, K> Vb1;
typedef CGAL::Triangulation_face_base_with_info_2<FaceInfo2, K>    Fbb;
typedef CGAL::Constrained_triangulation_face_base_2<K, Fbb>        Fb1;
typedef CGAL::Triangulation_data_structure_2<Vb1, Fb1>               TDS1;
typedef CGAL::Exact_predicates_tag                                Itag;
typedef CGAL::Constrained_Delaunay_triangulation_2<K, TDS1, Itag>  CDT1;

typedef CGAL::Delaunay_mesh_vertex_base_2<K>                Vb;
typedef CGAL::Delaunay_mesh_face_base_2<K>                  Fb;
typedef CGAL::Triangulation_data_structure_2<Vb, Fb>        TDS;
typedef CGAL::Constrained_Delaunay_triangulation_2<K, TDS>  CDT;

typedef CDT1::Point                                                Point;
typedef CGAL::Polygon_2<K>                                        Polygon2;
typedef CDT1::Face_handle                                          Face_handle;

typedef CGAL::Delaunay_mesh_size_criteria_2<CDT>            Criteria;
typedef CGAL::Delaunay_mesher_2<CDT, Criteria>              Mesher;

void mark_domains(CDT1& ct, Face_handle start, int index, std::list<CDT1::Edge>& border)
{
	if (start->info().nesting_level != -1) {
		return;
	}
	std::list<Face_handle> queue;
	queue.push_back(start);
	while (!queue.empty()) {
		Face_handle fh = queue.front();
		queue.pop_front();
		if (fh->info().nesting_level == -1) {
			fh->info().nesting_level = index;
			for (int i = 0; i < 3; i++) {
				CDT1::Edge e(fh, i);
				Face_handle n = fh->neighbor(i);
				if (n->info().nesting_level == -1) {
					if (ct.is_constrained(e)) border.push_back(e);
					else queue.push_back(n);
				}
			}
		}
	}
}
//explore set of facets connected with non constrained edges,
//and attribute to each such set a nesting level.
//We start from facets incident to the infinite vertex, with a nesting
//level of 0. Then we recursively consider the non-explored facets incident
//to constrained edges bounding the former set and increase the nesting level by 1.
//Facets in the domain are those with an odd nesting level.
void mark_domains(CDT1& cdt)
{
	for (CDT1::Face_handle f : cdt.all_face_handles()) {
		f->info().nesting_level = -1;
	}
	std::list<CDT1::Edge> border;
	mark_domains(cdt, cdt.infinite_face(), 0, border);
	while (!border.empty()) {
		CDT1::Edge e = border.front();
		border.pop_front();
		Face_handle n = e.first->neighbor(e.second);
		if (n->info().nesting_level == -1) {
			mark_domains(cdt, n, e.first->info().nesting_level + 1, border);
		}
	}
}

class MexFunction : public matlab::mex::Function {
public:
	void operator()(matlab::mex::ArgumentList outputs, matlab::mex::ArgumentList inputs) 
	{
		checkArguments(outputs, inputs);

		matlab::data::TypedArray<double> pnts = std::move(inputs[0]);
		double x = inputs[1][0];

		matlab::data::CellArray nodeMat = mexDT(pnts, x);
		outputs[0] = nodeMat;
	}

	void checkArguments(matlab::mex::ArgumentList outputs, matlab::mex::ArgumentList inputs) {
        std::shared_ptr<matlab::engine::MATLABEngine> matlabPtr = getEngine();
        matlab::data::ArrayFactory factory;

        if (inputs.size() != 2) {
            matlabPtr->feval(u"error", 
                0, std::vector<matlab::data::Array>({ factory.createScalar("Four inputs required") }));
        }

        if (inputs[0].getDimensions()[1] != 2) {
            matlabPtr->feval(u"error", 
                0, std::vector<matlab::data::Array>({ factory.createScalar("Points input 1 must be NX2 dimension") }));
        }
        if (inputs[0].getType() != matlab::data::ArrayType::DOUBLE ||
            inputs[0].getType() == matlab::data::ArrayType::COMPLEX_DOUBLE) {
            matlabPtr->feval(u"error", 
                0, std::vector<matlab::data::Array>({ factory.createScalar("Points Input 1 must be a noncomplex scalar double") }));
        }

		if (inputs[1].getDimensions()[0] != 1 || inputs[1].getDimensions()[1] != 1) {
            matlabPtr->feval(u"error", 0, 
				std::vector<matlab::data::Array>({ factory.createScalar("input 2 must be 1X1 dimension") }));
        }       
		if (inputs[1].getType() != matlab::data::ArrayType::DOUBLE ||
			inputs[1].getType() == matlab::data::ArrayType::COMPLEX_DOUBLE) {
			matlabPtr->feval(u"error", 
				0, std::vector<matlab::data::Array>({ factory.createScalar("input 2 must be a noncomplex scalar double") }));
		}
		if (outputs.size() > 1) {
            matlabPtr->feval(u"error", 0, 
                std::vector<matlab::data::Array>({ factory.createScalar("Only one output is returned") }));
        }
    }

	matlab::data::CellArray mexDT(matlab::data::TypedArray<double> pnts, double x)
	{
		// step.0 generate polygon
		Polygon2 pgon;
		for (int i = 0; i < pnts.getDimensions()[0]; ++i)
			pgon.push_back(Point(pnts[i][0], pnts[i][1]));
		 
		// step.1 lloyd opt
		CDT cdt;
		cdt.insert_constraint(pgon.vertices_begin(), pgon.vertices_end(), true);

		Mesher mesher(cdt);
		mesher.set_criteria(Criteria(0.125, x));
		mesher.refine_mesh();

		CGAL::lloyd_optimize_mesh_2(cdt, CGAL::parameters::max_iteration_number = 20);

		// step.2 transform to CDT1
		std::vector<std::pair<Point, unsigned>> points;
		int index = 0;
		for (CDT::Vertex_handle vh : cdt.finite_vertex_handles())
			points.push_back(std::make_pair(vh->point(), index++));

		CDT1 cdt1;
		cdt1.insert(points.begin(), points.end());
		cdt1.insert_constraint(pgon.vertices_begin(), pgon.vertices_end(), true);
		mark_domains(cdt1);

		// step.3 write
		Eigen::MatrixXd nodes = Eigen::MatrixXd(cdt1.number_of_vertices(), 2);
		std::vector<std::vector<int>> faces;

		for (CDT1::Finite_faces_iterator f = cdt1.finite_faces_begin(); f != cdt1.finite_faces_end(); f++)
		{
			// find the faces in polygon
			if (f->info().in_domain()) 
			{
				std::vector<int> face_i{0,0,0};
				for (int i = 0; i < 3; ++i)
				{
					Point p = f->vertex(i)->point();
					int idx = f->vertex(i)->info();
					nodes(idx, 0) = p.x();
					nodes(idx, 1) = p.y();
					face_i[i] = idx;
				}
				faces.push_back(face_i);
			}
		}

		matlab::data::ArrayFactory factory;
		matlab::data::TypedArray<double> nodesMat = factory.createArray<double>({ cdt1.number_of_vertices(), 2 });
		for (int i = 0; i < nodes.rows(); i++)
		{
			nodesMat[i][0] = nodes(i, 0);
			nodesMat[i][1] = nodes(i, 1);
		}

		matlab::data::TypedArray<double> facesMat = factory.createArray<double>({ faces.size(), 3 });
		for (int i = 0; i < faces.size(); ++i)
		{
			facesMat[i][0] = faces[i][0] + 1;
			facesMat[i][1] = faces[i][1] + 1;
			facesMat[i][2] = faces[i][2] + 1;
		}

		matlab::data::CellArray trimesh = factory.createCellArray({ 2, 1 });
		trimesh[0] = nodesMat;
		trimesh[1] = facesMat;
		return trimesh;
	}
};