平面点集三角化的示例

平面点集三角化的代码，虽然做了一些优化，但是效率还是非常低，在笔记本上跑，10K个点1秒左右，好处是代码好理解。源码如下：

 

package ImageProcess;

import java.awt.Color;
import java.awt.Graphics;
import java.awt.geom.Point2D;
import java.awt.geom.Point2D.Double;
import java.awt.image.BufferedImage;
import java.util.Arrays;
import java.util.Calendar;
import java.util.Comparator;

import javax.swing.JFrame;

public class Delaunay {
    public static long counter = 0;

    // 三角形的边：
    class IndexedEdge {
        public int p1;
        public int p2;
    }

    // 三角形：
    class IndexedTriangle {
        public int p1;
        public int p2;
        public int p3;
    }

    class Circle {
        public double xc;
        public double yc;
        public double rad;
    }

    // 顶点坐标：
    Point2D.Double[] points;
    int nPoint;

    // 三角形：
    IndexedTriangle[] triangles;
    int nTriangle;

    private static double EPS = 1.192092896e-07F;

    public static void main(String[] args) {
        window w = new window();
    }

    // 初始化：
    @SuppressWarnings("unchecked")
    public Delaunay(Point2D.Double[] verts) {
        int i;

        // 顶点数：
        nPoint = verts.length;

        // 总的三角数：
        nTriangle = 0;

        // 顶点数组：
        points = new Point2D.Double[nPoint + 3];
        for (i = 0; i < points.length; i++) {
            points[i] = new Point2D.Double();
            if (i < nPoint) {
                points[i].x = verts[i].x;
                points[i].y = verts[i].y;
            }
        }

        // 三角形数组：
        triangles = new IndexedTriangle[nPoint * 3];
        for (i = 0; i < triangles.length; i++) {
            triangles[i] = new IndexedTriangle();
        }

        // 顶点数组按横坐标排序：
        Arrays.sort(points, new Comparator<Point2D.Double>() {
            @Override
            public int compare(Double o1, Double o2) {
                return o1.x > o2.x ? 1 : -1;
            }
        });
    }

    /*
     * 检查点(px,py)是否在点(x1,y1),(x2,y2),(x3,y3)确定的圆内， 
     * 注意 : (px,py)点与其中一边张开的角小于该边的对角，就可以证明该点在三角形的外接圆的外面
     */
    public final boolean inCircle1(double px, double py, double x1, double y1, double x2, double y2, double x3,
            double y3) {

        double aa = (x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2);
        double bb = (x1 - x3) * (x1 - x3) + (y1 - y3) * (y1 - y3);
        double cc = (x2 - x3) * (x2 - x3) + (y2 - y3) * (y2 - y3);
        double AA = (px - x2) * (px - x2) + (py - y2) * (py - y2);
        double BB = (px - x3) * (px - x3) + (py - y3) * (py - y3);

        double a = Math.sqrt(aa);
        double b = Math.sqrt(bb);
        double A = Math.sqrt(AA);
        double B = Math.sqrt(BB);

        // 根据余弦定理计算两角的余弦值：
        double s = (aa + bb - cc) * (A * B);
        double r = (AA + BB - cc) * (a * b);

        // 如果不在同一边，则做一转换：
        double d = ((y1 - y3) * (x2 - x3) - (x1 - x3) * (y2 - y3)) * ((py - y3) * (x2 - x3) - (px - x3) * (y2 - y3));
        if (d < 0) {
            s *= -1;
        }

        // 如果点P对边的角度比该边的对角值大，则说明点在外接圆内：
        return s > r;
    }
    
    /**
     * 判断点(px,py)是否包含在(ax,ay),(bx,by),(cx,cy)三点确定圆内
     * @param px    p点的坐标
     * @param py
     * @param ax    三个点的坐标
     * @param ay
     * @param bx
     * @param by
     * @param cx
     * @param cy
     * @return
     */
    private boolean inCircle(double px, double py, double ax, double ay, double bx, double by, double cx, double cy) {
        double dx = ax - px;
        double dy = ay - py;
        double ex = bx - px;
        double ey = by - py;
        double fx = cx - px;
        double fy = cy - py;

        double ap = dx * dx + dy * dy;
        double bp = ex * ex + ey * ey;
        double cp = fx * fx + fy * fy;

        return dx * (ey * cp - bp * fy) -
               dy * (ex * cp - bp * fx) +
               ap * (ex * fy - ey * fx) < 0;
    }    

