From be997d3da34367a6bc5f52fe719f6b72dd8a6cd4 Mon Sep 17 00:00:00 2001 From: Krishna Vedala Date: Thu, 28 May 2020 13:39:58 -0400 Subject: [PATCH] updated documentation - NOTE TESTS FAIL @yanglbme this algorithm does not give right answers --- others/smallest-circle.cpp | 190 ++++++++++++++++++++++++++----------- 1 file changed, 132 insertions(+), 58 deletions(-) diff --git a/others/smallest-circle.cpp b/others/smallest-circle.cpp index 48437cd86..7eb913cd1 100644 --- a/others/smallest-circle.cpp +++ b/others/smallest-circle.cpp @@ -1,121 +1,195 @@ +/** + * @file + * @brief Get centre and radius of the + * [smallest circle](https://en.wikipedia.org/wiki/Smallest-circle_problem) + * that circumscribes given set of points. + * + * @see [other + * implementation](https://www.nayuki.io/page/smallest-enclosing-circle) + */ +#include #include #include -#include -using namespace std; +/** Define a point */ +struct Point { + double x, /**< abscissa */ + y; /**< ordinate */ -struct Point -{ - double x, y; - Point(double a = 0.0, double b = 0.0) - { + /** construct a point + * \param [in] a absicca (default = 0.0) + * \param [in] b ordinate (default = 0.0) + */ + 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))); +/** Compute the Euclidian distance between two points \f$A\equiv(x_1,y_1)\f$ and + * \f$B\equiv(x_2,y_2)\f$ using the formula: + * \f[d=\sqrt{\left(x_1-x_2\right)^2+\left(y_1-y_2\right)^2}\f] + * + * \param [in] A point A + * \param [in] B point B + * \return ditance + */ +double LenghtLine(const Point &A, const Point &B) { + double dx = B.x - A.x; + double dy = B.y - A.y; + return std::sqrt((dx * dx) + (dy * dy)); } -double TriangleArea(Point A, Point B, Point C) -{ +/** + * Compute the area of triangle formed by three points using [Heron's + * formula](https://en.wikipedia.org/wiki/Heron%27s_formula). + * If the lengths of the sides of the triangle are \f$a,\,b,\,c\f$ and + * \f$s=\displaystyle\frac{a+b+c}{2}\f$ is the semi-perimeter then the area is + * given by \f[A=\sqrt{s(s-a)(s-b)(s-c)}\f] + * \param [in] A vertex A + * \param [in] B vertex B + * \param [in] C vertex C + * \returns area of triangle + */ +double TriangleArea(const Point &A, const Point &B, const 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)); + return std::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++) - { +/** + * Check if a set of points lie within given circle. This is true if the + * distance of all the points from the centre of the circle is less than the + * radius of the circle + * \param [in] P set of points to check + * \param [in] Center coordinates to centre of the circle + * \param [in] R radius of the circle + * \returns True if P lies on or within the circle + * \returns False if P lies outside the circle + */ +bool PointInCircle(const std::vector &P, const 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) -{ +/** + * Find the centre and radius of a circle enclosing a set of points.\n + * The function returns the radius of the circle and prints the coordinated of + * the centre of the circle. + * \param [in] P vector of points + * \returns radius of the circle + */ +double circle(const std::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; + + /* This code is invalid and does not give correct result for TEST 3 + // for each point in the list + for (size_t i = 0; i < P.size() - 2; i++) + // for every subsequent point in the list + for (size_t j = i + 1; j < P.size(); j++) + // for every subsequent point in the list + for (size_t k = j + 1; k < P.size(); k++) { + // here, we now have picked three points from the given set + of + // points that we can use + // viz., P[i], P[j] and P[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; + } } - if (R <= minR) - { - minR = R; - minC = C; - } - } + */ + + // for each point in the list for (size_t i = 0; i < P.size() - 1; i++) - for (size_t j = i + 1; j < P.size(); j++) - { + // for every subsequent point in the list + 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)) - { + if (!PointInCircle(P, C, R)) { continue; } - if (R <= minR) - { + if (R <= minR) { minR = R; minC = C; } } - cout << minC.x << " " << minC.y << endl; + std::cout << minC.x << " " << minC.y << std::endl; return minR; } -void test() -{ - vector Pv(5); +/** Test case: result should be: + * \n Circle with + * \n radius 3.318493136080724 + * \n centre at (3.0454545454545454 1.3181818181818181) + */ +void test() { + std::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; + std::cout << circle(Pv) << std::endl; } -void test2() -{ - vector Pv(4); +/** Test case: result should be: + * \n Circle with + * \n radius 1.4142135623730951 + * \n centre at (1.0 1.0) + */ +void test2() { + std::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; + std::cout << circle(Pv) << std::endl; } -void test3() -{ - vector Pv(3); +/** Test case: result should be: + * \n Circle with + * \n radius 1.821078397711709 + * \n centre at (2.142857142857143 1.7857142857142856) + */ +void test3() { + std::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; + std::cout << circle(Pv) << std::endl; } -int main() -{ + +/** Main program */ +int main() { test(); - cout << endl; + std::cout << std::endl; test2(); - cout << endl; + std::cout << std::endl; test3(); return 0; }