计算几何

此文章是一些自己未接触过的计算几何知识，不是太明白。。。。。

zoj 1450 && hdu 3007     本题代码来自 himdd大牛  ：http://blog.himdd.com/archives/2666

/**
    最小圆覆盖 （随机增量法）
    求给你n个点，让你用一个最小的圆覆盖所有的点，求这个圆的坐标&&半径
**/
#include <iostream>
#include <cstdio>
#include <ctime>
#include <cmath>
#include <algorithm>
using namespace std;
const double eps = 1e-12;
const int len = 100010;
struct Point
{
    double x,y;
} p[len];
struct Line
{
    Point a,b;
};
int dbcmp(double n)
{
    return n < -eps ? -1 : n > eps;
}
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 intersection(Line u,Line v)
{
    Point ret=u.a;
    double t=((u.a.x-v.a.x)*(v.b.y-v.a.y)-(u.a.y-v.a.y)*(v.b.x-v.a.x))
             /((u.a.x-u.b.x)*(v.b.y-v.a.y)-(u.a.y-u.b.y)*(v.b.x-v.a.x));
    ret.x+=(u.b.x-u.a.x)*t;
    ret.y+=(u.b.y-u.a.y)*t;
    return ret;
}
//三角形外接圆圆心(外心)
Point center(Point a,Point b,Point c)
{
    Line u,v;
    u.a.x=(a.x+b.x)/2;
    u.a.y=(a.y+b.y)/2;
    u.b.x=u.a.x+(u.a.y-a.y);
    u.b.y=u.a.y-(u.a.x-a.x);
    v.a.x=(a.x+c.x)/2;
    v.a.y=(a.y+c.y)/2;
    v.b.x=v.a.x+(v.a.y-a.y);
    v.b.y=v.a.y-(v.a.x-a.x);
    return intersection(u,v);
}

void min_cir(Point * p, int n, Point & cir, double &r)
{
    random_shuffle(p, p + n);
    cir = p[0];
    r = 0;
    for(int i = 1; i < n; ++i)
    {
        if(dbcmp(dis(p[i],cir) -r) <= 0)continue;
        cir = p[i];
        r = 0;
        for(int j =0; j < i ; ++j)
        {
            if(dbcmp(dis(p[j], cir) -r) <= 0 )continue;
            cir.x = (p[i].x + p[j].x) /2;
            cir.y = (p[i].y + p[j].y) /2;
            r = dis( cir, p[j]);
            for(int k =0; k < j; ++k)
            {
                if(dbcmp( dis(p[k], cir) -r) <= 0) continue;
                cir = center(p[i], p[j], p[k]);
                r = dis( cir, p[k]);
            }

        }
    }
}
int main()
{
    int n;
    while( ~scanf("%d", &n), n)
    {
        for (int i = 0; i < n; ++i)
        {
            scanf("%lf%lf", &p[i].x, &p[i].y);
        }
        Point cir;
        double r;
        min_cir(p, n, cir, r);
        printf("%.2lf %.2lf %.2lf\n", cir.x, cir.y, r);
    }

}

POJ 3608 Bridge Across Islands

题意：求两凸包之间的距离；摘自 kuangbin 大牛：http://www.cnblogs.com/kuangbin/p/3221911.html

#include <stdio.h>
#include <algorithm>
#include <iostream>
#include <string.h>
#include <math.h>
using namespace std;

const double eps = 1e-8;
int sgn(double x)
{
    if(fabs(x) < eps)return 0;
    if(x < 0)return -1;
    else return 1;
}
struct Point{
    double x,y;
    Point(double _x = 0.0,double _y = 0.0){
        x = _x;y = _y;
    }
    Point operator -(const Point &b)const{
        return Point(x - b.x, y - b.y);
    }
    double operator ^(const Point &b)const{
        return x*b.y - y*b.x;
    }
    double operator *(const Point &b)const{
        return x*b.x + y*b.y;
    }
    void input(){
        scanf("%lf%lf",&x,&y);
    }
};
struct Line
{
    Point s,e;
    Line(){}
    Line(Point _s,Point _e){
        s = _s;
        e = _e;
    }
};
double dist(Point a,Point b){
    return sqrt((a-b)*(a-b));
}
Point NearestPointToLineSeg(Point P,Line L){
    Point result;
    double t = ((P-L.s)*(L.e-L.s))/((L.e-L.s)*(L.e-L.s));
    if(t >= 0 && t <= 1){
        result.x = L.s.x + (L.e.x - L.s.x)*t;
        result.y = L.s.y + (L.e.y - L.s.y)*t;
    }
    else{
        if(dist(P,L.s) < dist(P,L.e))
            result = L.s;
        else result = L.e;
    }
    return result;
}
/*
 * 求凸包，Graham算法
 * 点的编号0~n-1
 * 返回凸包结果Stack[0~top-1]为凸包的编号
 */