    /*
     * 检查p点是否在p1,p2,p3确定的圆内，并返回三点确定的圆的圆心和半径： 注意 : p点若在三角形的一条边上，则可确定其在圆内：
     */
    public boolean circumCircle(Point2D.Double p, Point2D.Double p1, Point2D.Double p2, Point2D.Double p3, Circle cir) {
        double m1, m2, mx1, mx2, my1, my2;
        double dx, dy, rsqr, drsqr;

        /* 检查是否共线 */
        if (Math.abs(p1.y - p2.y) < EPS && Math.abs(p2.y - p3.y) < EPS) {
            return false;
        }

        if (Math.abs(p2.y - p1.y) < EPS) {
            m2 = -(p3.x - p2.x) / (p3.y - p2.y);
            mx2 = (p2.x + p3.x) / 2.0;
            my2 = (p2.y + p3.y) / 2.0;
            cir.xc = (p2.x + p1.x) / 2.0;
            cir.yc = m2 * (cir.xc - mx2) + my2;
        } else if (Math.abs(p3.y - p2.y) < EPS) {
            m1 = -(p2.x - p1.x) / (p2.y - p1.y);
            mx1 = (p1.x + p2.x) / 2.0;
            my1 = (p1.y + p2.y) / 2.0;
            cir.xc = (p3.x + p2.x) / 2.0;
            cir.yc = m1 * (cir.xc - mx1) + my1;
        } else {
            m1 = -(p2.x - p1.x) / (p2.y - p1.y);
            m2 = -(p3.x - p2.x) / (p3.y - p2.y);
            mx1 = (p1.x + p2.x) / 2.0;
            mx2 = (p2.x + p3.x) / 2.0;
            my1 = (p1.y + p2.y) / 2.0;
            my2 = (p2.y + p3.y) / 2.0;
            cir.xc = (m1 * mx1 - m2 * mx2 + my2 - my1) / (m1 - m2);
            cir.yc = m1 * (cir.xc - mx1) + my1;
        }
        dx = p2.x - cir.xc;
        dy = p2.y - cir.yc;
        rsqr = dx * dx + dy * dy;
        cir.rad = Math.sqrt(rsqr);
        dx = p.x - cir.xc;
        dy = p.y - cir.yc;
        drsqr = dx * dx + dy * dy;
        return drsqr <= rsqr;
    }

    /**
     * 检查当前点的横坐标px与三个点确定的圆心的距离是否已大于圆的半径：
     * @param px    当前点的横坐标
     * @param p1x    三个顶点的坐标
     * @param p1y
     * @param p2x
     * @param p2y
     * @param p3x
     * @param p3y
     * @return
     */
    public final boolean isComplete(double px, double p1x, double p1y, double p2x, double p2y, double p3x, double p3y) {
        //counter++;
        double m1, m2, mx1, mx2, my1, my2;
        double dx, dy, rsqr;
        double xc, yc, rad;

        /* 检查是否共线 */
        if (Math.abs(p1y - p2y) < EPS && Math.abs(p2y - p3y) < EPS) {
            return false;
        }

        if (Math.abs(p2y - p1y) < EPS) {
            m2 = -(p3x - p2x) / (p3y - p2y);
            mx2 = (p2x + p3x) / 2.0;
            my2 = (p2y + p3y) / 2.0;
            xc = (p2x + p1x) / 2.0;
            yc = m2 * (xc - mx2) + my2;
        } else if (Math.abs(p3y - p2y) < EPS) {
            m1 = -(p2x - p1x) / (p2y - p1y);
            mx1 = (p1x + p2x) / 2.0;
            my1 = (p1y + p2y) / 2.0;
            xc = (p3x + p2x) / 2.0;
            yc = m1 * (xc - mx1) + my1;
        } else {
            m1 = -(p2x - p1x) / (p2y - p1y);
            m2 = -(p3x - p2x) / (p3y - p2y);
            mx1 = (p1x + p2x) / 2.0;
            mx2 = (p2x + p3x) / 2.0;
            my1 = (p1y + p2y) / 2.0;
            my2 = (p2y + p3y) / 2.0;
            xc = (m1 * mx1 - m2 * mx2 + my2 - my1) / (m1 - m2);
            yc = m1 * (xc - mx1) + my1;
        }
        dx = p2x - xc;
        dy = p2y - yc;
        rsqr = dx * dx + dy * dy;
        rad = Math.sqrt(rsqr);
        return px > rad + xc + EPS;
    }

