寫在前面：剛學專業課的時候，記得有天突發奇想，心說高三數學的壓軸題能不能寫個程序跑出答案，這樣豈不是解放了數萬苦逼高三生的雙手？但是當時也僅僅是停留在想法上面，因為高中的解析幾何雖然步驟程序化，但是有時候需要靈巧的因式分解，感覺以目前的編程水平還是寫不出來，但是了解到數學有一個分支——計算幾何，專門利用計算機來進行幾何計算的一門科學，並且還與計算機圖形學、計算機視覺和圖像處理、機器人、計算機輔助設計和制造等高深學科有着聯系(摘自《計算幾何與應用》導言)，所以今天懷着激動的心情開始了這個專題的學習。
數學的一大功用就是計算，而計算的前提便是把自然界中的數據抽象成數學語言，數軸、坐標系都是為了更好的計算而引入的工具，同樣，向量——也是這個作用。
眾所周知，向量有兩個必不可少的因素——大小、方向。而在實際問題中，向量即作為有大小、有方向的具有實際物理意義的載體，也是進行代數運算的重要載體，而計算幾何中則充分體現了向量的這一特點。
我們就一個題目來初涉計算幾何這個領域：關於n邊形的面積求解的問題。
如果給出你n邊形的n個頂點的坐標，你能求出這個n邊形的面積么？
如果想要編程實現，顯然我們需要讓這個求解的過程變得自動化、規律化。我們可以考慮，將其依次分割成多個三角形，然后依次求出處三角形的面積。而如何分割呢？我們考慮，從某個點出發，然后依次連接多邊形的兩個頂點，然后就可以把圖形分割成多個三角形。
這個點可以在多邊形內部，可以是多邊形的頂點，可以某個頂點。
像下面這幾種分割方法。
我們來看圖3.3，顯然，現在我們的目標是依次求出所有的三角形然后加和即可，這在編程上也很好實現，但是現在問題是，給出的是點的坐標，我們如何更方便快捷的計算出三角形的面積呢？
這里的方法涉及向量叉積的概念了。
向量的叉積：
我們知道，向量是既有大小又有方向的，對於確定的兩個向量a，向量b，通過平移，是可以組成一個確定的三角形的(SAS，其詳細證明筆者認為《幾何原本》可能會給出)。而我們再根據三角形的面積計算公式s = 1/2*a*b*sinθ(θ是向量a、b的夾角，這里公式不予證明)，以及高中學到的向量的點積公式，我們可以做出如下推導。
那么我們可以得出結論，三角形的面積實際可以通過任意兩個邊所形成的向量a、b通過叉積運算得來。
然而我們看到，如果分割過程沒有將n邊形分成標准的三角形怎么辦？如圖3.4(b)，顯然上面的方法是存在不足的。
從叉積的定義我們可以看出，叉積運算也是有符號的。通過簡單的證明，我們可以找到叉積運算的符號規律。(這個探討好像偏離了問題，但是是對后面探索的一個鋪墊)
我們可以得出這樣一個規律——如果向量b的方向在向量a的逆時鍾旋轉π的范圍里，那么向量a×向量b的值是大於零的。
我們再回到剛剛的論證當中，顯然我們會考慮到，如果這多邊形是凹的，就可能出現不完全分割的情況，而在編程過程中我們給出的算法顯然要呈現出一般化的規律，因此無法保證能否找到一個合適的點進行分割，所以我們猜測，是否完全分割和不完全分割對於面積的計算是遵循一個公式呢？
讓我換一個簡單的模型來探究一下 這種不完全分割是怎么來的，這里有個誤區——認為凹多邊形一定會不完全分割，其實不然，圖3.4的a、b的對比就佐證了這個觀點。
這個模型是這樣來的，我們規定了AB的斜率不動，然后讓頂點o可以隨意移動。
在情況1中，我們看到，S△AOB被“擠”成了一條直線。
在情況2中，AO的斜率是小於AB的，此時可以完全分割，計算面積只需依次就出三角形的面積再加和即可。而S△AOB
= 向量AB × 向量OA恰好是正值，這也呼應了前面的過程。
在情況3中，情況則大不相同，這里S△AOB一部分是被掏空的，一部分被計算了兩次，因此在面積求和的時候應該減掉一個S△AOB，而此時向量AB × 向量OA
也恰好是符號，再次呼應，這也為給出不規則多邊形面積公式奠定了基礎。
這是一個最簡單的模型，推廣以后(精髓所在)，便可以得到上圖給出的公式。
基於以上的證明，我們來做一道簡單的題目來完成算法的代碼實現。(Problem source : hdu 2036)
話說部分學生心態極好，每天就知道游戲，這次考試如此簡單的題目，也是雲里霧里，而且，還竟然來這么幾句打油詩。 好呀，老師的責任就是幫你解決問題，既然想種田，那就分你一塊。 這塊田位於浙江省溫州市蒼南縣靈溪鎮林家鋪子村，多邊形形狀的一塊地，原本是linle 的，現在就准備送給你了。不過，任何事情都沒有那么簡單，你必須首先告訴我這塊地到底有多少面積，如果回答正確才能真正得到這塊地。 發愁了吧？就是要讓你知道，種地也是需要AC知識的！以后還是好好練吧...
很明顯的多邊形求面積的問題，在編程實現的時候，基於上文計算公式的得出，為了方便計算，我們取原點為分割凸多邊形的那個點，這樣我們要遍歷到多邊形的n個頂點，並通過叉積來計算小三角形的面積。遍歷到第n個頂點的時候，第n個頂點會和第1個頂點構成三角形，因此我們需要將記錄點集的數組的末尾再添加一個元素——即第1個頂點，這才計算幾何問題中遍歷點集進行運算是很常見的操作。
參考代碼如下。
#include<iostream> #include<stdio.h> using namespace std; struct point { double x,y; } P[101]; double cross(point a,point b) { return a.x*b.y-a.y*b.x; } double area(int n) { double tmp=0; for(int i=0;i<n;i++) { tmp+=cross(P[i],P[i+1]); } return tmp/2.00; } int main() { int i,n; while(scanf("%d",&n),n) { for(i=0;i<n;i++) scanf("%lf%lf",&P[i].x,&P[i].y); P[n].x=P[0].x;P[n].y=P[0].y; printf("%.1lf\n",area(n)); } return 0; }
有關三角形和圓的結合。(Problem source : 1374)
You are given the cartesian coordinates of three non-collinear points in the plane. Your job is to calculate the circumference of the unique circle that intersects all three points.
題目大意，讓你給出三角形的三個頂點，然后讓你計算出該三角形的外接圓的周長。
這里我們需要推導一下計算外接圓半徑的公式。根據正弦定理我們知道，c/sinc =2r，而s = 1/2a*b*sinc，我們可以得到r =a*b*c / 4s。
這里求三角形的面積s的時候，我們常用s = sqrt(p*(p-a)*(p-b)*(p-c))這個公式，其中p = 0.5(a + b + c)。
有了明確的公式，算法是實現上非常容易。
代碼如下。
#include<stdio.h> #include<math.h> #define PI 3.141592653589793 typedef struct point { double x , y; }Point; double distance(Point a ,Point b) { return sqrt((a.x - b.x)*(a.x - b.x) + (a.y - b.y)*(a.y - b.y)); } int main() { double a , b , c , p , r , s , c_of_circle; Point v[3]; while(scanf("%lf %lf %lf %lf %lf %lf",&v[0].x , &v[0].y , &v[1].x ,&v[1].y , &v[2].x , &v[2].y) != EOF) { a = distance(v[0] , v[1]); b = distance(v[1] , v[2]); c = distance(v[2] , v[0]); p = (a + b + c)/2.0; s = sqrt(p * ( p - a ) * ( p - b ) * ( p - c )); r = (a*b*c) / (4.0*s); c_of_circle = 2*PI*r; printf("%.2lf\n",c_of_circle); } }
我們再來看一道關於求解三角形外接圓的題目。(Problem soure : poj 1329)
Description
(x - h)^2 + (y - k)^2 = r^2 (1)and an equation of the form
x^2 + y^2 + cx + dy - e = 0 (2)
Input
Output
題目大意：基礎三個點的坐標，讓你輸出對應外接圓的頂點式和一般式。
數理分析：這里我們假設給出的三點坐標為(x1,y1)，(x2,y3)，(x3,y3)根據外心的定義，我們需要找到兩條中垂線的交點。
拿出兩個點p1(x1,y1)，p2(x2,y2)來說，中垂線過p1p2的中點，並且斜率與p1p2所成直線的斜率的乘積是-1，通過簡單的解析幾何的知識，我們可得到一條中垂線的參數方程。
(x-(x1+x2)/2)(x1-x2) = -(y1-y2)(y-(y1+y2)/2)。
同理再得到一條這樣的方程，兩個方程進行連理，既可以一個三角形外心通用的計算公式。
知道了圓心的坐標，我們只需計算圓心到某個頂點的距離，即可得到半徑r。
編程實現：題目給出的格式比較繁瑣，這里注意分情況討論一下處理一下符號即可。
參考代碼如下。
#include<stdio.h> #include<cmath> using namespace std; struct point { double x , y; }p[3],o; double dis(point a , point b) { return sqrt((a.x-b.x)*(a.x-b.x) + (a.y-b.y)*(a.y-b.y)); } void get_C(point p1,point p2,point p3)//求圓心 { double t1,t2,t3; t1=p1.x*p1.x+p1.y*p1.y-p2.x*p2.x-p2.y*p2.y; t2=p1.x*p1.x+p1.y*p1.y-p3.x*p3.x-p3.y*p3.y; t3=2*(p1.x-p2.x)*(p1.y-p3.y)-2*(p1.x-p3.x)*(p1.y-p2.y); o.x=((p1.y-p3.y)*t1-(p1.y-p2.y)*t2)/t3; o.y=-((p1.x-p3.x)*t1-(p1.x-p2.x)*t2)/t3; } int main() { while(scanf("%lf %lf %lf %lf %lf %lf",&p[0].x , &p[0].y , &p[1].x ,&p[1].y , &p[2].x , &p[2].y) != EOF) { get_C(p[0],p[1],p[2]); double r = dis(p[1],o); double t=r*r-o.x*o.x-o.y*o.y; if(o.x<0&&o.y<0) { printf("(x + %.3lf)^2 + (y + %.3lf)^2 = %.3lf^2\n",-o.x,-o.y,r); if(t>=0) printf("x^2 + y^2 + %.3lfx + %.3lfy - %.3lf = 0\n",-2*o.x,-2*o.y,t); else printf("x^2 + y^2 + %.3lfx + %.3lfy + %.3lf = 0\n",-2*o.x,-2*o.y,-t); } else if(o.x<0) { printf("(x + %.3lf)^2 + (y - %.3lf)^2 = %.3lf^2\n",-o.x,o.y,r); if(t>=0) printf("x^2 + y^2 + %.3lfx - %.3lfy - %.3lf = 0\n",-2*o.x,2*o.y,t); else printf("x^2 + y^2 + %.3lfx - %.3lfy + %.3lf = 0\n",-2*o.x,2*o.y,-t); } else if(o.y<0) { printf("(x - %.3lf)^2 + (y + %.3lf)^2 = %.3lf^2\n",o.x,-o.y,r); if(t>=0) printf("x^2 + y^2 - %.3lfx + %.3lfy - %.3lf = 0\n",2*o.x,-2*o.y,t); else printf("x^2 + y^2 - %.3lfx + %.3lfy + %.3lf = 0\n",2*o.x,-2*o.y,-t); } else { printf("(x - %.3lf)^2 + (y - %.3lf)^2 = %.3lf^2\n",o.x,o.y,r); if(t>=0) printf("x^2 + y^2 - %.3lfx - %.3lfy - %.3lf = 0\n",2*o.x,2*o.y,t); else printf("x^2 + y^2 - %.3lfx - %.3lfy + %.3lf = 0\n",2*o.x,2*o.y,-t); } printf("\n"); } }
我們再來看一些關於最小面積覆蓋的問題。(Problem source : hdu 1859)
讀完題后其實不難發現，這道題目的本質其實就是找到這些點的最大的橫坐標和縱坐標(右上角的點)，和最小的和坐標和縱坐標(左下角的點)，編程實現上也基本沒有什么難度。
這道題似乎有點簡答的弱智了，然是在簡短的代碼中其實還是體現了計算幾何問題在編程實現時代碼的規整之美，計算幾何其實就是基於計算機自動化、程序化的特點來解決一些繁瑣計算的幾何問題，而自動化的背后往往呈現出簡單的原理，掌握到背后簡單的原理，解決計算幾何問題是將會更加游刃有余。
由於其實規整的矩形覆蓋，所以問題非常好解決，但是如果是圓形呢？那其實就是計算幾何里面的一類問題——最小圓覆蓋。在后面的文章中我們將繼續詳細的介紹。
參考代碼如下:
#include<stdio.h> using namespace std; int main() { int maxx , maxy , minx , miny; int x , y; while(scanf("%d %d" , &x , &y) != EOF) { if(x == 0 && y == 0) break; maxx = x , minx = x , miny = y , maxy = y; while(1) { scanf("%d %d" , &x , &y); if(x == 0 && y == 0) { printf("%d %d %d %d\n" , minx , miny , maxx , maxy); break; } if(x > maxx) maxx = x; if(x < minx) minx = x; if(y > maxy) maxy = y; if(y < miny) miny = y; } } }
我們再來看一道平面解析幾何中最基礎的問題，關於點的對稱、直線交點的計算。(Problem source : hdu2857)
Now, our problem is that, if a branch of light goes into a large and infinite mirror, of course，it will reflect, and leave away the mirror in another direction. Giving you the position of mirror and the two points the light goes in before and after the reflection, calculate the reflection point of the light on the mirror. You can assume the mirror is a straight line and the given two points can’t be on the different sizes of the mirror.
(X1,Y1),(X2,Y2) mean the different points on the mirror, and (Xs,Ys) means the point the light travel in before the reflection, and (Xe,Ye) is the point the light go after the reflection.
The eight real number all are rounded to three digits after the decimal point, and the absolute values are no larger than 10000.0.
題目大意：這里先給出兩個點得到一條直線，當做鏡子。然后給定入射點和出射點，讓你計算出這條光線與鏡子的交點。
數理分析：很基礎的幾何題目，只需做出入射點關於鏡子的對稱點，對稱點和出射點構成的直線與鏡子的交點便是所求作的點。
編程實現：數理思維雖然不難，但是體現到編程上還是需要解決很多問題的。就比如：
①給定直線的兩點，我們如何得到一條直線？
②給定兩條直線，我們又如何求其交點？
③這里的對稱點怎么找？
這些問題雖然不難，但是還是需要比較巧的推理技巧。
關於給定兩點求直線，這里我們借助解析幾何中常見的斜截式y = kx + b，然后將其轉化成一般式ax + by+ c = 0，簡單的推導如下。
關於給定兩條直求其交點，我們采取一樣的方法。
而對於找對稱點，我們先找到互為對稱點的這兩個點組成的線段的中點，然后再進行計算。而這個中點恰好又是過入射點且垂直於鏡子的直線與鏡子的交點，這里的計算也很簡單。簡單的推導如下。
編程實現：有了以上推導，編程實現上就顯得容易很多。在編程實現的時候需要注意，我們求出交點，顯然有多個返回值，所以在將各種操作(對稱、求交點)寫成函數的時候，可以用到指針或者引用。
參考程序如下。(注：在推導過程中一般式按照ax+by+c=0推導，在編程的時候按照ax+by=c，這體現在代碼上可能有正負號的微小差異，本質是相同的)
#include<stdio.h> void Intersect(double a1,double b1,double c1,double a2,double b2,double c2,double &x,double &y) { x = (b2*c1 - b1*c2)/(a1*b2 - a2*b1); y = (a2*c1 - a1*c2)/(a2*b1 - a1*b2); //計算兩直線的交點 } void line(double x1,double y1,double x2,double y2,double &a,double &b,double &c) { a = y1 - y2; //計算兩點構成直線的參數a，b，c b = x2 - x1; c = x2*y1 - x1*y2; } int main() { int T; double x1,x2,y1,y2,x0,y0,x3,y3,a1,a2,b1,b2,c1,c2,x,y,a3,b3,c3; scanf("%d",&T); while(T--) { scanf("%lf%lf%lf%lf%lf%lf%lf%lf",&x1,&y1,&x2,&y2,&x0,&y0,&x3,&y3); a1 = x1 - x2; b1 = y1 - y2; c1 = x0*x1 - x0*x2 + y0*y1 - y0*y2; line(x1,y1,x2,y2,a2,b2,c2); Intersect(a1,b1,c1,a2,b2,c2,x,y); x += x - x0; y += y - y0; line(x,y,x3,y3,a3,b3,c3); Intersect(a3,b3,c3,a2,b2,c2,x,y); printf("%.3lf %.3lf\n",x,y); } }
我們再來看一道關於求解多邊形重心的問題。(Problem source:hdu1115)
題目大意：順次給你n個頂點的坐標，讓你找出n多邊形的重心。
數理分析：有了上面關於求解多邊形面積的思維鋪墊，這道題目將會輕松很多。
我們知道，重心是來源於物理學中的概念。這里我們給出求幾何體重心普適性的公式：
①如果質量集中在頂點上， n個頂點坐標為(xi,yi)，質量為mi，則重心
X = ∑( xi×mi ) / ∑mi
Y = ∑( yi×mi ) / ∑mi
而特殊地，若每個點的質量相同，則
X = ∑xi / n
Y = ∑yi / n
②質量分布均勻(這里沒有給出普適性的公式，因為在此算法中不需要)
特殊地，質量均勻的三角形重心：
X = ( x0 + x1 + x2 ) / 3
Y = ( y0 + y1 + y2 ) / 3
有了這些公式，我們下面要做的就是找到公式中字母在實際問題中代表的意義。
先從整體的思路開始，秉承計算幾何慣有的思維，這里還是要將n多邊形分割成n - 2個三角形，基於三角形求重心的簡便性，我們可以套用質量均勻分布情況下的三角形重心求解公式，此時，n多邊形就轉化成了質量集中在n - 2個三角形重心的圖形了，此時我們再套用質量集中在頂點情況下的普適性公式。
那么公式中的m代表什么呢？它是三角形的面積。我們立體化的看這些平面幾何圖形，假設他們有密度、厚度，肯定是一樣的，那么這里質量就與平面中圖形的面積成正比了。而對於三角形面積的計算，在上面也有所推導證明。
編程實現：有了以上數理思維做鋪墊，在編程實現的時候只需設置窮舉正確的將n多邊形分割成三角形即可。
參考代碼如下。
#include<stdio.h> using namespace std; struct point { double x , y; }; double area(point p1 , point p2 , point p3) { return ((p2.x - p1.x) * (p3.y - p1.y) - (p3.x - p1.x)*(p2.y - p1.y))/2; } int main() { point p1,p2,p3; int t , n; double temp , sumarea , xg , yg; scanf("%d",&t); while(t--) { scanf("%d",&n); sumarea = xg = yg = 0; scanf("%lf%lf%lf%lf",&p1.x , &p1.y , &p2.x , &p2.y); for(int i = 2;i < n;i++) { scanf("%lf%lf",&p3.x , &p3.y); temp = area(p1 , p2 , p3); xg += (p1.x + p2.x + p3.x) * temp; yg += (p1.y + p2.y + p3.y) * temp; sumarea += temp; p2 = p3; } xg = xg / sumarea / 3.0; yg = yg / sumarea / 3.0; printf("%.2lf %.2lf\n" , xg , yg); } return 0; }
再來看一道有關求多邊形重心的問題(Problem source :pku 3855)
Description
Let's have a more formal definition for center of mass (COM). The center of mass for a square, (also circle, and other symmetric shapes) is its center point. And, if a simple shape C is partitioned into two simple shapes A and B with areas SA and SB, then COM(C) (as a vector) can be calculated by
Input
Output
題目大意：依舊是給出頂點數n，然后依次給出n個頂點的坐標，讓你求解n邊形的重心。
編程實現：雖然方法是和上面一題一模一樣的，但是在答案輸出的時卻對精度有着不同的要求。這里如果我們的精度是小數點后6位，而某個答案是-0.0000001，此時如果不作處理會輸出-0.000000，而正確答案應該是0.000000，雖然大小沒變，但是oj返回一個錯誤結果。
所以我們在這里要隨着精度的變化對代碼做出相應的改動。這里要求小數點后6位的精度，我們設置一個常量為0.000001，如果|x|<0.000001，我們就令x = 0.0，這樣就會輸出0.000000而不是-0.000000。
參考代碼如下。
#include<cstdio> #include<cmath> using namespace std; const double eps = 1e-6; struct point { double x , y; }; double area(point p1 , point p2 , point p3) { return ((p2.x - p1.x) * (p3.y - p1.y) - (p3.x - p1.x)*(p2.y - p1.y))/2; } int main() { point p1,p2,p3; int t , n; double temp , sumarea , xg , yg; int tt = 1; while(scanf("%d",&n)!=EOF && n) { sumarea = xg = yg = 0; scanf("%lf%lf%lf%lf",&p1.x , &p1.y , &p2.x , &p2.y); for(int i = 2;i < n;i++) { scanf("%lf%lf",&p3.x , &p3.y); temp = area(p1 , p2 , p3); xg += (p1.x + p2.x + p3.x) * temp; yg += (p1.y + p2.y + p3.y) * temp; sumarea += temp; p2 = p3; } xg = xg / sumarea / 3.0; yg = yg / sumarea / 3.0; if(fabs(xg) < eps) xg = 0.0; if(fabs(yg) < eps) yg = 0.0; printf("Stage #%d: %.6f %.6f\n" , tt++ , xg , yg); } return 0; }
考慮這樣一個問題，給定含有有限個元素的點集，我們如何變成計算能夠覆蓋所有點的最小的圓的圓心和半徑，這也就是所謂的最小圓覆蓋問題。(Problem source: zoj 1450)
Input
The input contains several problems. The first line of each problem is a line containing only one integer N which indicates the number of points to be covered. The next N lines contain N points. Each point is represented by x and y coordinates separated by a space. After the last problem, there will be a line contains only a zero.
Output
For each input problem, you should give a one-line answer which contains three numbers separated by spaces. The first two numbers indicate the x and y coordinates of the result circle, and the third number is the radius of the circle. (use escape sequence %.2f)
編程實現：對於計算幾何的問題，難點多在於編程實現自動化計算，在數理思維上不是很困難。
實現最校園覆蓋的算法有很多，這里介紹一種隨機增量法。
我們將設置三層遍歷，依次遍歷點集。
第一層遍歷，用於我們判斷該點是否在目前構造的圓內，否則的話，此點作為新的圓心。
第二層遍歷，用於我們判斷該點是否在目前構造的圓內，否則的話，此時選出的兩點形成一條線段，中點即為圓心。
第三層遍歷，用於我們判斷該點是否在目前構造的圓內，否則的話，此時遍歷選出的三點的形成的三角形的外心是圓心。
可能有人會疑惑這樣編程的順序，能否實現“最小”圓這一要求。
進行第一層比較，兩個點形成的圓必然比三個點(這三個點中包含前兩個點)形成的圓要小，這是顯而易見的。
進行第二層比較，如果兩個點形成了一個符合要求圓，那它是否是最小的呢？答案也是顯然的，如果我們假設還存在另外一個更小的圓，那么此時開始的兩個點將作為一條弦出現，如果這條弦是直徑，那么這是個等大的圓；如果不是直徑，那么這個圓的直徑將大於這條弦，這與我們初始的假設相悖。
所以我們得到結論，按照上述的過程，一旦構造出符合要求的圓，那么一定是最小的覆蓋圓。
在真正編程實現的時候還需要考慮精度問題，其具體的方法就是盡量減少* /運算的出現，因為其會產生誤差，雖然看似極小，但是數據量一大之后，聚會產生較為明顯差值從而導致錯誤。
這里對於中垂線段的確定，我們給出簡略的證明(針對遍歷三個點是確定外接圓圓心的情況)。
如圖所示，我們在線段ab及其重點ua的基礎上構造出矩形，我們可通過全等證垂直，這樣，中垂線段就用+ - 的形式表示出來，隨后只需求兩條線斷的交點即可求出三點確定的圓的圓心。這樣做相比於通過確立直線方程坐標然后聯立得到方程組，很好的減少了出現精度差的運算次數。
參考代碼如下。
#include<stdio.h> #include<cmath> #include<iostream> using namespace std; const int maxn = 505; const double eps = 1e-6; int n; double r; struct point { double x , y; }p[maxn] , o; double dis(point a , point b) { return sqrt((a.x-b.x)*(a.x-b.x) + (a.y-b.y)*(a.y-b.y)); } point Interesection(point u1,point u2,point v1,point v2) { point ans = u1; double t = ((u1.x - v1.x) * (v1.y - v2.y) - (u1.y - v1.y) * (v1.x - v2.x)) / ((u1.x - u2.x) * (v1.y - v2.y) - (u1.y - u2.y) * (v1.x - v2.x)); ans.x += (u2.x - u1.x) * t; ans.y += (u2.y - u1.y) * t; return ans; //高精度求線段交點，即為外接圓的圓心。 } point Circum(point a , point b , point c) { point ua,ub,va,vb; //求出兩條中垂線段 ua.x = ( a.x + b.x ) / 2; ua.y = ( a.y + b.y ) / 2; ub.x = ua.x - a.y + b.y; ub.y = ua.y + a.x - b.x; va.x = ( a.x + c.x ) / 2; va.y = ( a.y + c.y ) / 2; vb.x = va.x - a.y + c.y; vb.y = va.y + a.x - c.x; return Interesection(ua,ub,va,vb); } point min_center() { int i , j , k; o = p[0]; r = 0; for(i = 1;i < n;i++) //三層循環，因為3點將確定一個圓 { if(dis(o,p[i]) - r > eps ) { o = p[i]; r = 0; for(j = 0;j < i;j++) { if(dis(p[j],o) - r > eps) { o.x = (p[i].x + p[j].x)/2.0; o.y = (p[i].y + p[j].y)/2.0; r = dis(o,p[i]); for(k = 0;k < j;k++) { if(dis(p[k],o) - r > eps) { o = Circum(p[i],p[j],p[k]); r = dis(p[i],o); } } } } } } return o; } int main() { while(scanf("%d",&n) , n) { for(int i = 0;i < n;i++) scanf("%lf%lf",&p[i].x , &p[i].y); min_center(); if(fabs(o.x) < 0.01) o.x = 0.00; if(fabs(o.y) < 0.01) o.y = 0.00; printf("%.2lf %.2lf %.2lf\n",o.x,o.y,r); } }
那么基於對最小圓覆蓋模型的認識，我們再來看一道與之有關的題目(Problem source : hdu4720)
題目大意：給你四個點，用前三個點來確定一個最小覆蓋圓，然后判斷第四個點是否在這個圓內。
數理分析：這里其實是最小圓覆蓋模型的超級簡單版，把覆蓋的n多邊形變成了三角形。這里我們很容易看到，對於直角三角形和銳角三角形，我們要求的圓是外接圓，而對於鈍角三角形，我們要求的則是一個覆蓋圓。但是綜合起來，我們都可以講它們看成求三個點的最小覆蓋圓。
編程實現上沒有什么難度，可以直接使用上面我們學習過的關於最小圓覆蓋的代碼(雖然有點大材小用)。
參考代碼如下。
#include<stdio.h> #include<cmath> #include<iostream> using namespace std; const int maxn = 4; const double eps = 1e-6; int n; double r; struct point { double x , y; }p[maxn] , o; double dis(point a , point b) { return sqrt((a.x-b.x)*(a.x-b.x) + (a.y-b.y)*(a.y-b.y)); } point Interesection(point u1,point u2,point v1,point v2) { point ans = u1; double t = ((u1.x - v1.x) * (v1.y - v2.y) - (u1.y - v1.y) * (v1.x - v2.x)) / ((u1.x - u2.x) * (v1.y - v2.y) - (u1.y - u2.y) * (v1.x - v2.x)); ans.x += (u2.x - u1.x) * t; ans.y += (u2.y - u1.y) * t; return ans; //高精度求線段交點，即為外接圓的圓心。 } point Circum(point a , point b , point c) { point ua,ub,va,vb; //求出兩條中垂線段 ua.x = ( a.x + b.x ) / 2; ua.y = ( a.y + b.y ) / 2; ub.x = ua.x - a.y + b.y; ub.y = ua.y + a.x - b.x; va.x = ( a.x + c.x ) / 2; va.y = ( a.y + c.y ) / 2; vb.x = va.x - a.y + c.y; vb.y = va.y + a.x - c.x; return Interesection(ua,ub,va,vb); } point min_center() { int i , j , k; o = p[0]; r = 0; for(i = 1;i < n;i++) //三層循環，因為3點將確定一個圓 { if(dis(o,p[i]) - r > eps ) { o = p[i]; r = 0; for(j = 0;j < i;j++) { if(dis(p[j],o) - r > eps) { o.x = (p[i].x + p[j].x)/2.0; o.y = (p[i].y + p[j].y)/2.0; r = dis(o,p[i]); for(k = 0;k < j;k++) { if(dis(p[k],o) - r > eps) { o = Circum(p[i],p[j],p[k]); r = dis(p[i],o); } } } } } } return o; } int main() { n = 3; int t; while(~scanf("%d",&t)) { int tt = 1; while(t--) { for(int i = 0;i < n;i++) scanf("%lf%lf",&p[i].x , &p[i].y); scanf("%lf%lf",&p[3].x,&p[3].y); min_center(); printf("Case #%d: ",tt++); if(dis(p[3],o) > r) printf("Safe\n"); else printf("Danger\n"); } } }
基於上面對最小圓覆蓋問題的學習，我們再來看一道應用它的問題。(Problem source : hdu 3932)
題目大意：給出你n個點的坐標，你任取一個點，這個點到到這n個點中的某個點的距離都是存在最大值的，那么現在需要你找到這個最大值最小時的情況。輸出這個點的坐標，並輸出最小的最大距離。
數理分析：我們很明顯的能夠分析到這里需要用到最小圓覆蓋的模型。我們定性的來分析這種最值問題。如果我們n只有2個，那么顯然要求這兩個點構成線段的中點。如果n=3，我們則取這三個點構成三角形的外心。如果這n個點組成的n多邊形的外接圓是存在的，那么這個外接圓的圓心就是我們要找的點。我們看到，為了使這個最大距離盡量小，就使n個點到這個點的距離的差值盡量小些。因為我們發現，如果我們移動該點，使得它靠近了i點，那么相應的會遠離j點(也不盡然，這里是定性的分析)，為了使max(i,j)最小，我們就要一個“均衡受力”的點，所以，這個點應該是外接圓的圓心，因為它到各個點的距離都是一樣的。
那么如果不存在外接圓呢？那我們就要去相應的找最小覆蓋圓了。
有了這層數理分析，我們只要基於求解最小覆蓋圓的模板，就可以很好的解決這個問題了。
參考代碼如下。(Ps:此題目對數據的要求是精確到小數點后一位，四舍五入，但是實際測評中並沒有四舍五入，這里在代碼中就沒有體現四舍五入。)
#include<stdio.h> #include<cmath> #include<iostream> using namespace std; const int maxn = 1005; const double eps = 1e-6; int n , x , y; double r; struct point { double x , y; }p[maxn] , o; double dis(point a , point b) { return sqrt((a.x-b.x)*(a.x-b.x) + (a.y-b.y)*(a.y-b.y)); } point Interesection(point u1,point u2,point v1,point v2) { point ans = u1; double t = ((u1.x - v1.x) * (v1.y - v2.y) - (u1.y - v1.y) * (v1.x - v2.x)) / ((u1.x - u2.x) * (v1.y - v2.y) - (u1.y - u2.y) * (v1.x - v2.x)); ans.x += (u2.x - u1.x) * t; ans.y += (u2.y - u1.y) * t; return ans; //高精度求線段交點，即為外接圓的圓心。 } point Circum(point a , point b , point c) { point ua,ub,va,vb; //求出兩條中垂線段 ua.x = ( a.x + b.x ) / 2; ua.y = ( a.y + b.y ) / 2; ub.x = ua.x - a.y + b.y; ub.y = ua.y + a.x - b.x; va.x = ( a.x + c.x ) / 2; va.y = ( a.y + c.y ) / 2; vb.x = va.x - a.y + c.y; vb.y = va.y + a.x - c.x; return Interesection(ua,ub,va,vb); } point min_center() { int i , j , k; o = p[0]; r = 0; for(i = 1;i < n;i++) //三層循環，因為3點將確定一個圓 { if(dis(o,p[i]) - r > eps ) { o = p[i]; r = 0; for(j = 0;j < i;j++) { if(dis(p[j],o) - r > eps) { o.x = (p[i].x + p[j].x)/2.0; o.y = (p[i].y + p[j].y)/2.0; r = dis(o,p[i]); for(k = 0;k < j;k++) { if(dis(p[k],o) - r > eps) { o = Circum(p[i],p[j],p[k]); r = dis(p[i],o); } } } } } } return o; } int main() { while(~scanf("%d%d%d",&x,&y,&n)) { for(int i = 0;i < n;i++) scanf("%lf%lf",&p[i].x , &p[i].y); min_center(); printf("(%.1lf,%.1lf).\n",o.x,o.y); printf("%.1lf\n",r); } }
在了解了最小矩陣覆蓋、凸包、最小圓覆蓋等覆蓋類問題后，我們在這里再來看一道有關圖形覆蓋類問題。(Problem source 1077)
題目大意：給出你n個平面內的點的坐標，然后讓你用半徑為1的圓去圈點，讓你求出能夠圈出的最大的點。
數理分析：這道題在大的思路上肯定是要暴力窮舉，這不難看出來，剩下我們要解決的就是針對每次窮舉我們應該怎樣處理。
我們窮舉要選取的是兩個點，為什么是兩個點呢？因為我們通過這兩個點，再加上圓的半徑，就可以確定單位圓的位置。什么？你說可能有兩邊的？那么你想對了，這的確是這道題目需要着力思考的一個地方。
我們先不糾結於細節，單純的拿出一個圖來進行分析。
通過畫這一個圖，我們可以簡單的列出每次窮舉的時候需要用到的計算公式，原理也十分簡單，是最淺顯的幾何知識。
那么現在的問題就是選取哪個半圓的問題。選取了半圓，我們用正余弦計算正負號的問題。這里做一個簡單的思維鋪墊，我們假定選取的兩個點是有順序的，即取的兩個點A、B是有方向的,A->B，這樣的話，我們在選圓的時候只需要選一個即可，因為在后面遍歷的時候還會出現B->A，此時選出的圓將會對應着上次選擇沒有選到的那個圓，這要就做到了不重不漏。
既然這樣，我們開始基於一個點分析，加入我們選取的第一個點是A，那么就確立了關系A -> ? ，其中'?'就是我們下個要選取的字母。此時，其他字母相對於A有如下四種情況(之所以要分析這么詳細是選取圓時，為了計算圓心，基於選取兩點中點進行正余弦運算的正負號是不確定的)
其實細分到這，比較懶的一種做法就是四種情況分情況，對應四段代碼，但是通常情況下呈現出的符號是有歸納性的規律的。
這里答案其實呈現出多樣性，我只就一種情況進行分析。
在這四種情況中，dy / dx的正負我已標在圖中。對於②③圖，如果我們都選擇直線AB左半平面的圓，此時取ang = arctan(-dy/dx)，那么這兩個圖呈現出-x，+y的運算規律，考慮到第四象限角的正余弦值，顯然我們基於AB中點o，要進行-操作。
同理，我們分析①④，選取右平面，就會得到同上的結論。
這樣我們就很好的將四種情況統一到了一起，並且不會漏不會重復。
有了如上的數理思維，編程就很容易了。這里還有一個容易忽視的細節就是n = 1的臨界情況分析，顯然這種情況怎樣都是可以圈出一個的，所以我們在進行暴力枚舉的時候，應該注意比較大小的初始變量ans應該賦為1.
參考代碼如下。
#include<stdio.h> #include<cmath> using namespace std; const int maxn = 305; const double r = 1.0; const double eps = 1e-6; struct point { double x ,y; }p[maxn]; int n; double dis(int a , int b) { return (p[a].x-p[b].x)*(p[a].x-p[b].x) + (p[a].y-p[b].y)*(p[a].y-p[b].y); } void get_centre(int a,int b) { double ox,oy,dx,dy,r,temp; ox=(p[a].x+p[b].x)/2; oy=(p[a].y+p[b].y)/2; dx=p[b].x-p[a].x; dy=p[b].y-p[a].y; p[n].x=ox;p[n].y=oy; temp=dis(n,b); r=sqrt(1.0-temp); if(fabs(dx) < eps)//豎直的情況 p[n].x-=r; else { double ang=atan(-dy/dx); p[n].x-=r*sin(ang); //這里給出的情況呼應上文的分析 p[n].y-=r*cos(ang); } } int main() { int T; int i , j , k ; scanf("%d",&T); while(T--) { scanf("%d",&n); for(i = 0;i < n;i++) scanf("%lf%lf",&p[i].x,&p[i].y); int Max , ans = 1; for(i = 0;i < n;i++) { for(j = i + 1;j < n;j++) { if(dis(i,j) >= 4) continue; get_centre(i , j); for(k = 0,Max = 0;k < n;k++) { if(n - k + Max <= ans) break; double temp = sqrt(dis(k,n)); if(temp <= 1.0001) Max++; } if(Max > ans) ans = Max; } } printf("%d\n",ans); } }
基於上面對於單位圓點覆蓋的問題的探討，我們再來看一道類似的題目。(Problem source : pku 1379)
題目闡述的很直接，就是給定一個含n個元素的點集，讓你求找到用單位圓可以覆蓋的最多的點是多少。
基於對上題模型的思考，這里就可以簡單的修改一下上面代碼來實現了，這里就不再累述。
參考系：《計算幾何及應用》 金博 郭立