65.9K
CodeProject 正在变化。 阅读更多。
Home

创建 Voronoi 图 1/2

Kenneth Haugland
starIconstarIconstarIconstarIcon
emptyStarIcon
starIcon

4.93/5 (31投票s)

2012 年 6 月 29 日

CPOL

12分钟阅读

viewsIcon

120317

downloadIcon

2205

创建 Voronoi 图的详细描述，基本功能回顾

引言

Voronoi 图是计算几何中非常有用的工具，用途广泛，例如计算森林中每棵树的面积，或根据受害者地址推断城市中被污染的水井位置等等。总的来说，它有助于找出“谁离谁最近”。下面列出了一些使用 Voronoi 图的问题：

  • 碰撞检测
  • 模式识别
  • 地理优化
  • 几何聚类
  • 最近点对算法
  • k-最近邻查询

如果您对更多示例感兴趣，可以访问 此处 网站。

Voronoi 图的历史相当复杂，因为它被多次重新发现。Voronoi 图的第一个书面描述是由笛卡尔在 17 世纪中叶完成的。Dirichlet 在 1850 年首次构思了这个通用概念，但直到 Voronoi 的 1908 年文章才对其进行了严格的数学处理。有关其完整历史，您应该阅读 Aurenhammer 的文章。重要的一点是，Delaunay (1934) 证明了 Voronoi 图有一个姐妹（或者说是兄弟），称为 Delaunay 三角剖分，它（通常）可以从 Voronoi 图构造，反之亦然。（有趣的事实：Delaunay 的博士论文是在 Voronoi 的指导下完成的。） 

在最初发现后的几十年里，出现了几种构造 Voronoi 图的方法，但大多数方法通常相当慢 O(n^2)。有一种算法可以在 O(n*log n) 时间内完成，但实现起来非常困难，因此很少使用。直到 1985 年 Fortune 发现了扫描线算法（发表于 1987 年）——它也运行在 O(n*log n) 时间内——一种被普遍接受的构造 Voronoi 图的方法才得以实现。在本文中，我将尝试解释该算法的基本功能。然而，下面简要概述了构建 Voronoi 图的不同算法：

  • 扫描线
  • 分治法
  • 螺旋搜索
  • 三维凸包
  • 增量法

撰写本文介绍 Fortune 扫描线算法的原因是澄清 Voronoi 图的实现。这主要是为了在我尝试实现该算法时理清思路，并希望对他人也能有所帮助，从而提高您对基本原理的理解。是的，我熟悉 CodeProject 上 Ben 的文章。本文使用的一些代码就来自他的实现。

背景

Voronoi 图究竟是什么？嗯，有几种不同类型的 Voronoi 图，但它们都具有一些共同的、风格非常普遍的属性，可以通过简单的比喻来解释。它通常被描述为围绕一个点的区域，其中中心点是周围多边形中最近的点，如下所示：

存在一个通用的潜在模式，可以通过查找两点之间的中线来例证。我们通过构造两个半径大于两点之间距离一半的圆来找到与两点等距的线。两点 A 和 B 之间的中线图示： 

如果我们只对找出 A 和 B 两点之间连线的中点感兴趣，还有一种计算该点的方法。如果您看一下下面将绳子系在 A 和 B 两点的图示，您会立即看到中点是绳子最低的点。因此，这是一种计算两点之间中距离的补充方法，就像上面两个圆的图示一样。

如果我们用毯子代替绳子，并竖起同样高度的木棍，您可以想象，所有山谷都会给出图中的所有 Voronoi 线，前提是木棍足够高（或毯子的张力足够高），这样毯子就不会接触地面。

希望上述图像足够简单，每个人都能理解，但从这个类比中构造 Voronoi 图并不方便。用 Steven Fortune 的话说：“要实际实现一个抽象描述的几何算法，极其困难。”

一种更数学上严谨的看待 Voronoi 图的方法是放置圆锥体（或者说尖帽子——见下图），侧面的斜率为 45 度，以中心点为最高点。这将产生相同的重叠侧面，从而在 Voronoi 图的点之间产生直线，如第一张图所示。

这实际上意味着您总是在最近点之间的中点处找到重叠线。圆锥体是由许多不同半径的圆组成的，因此与中线的类比仍然有效。上面示例图的结果给出了下面的交线：

