The difference between Matlab delaunayn and Scipy Delaunay

Question

I am trying to replicate the three-dimensional Delaunay triangulation performed by the Matlab delaunayn function in Python using the scipy.spatial.Delaunay function. However, although the Matlab function gives me the result that I want and expect, scipy gives me something else. I find this strange given that both are wrappers in the QHull library. I assume that Matlab implicitly sets different parameters in its call. The situation that I am trying to reproduce between the two of them is in the Matlab documentation .

The setup is to have a cube with a dot in the center, as shown below. The blue lines that I provided to render the figure, but they have no purpose or meaning for this problem.

The triangulation that I expect from this leads to 12 simplexes (indicated in the Matlab example) and looks like this.

However, this python equivalent creates superfluous simplexes.

x = np.array([[-1,-1,-1],[-1,-1,1],[-1,1,-1],[1,-1,-1],[1,1,1],[1,1,-1],[1,-1,1],[-1,1,1],[0,0,0]]) simp = scipy.spatial.Delaunay(x).simplices 

The returned variable simp should be an array M x N, where M is the number of simplexes found (should be 12 for my case), and N is the number of points in the simplex. In this case, each simplex should be a tetrahedron, meaning that N is 4.

What I find is that M is actually 18 and that the additional 6 simplexes are not tetrahedrons, but rather 6 faces of the cube.

What's going on here? How can I limit returned simplexes to only tetrahedra? I used this simple case to demonstrate the problem, so I would like the solution not to be adapted to this problem.

EDIT

Thanks to Amro's answer, I was able to figure this out, and I can get a match in sympathies between Matlab and Scipy. There were two factors in the game. First, as stated, Matlab and Scipy use different QHull options. Secondly, QHull returns simplifications with zero volume. Matlab removes them, Scipy does not. This was obvious in the example above, because all 6 additional simplexes were coplanar faces of the zero volume of the cube. They can be deleted in N dimensions with the next bit of code.

 N = 3 # The dimensions of our points options = 'Qt Qbb Qc' if N <= 3 else 'Qt Qbb Qc Qx' # Set the QHull options tri = scipy.spatial.Delaunay(points, qhull_options = options).simplices keep = np.ones(len(tri), dtype = bool) for i, t in enumerate(tri): if abs(np.linalg.det(np.hstack((points[t], np.ones([1,N+1]).T)))) < 1E-15: keep[i] = False # Point is coplanar, we don't want to keep it tri = tri[keep] 

I believe that other conditions should be resolved, but I am guaranteed that my points do not contain duplicates, and the orientation condition does not affect the output that I can distinguish.

Answer 1

Some notes comparing MATLAB and SciPy:

• According to docs, MATLAB by default uses Qt Qbb Qc Qhull options for 3D input, while SciPy uses Qt Qbb Qc Qz .

• not sure if that matters, but your NumPy array is not in the same order as the points created with ndgrid in MATLAB.

In fact, if you look at the MATLAB code in edit delaunayn.m , you will see three additional steps:

• first it merges the repeating points of mergeDuplicatePoints (this is not a problem in your case)
• then it provides alignment of orientation for points (see code)

• finally, getting the result from Qhull (implemented as the qhullmx MEX function), there is the following comment above a few lines of code:

  Separate the simplified zero-volume styles that might have been created due to degeneracy.

Since the file is copyrighted, I will not post the code here, but you can check it at your end.