Note firsts:
The problem will take more then 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...

So I would please you to read through my code and tell me, where I can optimize things (and I think, that there are a lot of things to optimize)

First I think, that I'll briefly explain the 15 steps, which are listed in the paper, because I can understand, that no one wants to read the complete paper... (If I explained one step unclear then just take a quick look at the paper please):

The Steps: (how I understood them)

1. Store all edge-pixels (pixels in black) in an array
2. clear the accumulator-array (you'll see the use of it later)
3. Loop through each array-entry in the "edge-pixels-array"
4. 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
5. 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)
6. 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.
7. Calculate the length of minor axis using equations (if you need to look them up on the paper)
8. Add the assumed parameters of the ellipse to the accumulator-array (center, x-Width, y-Width, orientation)
9. Wait for the inner (3.) for-loop to finish
10. Average all values in the accumulator-array to find the average parameters for the investigated ellipse
11. Save the average ellipse parameters
12. Remove the pixels within the detected ellipse (so you don't check the pixels twice...)
13. Clear accumulator-Array
14. Wait for the outer two for-loops to finish
15. Output the found ellipses

And here is my implementation with Python:

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

I would be glad, if someone could help me speedin' up the process because it takes really, really, really, really long.. :D