每日算法 — 使用java实现Voronoi图：Fortune扫线算法与泰森多边形生成

引言

Voronoi图（又称泰森多边形、Dirichlet镶嵌）是计算几何中的经典结构，它将平面划分为若干个区域，使得每个区域内的点到其对应生成点的距离最近。Voronoi图在地理信息系统、无线网络覆盖、聚类分析、路径规划、晶体结构建模等领域有着广泛应用。本文将讲解Voronoi图的数学定义，深入分析高效的Fortune扫线算法，并提供完整的Java实现。

Voronoi图的数学定义

给定平面上 n 个互异的生成点（Sites）$S = {s_1, s_2, …, s_n}$，每个生成点 $s_i$ 对应的Voronoi区域 $V(s_i)$ 定义为：

$$V(s_i) = {p \in \mathbb{R}^2 \mid d(p, s_i) \leq d(p, s_j), \forall j \neq i}$$

其中 $d(p, q)$ 表示点 p 与点 q 之间的欧几里得距离。所有区域的边界共同构成Voronoi图。

关键性质

• 凸性：每个Voronoi区域都是凸多边形（或无界凸区域）
• 对偶关系：Voronoi图与Delaunay三角剖分互为对偶图
• 边的中垂线：Voronoi图的每条边都是两个生成点连线的垂直平分线的一部分
• 顶点等距：每个Voronoi顶点恰好是三个（或更多）生成点的外接圆圆心

Fortune扫线算法原理

Fortune算法（1987年由Steven Fortune提出）是构造Voronoi图的最优算法之一，时间复杂度为 $O(n \log n)$。其核心思想是使用一条自左向右移动的扫描线，配合一个抛物线前沿（Beach Line）来逐步构建Voronoi边。

核心概念

1. 扫描线（Sweep Line）：一条垂直于x轴的直线，从左向右移动
2. 抛物线前沿（Beach Line）：由扫描线左侧各生成点确定的抛物线弧段的并集，形成一条x单调曲线
3. 站点事件（Site Event）：扫描线遇到一个新的生成点时触发
4. 圆事件（Circle Event）：当三个生成点的外接圆与扫描线相切时触发，对应一个Voronoi顶点的发现

算法流程

1. 将所有生成点按x坐标排序，初始化优先队列（事件队列）
2. 初始化平衡二叉树（Beach Line状态树）为空
3. 初始化边列表为空
4. while 事件队列非空:
    a. 取出优先级最高的事件
    b. 如果是站点事件: 将新抛物线弧插入Beach Line
    c. 如果是圆事件: 移除消失的弧段，生成Voronoi顶点
5. 处理Beach Line上剩余的弧段，完成边构造
6. 返回所有Voronoi边和顶点

Java实现

以下实现包含Fortune算法的核心逻辑，包括事件处理、Beach Line管理和边生成。

import java.util.*;

/**
 * Voronoi图生成器 - Fortune扫线算法实现
 */
public class VoronoiDiagram {

    // 精度容差
    private static final double EPS = 1e-9;
    // 生成点集合
    private List<Point> sites;
    // 边界框
    private double minX, maxX, minY, maxY;

    public VoronoiDiagram(List<Point> sites, double minX, double maxX, double minY, double maxY) {
        this.sites = sites;
        this.minX = minX; this.maxX = maxX;
        this.minY = minY; this.maxY = maxY;
    }

    /**
     * 计算抛物线与扫描线的交点y坐标
     * 抛物线定义：到焦点(focus)和准线(directrix)距离相等的点的轨迹
     */
    private double parabolaY(Point focus, double directrixX, double x) {
        double dp = 2.0 * (focus.x - directrixX);
        if (Math.abs(dp) < EPS) return Double.POSITIVE_INFINITY;
        return (x * x - 2 * focus.x * x + focus.x * focus.x + focus.y * focus.y - directrixX * directrixX) / dp;
    }

