-
-
Notifications
You must be signed in to change notification settings - Fork 7.3k
/
smallest_circle.cpp
205 lines (193 loc) · 6.68 KB
/
smallest_circle.cpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
/**
* @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 <vector>
/** Define a point */
struct Point {
double x, /**< abscissa */
y; /**< ordinate */
/** construct a point
* \param [in] a absicca (default = 0.0)
* \param [in] b ordinate (default = 0.0)
*/
explicit Point(double a = 0.f, double b = 0.f) {
x = a;
y = 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:
* \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));
}
/**
* 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 std::sqrt(p * (p - a) * (p - b) * (p - c));
}
/**
* 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<Point> &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;
}
/**
* 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 = INFINITY;
double R;
Point C;
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 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;
}
}
// for each point in the list
for (size_t i = 0; i < P.size() - 1; i++)
// for every subsequent point in the list
for (size_t j = i + 1; j < P.size(); j++) {
// check for diameterically opposite points
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;
}
}
std::cout << minC.x << " " << minC.y << std::endl;
return minR;
}
/** Test case: result should be:
* \n Circle with
* \n radius 3.318493136080724
* \n centre at (3.0454545454545454, 1.3181818181818181)
*/
void test() {
std::vector<Point> Pv;
Pv.push_back(Point(0, 0));
Pv.push_back(Point(5, 4));
Pv.push_back(Point(1, 3));
Pv.push_back(Point(4, 1));
Pv.push_back(Point(3, -2));
std::cout << circle(Pv) << std::endl;
}
/** Test case: result should be:
* \n Circle with
* \n radius 1.4142135623730951
* \n centre at (1.0, 1.0)
*/
void test2() {
std::vector<Point> Pv;
Pv.push_back(Point(0, 0));
Pv.push_back(Point(0, 2));
Pv.push_back(Point(2, 2));
Pv.push_back(Point(2, 0));
std::cout << circle(Pv) << std::endl;
}
/** Test case: result should be:
* \n Circle with
* \n radius 1.821078397711709
* \n centre at (2.142857142857143, 1.7857142857142856)
* @todo This test fails
*/
void test3() {
std::vector<Point> Pv;
Pv.push_back(Point(0.5, 1));
Pv.push_back(Point(3.5, 3));
Pv.push_back(Point(2.5, 0));
Pv.push_back(Point(2, 1.5));
std::cout << circle(Pv) << std::endl;
}
/** Main program */
int main() {
test();
std::cout << std::endl;
test2();
std::cout << std::endl;
test3();
return 0;
}