在Python中提高重心坐标计算的效率

Question

背景：我正在尝试将一个不同形状的面孔变形成另一个面孔。

为了将一幅图像变形成另一幅图像，我使用了面部标记的德劳内三角剖分，并将一个肖像的三角形变形为第二个肖像的相应三角形。我使用重心坐标系将三角形内的点映射到其在另一个三角形上对应的变形位置。

我的第一种方法是使用逆乘法解决系统Ax = b，其中A由三角形的三个角组成，b表示当前点，x表示该点的重心坐标（alpha、beta和gamma）。我每个三角形只找到一次矩阵A的逆矩阵，然后对于该三角形内的每个点，通过找到A^-1和点b的点积来计算重心坐标。我发现这非常慢（函数需要36秒才能完成）。

根据其他帖子的建议，我尝试使用最小二乘解来提高此过程的效率。然而，当我使用numpy的lsq方法时，时间增加到了154秒。我认为这是因为每次内部循环运行时都要对A矩阵进行因式分解，而之前我能够在两个循环开始之前仅找到一次逆矩阵。

我的问题是，如何提高这个函数的效率？有没有一种方法可以存储A的分解，以便每次为新点计算最小二乘解时，不会重复相同的工作？

该函数的伪代码：

# 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

- keygrip_chris

---

提高效率的关键是消除循环并利用矢量化可能性。除此之外，我认为你可以让scipy.spatial为你承担更多的重负。你能否添加一个更高级别的预期输入输出类型描述来解释你的问题？ - Eelco Hoogendoorn

一个简单的优化方法是首先对x和y的循环进行向量化。只需使用np.mgrid创建一个点列表，将所有这些点转换为重心空间，并过滤掉所有没有完全正坐标的点。这应该可以轻松地提高一个数量级的性能。但是，如果我们从更高层次的角度看问题，我认为这个解决方案仍然不是最优的。 - Eelco Hoogendoorn

顺便提一下，lstsq本身不计算重心坐标。考虑一个以原点为中心的三角形。使用lstsq计算[0,0]的“重心坐标”将得到[0,0,0]，它们并不等于1。找到转换到重心坐标的方法应该基本上包括满足总和为1的约束条件。 - Eelco Hoogendoorn

barytransform = np.linalg.inv([[ax,bx,cx], [ay,by,cy], [1,1,1]])这将给你一个转换到重心空间的矩阵；然后获取重心坐标：barycoord = np.dot(barytransform, [x,y,1]) - Eelco Hoogendoorn

Answer 1

以下是我的建议，以您的伪代码表达。请注意，将三角形循环向量化也不应该更难。

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

- Eelco Hoogendoorn

---

在计算三角形的重心坐标后，如何在不使用循环的情况下执行实际变换是最有效的方法？ - keygrip_chris

我不太明白你试图解决的更高级别问题，无法具体回答那个问题。就我个人而言，我会猜测你的问题最好使用现有的更高级别功能来解决，比如scipy.spatial.Delaunay.find_simplex。 - Eelco Hoogendoorn

Answer 2

举个具体的例子，受这篇文章的启发，我们可以通过使用scipy.spatial的Delaunay()函数获取三角形（单纯形）并计算点的重心坐标来迭代遍历它们，如下面的代码片段/输出所示。由于需要遍历单纯形，因此在单纯形数量较少时速度应该很快。

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)

- Sandipan Dey