scipy.spatial.Delaunay: Finding points inside 3d poly yields points outside the 3d poly

Question

I want to find N random points which are located within a 3d poly e.g. a square-based pyramid. As speed will become a critical parameter later on I'd like to use scipy.spatial.Delaunay as proposed in Finding if point is in 3D poly in python.

Here is how I implemented it for N=10:

import numpy as np
import matplotlib.pyplot as plt
from scipy.spatial import Delaunay
from mpl_toolkits.mplot3d.art3d import Poly3DCollection

# -- SET VERTICES OF PYRAMID
a, h = 6,12
L = [h, a, a]
v = np.array([[-1*a/2,-1*a/2,0], [1*a/2,-1*a/2,0], [1*a/2,1*a/2,0], [-1*a/2,1*a/2,0], [0,0,1*h/2]]) # Pyramid base is centered in (0,0,0)

# --  PLOT PYRAMID
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.scatter3D(v[:, 0], v[:, 1], v[:, 2], marker='.', alpha = 1, c='k')
verts = [[v[0],v[1],v[4]], 
         [v[0],v[3],v[4]], 
         [v[2],v[1],v[4]],
         [v[2],v[3],v[4]]]
ax.add_collection3d(Poly3DCollection(verts, 
     facecolors='grey', linewidths=1, edgecolors='k', alpha=.25))

# -- FIND N RANDOM POINTS WITHIN THE 3D-PYRAMID
N = 10
pos = []
while len(pos) < N:
    test = False
    while test == False:
        testPos = L*(2*np.random.rand(3)-1)        
        test = Delaunay(v).find_simplex(testPos) >= 0  # True if point lies within poly

    pos.append(testPos)
    ax.scatter3D(testPos[1],testPos[2],testPos[0], marker='.', c='r')

However, when executing this code some of the points are always outside of the pyramid, as depicted in this examplary figure.

Is there a problem with the Delaunay triangulation or the pyramid? What am I doing wrong?

Answer 1

Apparently, I found the answer to my question:

I rearanged the vector L such that

L = [a, a, h]

and now it works. It somehow seems logical but to be honest, I don't know the real reason why it works like this and not the other way around. If somebody more experienced can explain it, I'd be interested!