地形中的大面积平面区域计算及金属共聚焦激光扫描显微镜技术研究

133理论计算机科学电子笔记46（2001）网址：http://www.elsevier.nl/locate/entcs/volume46.html13页地形中大面积平面区域的计算Katharina Lange1，2，Rahul Ray1，3，Pastel Smid1，4和Ulrich Wendt1，5材料科学Otto-von-Guericke University ofMagdeburg D-39106 Magdeburg，Germany摘要我们考虑的问题计算的最大区域的地形，是approximately包含在一些二维平面。我们把这个问题简化为下面的问题。 给定3度图G在单位球面上的一个嵌入S2，其顶点是加权的，计算包含在某个固定半径的球形圆盘中的最大权连通子图。我们给出了一个算法，解决这个问题的O（n2logn（log logn）3）时间，其中n表示G的顶点数，或者，地形的面数。 我们也给一种启发式算法，可用于计算近似平面的地形中的足够大的区域。我们讨论了一个基于网络的实现这种启发式，并显示了一些结果的地形表示三维（地形）图像的断裂表面的金属共聚焦激光扫描显微镜。1介绍地形是R3中的曲面，由两个实变量的实值函数描述。如果该函数是分段线性的，并且表面由三角形的集合组成，则地形称为不规则三角网（TIN）。我们考虑的问题，计算最大的连接区域，近似平面的三角网。更精确地说，给定TINT，我们要计算T的三角形的子集T，使得（i）对于每个三角形t∈T，存在另一个三角形tJ∈T，使得t和tJ共享一条边，1作者得到了德国研究共同体的支持，赠款WE 2301/3- 2和SM 57/4-12 电子邮件：Katharina.Lange@mb.Uni-Magdeburg.DE3电邮地址： rahul@isg.cs.uni-magdeburg.de4电邮地址： michiel@isg.cs.uni-magdeburg.de5电子邮件：Ulrich.Wendt@mb.Uni-Magdeburg.DE2000年1月，出版社dbyElsevierScienceB。 V. 在CCBY-NC-ND许可下开放访问。兰格、雷伊、斯米德和温特134+++不不不++(ii)存在一个二维平面，该平面近似包含T的所有三角形，并且（iii）T的三角形的总面积尽可能大为了严格地定义这个问题，我们必须指定近似包含在平面中的概念。在本文中，我们将这个概念定义如下。设T是一个小的正实数，T是一组满足上述（i）的我们说T是n-平面的，如果存在一个向量c使得T中所有三角形的法向量与c的夹角至多为n。现在我们展示如何用这个概念来重新表述我们的问题。单位球，即，以原点为中心并具有半径1的三维球的边界由S2表示。上半球定义为S2：=S2<${（x，y，z）∈R3：z> 0}。请注意，我们可以将R3中任何非垂直三角形的法向量作为S2上的一点。设T是一个三角网，设SS2是三角网的法向量集三角形在 设G=（S，E）是具有顶点集的无向图其中任意两个顶点由一条边连接当且仅当中的相应三角形共享一条边。（实际上，S是一个多重集，因为不同的三角形可以有相同的法线。 相等的法线被视为G中的不同顶点。）注意G的每个顶点的度最多为3。我们给G的每个顶点p一个权重wt（p），它等于G的面积在三角形中，S的任意子集C的重量（C）为定义为wt（C）：=p∈Cwt（p）。对任意点x∈S2，设Dx表示以x为中心的半径为λ的球面盘. 也就是说，Dx是所有点y ∈ S2的集合，使得向量x和y小于或等于x。进一步，设Gx表示G的一个子图，其顶点集为S<$Dx，其边集为所有边（p，q）∈E的集合，其中p和q都包含在Dx中.最后，我们定义Wx为任意连通分量的最大权重，图Gx.使用这个术语，我们的问题可以正式表述如下。问题1.1给定一个有n个顶点和一个实常数的图Gn>0，计算一个点x∈S2使得Wx最大。设（p，q）是图G的任意边.如果向量p和q之间的夹角大于2，那么很明显，在解决问题1.1时，可以忽略（p，q）。因此，我们可以不失一般性地假设G的任何边的端点（当被视为向量时）之间的角度至多为2π。1.1动机本文中考虑的问题出现在马格德堡大学的计算机科学和材料科学系之间的合作。该项目的目标是设计和实现算法，可用于定量的断裂表面形貌，给定的兰格、雷伊、斯米德和温特135×+作为通过共聚焦激光扫描显微镜获得的三维图像。我们使用的图像取自金属的断裂表面。 知道- 材料断口边缘形貌可以用来解释材料的力学性能和断裂机理。这些知识可以帮助检测材料微观结构中的弱点。这些关系可以作为改进和修改已知材料以及新的材料设计标准这适用于所有类型的材料，如金属、聚合物、陶瓷和复合材料。这些图像被给定为512 512像素阵列，其中存储在每个条目处的值等于对应像素的高度。我们研究的一个目标是在这幅图像中找到近似平面的大的连通区域这些平面区域的法向量和面积给出了关于断裂面生成过程的有用信息。在第4节中，我们将看到如何将像素数组转换为TIN。1.2我们的结果在第2节中，我们将展示如何计算几何和动态图算法可以用来解决问题1.1在O（n2logn（log logn）3）时间。问题1.1不太可能在次二次时间内解决在事实上，似乎即使是计算点y∈S2的问题，Wy近似于最优解，不能在次二次时间内求解。因此，在第3节中，我们描述了一个简单的基于网格的启发式算法。我们已经实现了这种启发式，以便它可以通过万维网使用我们在第4节中给出了有关此实现的一些细节。在第5节中，我们提出了一些实验结果的图像上获得的断裂表面的共聚焦激光扫描显微镜。这些结果表明，该启发式算法能够找到面积足够大的非平面区域1.3相关工作本文考虑的问题与地形简化问题有关。在这个问题中，我们想用一个“更小”的地形来近似一个多面体地形，即，其中一个具有最小数目的顶点。虽然这个问题已经在例如计算机图形社区中进行了研究，但我们知道从复杂性角度考虑这个问题的主要参考文献是Agarwal和Suri [2]。他们通过证明一个强相关问题是NP难的，证明了地形简化问题是困难的。他们还给出了一个多项式时间算法近似的最小地形。如果我们考虑问题1.1的情况下，当图G是完全的，所有的顶点都有相等的权重，那么我们就得到了计算的问题，一个球面盘的放置，它包含给定的S2上n个点的集合S的最大子集。Chazelle已经考虑过这个问题，兰格、雷伊、斯米德和温特136+李[4]的情况下，当S是一个点集在欧几里德平面。这个问题的时间复杂度为O（n）。 相关问题计算半平面排列中的最深点是3SUM-困难的，参见Gajentaan和Overmars [6]。这表明在次二次时间内解决磁盘布局问题可能是非常困难的。最近，Agarwal等人 [1]给出了求解最优磁盘布局问题的另一种O（n2）时间算法，以及运行时间接近线性的近似算法。我们在第2节中求解问题1.1的算法是受[1]中O（n2）时间算法的启发2解决问题1.1在本节中，我们给出一个解决问题1.1的算法。考虑图G=（S，E），并考虑以点为中心的球形圆盘Dpp，S。设A是S2上由圆盘Dp的边界定义的排列，其中p∈S。也就是说，A是将S2细分为顶点、边和面，这些顶点、边和面由圆盘Dp，p∈S的边界的重叠来定义。由于p∈Dx当且仅当x∈Dp，我们有Wx=Wy，对于任意两个点x和y，它们在A的同一个面f的内部。此外，对于f的每个顶点z，我们有Wz≥Wx。这证明了下面的引理。引理2.1为了解决问题1.1，假设点x∈S2，是排列A的顶点。在这一节中，我们对集合S作如下的一般位置描述。我们假设S的元素是两两不同的。此外，我们还假设对于S的任意两个不同的点p和q，球面圆盘Dp和Dq要么不相交，要么有正面积的交集（因此，Dp和Dq互不接触）.最后，对于S的任意三个不同的点p、q和r，球面圆盘Dp、Dq和Dr不相交于一个点。我们做这个假设只是为了简化我们算法的描述。该算法可以很容易地扩展到处理任意的点集。上面的讨论导致下面的初步算法解决问题1.1。步骤1：计算排列A。步骤2：设W：= 0。 对于A的每个顶点x，执行以下操作。• 计算图Gx。• 计算Gx的连通分量，并为每个分量计算其权重。这些连通分量的最大权给出了Wx的值。• 设置W：= max（W，Wx）。第三步：返回W。兰格、雷伊、斯米德和温特137一一+一显然，该算法正确地解决了问题1.1。让我们分析一下它的运行时间。回想一下，n表示点集S的元素数。 由于图G的每个顶点的度至多为3，因此G至多有3n/ 2条边。算法的时间复杂度为O（n2），所以算法的复杂度为O（n2Amato等人的出租法[3]。考虑A的任意顶点x。图Gx、它的连通分支和值Wx可以在O（n）时间内计算以来时间复杂度为O（n2）。因此，整个算法需要O（n3）时间。我们现在展示如何显著提高运行时间。请注意，前面算法的瓶颈是步骤2。改进的想法算法是遍历布置A并保持连接的com。在一个数据结构中的Gx的分量。考虑一下当我们沿着A的一条边从一个顶点x走到相邻的顶点y时会发生什么。此时，我们知道图Gx的连通分支，并且想要计算图G x的连通分支。Gy的连通分量从x移动到y意味着我们沿着的边缘移动球面圆盘Dx到位置Dy。在这个运动过程中，S的最多一个点可以进入或离开球面盘。(Here我们使用一般位置假设。）由于图G的度为3，因此可以通过执行至多常数数量的边插入和删除来从Gx获得图Gy我们假设我们有一个数据结构CC，它存储图的连接组件及其权重，并且支持边插入和删除，以及形式为“报告任何连接组件的最大权重”的查询。 (We稍后将详细说明此数据结构。为现在，我们把它当作一个黑盒子。对于任意点x∈S2，我们记为：CCx是图Gx的这个数据结构的实例。我们的改进算法做了以下工作。步骤1：计算排列A。步骤2：设x是A的任意顶点。• 计算图Gx。• 计算Gx的连通分量，并为每个分量计算其权重。计算Wx作为Gx的任何连通分支的最大权重。• 设置W：=Wx。• 构造数据结构CCx。步骤3：从x开始，遍历排列的顶点，例如，在深度-第一顺序。在一般步骤中，我们从顶点y走到相邻顶点z。此时，我们有数据结构CCy，存储图Gy的连接分量以及它们的权重。图Gz可以通过插入和删除至多一个常数个兰格、雷伊、斯米德和温特138++≤≤当前图Gy中的边。因此，我们通过在数据结构CCy中执行这些更新来获得数据结构CCz。然后，我们查询CCz以找到Wz的值，并设置W：= max（W，Wz）。步骤4：返回W。该算法的正确性是显而易见的。时间复杂度为O（n2），时间复杂度为O（1）。步骤2和步骤3的时间取决于数据结构CC。设P（n）、U（n）和Q（n）表示预处理时间，更新时间和查询时间。时间复杂度为O（n+p）。 在步骤3中，我们花费O（U（n）+Q（n））时间来计算每个A的顶点 由于这种安排有O（n2）个顶点，因此可以得出，算法的总运行时间为.2ΣO P（n）+n（U（n）+Q（n）.它仍然指定数据结构CC。在[9]中，Thorup给出了一种数据结构，它能在每次更新的平均时间为O（logn（log logn）3）的情况下，在边的插入和删除下，保持图的生成森林，其中n表示图的顶点数.给定图的任意两个顶点，可以在O（logn/ log log logn）时间内确定它们是否在同一连通分支中（一个更简单但理论上效率稍低的数据结构由Holm等人给出 [7]。该数据结构可以很容易地扩展，以便在同一时间范围内保持所有连接组件的权重。如果我们将这些权重存储在堆中，那么我们可以在O（1）时间内提取最大连通分量的权重。而且，这个堆可以在每个操作的O（logn）时间内更新。 数据结构可以通过将所有边连续插入到最初为空的图中来构建。因此，我们有P（n）=O（nlogn（log logn）3），U（n）=O（logn（log logn）3），并且Q（n）=O（1）. 因此，我们证明了以下结果。因此，时间复杂度为O（n2log n（log log n）3）。在1.3节中，我们论证了问题1.1不太可能在次二次时间内得到解。相反，人们可以询问近似最优解的时间复杂度。也就是说，设δ是一个实数，0<δ 1，设x∈S2计算一个点y∈S2是一个点，其中Wx是最大值。我们的目标是使得Wy≥（1 −δ）Wx。考虑下面的例子。设D是半径为λ的球面圆盘，设p和q是D的边界上两个径向相对的点。设m是一个大整数。对于每个i，1≤i≤m，设ai：=p和bi：=q。考虑边（ai，bi），1≤i≤m，和（bi，ai+1），1≤im。请注意，这些边缘形成在点p和点q之间交替的a1和bm之间的路径。设G是一个包含点ai和bi的图，1我m作为顶点，以及上面的边缘。这个图包含更多的顶点和边，但是所有这些顶点都 我们假设G的所有顶点都有单位权。 (It很容易以构造一个TIN，其中G是对应的图。）对于这张图，兰格、雷伊、斯米德和温特139不不++∈不不--D的中心给出问题1.1的最优解。如果我们移动圆盘D，则点p或点q离开圆盘，因此，包含在圆盘中的G的每个连通子图都由一个顶点组成。这表明问题1.1的任何近似算法都必须返回D的中心。正因为如此，我们认为，这是困难的，以解决问题1.1的近似版本在次二次时间。3一种寻找大平面区域的启发式算法在本节中，我们给出了一个简单的启发式方法来计算三角网中的大连通区域，这些区域是平面的。(Di不同的区域不需要近似包含在同一个二维平面中。）我们不能证明其输出质量的任何非平凡的界限，但实验表明，输出在实践中是好的，启发式是快速的。设是一个由n个三角形组成的TIN，令S>0，令S是这些三角形的法向量的集合，令G=（S，E）是定义的图第1节。回想一下，S的任何元素p的权重wt（p）等于法向量为p的三角形的面积。此外，还记得S实际上是一个多重集。我们的目标是找到所有面积至少为A的平面区域，其中A是某个给定的正实数。我们在上半球S2使用经度线，纬度，使得每个网格单元包含在半径- 是的更精确地说，我们选择一个合适的实数δ，δ = Θ（π），并使用X：=[π/δ|等间隔的经度线lo i，0 ≤i< X，且Y：=[π/（2 δ）]|等间距的纬线la j，0 ≤j< Y。对于任意两个指数i和j，其中0≤i X和0≤j Y，我们称之为pair（i，j）由loi、loi+1、laj和laj+1界定的网格单元的索引。为在O（1）时间内计算出包含p的网格单元的索引。注意，每个网格单元最多与四个其他单元相邻除了那些发生在北极的。我们的启发式算法将TIN、正实数A和A以及图G=（S，E）作为输入。它首先计算所有面积至少为A/4的平面区域的子集。（下面，我们将清楚为什么我们选择A/4而不是A。事实上，因子1/ 4可以用0到1之间的任何常数代替。）第一步：初始化一个数组C [0.. X1，0... Y1]，并与每个条目一起存储空的点列表和空的边列表步骤2：对于每个点p∈S，计算包含p的网格单元的索引（ip，jp），并将p添加到C[ip，jp]的点列表中步骤3：对于每个边（p，q）∈E，其中（ip，jp）=（iq，jq），将（p，q）添加到C[ip，jp]的边列表中。第四步：初始化一个空列表L。兰格、雷伊、斯米德和温特140不不不不不步骤5：对于每个i和j，其中0≤i X和0≤j Y，执行以下操作。• 计算以C[i，j]的点列表为顶点集、以C[i，j]的边列表为边集的图Gij• 对于Gij的每个连通分支，如果其权重至少为A/4，则将其添加到列表L步骤6：返回列表L。考虑在步骤6中返回的列表L。L的每个元素对应于G的一个连通子图，其权重至少为A/4，并且包含在一个网格单元中，即，它对应于TIN中面积至少为A/4的平面区域。可能发生的是，这样的区域R可以通过添加与R相邻并且其法线为在与产生R的网格单元相邻的网格单元中。 当然，放大区域应该仍然是平面的。下面，我们描述一个步骤7（边界校正步骤）：对存储在列表L中的G的每个连接子图执行该步骤.设GJ是任意一个这样的子图。我们首先计算包含所有的最小封闭球面圆盘DJGJ的顶点，使用所提出的线性时间算法，例如，在de Berg等人 [5]中。设c是DJ的中心，D是以c为中心的半径为π的球面盘。注意，GJ包含在D中，因为DJ包含在D中。现在我们回到三角网，标记与GJ的顶点对应的三角形。设T是所有有标记三角形的集合。 注意T中三角形的所有法线都包含在D中。现在我们考虑每个未标记的三角形t的与至少一个标记三角形共享一条边，如果其法线包含在D中，则标记t。 设TJ是在我们考虑了所有这样的三角形t之后的所有有标记三角形的集合。 如果T=TJ，那么我们停止边界对这个子图GJ的校正步骤。否则，我们重复这个过程，考虑与至少一个标记的三角形共享一条边的未标记的三角形，三角步骤8：在完成L中每个元素的边界校正步骤之后，我们有一个平面区域的集合，每个平面区域的面积至少为A/4。由于这些区域可能重叠，我们继续如下。我们根据它们的面积对所有这些放大的区域进行排序。如果最大的区域的面积至少为A，那么我们报告它，并标记它的所有三角形。然后我们考虑第二大区域。如果它的面积至少是A，并且如果它的三角形还没有被标记，我们也报告它并标记它的所有三角形。我们继续这样做，直到我们报告了面积至少为A且不与任何先前报告的区域重叠的所有区域。关于产品的质量，我们能说些什么呢？设R为TIN中的任意平面区域 ，设GJ是G. 我们的算法报告区域R当且仅当GJ包含连通兰格、雷伊、斯米德和温特141××图1.一、 将像素数组转换为TIN。子图GJJ的重量至少A/4，完全包含在我们的网格细胞之一。我们可以通过运行算法来提高找到R使用网格的移动副本多次。4基于Web的实现在本节中，我们将讨论上一节的启发式算法的实现。我们称之为PLANE FINDER的程序以tif（标记图像文件）格式的k l像素数组作为输入。存储在每个条目处的值是相应表面点的高度，即，体素位置的z值。PLANE FINDER首先计算TIN，如图1所示。 TIN的顶点是（i）像素的中心，其中z坐标由像素的高度给出，以及（ii）所有2 × 2像素块的中心，其中z坐标由组成块的四个像素的平均高度给出。这些顶点由图1右边部分所示的 注意，三角形的总数，我们用n表示，等于n = 4（k −1）（l−1）。给定此TIN，PLANE FINDER如第3节所述进行。该程序是用C++编写的，并使用LEDA库[8]。我们已经添加了一个网络界面的飞机发现，使用户可以上传tif-图像和运行程序从任何平台与Netscape 4.x或更高。此接口是用Common Gateway Interface脚本编写的，该脚本在我们的Web服务器（Sun Ultra Sparc计算机）上运行。图2显示了该网页的屏幕截图，http://wwwisg.cs.uni-magdeburg.de/~rahul/proj.html目前支持以下操作。• 对于用户提供的某些值k和k，显示k个最大的• 显示所有面积至少为A（同样是通过启发式发现的）的平面区域，对于用户提供的某些值A和A兰格、雷伊、斯米德和温特142×图二、 PLANE FINDER网页的屏幕截图。输出是原始的tif-image，其中所有被检测到的区域（当且仅当与像素重叠的三角形中至少有一个属于该区域时，像素才被着色为白色。）5实验结果我们已经运行了PLANE FINDER对共聚焦激光扫描显微镜获得的断裂表面的各种图像。这些图像是tif格式，由512个 512像素。 因此，相应的TIN包括n=1044484个三角形。这些表面中的大多数是非常粗糙的，并且仅包含小的平面区域。对于一个典型的图像，PLANE FINDER在我们的Web服务器上大约需要35秒。当程序通过网络使用时，上传和发回输出的时间必须增加。对于作者来说，在马格德堡家中的电脑上处理一个查询大约需要65秒。图3中给出了一个示例，其中显示了原始图像（左侧）以及PLANEFINDER的输出（右侧）。程序运行时，θ等于5度，A等于1500平方单位。的兰格、雷伊、斯米德和温特143图3.第三章。 在左侧，显示了钢断裂表面的共焦激光扫描显微照片（形貌图像）。颜色代表像素的高度由PLANE FINDER计算的四个平面区域在图中被标记为白色。右图像。图像尺寸为100 µ m×100 µ m。右边部分用白色示出了所发现的面积至少为A= 1500平方单位的四个平面区域这里，面积是指空间中的面积，而不是投影面积。图4显示了一个图像示例，其中我们添加了四个大的倾斜三角形。其中一个是平面的，其中θ= 5度，而其他的是对于稍微大于5度的值的平面。PLANE FINDER在此图像上运行，A= 3000平方单位。6总结发言我们已经考虑了在一个连接的和近似平面的地形中计算大区域的问题。我们证明了计算最大区域的问题可以在大约是地形中三角形数量的二次方的时间内解决，并认为这不太可能更快地解决问题。我们留下了严格证明这一点的问题。我们还认为，甚至可能很难近似这个最大的区域。正式证明这一主张也是一个悬而未决的问题。在本文中，我们定义了一个三角形的连通集是“近似平面的”，如果它们的法线包含在一个固定的小半径的球面盘中。对于近似平面的其他概念，解决这个问题将是有趣的。例如，我们可以要求法线的角直径最大为100，或者三角形集包含在两个距离为100的平行平面之间。我们的软件还处于初级阶段，有很大的改进余地。我们计划通过在预处理步骤中应用均值图像滤波器这种过滤器将消除局部粗糙度。兰格、雷伊、斯米德和温特144见图4。在顶部，显示地形图像及其在xy平面上的平行投影。出于测试目的，插入了一个“近似平面”三角形和三个“不太近似平面”三角形。由PLANE FINDER计算的平面区域在底部图像中标记为白色。目前，我们正在进行内存优化。此外，我们将添加各种功能来改善用户交互。在PLANE FINDER的未来版本中，相邻的平面区域将以明显不同的颜色绘制，并为用户提供关于平面区域的更多信息。例如，通过点击一个区域，将显示该区域的此外，我们将添加一个函数，它提供图像的三维视图，并能够旋转它。引用[1] P. K.阿加瓦尔，T.哈格鲁普河Ray，M. Sharir，M. Smid和E.欢迎平移平面对象以最大化点包容：精确和近似算法。Mandarin pt，2001.[2] P. K. Agarwal和S.苏瑞曲面近似和几何分区。SIAM J. Comput. ，27：1016兰格、雷伊、斯米德和温特145[3] N. M. Amato，M. T. Goodrich和E. A.拉莫斯 计算安排 曲线段：通过采样的分治算法。在proc 第11届ACM-SIAM研讨会离散算法，第705-706页，2000年。[4] B. Chazelle和D. T.李你关于一个圆的放置问题。Computing，36：1-16，1986.[5] M. de Berg，M. van Kreveld，M. Overmars和O.施瓦茨科普夫计算几何：算法与应用。Springer-Verlag，柏林，1997年。[6] A. Gajentaan和M. H.奥弗马斯关于计算几何中的一类O（n2）问题。Comput.几何理论应用，5：165[7] J. Holm，K. de Lichtenberg和M.索鲁普连通性、最小生成树、2-边和双连通性的多对数确定性全动态算法在Proc.30th Annu. ACM研讨会。理论计算第79- 89页[8] K. 梅尔霍尔和S。 不是她。 LEDA：一个组合计算和几何计算的平台.剑桥大学出版社，英国剑桥，1999年[9] M. 索鲁普近似最优的全动态图连通性。在Proc.32ndAnnu. ACM研讨会。理论计算，2000年。