SciPy spatial Delaunay/ConvexHull confusion

I am trying to generate random convex polyhedra. I generate a set of random 3D coordinates and then find their convex hull (so far so good).

I then thought I'd use a Delaunay triangulation to give me a triangulation of the convex hulls. This is where my basic understanding started to show!

Here's the code

import numpy as np
from scipy.spatial import ConvexHull
from scipy.spatial import Delaunay
import matplotlib as mpl
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# Generate random points & convex hull
points = np.random.rand(20,3)
hull = ConvexHull(points)

fig = plt.figure()
ax = fig.gca(projection = '3d')

# Plot hull's vertices
for vert in hull.vertices:    
    ax.scatter(points[vert,0], points[vert,1], zs=points[vert,2])#, 'ro')

# Calculate Delaunay triangulation & plot
tri = Delaunay(points[hull.vertices])
for simplex in tri.simplices:
    vert1 = [points[simplex[0],0], points[simplex[0],1], points[simplex[0],2]]
    vert2 = [points[simplex[1],0], points[simplex[1],1], points[simplex[1],2]]
    vert3 = [points[simplex[2],0], points[simplex[2],1], points[simplex[2],2]]
    vert4 = [points[simplex[3],0], points[simplex[3],1], points[simplex[3],2]]
    ax.plot([vert1[0], vert2[0]], [vert1[1], vert2[1]], zs = [vert1[2], vert2[2]])
    ax.plot([vert2[0], vert3[0]], [vert2[1], vert3[1]], zs = [vert2[2], vert3[2]])
    ax.plot([vert3[0], vert4[0]], [vert3[1], vert4[1]], zs = [vert3[2], vert4[2]])
    ax.plot([vert4[0], vert1[0]], [vert4[1], vert1[1]], zs = [vert4[2], vert1[2]])

plt.show()

A couple of things are concerning me, the plot sometimes misses out points on the hull & this seems to be Delaunay tetrahedralisation, which I suppose I shouldn't be that surprised about buts not what I'm after.

I would like a triangulation of just the hull surface,so I guess a simplex containing surface facets? Is this possible?

Thanks

EDIT: After pv's revelatory post below I have amended the code as follows;

import numpy as np
import pylab as pl
import scipy as sp
from scipy.spatial import ConvexHull
from scipy.spatial.distance import euclidean
import matplotlib.pyplot as plt
import mpl_toolkits.mplot3d as a3

aspect = 0
while aspect == 0:

    # Generate random points & convex hull
    points = np.random.rand(20,3)
    hull = ConvexHull(points)

    # Check aspect ratios of surface facets
    aspectRatio = []
    for simplex in hull.simplices:
        a = euclidean(points[simplex[0],:], points[simplex[1],:])
        b = euclidean(points[simplex[1],:], points[simplex[2],:])
        c = euclidean(points[simplex[2],:], points[simplex[0],:])
        circRad = (a*b*c)/(np.sqrt((a+b+c)*(b+c-a)*(c+a-b)*(a+b-c)))
        inRad = 0.5*np.sqrt(((b+c-a)*(c+a-b)*(a+b-c))/(a+b+c))
        aspectRatio.append(inRad/circRad)

    # Threshold for minium allowable aspect raio of surface facets
    if np.amin(aspectRatio) > 0.3:
        aspect = 1

ax = a3.Axes3D(pl.figure())
facetCol = sp.rand(3) #[0.0, 1.0, 0.0]

# Plot hull's vertices
#for vert in hull.vertices:    
#    ax.scatter(points[vert,0], points[vert,1], zs=points[vert,2])

# Plot surface traingulation
for simplex in hull.simplices:
    vtx = [points[simplex[0],:], points[simplex[1],:], points[simplex[2],:]]
    tri = a3.art3d.Poly3DCollection([vtx], linewidths = 2, alpha = 0.8)
    tri.set_color(facetCol)
    tri.set_edgecolor('k')
    ax.add_collection3d(tri)

plt.axis('off')
plt.show()

Now everything's working as I hoped. I added in an aspect ratio threshold to ensure a nicer triangulation.

Solution

Some things:

You give points[hull.vertices] as an argument to Delaunay, so the integers in tri.simplices are indices into points[hull.vertices], not into points, so that you end up plotting the wrong points
Tetrahedra have 6 ridges, but you are only plotting 4
If you need just the triangulation of the convex hull surface, that is available as hull.simplices

I.e.,

for simplex in hull.simplices:
    xs, ys, zs = points[simplex].T
    xs = np.r_[xs, xs[0]] # close polygons
    ys = np.r_[ys, ys[0]]
    zs = np.r_[zs, zs[0]]
    ax.plot(xs, ys, zs)

Or just:

ax.plot_trisurf(points[:,0], points[:,1], points[:,2],
                triangles=hull.simplices)