mirror of
https://hub.njuu.cf/TheAlgorithms/C-Plus-Plus.git
synced 2023-10-11 13:05:55 +08:00
updated documentation - NOTE TESTS FAIL
@yanglbme this algorithm does not give right answers
This commit is contained in:
Krishna Vedala
parent
00c8c0bfdb
commit
be997d3da3
@ -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 <cmath>
|
||||||
#include <iostream>
|
#include <iostream>
|
||||||
#include <vector>
|
#include <vector>
|
||||||
#include <math.h>
|
|
||||||
|
|
||||||
using namespace std;
|
/** Define a point */
|
||||||
|
struct Point {
|
||||||
|
double x, /**< abscissa */
|
||||||
|
y; /**< ordinate */
|
||||||
|
|
||||||
struct Point
|
/** construct a point
|
||||||
{
|
* \param [in] a absicca (default = 0.0)
|
||||||
double x, y;
|
* \param [in] b ordinate (default = 0.0)
|
||||||
Point(double a = 0.0, double b = 0.0)
|
*/
|
||||||
{
|
Point(double a = 0.0, double b = 0.0) {
|
||||||
x = a;
|
x = a;
|
||||||
y = b;
|
y = b;
|
||||||
}
|
}
|
||||||
};
|
};
|
||||||
|
|
||||||
double LenghtLine(Point A, Point B)
|
/** 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:
|
||||||
return sqrt(abs((B.x - A.x) * (B.x - A.x)) + abs((B.y - A.y) * (B.y - A.y)));
|
* \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 a = LenghtLine(A, B);
|
||||||
double b = LenghtLine(B, C);
|
double b = LenghtLine(B, C);
|
||||||
double c = LenghtLine(C, A);
|
double c = LenghtLine(C, A);
|
||||||
double p = (a + b + c) / 2;
|
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<Point> &P, Point Center, double R)
|
/**
|
||||||
{
|
* Check if a set of points lie within given circle. This is true if the
|
||||||
for (size_t i = 0; i < P.size(); i++)
|
* 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<Point> &P, const Point &Center, double R) {
|
||||||
|
for (size_t i = 0; i < P.size(); i++) {
|
||||||
if (LenghtLine(P[i], Center) > R)
|
if (LenghtLine(P[i], Center) > R)
|
||||||
return false;
|
return false;
|
||||||
}
|
}
|
||||||
return true;
|
return true;
|
||||||
}
|
}
|
||||||
|
|
||||||
double circle(vector<Point> 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<Point> &P) {
|
||||||
double minR = INT8_MAX;
|
double minR = INT8_MAX;
|
||||||
double R;
|
double R;
|
||||||
Point C;
|
Point C;
|
||||||
Point minC;
|
Point minC;
|
||||||
|
|
||||||
|
/* 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 (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 (size_t j = i + 1; j < P.size(); j++)
|
||||||
for (size_t k = j + 1; k < P.size(); k++)
|
// for every subsequent point in the list
|
||||||
{
|
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)));
|
// here, we now have picked three points from the given set
|
||||||
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)));
|
of
|
||||||
R = (LenghtLine(P[i], P[j]) * LenghtLine(P[j], P[k]) * LenghtLine(P[k], P[i])) / (4 * TriangleArea(P[i], P[j], P[k]));
|
// points that we can use
|
||||||
if (!PointInCircle(P, C, R))
|
// 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 -
|
||||||
continue;
|
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)
|
if (R <= minR) {
|
||||||
{
|
|
||||||
minR = R;
|
minR = R;
|
||||||
minC = C;
|
minC = C;
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
*/
|
||||||
|
|
||||||
|
// for each point in the list
|
||||||
for (size_t i = 0; i < P.size() - 1; i++)
|
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.x = (P[i].x + P[j].x) / 2;
|
||||||
C.y = (P[i].y + P[j].y) / 2;
|
C.y = (P[i].y + P[j].y) / 2;
|
||||||
R = LenghtLine(C, P[i]);
|
R = LenghtLine(C, P[i]);
|
||||||
if (!PointInCircle(P, C, R))
|
if (!PointInCircle(P, C, R)) {
|
||||||
{
|
|
||||||
continue;
|
continue;
|
||||||
}
|
}
|
||||||
if (R <= minR)
|
if (R <= minR) {
|
||||||
{
|
|
||||||
minR = R;
|
minR = R;
|
||||||
minC = C;
|
minC = C;
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
cout << minC.x << " " << minC.y << endl;
|
std::cout << minC.x << " " << minC.y << std::endl;
|
||||||
return minR;
|
return minR;
|
||||||
}
|
}
|
||||||
|
|
||||||
void test()
|
/** Test case: result should be:
|
||||||
{
|
* \n Circle with
|
||||||
vector<Point> Pv(5);
|
* \n radius 3.318493136080724
|
||||||
|
* \n centre at (3.0454545454545454 1.3181818181818181)
|
||||||
|
*/
|
||||||
|
void test() {
|
||||||
|
std::vector<Point> Pv(5);
|
||||||
Pv.push_back(Point(0, 0));
|
Pv.push_back(Point(0, 0));
|
||||||
Pv.push_back(Point(1, 3));
|
Pv.push_back(Point(1, 3));
|
||||||
Pv.push_back(Point(4, 1));
|
Pv.push_back(Point(4, 1));
|
||||||
Pv.push_back(Point(5, 4));
|
Pv.push_back(Point(5, 4));
|
||||||
Pv.push_back(Point(3, -2));
|
Pv.push_back(Point(3, -2));
|
||||||
cout << circle(Pv) << endl;
|
std::cout << circle(Pv) << std::endl;
|
||||||
}
|
}
|
||||||
|
|
||||||
void test2()
|
/** Test case: result should be:
|
||||||
{
|
* \n Circle with
|
||||||
vector<Point> Pv(4);
|
* \n radius 1.4142135623730951
|
||||||
|
* \n centre at (1.0 1.0)
|
||||||
|
*/
|
||||||
|
void test2() {
|
||||||
|
std::vector<Point> Pv(4);
|
||||||
Pv.push_back(Point(0, 0));
|
Pv.push_back(Point(0, 0));
|
||||||
Pv.push_back(Point(0, 2));
|
Pv.push_back(Point(0, 2));
|
||||||
Pv.push_back(Point(2, 2));
|
Pv.push_back(Point(2, 2));
|
||||||
Pv.push_back(Point(2, 0));
|
Pv.push_back(Point(2, 0));
|
||||||
cout << circle(Pv) << endl;
|
std::cout << circle(Pv) << std::endl;
|
||||||
}
|
}
|
||||||
|
|
||||||
void test3()
|
/** Test case: result should be:
|
||||||
{
|
* \n Circle with
|
||||||
vector<Point> Pv(3);
|
* \n radius 1.821078397711709
|
||||||
|
* \n centre at (2.142857142857143 1.7857142857142856)
|
||||||
|
*/
|
||||||
|
void test3() {
|
||||||
|
std::vector<Point> Pv(3);
|
||||||
Pv.push_back(Point(0.5, 1));
|
Pv.push_back(Point(0.5, 1));
|
||||||
Pv.push_back(Point(3.5, 3));
|
Pv.push_back(Point(3.5, 3));
|
||||||
Pv.push_back(Point(2.5, 0));
|
Pv.push_back(Point(2.5, 0));
|
||||||
cout << circle(Pv) << endl;
|
std::cout << circle(Pv) << std::endl;
|
||||||
}
|
}
|
||||||
int main()
|
|
||||||
{
|
/** Main program */
|
||||||
|
int main() {
|
||||||
test();
|
test();
|
||||||
cout << endl;
|
std::cout << std::endl;
|
||||||
test2();
|
test2();
|
||||||
cout << endl;
|
std::cout << std::endl;
|
||||||
test3();
|
test3();
|
||||||
return 0;
|
return 0;
|
||||||
}
|
}
|
||||||
|
|||||||
Loading…
Reference in New Issue
Block a user