using scipy.spatial.delaunay in place of matplotlib.tri.triangulation's built-in version

So you need to pass both the points and the triangles from qhull to the Triangulation constructor:

import numpy as np
import scipy.spatial
import matplotlib
import math

import matplotlib.tri as mtri
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d
import Axes3D

# First create the x and y coordinates of the points.
n_angles = 20
n_radii = 10
min_radius = 0.15
radii = np.linspace(min_radius, 0.95, n_radii)
angles = np.linspace(0, 2 * math.pi, n_angles, endpoint = False)
angles = np.repeat(angles[..., np.newaxis], n_radii, axis = 1)
angles[: , 1::2] += math.pi / n_angles
x = (radii * np.cos(angles)).flatten()
y = (radii * np.sin(angles)).flatten()

# Create the Delaunay tessalation using scipy.spatial
pts = np.vstack([x, y]).T
tess = scipy.spatial.Delaunay(pts)

# Create the matplotlib Triangulation object
x = tess.points[: , 0]
y = tess.points[: , 1]
tri = tess.vertices # or tess.simplices depending on scipy version
triang = mtri.Triangulation(x = pts[: , 0], y = pts[: , 1], triangles = tri)

# Plotting
z = x * x + y * y
fig = plt.figure()
ax = fig.gca(projection = '3d')
ax.plot_trisurf(triang, z)
plt.show()

Suggestion : 2

Indices of coplanar points and the corresponding indices of the nearest facet and the nearest vertex. Coplanar points are input points which were not included in the triangulation due to numerical precision issues.,The coordinates for the first point are all positive, meaning it is indeed inside the triangle. The third point is on a vertex, hence its null third coordinate.,Unless you pass in the Qhull option “QJ”, Qhull does not guarantee that each input point appears as a vertex in the Delaunay triangulation. Omitted points are listed in the coplanar attribute.,Indices of the points forming the simplices in the triangulation. For 2-D, the points are oriented counterclockwise.

>>> points = np.array([
      [0, 0],
      [0, 1.1],
      [1, 0],
      [1, 1]
   ]) >>>
   from scipy.spatial
import Delaunay
   >>>
   tri = Delaunay(points)

>>>
import matplotlib.pyplot as plt >>>
   plt.triplot(points[: , 0], points[: , 1], tri.simplices) >>>
   plt.plot(points[: , 0], points[: , 1], 'o') >>>
   plt.show()

>>> tri.simplices
array([
   [2, 3, 0], # may vary[3, 1, 0]
], dtype = int32)

>>> points[tri.simplices]
array([
   [
      [1., 0.], # may vary[1., 1.],
      [0., 0.]
   ],
   [
      [1., 1.],
      [0., 1.1],
      [0., 0.]
   ]
])

>>> tri.neighbors[1]
array([-1, 0, -1], dtype = int32) >>>
   points[tri.simplices[1, 1]]
array([0., 1.1])

>>> p = np.array([(0.1, 0.2), (1.5, 0.5), (0.5, 1.05)]) >>>
   tri.find_simplex(p)
array([1, -1, 1], dtype = int32)

Suggestion : 3

Catalan Number: the number of triangulations of a convex polygon with n + 2 vertices,The flip graph of any point set in the plane is connected (you can convert between different triangulations with edge flipping), which was proven by Charles Lawson in 1971.,Convex Hull algorithms are fundamental to computational geometry and it is useful to have a working knowledge of the terminology used to describe their performance.,Choose an interior point and draw edges to the three vertices of the triangle that contains it

1._

In[1]:

   %
   load_ext autoreload %
   autoreload 2

import itertools
import circumcircle
import triangle
import numpy as np
import scipy
import scipy.spatial
import scipy.optimize
import matplotlib
import matplotlib.pyplot as plt
from matplotlib.patches
import Polygon
   %
   matplotlib inline
