计算几何在MATLAB R2009a，第一部分

发布的罗兰美国舒尔，二零零九年七月十五日

14次观看(过去30天)| 0喜欢| 15个评论

我很高兴向大家介绍达米安·希伊作为本周的客座博主。Damian是The MathWorks的几何开发人员，他的背景是几何建模和网格生成领域。

内容

一个叫爱扑塞隆的生物

当我开始写这篇文章时，我回顾了罗兰过去的一些博客，以刷新我对过去的看法。当我看了她的博客浮点精度一瞥，浮点数的另一课,共线性我知道我来到了熟悉的领域。很有可能，您在过去的编程过程中也遇到过。

我第一次遇到的ulp(单位在最后的地方)和他的同伙每股收益90年代初我还是个研究生的时候我的导师和我在研究一个算法这个算法是基于德劳内三角．三角测量算法经常会因为数值问题而“崩溃”。为了让我们的数据集成功地三角测量，我浪费了精力去“微调”公差。如你所知，这就像试图把一个皱纹从一个完整的地毯。果不其然，不久之后，我的导师就会展示一个来自研究合作者的新数据集，我就可以重新开始工作了。令人沮丧的部分是缺乏关于如何处理这些故障的信息，因为它们很少在文本或技术出版物中得到解决。当我从研究生院毕业后，情况开始好转。处理几何计算的技术期刊和会议开始征求关于鲁棒性的研究论文。这些年来发表的论文令人鼓舞，我有一种宽慰的感觉，帮助正在到来。

为什么几何计算容易受到鲁棒性问题的影响

在罗兰的博客上共线性，我们意识到三个点的共线性检验可以看作是计算由三个点组成的三角形的面积，然后检查零面积。这将几何检验简化为计算行列式。事实上，这个几何测试还可以用来判断一个查询点是位于由两点定义的有向线的左边还是右边。来，试试吧…

P1 = [1 1];P2 = [5 5];P3 = [2 4];plotPointLine (p1, p2, p3);轴平等的；pointLineOrientation (p1, p2, p3)

ans = LEFT

Q1 = [1 1];Q2 = [5 5];Q3 = [4 2];pointLineOrientation (q1、q2、q3)

ans = RIGHT

dbtypepointLineOrientation

2 mat = [p1(1)-p3(1) p1(2)-p3(2);．.．3 p2(1)-p3(1) p2(2)-p3(2)];4 detmat = det(mat);5 if(detmat == 0) 6 orient = 'ON';如果(detmat < 0) 8 orient = 'RIGHT';9 else orient = 'LEFT';10结束11结束12

dbtypeplotPointLine