Fortune 的扫描线算法基于这些原理，但不是将圆锥体顶端朝上放置，而是将圆锥体倒置。圆锥体的侧面倾斜 45 度，而交替扫描线（下图中的平面称为 pi）与倾斜的圆锥体呈相反的 45 度角。我希望您也会注意到，通过将圆锥体倒置，就没有高度限制。所有角度都相对于 z = 0 平面给出。这种方法如下图所示，扫描线不断地从 0 移动到 x = 无穷大（或者更确切地说，是能够完成我们感兴趣区域的 Voronoi 图中所有线的边界矩形）。

对于熟悉不同几何曲线的人来说（请参阅维基百科关于圆锥曲线的文章），扫描平面 pi 和圆锥体之间的交线会形成抛物线方程，这是很清楚的。这条交线取决于圆锥体中心 P(x,y) 的相对位置以及扫描线 L 的位置。从现在开始，我们将把由不同弧生成的弯曲交线称为“海滩线”。

该算法必须始终维护扫描线的当前位置（y 坐标）。它按从左到右（或从上到下，因为您可以在 X 和 Y 方向上进行扫描）的顺序存储定义海滩线的点集（圆锥体的中心点）。然而，重要的是要知道，算法永远不需要存储海滩线的所有抛物线弧。我们将只维护扫描线状态，并对“发生重大事件”的离散事件点感兴趣。也就是说，任何改变 Voronoi 图和海滩线拓扑结构的事件。

有两种事件可能改变海滩线。第一种是点事件，第二种是圆事件。

点事件  

当扫描线扫过一个新点时，一个新的弧将被插入到海滩线中。

对于由 Fortune 扫描线算法给出的任何抛物线，其方程如下：Pjx 和 Pjy 是扫描线上圆锥体中心的 x 和 y 位置，Ly 是扫描线的位置，X 是扫描线上新交点的位置。

必须对扫描线经过的所有点使用该公式，才能找到离该点到扫描线距离最短的点。现在，您可以通过计算扫描线下落时所有点的 beta j 来轻松生成所有弧。然而，当您尝试计算 Voronoi 图时，第一个交点并不那么有趣。有趣的是找到上图中左边的两个点，海滩线在那里形成两个交点，它们将构成 Voronoi 图中新线的起点。

因此，我们假设有两个点形成两条抛物线，有两个交点，并且每条单独的抛物线都受上述抛物线方程的控制。我们现在可以令两个抛物线方程相等，并通过使用二次公式来计算两个交点发生的位置。结果是两个不同的 x 值，您可以将它们代入任一抛物线方程以获得结果的 y 值，从而得到两个交点。方程如下，其中 Pi 和 Pj 是圆锥体的两个中心点：

我不想简化这个方程，因为它相当长，但想法是将其写成如下形式的方程：

通过以下方程找到 X 的两个解：

但方程中有一个相当大的问题，因为点事件，您什么时候应该猜测 Voronoi 图中的线条？如果您在计算前稍微增加一点，线条就会产生很大的误差；如果您等待太久，就会变得一团糟，因为您必须缩短两点之间的距离，因为新的交点会破坏计算出的线条。

您可能会问，对于给定的点，您需要检查抛物线多久？答案是，只要中心点不被线包围，即没有更下游的事件会改变围绕该点的任何线的边界，您就必须继续检查该点的抛物线。这通常不是问题，因为检查是由两个条件驱动的，并且当发生圆事件（也称为顶点事件）时，抛物线弧会终止。

圆事件  

仅当两条抛物线“吞没”第三条抛物线时，才会发生此事件，这意味着我们只需要检查三个（或更多）点是否会形成一个不包含其他点的圆。 

您可能会问：“哪个弧会消失？”答案是，如果您有三个点形成像上图那样的圆，那么中间的弧会消失。

这里有一个计算圆的中心和半径的公式：在此处。如果您想知道另一个点是否在圆内，您只需计算圆心和潜在点之间的距离。如果距离小于半径，则该点在圆内。否则，圆心将为您提供另一个顶点点，并且一条抛物线将消失。计算是通过使用圆的一般方程完成的：

(X - Xc)*2 + (Y - Yc)*2 = R*2

并将三个点代入 X 和 Y，您将得到三个未知变量的三个方程，然后可以解出 Xc、Yc 和 R。还有其他（更快的）计算中心点的方法，但这部分留给您自行研究。 

代码中使用的基本函数

请注意，此代码使用了 .NET 的原生库，对于任意点集可能有一些限制。这样做是为了不影响代码的速度。 

对于任意点集，我们需要做的第一件事是按照扫描线方向对它们进行排序。一种通用的算法是以词典顺序对点进行排序（无论是在 x 方向还是 y 方向），这意味着如果您按 y 值排序，并得到相等的 y 值，则在 y 检查后，您将按 x 方向对这两个点进行排序。代码如下： 

