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 < a and 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]
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))
This is just wrong. The specification calls for you to find the most common minor axis (presumably with some margin of error, although that's not stunningly clear) and to take just the points which correspond to that minor axis as the points on the ellipse.
The details are a bit tricky. You're casting to int when producing the output of the parameters, so do you assume that all minor axes should be integers? If so, that makes the grouping relatively easy. Otherwise the easiest implementation would be to add a parameter bMargin, sort the values of b and scan to find the value of b which has most points in the range [b, b + bMargin].
Additional points
tmp2 = math.sqrt(abs(math.pow((EdgePixel[i][0] - EdgePixel[j][0]),2)) +
abs(math.pow((EdgePixel[i][1] - EdgePixel[j][1]),2)))
This is another bug. It should be using k rather than i.
Why accumulate values of a and orientation and then take their averages? You already know that aAverage will be a and orientationAverage will be orientation.
By looking at how the results of sqrt were used and observing that they're nearly all squared again, it's possible to simplify a lot. This is steps 3 to 9 (admittedly with too few comments and needing a better name than 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([cx, cy, minorSemiAxis, k])
Higher level suggestions
Looking at the filters it's become clear what the biggest speed improvement would be. There's a maximum permitted major axis; once you're looking for possible edge points, there's a constraint that every point must be inside the circle which has EdgePoints[i] and EdgePoints[j] as its diameter. Rather than scanning every single edge pixel to see whether it meets those criteria, you should find or implement a geometric lookup data structure. It doesn't necessarily have to support circular or semicircular lookups: bounding boxes should be good enough to reduce the points tested considerably. As a first implementation, unless there's suitable library support, I'd try a simple quadtree.