delaunay2d/delaunay.h

/**
 * @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"

#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;
}

void circumcircle(vertex2dc *v0, vertex2dc *v1, vertex2dc *v2, double &cx, double &cy, double &cr) // calculate the circumcircle from three points
{
	double s = 0.5 / ((v1->x - v0->x) * (v2->y - v0->y) - (v1->y - v0->y) * (v2->x - v0->x));
	double m = v1->x*v1->x - v0->x*v0->x + v1->y*v1->y - v0->y*v0->y;
	double u = v2->x*v2->x - v0->x*v0->x + v2->y*v2->y - v0->y*v0->y;

	cx = ((v2->y - v0->y)*m + (v0->y - v1->y)*u)*s;
	cy = ((v0->x - v2->x)*m + (v1->x - v0->x)*u)*s;
	cr = (v0->x - cx)*(v0->x - cx) + (v0->y - cy)*(v0->y - cy); // not need to calculate the squared root here
	return;
}
// End vertex definition

// Start triangle definition
struct triangle
{
	int id;
	vertex2dc *vert[3]; // vertex of the triangle
	triangle *neigh[3]; // neighbors 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; neigh[0] = neigh[1] = neigh[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;
		neigh[0] = neigh[1] = neigh[2] = nullptr;

		circumcircle(vert[0], vert[1], vert[2], cx, cy, cr);
		return;
	}

	void set_neighbor(triangle *n0ptr, triangle *n1ptr, triangle *n2ptr)
	{
		neigh[0] = n0ptr; neigh[1] = n1ptr; neigh[2] = n2ptr;
		return;
	}

	bool bound_location(double inx, double iny) // Test if the location is inside the triangle
	{
		double l1x, l1y, l2x, l2y;
		for (int i = 0; i < 3; i++)
		{
			l1x = vert[(i+1)%3]->x - vert[i]->x;
			l1y = vert[(i+1)%3]->y - vert[i]->y;
			l2x = inx - vert[i]->x;
			l2y = iny - vert[i]->y;

			if ((l1x*l2y - l1y*l2x) < 0) // This condition includes points on the triangle's edge
			{
				return false;
			}
		}
		return true;
	}
};
// 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)
{
	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); // lower left corner
	assit_vert.push_back(tmp_vert);

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

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

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

	triangle *old_tri = nullptr, *tmp_tri = nullptr;
	triangle *cnst_tri[3], *old_neigh[6];
	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
	out_tris.push_back(tmp_tri);

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

	out_tris[0]->set_neighbor(nullptr, nullptr, out_tris[1]);
	out_tris[1]->set_neighbor(out_tris[0], nullptr, nullptr);

	// loop all input vertice
	double dist;
	for (int i = 0; i < in_verts.size(); ++i)
	{
		// determine the triangle that includes the new vertex and remove it from out_tris
		for (t_iter = out_tris.begin(); t_iter != out_tris.end(); )
		{
			old_tri = *t_iter;
			if (old_tri->bound_location(in_verts[i].x, in_verts[i].y))
			{
				t_iter = out_tris.erase(t_iter);
				break;
			}
			else t_iter++;
		}

		// build three new triangles
		for (int n = 0; n < 3; ++n)
		{
			tmp_tri = new triangle(old_tri->vert[n], old_tri->vert[(n+1)%3], &in_verts[i]);
			cnst_tri[n] = tmp_tri;
			out_tris.push_back(tmp_tri);
		}

		// sort neighbors
		for (int n = 0; n < 3; ++n)
		{
			if (old_tri->neigh[n] == nullptr)
			{
				cnst_tri[n]->set_neighbor(nullptr, cnst_tri[(n+1)%3], cnst_tri[(n+2)%3]);
			}
			else
			{
				cnst_tri[n]->set_neighbor(old_tri->neigh[n], cnst_tri[(n+1)%3], cnst_tri[(n+2)%3]);
				for (int k = 0; k < 3; ++k) // replace neighbor for the oppositing triangle 
				{
					if (old_tri->neigh[n]->neigh[k] == old_tri)
					{
						old_tri->neigh[n]->neigh[k] = cnst_tri[n];
						break;
					}
				}
			}
		}

		// delete the old triangle
		delete old_tri; old_tri = nullptr;

		// test if the cnst_tri need to be flipped
		for (int n = 0; n < 3; ++n)
		{
			if (cnst_tri[n]->neigh[0] != nullptr) // must has neighbor on this side
			{
				old_tri = cnst_tri[n]->neigh[0];
				for (int v = 0; v < 3; ++v)
				{
					tmp_vert = old_tri->vert[v];
					if (tmp_vert != cnst_tri[n]->vert[0] && tmp_vert != cnst_tri[n]->vert[1]) // find the opposite vertex
					{
						//dist = (cnst_tri[n]->cx - tmp_vert->x) * (cnst_tri[n]->cx - tmp_vert->x) + 
						//	(cnst_tri[n]->cy - tmp_vert->y) * (cnst_tri[n]->cy - tmp_vert->y);

						//if ((dist - cnst_tri[n]->cr) <= ZERO) // need to be flipped

						dist = (old_tri->cx - cnst_tri[n]->vert[2]->x) * (old_tri->cx - cnst_tri[n]->vert[2]->x) + 
							(old_tri->cy - cnst_tri[n]->vert[2]->y) * (old_tri->cy - cnst_tri[n]->vert[2]->y);
						
						if ((dist - old_tri->cr) <= ZERO) // need to be flipped
						{
							// record the original neighbors
							old_neigh[0] = cnst_tri[n]->neigh[0];
							old_neigh[1] = cnst_tri[n]->neigh[1];
							old_neigh[2] = cnst_tri[n]->neigh[2];
							old_neigh[3] = old_tri->neigh[0];
							old_neigh[4] = old_tri->neigh[1];
							old_neigh[5] = old_tri->neigh[2];

							cnst_tri[n]->set(cnst_tri[n]->vert[0], tmp_vert, cnst_tri[n]->vert[2]); // flip
						
							if (v == 0)
							{
								old_tri->set(old_tri->vert[0], old_tri->vert[1], cnst_tri[n]->vert[2]); //flip

								// Sort neighbors
								cnst_tri[n]->set_neighbor(old_neigh[5], old_tri, old_neigh[2]); 
								old_tri->set_neighbor(old_neigh[3], old_neigh[1], cnst_tri[n]);

								for (int k = 0; k < 3; ++k)
								{
									if (cnst_tri[n]->neigh[0] != nullptr && cnst_tri[n]->neigh[0]->neigh[k] == old_tri)
									{
										cnst_tri[n]->neigh[0]->neigh[k] = cnst_tri[n];
										break;
									}
								}
								
								for (int k = 0; k < 3; ++k)
								{
									if (old_tri->neigh[1] != nullptr && old_tri->neigh[1]->neigh[k] == cnst_tri[n])
									{
										old_tri->neigh[1]->neigh[k] = old_tri;
										break;
									}
								}
							}
							else if (v == 1)
							{
								old_tri->set(cnst_tri[n]->vert[2], old_tri->vert[1], old_tri->vert[2]); //flip

								// Sort neighbors
								cnst_tri[n]->set_neighbor(old_neigh[3], old_tri, old_neigh[2]); 
								old_tri->set_neighbor(cnst_tri[n], old_neigh[4], old_neigh[1]);

								for (int k = 0; k < 3; ++k)
								{
									if (cnst_tri[n]->neigh[0] != nullptr && cnst_tri[n]->neigh[0]->neigh[k] == old_tri)
									{
										cnst_tri[n]->neigh[0]->neigh[k] = cnst_tri[n];
										break;
									}
								}
								
								for (int k = 0; k < 3; ++k)
								{
									if (old_tri->neigh[2] != nullptr && old_tri->neigh[2]->neigh[k] == cnst_tri[n])
									{
										old_tri->neigh[2]->neigh[k] = old_tri;
										break;
									}
								}
							}
							else
							{
								old_tri->set(old_tri->vert[0], cnst_tri[n]->vert[2], old_tri->vert[2]); //flip

								// Sort neighbors
								cnst_tri[n]->set_neighbor(old_neigh[4], old_tri, old_neigh[2]); 
								old_tri->set_neighbor(old_neigh[1], cnst_tri[n], old_neigh[5]);

								for (int k = 0; k < 3; ++k)
								{
									if (cnst_tri[n]->neigh[0] != nullptr && cnst_tri[n]->neigh[0]->neigh[k] == old_tri)
									{
										cnst_tri[n]->neigh[0]->neigh[k] = cnst_tri[n];
										break;
									}
								}
								
								for (int k = 0; k < 3; ++k)
								{
									if (old_tri->neigh[0] != nullptr && old_tri->neigh[0]->neigh[k] == cnst_tri[n])
									{
										old_tri->neigh[0]->neigh[k] = old_tri;
										break;
									}
								}
							}
						}
						break;
					}
				}
			}
		}
	}

	// remove any triangles has an assistant vertex from out_tris
	for (t_iter = out_tris.begin(); t_iter != out_tris.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])
		{
			for (int k = 0; k < 3; ++k)
			{
				if (tmp_tri->neigh[1] != nullptr && tmp_tri->neigh[1]->neigh[k] == tmp_tri)
				{
					tmp_tri->neigh[1]->neigh[k] = nullptr;
					break;
				}
			}

			// destroy the memories located and remove from the vector
			t_iter = out_tris.erase(t_iter);
			delete tmp_tri; tmp_tri = nullptr;
		}
		else if (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])
		{
			for (int k = 0; k < 3; ++k)
			{
				if (tmp_tri->neigh[2] != nullptr && tmp_tri->neigh[2]->neigh[k] == tmp_tri)
				{
					tmp_tri->neigh[2]->neigh[k] = nullptr;
					break;
				}
			}

			// destroy the memories located and remove from the vector
			t_iter = out_tris.erase(t_iter);
			delete tmp_tri; tmp_tri = nullptr;
		}
		else if (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])
		{
			for (int k = 0; k < 3; ++k)
			{
				if (tmp_tri->neigh[0] != nullptr && tmp_tri->neigh[0]->neigh[k] == tmp_tri)
				{
					tmp_tri->neigh[0]->neigh[k] = nullptr;
					break;
				}
			}

			// destroy the memories located and remove from the vector
			t_iter = out_tris.erase(t_iter);
			delete tmp_tri; tmp_tri = nullptr;
		}
		else t_iter++;
	}

	// assign triangles index
	for (int i = 0; i < out_tris.size(); i++)
	{
		out_tris[i]->id = i;
	}
	
	// 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])
			{
				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 false;

	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