diff --git a/smallest-circle.cpp b/smallest-circle.cpp new file mode 100644 index 000000000..f48af5e63 --- /dev/null +++ b/smallest-circle.cpp @@ -0,0 +1,121 @@ +#include +#include +#include + +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 &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 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 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 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 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; +}