HDU 3694 Fermat Point in Quadrangle(四邊形的費馬點)

Fermat Point in Quadrangle

Time Limit: 2000/1000 MS (Java/Others)    Memory Limit: 65536/32768 K (Java/Others)
Total Submission(s): 1460    Accepted Submission(s): 241


Problem DescriptionIn geometry the Fermat point of a triangle, also called Torricelli point, is a point such that the total distance from the three vertices of the triangle to the point is the minimum. It is so named because this problem is first raised by Fermat in a private
letter. In the following picture, P0 is the Fermat point. You may have already known the property that:

Alice and Bob are learning geometry. Recently they are studying about the Fermat Point. 

Alice: I wonder whether there is a similar point for quadrangle.

Bob: I think there must exist one.

Alice: Then how to know where it is? How to prove?

Bob: I don’t know. Wait… the point may hold the similar property as the case in triangle. 

Alice: It sounds reasonable. Why not use our computer to solve the problem? Find the Fermat point, and then verify your assumption.

Bob: A good idea.

So they ask you, the best programmer, to solve it. Find the Fermat point for a quadrangle, i.e. find a point such that the total distance from the four vertices of the quadrangle to that point is the minimum.

 

InputThe input contains no more than 1000 test cases.

Each test case is a single line which contains eight float numbers, and it is formatted as below:

x1 y1 x2 y2 x3 y3 x4 y4

xi, yi are the x- and y-coordinates of the ith vertices of a quadrangle. They are float numbers and satisfy 0 ≤ xi ≤ 1000 and 0 ≤ yi ≤ 1000 (i = 1, …, 4).

The input is ended by eight -1.

 

OutputFor each test case, find the Fermat point, and output the total distance from the four vertices to that point. The result should be rounded to four digits after the decimal point. 

Sample Input

0 0 1 1 1 0 0 11 1 1 1 1 1 1 1-1 -1 -1 -1 -1 -1 -1 -1

 

Sample Output

2.82840.0000

 

Source2010 Asia Fuzhou Regional Contest 

Recommendaxubiao 這個題目wa的我有點傷心，偶爾犯一些小錯誤，在別的類型題目可能是容易發現的，但是在計算幾何裡面不但不容易發現而且是致命的拿到這個題目想都沒想類比退火的思路，很快寫出代碼（下下面會有） 老是過不了，不知道為什麼四變形就不可以！然後就想這樣子做首先按照點排序，和凸包一樣的排序，判斷四點是否共線，不成立後判斷是否是凸包！我的解決方案有點極端：首先四個點枚舉，每一個點作為費馬點一次求出距離，這個是假設不是凸包的情況下然後在求是凸包的情況下，防止錯誤兩者取最小的保證正確也就是說費馬點不是對角線的交點就是四個點中的某一個，嘗試一下取最下至於下面代碼為什麼這麼長，因為老是wa老是改，就成這樣了！

