/**
 * @defgroup   DELAUNAY
 *
 * @brief      An implementation of the 2D Delaunay triangulation using the Bowyer-Watson algorithm.
 *
 * @author     Yi Zhang
 * @date       2021-09-12
 */

#ifndef _BW_2D_DELAUNAY_H
#define _BW_2D_DELAUNAY_H
#include "cmath"
#include "vector"

#include "iostream"
#include "fstream"
#include "iomanip"

#define ZERO 1e-5

// Start vertex definition
struct vertex2dc
{
	unsigned int id; // index of the vertex
	double x, y; // position of the vertex

	vertex2dc() : x(NAN), y(NAN), id(0) {}
	vertex2dc(double inx, double iny, unsigned int inid = 0) {set(inx, iny, inid);}
	void set(double inx, double iny, unsigned int inid = 0)
	{
		x = inx; y = iny; id = inid;
		return;
	}
};

bool operator ==(const vertex2dc &a, const vertex2dc &b) // overload the == operator for vertex2dc type
{
	if(fabs(a.x - b.x) <= ZERO && fabs(a.y - b.y) <= ZERO)
	{
		return true;
	}
	return false;
}
// End vertex definition

// Start edge definition
struct edge
{
	vertex2dc *vert[2]; // vertex of the edge

	edge() {vert[0] = vert[1] = nullptr;}
	edge(vertex2dc *v0ptr, vertex2dc *v1ptr) {set(v0ptr, v1ptr);}
	void set(vertex2dc *v0ptr, vertex2dc *v1ptr)
	{
		vert[0] = v0ptr; vert[1] = v1ptr;
		return;
	}
};

bool operator ==(const edge &a, const edge &b) // overload the == operator for edge type
{
	if((a.vert[0] == b.vert[0] && a.vert[1] == b.vert[1]) || 
		(a.vert[0] == b.vert[1] && a.vert[1] == b.vert[0]))
	{
		return true;
	}
	return false;
}
// End edge definition

// Start triangle definition
struct triangle
{
	vertex2dc *vert[3]; // vertex of the triangle
	double cx, cy; // center of the triangle's circumcircle
	double cr; // radius of the circumcircle

	triangle() {vert[0] = vert[1] = vert[2] = nullptr;}
	triangle(vertex2dc *v0ptr, vertex2dc *v1ptr, vertex2dc *v2ptr) {set(v0ptr, v1ptr, v2ptr);}
	void set(vertex2dc *v0ptr, vertex2dc *v1ptr, vertex2dc *v2ptr)
	{
		vert[0] = v0ptr; vert[1] = v1ptr; vert[2] = v2ptr;

		double s = 0.5 / ((vert[1]->x - vert[0]->x) * (vert[2]->y - vert[0]->y) - (vert[1]->y - vert[0]->y) * (vert[2]->x - vert[0]->x));
		double m = vert[1]->x * vert[1]->x - vert[0]->x * vert[0]->x + vert[1]->y * vert[1]->y - vert[0]->y * vert[0]->y;
		double u = vert[2]->x * vert[2]->x - vert[0]->x * vert[0]->x + vert[2]->y * vert[2]->y - vert[0]->y * vert[0]->y;

		cx = ((vert[2]->y - vert[0]->y) * m + (vert[0]->y - vert[1]->y) * u) * s;
		cy = ((vert[0]->x - vert[2]->x) * m + (vert[1]->x - vert[0]->x) * u) * s;
		cr = (vert[0]->x - cx) * (vert[0]->x - cx) + (vert[0]->y - cy) * (vert[0]->y - cy); // not need to sqrt() here
		return;
	}
};
// End triangle definition

/**
 * @brief      2D Delaunay triangulation of some given points using the Bowyer-Watson algorithm.
 *
 * @param      in_verts  Input vertexes. Defined by the user.
 * @param      out_tris  Output triangles. Compute by the function.
 */
