3D Delaunay triangulation: bad output (extra simplices appearing)

Question

I am using python3.11 to create the Delaunay triangulation of a point cloud with script.Delaunay and it is misbehaving by creating some extra faces. In the image below you can see a 3D scatter plot of the points.

The image was created using the next very few lines of code:

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()

The data puntos can be downloaded as a csv file in this link. Now, as I said, I am interested in obtaining the Delaunay triangulation of that point could, for which the following piece of code is used.

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()

This produces the following image (point cloud images and triangulation image do not have exactly same aspect ratio, but I find it sufficient this way):

As you might see, the triangulation is creating some extra faces in the boundary of the surface, and that is repeated along the four sides of the boundary. Anyone knows why this happens and how can I solve it?

Thank you in advance!

Answer 1

Referring to my comments for your original question... The extra simplices occur because some of your vertices are positioned inside, but near, the convex hull of the Delaunay triangulation. The code below adjusts the position of the problematic vertices by finding the nearest point on the convex hull. For this example, I use the Tinfour Software Library which is written in Java. But you should be able to adapt the ideas to Python if you wish to do so.

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());
    }
  }

}

And here are the modified vertices. Some of the z values may be a little different because Tinfour only supports single-precision floating point values for its z elements.

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