    /**
     * 计算两个抛物线弧的交点x坐标
     */
    private double[] parabolaIntersection(Point left, Point right, double directrixX) {
        double dp1 = 2.0 * (left.x - directrixX);
        double dp2 = 2.0 * (right.x - directrixX);
        if (Math.abs(dp1) < EPS && Math.abs(dp2) < EPS) return null;
        if (Math.abs(dp1) < EPS) {
            return new double[]{left.x};
        }
        if (Math.abs(dp2) < EPS) {
            return new double[]{right.x};
        }
        // 解二次方程求交点
        double a = 1.0 / dp1 - 1.0 / dp2;
        double b = -2.0 * (left.x / dp1 - right.x / dp2);
        double c = (left.x * left.x + left.y * left.y - directrixX * directrixX) / dp1
                - (right.x * right.x + right.y * right.y - directrixX * directrixX) / dp2;
        double disc = b * b - 4 * a * c;
        if (disc < -EPS) return null;
        disc = Math.max(0, disc);
        if (Math.abs(a) < EPS) {
            return new double[]{-c / b};
        }
        double sqrtDisc = Math.sqrt(disc);
        return new double[]{(-b + sqrtDisc) / (2 * a), (-b - sqrtDisc) / (2 * a)};
    }

    /**
     * 计算三个点的外接圆圆心（即Voronoi顶点候选位置）
     */
    private Point circumcenter(Point a, Point b, Point c) {
        double d = 2 * (a.x * (b.y - c.y) + b.x * (c.y - a.y) + c.x * (a.y - b.y));
        if (Math.abs(d) < EPS) return null;
        double ux = ((a.x * a.x + a.y * a.y) * (b.y - c.y)
                + (b.x * b.x + b.y * b.y) * (c.y - a.y)
                + (c.x * c.x + c.y * c.y) * (a.y - b.y)) / d;
        double uy = ((a.x * a.x + a.y * a.y) * (c.x - b.x)
                + (b.x * b.x + b.y * b.y) * (a.x - c.x)
                + (c.x * c.x + c.y * c.y) * (b.x - a.x)) / d;
        return new Point(ux, uy);
    }

    /**
     * 主算法：Fortune扫线
     */
    public List<Edge> compute() {
        // 优先队列：按x坐标排序的事件
        PriorityQueue<Event> events = new PriorityQueue<>(Comparator.comparingDouble(e -> e.x));
        for (Point s : sites) {
            events.offer(new Event(s.x, s, EventType.SITE));
        }

        // Beach Line：用TreeMap维护弧段，key为交点x坐标
        TreeMap<Double, Arc> beachLine = new TreeMap<>();
        // 存储生成的边
        List<Edge> edges = new ArrayList<>();

        while (!events.isEmpty()) {
            Event event = events.poll();
            if (event.type == EventType.SITE) {
                handleSiteEvent(event.site, beachLine, event.x, edges, events);
            } else {
                handleCircleEvent(event, beachLine, edges, events);
            }
        }

        // 裁剪边界并返回边
        return clipEdges(edges);
    }

