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

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