Return to Answer

1. The first thing to say is that you have made this question unnecessarily difficult for us because your functions have no documentation. We have to read and reverse-engineer your code to try to figure out what the purpose of each function is. Python allows you to write a "docstring" for each function in which you explain what the function does, what arguments it takes, and what values it returns.

def containing_tet(tg, point):
 """If point is inside the i'th tetrahedron in tg, return the pair
 (i, bcoords) where bcoords are the barycentric coordinates of
 point within that tetrahedron. Otherwise, return (None, None).
 """

and so on for your other functions. Writing this kind of documentation will help your colleagues (and you in a few months when you have forgotten all the details).

1. Why is containing_tet slow? Well, you say that you are calling it many times, and each time it has to repeat some work. First, it collects the vertices of each tetrahedron:

    verts = [tg.points[j] for j in tet]
   

T = (np.array(vertices[:-1])-vertices[-1]).T
... la.inv(T) ...

1. You are calling barycentric_coords for many points. So you should vectorize this code: that is, read all the targets into a NumPy array, and then compute the barycentric coordinates for all the targets at once.

2. Having written the above, however, I am wondering why you did not use scipy.spatial.Delaunay ? If I had to write this code I would do the triangulations like this:

    >>> import numpy as np
    >>> import scipy.spatial
    >>> points = np.array([(0,0,0),(0,0,1),(0,1,0),(0,1,1),(1,0,0),(1,0,1),(1,1,0),(1,1,1)])
    >>> tri = scipy.spatial.Delaunay(points)
    >>> tri.simplices
    array([[3, 2, 4, 0],
    [3, 1, 4, 0],
    [3, 6, 2, 4],
    [3, 6, 7, 4],
    [3, 5, 1, 4],
    [3, 5, 7, 4]], dtype=int32)
   

>>> targets = np.array([[.1,.1,.1], [.9,.9,.9], [.1,.6,.7], [.4,.9,.1]])

>>> tetrahedra = tri.find_simplex(targets)
>>> tetrahedra
array([0, 3, 1, 2], dtype=int32)

>>> X = tri.transform[tetrahedra,:3]
>>> Y = targets - tri.transform[tetrahedra,3]
>>> b = np.einsum('ijk,ik->ij', X, Y)
>>> bcoords = np.c_[b, 1 - b.sum(axis=1)]
>>> bcoords
array([[ 0.1, 0. , 0.1, 0.8],
 [ 0.1, 0. , 0.8, 0.1],
 [ 0.6, 0.1, 0.1, 0.2],
 [ 0.1, 0.3, 0.5, 0.1]])

>>> b = np.array([x.dot(y) for x, y in zip(X, Y)])

Answers to questions

1. How to structure your code? I would expect you to carry out the following sequence of operations. For maximum speed, each operation should be vectorized, that is, it should compute a NumPy array of results using a handful of calls to NumPy functions, with no Python for loops. (The equivalent outputs in my §4 above, if any, are in parentheses.)

   1. Load the array of vertices to triangulate (points).
   2. Perform the Delaunay triangulation, getting an array of arrays of indices of vertices (tri.simplices).
   3. Find the vertices for each tetrahedron, getting an array of arrays of vertices (points[tri.simplices]).
   4. Compute the affine transform for each tetrahedron, getting an array of inverse matrices (tri.transform[:,:3]) and an array of offset vectors (tri.transform[:,3]).
   5. Load the array of target points (targets).
   6. Compute the barycentric coordinates of each target point in each tetrahedron. (No exact equivalent in my code: I used the output of find_simplex to compute the barycentric coordinates only for the tetrahedron it was found in, but the algorithm should be similar, with suitable adjustment of the indices in the np.einsum call.)
   7. Identify the coordinates that are non-negative and so find which tetrahedron each target was found in (tetrahedra).

