Geekly articles weekly

财富算法，实现细节

在过去的几周中，我一直在努力用C ++实现Fortune的算法 。 该算法需要很多2D点，并根据它们构建Voronoi图 。 如果您不知道什么是Voronoi图，请看一下该图：


对于每个称为“站点”的入口点，我们需要找到比这个地方更近的地方。 这样的点集形成了上图所示的单元格。

在《财富》算法中，他及时地建立了这样的图表是很了不起的 O（n logn）（这是比较算法的最佳选择），其中 n是个地方。

我写这篇文章是因为我认为该算法的实现非常困难。 目前，这是我必须实现的最复杂的算法。 因此，我想分享我遇到的问题以及如何解决它们。

像往常一样，代码被发布在github上 ，并且我使用的所有参考资料都列在本文的结尾。

财富算法说明


我不会解释该算法的工作原理，因为其他人已经做得很好。 我可以建议您学习这两篇文章： here和here 。 第二个非常有趣-作者用Javascript编写了一个交互式演示，对于理解算法的操作很有用。 如果您需要所有证据都更正规的方法，建议阅读《 计算几何》第3版第7章。

此外，我更喜欢处理没有详细记录的实现细节。 正是他们使算法的实现变得如此复杂。 特别是，我将重点介绍所使用的数据结构。

我只是写了该算法的伪代码，以便您了解其全局结构：

 将地点事件添加到每个地点的事件队列中
直到事件队列为空
    检索顶级事件
    如果事件是地方事件
        在海岸线上插入新弧线
        检查新的圈子活动
    否则
        在图中创建一个顶点
        我们从海岸线上去除了拉紧的弧线
        删除无效事件
        检查新的圈子活动

图表数据结构


我遇到的第一个问题是选择存储Voronoi图的方式。

我决定使用一种广泛用于计算几何的数据结构，称为双连接边列表（DCEL）。

