334 lines
8.7 KiB
C++
334 lines
8.7 KiB
C++
/*
|
||
* © Copyright 2014 – Urheberrechtshinweis
|
||
* Alle Rechte vorbehalten / All Rights Reserved
|
||
*
|
||
* Programmcode ist urheberrechtlich geschuetzt.
|
||
* Das Urheberrecht liegt, soweit nicht ausdruecklich anders gekennzeichnet, bei Frank Ebner.
|
||
* Keine Verwendung ohne explizite Genehmigung.
|
||
* (vgl. § 106 ff UrhG / § 97 UrhG)
|
||
*/
|
||
|
||
#ifndef CIRCLE2_H
|
||
#define CIRCLE2_H
|
||
|
||
#include <vector>
|
||
|
||
#include "Point2.h"
|
||
#include "Ray2.h"
|
||
#include "Line2.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;
|
||
|
||
}
|
||
|
||
/** does this circle intersect with the given ray? */
|
||
bool intersects(const Line2 line) const {
|
||
|
||
// https://math.stackexchange.com/questions/311921/get-location-of-vector-circle-intersection
|
||
Point2 dir = line.p2 - line.p1;
|
||
Point2 start = line.p1;
|
||
const float a = dir.x*dir.x + dir.y*dir.y;
|
||
const float b = 2 * dir.x * (start.x-center.x) + 2 * dir.y * (start.y - center.y);
|
||
const float c = (start.x-center.x) * (start.x-center.x) + (start.y - center.y)*(start.y - center.y) - radius*radius;
|
||
const float discr = b*b - 4*a*c;
|
||
|
||
if (discr < 0) {return false;}
|
||
|
||
const float t = (2*c) / (-b + std::sqrt(discr));
|
||
return (t <= 1) && (t >= 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
|