Полностью сделал задачу 6 на вычислительную геометрию
This commit is contained in:
AZEN-SGG
parent
a08f95aef8
commit
3efd67c703
15 changed files with 570 additions and 60 deletions
175
ComputationalGeometry/6Ex/PythonVersion/main.py
Normal file
175
ComputationalGeometry/6Ex/PythonVersion/main.py
Normal file
|
|
@ -0,0 +1,175 @@
|
|||
# Python3 program to find the minimum enclosing
|
||||
# circle for N integer points in a 2-D plane
|
||||
from math import sqrt
|
||||
from random import randint, shuffle
|
||||
|
||||
# Defining infinity
|
||||
INF = 1e18
|
||||
MAX_POINT_COORD = 100
|
||||
|
||||
|
||||
# Structure to represent a 2D point
|
||||
class Point:
|
||||
def __init__(self, X=0, Y=0) -> None:
|
||||
self.X = X
|
||||
self.Y = Y
|
||||
|
||||
|
||||
# Structure to represent a 2D circle
|
||||
class Circle:
|
||||
def __init__(self, c=Point(), r=0) -> None:
|
||||
self.C = c
|
||||
self.R = r
|
||||
|
||||
def print(self) -> None:
|
||||
print(f"(x - {self.C.X})^2 + (y - {self.C.Y})^2 = {self.R}^2")
|
||||
|
||||
|
||||
# Function to return the euclidean distance
|
||||
# between two points
|
||||
def dist(a, b):
|
||||
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
|
||||
def is_inside(c, p):
|
||||
return dist(c.C, p) <= c.R
|
||||
|
||||
|
||||
# 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
|
||||
def get_circle_center(bx, by,
|
||||
cx, cy):
|
||||
B = bx * bx + by * by
|
||||
C = cx * cx + cy * cy
|
||||
D = bx * cy - by * cx
|
||||
return Point((cy * B - by * C) / (2 * D),
|
||||
(bx * C - cx * B) / (2 * D))
|
||||
|
||||
|
||||
# Function to return the smallest circle
|
||||
# that intersects 2 points
|
||||
def circle_from1(A, B):
|
||||
# Set the center to be the midpoint of A and B
|
||||
C = Point((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)
|
||||
|
||||
|
||||
# Function to return a unique circle that
|
||||
# intersects three points
|
||||
def circle_from2(A, B, C):
|
||||
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))
|
||||
|
||||
|
||||
# Function to check whether a circle
|
||||
# encloses the given points
|
||||
def is_valid_circle(c, P):
|
||||
# Iterating through all the points
|
||||
# to check whether the points
|
||||
# lie inside the circle or not
|
||||
for p in P:
|
||||
if (not is_inside(c, p)):
|
||||
return False
|
||||
return True
|
||||
|
||||
|
||||
# Function to return the minimum enclosing
|
||||
# circle for N <= 3
|
||||
def min_circle_trivial(P):
|
||||
assert (len(P) <= 3)
|
||||
if not P:
|
||||
return Circle()
|
||||
|
||||
elif (len(P) == 1):
|
||||
return Circle(P[0], 0)
|
||||
|
||||
elif (len(P) == 2):
|
||||
return circle_from1(P[0], P[1])
|
||||
|
||||
# To check if MEC can be determined
|
||||
# by 2 points only
|
||||
for i in range(3):
|
||||
for j in range(i + 1, 3):
|
||||
|
||||
c = circle_from1(P[i], P[j])
|
||||
if (is_valid_circle(c, P)):
|
||||
return c
|
||||
|
||||
return circle_from2(P[0], P[1], P[2])
|
||||
|
||||
|
||||
# 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.
|
||||
def welzl_helper(P, R, n):
|
||||
# Base case when all points processed or |R| = 3
|
||||
if (n == 0 or len(R) == 3):
|
||||
return min_circle_trivial(R)
|
||||
|
||||
# Pick a random point randomly
|
||||
idx = randint(0, n - 1)
|
||||
p = P[idx]
|
||||
|
||||
# Put the picked point at the end of P
|
||||
# since it's more efficient than
|
||||
# deleting from the middle of the vector
|
||||
P[idx], P[n - 1] = P[n - 1], P[idx]
|
||||
|
||||
# Get the MEC circle d from the
|
||||
# set of points P - :p
|
||||
d = welzl_helper(P, R.copy(), n - 1)
|
||||
|
||||
# If d contains p, return d
|
||||
if (is_inside(d, p)):
|
||||
return d
|
||||
|
||||
# Otherwise, must be on the boundary of the MEC
|
||||
R.append(p)
|
||||
|
||||
# Return the MEC for P - :p and R U :p
|
||||
return welzl_helper(P, R.copy(), n - 1)
|
||||
|
||||
|
||||
def welzl(P):
|
||||
P_copy = P.copy()
|
||||
shuffle(P_copy)
|
||||
return welzl_helper(P_copy, [], len(P_copy))
|
||||
|
||||
|
||||
def generate(num: int) -> list[Point]:
|
||||
data: list[Point] = list()
|
||||
|
||||
for i in range(num):
|
||||
data.append(Point(randint(-MAX_POINT_COORD, MAX_POINT_COORD), randint(-MAX_POINT_COORD, MAX_POINT_COORD)))
|
||||
|
||||
return data
|
||||
|
||||
|
||||
def printPoints(data: list[Point]) -> None:
|
||||
for i in range(len(data)):
|
||||
print(f"{chr(i + 65)} = ({data[i].X}, {data[i].Y})")
|
||||
|
||||
|
||||
# Driver code
|
||||
if __name__ == '__main__':
|
||||
num = int(input("Enter number of points: "))
|
||||
data = generate(num)
|
||||
|
||||
printPoints(data)
|
||||
|
||||
mec = welzl(data)
|
||||
mec.print()
|
||||
|
|
@ -24,21 +24,23 @@ double powd(double number) {
|
|||
|
||||
circle primitive(point * N, int nlen) {
|
||||
if (nlen == 3) {
|
||||
point temp;
|
||||
circle crcl;
|
||||
|
||||
for (int i = 0; i < nlen; ++i) {
|
||||
temp = N[2];
|
||||
N[2] = N[i];
|
||||
N[i] = temp;
|
||||
|
||||
crcl = centermass(N[0], N[1]);
|
||||
if (belongs(crcl, N[2])) return crcl;
|
||||
|
||||
N[i] = N[2];
|
||||
N[2] = temp;
|
||||
}
|
||||
for (int i = 0; i < 3; ++i) {
|
||||
for (int j = i + 1; j < 3; ++j) {
|
||||
circle crcl = centermass(N[i], N[j]);
|
||||
bool all_inside = true;
|
||||
|
||||
for (int k = 0; k < 3; ++k) {
|
||||
if (!belongs(crcl, N[k])) {
|
||||
all_inside = false;
|
||||
break;
|
||||
}
|
||||
}
|
||||
if (all_inside) {
|
||||
return crcl;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
return byThreePoints(N);
|
||||
} else if (nlen == 2) {
|
||||
return centermass(N[0], N[1]);
|
||||
|
|
@ -54,28 +56,28 @@ circle centermass(point p1, point p2) {
|
|||
}
|
||||
|
||||
circle byThreePoints(point * warp) {
|
||||
point center;
|
||||
double radius, x, y;
|
||||
double ang_a;
|
||||
double ang_b;
|
||||
|
||||
if (fabs(warp[1].x - warp[0].x) < exp || fabs(warp[2].x - warp[1].x) < exp) return (circle){.center=(point){0, 0}, .radius=-1};
|
||||
|
||||
ang_a = straightAngle(warp[1], warp[0]);
|
||||
ang_b = straightAngle(warp[2], warp[1]);
|
||||
|
||||
if (fabs(ang_b) < exp || fabs(ang_b - ang_a) < exp) {
|
||||
return (circle){.center=(point){0, 0}, .radius=-1};
|
||||
}
|
||||
|
||||
x = (ang_a * ang_b * (warp[0].y - warp[2].y) + ang_b * (warp[0].x + warp[1].x) - ang_a * (warp[1].x + warp[2].x)) / (2 * (ang_b - ang_a));
|
||||
y = (-(1/ang_b) * (x - (warp[1].x + warp[2].x) / 2) + (warp[1].y + warp[2].y) / 2);
|
||||
center = (point){x, y};
|
||||
point center = getCenter(warp);
|
||||
double radius;
|
||||
|
||||
center.x += warp[0].x, center.y += warp[0].y;
|
||||
radius = distance(center, warp[0]);
|
||||
|
||||
return (circle){center, radius};
|
||||
}
|
||||
|
||||
point getCenter(point * warp) {
|
||||
double sqrt1, sqrt2, scalar;
|
||||
point vec1 = (point){warp[1].x - warp[0].x, warp[1].y - warp[0].y};
|
||||
point vec2 = (point){warp[2].x - warp[0].x, warp[2].y - warp[0].y};
|
||||
|
||||
sqrt1 = vec1.x * vec1.x + vec1.y * vec1.y;
|
||||
sqrt2 = vec2.x * vec2.x + vec2.y * vec2.y;
|
||||
scalar = (vec1.x * vec2.y - vec2.x * vec1.y) * 2;
|
||||
|
||||
return (point){.x=((vec2.y * sqrt1 - vec1.y * sqrt2) / scalar), \
|
||||
.y=((vec1.x * sqrt2 - vec2.x * sqrt1) / scalar)};
|
||||
}
|
||||
|
||||
double straightAngle(point p1, point p2) {
|
||||
return (p1.y - p2.y) / (p1.x - p2.x);
|
||||
}
|
||||
|
|
@ -85,7 +87,14 @@ double distance(point p1, point p2) {
|
|||
}
|
||||
|
||||
void printCircle(circle crcl) {
|
||||
printf("(x - %.4lf)^2 + (y - %.4lf)^2 = %.4lf^2\n", crcl.center.x, crcl.center.y, crcl.radius);
|
||||
// printf("Center of circle at point (%.2lf, %.2lf)\nRadius is %.2lf\n", crcl.center.x, crcl.center.y, crcl.radius);
|
||||
printf("(x - %.4lf)^2 + (y - %.4lf)^2 = %.4lf^2\n\n", crcl.center.x, crcl.center.y, crcl.radius);
|
||||
printf("Center of circle at point (%.2lf, %.2lf)\nRadius is %.2lf\n", crcl.center.x, crcl.center.y, crcl.radius);
|
||||
}
|
||||
|
||||
int isCover(circle crcl, points pts) {
|
||||
for (int i = 0; i < pts.length; ++i) {
|
||||
if (!belongs(crcl, pts.array[i])) return i;
|
||||
}
|
||||
|
||||
return pts.length;
|
||||
}
|
||||
|
|
@ -14,8 +14,10 @@ double powd(double number);
|
|||
circle primitive(point * N, int nlen);
|
||||
circle centermass(point p1, point p2);
|
||||
circle byThreePoints(point * warp);
|
||||
point getCenter(point * warp);
|
||||
double straightAngle(point p1, point p2);
|
||||
double distance(point p1, point p2);
|
||||
void printCircle(circle crcl);
|
||||
int isCover(circle crcl, points pts);
|
||||
|
||||
#endif
|
||||
|
|
|
|||
|
|
@ -2,49 +2,46 @@
|
|||
|
||||
circle hope(point * ps, int length) {
|
||||
circle crcl;
|
||||
circle minimum = (circle){(point){0, 0}, 0};
|
||||
|
||||
if (length == 1) {
|
||||
return (circle){ps[0], 0};
|
||||
} else if (length < 1) {
|
||||
return minimum;
|
||||
}
|
||||
|
||||
for (int i = 0; i < length; ++i) {
|
||||
point temp = ps[i];
|
||||
ps[i] = ps[length - 1];
|
||||
point temp = ps[0];
|
||||
|
||||
for (int i = 1; i < length; ++i) ps[i - 1] = ps[i];
|
||||
ps[length - 1] = temp;
|
||||
|
||||
for (int j = 0; j < length - 1; ++j) {
|
||||
crcl = centermass(temp, ps[i]);
|
||||
if (isSuit(crcl, ps, length - 1)) return crcl;
|
||||
crcl = centermass(temp, ps[j]);
|
||||
if (isCover(crcl, (points){ps, length}) == length) {
|
||||
if (minimum.radius < exp || (minimum.radius - crcl.radius > exp)) {
|
||||
minimum = crcl;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
ps[length - 1] = ps[i];
|
||||
ps[i] = temp;
|
||||
}
|
||||
|
||||
for (int i = 0; i < length; ++i) {
|
||||
point temp = ps[i];
|
||||
ps[i] = ps[length - 1];
|
||||
ps[length - 1] = temp;
|
||||
|
||||
for (int j = 0; j < length - 1; ++j) {
|
||||
temp = ps[j];
|
||||
ps[j] = ps[length - 2];
|
||||
ps[length - 2] = temp;
|
||||
|
||||
for (int k = 0; k < length - 2; ++k) {
|
||||
point warp[] = {ps[length - 1], ps[length - 2], ps[k]};
|
||||
point warp[] = {ps[i], ps[j], ps[k]};
|
||||
crcl = byThreePoints(warp);
|
||||
|
||||
if (fabs(crcl.radius + 1) > exp && isSuit(crcl, ps, length - 2)) return crcl;
|
||||
if ((fabs(crcl.radius + 1) > exp) && (isCover(crcl, (points){ps, length}) == length)) {
|
||||
if (minimum.radius < exp || (minimum.radius - crcl.radius > exp)) {
|
||||
minimum = crcl;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
temp = ps[j];
|
||||
ps[j] = ps[length - 2];
|
||||
ps[length - 2] = temp;
|
||||
}
|
||||
|
||||
temp = ps[i];
|
||||
ps[i] = ps[length - 1];
|
||||
ps[length - 1] = temp;
|
||||
}
|
||||
|
||||
return (circle){.center=(point){0, 0}, .radius=0};
|
||||
return minimum;
|
||||
}
|
||||
|
||||
bool isSuit(circle crcl, point * ps, int length) {
|
||||
|
|
|
|||
|
|
@ -1 +1 @@
|
|||
1 0 -1 0
|
||||
36.8400 -33.0600 -11.2900 95.9600 -56.8500 -77.7900 61.9600 59.0600 -9.8400 -14.3700 1.0100 85.6400
|
||||
|
|
@ -6,6 +6,7 @@ int main(void) {
|
|||
points pts;
|
||||
point N[3];
|
||||
circle crcl;
|
||||
int numLose;
|
||||
|
||||
pts = getPoints();
|
||||
if (pts.array == NULL) return -1;
|
||||
|
|
@ -16,9 +17,15 @@ int main(void) {
|
|||
crcl = MEC(pts.array, pts.length, N, 0);
|
||||
printCircle(crcl);
|
||||
|
||||
if ((numLose = isCover(crcl, pts)) == pts.length) printf("\nThe circle covers all the points!\n");
|
||||
else printf("\nThe circle misses a maximum of %d points", pts.length - numLose);
|
||||
|
||||
printf("\nReliable algorithm:\n");
|
||||
crcl = hope(pts.array, pts.length);
|
||||
printCircle(crcl);
|
||||
|
||||
if (isCover(crcl, pts) == pts.length) printf("\nThe circle covers all the points!\n");
|
||||
else printf("\nThe circle misses a maximum of %d points", pts.length - numLose);
|
||||
|
||||
free(pts.array);
|
||||
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue