smallest-circle Algorithm
The smallest-circle algorithm, also known as the minimum enclosing circle algorithm, is a computational geometry technique used to find the smallest circle that completely encloses a given set of points in a two-dimensional plane. This minimal bounding circle is essential in various applications, such as clustering analysis, data compression, visualization, and collision detection. The algorithm's primary goal is to find the center point and radius of the smallest circle that can enclose all the given points without any assumptions about the data distribution.
Several smallest-circle algorithms have been proposed, each differing in their approach and computational complexity. One popular method is Welzl's algorithm, which is a randomized, recursive, and efficient algorithm that works based on the idea of linear programming. Other approaches include the Skyum algorithm, which is a simple and elegant incremental algorithm, and the Megiddo algorithm, which is a linear-time deterministic algorithm. In general, these algorithms use techniques such as iterative improvements, divide and conquer, or convex hulls to determine the smallest possible enclosing circle efficiently. The performance of these algorithms varies depending on factors such as the number of input points, the distribution of the points, and the required precision for the final circle.
#include <iostream>
#include <vector>
#include <math.h>
using namespace std;
struct Point
{
double x, y;
Point(double a = 0.0, double b = 0.0)
{
x = a;
y = b;
}
};
double LenghtLine(Point A, Point B)
{
return sqrt(abs((B.x - A.x) * (B.x - A.x)) + abs((B.y - A.y) * (B.y - A.y)));
}
double TriangleArea(Point A, Point B, Point C)
{
double a = LenghtLine(A, B);
double b = LenghtLine(B, C);
double c = LenghtLine(C, A);
double p = (a + b + c) / 2;
return sqrt(p * (p - a) * (p - b) * (p - c));
}
bool PointInCircle(vector<Point> &P, Point Center, double R)
{
for (size_t i = 0; i < P.size(); i++)
{
if (LenghtLine(P[i], Center) > R)
return false;
}
return true;
}
double circle(vector<Point> P)
{
double minR = INT8_MAX;
double R;
Point C;
Point minC;
for (size_t i = 0; i < P.size() - 2; i++)
for (size_t j = i + 1; j < P.size(); j++)
for (size_t k = j + 1; k < P.size(); k++)
{
C.x = -0.5 * ((P[i].y * (P[j].x * P[j].x + P[j].y * P[j].y - P[k].x * P[k].x - P[k].y * P[k].y) + P[j].y * (P[k].x * P[k].x + P[k].y * P[k].y - P[i].x * P[i].x - P[i].y * P[i].y) + P[k].y * (P[i].x * P[i].x + P[i].y * P[i].y - P[j].x * P[j].x - P[j].y * P[j].y)) / (P[i].x * (P[j].y - P[k].y) + P[j].x * (P[k].y - P[i].y) + P[k].x * (P[i].y - P[j].y)));
C.y = 0.5 * ((P[i].x * (P[j].x * P[j].x + P[j].y * P[j].y - P[k].x * P[k].x - P[k].y * P[k].y) + P[j].x * (P[k].x * P[k].x + P[k].y * P[k].y - P[i].x * P[i].x - P[i].y * P[i].y) + P[k].x * (P[i].x * P[i].x + P[i].y * P[i].y - P[j].x * P[j].x - P[j].y * P[j].y)) / (P[i].x * (P[j].y - P[k].y) + P[j].x * (P[k].y - P[i].y) + P[k].x * (P[i].y - P[j].y)));
R = (LenghtLine(P[i], P[j]) * LenghtLine(P[j], P[k]) * LenghtLine(P[k], P[i])) / (4 * TriangleArea(P[i], P[j], P[k]));
if (!PointInCircle(P, C, R))
{
continue;
}
if (R <= minR)
{
minR = R;
minC = C;
}
}
for (size_t i = 0; i < P.size() - 1; i++)
for (size_t j = i + 1; j < P.size(); j++)
{
C.x = (P[i].x + P[j].x) / 2;
C.y = (P[i].y + P[j].y) / 2;
R = LenghtLine(C, P[i]);
if (!PointInCircle(P, C, R))
{
continue;
}
if (R <= minR)
{
minR = R;
minC = C;
}
}
cout << minC.x << " " << minC.y << endl;
return minR;
}
void test()
{
vector<Point> Pv(5);
Pv.push_back(Point(0, 0));
Pv.push_back(Point(1, 3));
Pv.push_back(Point(4, 1));
Pv.push_back(Point(5, 4));
Pv.push_back(Point(3, -2));
cout << circle(Pv) << endl;
}
void test2()
{
vector<Point> Pv(4);
Pv.push_back(Point(0, 0));
Pv.push_back(Point(0, 2));
Pv.push_back(Point(2, 2));
Pv.push_back(Point(2, 0));
cout << circle(Pv) << endl;
}
void test3()
{
vector<Point> Pv(3);
Pv.push_back(Point(0.5, 1));
Pv.push_back(Point(3.5, 3));
Pv.push_back(Point(2.5, 0));
cout << circle(Pv) << endl;
}
int main()
{
test();
cout << endl;
test2();
cout << endl;
test3();
return 0;
}
