Faster computation of barycentric coordinates for many points

I'm just starting to understand the Python syntax and I created a module that does what I wanted, but really slow.

Here are the stats of cProfile, top 10 ordered by internal time:

ncalls  tottime  percall  cumtime  percall filename:lineno(function)
  1291716   45.576    0.000  171.672    0.000 geometry.py:10(barycentric_coords)

  6460649   31.617    0.000   31.617    0.000 {numpy.core.multiarray.array}
  2583432   15.979    0.000   15.979    0.000 {method 'reduce' of 'numpy.ufunc'
objects}
     2031   12.032    0.006  193.333    0.095 geometry.py:26(containing_tet)
  1291716   10.944    0.000   58.323    0.000 linalg.py:244(solve)
  1291716    7.075    0.000    7.075    0.000 {numpy.linalg.lapack_lite.dgesv}
  1291716    5.750    0.000    9.865    0.000 linalg.py:99(_commonType)
  2583432    5.659    0.000    5.659    0.000 {numpy.core.multiarray._fastCopyAn
dTranspose}
  1291716    5.526    0.000    7.299    0.000 twodim_base.py:169(eye)
  1291716    5.492    0.000   12.791    0.000 numeric.py:1884(identity)

To process data and create a 24 bitmap of 300*300 pixels it would take about 1 1/2 hour on my laptop with the latest intel i5, a SSD disk and 12gb RAM!

I am not sure if it's a good idea to post all the code, but if not how can you get the entire picture?

The first thing colors_PPM.py does is to access a function called geometry.

geometry.py itself calls tetgen, tetgen.py calls tetgen.exe.

The module colors_PPM.py reads a list in targets.csv and a list in XYZcolorlist_D65.csv, it excludes some elements of XYZcolorlist_D65.csv and then reads one by one the rows of targets.csv, performs a delaunay triangulation via tetgen and returns 4 names[] and 4 bcoords[].

Then random is used to choose one name by a series of if elif else tests.

Finally the result is exported in a bitmap file epson_gamut.pbm.

Do you see any way this could run faster? I know it should seem quite a mess.

.csv files structures exampl

1. Introduction

Thanks for running the profiler. As you can see from the output, most of the runtime is being spent in your containing_tet function.

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.

It looks to me as though containing_tet tries to find a tetrahedron containing point, and it does so by computing the barycentric coordinates of point within each tetrahedron and checking that all the coordinates are non-negative.

So you could have written:

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. 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]

and then in barycentric_coords it computes the transform matrix:

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

but these computations will be the same every time (they do not depend on point). It would be best to compute them just once.

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.

1. 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)

and then if I have an array of targets to query:

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

I can find which tetrahedron each point belongs to by calling scipy.spatial.Delaunay.find_simplex:

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

And then I can find the barycentric coordinates of each point within its tetrahedron like this:

>>> 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]])

This is essentially following the recipe in the scipy.spatial.Delaunay documentation, except that I transform each point using the affine transformation for the tetrahedron it was found in. Note that in the final result, all the barycentric coordinates are in the range 0–1 as you would expect.

The bit of the computation that is tricky to figure out how to vectorize is the computation of b. numpy.dot computes the dot product of two arrays, but here I need an array of dot products. I could loop over the elements of X and Y in Python, like this:

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

but that would be using slow Python iteration rather than fast NumPy vector operations. Hence the rather hairy use of numpy.einsum.

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.)

1. Answers to questions

In comments you asked:

-
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)
2023650225000

StackExchange Code Review Q#41024, answer score: 13