psemek/libs/cg/include/psemek/cg/voronoi.hpp

#pragma once

#include <psemek/geom/gauss.hpp>
#include <psemek/cg/triangulation/delaunay.hpp>
#include <psemek/cg/dual.hpp>
#include <psemek/cg/bbox.hpp>

namespace psemek::cg
{

	// Point #0 of the resulting dcel corresponds to the outer face,
	// i.e. the point at infinity, and has unspecified coordinates.
	template <typename Index = std::size_t, typename Iterator>
	auto voronoi(Iterator begin, Iterator end)
	{
		using point_type = std::decay_t<decltype(*begin)>;

		auto del = dual(delaunay<Index>(begin, end));

		dcel<point_type, util::empty, Iterator, Index> result;

		result.points.reserve(del.points.size());
		for (std::size_t i = 0; i < del.points.size(); ++i)
		{
			auto p = del.point(i);

			point_type q[3];

			auto e = p.edge();
			q[0] = *e.face().data();
			e = e.twin().next();
			q[1] = *e.face().data();
			e = e.twin().next();
			q[2] = *e.face().data();

			using T = typename point_type::scalar_type;

			geom::matrix<T, 3, 3> m;
			geom::vector<T, 3> b;

			for (std::size_t i = 0; i < 3; ++i)
			{
				m[i][0] = q[i][0];
				m[i][1] = q[i][1];
				m[i][2] = T{1};

				b[i] = q[i][0] * q[i][0] + q[i][1] * q[i][1];
			}

			geom::gauss(m, b);

			auto newp = result.push_point(point_type{b[0] / 2, b[1] / 2});
			newp.edge(result.edge(p.edge().index()));
		}

		result.edges = std::move(del.edges);
		result.faces = std::move(del.faces);

		return result;
	}

	// Replaces the #0 infinite point of a dcel with a proper outer face
	template <typename Point, typename Face, typename Edge, typename Index>
	void remove_infinite_point(dcel<Point, Face, Edge, Index> & dcel)
	{
		// insert outer face at position #0
		auto outer_face = dcel.insert_face(0);

		std::vector<std::size_t> delete_edges;
		std::vector<std::size_t> delete_faces;

		auto infinite_point = dcel.point(0);

		// iterate edges of infinite vertex #0 and fix hull
		for (auto e = infinite_point.edge();;)
		{
			for (auto n = e.next(); n.twin().origin() != infinite_point; n = n.next())
			{
				n.face(outer_face);
				n.origin().edge(n);
				if (!outer_face.edge())
					outer_face.edge(n);
			}

			auto n = e;
			for (; n.next() != e.twin(); n = n.next().twin());

			if (n.origin() != infinite_point)
			{
				auto t = e.next();
				for (; t.twin().origin() == infinite_point; t = t.twin().next());

				n.next(t);
			}

			delete_faces.push_back(e.face().index());
			delete_edges.push_back(e.index());
			e = e.twin();
			delete_edges.push_back(e.index());
			e = e.next();

			if (e == infinite_point.edge()) break;
		}

		// remove edges incident to infinite vertex #0
		std::sort(delete_edges.begin(), delete_edges.end(), std::greater<>{});
		for (auto e : delete_edges)
			dcel.remove_edge(dcel.edge(e));

		// remove vertex #0
		dcel.remove_point(infinite_point);

		std::sort(delete_faces.begin(), delete_faces.end(), std::greater<>{});
		for (auto f : delete_faces)
			dcel.remove_face(dcel.face(f));
	}

	// iteratively removes degenerate edges
	template <typename Point, typename Face, typename Edge, typename Index>
	void remove_degenerate_edges(dcel<Point, Face, Edge, Index> & dcel)
	{
		std::vector<std::size_t> delete_points;
		std::vector<std::size_t> delete_edges;

		for (std::size_t i = 0; i < dcel.points.size(); ++i)
		{
			auto p = dcel.point(i);
			while (p.edge().twin().next() == p.edge())
			{
				delete_points.push_back(p.index());
				delete_edges.push_back(p.edge().index());
				delete_edges.push_back(p.edge().twin().index());

				p.edge().next().origin().edge(p.edge().next());
				p.edge().face().edge(p.edge().next());

				auto n = p.edge();
				while (true)
				{
					auto m = n.next().twin();
					if (m == p.edge()) break;
					n = m;
				}

				n.next(p.edge().next());

				p = p.edge().next().origin();
			}
		}

		std::sort(delete_points.begin(), delete_points.end(), std::greater<>{});
		for (auto p : delete_points)
			dcel.remove_point(dcel.point(p));

		std::sort(delete_edges.begin(), delete_edges.end(), std::greater<>{});
		for (auto e : delete_edges)
			dcel.remove_edge(dcel.edge(e));
	}

	// Outputs 4 points such that the finite faces of the voronoi diagram of (anything
	// inside the convex hull of input + the output points) covers the convex hull of input
	template <typename InputIterator, typename OutputIterator>
	OutputIterator bounded_voronoi_extra_points(InputIterator begin, InputIterator end, OutputIterator out)
	{
		using point_type = std::decay_t<decltype(*begin)>;

		static_assert(point_type::dimension == 2);

		auto box = bbox(begin, end);

		point_type const center {(box[0].min + box[0].max) / 2, (box[1].min + box[1].max) / 2};

		auto const R = std::sqrt(box[0].length() * box[0].length() + box[1].length() * box[1].length());

		*out++ = point_type{center[0] + R, center[1]};
		*out++ = point_type{center[0], center[1] + R};
		*out++ = point_type{center[0] - R, center[1]};
		*out++ = point_type{center[0], center[1] - R};

		return out;
	}

	// Compute a clipped voronoi tessellation of (begin, end - 4)
	// Last 4 points will be used for temporary data
	template <typename InputIterator>
	auto bounded_voronoi(InputIterator begin, InputIterator end)
	{
		bounded_voronoi_extra_points(begin, end - 4, end - 4);
		auto dcel = cg::voronoi(begin, end);
		cg::remove_infinite_point(dcel);
		return dcel;
	}

}