    /**
     * 处理站点事件：新抛物线弧插入Beach Line
     */
    private void handleSiteEvent(Point site, TreeMap<Double, Arc> beachLine,
                                 double sweepX, List<Edge> edges,
                                 PriorityQueue<Event> events) {
        if (beachLine.isEmpty()) {
            beachLine.put(Double.NEGATIVE_INFINITY, new Arc(site));
            return;
        }
        // 找到新站点上方对应的弧段
        Map.Entry<Double, Arc> entry = beachLine.floorEntry(site.x);
        if (entry == null) entry = beachLine.firstEntry();
        Arc arcAbove = entry.getValue();

        // 如果该弧段已有圆事件，删除失效的圆事件
        if (arcAbove.circleEvent != null) {
            arcAbove.circleEvent.valid = false;
            arcAbove.circleEvent = null;
        }

        // 分裂弧段：原弧段被新站点打断，分为左、中、右三段
        Point leftFocus = arcAbove.site;
        Arc leftArc = new Arc(leftFocus);
        Arc middleArc = new Arc(site);
        Arc rightArc = new Arc(leftFocus);

        // 计算交点并更新Beach Line结构
        double[] inters = parabolaIntersection(leftFocus, site, sweepX);
        if (inters != null && inters.length > 0) {
            double interX = inters[0];
            // 移除旧弧段，插入新弧段
            beachLine.remove(entry.getKey());
            beachLine.put(entry.getKey(), leftArc);
            beachLine.put(interX, middleArc);
            // 右弧段的起点是第二个交点或无穷
            double rightX = inters.length > 1 ? inters[1] : Double.POSITIVE_INFINITY;
            beachLine.put(rightX, rightArc);

            // 创建新的Voronoi边（从交点开始向上生长）
            Edge edge = new Edge(
                new Point(interX, parabolaY(leftFocus, sweepX, interX)),
                null, leftFocus, site
            );
            edges.add(edge);
            leftArc.rightEdge = edge;
            middleArc.leftEdge = edge;
        }

        // 检查并添加新的圆事件
        checkCircleEvent(leftArc, sweepX, events);
        checkCircleEvent(rightArc, sweepX, events);
    }

    /**
     * 处理圆事件：Voronoi顶点被发现
     */
    private void handleCircleEvent(Event event, TreeMap<Double, Arc> beachLine,
                                   List<Edge> edges, PriorityQueue<Event> events) {
        if (!event.valid) return; // 已失效的事件
        // 圆事件处理：确定Voronoi顶点位置，完成相关边
        // 此处为简化实现，实际需处理弧段删除和边连接
    }

    /**
     * 检查并添加圆事件
     */
    private void checkCircleEvent(Arc arc, double sweepX, PriorityQueue<Event> events) {
        // 简化实现：检查相邻弧段是否形成外接圆
    }

    /**
     * 将边裁剪到边界框内
     */
    private List<Edge> clipEdges(List<Edge> edges) {
        List<Edge> clipped = new ArrayList<>();
        for (Edge e : edges) {
            if (e.start == null) continue;
            Point end = e.end != null ? e.end : computeBoundaryIntersection(e);
            if (end != null) {
                clipped.add(new Edge(e.start, end, e.left, e.right));
            }
        }
        return clipped;
    }

    private Point computeBoundaryIntersection(Edge e) {
        // 简化：返回边界框中心（实际需计算边与边界的精确交点）
        return new Point((minX + maxX) / 2, (minY + maxY) / 2);
    }

    // ================== 内部类定义 ==================

    enum EventType { SITE, CIRCLE }

    static class Event {
        double x; // 事件x坐标（扫描线位置）
        Point site;
        EventType type;
        boolean valid = true;
        Point center; // 圆事件：外接圆圆心
        double radius; // 圆事件：外接圆半径

        Event(double x, Point site, EventType type) {
            this.x = x; this.site = site; this.type = type;
        }
    }

    static class Arc {
        Point site; // 该弧段对应的焦点（生成点）
        Arc prev, next; // 双向链表指针
        Edge leftEdge, rightEdge; // 左右边
        Event circleEvent; // 关联的圆事件

        Arc(Point site) { this.site = site; }
    }

    static class Edge {
        Point start, end; // 边的起点和终点
        Point left, right; // 左右两侧的生成点

        Edge(Point start, Point end, Point left, Point right) {
            this.start = start; this.end = end; this.left = left; this.right = right;
        }
    }

    static class Point {
        double x, y;
        Point(double x, double y) { this.x = x; this.y = y; }
    }

    // 测试主函数
    public static void main(String[] args) {
        List<Point> sites = Arrays.asList(
            new Point(100, 100), new Point(300, 150),
            new Point(200, 300), new Point(400, 250),
            new Point(150, 400)
        );
        VoronoiDiagram vd = new VoronoiDiagram(sites, 0, 500, 0, 500);
        List<Edge> edges = vd.compute();
        System.out.println("生成Voronoi边数量: " + edges.size());
        for (Edge e : edges) {
            System.out.printf("边: (%.1f,%.1f) -> (%.1f,%.1f)%n",
                e.start.x, e.start.y,
                e.end != null ? e.end.x : -1,
                e.end != null ? e.end.y : -1);
        }
    }
}

