椭圆检测算法

Question

这个问题需要超过2分钟，所以如果你没有太多的时间，我就不会建议你去处理这个问题。

我找到了关于如何编写“一种新的高效椭圆检测方法”的这论文，所以我想为什么不尝试一下呢？我用Python编写了它，但与论文的标题相反，我的实现非常慢，需要一个325x325图像(6000个黑色像素)，就像这个这里，有多个圆圈/椭圆，大约30分钟.

请阅读我的代码，并告诉我哪里可以优化(我认为，有很多事情要优化)。

我将简要解释论文中列出的15个步骤(如果我解释了一个不清楚的步骤，那么就简单地看一下论文)：

步骤：(我如何理解它们)

1. 将所有边缘像素(黑色像素)存储在数组中。
2. 清除累加器数组(稍后您将看到它的使用)
3. 循环遍历“边缘像素数组”中的每个数组条目。
4. 再次循环每个数组入口，检查两个坐标(从步骤3+4)之间的距离是否在我的椭圆的最小半径和最大半径之间(最小半径和最大半径只是函数参数)，如果这是真的，那么继续步骤5-14。
5. 计算主轴的中心、方向和假定长度(你可以在纸上看到方程，但我不认为它们是必要的)
6. 第三次循环每个数组入口，检查这个坐标和假定的点中心之间的距离是否在最小半径和最大半径之间。这是正确的，然后继续执行步骤7-9。
7. 用方程式计算长轴的长度(如果你需要在纸上查找的话)
8. 将椭圆的假定参数添加到累加器数组(中心、x宽度、y宽度、方向)。
9. 等待内部(3.) for循环完成
10. 累加器数组中所有值的平均值，以查找所调查椭圆的平均参数。
11. 保存平均椭圆参数
12. 删除检测到的椭圆中的像素(这样您就不会对像素进行两次检查……)
13. 清除累加器-阵列
14. 等待外部两个for循环完成
15. 输出所发现的椭圆

如果有人能帮我加速这个过程，我会很高兴，因为这需要很长时间。

import cv2
from PIL import Image
import math

def detectEllipse(filePath, minR, maxR, minAmountOfEdges):
    image = cv2.imread(outfile) # proceed with lower res.
    w, h = len(image[0]), len(image)

    # Ellipse Detection
    output = image.copy()
    gray = cv2.cvtColor(image, cv2.COLOR_BGR2GRAY)
    EdgePixel = []

    # Step 1 - Save all Edge-Pixels in an array
    for x in range(gray.shape[0]):
        for y in range(gray.shape[1]):
            if gray[x, y] == 0:
            EdgePixel.append([x,y])


    # Step 2 - Initialize AccumulatorArray and EllipsesFound-Array
    AccumulatorArray = []
    EllipsesFound = []

    # Step 3 - Loop through all pixels
    ignore_indices = set()
    for i in range(len(EdgePixel)):
        if i in ignore_indices:
            continue
        # Step 4 - Loop through all pixels a second time
        for j in range(len(EdgePixel)):
            if j in ignore_indices:
                continue
            if i != j:
            xAverage, yAverage, aAverage, bAverage, orientationAverage = 0, 0, 0, 0, 0
            # (Step 12, clear Accumulator-Array)
            AccumulatorArray = []
            
            tmp = math.sqrt(abs(math.pow(EdgePixel[i][0]-EdgePixel[j][0], 2)) + 
                  abs(math.pow(EdgePixel[i][1]-EdgePixel[j][1], 2)))
            distance = int(tmp)
            # Step 4.1 - Check if the distance is "ok"
            if(distance / 2 > minR and distance / 2 < maxR):
                # Step 5 - Calculate 4 parameters of the ellipse
                xMiddle     = (EdgePixel[i][0] + EdgePixel[j][0]) / 2
                yMiddle     = (EdgePixel[i][1] + EdgePixel[j][1]) / 2
                a           = tmp / 2
                if(EdgePixel[j][0] != EdgePixel[i][0]): # To prevent division by 0
                    orientation = math.atan((EdgePixel[j][1]-EdgePixel[i][1])/
                                            (EdgePixel[j][0]-EdgePixel[i][0]))
                else:
                    orientation = 0
                
                # Step 6 - Lop through all pixels a third time
                for k in range(len(EdgePixel)):
                    if k in ignore_indices:
                        continue
                    if len(AccumulatorArray) > minAmoutOfEdges:
                        continue
                    if k != i and k != j:
                        # Distance from x,y to the middlepoint
                        innerDistance = math.sqrt(abs(math.pow((xMiddle - EdgePixel[k][0]),2)) + 
                                        abs(math.pow((yMiddle - EdgePixel[k][1]),2)))
                        # Distance from x,y to x2,y2
                        tmp2 = math.sqrt(abs(math.pow((EdgePixel[i][0] - EdgePixel[j][0]),2)) + 
                                         abs(math.pow((EdgePixel[i][1] - EdgePixel[j][1]),2)))
                        
                        # Distance from x,y and x0,y0 has to be smaller then the distance from x1,y1 to x0,y0
                        if(innerDistance < a and innerDistance > minR and innerDistance < maxR):
                            # Step 7 - Calculate length of minor axis
                            
                            # calculate cos^2(tau):
                            tau = math.pow(((math.pow(a,2)+math.pow(innerDistance,2)-math.pow(tmp2,2))/(2*a*innerDistance)),2)  
                            bSquared = (math.pow(a, 2)*math.pow(innerDistance, 2)*(1-tau))/(math.pow(a, 2)-math.pow(innerDistance, 2)*tau)
                            # It follows:
                            b = math.sqrt(bSquared) # this is the minor axis length
                            
                            # Step 8 - Add the parameters to the AccumulatorArray
                            Data = [xMiddle, yMiddle, a, b, orientation]
                            AccumulatorArray.append(Data)
                # Step 9 (repeat from Step 6 till here)
                # Step 10 - Check if the algorithm has detected enough Edge-Points and then average the results
                if len(AccumulatorArray) > minAmoutOfEdges: 
                    # Averageing
                    for k in range(len(AccumulatorArray)):
                        tmpData = AccumulatorArray[k]
                        xAverage            += tmpData[0]
                        yAverage            += tmpData[1]
                        aAverage            += tmpData[2]
                        bAverage            += tmpData[3]
                        orientationAverage  += tmpData[4]
                    xAverage            = int(xAverage / len(AccumulatorArray))
                    yAverage            = int(yAverage / len(AccumulatorArray))
                    aAverage            = int(aAverage / len(AccumulatorArray))
                    bAverage            = int(bAverage / len(AccumulatorArray))
                    orientationAverage  = int(orientationAverage / len(AccumulatorArray))

                    # Step 11 - Save the found Ellipse-Parameters
                    EllipsesFound.append([xAverage, yAverage, aAverage, bAverage, orientationAverage])
                    
                    # Step 12 - Remove the Pixels on the detected ellipse from the Edge-Array
                    for k in range(len(EdgePixel)):
                        if ((math.pow(EdgePixel[k][0] - xAverage, 2) / math.pow((aAverage+5), 2)) + 
                               ((math.pow(EdgePixel[k][1] - yAverage, 2) / math.pow((bAverage+5), 2)))) <= 1:
                            ignore_indices.add(k)
    return

detectEllipses("LinkToImage", 5, 15, 100)

下面是该程序的概要文件，大部分时间用于pow2 (简单地乘以2个值)和math.sqrt：

65082750次函数调用(65080245次原始调用)，在22.924秒内按:内部时间n次调用时每次调用和时间调用文件名:lineno(函数)1 17.666 17.666 22.724 22.724 ElipseDetection.py:16(detectEllipses) 34239900 2.410 0.000 0.000 0.000 ElipseDetection.py:162(pow2) 15660662 1.806 0.000 1.8060.000{内置方法math.sqrt} 14612430/14612383 0.699 0.000 0.699 0.000 {内置方法builtins.len} 488000 0.062 0.000 0.062 0.000 {“列表”对象的“附加”方法}

Answer 1

原始点

步骤2-初始化AccumulatorArray和椭圆数组AccumulatorArray = []

在这个作用域中不需要它，在使用它的范围中，发生的第一件事就是重新初始化它，所以这一行完全没有意义。

对于范围内的I(len(EdgePixel))：如果我在ignore_indices中:继续#步骤4-循环第二次遍历范围中的j的所有像素(len(EdgePixel))：如果j在ignore_indices:继续如果i != j：

这应该是if i < j:或if j < i:中的一个，因为否则任何被拒绝为主轴的点都将被第二次测试，而指数将被反转，使工作负荷加倍。

xAverage, yAverage, aAverage, bAverage, orientationAverage = 0, 0, 0, 0, 0

这些变量的范围是什么？

tmp = math.sqrt(abs(math.pow(EdgePixel[i][0]-EdgePixel[j][0], 2)) + abs(math.pow(EdgePixel[i][1]-EdgePixel[j][1], 2))) distance = int(tmp) # Step 4.1 - Check if the distance is "ok" if(distance / 2 > minR and distance / 2 < maxR):

那些abs是做什么用的？看起来，您正在试图绕过标准库中的一个bug，但在这种情况下，您应该清楚地记录该bug。

为什么要拿sqrt？据我所知，distance不在其他地方使用。由于sqrt是一个在感兴趣范围内单调增加的函数，所以您可以将minR和maxR相乘一次，然后比较平方距离。

在我看来，FWIW distance也不是一个有用的名字。distanceIJ可以解释它是多个距离中的哪一个，但最好还是majorAxis。