import math
from ipywidgets
import interact, interactive, fixed
import ipywidgets as widgets
import cPickle as pickle
import shapefile #Python pyshp library
from mpl_toolkits.mplot3d
import Axes3D
from mpl_toolkits.mplot3d.art3d
import Poly3DCollection
from matplotlib.collections
import PolyCollection
from matplotlib.collections
import PatchCollection
from matplotlib
import colors
import time

2._

% load_ext autoreload
   %
   autoreload 2

import itertools
import circumcircle
import triangle
import numpy as np
import scipy
import scipy.spatial
import scipy.optimize
import matplotlib
import matplotlib.pyplot as plt
from matplotlib.patches
import Polygon
   %
   matplotlib inline
import math
from ipywidgets
import interact, interactive, fixed
import ipywidgets as widgets
import cPickle as pickle
import shapefile #Python pyshp library
from mpl_toolkits.mplot3d
import Axes3D
from mpl_toolkits.mplot3d.art3d
import Poly3DCollection
from matplotlib.collections
import PolyCollection
from matplotlib.collections
import PatchCollection
from matplotlib
import colors
import time

/Users/treddy / anaconda / lib / python2 .7 / site - packages / matplotlib / font_manager.py: 273: UserWarning: Matplotlib is building the font cache using fc - list.This may take a moment.
warnings.warn('Matplotlib is building the font cache using fc-list. This may take a moment.')

/Users/treddy / anaconda / lib / python2 .7 / site - packages / matplotlib / font_manager.py: 273: UserWarning: Matplotlib is building the font cache using fc - list.This may take a moment.
warnings.warn('Matplotlib is building the font cache using fc-list. This may take a moment.')

1._

In[2]:

   #plot an example of a polygon and an object that fails to be a polygon

fig_2_1 = matplotlib.pyplot.figure()

polygon_vertices = np.array([
   [0, 2], #top left[2, 2], #top right[2, 0], #bottom right[0, 0]
]) #bottom left

self_intersection_vertices = np.array([
   [0, 2], #top left[2, 2], #top right[0, 0], #bottom left[2, 0]
]) #bottom right

ax = fig_2_1.add_subplot(121)
ax.scatter(polygon_vertices[..., 0], polygon_vertices[..., 1], color = 'black', s = 200, label = 'vertices')
polygon = Polygon(polygon_vertices, alpha = 0.5)
ax.add_patch(polygon)
ax.set_title('A Polygon', fontsize = 18)

ax2 = fig_2_1.add_subplot(122)
ax2.scatter(self_intersection_vertices[..., 0], self_intersection_vertices[..., 1], color = 'black', s = 200, label = 'vertices')
polygon = Polygon(self_intersection_vertices, alpha = 0.5)
ax2.add_patch(polygon)
ax2.set_title('Self-Intersection = Not a Polygon', fontsize = 18)

for axis in [ax, ax2]:
   axis.set_xlabel('x', fontsize = 16)
axis.set_ylabel('y', fontsize = 16)

#some text labels
for the polygon
case
ax.text(0.85, 2.1, 'Edge', fontsize = 16)
ax.text(0.85, -0.2, 'Edge', fontsize = 16)
ax.text(-0.25, 1.0, 'Edge', fontsize = 16, rotation = 90)
ax.text(2.1, 1.0, 'Edge', fontsize = 16, rotation = 90)
ax.text(-0.4, 2.3, 'Vertex', fontsize = 16, rotation = -45)
ax.text(2.0, 2.3, 'Vertex', fontsize = 16, rotation = 45)
ax.text(-0.4, -0.1, 'Vertex', fontsize = 16, rotation = 45)
ax.text(2.0, -0.1, 'Vertex', fontsize = 16, rotation = -45)

fig_2_1.set_size_inches(17, 8.5)

2._

#plot an example of a polygon and an object that fails to be a polygon

fig_2_1 = matplotlib.pyplot.figure()

polygon_vertices = np.array([
   [0, 2], #top left[2, 2], #top right[2, 0], #bottom right[0, 0]
]) #bottom left

self_intersection_vertices = np.array([
   [0, 2], #top left[2, 2], #top right[0, 0], #bottom left[2, 0]
]) #bottom right

ax = fig_2_1.add_subplot(121)
ax.scatter(polygon_vertices[..., 0], polygon_vertices[..., 1], color = 'black', s = 200, label = 'vertices')
polygon = Polygon(polygon_vertices, alpha = 0.5)
ax.add_patch(polygon)
ax.set_title('A Polygon', fontsize = 18)

ax2 = fig_2_1.add_subplot(122)
ax2.scatter(self_intersection_vertices[..., 0], self_intersection_vertices[..., 1], color = 'black', s = 200, label = 'vertices')
polygon = Polygon(self_intersection_vertices, alpha = 0.5)
ax2.add_patch(polygon)
ax2.set_title('Self-Intersection = Not a Polygon', fontsize = 18)

for axis in [ax, ax2]:
   axis.set_xlabel('x', fontsize = 16)
axis.set_ylabel('y', fontsize = 16)

#some text labels
for the polygon
case
ax.text(0.85, 2.1, 'Edge', fontsize = 16)
ax.text(0.85, -0.2, 'Edge', fontsize = 16)
ax.text(-0.25, 1.0, 'Edge', fontsize = 16, rotation = 90)
ax.text(2.1, 1.0, 'Edge', fontsize = 16, rotation = 90)
ax.text(-0.4, 2.3, 'Vertex', fontsize = 16, rotation = -45)
ax.text(2.0, 2.3, 'Vertex', fontsize = 16, rotation = 45)
ax.text(-0.4, -0.1, 'Vertex', fontsize = 16, rotation = 45)
ax.text(2.0, -0.1, 'Vertex', fontsize = 16, rotation = -45)

fig_2_1.set_size_inches(17, 8.5)

1._

In[3]:

   @interact(y_value_p1 = widgets.FloatSlider(min = -1, max = 2.1, step = 0.1, value = -1),
      y_value_p2 = widgets.FloatSlider(min = -1, max = 2.1, step = 0.1, value = -1),
      new_diagonal_p2 = False)
def polygon_sweep(y_value_p1, y_value_p2, new_diagonal_p2):
   fig_2_2 = matplotlib.pyplot.figure()
fig_2_2.set_size_inches(17, 8.5)

polygon_1 = np.array([
   [0, 2], #top[-1, 1], #left[0, -1], #bottom[1, 1]
]) #right

polygon_2 = np.array([
   [0, 0.7], #top[-1, 2], #left[0, -1], #bottom[1, 2]
]) #right

ax1 = fig_2_2.add_subplot(121)
p1 = Polygon(polygon_1, color = 'blue', alpha = 0.2)
ax1.add_patch(p1)
ax1.scatter(polygon_1[..., 0], polygon_1[..., 1], s = 80)
ax2 = fig_2_2.add_subplot(122)
p2 = Polygon(polygon_2, color = 'blue', alpha = 0.2)
ax2.add_patch(p2)
ax2.scatter(polygon_2[..., 0], polygon_2[..., 1], s = 80)

for axis in [ax1, ax2]:
   axis.set_xlim(-1.2, 1.2)
axis.set_ylim(-1.1, 2.1)
ax1.axhline(y = y_value_p1, lw = 4, color = 'black')
ax2.axhline(y = y_value_p2, lw = 4, color = 'black')
if new_diagonal_p2:
   ax2.axvline(x = 0, ymin = 0.03, ymax = 0.55, color = 'green', lw = 4)

2._

@interact(y_value_p1 = widgets.FloatSlider(min = -1, max = 2.1, step = 0.1, value = -1),
   y_value_p2 = widgets.FloatSlider(min = -1, max = 2.1, step = 0.1, value = -1),
   new_diagonal_p2 = False)
def polygon_sweep(y_value_p1, y_value_p2, new_diagonal_p2):
   fig_2_2 = matplotlib.pyplot.figure()
fig_2_2.set_size_inches(17, 8.5)

polygon_1 = np.array([
   [0, 2], #top[-1, 1], #left[0, -1], #bottom[1, 1]
]) #right

polygon_2 = np.array([
   [0, 0.7], #top[-1, 2], #left[0, -1], #bottom[1, 2]
]) #right

ax1 = fig_2_2.add_subplot(121)
p1 = Polygon(polygon_1, color = 'blue', alpha = 0.2)
ax1.add_patch(p1)
ax1.scatter(polygon_1[..., 0], polygon_1[..., 1], s = 80)
ax2 = fig_2_2.add_subplot(122)
p2 = Polygon(polygon_2, color = 'blue', alpha = 0.2)
ax2.add_patch(p2)
ax2.scatter(polygon_2[..., 0], polygon_2[..., 1], s = 80)

for axis in [ax1, ax2]:
   axis.set_xlim(-1.2, 1.2)
axis.set_ylim(-1.1, 2.1)
ax1.axhline(y = y_value_p1, lw = 4, color = 'black')
ax2.axhline(y = y_value_p2, lw = 4, color = 'black')
if new_diagonal_p2:
   ax2.axvline(x = 0, ymin = 0.03, ymax = 0.55, color = 'green', lw = 4)

Suggestion : 4

If the planar region $D$ ($U$) is rectangular, then one defines a meshgrid on it, and the points of the grid are the input points for the scipy.spatial.Delaunay function that defines the planar triangulation of $D$, respectively $U$.,The Delaunay triangulation of a planar region is defined and illustrated in a Python Plotly tutorial posted here.,The images of the points2D through the surface parameterization are 3D points. The same simplices define the triangles on the surface.,We consider polar coordinates on the disk, $D(0, 1)$, centered at origin and of radius 1, and define a meshgrid on the set of points $(r, \theta)$, with $r\in[0,1]$ and $\theta\in[0,2\pi]$:

import plotly.plotly as py
import plotly.graph_objs as go

import numpy as np
import matplotlib.cm as cm
from scipy.spatial
import Delaunay

u = np.linspace(0, 2 * np.pi, 24)
v = np.linspace(-1, 1, 8)
u, v = np.meshgrid(u, v)
u = u.flatten()
v = v.flatten()

#evaluate the parameterization at the flattened u and v
tp = 1 + 0.5 * v * np.cos(u / 2.)
x = tp * np.cos(u)
y = tp * np.sin(u)
z = 0.5 * v * np.sin(u / 2.)

#define 2 D points, as input data
for the Delaunay triangulation of U
points2D = np.vstack([u, v]).T
tri = Delaunay(points2D) #triangulate the rectangle U

print tri.simplices.shape, '\n', tri.simplices[0]

(322, 3)[73 96 72]

def map_z2color(zval, colormap, vmin, vmax):
   #map the normalized value zval to a corresponding color in the colormap

if vmin > vmax:
   raise ValueError('incorrect relation between vmin and vmax')
t = (zval - vmin) / float((vmax - vmin)) #normalize val
R, G, B, alpha = colormap(t)
return 'rgb(' + '{:d}'.format(int(R * 255 + 0.5)) + ',' + '{:d}'.format(int(G * 255 + 0.5)) + \
   ',' + '{:d}'.format(int(B * 255 + 0.5)) + ')'

def tri_indices(simplices):
   #simplices is a numpy array defining the simplices of the triangularization
#returns the lists of indices i, j, k

return ([triplet[c]
      for triplet in simplices
   ]
   for c in range(3))

def plotly_trisurf(x, y, z, simplices, colormap = cm.RdBu, plot_edges = None):
   #x, y, z are lists of coordinates of the triangle vertices