我的VoronoiDiagram类使用四个容器作为字段：

 class VoronoiDiagram { public: // ... private: std::vector<Site> mSites; std::vector<Face> mFaces; std::list<Vertex> mVertices; std::list<HalfEdge> mHalfEdges; }

我将详细讨论它们中的每一个。

Site类描述了入口点。 每个位置都有一个索引，这对于将其放入队列，坐标和指向单元格（ face ）的指针很有用：

 struct Site { std::size_t index; Vector2 point; Face* face; };

单元的Vertex由Vertex类表示，它们只有一个坐标字段：

 struct Vertex { Vector2 point; };

这是半边的实现：

 struct HalfEdge { Vertex* origin; Vertex* destination; HalfEdge* twin; Face* incidentFace; HalfEdge* prev; HalfEdge* next; };

您可能想知道，什么是半肋骨？ Voronoi图中的一条边对于两个相邻的单元是公用的。 在DCEL数据结构中，我们将这些边分成两个半边，每个边对应一个半边，并由twin指针链接。 而且，半边缘具有起点和终点。 incidentFace字段指示半边线所属的面。 DCEL中的单元被实现为半边的循环双向链接列表，其中相邻的半边连接在一起。 因此，上一个和next字段指示单元格中的上一个和下一个半边缘。

下图显示了红色半边的所有这些字段：


最后， Face类定义了单元格。 它仅包含一个指向其位置的指针，另一个指向其半肋的指针。 选择哪个半边都无所谓，因为该单元格是一个封闭的多边形。 因此，在遍历循环链表时，我们可以访问所有半边。

 struct Face { Site* site; HalfEdge* outerComponent; };

事件队列


实现事件队列的标准方法是使用优先级队列。 在处理地点和圈子事件的过程中，我们可能需要从队列中删除圈子事件，因为它们不再有效。 但是大多数标准优先级队列实现不允许您删除不在最前面的项目。 这尤其适用于std::priority_queue 。

有两种方法可以解决此问题。 第一个更简单的方法是向事件添加valid标志。 valid最初设置为true 。 然后，不必从队列中删除circle事件，我们只需将其标志设置为false 。 最后，在处理主循环中的所有事件时，我们检查事件的valid标志值是否为false ，如果是，则直接跳过该事件并处理下一个事件。

我应用的第二种方法是不使用std::priority_queue 。 相反，我实现了自己的优先级队列，该队列支持删除其中包含的任何元素。 这样的队列的实现非常简单。 我选择此方法是因为它使算法代码更整洁。

海岸线


海岸线数据结构是算法的复杂部分。 如果实现不正确，则不能保证算法将在 O（n logn）。 获得这种时间复杂性的关键是使用自平衡树。 但是说起来容易做起来难！

建议我研究的大多数资源（上文提到的两篇文章和《 计算几何》一书）将海岸线实现为一棵树，其中内部节点表示断点，叶子表示弧形。 但是他们没有说如何平衡一棵树。 我认为这样的模型不是最好的，原因如下：

  • 其中包含冗余信息：我们知道两个相邻弧之间有一个断点，因此没有必要将这些点表示为节点
  • 它不足以实现自我平衡：只有断点形成的子树才能达到平衡。 我们确实无法平衡整个树，因为否则弧会成为内部节点和断点的叶子。 编写仅平衡内部节点形成的子树的算法对我来说就像一场噩梦。

因此，我决定以不同的方式呈现海岸线。 在我的实现中，海岸线也是一棵树，但是所有节点都是弧形的。 这样的模型没有任何列出的缺点。

这是我的实现中Arc arc的定义：

 struct Arc { enum class Color{RED, BLACK}; // Hierarchy Arc* parent; Arc* left; Arc* right; // Diagram VoronoiDiagram::Site* site; VoronoiDiagram::HalfEdge* leftHalfEdge; VoronoiDiagram::HalfEdge* rightHalfEdge; Event* event; // Optimizations Arc* prev; Arc* next; // Only for balancing Color color; };

前三个字段用于构造树。 leftHalfEdge字段指示由弧的最左点绘制的半边。 rightHalfEdge位于最右点绘制的半边上。 两个指针prev和next用于直接访问海岸线的上一个弧线和下一个弧线。 此外，它们还允许您绕过海岸线，成为双向链接列表。 最后，每个弧线都有一种用于平衡海岸线的颜色。

为了平衡海岸线，我决定采用红黑方案 。 编写代码时，我受到《 算法简介 》一书的启发。 第13章介绍了两种有趣的算法， insertFixup和deleteFixup ，它们在插入或删除后平衡树。

但是，我无法使用书中所示的insert方法，因为键用于查找节点的插入点。 财富算法中没有键，我们只知道我们需要在海岸线的另一弧之前或之后插入弧。 为了实现这一点，我创建了insertBefore和insertAfter ：

 void Beachline::insertBefore(Arc* x, Arc* y) { // Find the right place if (isNil(x->left)) { x->left = y; y->parent = x; } else { x->prev->right = y; y->parent = x->prev; } // Set the pointers y->prev = x->prev; if (!isNil(y->prev)) y->prev->next = y; y->next = x; x->prev = y; // Balance the tree insertFixup(y); }

通过三个步骤在x之前插入y ：

  1. 寻找一个插入新节点的地方。 为此，我使用了以下观察：左子项x或右子项x->prev是Nil ，而Nil的子项在x之前和x->prev 。
  2. 在海岸线内部，我们保持双向链表的结构，因此必须相应地更新元素x->prev ， y和x的prev和next指针。
  3. 最后，我们只需调用本书中介绍的insertFixup方法来平衡树。

insertAfter的实现类似。

从书中删除的方法可以直接实施。

图表限制


这是上述财富算法的输出：


图像边界上的单元格的某些边缘存在一个小问题：由于它们是无限的，因此未绘制它们。

更糟糕的是，一个单元可能不是单个片段。 例如，如果我们在同一条线上取三个点，那么中点将具有两个没有连接在一起的无限半边。 这不太适合我们，因为我们将无法访问半边之一，因为该单元格是边的链接列表。

为了解决这些问题，我们将限制图表。 我的意思是，我们将限制图中的每个单元，以使它们不再具有无限的边缘，并且每个单元都是封闭的多边形。

幸运的是，Fortune的算法使我们能够快速找到无尽的边缘：在算法结束时，它们对应于仍在海岸线中的半边缘。

我的限制算法收到一个框作为输入，包括三个步骤：

  1. 它提供了图的每个顶点在矩形内的位置。
  2. 切掉每一个无限的边缘。
  3. 关闭单元格。

第1阶段很简单-如果矩形不包含顶点，我​​们只需要对其进行扩展即可。

第2阶段也非常简单-它包括计算射线与矩形之间的交点。

第三阶段也不是很复杂，只需要注意即可。 我分两个阶段执行。 首先，我将矩形的角点添加到应位于其顶点处的单元格中。 然后，确保单元的所有顶点都通过半边连接。

我建议您研究代码并询问是否需要有关此部分的详细信息。

这是边界算法的输出图：


现在我们看到所有边缘都被绘制了。 如果缩小，则可以确保所有单元格都已关闭：


矩形相交


太好了！ 但是从文章开始的第一张图片更好，对吧？

在许多应用中，使Voronoi图与矩形相交是很有用的，如第一幅图所示。

好消息是，在限制图表之后，这容易得多。 坏消息是，尽管算法不是很复杂，但我们仍需小心。

这个想法是这样的：我们绕过每个像元的半边并检查半边与矩形之间的交点。 可能有五种情况：

  1. 半肋完全位于矩形内部：我们保存了一个半肋
  2. 半肋完全位于矩形的外部：我们丢弃了这样一个半肋
  3. 半肋超出矩形的范围：我们截断半肋并将其保存为最后出现的半肋 。
  4. 半肋进入矩形内部：我们将半肋截断以将其与最后一个出去的半肋连接起来（在3或5的情况下将其保存）
  5. 半肋穿过矩形两次：我们截断半肋，并添加一个半肋以使其与最后一个出的半肋关联，然后将其保存为新的最后一个半肋 。

是的，有很多情况。 我创建了一张图片以显示所有图片：


橙色多边形是原始单元格，红色是截断的单元格。 截断的半肋标记为红色。 添加了绿色肋骨，以连接进入矩形的半肋与半肋出来。

将此算法应用于有界图，我们可以得到预期的结果：


结论


原来这篇文章很长。 而且我敢肯定，您仍然不清楚很多方面。 但是，我希望它对您有用。 检查代码以获取详细信息。

总结并确保我们不会白白浪费时间，我在（便宜的）笔记本电脑上测量了计算不同地方的Voronoi图的时间：

  • n=1000：2毫秒
  • n=$10,00：33毫秒
  • n=$100,00：450毫秒
  • n=$1,000,00：6600毫秒

我没有什么可与这些指标进行比较，但是看来它是如此之快！

参考文献


Source: https://habr.com/ru/post/zh-CN430628/

More articles:


All Articles