计算几何之模拟退火
1.POJ1379
HDU 1109题面相同
题意：给定0——X,0——Y的矩形，给定n个点的坐标，在矩形中求得一个点是的该点到所有点的最短距离最大。(保留小数 )
分析：模拟退火(随机化)算法，可求解精度要求较小的几何寻点问题.
随机选取多个(20个)初始点，进行多次随机坐标变换，并且根据降温概率进行选择是否转移次优解。
随机坐标变换：随机半径，随机角度进行改变坐标.
#include<cstdio>
#include<algorithm>
#include<cstring>
#include<cmath>
#include<ctime>
using namespace std;
typedef double db;
const int maxn=1e3+10;
const int maxm=20;
const db PI=acos(-1.0);
db X,Y;
int n;
db dis[100];
struct point{
db x,y;
point(){}
point( db _x,db _y ):x(_x),y(_y){ }
point operator - ( point p ) { return point(x-p.x,y-p.y) ; }
db dis() { return sqrt(x*x+y*y); }
}a[maxn],ans[maxn];
db getMinDis( point p )
{
double res=1e9;
for( int i=0;i<n;i++ ) res=min(res,(p-a[i]).dis() );
return res;
}
point SA()
{
for( int i=0;i<maxm;i++ )
{
ans[i].x=( rand()%1000+1 )/1000.0*X;
ans[i].y=( rand()%1000+1 )/1000.0*Y;
dis[i]=getMinDis(ans[i]);
}
db delta=max(X,Y)/sqrt(n*1.0),d=0.9;
while( delta>0.01 )
{
for( int i=0;i<maxm;i++ )
{
for( int j=0;j<maxm;j++ )
{
point t=ans[i];
db angle=rand()%1000/1000.0*2*PI;
db dx=delta*cos(angle)*(rand()%1000/1000.0);
db dy=delta*sin(angle)*(rand()%1000/1000.0);
t.x+=dx;t.y+=dy;
if( t.x<0 || t.x>X || t.y<0 || t.y>Y ) continue;
db tmp=getMinDis(t);
db delta2=tmp-dis[i];
if( delta2>0 ) dis[i]=tmp,ans[i]=t;
else if( tmp>dis[i] && exp( -delta2/delta )*RAND_MAX>rand() );
else t.x-=dx,t.y-=dy;
}
}
delta*=d;
}
db dd=0;int pp=0;
for( int i=0;i<maxm;i++ )
{
if( dis[i]>dd ) dd=dis[i],pp=i;
}
return ans[pp];
}
int main()
{
srand(time(0));
int T;
scanf("%d",&T);
while( T-- )
{
scanf("%lf%lf%d",&X,&Y,&n);
for( int i=0;i<n;i++ ) scanf("%lf%lf",&a[i].x,&a[i].y);
point p=SA();
printf("The safest point is (%.1f, %.1f).\n", p.x, p.y);
}
}2.POJ 2420
题意：给定n个点的坐标，求是否存在一点到n个点的距离之和最小.输出最小距离之和.(n<=100,0<=xi、yi<=1000 )
分析：该点一定是在n个点环绕的封闭区域中，这道题我们只用取横纵坐标的最大最小值来限制随机点的坐标,然后根据模拟退火进行随机搜索.
注:这道题是多组读入，降温系数d不可太大容易TLE.
#include<cstdio>
#include<algorithm>
#include<cstring>
#include<cmath>
#include<ctime>
using namespace std;
typedef double db;
const int maxn=1e3+10;
const int maxm=30;
const db PI=acos(-1.0);
db X,Y,X1,Y1;
int n;
db dis[100];
struct point{
db x,y;
point(){}
point( db _x,db _y ):x(_x),y(_y){ }
point operator - ( point p ) { return point(x-p.x,y-p.y) ; }
db dis() { return sqrt(x*x+y*y); }
}a[maxn],ans[maxn];
db getMinDis( point p )
{
double res=0;
for( int i=0;i<n;i++ ) res+=(p-a[i]).dis();
return res;
}
db SA()
{
for( int i=0;i<maxm;i++ )
{
ans[i].x=( rand()%1000+1 )/1000.0*X;
ans[i].y=( rand()%1000+1 )/1000.0*Y;
dis[i]=getMinDis(ans[i]);
}
db delta=max(X-X1,Y-Y1),d=0.5;
while( delta>0.01 )
{
for( int i=0;i<maxm;i++ )
{
for( int j=0;j<maxm;j++ )
{
point t=ans[i];
db angle=rand()%1000/1000.0*2*PI;
db dx=delta*cos(angle)*(rand()%1000/1000.0);
db dy=delta*sin(angle)*(rand()%1000/1000.0);
t.x+=dx;t.y+=dy;
if( t.x<X1 || t.x>X || t.y<Y1 || t.y>Y ) continue;
db tmp=getMinDis(t);
db delta2=tmp-dis[i];
if( delta2<0 ) dis[i]=tmp,ans[i]=t;
else if( tmp>dis[i] && exp( -delta2/delta )*RAND_MAX>rand() );
else t.x-=dx,t.y-=dy;
}
}
delta*=d;
}
db res=1e12;
for( int i=0;i<maxm;i++ ) res=min(res,dis[i]);
return res;
}
int main()
{
srand(time(0));
while( ~scanf("%d",&n) )
{
X1=1e9,Y1=1e9;
for( int i=0;i<n;i++ )
{
scanf("%lf%lf",&a[i].x,&a[i].y);
X=max(X,a[i].x);Y=max(Y,a[i].y);
X1=min(X1,a[i].x);Y1=min(Y1,a[i].y);
}
printf("%.0lf\n",SA());
}
}
3.POJ3466
题意：给定n个顶点围成的多边形，求在多边形内部是否能放下半径为x的圆.
#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <algorithm>
using namespace std;
const int maxn = 100;
const double eps = 1e-3;
const int times = 5;
const double INF = 99999999;
int n;
struct point
{
double x, y, r;
} p[maxn], test[maxn];
inline int dblcmp( double x )
{
if( fabs(x) < eps )
return 0;
return x > 0 ? 1 : -1;
}
inline double sqr( double x )
{
return x*x;
}
double dis( point& aa, point& bb )
{
return sqrt(sqr(aa.x-bb.x)+sqr(aa.y-bb.y));
}
double cross( point& k, point& aa, point& bb )
{
return (k.x-aa.x)*(k.y-bb.y)-(k.y-aa.y)*(k.x-bb.x);
}
//线段到点的距离
double seg_point_dis( point& l1, point& l2, point& k )
{
double a, b, c;
a = dis(l1, k);
b = dis(l2, k);
c = dis(l1, l2);
//这里的判断是由余弦定理推出来的
if( dblcmp(a*a+c*c-b*b) < 0 )
return a;
else if( dblcmp(b*b+c*c-a*a) < 0 )
return b;
else
return fabs(cross(k, l1, l2)/c);
}
//判断点是否在多边形内
//思想：过点p做一条射线，判断交点数
bool point_inside( point& aa )
{
int i, cnt = 0;
double t;
//确保p[n] = p[0]
for( i = 0; i < n; ++i )
{
if( (p[i].y <= aa.y && p[i+1].y > aa.y) ||
(p[i+1].y <= aa.y && p[i].y > aa.y) )
{
if( !dblcmp(p[i].y-p[i+1].y) )
{
if( dblcmp(p[i].y-aa.y) == 0 )
cnt++;
t = -INF;
}
else
t = p[i+1].x - (p[i+1].x-p[i].x)*(p[i+1].y-aa.y)/(p[i+1].y-p[i].y);
if( dblcmp( t - aa.x ) >= 0 )
cnt++;
}
}
return cnt%2;
}
void cal( point& aa )
{
double t;
aa.r = INF;
for( int i = 0; i < n; ++i )
{
t = seg_point_dis(p[i], p[i+1], aa);
aa.r = min(aa.r, t);
}
}
int main()
{
int i, j;
double R, delte, maxx, maxy, minx, miny, ang;
point temp;
bool ok;
while( scanf("%d", &n), n )
{
ok = 0;
maxx = maxy = 0;
minx = miny = INF;
for( i = 0; i < n; ++i )
{
scanf("%lf %lf", &p[i].x, &p[i].y);
maxx = max(maxx, p[i].x);
maxy = max(maxy, p[i].y);
minx = min(minx, p[i].x);
miny = min(miny, p[i].y);
}
p[n] = p[0];
for( i = 0; i < n; ++i )
{
test[i].x = (p[i].x+p[i+1].x)/2;
test[i].y = (p[i].y+p[i+1].y)/2;
test[i].r = 0;
}
scanf("%lf", &R);
maxx -= minx; maxy -= miny;
delte = sqrt(maxx*maxx+maxy*maxy)/2;
while( !ok && delte > eps )
{
for( i = 0; !ok && i < n; ++i )
{
for( j = 0; !ok && j < times; ++j )
{
ang = rand();
temp.x = test[i].x + delte*cos(ang);
temp.y = test[i].y + delte*sin(ang);
if( point_inside(temp) )
{
cal(temp);
if( temp.r > test[i].r )
{
test[i] = temp;
if( dblcmp(temp.r-R) >= 0 )
ok = 1;
}
}
}
}
delte *= 0.90;
}
if( ok )
printf("Yes\n");
else
printf("No\n");
}
return 0;
}搁着，暂时补不动....
