凸包を求める

前回のようなランダムな点があったとき、その点を全部含んでへこみのない多角形を作成するという問題。


アルゴリズムとしては、「他の点でできる三角形に含まれてたら凸包上の点じゃない」という考え方で、全部の組み合わせの三角形について全部の点を走査します。最適なアルゴリズムではないですが、とりあえず今回はこれで。


で、点Pが三角形ABCに含まれるかどうかは、三角形ABCが時計周りだとしたら、三角形ABP、BCP、CAPも時計回りだということを使います。


点Pが三角形ABCの外にある場合は、三角形ABP、BCP、CAPのどれかが三角形ABCと逆周りになります。


三角形の向きは符号付面積を使うのですが、三角形の座標をそれぞれ(x1, y1)(x2, y2)(x3, y3)とすると(x1 - x2)(y1-y3)-(x1-x3)(y1-y2)で得られます。


そうやって凸包上に乗らない点を削除したら、隣り合う2点と重心とでできる三角形が時計周りになるように、それぞれの点を並べなおします。これにも符号付三角形を使います。

public class ConvexHull extends JComponent{
    static List<Point> points = new ArrayList<Point>();
    static List<Point> hull;
    
    //凸包を求める
    static void createConvexHull(){
        // 凸包上の点を得る
        hull = (ArrayList<Point>)(((ArrayList<Point>)points).clone());
        for(int i = 0; i < hull.size() - 2; ++i){
            for(int j = i + 1; j < hull.size() - 1; ++j){
                for(int k = j + 1; k < hull.size(); ++k){
                    for(int l = 0; l < hull.size(); ++l){
                        if(i == l || j == l || k == l) continue;
                        //点lが三角形ijkに含まれていたら、凸包の点ではない
                        if(isContained(
                                hull.get(i), hull.get(j), hull.get(k), hull.get(l)))
                        {
                            hull.remove(l);
                            if(i > l) --i;
                            if(j > l) --j;
                            if(k > l) --k;
                            --l;
                        }
                    }
                }
            }
        }
        
        //重心を求める
        int tempx = 0, tempy = 0;
        for(Point p : hull){
            tempx += p.x;
            tempy += p.y;
        }
        final Point centroid = new Point(tempx / hull.size(), tempy / hull.size());
        
        //重心からの角度によって並べ替え
        Collections.sort(hull, new Comparator() {
            public int compare(Object o1, Object o2) {
                Point p1 = (Point) o1;
                Point p2 = (Point) o2;
                int xx1 = p1.x - centroid.x;
                int xx2 = p2.x - centroid.x;
                if(xx1 * xx2 < 0) return xx2 - xx1;
                return getSignedArea(p1, p2, centroid);
            }
        });
    }

    public static void main(String[] args){
        final int n = 25;//点の個数
        final int width = 220;//領域の大きさ
        //ランダムな座標を得る
        final int r = width / n * 2 - 1;
        int x = 0;
        int y = 0;
        List<Integer> xs = new ArrayList<Integer>();
        List<Integer> ys = new ArrayList<Integer>();
        for(int i = 0; i < n; ++i){
            x+= Math.random() * r + 1;
            xs.add(x);
            y += Math.random() * r + 1;
            ys.add(y);
        }
        
        Collections.shuffle(xs);
        Collections.shuffle(ys);
        for(int i = 0; i < n; ++i){
            points.add(new Point(xs.get(i), ys.get(i)));
        }
        
        createConvexHull();
        
        JFrame f = new JFrame("凸包");
        f.setDefaultCloseOperation(JFrame.EXIT_ON_CLOSE);
        f.add(new ConvexHull());
        f.setSize(400, 300);
        f.setVisible(true);
    }
    
    @Override
    protected void paintComponent(Graphics g) {
        g.setColor(Color.WHITE);
        g.fillRect(0, 0, 400, 350);
        g.setColor(Color.BLACK);
        for(Point p : points){
            g.drawOval(50 + p.x, 10 + p.y, 3, 3);
        }
        Polygon pl = new Polygon();
        for(Point p : hull){
            g.fillOval(50 + p.x, 10 + p.y, 3, 3);
            pl.addPoint(50 + p.x, 10 + p.y);
        }
        g.drawPolygon(pl);
    }

    /** 符号付面積を求める */
    static int getSignedArea(Point p1, Point p2, Point p3){
        return (p2.x - p1.x) * (p3.y - p1.y) 
                - (p3.x - p1.x) * (p2.y - p1.y);
    }

    /** 点の三角形の内外判定 */
    static boolean isContained(Point p1, Point p2, Point p3, Point p){
        int a = getSignedArea(p1, p2, p3);
        if(a == 0) return false;//3点が直線に並んでる
        int b = getSignedArea(p1, p2, p);
        if(a * b <= 0) return false;//aとbの符号が違うということは3角形の外
        b = getSignedArea(p2, p3, p);
        if(a * b <= 0) return false;
        b = getSignedArea(p3, p1, p);
        if(a * b <= 0) return false;
        return true;
    }
}