Increasing efficiency of barycentric coordinate calculation in python

Question

Background: I'm attempting to warp one face to another of a different shape.

In order to warp one image to another, I'm using a delaunay triangulation of facial landmarks and warping the triangles of one portrait to the corresponding triangles of the second portrait. I'm using a barycentric coordinate system to map a point within a triangle to its corresponding warped location on the other triangle.

My first approach was to solve the system Ax = b with the inverse multiplication method, where A consists of the three corners of the triangle, b represents the current point, and x represents the barycentric coordinates of this point (alpha, beta, and gamma). I found the inverse of matrix A once per triangle, and then for every point within that triangle calculated the barycentric coordinates by finding the dot product of A^-1 and the point b. I found this to be very slow (the function takes 36 seconds to complete).

Following the recommendation of other posts, I attempted to use a least squares solution to improve the efficiency of this process. However, the time increased to 154 seconds when I used numpy's lsq method. I believe this is due to the fact that the A matrix is factored every single time the inside loop runs, while before I was able to find the inverse only one time, before the two loops begin.

My question is, how can I improve the efficiency of this function? Is there a way to store the factorization of A so that each time the least squares solution is calculated for a new point, it isn't repeating the same work?

Pseudocode for this function:

# Iterate through each triangle (and get corresponding warp triangle)
for triangle in triangulation:

    # Extract corners of the unwarped triangle
    a = firstCornerUW
    b = secondCornerUW
    c = thirdCornerUW

    # Extract corners of the warp triangle
    a_prime = firstCornerW
    b_prime = secondCornerW
    c_prime = thirdCornerW

    # This matrix will be the same for all points within the triangle
    triMatrix = matrix of a, b, and c

    # Bounding box of the triangle
    xleft = min(ax, bx, cx)
    xright = max(ax, bx, cx)
    ytop = min(ay, by, cy)
    ybottom = max(ay, by, cy)

    for x in range(xleft, xright):

        for y in range(ytop, ybottom):

            # Store the current point as a matrix
            p = np.array([[x], [y], [1]])

            # Solve for least squares solution to get barycentric coordinates
            barycoor = np.linalg.lstsq(triMatrix, p)

            # Pull individual coordinates from the array
            alpha = barycoor[0]
            beta = barycoor[1]
            gamma = barycoor[2]

            # If any of these conditions are not met, the point is not inside the triangle
            if alpha, beta, gamma > 0 and alpha + beta + gamma <= 1:

                # Now calculate the warped point by multiplying by alpha, beta, and gamma
                # Warp the point from image to warped image

Amhara

The essence to increasing efficiency here is to eliminate the loops, and exploit vectorization possibilities. On top of that, I get the impression you could make scipy.spatial do a lot of the heavy lifting for you. Could you add a more high-level expected-input-ouput type description of your problem? – Mendez

A simple optimization would be to first vectorize the loops over x and y. Just create a list of points using np.mgrid, transform all these points into barycentric space, and filter out all points that do not have exclusively positive coordinates. That should easily yield an order of magnitude in performance. But still, I don't think this solution is optimal if we take a higher-level perspective of the problem. – Mendez

btw; lstsq does not compute barycentric coordinates per se. Consider a triangle with its COM at the origin. The 'barycentric coordinates' for [0,0] as computed with lstsq are [0,0,0]; which do not sum to one. Finding the transform to barycentric cords should fundamentally include the sum-to-one constraint. – Mendez

barytransform = np.linalg.inv([[ax,bx,cx], [ay,by,cy], [1,1,1]]) That should give you a matrix transforming to barycentric space; then to get the barycentric coords: barycoord = np.dot(barytransform, [x,y,1]) – Mendez

Answer 1

Here are my suggestions, expressed in your pseudocode. Note that vectorizing the loop over the triangles should not be much harder either.

# Iterate through each triangle (and get corresponding warp triangle)
for triangle in triangulation:

    # Extract corners of the unwarped triangle
    a = firstCornerUW
    b = secondCornerUW
    c = thirdCornerUW

    # Bounding box of the triangle
    xleft = min(ax, bx, cx)
    xright = max(ax, bx, cx)
    ytop = min(ay, by, cy)
    ybottom = max(ay, by, cy)

    barytransform = np.linalg.inv([[ax,bx,cx], [ay,by,cy], [1,1,1]])     

    grid = np.mgrid[xleft:xright, ytop:ybottom].reshape(2,-1)
    grid = np.vstack((grid, np.ones((1, grid.shape[1]))))

    barycoords = np.dot(barytransform, grid)
    barycoords = barycoords[:,np.all(barycoords>=0, axis=0)]

Mendez

After computing all of the barycentric coordinates for the triangle, what would be the most efficient method of performing the actual transformation without any use of loops? – Amhara

I do not understand the higher-level problem you are trying to address well enough to answer that question with any specificity. Personally, I would guess your problem is best tackled by existing higher-level functionality, like scipy.spatial.Delaunay.find_simplex – Mendez

Answer 2

Just adding a concrete example, motivated by this post, we can iterate over the triangles (simplices) by obtaining them using scipy.spatial's Delaunay() function and computing the barycentric coordinates of the points, as shown in the next code snippet / output. Since it requires to iterate over the simplices, it should be fast when the number of simplices is small.

from scipy.spatial import Delaunay
import numpy as np
from time import time

t = time()

a, b, c = [250, 100], [100, 400], [400, 400]

tri = Delaunay(np.array([a, b, c]))

# bounding box of the triangle
xleft, xright = min(a[0], b[0], c[0]), max(a[0], b[0], c[0])
ytop, ybottom = min(a[1], b[1], c[1]), max(a[1], b[1], c[1])

xv, yv = np.meshgrid(range(xleft, xright), range(ytop, ybottom))
xv, yv = xv.flatten(), yv.flatten()

pp = np.vstack((xv, yv)).T
ss = tri.find_simplex(pp)
ndim = tri.transform.shape[2]
print(len(np.unique(ss)))
# 2

out = np.zeros((450,450,3), dtype=np.uint8)
for i in np.unique(ss): # for all simplices (triangles)
    p = pp[ss == i] # all points in the simplex
    # compute the barycentric coordinates of the points
    b = tri.transform[i, :ndim].dot(np.transpose(p) - tri.transform[i, ndim].reshape(-1,1))
    αβγ = np.c_[np.transpose(b), 1 - b.sum(axis=0)] 
    indices = np.where(np.all(αβγ>=0, axis=1))
    out[p[indices,0], p[indices,1]] = αβγ[indices]@np.array([[0,0,255], [0,255,0], [255,0,0]])

print(f'time: {time() - t} sec')
# time: 0.03899240493774414 sec

plt.imshow(out)

Meares