Hough Circular Transform

Question

I'm writing a very basic crater-detection code using Hough Circular Transforms. By far, the longest part of my code is actually doing the transform, which involves numerous nested loops/conditions. I've tried several things to try to speed it up, succeeding with some but not others, but I'm also not great with writing Python-ese (all those tricks with addressing vectors that speed things up).

In the code block below, I have included a lot of commented out code so you can see what I've tried already, plus my comments that include further notes about speed.

Variables not self-contained in this block:

• M_ImageEdges is a 2D array that is the Sobel edge-detected original (slightly blurred) image.
• lon_size is a variable containing the axis-0 dimension of that image.
• lat_size is a variable containing the axis-1 dimension of that image.
• i_counter_crater is what test location this particular iteration is looking at.
• DIAM_GUESS_SEED is a 1D array of possible crater diameters to look for (would be better just to pass the one for this loop, but that's a trivial speed improvement).
• reduction_factor is an overall scaling of the image to reduce it in size for faster computation, but I need to know where things are (and how big) on the original so this variable is propagated through.

    #This program is really slow due to all the nested for() loops.  Try to decrease
    # the parameter space.
    scale_r     =  1 #legacy, but left in case you want it; basically skip every scale_r-1 values in the Accumulation Matrix
    scale_theta = 15 #skip every scale_theta-1 values when adding to the Accumulation Matrix

    #Create the Hough accumulation matrix.
    # Note: Radius can't be >50% of the window!
    Hough = np.zeros((lon_size,lat_size,DIAM_GUESS_SEED[i_counter_crater])) #initialize to zero (remember: we're looking for a maximum)

    #Create the angles we search around each point, and for speed since we
    # always calculate the same angle values, cache the sine and cosine results.
    theta = [x*scale_theta for x in range(0,int(360/scale_theta))]
    sin_cache = [math.sin(theta[x] * math.pi/180.) for x in range(0,len(theta))]
    cos_cache = [math.cos(theta[x] * math.pi/180.) for x in range(0,len(theta))]


    #Try to constrain the radius parameter space.  This will be different depending
    # on how well we know or don't know the crater size.  Since this particular
    # version is ONLY being used for a seeded-crater-detection, we are going to
    # go between ±*sqrt(2).
    radius_guess_min = round(DIAM_GUESS_SEED[i_counter_crater]/2/reduction_factor /2)
    radius_guess_max = round(DIAM_GUESS_SEED[i_counter_crater]/2/reduction_factor *2)
    if (f_outputtiming == 1):
        print("Setting up variables and waves to proceed took %.3f seconds." % (time.time()-timer_dummy))
        timer_dummy = time.time()


    #The primary Hough Transformation loop.

    #For this version, allow the crater center to be anywhere in the cell other
    # than the min/max pixels along the edge.
    for counter_cell_x in range(1,lon_size,1):
        for counter_cell_y in range(1,lat_size,1):

            #The whole point of the edge detection is to only include edge pixels in
            # this calculation, to use them in "voting" on the best circle.
            if (M_ImageEdges[counter_cell_x,counter_cell_y] >= threshold):

                #Want to test for possible circles with a radius r.
                for counter_test_r in range(radius_guess_min,radius_guess_max,scale_r):

                    #Here's the magic: If we're on a >threshold edge point, calculate where the (h,k)
                    # center of the sample cirlce with this test radius would be.  Then, increment the
                    # counter in the accumulation matrix at that circle center location.  If there are
                    # a lot of circle centers that converge on a single spot in the accumulation matrix
                    # then, voilà, you have a best-candidate circle.
                    #Python speed: Using int() is faster than round()
                    #Python speed: Using the if() statement was faster
                    # than doing min(max()) in hh and kk definitions.
                    #Python speed: math.sin/cos is much faster than np.
                    #Python speed: It is a bit faster to cache the sin/cos.
#                    for counter_cell_theta in range(0,360,scale_theta):
#                        hh = int(counter_cell_x - counter_test_r * math.cos(counter_cell_theta * math.pi/180.))
#                        kk = int(counter_cell_y - counter_test_r * math.sin(counter_cell_theta * math.pi/180.))
#                        if ( (1 <= hh) and (hh < lon_size-1) and (1 <= kk) and (kk < lat_size-1) ):
#                            Hough[hh,kk,counter_test_r] += M_ImageEdges[counter_cell_x,counter_cell_y]
                    for counter_cell_theta in range(0,len(theta)):
                        hh = int(counter_cell_x - counter_test_r * cos_cache[counter_cell_theta])
                        kk = int(counter_cell_y - counter_test_r * sin_cache[counter_cell_theta])
                        if ( (1 <= hh) and (hh < lon_size-1) and (1 <= kk) and (kk < lat_size-1) ):
                            Hough[hh,kk,counter_test_r] += M_ImageEdges[counter_cell_x,counter_cell_y]
#                    hh_cache = [int(counter_cell_x - counter_test_r * cos_cache[counter_cell_theta]) for counter_cell_theta in range(0,len(theta))]
#                    kk_cache = [int(counter_cell_y - counter_test_r * sin_cache[counter_cell_theta]) for counter_cell_theta in range(0,len(theta))]
#                    for counter_cell_theta in range(0,len(theta)):
#                        if ( (1 <= hh_cache[counter_cell_theta]) and (hh_cache[counter_cell_theta] < lon_size-1) and (1 <= kk_cache[counter_cell_theta]) and (kk_cache[counter_cell_theta] < lat_size-1) ):
#                            Hough[hh_cache[counter_cell_theta],kk_cache[counter_cell_theta],counter_test_r] += M_ImageEdges[counter_cell_x,counter_cell_y]

                            
#    #Use NumPy's argwhere to avoid the two nested for() loops.  Note: I have
#    # this commented out because I found it to be comparable speed to the loops
#    # and it did not always yield the same results (the number of points above
#    # the threshold here varied from the above, so, not sure what's going on.
#    cheatsheet = np.argwhere(M_ImageEdges >= threshold)
#
#    #Use the cheatsheet of where the values are that we can use.
#    for locale in cheatsheet:
#        counter_cell_x = int(locale[0]) #not setting these to int() makes the math take much longer, even though they seem to be integers already
#        counter_cell_y = int(locale[1]) #not setting these to int() makes the math take much longer, even though they seem to be integers already
#
#        #Want to test for possible circles with a radius r.
#        for counter_test_r in range(radius_guess_min,radius_guess_max+1,scale_r):
#
#            #Here's the magic: If we're on a >threshold edge point, calculate where the (h,k)
#            # center of the sample cirlce with this test radius would be.  Then, increment the
#            # counter in the accumulation matrix at that circle center location.  If there are
#            # a lot of circle centers that converge on a single spot in the accumulation matrix
#            # then, voilà, you have a best-candidate circle.
#            #Python speed issue: Using int() is faster than round()
#            for counter_cell_theta in range(0,360,scale_theta):
#                hh = int(counter_cell_x - counter_test_r * math.cos(counter_cell_theta * math.pi/180.))
#                kk = int(counter_cell_y - counter_test_r * math.sin(counter_cell_theta * math.pi/180.))
#                if ( (1 <= hh) and (hh < lon_size-1) and (1 <= kk) and (kk < lat_size-1) ):
#                    Hough[hh,kk,counter_test_r] += M_ImageEdges[counter_cell_x,counter_cell_y]
    
    #Normalize the Accumulation Matrix (outside the loop for speed).
    Hough /= 255.0

In sum and substance I have tried (and in some cases implemented):

• Creating a cache of commonly calculated values to save time.
• Tried using NumPy's version of some calculations instead of Math's.
• Tried using NumPy to figure out where values I want to look at are to avoid one if statement.
• Cheated by removing some fidelity in the calculation (scale_theta).
• Looked at different ways to get integers (round vs int).
• Trying a Python-esque version of calculating two variables (hh, kk) as vectors and then referencing those, but it came out to be sometimes faster and sometimes slower.
• Moved a division outside the loops.

So, that's what I've tried, and that's where I am. I feel like the hh/kk/if statement block there should be able to be Python-ized, but I haven't figured out how.

P.S. I know there's an OpenCV implementation. I am not using it for several reasons, so please let me know if my code block can be sped up, not tell me to try OpenCV.

Answer 1

Better algorithms

I'm not quite a specialist in this area; so I assume you have selected the best possible algorithm and the best possible parameters, so microoptimizations are the only thing left to do here.

Python is not fast

I don't like to say this; usually Python is fast enough if you choose the right tools. OpenCV has HCT implementation; have you tried it? If it's not an option too, maybe you can switch to C? Ok, I'll assume that's not an option too.

Data types

int is faster than float. Can using ints be an option here? You have very much rounding, and HCT description in Wiki says voting should be A[x,y,r] += 1, not += M_ImageEdges[...]. Ok, still assuming you know what you're doing.

Precalculations

You've taken sin and cos out of the loop. The same can be done to precalculate whole index offsets:

offsets = [[int(counter_test_r * math.cos(counter_cell_theta * math.pi/180.)),
            int(counter_test_r * math.sin(counter_cell_theta * math.pi/180.)),
            counter_test_r]  for counter_test_r in range(radius_guess_min,radius_guess_max,scale_r)]

#instead of loops over theta and radius:
for offset_x, offset_y, counter_test_r in offsets:
    hh = counter_cell_x - offset_x
    kk = counter_cell_y - offset_y
    if (1 <= hh < lon_size-1) and (1 <= kk < lat_size-1 ):
                        Hough[hh,kk,counter_test_r] += M_ImageEdges[counter_cell_x,counter_cell_y]

Use numpy to the max

Python is slow because of type checks. Numpy guarantees the array has only one data type, so it can process operations without checks in loops. If offsets is numpy.ndarray, we can add/subtract it as array, making calculations even faster, like:

offsets = numpy.array([[-int(counter_test_r * math.cos(counter_cell_theta * math.pi/180.)),
            -int(counter_test_r * math.sin(counter_cell_theta * math.pi/180.)),
            counter_test_r]  for counter_test_r in range(radius_guess_min,radius_guess_max,scale_r)]
...
pts = offsets + [counter_cell_x, counter_cell_y, 0]
for hh, kk, counter_test_r in pts:
    ...

But this won't make things much better. Once again, according to Wiki, "since the parameter space of the CHT is three dimensional, it may require lots of storage and computation".