const int MAXN = 10010;
Point list[MAXN];
int Stack[MAXN],top;
//相对于list[0]的极角排序
bool _cmp(Point p1,Point p2){
    double tmp = (p1-list[0])^(p2-list[0]);
    if(sgn(tmp) > 0)return true;
    else if(sgn(tmp) == 0 && sgn(dist(p1,list[0]) - dist(p2,list[0])) <= 0)
        return true;
    else return false;
}
void Graham(int n){
    Point p0;
    int k = 0;
    p0 = list[0];
    //找最下边的一个点
    for(int i = 1;i < n;i++){
        if( (p0.y > list[i].y) || (p0.y == list[i].y && p0.x > list[i].x) ){
            p0 = list[i];
            k = i;
        }
    }
    swap(list[k],list[0]);
    sort(list+1,list+n,_cmp);
    if(n == 1){
        top = 1;
        Stack[0] = 0;
        return;
    }
    if(n == 2){
        top = 2;
        Stack[0] = 0;
        Stack[1] = 1;
        return ;
    }
    Stack[0] = 0;
    Stack[1] = 1;
    top = 2;
    for(int i = 2;i < n;i++){
        while(top > 1 && sgn((list[Stack[top-1]]-list[Stack[top-2]])^(list[i]-list[Stack[top-2]])) <= 0)
            top--;
        Stack[top++] = i;
    }
}
//点p0到线段p1p2的距离
double pointtoseg(Point p0,Point p1,Point p2){
    return dist(p0,NearestPointToLineSeg(p0,Line(p1,p2)));
}
//平行线段p0p1和p2p3的距离
double dispallseg(Point p0,Point p1,Point p2,Point p3){
    double ans1 = min(pointtoseg(p0,p2,p3),pointtoseg(p1,p2,p3));
    double ans2 = min(pointtoseg(p2,p0,p1),pointtoseg(p3,p0,p1));
    return min(ans1,ans2);
}
//得到向量a1a2和b1b2的位置关系
double Get_angle(Point a1,Point a2,Point b1,Point b2){
    Point t = b1 - ( b2 - a1 );
    return (a2-a1)^(t-a1);
}
//旋转卡壳，求两个凸包的最小距离
double rotating_calipers(Point p[],int np,Point q[],int nq){
    int sp = 0, sq = 0;
    for(int i = 0;i < np;i++)
        if(sgn(p[i].y - p[sp].y) < 0)
            sp = i;
    for(int i = 0;i < nq;i++)
        if(sgn(q[i].y - q[sq].y) > 0)
            sq = i;
    double tmp;
    double ans = 1e99;
    for(int i = 0;i < np;i++){
        while(sgn(tmp = Get_angle(p[sp],p[(sp+1)%np],q[sq],q[(sq+1)%nq])) < 0 )
            sq = (sq + 1)%nq;
        if(sgn(tmp) == 0)
            ans = min(ans,dispallseg(p[sp],p[(sp+1)%np],q[sq],q[(sq+1)%nq]));
        else ans = min(ans,pointtoseg(q[sq],p[sp],p[(sp+1)%np]));
        sp = (sp+1)%np;
    }
    return ans;
}

double solve(Point p[],int n,Point q[],int m){
    return min(rotating_calipers(p,n,q,m),rotating_calipers(q,m,p,n));
}
Point p[MAXN],q[MAXN];
int main()
{
    int n,m;
    while(scanf("%d%d",&n,&m)==2)
    {
        if(n == 0 && m == 0)break;
        for(int i = 0;i < n;i++)
            list[i].input();
        Graham(n);
        n = top;
        for(int i = 0;i < n;i++)
            p[i] = list[Stack[i]];
        for(int i = 0;i < m;i++)
            list[i].input();
        Graham(m);
        m = top;
        for(int i = 0;i < m;i++)
            q[i] = list[Stack[i]];
        printf("%.5lf\n",solve(p,n,q,m));
    }
    return 0;
}