2. How to read the names from the colorlist CSV? In NumPy each array must have a single data type (that's one the requirements for NumPy to be able to process them quickly). So you must read the names into one array and the points into another array. An easy way to do this is to call numpy.loadtxt twice:

    >>> names = np.loadtxt('XYZcolorlist_D65.csv', usecols=(0,), delimiter=',', dtype=str) # use dtype=bytes in Python 3
    >>> points = np.loadtxt('XYZcolorlist_D65.csv', usecols=(1,2,3), delimiter=',')
    >>> next(zip(names, points))
    ('255 63 127', array([ 35.53443021, 21.38072197, 20.3661096 ]))
   

1. Introduction

The first thing to say is that you have made this question unnecessarily difficult for us because your functions have no documentation. We have to read and reverse-engineer your code to try to figure out what the purpose of each function is. Python allows you to write a "docstring" for each function in which you explain what the function does, what arguments it takes, and what values it returns.

def containing_tet(tg, point):
 """If point is inside the i'th tetrahedron in tg, return the pair
 (i, bcoords) where bcoords are the barycentric coordinates of
 point within that tetrahedron. Otherwise, return (None, None).
 """

and so on for your other functions. Writing this kind of documentation will help your colleagues (and you in a few months when you have forgotten all the details).

2. Diagnosis

Why is containing_tet slow? Well, you say that you are calling it many times, and each time it has to repeat some work. First, it collects the vertices of each tetrahedron:

verts = [tg.points[j] for j in tet]

T = (np.array(vertices[:-1])-vertices[-1]).T
... la.inv(T) ...

Then you need to reorganize your code so that each operation is vectorized: that is, you should read all the targets into a NumPy array, and then compute the barycentric coordinates for all the targets at once.

3. Using scipy.spatial

Having written the above, however, I am wondering why you did not use scipy.spatial.Delaunay ? Is it because TetGen does a better job?

Since I don't have a copy of TetGen handy, if I had to write this code I would generate the triangulations like this:

>>> import numpy as np
>>> import scipy.spatial
>>> points = np.array([(0,0,0),(0,0,1),(0,1,0),(0,1,1),(1,0,0),(1,0,1),(1,1,0),(1,1,1)])
>>> tri = scipy.spatial.Delaunay(points)
>>> tri.simplices
array([[3, 2, 4, 0],
 [3, 1, 4, 0],
 [3, 6, 2, 4],
 [3, 6, 7, 4],
 [3, 5, 1, 4],
 [3, 5, 7, 4]], dtype=int32)

>>> targets = np.array([[.1,.1,.1], [.9,.9,.9], [.1,.6,.7], [.4,.9,.1]])

>>> tetrahedra = tri.find_simplex(targets)
>>> tetrahedra
array([0, 3, 1, 2], dtype=int32)

>>> X = tri.transform[tetrahedra,:3]
>>> Y = targets - tri.transform[tetrahedra,3]
>>> b = np.einsum('ijk,ik->ij', X, Y)
>>> bcoords = np.c_[b, 1 - b.sum(axis=1)]
>>> bcoords
array([[ 0.1, 0. , 0.1, 0.8],
 [ 0.1, 0. , 0.8, 0.1],
 [ 0.6, 0.1, 0.1, 0.2],
 [ 0.1, 0.3, 0.5, 0.1]])

>>> b = np.array([x.dot(y) for x, y in zip(X, Y)])

4. Answers to questions

1. How to structure your code? See §5 below. Notice how using NumPy efficiently requires us to turn the code "inside out". Instead of looping over the targets at top level, and in each loop iteration performing a series of operations on that one target, we lift the series of operations up to top level and make each operation run over all elements of a NumPy array.

2. How to read the colors from the colorlist CSV? In NumPy each array must have a single data type (that's one the requirements for NumPy to be able to process them quickly). So you must read the names into one array and the points into another array. An easy way to do this is to call numpy.loadtxt twice:

    >>> colors = np.loadtxt('XYZcolorlist_D65.csv', usecols=(0,), delimiter=',', 
    ... converters={0:lambda s:s.split()}, dtype=np.uint8)
    >>> colors[:5]
    array([[255, 63, 127],
    [255, 95, 127],
    [255, 127, 127],
    [255, 159, 127],
    [255, 223, 159]], dtype=uint8)
    >>> points = np.loadtxt('XYZcolorlist_D65.csv', usecols=(1,2,3), delimiter=',')
    >>> points[:5]
    array([[ 35.53443021, 21.38072197, 20.3661096 ],
    [ 40.20749455, 26.52829494, 22.70942844],
    [ 43.66474384, 32.34826255, 23.61818015],
    [ 47.12256284, 39.17809444, 22.9366615 ],
    [ 61.73791496, 62.83876017, 32.39362009]])
   

5. Revised code.

import numpy as np
import scipy.spatial
# Configuration.
POINTS_FILENAME = 'XYZcolorlist_D65.csv'
TARGETS_FILENAME = 'targets.csv'
OUTPUT_FILENAME = 'gamut.ppm'
DEFAULT_COLOR = np.array([[255, 255, 255]], dtype=np.uint8)
# Load colors
colors = np.loadtxt(POINTS_FILENAME, usecols=(0,), delimiter=',', 
 converters={0:lambda s:s.split()}, dtype=np.uint8)
# Load points
points = np.loadtxt(POINTS_FILENAME, usecols=(1, 2, 3), delimiter=',')
# Load targets
targets = np.loadtxt(TARGETS_FILENAME, delimiter=',')
ntargets = len(targets)
# Compute Delaunay triangulation of points.
tri = scipy.spatial.Delaunay(points)
# Find the tetrahedron containing each target (or -1 if not found)
tetrahedra = tri.find_simplex(targets)
# Affine transformation for tetrahedron containing each target
X = tri.transform[tetrahedra, :3]
# Offset of each target from the origin of its containing tetrahedron
Y = targets - tri.transform[tetrahedra, 3]
# First three barycentric coordinates of each target in its tetrahedron.
# The fourth coordinate would be 1 - b.sum(axis=1), but we don't need it.
b = np.einsum('...jk,...k->...j', X, Y)
# Cumulative sum of barycentric coordinates of each target.
bsum = np.c_[b.cumsum(axis=1), np.ones(ntargets)]
# A uniform random number in [0, 1] for each target.
R = np.random.uniform(0, 1, size=(ntargets, 1))
# Randomly choose one of the tetrahedron vertices for each target,
# weighted according to its barycentric coordinates, and get its
# color.
C = colors[tri.simplices[tetrahedra, np.argmax(R <= bsum, axis=1)]]
# Mask out the targets where we failed to find a tetrahedron.
C[tetrahedra == -1] = DEFAULT_COLOR
# Determine width and height of image.
# (Since I don't have your TIFF, this is the best I can do!)
width, height = 1, ntargets
# Write output as image in PPM format.
with open(OUTPUT_FILENAME, 'wb') as ppm:
 ppm.write("P3\n{} {}\n255\n".format(width, height).encode('ascii'))
 np.savetxt(ppm, C, fmt='%d')

1. How to structure your code? I would expect you to carry out the following sequence of operations. For maximum speed, each operation should be vectorized, that is, it should compute a NumPy array of results using a handful of calls to NumPy functions, with no Python for loops. (The equivalent outputs in my §4 above, if any, are in parentheses.)

   1. Load the array of vertices to triangulate (points).
   2. Perform the Delaunay triangulation, getting an array of arrays of indices of vertices (tri.simplices).
   3. Find the vertices for each tetrahedron, getting an array of arrays of vertices (points[tri.simplices]).
   4. Compute the affine transform for each tetrahedron, getting an array of inverse matrices (tri.transform[:,:3]) and an array of offset vectors (tri.transform[:,3]).
   5. Load the array of target points (targets).
   6. Compute the barycentric coordinates of each target point in each tetrahedron. (No exact equivalent in my code: I used the output of find_simplex to compute the barycentric coordinates only for the tetrahedron it was found in, but the algorithm should be similar, with suitable adjustment of the indices in the np.einsum call.)
   7. Identify the coordinates that are non-negative and so find which tetrahedron each target was found in (tetrahedra).

1. How to structure your code? I would expect you to carry out the following sequence of operations. For maximum speed, each operation should be vectorized, that is, it should compute a NumPy array of results using a handful of calls to NumPy functions, with no Python for loops. (The equivalent outputs in my §4 above, if any, are in parentheses.)

   1. Load the array of vertices to triangulate (points).
   2. Perform the Delaunay triangulation, getting an array of arrays of indices of vertices (tri.simplices).
   3. Find the vertices for each tetrahedron, getting an array of arrays of vertices (points[tri.simplices]).
   4. Compute the affine transform for each tetrahedron, getting an array of inverse matrices (tri.transform[:,:3]) and an array of offset vectors (tri.transform[:,3]).
   5. Load the array of target points (targets).
   6. Compute the barycentric coordinates of each target point in each tetrahedron. (No exact equivalent in my code: I used the output of find_simplex to compute the barycentric coordinates only for the tetrahedron it was found in, but the algorithm should be similar, with suitable adjustment of the indices in the np.einsum call.)
   7. Identify the coordinates that are non-negative and so find which tetrahedron each target was found in (tetrahedra).

2. How to read the names from the colorlist CSV? In NumPy each array must have a single data type (that's one the requirements for NumPy to be able to process them quickly). So you must read the names into one array and the points into another array. An easy way to do this is to call numpy.loadtxt twice:

    >>> names = np.loadtxt('XYZcolorlist_D65.csv', usecols=(0,), delimiter=',', dtype=str) # use dtype=bytes in Python 3
    >>> points = np.loadtxt('XYZcolorlist_D65.csv', usecols=(1,2,3), delimiter=',')
    >>> next(zip(names, points))
    ('255 63 127', array([ 35.53443021, 21.38072197, 20.3661096 ]))
   

 >>> X = tri.transform[tetrahedra,:3]
 >>> Y = targets - tri.transform[tetrahedra,3]
 >>> b = np.einsum('ijk,ik->ij', X, Y)
 >>> np.c_[b, 1 - b.sum(axis=1)]
 array([[ 0.1, 0. , 0.1, 0.8],
 [ 0.1, 0. , 0.8, 0.1],
 [ 0.6, 0.1, 0.1, 0.2],
 [ 0.1, 0.3, 0.5, 0.1]])

Note that you'll have to do something about the points that were not found in any tetrahedron. scipy.spatial.Delaunay.find_simplex returns -1 for these points. You could mask out the points you need, as described here, or you could try using a masked array. (Or try both and see which is faster.)

 >>> X = tri.transform[tetrahedra,:3]
 >>> Y = targets - tri.transform[tetrahedra,3]
 >>> b = np.einsum('ijk,ik->ij', X, Y)
 >>> bcoords = np.c_[b, 1 - b.sum(axis=1)]
 >>> bcoords
  array([[ 0.1, 0. , 0.1, 0.8],
 [ 0.1, 0. , 0.8, 0.1],
 [ 0.6, 0.1, 0.1, 0.2],
 [ 0.1, 0.3, 0.5, 0.1]])

Note that you'll have to do something about the points that were not found in any tetrahedron. scipy.spatial.Delaunay.find_simplex returns -1 for these points. You could mask out the points you need, as described here, or you could try using a masked array. (Or try both and see which is faster.)

Answers to questions

In comments you asked:

1. What is vectorization? See the "What is NumPy? " section of the NumPy documentation, in particular the section starting:

   Vectorization describes the absence of any explicit looping, indexing, etc., in the code - these things are taking place, of course, just "behind the scenes" (in optimized, pre-compiled C code).

Taking advantage of NumPy vector operations is the whole point for using NumPy. For example, this program computes 9 million multiplications and additions in Python:

 >>> v = np.arange(3000)
 >>> sum(i * j for i in v for j in v)
 20236502250000

It takes about ten seconds on this computer. But the corresponding operation in NumPy:

 >>> np.sum(np.outer(v, v))
 20236502250000

is about 100 times faster, because the loops involved in this computation run in fast machine code and not in slow Python code. So when you are working on a program in NumPy, you need to scrutinize every Python loop — that is, every for statement — to see if you can turn it into a single NumPy operation.

1. How to structure your code? I would expect you to carry out the following sequence of operations. For maximum speed, each operation should be vectorized, that is, it should compute a NumPy array of results using a handful of calls to NumPy functions, with no Python for loops. (The equivalent outputs in my §4 above, if any, are in parentheses.)

   1. Load the array of vertices to triangulate (points).
   2. Perform the Delaunay triangulation, getting an array of arrays of indices of vertices (tri.simplices).
   3. Find the vertices for each tetrahedron, getting an array of arrays of vertices (points[tri.simplices]).
   4. Compute the affine transform for each tetrahedron, getting an array of inverse matrices (tri.transform[:,:3]) and an array of offset vectors (tri.transform[:,3]).
   5. Load the array of target points (targets).
   6. Compute the barycentric coordinates of each target point in each tetrahedron. (No exact equivalent in my code: I used the output of find_simplex to compute the barycentric coordinates only for the tetrahedron it was found in, but the algorithm should be similar, with suitable adjustment of the indices in the np.einsum call.)
   7. Identify the coordinates that are non-negative and so find which tetrahedron each target was found in (tetrahedra).