kibalt? 10/12 12:57

minCircle.cpp

@@ -1,122 +1,166 @@
-//
-// Created by yarintz33 on 12/9/21.
-//
-#include <math.h>
+
 #include "minCircle.h"
+// circle for N integer points in a 2-D plane
+#include <algorithm>
+#include <assert.h>
+#include <iostream>
+#include <math.h>
+#include <vector>
+using namespace std;
 
 // Defining infinity
-const float INF = 91000000;
+//const double INF = 1e18;
 
+// Structure to represent a 2D point
 // Function to return the euclidean distance
 // between two points
-float dist(const Point& a, const Point& b)
+float dist(const Point& a,const Point& b)
 {
-    return sqrt(pow(a.x - b.x, 2) + pow(a.y - b.y, 2));
+    return sqrt(pow(a.x - b.x, 2)
+                + pow(a.y - b.y, 2));
 }
 
 // Function to check whether a point lies inside
 // or on the boundaries of the circle
-bool is_inside(Circle& c, Point& p)
+bool is_inside(Circle& c, Point*& p)
 {
-    return dist(c.center, p) <= c.radius;
+    return dist(c.center, *p) <= c.radius;
 }
 
-
-// The following two functions are the functions used
-// To find the equation of the circle when three
-// points are given.
+// The following two functions are used
+// To find the equation of the circle when
+// three points are given.
 
 // Helper method to get a circle defined by 3 points
-Point get_circle_center(float bx, float by,
-                        float cx, float cy)
+Point get_circle_center(double bx, double by,
+                        double cx, double cy)
 {
-    float B = bx * bx + by * by;
-    float C = cx * cx + cy * cy;
-    float D = bx * cy - by * cx;
-    return Point((cy * B - by * C) / (2 * D),
-                 (bx * C - cx * B) / (2 * D) );
+    double B = bx * bx + by * by;
+    double C = cx * cx + cy * cy;
+    double D = bx * cy - by * cx;
+    return Point( (cy * B - by * C) / (2 * D),
+                  (bx * C - cx * B) / (2 * D));
 }
 
-// Function to return a unique circle that intersects
-// three points
-Circle circle_from(Point& A, Point& B,
-                    Point& C)
+// Function to return a unique circle that
+// intersects three points
+Circle circle_from(const Point& A, const Point& B,
+                    const Point& C)
 {
     Point I = get_circle_center(B.x - A.x, B.y - A.y,
                                 C.x - A.x, C.y - A.y);
+
     I.x += A.x;
     I.y += A.y;
-    return Circle( I, dist(I, A)) ;
+    return Circle( I, dist(I, A));
 }
 
 // Function to return the smallest circle
 // that intersects 2 points
-Circle circle_from(Point& A, Point& B)
+Circle circle_from(const Point& A, const Point& B)
 {
     // Set the center to be the midpoint of A and B
-    Point C ( (A.x + B.x) / 2.0, (A.y + B.y) / 2.0);
+    Point C( (A.x + B.x) / 2.0, (A.y + B.y) / 2.0 );
 
     // Set the radius to be half the distance AB
-
-    return Circle(C, dist(A, B) / 2.0); //C, dist(A, B) / 2.0 };
+    return Circle( C, dist(A, B) / 2.0 );
 }
 
-// Function to check whether a circle encloses the given points
-bool is_valid_circle(Circle& c,  Point** points, size_t t)
+// Function to check whether a circle
+// encloses the given points
+bool is_valid_circle(Circle& c,
+                     vector<Point*>& P)
 {
 
-    // Iterating through all the points to check
-    // whether the points lie inside the circle or not
-    for (int i = 0; i < t; i++)
-        if (!is_inside(c, *points[i]))
+    // Iterating through all the points
+    // to check  whether the points
+    // lie inside the circle or not
+    for (Point*& p : P)
+        if (!is_inside(c, p))
             return false;
     return true;
 }
 
+// Function to return the minimum enclosing
+// circle for N <= 3
+Circle min_circle_trivial(vector<Point*>& P)
+{
+    assert(P.size() <= 3);
+    if (P.empty()) {
+        return Circle(Point(0,0), 0 );
+    }
+    else if (P.size() == 1) {
+        return Circle(*P[0] , 0);
+    }
+    else if (P.size() == 2) {
+        return circle_from(*P[0], *P[1]);
+    }
 
-Circle findMinCircle(Point** points,size_t size){
+    // To check if MEC can be determined
+    // by 2 points only
+    for (int i = 0; i < 3; i++) {
+        for (int j = i + 1; j < 3; j++) {
 
-    // To find the number of points
-    int n = size;
+            Circle c = circle_from(*P[i], *P[j]);
+            if (is_valid_circle(c, P))
+                return c;
+        }
+    }
+    return circle_from(*P[0], *P[1], *P[2]);
+}
 
-    if (n == 0)
-        return Circle( Point(0,0), 0 );
-    if (n == 1)
-        return Circle( *points[0], 0 );
+// Returns the MEC using Welzl's algorithm
+// Takes a set of input points P and a set R
+// points on the circle boundary.
+// n represents the number of points in P
+// that are not yet processed.
+Circle welzl_helper(vector<Point*>& P,
+                     vector<Point*> R, int n)
+{
+    // Base case when all points processed or |R| = 3
+    if (n == 0 || R.size() == 3) {
+        return min_circle_trivial(R);
+    }
 
-    // Set initial MEC to have infinity radius
-    float infimum = INF;
-    Circle mec(Point(0,0), infimum);
+    // Pick a random point randomly
+    int idx = rand() % n;
+    Point* p = P[idx];
 
-    // Go over all pair of points
-    for (int i = 0; i < n; i++) {
-        for (int j = i + 1; j < n; j++) {
+    // Put the picked point at the end of P
+    // since it's more efficient than
+    // deleting from the middle of the vector
+    swap(P[idx], P[n - 1]);
 
-            // Get the smallest circle that intersects P[i] and P[j]
-            Circle tmp = circle_from(*points[i], *points[j]); // points[i], points[j]
+    // Get the MEC circle d from the
+    // set of points P - {p}
+    Circle d = welzl_helper(P, R, n - 1);
 
-            // Update MEC if tmp encloses all points
-            // and has a smaller radius
-            if (tmp.radius < mec.radius && is_valid_circle(tmp, points, size))
-                mec = tmp;
-        }
+    // If d contains p, return d
+    if (is_inside(d, p)) {
+        return d;
     }
 
-    // Go over all triples of points
-    for (int i = 0; i < n; i++) {
-        for (int j = i + 1; j < n; j++) {
-            for (int k = j + 1; k < n; k++) {
+    // Otherwise, must be on the boundary of the MEC
+    R.push_back(p);
 
-                // Get the circle that intersects P[i], P[j], P[k]
-                Circle tmp = circle_from(*points[i], *points[j], *points[k]);
+    // Return the MEC for P - {p} and R U {p}
+    return welzl_helper(P, R, n - 1);
+}
 
-                // Update MEC if tmp encloses all points
-                // and has smaller radius
-                if (tmp.radius < mec.radius && is_valid_circle(tmp, points,size))
-                    mec = tmp;
-            }
-        }
+Circle welzl(vector<Point*>& P)
+{
+    vector<Point*> P_copy = P;
+    random_shuffle(P_copy.begin(), P_copy.end());
+    return welzl_helper(P_copy, {}, P_copy.size());
+}
+
+
+Circle findMinCircle(Point** points,size_t size){
+    vector<Point*> p ={};
+    for(int i = 0 ; i<5; i++){
+        p.push_back(points[i]);
     }
-    return mec;
+
+    return welzl(p);
 
 }