Indoor/geo/Circle2.h

#ifndef CIRCLE2_H
#define CIRCLE2_H

#include <vector>

#include "Point2.h"
#include "Ray2.h"

#include "../Assertions.h"

//#include <KLib/misc/gnuplot/Gnuplot.h>
//#include <KLib/misc/gnuplot/GnuplotPlot.h>
//#include <KLib/misc/gnuplot/GnuplotPlotElementLines.h>
//#include <KLib/misc/gnuplot/GnuplotPlotElementPoints.h>

struct Circle2 {

	Point2 center;

	float radius;

public:

	/** empty ctor */
	Circle2() : center(), radius(0) {;}

	/** ctor */
	Circle2(const Point2 center, const float radius) : center(center), radius(radius) {;}


	/** does this circle contain the given point? */
	bool contains(const Point2 p) const {
		return center.getDistance(p) <= radius;
	}

	/** does this circle contain the given point? */
	bool containsAll(const std::vector<Point2>& pts) const {
		for (const Point2& p : pts) {
			if (!contains(p)) {return false;}
		}
		return true;
	}

	/** get a point on the circle for the given radians */
	Point2 getPointAt(const float rad) const {
		return center + Point2(std::cos(rad), std::sin(rad)) * radius;
	}


	/** does this circle contain the given circle? */
	bool contains(const Circle2 c) const {
		return (center.getDistance(c.center)+c.radius) <= radius;
	}

	/** does this circle intersect with the given ray? */
	bool intersects(const Ray2 ray) const {

		// https://math.stackexchange.com/questions/311921/get-location-of-vector-circle-intersection
		const float a = ray.dir.x*ray.dir.x + ray.dir.y*ray.dir.y;
		const float b = 2 * ray.dir.x * (ray.start.x-center.x) + 2 * ray.dir.y * (ray.start.y - center.y);
		const float c = (ray.start.x-center.x) * (ray.start.x-center.x) + (ray.start.y - center.y)*(ray.start.y - center.y) - radius*radius;
		const float discr = b*b - 4*a*c;

		return discr >= 0;

	}


	/** configure this sphere to contain the given point-set */
	void adjustToPointSet(const std::vector<Point2>& lst) {

		//adjustToPointSetAPX(lst);
		adjustToPointSetExact(lst);

		// validate
		for (const Point2& p : lst) {
			Assert::isTrue(this->contains(p), "calculated circle seems incorrect");
		}

	}

	/** combine two circles into a new one containing both */
	static Circle2 join(const Circle2& a, const Circle2& b) {

		// if one already contains the other, just return it as-is
		if (a.contains(b)) {return a;}
		if (b.contains(a)) {return b;}

		// create both maximum ends
		const Point2 out1 = a.center + (a.center-b.center).normalized() * a.radius;
		const Point2 out2 = b.center + (b.center-a.center).normalized() * b.radius;

		// center is within both ends, so is the radius
		Point2 newCen = (out1 + out2) / 2;
		float newRad = out1.getDistance(out2) / 2 * 1.0001;	// slightly larger to prevent rounding issues
		Circle2 res(newCen, newRad);

		// check
		Assert::isTrue(res.contains(a), "sphere joining error. base-spheres are not contained.");
		Assert::isTrue(res.contains(b), "sphere joining error. base-spheres are not contained.");

		return res;

	}

	float getArea() const  {
		return M_PI * (radius*radius);
	}

private:

	/*
	void show(std::vector<Point2> pts, const Point2 P, const Point2 Q, const Point2 R = Point2(NAN, NAN)) {

		static K::Gnuplot gp;

		K::GnuplotPlot plot;
		K::GnuplotPlotElementPoints gPoints; plot.add(&gPoints); gPoints.setColor(K::GnuplotColor::fromHexStr("#0000ff")); gPoints.setPointSize(1);
		K::GnuplotPlotElementLines gLines; plot.add(&gLines);

		for (const Point2 p : pts) {
			K::GnuplotPoint2 p2(p.x, p.y);
			gPoints.add(p2);
		}

		K::GnuplotPoint2 gP(P.x, P.y);
		K::GnuplotPoint2 gQ(Q.x, Q.y);
		gLines.addSegment(gP, gQ);

		K::GnuplotPoint2 gR(R.x, R.y);
		//gLines.addSegment(gP, gR);
		gLines.addSegment(gQ, gR);

		for (float f = 0; f < M_PI*2; f += 0.1) {
			Point2 p = center + Point2(std::cos(f), std::sin(f)) * radius;
			K::GnuplotPoint2 gp (p.x, p.y);
			gLines.add(gp);
		}

		gp << "set size ratio -1\n";
		gp.draw(plot);
		gp.flush();

		int i = 0; (void) i;

	}
	*/