#include <iostream>#include <string.h>#include <stdio.h>#include <math.h>#include <algorithm>using namespace std;struct point{double x;double y;}po[100],my_p;struct line{point a;point b;};double dis(point &a,point &b){return sqrt((a.x-b.x)*(a.x-b.x) + (a.y-b.y)*(a.y-b.y));}double across(point &a,point &b,point &c){return (b.x-a.x)*(c.y-a.y)-(b.y-a.y)*(c.x-a.x);}int cmp(const void *a,const void *b){return across(po[0],*(point*)a,*(point*)b) > 1e-8 ? -1 : 1;}int is_convex(int num,point *p){int i;p[num]=p[0];num++;p[num]=po[1];num++;int count=0;for(i=2;i < num;i++){if(across(p[i-2],p[i-1],p[i]) < 1e-8){my_p=p[i-1];return 0;}}return 1;}point intersection(line &u,line &v){point ret=u.a;double t=((u.a.x-v.a.x)*(v.a.y-v.b.y) - (u.a.y-v.a.y)*(v.a.x-v.b.x))/((u.a.x-u.b.x)*(v.a.y-v.b.y)-(u.a.y-u.b.y)*(v.a.x-v.b.x));ret.x+=(u.b.x-u.a.x)*t;ret.y+=(u.b.y-u.a.y)*t;return ret;}bool zero(double a)//判斷結果是否為0{return fabs(a) <= 1e-8;}bool parallel(point &u1,point &u2,point &v1,point &v2)//判斷直線u1 u2 and v1 v2 是否平行{return zero((u1.x-u2.x)*(v1.y-v2.y)-(v1.x-v2.x)*(u1.y-u2.y));}bool equal(point &a,point &b){return fabs(a.x-b.x) < 1e-8 && fabs(a.y-b.y) < 1e-8;}double MIN(double a,double b){return a < b ? a : b;}int main(){int i,j,k;point temp;point my_temp;line a,b;double ans,ans1,tt;int pos;while(1){scanf("%lf%lf",&po[0].x,&po[0].y);temp=po[0];pos=0;for(i=1;i<4;i++){scanf("%lf%lf",&po[i].x,&po[i].y);if(temp.y > po[i].y)temp=po[i],pos=i;}if(po[0].x==-1 && po[0].y==-1 && po[1].x==-1 && po[1].y==-1 && po[2].x==-1 && po[2].y==-1 && po[3].x==-1 && po[3].y==-1)return 0;my_temp=po[0];po[0]=po[pos];po[pos]=my_temp;qsort(po+1,3,sizeof(po[0]),cmp);ans=ans1=100000000;tt=0;for(i=0;i<4;i++){tt=0;for(j=0;j<4;j++)tt+=dis(po[i],po[j]);ans1=MIN(tt,ans1);}if(parallel(po[0],po[1],po[1],po[2]) && parallel(po[1],po[2],po[2],po[3]))NULL;else if(is_convex(4,po)){a.a=po[0],a.b=po[2];b.a=po[1],b.b=po[3];ans=0;my_p=intersection(a,b);for(i=0;i<4;i++)ans+=dis(my_p,po[i]);//printf("%lf %lf \n",my_p.x,my_p.y);}printf("%.4lf\n",MIN(ans,ans1));}return 0;}

類比退火代碼：

#include <iostream>#include <stdio.h>#include <cmath>#include <string.h>using namespace std;#define inf 1e-10struct point{double x;double y;}po[4];double dis(point &a,point &b){return (a.x-b.x)*(a.x-b.x) + (a.y-b.y)*(a.y-b.y);}double find_ans(){int i;point temp,p,q;double ans=0;temp=p=q=po[0];for(i=1;i<4;i++)ans+=dis(temp,po[i]);double step=1500,rec;while(step > inf){q=p;temp=q;rec=0;temp.x+=step;for(i=0;i<4;i++)rec+=dis(temp,po[i]);if(rec < ans){p=temp;ans=rec;}temp=q;rec=0;temp.x-=step;for(i=0;i<4;i++)rec+=dis(temp,po[i]);if(rec < ans){p=temp;ans=rec;}temp=q;rec=0;temp.y+=step;for(i=0;i<4;i++)rec+=dis(temp,po[i]);if(rec < ans){p=temp;ans=rec;}temp=q;rec=0;temp.y-=step;for(i=0;i<4;i++)rec+=dis(temp,po[i]);if(rec < ans){p=temp;ans=rec;}step/=2;}return ans;}int main(){int i,j,k;while(1){for(i=0;i<4;i++)scanf("%lf%lf",&po[i].x,&po[i].y);if(po[0].x==-1 && po[0].y==-1 && po[1].x==-1 && po[1].y==-1 && po[2].x==-1 && po[2].y==-1 && po[3].x==-1 && po[3].y==-1)return 0;printf("%.4lf\n",find_ans());}return 0;