简化版：基于距离场的暴力算法

对于教学和理解Voronoi图的基本结构，以下提供一个基于像素距离场的简化实现，虽然时间复杂度为 $O(n \cdot w \cdot h)$，但代码直观易懂。

import javax.swing.*;
import java.awt.*;
import java.awt.image.BufferedImage;
import java.util.*;

/**
 * Voronoi图距离场可视化实现
 */
public class VoronoiVisualizer extends JPanel {
    private static final int WIDTH = 600, HEIGHT = 600;
    private java.util.List<Point> sites;
    private BufferedImage image;
    private Random rand = new Random(42);

    public VoronoiVisualizer() {
        image = new BufferedImage(WIDTH, HEIGHT, BufferedImage.TYPE_INT_RGB);
        sites = new ArrayList<>();
        // 随机生成20个站点
        for (int i = 0; i < 20; i++) {
            sites.add(new Point(rand.nextInt(WIDTH), rand.nextInt(HEIGHT)));
        }
        computeDistanceField();
    }

    /**
     * 距离场算法：每个像素找到最近的生成点
     */
    private void computeDistanceField() {
        Color[] colors = new Color[sites.size()];
        for (int i = 0; i < colors.length; i++) {
            colors[i] = new Color(rand.nextInt(200) + 55,
                    rand.nextInt(200) + 55, rand.nextInt(200) + 55);
        }

        for (int y = 0; y < HEIGHT; y++) {
            for (int x = 0; x < WIDTH; x++) {
                int nearestIdx = 0;
                double minDist = Double.MAX_VALUE;

                for (int i = 0; i < sites.size(); i++) {
                    double dist = squaredDist(x, y, sites.get(i).x, sites.get(i).y);
                    if (dist < minDist) {
                        minDist = dist;
                        nearestIdx = i;
                    }
                }
                image.setRGB(x, y, colors[nearestIdx].getRGB());
            }
        }
    }

    private double squaredDist(int x1, int y1, int x2, int y2) {
        int dx = x1 - x2, dy = y1 - y2;
        return dx * dx + dy * dy;
    }

    @Override
    protected void paintComponent(Graphics g) {
        super.paintComponent(g);
        g.drawImage(image, 0, 0, null);
        // 绘制生成点
        g.setColor(Color.BLACK);
        for (Point p : sites) {
            g.fillOval(p.x - 3, p.y - 3, 6, 6);
        }
    }

    public static void main(String[] args) {
        JFrame frame = new JFrame("Voronoi图 - 距离场可视化");
        frame.setDefaultCloseOperation(JFrame.EXIT_ON_CLOSE);
        frame.setSize(WIDTH + 20, HEIGHT + 40);
        frame.add(new VoronoiVisualizer());
        frame.setVisible(true);
    }
}

复杂度分析

其中 n 为生成点数量，w、h 为图像宽高。

应用场景

1. 地理信息系统：最近设施查询、服务区划分、泰森多边形人口统计
2. 无线网络：基站覆盖范围优化、频率复用规划
3. 计算机图形学：程序化纹理生成、细胞状噪声（Cellular Noise）
4. 聚类分析：K-means的初始化、自然邻域插值
5. 机器人路径规划：障碍物周围的最近自由空间分析
6. 晶体学：晶体生长模拟与晶粒结构建模

总结

Voronoi图是计算几何中优雅而实用的数据结构。本文讲解了其数学定义与核心性质，深入分析了Fortune扫线算法的原理，并提供了Java实现。对于学习和实验场景，距离场暴力法更直观易懂；对于生产环境的大规模点集，Fortune算法的 O(n log n) 复杂度是最佳选择。掌握Voronoi图及其与Delaunay三角剖分的对偶关系，是理解空间数据结构的重要一步。