	// Graphic Gems 2 - Jon Rokne
	void adjustToPointSetExact(const std::vector<Point2>& _lst) {

		if (_lst.size() == 2) {
			this->center = (_lst[0] + _lst[1]) / 2;
			this->radius = _lst[0].getDistance(_lst[1]) / 2 * 1.0001f;
			return;
		}

		std::vector<Point2> lst = _lst;

		// find point P having min(p.y)
		// NOTE: seems like the original work uses another coordinate system. so we search for max(p.y) instead!
		auto compMinY = [] (const Point2 p1, const Point2 p2) {return p1.y < p2.y;};
		auto it1 = std::max_element(lst.begin(), lst.end(), compMinY);
		Point2 P = *it1;
		lst.erase(it1);

		// find a point Q such that the angle of the line segment
		// PQ with the x axis is minimal
		auto compMinAngleX = [&] (const Point2 p1, const Point2 p2) {

			const Point2 PQ1 = p1 - P;
			const Point2 PQ2 = p2 - P;

			const float angle1 = dot(PQ1.normalized(), Point2(0,1));
			const float angle2 = dot(PQ2.normalized(), Point2(0,1));

			return std::acos(angle1) < std::acos(angle2);

		};
		auto it2 = std::min_element(lst.begin(), lst.end(), compMinAngleX);
		Point2 Q = *it2;
		lst.erase(it2);

		// get the angle abc which is the angle at "b"
		auto angle = [] (const Point2 a, const Point2 b, const Point2 c) {
			const Point2 ba = a - b;
			const Point2 bc = c - b;
			return std::acos(dot(ba.normalized(), bc.normalized()));
		};

		// TODO: how many loops?
		for (size_t xx = 0; xx < lst.size()*10; ++xx) {

			auto compMinAnglePRQ = [P,Q,angle]  (const Point2 p1, const Point2 p2) {
				return std::abs(angle(P,p1,Q)) < std::abs(angle(P,p2,Q));
			};

			auto it3 = std::min_element(lst.begin(), lst.end(), compMinAnglePRQ);

			Point2 R = *it3;
			lst.erase(it3);

			const float anglePRQ = angle(P,R,Q);
			const float angleRPQ = angle(R,P,Q);
			const float anglePQR = angle(P,Q,R);

			//check for case 1 (angle PRQ is obtuse), the circle is determined
			//by two points, P and Q. radius = |(P-Q)/2|, center = (P+Q)/2
			if (anglePRQ > M_PI/2) {

				this->radius = P.getDistance(Q) / 2 * 1.001f;
				this->center = (P+Q)/2;
				//if (!containsAll(_lst)) {show(_lst, P, Q, R);}
				return;

			}

			if (angleRPQ > M_PI/2) {
				lst.push_back(P);
				P = R;
				continue;
			}

			if (anglePQR > M_PI/2) {
				lst.push_back(Q);
				Q = R;
				continue;
			}

			const Point2 mPQ = (P+Q)/2;
			const Point2 mQR = (Q+R)/2;

			const float numer = -(-mPQ.y * R.y + mPQ.y * Q.y + mQR.y * R.y - mQR.y * Q.y - mPQ.x * R.x + mPQ.x * Q.x + mQR.x * R.x - mQR.x * Q.x);
			const float denom =  (-Q.x * R.y + P.x * R.y - P.x * Q.y + Q.y * R.x - P.y * R.x + P.y * Q.x);
			const float t = numer / denom;

			const float cx  = -t * (Q.y - P.y) + mPQ.x;
			const float cy  =  t * (Q.x - P.x) + mPQ.y;
			this->center = Point2(cx, cy);
			this->radius = this->center.getDistance(P) * 1.001f;

			//if (!containsAll(_lst)) {show(_lst, P, Q, R);}

			return;

		}

		throw Exception("should not happen");

	}

	/*
	void adjustToPointSetRefine(const std::vector<Point2>& lst) {

		float bestArea = 99999999;

		for (size_t i = 0; i < lst.size(); ++i) {
			for (size_t j = 0; j < lst.size(); ++j) {
				if (i == j) {continue;}

				const Point2 center = (lst[i] + lst[j]) / 3;
				const float radius = lst[i].getDistance(lst[j]);
				const Circle2 test(center, radius);

				if (test.containsAll(lst)) {
					if (test.getArea() < bestArea) {
						bestArea = test.getArea();
						this->radius = test.radius;
						this->center = test.center;
					}
				}
			}
		}

	}
*/
	void adjustToPointSetAPX(const std::vector<Point2>& lst) {

		// NOT OPTIMAL but fast

		// calculate the point set's center
		Point2 sum(0,0);
		for (const Point2 p : lst) {sum += p;}
		const Point2 center = sum / lst.size();

		// calculate the sphere's radius (maximum distance from center
		float radius = 0;
		for (const Point2 p : lst) {
			const float dist = center.getDistance(p);
			if (dist > radius) {radius = dist;}
		}

		this->radius = radius;
		this->center = center;




	}

};

#endif // CIRCLE2_H