    int Triangulate() {
        boolean[] complete;
        IndexedEdge[] edges;
        int nedge = 0; // 边的个数：
        int trimax; // 三角最大数
        int emax; // 边的最大数
        int i, j, k;
        double xp, yp, x1, y1, x2, y2, x3, y3, xc, yc, r;
        double xmin, xmax, ymin, ymax, xmid, ymid;
        double dx, dy, dmax;

        /* 为各变量分配内存： */
        trimax = 2 * (nPoint + 3) - 5;
        complete = new boolean[trimax];
        emax = 3 * (nPoint + 3) - 6;
        edges = new IndexedEdge[emax];
        for (i = 0; i < edges.length; i++) {
            edges[i] = new IndexedEdge();
        }
        /*
         * 寻找顶点的范围并构造超级三角形：
         */
        xmin = points[0].x;
        ymin = points[0].y;
        xmax = xmin;
        ymax = ymin;
        for (i = 1; i < nPoint; i++) {
            if (points[i].x < xmin) {
                xmin = points[i].x;
            }
            if (points[i].x > xmax) {
                xmax = points[i].x;
            }
            if (points[i].y < ymin) {
                ymin = points[i].y;
            }
            if (points[i].y > ymax) {
                ymax = points[i].y;
            }
        }
        dx = xmax - xmin;
        dy = ymax - ymin;
        dmax = (dx > dy) ? dx : dy;
        xmid = (xmax + xmin) / 2.0;
        ymid = (ymax + ymin) / 2.0;
        /*
         * 构造超级三角形，将该三角形的顶点添加到顶点数组的末尾， 并将该三角形放在三角形数组的起始位置：
         */
        points[nPoint + 0].x = xmid - 20 * dmax;
        points[nPoint + 0].y = ymid - dmax;
        points[nPoint + 1].x = xmid;
        points[nPoint + 1].y = ymid + 20 * dmax;
        points[nPoint + 2].x = xmid + 20 * dmax;
        points[nPoint + 2].y = ymid - dmax;

        //顺时针方向添加超级三角形：
        triangles[0].p1 = nPoint;
        triangles[0].p2 = nPoint + 1;
        triangles[0].p3 = nPoint + 2;
        complete[0] = false;
        nTriangle = 1;
        
        /*
         * 将所有点依次加入现有网格中：
         */
        for (i = 0; i < nPoint; i++) {
            // 取一个点：
            xp = points[i].x;
            yp = points[i].y;
            nedge = 0;

            /*
             * 建立边的缓存。如果有一点P落在三角形的外接圆内，那么将三角形的三条边放到缓存中，并移除这个三角形：
             */
            for (j = 0; j < nTriangle; j++) {
                if (complete[j]) {
                    continue;
                }

                IndexedTriangle curTri = triangles[j];
                x1 = points[curTri.p1].x;
                y1 = points[curTri.p1].y;
                x2 = points[curTri.p2].x;
                y2 = points[curTri.p2].y;
                x3 = points[curTri.p3].x;
                y3 = points[curTri.p3].y;

                // 如果点(xp,yp)在三个点确定的圆内，则重构三角形：
                if (inCircle(xp, yp, x1, y1, x2, y2, x3, y3)) {
                    // 缓存3条边：
                    edges[nedge + 0].p1 = curTri.p1;
                    edges[nedge + 0].p2 = curTri.p2;
                    edges[nedge + 1].p1 = curTri.p2;
                    edges[nedge + 1].p2 = curTri.p3;
                    edges[nedge + 2].p1 = curTri.p3;
                    edges[nedge + 2].p2 = curTri.p1;
                    nedge += 3;

                    // 还原三个点完成标记：
                    curTri.p1 = triangles[nTriangle - 1].p1;
                    curTri.p2 = triangles[nTriangle - 1].p2;
                    curTri.p3 = triangles[nTriangle - 1].p3;
                    complete[j] = false;

                    // 减去总数：
                    nTriangle--;
                    j--;
                } else {
                    //如果当前点到三角形的外接圆心的距离已经开始大于半径了，则后续点则无需比较：
                    complete[j] = this.isComplete(xp, x1, y1, x2, y2, x3, y3);
                }
            }
            
            /*
             * 标记公共边。如果两个共边的三角形的外接圆均包含该点，则此公共边在重构三角形时应该去掉，在这里先做上标记：
             */
            for (j = 0; j < nedge - 1; j++) {
                for (k = j + 1; k < nedge; k++) {
                    if ((edges[j].p1 == edges[k].p2) && (edges[j].p2 == edges[k].p1)) {
                        edges[j].p1 = -1;
                        edges[j].p2 = -1;
                        edges[k].p1 = -1;
                        edges[k].p2 = -1;
                    }
                }
            }
            
            /*
             * 跳过已经做了标记的边，构建三角形。所有的边按顺时针方向排列：
             */
            for (j = 0; j < nedge; j++) {
                if (edges[j].p1 < 0 || edges[j].p2 < 0)
                    continue;

                triangles[nTriangle].p1 = edges[j].p1;
                triangles[nTriangle].p2 = edges[j].p2;
                triangles[nTriangle].p3 = i;
                complete[nTriangle] = false;
                nTriangle++;
            }
        }
        /*
         * 移除和超级三角形三个顶点相连的三角形的边：
         */
        for (i = 0; i < nTriangle; i++) {
            if (triangles[i].p1 >= nPoint || triangles[i].p2 >= nPoint || triangles[i].p3 >= nPoint) {
                triangles[i].p1 = triangles[nTriangle - 1].p1;
                triangles[i].p2 = triangles[nTriangle - 1].p2;
                triangles[i].p3 = triangles[nTriangle - 1].p3;
                nTriangle--;
                i--;
            }
        }

        return 0;
    }