1函数plotPointLine(p1,p2,p3) 2 P = [p1;p2;p3];3情节(P (1:2, 1), P (1:2, 2));4 .坚持;5图(p3 (1), p3 (2) * r) 6 ptlabels = arrayfun (@ (n) {sprintf (% d, n)}, (1:3) ');7 Hpl =文本(P (: 1), P (:, 2), ptlabels,‘FontWeight’,‘大胆’);

如果点在直线的左边，行列式的符号为正，对应于一个面积为正的三角形。相反，如果这个点在另一边。正如我们之前看到的，如果查询点在直线上，所有三个点都是共线的，行列式为零。几何计算中的许多问题可归结为求行列式的符号。二维德劳内三角剖分的构造基于两个几何测试。我们刚刚看到的点线方向测试，以及用于确定查询点是否在由三角形的三个顶点定义的圆周内的圆内的点内测试。令人惊讶的是，你只需要在高中学过的两个简单的几何测试和一些编程知识就可以编写一个2D Delaunay三角剖分算法。然而，如果你真的尝试这样做，你会发现事情有时会出错——非常错误。

鲁棒Delaunay三角测量

Delaunay三角剖分算法的完整性取决于几何测试中使用的行列式的正确评估。当行列式模糊地接近于零时，问题就出现了。由于数值舍入的存在，方向测试可能会返回不正确的结果，然后算法可能会失败。

20世纪90年代初开发的Qhull计算几何库在计算过程中发现了这些数值精度问题。Qhull通过警告或错误消息向用户传达了问题，并提供了一些有用的交互选项来尝试解决问题。

例如，如果我们使用基于qhull的delaunayn方法计算单位正方形角点的Delaunay三角剖分，则会发出数值精度警告;这是因为这四个点都是共圆的。Qhull提供了一个“解决”这些问题的选项;“QJ”选项为这些点添加了噪音，以避免这种退化的情况。

X = [0 0;1 0;1 1;0 1] tri = delaunayn(X， {“QJ”})

X = 00 1 0 1 1 0 1警告:qhull精度警告:初始船体狭窄(余弦最小角度为1.0000000000000000)。共面点可导致宽面。选项'QbB'(缩放到单位框)或'QbB'(缩放到最后一个坐标)可能会删除此警告。使用“Pp”可以跳过此警告。参见http://www.qhull.org/html/qh-impre.htm tri = 1 2 4 4 2 3中的“限制”

虽然这种解决方法是有帮助的，但不能保证是健壮的，并且在实践中仍然可能出现失败。幸运的是，自20世纪90年代初以来，技术已经取得了进步。

在Loren早先关于共线性的博客中，后续评论引用了与行列式计算相关的数值不准确性，并推荐使用rand/SVD，因为它在数值上更可靠。但是，减少数值舍入也不能保证几何测试的正确结果。该实现仍然容易受到数值失败的影响。在过去的十年中，研究人员将注意力集中在诸如此类的稳健性问题上，特别是在Delaunay三角测量的稳健性计算上。对于Delaunay三角剖分，普遍接受的解决方案是基于精确几何计算(EGC)的概念。

罗兰的一个回复是这样的共线性博客强调了这个解决方案，但并没有引起人们的关注。这个想法是基于使用精确的算术计算行列式。由于这在计算上是昂贵的，所以应用了两步过程;行列式通常采用浮点算法计算，并在计算过程中计算出相应的误差范围。使用误差界限过滤掉模糊情况，然后对这些情况应用EGC以确保得到正确的结果。

在R2009a中，我们采用了来自计算几何算法库的2D和3D Delaunay三角测量(CGAL)在MATLAB中提供更健壮、更快和内存高效的解决方案。金宝搏官方网站CGAL采用EGC和浮点滤波器来保证数值的鲁棒性。让我们看看EGC在共线测试中的表现。我们将使用新的类DelaunayTri，这是基于CGAL的一个小实验。

共线检验重述

让我们再来看看洛伦共线检验，首先选择3个已知共线的点。如果我们试图构造德劳内三角测量使用DelaunayTri，则不应形成三角形，因为这三点共线。

格式长P1 = [1 1];P2 = [7500 7500];P3 = [750000 750000];

下面是基于行列式计算的测试

共线= (det([p2-p1;P3-p1]) == 0)

共线= 1

而测试则基于等级计算。

共线= (rank ([p1;p2;P3]) < 2)

共线= 1

P = [p1;p2;p3];

这是基于EGC的测试。

dt = DelaunayTri(P)

dt = DelaunayTri属性:约束:[]X: [3x2 double]三角测量:[]

三角测量是空的［］因为这些点是共线的，这就是我们所期望的。现在移动点p1使这三个点不共线。

P1 = [1 1+eps(5)];

下面是基于行列式计算的测试

共线= (det([p2-p1;P3-p1]) == 0)给出不正确的结果

共线= 1

而测试则基于等级计算。

共线= (rank ([p1;p2;P3]) < 2)给出不正确的结果

共线= 1

最后，这是基于egc的测试。

P = [p1;p2;p3];dt = DelaunayTri(P);共线= isempty(dt(:，:))给出正确的结果

共线= 0

在这种情况下，一个三角形被构造，这表明三个点不被认为是共线的。我们知道这是正确的评估。

我们能从中得出什么结论?

EGC在Delaunay三角剖分、Voronoi图和凸包等算法的鲁棒实现中发挥着非常重要的作用。长期以来，这些算法都很容易受到数值精度问题的影响。

这引发了一个有趣的问题;我们是否应该在MATLAB中采用精确的几何测试来开发通用算法?例如;我们是否应该使用精确共线测试和精确点线方向测试之类的?

一般来说，我们不应该这样做，因为我们工作的环境是有限精度的。如果我们广泛地使用EGC，我们就会意识到，虽然它在数值上是正确的，但它可能无法捕捉到我们的编程意图。在算法的流程中，我们从计算的一个步骤中获取输出，并将其作为输入传递给下一个步骤。

如果我们将不精确的输入传递给精确的几何测试会发生什么?假设我们从三个共线点开始计算，接下来我们应用旋转矩阵以使这些点在特定方向上重新定向。在这样做的过程中，我们引入了不准确性;精确的几何测试能纠正这些错误吗?不，不会;你自己看吧。

P1 = [1 1];P2 = [7500 7500];P3 = [750000 750000];P = [p1;p2;p3];T = [cos(/3) sin(/3);sin(π/ 3)因为(π/ 3)];P2 = p * t;dt = DelaunayTri(P2);共线= isempty(dt(:，:))

共线= 0

由点{p1, p2, p3}构成三角形，因此根据EGC，这三个点被认为是非共线的。

虽然结果在技术上是正确的，但可能不是你所期望的。在Delaunay三角测量及其应用的背景下，不精确输入的后果是相对良性的。EGC是解决该领域长期存在的健壮性问题的强大工具。

在我们有限精度的编程世界中，明智地选择公差通常足以捕获正确的行为。这就引出了一个问题;什么是宽容的明智选择?例如，在共线性测试中使用合适的公差是多少?好的，在这种情况下，一个好的容差可以在动态或相对意义上捕捉行列式计算中的损失。当点的坐标很小时，它应该工作，当点的坐标很大时，它应该工作。预定义的固定值不能在很大的坐标值范围内工作。当我计算用于几何测试的公差敏感行列式时，我检查了膨胀过程中产生的产品。下载188bet金宝搏然后我选择的产品与最大的幅度，并使用它来衡量eps。

例如，当计算矩阵的行列式时[b;c d]，我选择大小，

Mag = max(abs(a*d)， abs(b*c))

然后计算相对公差，

relativeTol = eps(mag);

我应用容差如下。

共线= (abs(det([p2-p1;p3-p1]) <= relativeTol))

(为简单起见，产品下载188bet金宝搏* d而且b * c不是重用来计算的依据．)

底线是什么?

那么，讲精确几何计算有什么意义呢?哦，我希望这个概述能让你对几何计算中的数值问题有一些见解。但对用户来说，重要的一点是，在MATLAB中构建2D/3D Delaunay三角剖分、Voronoi图或凸包时，您不再需要担心数值鲁棒性问题。由于其他应用程序(如最近邻居查询和分散数据插值)都是基于这些特性构建的，因此您可以期望在这些特性上实现健壮性、性能改进和内存效率。如果您想了解更多，请查看视频以及在R2009a中突出新的计算几何特征的Delaunay三角测量演示。以下链接是新功能的参考页面: