title: Difference between Matlab delaunayn and Scipy Delaunay
body: <p>I'm trying to replicate an N dimensional Delaunay triangulation that is performed by the Matlab delaunayn function in Python using the scipy.spatial.Delaunay function. However, while the Matlab function gives me the result I want and expect, scipy is giving me something different. I find this odd considering both are wrappers of the QHull library. I assume Matlab is implicitly setting different parameters in its call. The situation I'm trying to replicate between the two of them is found in <a href="http://www.mathworks.com/help/matlab/ref/delaunayn.html?requestedDomain=www.mathworks.com#examples" rel="nofollow noreferrer">Matlab's documentation</a>.</p>
  
  <p>The set up is to have a cube with a point in the center as below. The blue lines I provided to help visualize the shape, but they serve no purpose or meaning for this problem.</p>
  
  <p><a href="https://i.stack.imgur.com/HHByK.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/HHByK.png" alt="A cube with a point in the center"></a></p>
  
  <p>The triangulation I expect from this results in 12 simplices (listed in the Matlab example) and looks like the following.</p>
  
  <p><a href="https://i.stack.imgur.com/uUSMk.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/uUSMk.png" alt="Matlab&#39;s triangulation"></a></p>
  
  <p>However this python equivalent produces "extra" simplices.</p>
  
  <pre><code>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
  </code></pre>
  
  <p>The returned variable <code>simp</code> should be an M x N array where M is the number of simplices 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 N is 4.</p>
  
  <p>What I'm finding though is that M is actually 18 and that the extra 6 simplices are not tetrahedrons, but rather the 6 faces of the cube.</p>
  
  <p>What's going on here? How can I limit the returned simplices to only be tetrahedrons? I used this simple case to demonstrate the problem so I'd like a solution that isn't tailored to this problem.</p>
  
  <h2>EDIT</h2>
  
  <p>Thanks to an answer by Amro, I was able to figure this out and I can get a match in simplices between Matlab and Scipy. There were two factors in play. First, as pointed out, Matlab and Scipy use different QHull options. Second, QHull returns simplices with zero volume. Matlab removes these, Scipy doesn't. That was obvious in the example above because all 6 extra simplices were the zero-volume coplanar faces of the cube. These can be removed, in N dimensions, with the following bit of code.</p>
  
  <pre><code>N = 3 # The dimensions of our points
  options = 'Qt Qbb Qc' if N &lt;= 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)))) &lt; 1E-15:
          keep[i] = False # Point is coplanar, we don't want to keep it
  tri = tri[keep]
  </code></pre>
  
  <p>I suppose the other conditions should be addressed, but I'm guaranteed that my points contain no duplicates already, and the orientation condition appears to have no affect on the outputs that I can discern.</p>