# Step 5 - Calculate 4 parameters of the ellipse xMiddle = (EdgePixel[i][0] + EdgePixel[j][0]) / 2 yMiddle = (EdgePixel[i][1] + EdgePixel[j][1]) / 2 a = tmp / 2 if(EdgePixel[j][0] != EdgePixel[i][0]): # To prevent division by 0 orientation = math.atan((EdgePixel[j][1]-EdgePixel[i][1])/ (EdgePixel[j][0]-EdgePixel[i][0])) else: orientation = 0

首先，你的压痕在这里好像断了。在将代码粘贴到任何StackExchange站点之前，您应该确保它没有制表符，或者您的制表符是4个空格，因为缩进将强制转换为使用该制表符的空格。

EdgePixel[i][0]等的名称可能是ax, ay, bx, by、xp, yp, xq, yq或任何简单和一致的名称，这将更容易读懂。

避免除以零的正确方法是使用atan2而不是atan。

# Step 6 - Lop through all pixels a third time for k in range(len(EdgePixel)): if k in ignore\_indices: continue if len(AccumulatorArray) > minAmoutOfEdges: continue if k != i and k != j:

检查拼写。

第二个跳过的情况是什么？我在算法描述中没有看到这一点，它似乎降低了区分多个候选椭圆的能力。

为什么在两个条件下不一致地使用continue，第三个条件则使用大量嵌套的if块？写起来似乎更一致

            for k in range(len(EdgePixel)):
                if k in ignore_indices or k == i or k == j:
                    continue

然后在循环的主要内容上丢失一个缩进级别。

# Distance from x,y to the middlepoint innerDistance = math.sqrt(abs(math.pow((xMiddle - EdgePixel[k][0]),2)) + abs(math.pow((yMiddle - EdgePixel[k][1]),2))) # Distance from x,y to x2,y2 tmp2 = math.sqrt(abs(math.pow((EdgePixel[i][0] - EdgePixel[j][0]),2)) + abs(math.pow((EdgePixel[i][1] - EdgePixel[j][1]),2))) # Distance from x,y and x0,y0 has to be smaller then the distance from x1,y1 to x0,y0 if(innerDistance < a and innerDistance > minR and innerDistance < maxR):

以前关于abs和sqrt的评论也适用于这里。

# Step 8 - Add the parameters to the AccumulatorArray Data = [xMiddle, yMiddle, a, b, orientation] AccumulatorArray.append(Data)

那k呢？如果存储的话，可以大大简化对ignore_indices的更新，因为这只是从累加器中添加所选参数的k，而不是重复计算以确定椭圆上有哪些点。

# Step 10 - Check if the algorithm has detected enough Edge-Points and then average the results if len(AccumulatorArray) > minAmoutOfEdges: # Averageing for k in range(len(AccumulatorArray)): tmpData = AccumulatorArray[k] xAverage += tmpData[0] yAverage += tmpData[1] aAverage += tmpData[2] bAverage += tmpData[3] orientationAverage += tmpData[4] xAverage = int(xAverage / len(AccumulatorArray)) yAverage = int(yAverage / len(AccumulatorArray)) aAverage = int(aAverage / len(AccumulatorArray)) bAverage = int(bAverage / len(AccumulatorArray)) orientationAverage = int(orientationAverage / len(AccumulatorArray))

这是不对的。该规范要求您找到最常见的短轴(可能存在一定的误差，尽管这并不十分清楚)，并将与该小轴对应的点作为椭圆上的点。

细节有点棘手。在生成参数的输出时，您将转换为int，那么您是否假定所有的小轴都应该是整数？如果是这样的话，这使得分组相对容易。否则，最简单的实现将是添加一个参数bMargin，对b值进行排序，并扫描以查找b的值，该值在范围[b, b + bMargin]中的点数最多。然后，如果该集群有minAmountOfEdges点，则取集群中的点，平均它们的参数以获得EllipsesFound的值，并将它们的k值添加到ignore_indices中。

附加点

tmp2 = math.sqrt(abs(math.pow((EdgePixel[i][0] - EdgePixel[j][0]),2)) + abs(math.pow((EdgePixel[i][1] - EdgePixel[j][1]),2)))

这是另一个窃听器。它应该使用k而不是i。

为什么要累积xMiddle、yMiddle、a和orientation的值，然后取它们的平均值？Data中唯一真正依赖于k的元素是b，因此唯一值得积累和平均的元素。

通过查看如何使用sqrt的结果，并观察到它们几乎都是平方的，就可以简化很多。这是步骤3至9(无可否认，评论太少，需要一个比foo更好的名称)：

