Minimum coverage Circle Algorithm

The smallest circle coverage is a classic problem. The question is about N points on the plane. Finding the circle with the smallest radius can cover all points.

 

The algorithm is a bit difficult, so I want to explain my understanding.

If a minimum covering circle is required, the circle must be determined by at least three points. An algorithm is to take any three points as the circle, and then judge whether the point farthest from the center of the circle is in the circle. If it is in the circle, the circle is completed. If it is not in the circle, the circle is updated with the farthest point.

The algorithm introduced here is to first select any two points and use the line of the two points as the diameter for the circle. Then we can judge the remaining points to see if they are all in the circle (or on the circle). If they are all in, it indicates that the circle has been found. If none of them exist: assume that the first two vertices we use are P [1], p [2], and the first vertex we find is not in the circle (or on the circle) the point is P [I], so we use this point P [I] To find the smallest coverage circle covering P [1] to P [I-1.

So, how should we find the minimum coverage circle of the point P [I] from P [1] to P [I-1?

We first use P [1] and P [I] to do the circle, and then from 2 to the I-1 to determine whether a bit is not on this circle, if are in, it means that the circle that overwrites 1 to the I-1 has been found. If none of them exist: Suppose we find the first vertex not on this circle as P [J], so we use two known Points P [J] and P [I] To find the smallest covering circle covering 1 to J-1.

For the two known vertices P [J] and P [I] To find the smallest covering circle, as long as from 1 to the J-1, the K point is obtained through P [K], P [J], p [I] the circle of the three points, and then judge whether k + 1 to the J-1 are on the circle, if are in, it indicates to find the circle; if there is any other, use the new P [k] to update the circle.

As a result, this problem is transformed into several subproblems to solve.

Since the three points determine a circle, what we do in the process is from no definite point to a definite point, and then to two definite points, then there are three definite points to find the circle.

About the correctness of the proof and the complexity of the calculation here is not introduced, you can see the complete algorithm Introduction: http://wenku.baidu.com/view/162699d63186bceb19e8bbe6.html

Well. Details.

A. How can we find the circle through three points?

Calculate the cross product first.

If the cross product is 0, that is, the three points are in the same straight line, find the farthest pair of points, and use their line diameter as the circle;

If the cross product is not 0, that is, the three vertices are not collocated, it is the second problem. How can we find the outer circle of the triangle?

B. How do I find the outer circle of a triangle?

Suppose three points (x1, Y1), (X2, Y2), (X3, Y3 );

If the linear L1 equation of (x1, Y1) and (X2, Y2) is AX + by = C, its midpoint is (midx, midy) = (X1 + x2) /2, (Y1 + y2)/2), L1 vertical line equation is a1x + b1y = C1; then in its vertical line equation a1 =-B = x2-x1, b1 = A = y2-y1, C1 =-B * midx + A * midy = (X2 ^ 2-x1 ^ 2) + (Y2 ^ 2-y1 ^ 2)/2;

Likewise, we can know the vertical line equations of (x1, Y1) and (X3, Y3.

The intersection of the two vertical lines is the center of the center.

C. How to calculate the intersection of two straight lines?

Set the two straight lines to a1x + b1y = C1 and a2x + b2y = c2.

Set a variable det = A1 * B2-A2 * B1;

If det = 0, the two straight lines are parallel. if not equal to 0, the intersection is X = (B2 * C1-B1 * C2)/Det, y = (A1 * C2-A2 * C1)/det;

D. There is wood ..

1 # include <stdio. h> 2 # include <math. h> 3 struct tpoint 4 {5 double X, Y; 6}; 7 tpoint A [1005], D; 8 Double R; 9 10 double distance (tpoint P1, tpoint P2) // The distance between two points is 11 {12 Return (SQRT (p1.x-p2.x) * (p1.x-p2.x) + (p1.y-p2.y) * (p1.y-p2.y ))); 13} 14 double multiply (tpoint P1, tpoint P2, tpoint P0) 15 {16 return (p1.x-Snapshot X) * (p2.y-Snapshot y)-(p2.x-Snapshot X) * (p1.y-Policy); 17} 18 void minidiscwith2point (tpoint P, tpoin T q, int N) 19 {20 d. X = (P. X + q. x)/2.0; 21 d. y = (P. Y + q. y)/2.0; 22 r = distance (p, q)/2; 23 int K; 24 Double C1, C2, T1, T2, T3; 25 for (k = 1; k <= N; k ++) 26 {27 if (distance (D, a [k]) <= r) continue; 28 If (multiply (p, q, A [k])! = 0.0) 29 {30 C1 = (P. x * P. X + P. y * P. y-q.x * q. x-q.y * q. y)/2.0; 31 C2 = (P. x * P. X + P. y * P. y-A [K]. x * A [K]. x-A [K]. y * A [K]. y)/2.0; 32 33 d. X = (C1 * (p. y-A [K]. y)-C2 * (p. y-q.y)/(P. x-q.x) * (p. y-A [K]. y)-(P. x-A [K]. x) * (p. y-q.y); 34 D. y = (C1 * (p. x-A [K]. x)-C2 * (p. x-q.x)/(P. y-q.y) * (p. x-A [K]. x)-(P. y-A [K]. y) * (p. x-q.x); 35 r = distance (D, a [k]); 36} 37 else38 {39 T1 = distance (p, q); 40 t2 = distance (Q, A [k]); 41 T3 = distance (P, A [k]); 42 if (T1> = t2 & T1> = T3) 43 {d. X = (P. X + q. x)/2.0; D. y = (P. Y + q. y)/2.0; r = distance (p, q)/2.0;} 44 else if (T2> = T1 & T2> = T3) 45 {d. X = (a [K]. X + q. x)/2.0; D. y = (a [K]. Y + q. y)/2.0; r = distance (A [K], q)/2.0;} 46 else47 {d. X = (a [K]. X + P. x)/2.0; D. y = (a [K]. Y + P. y)/2.0; r = distance (A [K], p)/2.0;} 48} 49} 50} 51 52 void minidiscwithpoint (tpoint Pi, int N) 53 {54 d. X = (Pi. X + A [1]. x)/2.0; 55 d. y = (Pi. Y + A [1]. y)/2.0; 56 R = distance (PI, a [1])/2.0; 57 Int J; 58 for (j = 2; j <= N; j ++) 59 {60 if (distance (D, a [J]) <= r) continue; 61 else62 {63 minidiscwith2point (PI, a [J], J-1 ); 64} 65} 66} 67 int main () 68 {69 int I, n; 70 while (scanf ("% d", & N) 71 {72 for (I = 1; I <= N; I ++) 73 {74 scanf ("% lf", & A [I]. x, & A [I]. y); 75} 76 if (n = 1) 77 {printf ("%. 2lf %. 2lf 0.00 \ n ", a [1]. x, a [1]. y); continue;} 78 r = distance (A [1], a [2])/2.0; 79 d. X = (a [1]. X + A [2]. x)/2.0; 80 d. y = (a [1]. Y + A [2]. y)/2.0; 81 for (I = 3; I <= N; I ++) 82 {83 If (distance (D, a [I]) <= r) continue; 84 else85 minidiscwithpoint (A [I], I-1); 86} 87 printf ("%. 2lf %. 2lf %. 2lf \ n ", D. x, D. y, R); 88} 89 return 0; 90}

View code