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
vsint
). - 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.