ignore_indices = set()
for i in range(len(EdgePixel)):
    if i in ignore_indices:
        continue

    px, py = EdgePixel[i][0], EdgePixel[i][1]

    for j in range(len(EdgePixel)):
        if j <= i or j in ignore_indices:
            continue

        qx, qy = EdgePixel[j][0], EdgePixel[j][1]
        majorAxisSq = math.pow(px - qx, 2) + math.pow(py - qy, 2)

        if (majorAxisSq < minAxisSq or majorAxisSq > maxAxisSq):
            continue

        AccumulatorArray = []
        cx, cy = (px + qx) / 2, (py + qy) / 2

        for k in range(len(EdgePixel)):
            if k == i or k == j or k in ignore_indices:
                continue

            rx, ry = EdgePixel[k][0], EdgePixel[k][1]
            crSq = math.pow(cx - rx, 2) + math.pow(cy - ry, 2)

            # Require |C-R| < majorAxis and point not so close to centre that minor axis (< |C-R|) is too small
            if crSq > majorAxisSq or crSq < minAxisSq:
                continue

            # Calculate length of minor axis. TODO Numerical analysis
            var fSq = math.pow(rx - qx, 2) + math.pow(ry - qy, 2);
            foo = math.pow(majorAxisSq/4 + crSq - fSq, 2)
            bSquared = (majorAxisSq * crSq * (majorAxisSq * crSq - foo)) / (majorAxisSq - 4 * crSq * foo)
            minorSemiAxis = math.sqrt(bSquared)

            AccumulatorArray.append([minorSemiAxis, k])

高级建议

从过滤器的角度来看，最大的速度改进将是什么？有一个最大允许的长轴；一旦你要寻找可能的边缘点，就会有一个约束，那就是每个点都必须在圆周内，它的直径是EdgePoints[i]和EdgePoints[j]。与其扫描每个边缘像素以查看它是否符合这些条件，还不如找到或实现一个几何查找数据结构。它不一定要支持圆形或半圆形的查找:包围框应该足够好，以减少测试的点数。作为第一个实现，除非有合适的库支持，否则我会尝试一个简单的四叉树。

Answer 2

太多索引

Python不是一种编译语言。每次建立索引时，Python解释器都会重新执行索引查找。它不能确定一个值是否不变，也不需要重复索引。

只考虑以下几点：

    for i in range(len(EdgePixel)):
        ...
        for j in range(len(EdgePixel)):
            ...
            if ...:
                for k in range(len(EdgePixel)):
                    if k != i and k != j:
                        innerDistance = math.sqrt(abs(math.pow((xMiddle - EdgePixel[k][0]),2)) + 
                                        abs(math.pow((yMiddle - EdgePixel[k][1]),2)))
                        tmp2 = math.sqrt(abs(math.pow((EdgePixel[i][0] - EdgePixel[j][0]),2)) + 
                                         abs(math.pow((EdgePixel[i][1] - EdgePixel[j][1]),2)))

由于您有一个三重嵌套循环，内环执行大约N^3次数。在innerDistance计算中，EdgePixel[k]被计算了两次。第一次获取EdgePixel[k][0]值，然后第二次获取EdgePixel[k][1]值。这个索引不是一个简单的操作。每次EdgePixel首先在本地字典中查找，然后在本地字典中查找k，然后用这两个值调用PyObject_GetItem(PyObject *o, PyObject *key)来获取EdgePixel[k]的值！

由于循环中的值没有变化，所以在for语句本身中循环时查找它们一次：

    for i, (xi, yi) in enumerate(EdgePixel):
        ...
        for j, (xj, yj) in enumerate(EdgePixel):
            ...
            if ...:
                for k, (xk, yk) in enumerate(EdgePixel):
                    if k != i and k != j:
                        innerDistance = math.sqrt(abs(math.pow((xMiddle - xk),2)) + 
                                        abs(math.pow((yMiddle - yk),2)))
                        tmp2 = math.sqrt(abs(math.pow((xi - xj),2)) + 
                                         abs(math.pow((yi - yj),2)))

由于[0]和[1]已经被语义上更有意义的x和y命名的变量所取代，这不仅速度更快，而且更容易阅读。

将不变量移出循环

tmp2不依赖于k或最内部循环中的任何其他值。它可以在k循环之外进行计算。

    for i, (xi, yi) in enumerate(EdgePixel):
        ...
        for j, (xj, yj) in enumerate(EdgePixel):
            ...
            if ...:
                tmp2 = math.sqrt(abs(math.pow((xi - xj),2)) + 
                                 abs(math.pow((yi - yj),2)))
                for k, (xk, yk) in enumerate(EdgePixel):
                    if k != i and k != j:
                        innerDistance = math.sqrt(abs(math.pow((xMiddle - xk),2)) + 
                                        abs(math.pow((yMiddle - yk),2)))

当然，除非计算是不正确的，并且应该依赖于k，如Peter所指出的。在这种情况下，这个运动只是在优化错误的代码。;-)