RIA算法解决最小覆盖圆问题

一.概念引入

        最小包围圆问题：对于给定的平面上甩个点所组成的一个集合P，求出P的最小包围圆，即包含P中所有点、半径最小的那个圆。也就是求出这个最小
包围圆的圆心位置和半径。

        下面是若干性质。

• 有限点集P的最小包围圆是唯一的。这里约定，若P中只有一个点v，则最小包围圆是退化的，其半径为0，圆心为点v。
• 非退化最小包围圆可以由2个或者3个边界点定义。边界上只有两个点，则必定是直径两端，其它点都在圆内部，这个咱就不证明了。
• 点集P中，距离最大的2个点A、B不一定都在边界上，但是必有d≥|AB|，笔者认为这点很重要。

• 直角三角形或钝角三角形的3个顶点的最小包围圆是以最长边为直径的圆；锐角三角形3个顶点的最小包围圆是三角形的外接圆。
• 新加入点一定在圆上

二.算法实现

           加上shuffle后，ZOJ1450第二组数据结果不稳定，不加的话全部正确，不过两者都还是WA，莫非求圆心也没错(水平竖直斜向都可以求出)，别人的C++的都AC（用了random_shuffle()函数了，原来判断点在圆心内写错了，网上找不到java代码，不管了……

package a;
import java.util.Random;
import java.util.Scanner;

/*
 * hdu只有500点，可以直接两点间最大距离暴力试试
 */
public class HDU3007 {

    public static void main(String[] args) {
        new RIA().go();
    }
}

class Point {
    double x;
    double y;
    
    public Point() {
        this.x = 0;
        this.y = 0;
    }
}

class Line {
    Point a;
    Point b;
    
    public Line() {
        this.a = new Point();
        this.b = new Point();
    }
    public Line(Point a,Point b) {
        this.a = a;
        this.b = b;
    }
    
    //求两直线的交点，斜率相同的话res=u.a
    Point intersection(Line u,Line v){
        Point res = 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));
        res.x += (u.b.x-u.a.x)*t;
        res.y += (u.b.y-u.a.y)*t;
        return res;
    }
    
    //三角形外接圆圆心(外心)
//    Point center(Point a,Point b,Point c) {
//        //加上这个才没有编译器提示未初始化，因为new所以也写了构造方法
//        Line u = new Line(),v = new Line();
//        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);
//    }
    Point center(Point a,Point b,Point c) {
        Point ret = new Point();
        double a1=b.x-a.x, b1=b.y-a.y, c1=(a1*a1+b1*b1)/2;
        double a2=c.x-a.x, b2=c.y-a.y, c2=(a2*a2+b2*b2)/2;
        double d = a1*b2 - a2*b1;
        ret.x = a.x + (c1*b2-c2*b1)/d;
        ret.y = a.y + (a1*c2-a2*c1)/d;
        return ret;
    }
}

class RIA {
    int n;
    double x;
    double y;
    
    public void go() {
        Scanner sc = new Scanner(System.in);
        while(true) {
            n = sc.nextInt();
            if(0==n) 
                break;
            Point point[] = new Point[n];
            for(int i=0; i<n; i++) {//不加的话空指针异常
                point[i] = new Point();
            }
            for(int i=0; i<n; i++) {
                x = sc.nextDouble();
                y = sc.nextDouble();
                point[i].x = x;
                point[i].y = y;
            }
            //shuffle(point);
            solve(point);
        }
    }
    
    private void shuffle(Point[] point) {
        
        for(int i=0; i<point.length; i++) {
            //Random r = new Random();
            //int j = r.nextInt(point.length);
            int j = (int)(Math.random()*point.length);
            if(i!=j) {
                Point temp = point[i];
                point[i] = point[j];
                point[j] = temp;
            }
        }
    }

    private void solve(Point[] point) {
        
        Point circle = point[0];
        double r = 0;
        
        for(int i=1; i<n; i++) {
            double dis = distance(circle, point[i]);
            if(Double.compare(dis, r)<=0) {
                continue;
            }
            circle = point[i];
            r = 0;
            for(int j=0; j<i; j++) {
                dis = distance(circle, point[j]);
                if(Double.compare(dis, r)<=0) {
                    continue;
                }
                circle.x = (point[j].x + point[i].x)/2;
                circle.y = (point[j].y + point[i].y)/2;
                r = distance(circle, point[j]);
                
                for(int k=0; k<j; k++) {
                    dis = distance(circle, point[k]);
                    if(Double.compare(dis, r)<=0) {
                        continue;
                    }
                    Line line = new Line();
                    circle = line.center(point[i],point[j],point[k]);
                    r = distance(point[k], circle);
                }
                
            }
                
        }
        //没有lf只说
        System.out.println(String.format("%.2f", circle.x) + 
                " "+String.format("%.2f", circle.y)+
                " "+String.format("%.2f", r));
    //这样不行，若是初试不足三位，那么输出就不够三位
//        System.out.println((double)Math.round(circle.x*100)/100 + 
//                " "+(double)Math.round(circle.y*100)/100+
//                " "+(double)Math.round(r*100)/100);
        
    }

    public double distance(Point p1, Point p2) {
        return (Math.hypot((p1.x - p2.x), (p1.y - p2.y)));
    }
}

三.若干思考

        RIA算法叫随机增量法，加入随机性后复杂度是线性的(表示目前不太理解)，昨晚又想了想，第一层循环是产生新加入的点，由性质知该点必须在圆上，所以三层循环里每层都有point[i]去组成圆(第一层中是退化的)；

        第一层中为什么半径是0呢？和圆心是point[i]一样，笔者认为主要是为让第二层一定进行下去(略过if判断)，或者就认为此时只有一个点是退化圆。

        如何保证最小？因为每次都是最小的(看倒数第二条性质)，所以结果是最小的。

四.浮点数

        Double.compare(p,q)，若是和0比，下面也可以：

double exp = 1e-10;
if (Math.abs(val1 - val2)>-1*exp && Math.abs(val1 - val2)<exp) {
 //do things
}