''' <summary>
''' Returns a sorted pointcollection based on Lexografically sorting
''' </summary>
''' <param name="samplepoints"></param>
''' <param name="SortByXdirection"></param>
''' <returns></returns>
''' <remarks></remarks>
Private Function SortPoints(ByVal samplepoints As PointCollection, _
        ByVal SortByXdirection As Boolean) As PointCollection
    'Create another array so we can keep the original array out of order
    Dim copySamplePoints As Point() = New Point(samplepoints.Count - 1) {}
    samplepoints.CopyTo(copySamplePoints, 0)

    If SortByXdirection Then
        Array.Sort(copySamplePoints, New PointSort(PointSort.Mode.X))
    Else
        Array.Sort(copySamplePoints, New PointSort(PointSort.Mode.Y))
    End If

    Dim result As New PointCollection

    For Each p As Point In copySamplePoints
        result.Add(p)
    Next

    Return result
End Function

Private Class PointSort
    Implements IComparer
    Public Enum Mode
        X
        Y
    End Enum

    Private currentMode As Mode = Mode.X

    Public Sub New(ByVal mode As Mode)
        currentMode = mode
    End Sub

    'Comparing function
    'Returns one of three values - 0 (equal), 1 (greater than), -1 (less than)
    Private Function IComparer_Compare(ByVal a As Object, ByVal b As Object) _
                     As Integer Implements IComparer.Compare
        Dim point1 As Point = CType(a, Point)
        Dim point2 As Point = CType(b, Point)

        If currentMode = Mode.X Then
            'Compare X values
            If point1.X > point2.X Then
                Return 1
            ElseIf point1.X < point2.X Then
                Return -1
            Else
                If point1.Y > point2.Y Then
                    Return 1
                ElseIf point1.Y < point2.Y Then
                    Return -1
                Else 'Identical points
                    Return 0
                End If
            End If
        Else
            If point1.Y > point2.Y Then
                'Compare Y Values
                Return 1
            ElseIf point1.Y < point2.Y Then
                Return -1
            Else 'Y1 = Y2
                If point1.X > point2.X Then
                    'Compare Y Values
                    Return 1
                ElseIf point1.X < point2.X Then
                    Return -1
                Else 'Identical points
                    Return 0
                End If
            End If
        End If
    End Function
End Class

用于查找新点事件的交替抛物线函数： 

''' <summary>
''' Returns a polyline between the newpoint and the contructed arc of the oldpoint
''' </summary>
''' <param name="NewPoint">New occuring point on the sweep line</param>
''' <param name="OldPoint">Point (above the sweep line) that generates the prarabolic arc</param>
''' <param name="SweepLine">Position of the sweep line (Y) </param>
''' <returns>A straight line between the NewPoint to the Intersection point on the parabolic arc</returns>
''' <remarks></remarks>
Private Function Intersection(ByVal NewPoint As Point, ByVal OldPoint As Point, _
        ByVal SweepLine As Double) As Polyline
    Dim result As New Polyline
    result.Points.Add(NewPoint)

    Dim endpoint As New Point
    endpoint = NewPoint
    Dim res As Double
    res = (1 / (2 * ((OldPoint.Y) - SweepLine)))
    res *= (NewPoint.X ^ 2 - 2 * OldPoint.X * NewPoint.X + OldPoint.X ^ 2 + OldPoint.Y ^ 2 - SweepLine ^ 2)
    endpoint.Y = res
    result.Points.Add(endpoint)
    result.Stroke = Brushes.Red
    result.StrokeThickness = 2
    Return result
End Function

下一个函数用于计算与上一个函数相同的值，但它只返回旧点创建的弧中的 y 值。 

''' <summary>
''' Returns the y value construcked by the event from the newpoint and the contructed arc of the oldpoint
''' </summary>
''' <param name="NewPoint">This is the point on the sweep line</param>
''' <param name="OldPoint">This is one point above the sweep line</param>
''' <param name="SweepLine">Y position of the sweep line</param>
''' <returns>Calculates the distance fromn a point (NewPoint) to the Parabolic
''' arc generated by the OldPoint and the Sweep line</returns>
''' <remarks>The Function only gives you the Y value of the position
''' on the parabolic intersection given the X location</remarks>
Private Function Intersection2(ByVal NewPoint As Point, _
        ByVal OldPoint As Point, ByVal SweepLine As Double) As Double
    Dim res As Double
    res = (1 / (2 * ((OldPoint.Y) - SweepLine)))
    res *= (NewPoint.X ^ 2 - 2 * OldPoint.X * NewPoint.X + OldPoint.X ^ 2 + OldPoint.Y ^ 2 - SweepLine ^ 2)
    Return (res)
End Function

我们需要处理的下一个函数是所谓的点事件，或者我称之为抛物线切割（下面的代码中将扫描线设置在 ys + 500 的位置）。

''' <summary>
''' Calculates the line between two Parabolic intersections
''' </summary>
''' <param name="Point1">A point in the Voronoi diagram</param>
''' <param name="point2">A Point in the Voronoi diagram (different from point 1)</param>
''' <param name="ys">The position of the sweep line</param>
''' <returns>A staight line between the two intersection poitns</returns>
''' <remarks>It would give wrong values if the two points have the same or
''' close to the same Y coordinate, as the double has limited significant storage</remarks>
Private Function ParabolicCut(ByVal Point1 As Point, ByVal point2 As Point, ys As Double) As Polyline

    'Stores the values in double format, as I didnt bother to rewrite Benjamin Dittes quadratic equation.
    Dim x1, y1, x2, y2 As Double

    'Inizialize Point 1
    x1 = Point1.X
    y1 = Point1.Y

    'Inizialize Point 2 
    x2 = point2.X
    y2 = point2.Y

    'Sweep line, added 500 to make sure the two calculated points are well of the boundaries
    ys = ys + 500

    'Setting ut calculation of the quadratic equation
    Dim a1 As Double = 1 / (2 * (y1 - ys))
    Dim a2 As Double = 1 / (2 * (y2 - ys))

    'The two resulting values of x from the quadratic equation
    Dim xs1 As Double = 0.5 / (2 * a1 - 2 * a2) * (4 * a1 * x1 - 4 * a2 * x2 + 2 * _
        Math.Sqrt(-8 * a1 * x1 * a2 * x2 - 2 * a1 * y1 + 2 * a1 * y2 + 4 * a1 * a2 _
        * x2 * x2 + 2 * a2 * y1 + 4 * a2 * a1 * x1 * x1 - 2 * a2 * y2))
    Dim xs2 As Double = 0.5 / (2 * a1 - 2 * a2) * (4 * a1 * x1 - 4 * a2 * x2 - 2 * _
        Math.Sqrt(-8 * a1 * x1 * a2 * x2 - 2 * a1 * y1 + 2 * a1 * y2 + 4 * a1 * _
        a2 * x2 * x2 + 2 * a2 * y1 + 4 * a2 * a1 * x1 * x1 - 2 * a2 * y2))

    'Generate two points to store the two intersection points
    Dim p1, p2 As New Point
    p1.X = xs1
    p2.X = xs2
    p1.Y = 0
    p2.Y = 0

    'Find the resulting Y values calculated from the quadratic equation
    ' (It dosent matter that the Y is 0 as the function Intersection2
    ' only uses the X value for computation)
    p1.Y = Intersection2(p1, Point1, ys)
    p2.Y = Intersection2(p2, Point1, ys)

    'Make a polyline to store the result
    Dim result As New Polyline

    'Add the two calculated poitns
    result.Points.Add(p1)
    result.Points.Add(p2)

    'Making the line visible
    result.Stroke = Brushes.Black
    result.StrokeThickness = 2

    'Return the polyline
    Return result