    public void showPoints(Graphics g) {
        int i, x, y;

        g.translate(20, 50);

        g.setColor(Color.RED);
        for (i = 0; i < nPoint; i++) {
            x = (int) Math.round(points[i].x);
            y = (int) Math.round(points[i].y);

            g.drawOval(x - 2, y - 2, 4, 4);
        }
        g.setColor(Color.black);
    }

    public void ShowTriangles(Graphics g) {
        int i, x1, y1, x2, y2, x3, y3;

        for (i = 0; i < nTriangle; i++) {
            x1 = (int) Math.round(points[triangles[i].p1].x);
            y1 = (int) Math.round(points[triangles[i].p1].y);
            x2 = (int) Math.round(points[triangles[i].p2].x);
            y2 = (int) Math.round(points[triangles[i].p2].y);
            x3 = (int) Math.round(points[triangles[i].p3].x);
            y3 = (int) Math.round(points[triangles[i].p3].y);

            g.drawLine(x1, y1, x2, y2);
            g.drawLine(x2, y2, x3, y3);
            g.drawLine(x3, y3, x1, y1);
        }
    }
}

/**用于显示的类*/

class window extends JFrame {
    public window() {
        this.setBounds(0, 0, 1200, 900);

        this.setDefaultCloseOperation(JFrame.EXIT_ON_CLOSE);

        this.setVisible(true);
    }

    public void paint(Graphics g) {
        g.clearRect(0, 0, this.getWidth(), this.getHeight());

        BufferedImage bi = new BufferedImage(2000, 1000, BufferedImage.TYPE_3BYTE_BGR);
        Graphics gc = bi.getGraphics();
        gc.setColor(Color.white);
        gc.fillRect(0, 0, 2000, 1000);

        Point2D.Double[] lst = new Point2D.Double[10000];
        for (int i = 0; i < lst.length; i++) {
            lst[i] = new Point2D.Double(Math.random() * 800, Math.random() * 800);
        }

        Calendar t1 = Calendar.getInstance();
        Delaunay d = new Delaunay(lst);

        d.Triangulate();

        Calendar t2 = Calendar.getInstance();

        System.out.println("耗时：" + (t2.getTimeInMillis() - t1.getTimeInMillis()) + "ms");

        d.showPoints(gc);

        d.ShowTriangles(gc);

        g.drawImage(bi, 0, 0, null);

        System.out.println("size=" + d.nTriangle);

        System.out.println("判断" + Delaunay.counter + "次");
    }
}