void triangulation(std::vector<vertex2dc> &in_verts, std::vector<triangle> &out_tris, int valid_size)
{
	if (!out_tris.empty()) out_tris.clear();
	if (in_verts.size() < 3) return;

	// locate the surrounding box and initiate the staring two triangles
	double xmin = in_verts[0].x, xmax = in_verts[0].x;
	double ymin = in_verts[0].y, ymax = in_verts[0].y;
	for (int i = 0; i < in_verts.size(); ++i)
	{
		xmin = std::min(xmin, in_verts[i].x);
		xmax = std::max(xmax, in_verts[i].x);
		ymin = std::min(ymin, in_verts[i].y);
		ymax = std::max(ymax, in_verts[i].y);
	}

	double midx = 0.5*(xmin + xmax);
	double midy = 0.5*(ymin + ymax);
	double maxi_s = std::max(xmax - xmin, ymax - ymin); // use an four times bigger rectangle to include all points

	vertex2dc *tmp_vert = nullptr;
	std::vector<vertex2dc*> assit_vert;

	tmp_vert = new vertex2dc(midx - maxi_s, midy - maxi_s, 0); // lower left corner
	assit_vert.push_back(tmp_vert);

	tmp_vert = new vertex2dc(midx + maxi_s, midy - maxi_s, 1); // lower right corner
	assit_vert.push_back(tmp_vert);

	tmp_vert = new vertex2dc(midx + maxi_s, midy + maxi_s, 2); // upper right corner
	assit_vert.push_back(tmp_vert);

	tmp_vert = new vertex2dc(midx - maxi_s, midy + maxi_s, 3); // upper left corner
	assit_vert.push_back(tmp_vert);

	triangle *tmp_tri = nullptr;
	std::vector<triangle*>  exist_tri, cnst_tri;
	std::vector<triangle*>::iterator t_iter;

	tmp_tri = new triangle(assit_vert[0], assit_vert[1], assit_vert[2]); // order the vertex anti-clock wise
	exist_tri.push_back(tmp_tri);

	tmp_tri = new triangle(assit_vert[0], assit_vert[2], assit_vert[3]); // order the vertex anti-clock wise
	exist_tri.push_back(tmp_tri);

	// loop all input vertice
	bool removed;
	double dist;
	edge tmp_edge;
	std::vector<edge> cnst_edge;
	std::vector<edge>::iterator e_iter;
	//for (int i = 0; i < in_verts.size(); ++i)
	for (int i = 0; i < valid_size; ++i)
	{
		// determine triangles that include the point and add the triangle to the cnst_tri and remove the triangle from exist_tri
		// this is the part that could take a lot of time if we are working with a large amount of points. We will fix this later
		cnst_tri.clear();
		for (t_iter = exist_tri.begin(); t_iter != exist_tri.end(); )
		{
			tmp_tri = *t_iter;

			dist = (tmp_tri->cx - in_verts[i].x) * (tmp_tri->cx - in_verts[i].x) + 
				(tmp_tri->cy - in_verts[i].y) * (tmp_tri->cy - in_verts[i].y);

			if ((dist - tmp_tri->cr) <= ZERO) // Points on the circumcircle are also included
			{
				t_iter = exist_tri.erase(t_iter);
				cnst_tri.push_back(tmp_tri);
			}
			else t_iter++;
		}

		// loop to remove duplicate edges
		cnst_edge.clear();
		for (int c = 0; c < cnst_tri.size(); ++c)
		{
			for (int e = 0; e < 3; ++e)
			{
				tmp_edge.set(cnst_tri[c]->vert[e], cnst_tri[c]->vert[(e+1)%3]);

				removed = false;
				for (e_iter = cnst_edge.begin(); e_iter != cnst_edge.end(); )
				{
					if (tmp_edge == *e_iter) // duplicate edge, remove from cnst_edge
					{
						e_iter = cnst_edge.erase(e_iter);
						removed = true;
						break; // no need to search more
					}
					else e_iter++;
				}

				if (!removed) // not a duplicate edge, add to the cnst_edge
				{
					cnst_edge.push_back(tmp_edge);
				}
			}
		}

		// construct new triangles and add to exist_tri
		for (int c = 0; c < cnst_edge.size(); ++c)
		{
			tmp_tri = new triangle(cnst_edge[c].vert[0], cnst_edge[c].vert[1], &in_verts[i]); // order the vertex anti-clock wise
			exist_tri.push_back(tmp_tri);
		}

		// destroy memories used by cnst_edge
		for (int c = 0; c < cnst_tri.size(); ++c)
		{
			tmp_tri = cnst_tri[c];
			delete tmp_tri; tmp_tri = nullptr;
		}
	}

	// Start testing code
	std::ofstream outfile("test_backup.msh");
	outfile << "$MeshFormat" << std::endl << "2.2 0 8" << std::endl << "$EndMeshFormat "<<std::endl;
	outfile << "$Nodes" << std::endl << valid_size + 4 << std::endl;
	for (int i = 0; i < 4; i++)
	{
		outfile << i + 1 << " " << std::setprecision(16) 
			<< assit_vert[i]->x << " " << assit_vert[i]->y << " 0.0" << std::endl;
	}
	for (int i = 0; i < valid_size; i++)
	{
		outfile << i + 5 << " " << std::setprecision(16) 
			<< in_verts[i].x << " " << in_verts[i].y << " 0.0" << std::endl;
	}
	outfile<<"$EndNodes"<<std::endl;
	outfile << "$Elements" << std::endl << exist_tri.size() <<std::endl;
	for (int i = 0; i < exist_tri.size(); i++)
	{
		outfile << i + 1 << " 2 0";
		for (int j = 0; j < 3; j++)
		{
			outfile << " " << exist_tri[i]->vert[j]->id + 1;
		}
		outfile << std::endl;
	}
	outfile << "$EndElements"<< std::endl;
	outfile.close();
	// End testing code

	// remove any triangles has an assistant vertex from exist_tri
	for (t_iter = exist_tri.begin(); t_iter != exist_tri.end(); )
	{
		tmp_tri = *t_iter;
		if (tmp_tri->vert[0] == assit_vert[0] || tmp_tri->vert[0] == assit_vert[1] || tmp_tri->vert[0] == assit_vert[2] || tmp_tri->vert[0] == assit_vert[3] || 
			tmp_tri->vert[1] == assit_vert[0] || tmp_tri->vert[1] == assit_vert[1] || tmp_tri->vert[1] == assit_vert[2] || tmp_tri->vert[1] == assit_vert[3] || 
			tmp_tri->vert[2] == assit_vert[0] || tmp_tri->vert[2] == assit_vert[1] || tmp_tri->vert[2] == assit_vert[2] || tmp_tri->vert[2] == assit_vert[3])
		{
			// destroy the memories located and remove from the vector
			t_iter = exist_tri.erase(t_iter);
			delete tmp_tri; tmp_tri = nullptr;
		}
		else t_iter++;
	}

	// copy exist_tri to out_tris and destroy memories located
	out_tris.resize(exist_tri.size());
	for (int i = 0; i < exist_tri.size(); ++i)
	{
		out_tris[i].set(exist_tri[i]->vert[0], exist_tri[i]->vert[1], exist_tri[i]->vert[2]);
		delete exist_tri[i]; exist_tri[i] = nullptr;
	}

	// destroy memories located for assit_vert
	for (int i = 0; i < 4; ++i)
	{
		delete assit_vert[i]; assit_vert[i] = nullptr;
	}
	return;
}

