Ellipse-detection algorithm

The problem will take more than 2 minutes, so if you don't have much time I wouldn't recommend you to deal with the problem.

I found this paper on how to program "a new efficient ellipse detection method", so I thought why not just try it. I programmed it in Python but my implementation is - in contrast to the title of the paper - really slow and needs for an 325x325 image (with 6000 black pixels), like this one here, with multiple circles/ellipses on the order of 30 minutes...

Please read through my code and tell me where I can optimize things (and I think, that there are a lot of things to optimize).

I'll briefly explain the 15 steps, which are listed in the paper (if I explained one step unclear then just take a quick look at the paper):

The steps: (how I understood them)

-
Store all edge-pixels (pixels in black) in an array

-
clear the accumulator-array (you'll see the use of it later)

-
Loop through each array-entry in the "edge-pixels-array"

-
Loop through each array-entry again, check if the distance between the two coordinates (from step 3+4) is in between the min-radius and max-radius of my ellipse (min-radius and max-radius are just function parameters)

If this is true, then proceed with steps 5-14.

-
Calculate the center, orientation and the assumed length of the major axis (you can see the equations on the paper, but I don't think that they are necessary)

-
Loop through each array-entry a third time, check if the distance between this coordinate and the assumed center of the point is between the min-radius and the max-radius. It this is true, then proceed with steps 7.-9.

-
Calculate the length of minor axis using equations (if you need to look them up on the paper)

-
Add the assumed parameters of the ellipse to the accumulator-array (center, x-Width, y-Width, orientation)

-
Wait for the inner (3.) for loop to finish

-
Average all values in the accumulator-array to find the average parameters for the investigated ellipse

-

Original points

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

It's not needed in this scope, and in the scope where it is used the first thing that happens is that it's reinitialised, so this line is completely pointless.

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:

This should be one of if i < j: or if j < i: because otherwise any pair of points which is rejected as a major axis will be tested a second time with the indices reversed, doubling the workload.

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

What's the scope of these variables?

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

What are those abs for? It looks like you're trying to work around a bug in the standard library, but in that case you should document the bug clearly.

Why take the sqrt? As far as I can see distance isn't used anywhere else. Since sqrt is a monotonically increasing function over the range of interest you could instead square minR and maxR once and then compare the squared distances.

FWIW distance is also not a useful name in my opinion. distanceIJ would explain which of the many distances it is, but better still would be 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

Firstly, your indentation seems broken here. Before pasting code into any StackExchange site you should ensure that it has no tabs or that your tabstop is 4 spaces, because the indentation will be forceably converted to spaces using that tabstop.

This would be more readable with names for EdgePixel[i][0] etc. Could be ax, ay, bx, by, xp, yp, xq, yq, or anything simple and consistent.

The correct way to avoid division by zero is to use atan2 instead of 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:

Check the spelling.

What's the second skip condition about? I don't see that in the algorithm description, and it seems to reduce the ability to discriminate between multiple candidate ellipses.

Why the inconsistent use of continue for two conditions and then a massive nested if block for the third? It would seem more consistent to write

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

and then lose a level of indentation on the main content of the loop.

# 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  minR and innerDistance < maxR):

Previous comments about abs and sqrt apply also here.

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

What about k? If you store that you greatly simplify the update to ignore_indices, because it's simply a case of adding the k of the selected parameters from the accumulator rather than repeating the calculation to determine which points are on the ellipse.

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

StackExchange Code Review Q#154065, answer score: 12