#simplices are the simplices that define the triangularization;
#simplices is a numpy array of shape(no_triangles, 3)
#insert here the type check
for input data

points3D = np.vstack((x, y, z)).T
tri_vertices = map(lambda index: points3D[index], simplices) # vertices of the surface triangles
zmean = [np.mean(tri[: , 2]) for tri in tri_vertices] # mean values of z - coordinates of
   #triangle vertices
min_zmean = np.min(zmean)
max_zmean = np.max(zmean)
facecolor = [map_z2color(zz, colormap, min_zmean, max_zmean) for zz in zmean]
I, J, K = tri_indices(simplices)

triangles = go.Mesh3d(x = x,
   y = y,
   z = z,
   facecolor = facecolor,
   i = I,
   j = J,
   k = K,
   name = ''
)

if plot_edges is None: # the triangle sides are not plotted
return [triangles]
else:
   #define the lists Xe, Ye, Ze, of x, y, resp z coordinates of edge end points
for each triangle
#None separates data corresponding to two consecutive triangles
lists_coord = [
   [
      [T[k % 3][c]
         for k in range(4)
      ] + [None]
      for T in tri_vertices
   ]
   for c in range(3)
]
Xe, Ye, Ze = [reduce(lambda x, y: x + y, lists_coord[k]) for k in range(3)]

#define the lines to be plotted
lines = go.Scatter3d(x = Xe,
   y = Ye,
   z = Ze,
   mode = 'lines',
   line = dict(color = 'rgb(50,50,50)', width = 1.5)
)
return [triangles, lines]

data1 = plotly_trisurf(x, y, z, tri.simplices, colormap = cm.RdBu, plot_edges = True)

Suggestion : 5

Here are the examples of the python api scipy.spatial.Delaunay taken from open source projects. By voting up you can indicate which examples are most useful and appropriate.

def triangulate(vertices):
   ''
'Create the triangulation of a set of points, and output the resulting triangles'
''
from scipy.spatial
import Delaunay
d = Delaunay(vertices)
polys = d.points[d.vertices].tolist()
genericpolys = [GenericPoly.gen_from_point_lists(poly, []) for poly in polys]
return genericpolys

    @dec.skipif(not has_matplotlib, "Matplotlib not available")
    def test_delaunay(self):
       # Smoke test
    fig = plt.figure()
    obj = Delaunay(self.points)
    s_before = obj.simplices.copy()
    r = delaunay_plot_2d(obj, ax = fig.gca())
    assert_array_equal(obj.simplices, s_before) # shouldn 't modify
    assert_(r is fig)
    delaunay_plot_2d(obj, ax = fig.gca())

def triangulatePoints(points):
   points = toNumpy(points)
triangulation = sp.Delaunay(points)
triangles = []
for i in range(len(triangulation.simplices)):
   verts = map(lambda p: shapes.Point(p[0], p[1]),
      points[triangulation.simplices[i,: ]])
triangle = shapes.Triangle(
   verts[0], verts[1], verts[2])
triangles.append(triangle)
return triangles

def inside_convex_poly(pts):
   ""
"Returns a function that checks if inputs are inside the convex hull of polyhedron defined by pts

Alternative method to check is to get faces of the convex hull, then check
if each normal is pointed away from each point.
As it turns out, this is vastly slower than using qhull 's find_simplex, even though the simplex is not needed.
""
"
tri = Delaunay(pts)
return lambda x: tri.find_simplex(x) != -1

    def get_data(self, element, ranges):
       try:
       from scipy.spatial
    import Delaunay
    except:
       SkipRendering("SciPy not available, cannot plot Trisurface")
    x, y, z = (element.dimension_values(i) for i in range(3))
    points2D = np.vstack([x, y]).T
    tri = Delaunay(points2D)
    simplices = tri.simplices
    return (x, y, z, simplices, self.colorbar, 'black'), {}