End Function

我还创建了一个自己的圆类，其中包含所有必需的函数。

''' <summary>
''' A class for preforming circle tests
''' </summary>
''' <remarks></remarks>
Public Class Circle
    Private _CenterPoint As Point
    Public Property CenterPoint() As Point
        Get
            Return _CenterPoint
        End Get
        Set(ByVal value As Point)
            _CenterPoint = value
        End Set
    End Property

    Private _radius As Double
    Public Property Radius() As Double
        Get
            Return _radius
        End Get
        Set(ByVal value As Double)
            _radius = value
        End Set
    End Property

    'Calculates the Eucladian distence between two poitns
    Private Function Distance(ByVal p1 As Point, ByVal p2 As Point) As Double
        Return Math.Sqrt((p1.X - p2.X) ^ 2 + (p1.Y - p2.Y) ^ 2)
    End Function

    ''' <summary>
    ''' Checks if a Point1 is inside the circle
    ''' </summary>
    ''' <param name="p1">Test point</param>
    ''' <returns>If point is on the circle it is considered to be outside,
    ''' as three or more points is checked for</returns>
    ''' <remarks></remarks>
    Public Function DoesCircleContainPoint(ByVal p1 As Point) As Boolean
        Dim d As Double
        d = Distance(p1, CenterPoint)
        'Replace this with <= if you just want three points to form a circle
        If d < Radius Then
            Return True
        Else
            Return False
        End If
    End Function

    ''' <summary>
    ''' Returns the centerpoint of the created circle
    ''' </summary>
    ''' <param name="can"></param>
    ''' <remarks></remarks>
    Public Sub ReturnCenterPoint(ByVal can As Canvas)
        Dim p As New Ellipse
        p.Fill = Brushes.Blue
        p.StrokeThickness = 2
        p.Stroke = Brushes.Blue

        ' Set the width and height of the Ellipse.
        p.Width = 5
        p.Height = 5

        can.Children.Add(p)
        Canvas.SetLeft(p, CenterPoint.X - 2.5)
        Canvas.SetTop(p, CenterPoint.Y - 2.5)
    End Sub

    ''' <summary>
    ''' Returns the complete circle created from the three points
    ''' </summary>
    ''' <param name="can"></param>
    ''' <remarks></remarks>
    Public Sub ReturnEllipse(ByVal can As Canvas)
        Dim p As New Ellipse
        p.Fill = Brushes.Yellow
        p.StrokeThickness = 2
        p.Stroke = Brushes.Black

        ' Set the width and height of the Ellipse.
        p.Width = Radius * 2
        p.Height = Radius * 2

        can.Children.Add(p)
        Canvas.SetLeft(p, CenterPoint.X - Radius)
        Canvas.SetTop(p, CenterPoint.Y - Radius)
    End Sub

    Public Sub CreatCircleFromThreePoints(ByVal Point1 As Point, _
           ByVal Point2 As Point, ByVal Point3 As Point)
        Dim Pt(3) As Point
        Pt(0) = Point1
        Pt(1) = Point2
        Pt(2) = Point3

        GetCircleRectFromPoints(Pt)
    End Sub

    Private Sub GetCircleRectFromPoints(ByVal Pts() As Point) 'As RectangleF
        '(X - Cx)^2 + (Y - Cy)^2 = R^2
        '(R^2 - Cx^2 - Cy^2) + 2X*Cx + 2Y*Cy = X^2 +  Y^2
        'Cr + A*X + B*Y = X^2 + Y^2
        'Solve matrix using gaussian elimination.
        Dim N As Integer = Pts.Length - 1
        Dim X(N, N) As Double, Y(N) As Double
        For I As Integer = 0 To N
            X(I, 0) = 1 : X(I, 1) = Pts(I).X : X(I, 2) = Pts(I).Y
            Y(I) = X(I, 1) * X(I, 1) + X(I, 2) * X(I, 2)
        Next
        Dim MatInv As New Elimination
        MatInv.ComputeCoefficents(X, Y)
        Dim A As Single = CSng(Y(1) / 2)
        Dim B As Single = CSng(Y(2) / 2)
        Dim R As Single = CSng(Math.Sqrt(Y(0) + A * A + B * B))
        CenterPoint = New Point(A, B)
        Radius = R
    End Sub

    Public Class Elimination
        Sub ComputeCoefficents(ByVal X(,) As Double, ByVal Y() As Double)
            Dim I, J, K, K1, N As Integer
            N = Y.Length
            For K = 0 To N - 1
                K1 = K + 1
                For I = K To N - 1
                    If X(I, K) <> 0 Then
                        For J = K1 To N - 1
                            X(I, J) /= X(I, K)
                        Next J
                        Y(I) /= X(I, K)
                    End If
                Next I
                For I = K1 To N - 1
                    If X(I, K) <> 0 Then
                        For J = K1 To N - 1
                            X(I, J) -= X(K, J)
                        Next J
                        Y(I) -= Y(K)
                    End If
                Next I
            Next K
            For I = N - 2 To 0 Step -1
                For J = N - 1 To I + 1 Step -1
                    Y(I) -= X(I, J) * Y(J)
                Next J
            Next I
        End Sub
    End Class
End Class

主代码

在我撰写本文时，我很快意识到，如果我包含创建完整 Voronoi 图所需的所有必要步骤，它将会变得非常庞大。因此，我只完成了一个版本，该版本绝不包含创建完整 Voronoi 图的所有内容，而是让您有机会理解创建它所需的最基本功能。最终的完整实现和创建将在下一篇文章中完成。

历史

首次发布：2012 年 6 月 29 日。

第二次发布：2012 年 7 月 1 日，语言略有改动，在文章中添加了一些参考文献。

参考文献  

© . All rights reserved.