3D Delaunay 三角剖分：输出错误（出现额外的单纯形）

Question

我正在使用 python3.11 通过 script.Delaunay 创建点云的 Delaunay 三角剖分，它通过创建一些额外的面而表现不佳。在下图中，您可以看到点的 3D 散点图。

该图像是使用接下来的几行代码创建的：

import plotly.graph_objects as go
fig = go.Figure()
fig.add_scatter3d(x = puntos[:,0], y = puntos[:,1], z = puntos[:,2],mode='markers', marker=dict(
                size=1,
                color='rgb(0,0,0)',
                opacity=0.8
            ))
fig.update_layout(scene = dict(aspectmode = 'data'))

fig.show()

数据蓬托可以下载为 csv 文件in this link。现在，正如我所说，我有兴趣获得该点的 Delaunay 三角剖分，为此使用了以下代码段。

import numpy as np
import pandas as pd
from scipy.spatial import Delaunay
import plotly.figure_factory as ff

puntos = pd.read_csv('puntos.csv')
puntos = puntos[['0', '1', '2']]
tri = Delaunay(np.array([puntos[:,0], puntos[:,1]]).T)
simplices = tri.simplices
fig = ff.create_trisurf(x=puntos[:,0], y=puntos[:,1], z=puntos[:,2],
                         simplices=simplices, aspectratio=dict(x=1, y=1, z=0.3))
fig.show()

这会产生以下图像（点云图像和三角测量图像的纵横比不完全相同，但我觉得这样就足够了）：

正如您所看到的，三角剖分在曲面的边界中创建了一些额外的面，并且沿着边界的四个边重复进行。任何人都知道为什么会发生这种情况，我该如何解决？

先感谢您！

Answer 1

首先，我祝贺你有这样一个精美的例子。我建议您探索函数 Delaunay 的示例。 documentation 展示了您可能感兴趣的几个属性。

Answer 2

关于您的原始问题，请参考我的 cmets ... 出现额外的单纯形是因为您的某些顶点位于 Delaunay 三角剖分的凸包内部，但靠近凸包。下面的代码通过在凸包上找到最近的点来调整有问题的顶点的位置。对于此示例，我使用用 Java 编写的 Tinfour Software Library。但是如果你愿意的话，你应该能够将这些想法应用到 Python 中。

import java.io.File;
import java.io.IOException;
import java.util.ArrayList;
import java.util.Comparator;
import java.util.List;
import org.tinfour.common.IIncrementalTin;
import org.tinfour.common.IQuadEdge;
import org.tinfour.common.Vertex;
import org.tinfour.semivirtual.SemiVirtualIncrementalTin;
import org.tinfour.utils.loaders.VertexReaderText;

public class AdjustEdgePoints {
  public static void main(String[] args) throws IOException {
    File input = new File("puntos.csv");
    List<Vertex> vertices = null;
    try ( VertexReaderText vrt = new VertexReaderText(input)) {
      vertices = vrt.read(null);
    }

    double pointSpacing = 0.2;
    double edgeTooLong = 0.4;  // based on point spacng
    double pointTooClose = 0.021;

    // vertices are numbered 0 to n-1
    boolean[] doNotTest = new boolean[vertices.size()];
    boolean[] modified = new boolean[vertices.size()];
    List<Vertex> replacements = new ArrayList<>();

    IIncrementalTin tin = new SemiVirtualIncrementalTin(pointSpacing);
    tin.add(vertices, null);
    List<IQuadEdge> perimeter = tin.getPerimeter();  // the convex hull
    // mark all vertices on the perimeter as do-not-test
    for (IQuadEdge edge : perimeter) {
      Vertex A = edge.getA(); // vertices are for edge are A and B
      doNotTest[A.getIndex()] = true;
    }

    // For all excessively long edges, find the vertices that are too close
    // and move them to the edge.
    for (IQuadEdge edge : perimeter) {
      double eLength = edge.getLength();
      if (eLength < edgeTooLong) {
        continue; // no processing required
      }
      Vertex A = edge.getA();
      Vertex B = edge.getB();
      double eX = B.getX() - A.getX();  // vector in direction of edge
      double eY = B.getY() - A.getY();
      double pX = -eY / eLength; // unit vector perpendicular to edge
      double pY = eX / eLength;

      for (Vertex v : vertices) {
        if (doNotTest[v.getIndex()]) {
          continue;
        }
        double vX = v.getX() - A.getX();
        double vY = v.getY() - A.getY();
        // compute t, the parameter for a point on the line of the edge
        // closest to the vertex.  We are only interested in this point
        // if it falls between the two endpoints of the edge.
        // in that case, t will be in the range 0 < t < 1
        double t = (vX * eX + vY * eY) / (eLength * eLength);
        if (0 < t && t < 1) {
          double s = pX * vX + pY * vY; // distance of V from edge
          if (s < pointTooClose) {
            double x = A.getX() + t * eX; // point on edge
            double y = A.getY() + t * eY;
            Vertex X = new Vertex(x, y, v.getZ(), v.getIndex());
            modified[v.getIndex()] = true;
            doNotTest[v.getIndex()] = true;
            replacements.add(X);
          }
        }
      } // end of vertices loop
    } // end of perimeter loop

    System.out.println("i,x,y,z");
    for (Vertex v : vertices) {
      if (!modified[v.getIndex()]) {
        System.out.format("%d,%19.16f,%19.16f,%19.16f%n",
          v.getIndex(), v.getX(), v.getY(), v.getZ());
      }
    }

    replacements.sort(new Comparator<Vertex>() {
      @Override
      public int compare(Vertex arg0, Vertex arg1) {
        return Integer.compare(arg0.getIndex(), arg1.getIndex());
      }
    });
    System.out.println("");
    for (Vertex v : replacements) {
      System.out.format("%d,%19.16f,%19.16f,%19.16f%n",
        v.getIndex(), v.getX(), v.getY(), v.getZ());
    }
  }

}

