2

我编写了一些 Python 代码来完成一些图像处理工作，但是运行起来需要大量时间。我花了最后几个小时试图优化它，但我认为我已经达到了我的能力的极限。

查看分析器的输出，下面的函数占用了我的代码总时间的很大一部分。有什么办法可以加快速度吗？

def make_ellipse(x, x0, y, y0, theta, a, b):
    c = np.cos(theta)
    s = np.sin(theta)
    a2 = a**2
    b2 = b**2
    xnew = x - x0
    ynew = y - y0
    ellipse = (xnew * c + ynew * s)**2/a2 + (xnew * s - ynew * c)**2/b2 <= 1

    return ellipse

为了给出上下文，它被称为x和y作为np.meshgrid具有相当大的网格大小的输出，所有其他参数都是简单的整数值。

尽管该功能似乎花费了很多时间，但可能还有其他方法可以加快其余代码的速度。我已经把其余的代码放在了这个要点上。

任何想法将不胜感激。我试过使用 numba 和autojiting 主要功能，但这并没有多大帮助。

4

4 回答 4

3

让我们尝试结合它的调用者优化 make_ellipse。

首先，请注意a和b在许多调用中是相同的。由于 make_ellipse 每次都将它们平方，只需让调用者执行此操作即可。

其次，请注意这在我的系统np.cos(np.arctan(theta))上1 / np.sqrt(1 + theta**2)似乎稍快一些。可以使用类似的技巧从 theta 或从 cos(theta) 计算正弦（反之亦然）。

第三，不太具体，考虑短路一些最终的椭圆公式计算。例如，(xnew * c + ynew * s)**2/a2如果大于 1，则椭圆值必须为 False。如果这种情况经常发生，您可以在这些位置“屏蔽”椭圆（昂贵）计算的后半部分。我还没有彻底计划好，但请参阅 numpy.ma 了解一些可能的线索。

于 2013-07-03T15:01:51.013 回答
2

它不会在所有情况下都加快速度，但是如果您的椭圆不占用整个图像，您应该将椭圆内的点的搜索限制在其边界矩形内。我对数学很懒，所以我用谷歌搜索它并重用@JohnZwinck 的反正切技巧的整洁余弦来提出这个函数：

def ellipse_bounding_box(x0, y0, theta, a, b):
    x_tan_t = -b * np.tan(theta) /  a
    if np.isinf(x_tan_t) :
        x_cos_t = 0
        x_sin_t = np.sign(x_tan_t)
    else :
        x_cos_t = 1 / np.sqrt(1 + x_tan_t*x_tan_t)
        x_sin_t = x_tan_t * x_cos_t
    x = x0 + a*x_cos_t*np.cos(theta) - b*x_sin_t*np.sin(theta)

    y_tan_t = b / np.tan(theta) /  a
    if np.isinf(y_tan_t):
        y_cos_t = 0
        y_sin_t = np.sign(y_tan_t)
    else:
        y_cos_t = 1 / np.sqrt(1 + y_tan_t*y_tan_t)
        y_sin_t = y_tan_t * y_cos_t
    y = y0 + b*y_sin_t*np.cos(theta) + a*y_cos_t*np.sin(theta)

    return np.sort([-x, x]), np.sort([-y, y])

您现在可以将原始函数修改为以下内容：

def make_ellipse(x, x0, y, y0, theta, a, b):
    c = np.cos(theta)
    s = np.sin(theta)
    a2 = a**2
    b2 = b**2
    x_box, y_box = ellipse_bounding_box(x0, y0, theta, a, b)
    indices = ((x >= x_box[0]) & (x <= x_box[1]) & 
               (y >= y_box[0]) & (y <= y_box[1]))
    xnew = x[indices] - x0
    ynew = y[indices] - y0
    ellipse = np.zeros_like(x, dtype=np.bool)
    ellipse[indices] = ((xnew * c + ynew * s)**2/a2 +
                        (xnew * s - ynew * c)**2/b2 <= 1)
    return ellipse
于 2013-07-03T16:08:04.650 回答
2

由于除 x 和 y 之外的所有内容都是整数，因此您可以尝试最小化数组计算的数量。我想大部分时间都花在了这句话上：

ellipse = (xnew * c + ynew * s)**2/a2 + (xnew * s - ynew * c)**2/b2 <= 1

像这样的简单重写应该减少数组操作的数量：

a = float(a)
b = float(b)
ellipse = (xnew * (c/a) + ynew * (s/a))**2 + (xnew * (s/b) - ynew * (c/b))**2 <= 1

原来的 12 个数组操作现在是 10 个（加上 4 个标量操作）。我不确定 numba 的 jit 是否会尝试这个。它可能只是先做所有的广播，然后 jit 生成的操作。在这种情况下，重新排序以便立即完成常见操作应该会有所帮助。

进一步，您可以再次将其重写为

ellipse = ((xnew + ynew * (s/c)) * (c/a))**2 + ((xnew * (s/c) - ynew) * (c/b))**2 <= 1

或者

t = numpy.tan(theta)
ellipse = ((xnew + ynew * t) * (b/a))**2 + (xnew * t - ynew)**2 <= (b/c)**2

用一个标量替换另一个数组操作，并消除其他标量操作以获得 9 个数组操作和 2 个标量操作。

与往常一样，请注意输入的范围以避免舍入错误。

不幸的是，如果两个加数中的任何一个大于比较的右侧，就没有好方法可以提前计算和保释。这将是一个明显的加速，但你需要 cython（或 c/c++）来编码。

于 2013-07-03T18:01:36.470 回答
0

您可以通过使用 Cython 大大加快速度。关于如何执行此操作有一个非常好的文档。

于 2013-07-03T14:29:44.140 回答