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
}