这是修改后的顶点。某些 z 值可能略有不同，因为 Tinfo 仅支持其 z 元素的单精度浮点值。

1,-1.7362462292391640,-1.9574243190958676, 0.2769193053245544
2,-1.5328386393032090,-1.9691151455342581, 0.2039903998374939
3,-1.3356454135288550,-1.9804488021499242, 0.1507489830255508
7,-0.5545998965948742,-2.0099766379255533, 0.1641141176223755
8,-0.3452278908343442,-2.0128547506964680, 0.2061954736709595
9,-0.1280981557739663,-2.0158395043704496, 0.2454121708869934
10, 0.0960477228578226,-2.0189207048083720, 0.2699134647846222
11, 0.3249315805387950,-2.0220670354338774, 0.2698880732059479
12, 0.5546177976887393,-2.0252243956190710, 0.2401449084281921
13, 0.7809245195097682,-2.0283352998865780, 0.1818846762180328
14, 1.0012918968621154,-2.0313645595086890, 0.1022855937480927
15, 1.2151144438651920,-2.0343038512301570, 0.0123253259807825
20,-1.9574243190958676,-1.7362462292391640, 0.2769193053245544
40,-1.9691151455342581,-1.5328386393032090, 0.2039903998374939
60,-1.9804488021499242,-1.3356454135288550, 0.1507489830255508
99, 2.0343038512301570,-1.2151144438651922, 0.0123253259807825
119, 2.0313645595086890,-1.0012918968621158, 0.1022855937480927
139, 2.0283352998865780,-0.7809245195097688, 0.1818846762180328
140,-2.0099766379255533,-0.5545998965948735, 0.1641141176223755
159, 2.0252243956190710,-0.5546177976887400, 0.2401449084281921
160,-2.0128547506964680,-0.3452278908343438, 0.2061954736709595
179, 2.0220670354338774,-0.3249315805387953, 0.2698880732059479
180,-2.0158395043704496,-0.1280981557739660, 0.2454121708869934
199, 2.0189207048083720,-0.0960477228578232, 0.2699134647846222
200,-2.0189207048083720, 0.0960477228578229, 0.2699134647846222
219, 2.0158395043704496, 0.1280981557739658, 0.2454121708869934
220,-2.0220670354338774, 0.3249315805387953, 0.2698880732059479
239, 2.0128547506964680, 0.3452278908343438, 0.2061954736709595
240,-2.0252243956190710, 0.5546177976887398, 0.2401449084281921
259, 2.0099766379255533, 0.5545998965948733, 0.1641141176223755
260,-2.0283352998865780, 0.7809245195097683, 0.1818846762180328
280,-2.0313645595086890, 1.0012918968621158, 0.1022855937480927
300,-2.0343038512301570, 1.2151144438651922, 0.0123253259807825
339, 1.9804488021499242, 1.3356454135288547, 0.1507489830255508
359, 1.9691151455342581, 1.5328386393032085, 0.2039903998374939
379, 1.9574243190958676, 1.7362462292391640, 0.2769193053245544
384,-1.2151144438651920, 2.0343038512301570, 0.0123253259807825
385,-1.0012918968621156, 2.0313645595086890, 0.1022855937480927
386,-0.7809245195097685, 2.0283352998865780, 0.1818846762180328
387,-0.5546177976887396, 2.0252243956190710, 0.2401449084281921
388,-0.3249315805387951, 2.0220670354338774, 0.2698880732059479
389,-0.0960477228578228, 2.0189207048083720, 0.2699134647846222
390, 0.1280981557739660, 2.0158395043704496, 0.2454121708869934
391, 0.3452278908343442, 2.0128547506964680, 0.2061954736709595
392, 0.5545998965948740, 2.0099766379255533, 0.1641141176223755
396, 1.3356454135288547, 1.9804488021499242, 0.1507489830255508
397, 1.5328386393032085, 1.9691151455342581, 0.2039903998374939
398, 1.7362462292391640, 1.9574243190958676, 0.2769193053245544