/**
 * @brief      Check for duplicated vertex
 *
 * @param[in]  in_verts  Input vertexes
 *
 * @return     If there is duplicated vertex
 */
bool duplicated_vertex(const std::vector<vertex2dc> &in_verts)
{
	if (in_verts.empty()) return false;

	for (int i = 0; i < in_verts.size()-1; ++i)
	{
		for (int j = i+1; j < in_verts.size(); ++j)
		{
			if (in_verts[i] == in_verts[j] && in_verts[i].id != in_verts[j].id)
			{
				return true;
			}
		}
	}
	return false;
}

/**
 * @brief      Check to see if the triangulation is fully delaunay
 *
 * @param[in]  in_tris   Input triangles
 * @param[in]  in_verts  Input vertexes
 *
 * @return     If the triangulation is fully delaunay
 */
bool fully_delaunay(const std::vector<triangle> &in_tris, const std::vector<vertex2dc> &in_verts)
{
	if (in_tris.empty()) return true;

	int count;
	double dist;
	for (int i = 0; i < in_tris.size(); ++i)
	{
		count = 0;
		for (int j = 0; j < in_verts.size(); ++j)
		{
			dist = (in_tris[i].cx - in_verts[j].x) * (in_tris[i].cx - in_verts[j].x) + 
				(in_tris[i].cy - in_verts[j].y) * (in_tris[i].cy - in_verts[j].y);

			if ((dist - in_tris[i].cr) <= ZERO)
			{
				count++;
			}
		}

		if (count > 3)
		{
			return false;
		}
	}

	return true;
}

#endif // _BW_2D_DELAUNAY_H