From 7b1e723bef24191d4cac6772129c60f441ebe04b Mon Sep 17 00:00:00 2001
From: Bruno Freitas Tissei <bft15@inf.ufpr.br>
Date: Thu, 29 Aug 2019 13:56:45 -0300
Subject: [PATCH] Add circle
Signed-off-by: Bruno Freitas Tissei <bft15@inf.ufpr.br>
---
algorithms/geometry/circle.cpp | 47 +++++++++++++++++++++
algorithms/geometry/geometry_primitives.cpp | 5 +++
2 files changed, 52 insertions(+)
create mode 100644 algorithms/geometry/circle.cpp
diff --git a/algorithms/geometry/circle.cpp b/algorithms/geometry/circle.cpp
new file mode 100644
index 0000000..91f3555
--- /dev/null
+++ b/algorithms/geometry/circle.cpp
@@ -0,0 +1,47 @@
+/// Circle
+
+struct Circle {
+ Point<double> c;
+ double r;
+
+ Circle(Point<double> c, double r) : c(c), r(r) {}
+
+ // Circumcircle
+ Circle(Point<double> a, Point<double> b, Point<double> c) {
+ Point<double> u((b - a).y, -(b - a).x);
+ Point<double> v((c - a).y, -(c - a).x);
+ Point<double> n = (c - b)*0.5;
+
+ double t = u.cross(n) / v.cross(u);
+
+ this->c = (a + c)*0.5 + v*t;
+ this->r = dist(this->c, a);
+ }
+
+ bool contains(Point<double> p) {
+ return (dist(c, p) <= r + EPS);
+ }
+};
+
+// Minumum enclosing circle
+Circle min_enclosing(vector<Point<double>> p) {
+ random_shuffle(all(p));
+ Circle C(p[0], 0.0);
+
+ for (int i = 0; i < p.size(); ++i) {
+ if (C.contains(p[i])) continue;
+ C = Circle(p[i], 0.0);
+
+ for (int j = 0; j < i; ++j) {
+ if (C.contains(p[j])) continue;
+ C = Circle((p[j] + p[i])*0.5, 0.5*dist(p[j], p[i]));
+
+ for (int k = 0; k < j; ++k) {
+ if (C.contains(p[k])) continue;
+ C = Circle(p[j], p[i], p[k]);
+ }
+ }
+ }
+
+ return C;
+}
diff --git a/algorithms/geometry/geometry_primitives.cpp b/algorithms/geometry/geometry_primitives.cpp
index 590ee0a..041afc2 100644
--- a/algorithms/geometry/geometry_primitives.cpp
+++ b/algorithms/geometry/geometry_primitives.cpp
@@ -12,6 +12,7 @@ struct Point {
Point operator+(Point p) { return Point(x+p.x, y+p.y); }
Point operator-(Point p) { return Point(x-p.x, y-p.y); }
+ Point operator*(T s) { return Point(x*s, y*s); }
T dot(Point p) { return (x*p.x) + (y*p.y); }
T cross(Point p) { return (x*p.y) - (y*p.x); }
@@ -51,6 +52,10 @@ struct Point {
}
};
+double dist(Point<double> a, Point<double> b) {
+ return hypot(a.x - b.x, a.y - b.y);
+}
+
template <typename T>
struct Segment {
Point<T> a, b;
--
GitLab