Create a gist now

Instantly share code, notes, and snippets.

What would you like to do?
# usage: `python ~/Desktop/hough_by_hand.py 1994-654-12_v02.tif`
# output is a plot that pops up
import cv2
import numpy as np
from matplotlib import pyplot as plt
import sys
import math
import pdb
# Loading image
if len(sys.argv) == 2:
filename = sys.argv[1]
else:
print "No input image given! \n"
img_orig = cv2.imread(filename,)
img = img_orig[:,:,::-1] # color channel plotting mess http://stackoverflow.com/a/15074748/2256243
# http://nabinsharma.wordpress.com/2012/12/26/linear-hough-transform-using-python/
def hough_transform(img_bin, theta_res=1, rho_res=1):
nR,nC = img_bin.shape
theta = np.linspace(-90.0, 0.0, np.ceil(90.0/theta_res) + 1.0)
theta = np.concatenate((theta, -theta[len(theta)-2::-1]))
D = np.sqrt((nR - 1)**2 + (nC - 1)**2)
q = np.ceil(D/rho_res)
nrho = 2*q + 1
rho = np.linspace(-q*rho_res, q*rho_res, nrho)
H = np.zeros((len(rho), len(theta)))
for rowIdx in range(nR):
for colIdx in range(nC):
if img_bin[rowIdx, colIdx]:
for thIdx in range(len(theta)):
rhoVal = colIdx*np.cos(theta[thIdx]*np.pi/180.0) + \
rowIdx*np.sin(theta[thIdx]*np.pi/180)
rhoIdx = np.nonzero(np.abs(rho-rhoVal) == np.min(np.abs(rho-rhoVal)))[0]
H[rhoIdx[0], thIdx] += 1
return rho, theta, H
def top_n_rho_theta_pairs(ht_acc_matrix, n, rhos, thetas):
'''
@param hough transform accumulator matrix H (rho by theta)
@param n pairs of rho and thetas desired
@param ordered array of rhos represented by rows in H
@param ordered array of thetas represented by columns in H
@return top n rho theta pairs in H by accumulator value
@return x,y indexes in H of top n rho theta pairs
'''
flat = list(set(np.hstack(ht_acc_matrix)))
flat_sorted = sorted(flat, key = lambda n: -n)
coords_sorted = [(np.argwhere(ht_acc_matrix == acc_value)) for acc_value in flat_sorted[0:n]]
rho_theta = []
x_y = []
for coords_for_val_idx in range(0, len(coords_sorted), 1):
coords_for_val = coords_sorted[coords_for_val_idx]
for i in range(0, len(coords_for_val), 1):
n,m = coords_for_val[i] # n by m matrix
rho = rhos[n]
theta = thetas[m]
rho_theta.append([rho, theta])
x_y.append([m, n]) # just to unnest and reorder coords_sorted
return [rho_theta[0:n], x_y]
def valid_point(pt, ymax, xmax):
'''
@return True/False if pt is with bounds for an xmax by ymax image
'''
x, y = pt
if x <= xmax and x >= 0 and y <= ymax and y >= 0:
return True
else:
return False
def round_tup(tup):
'''
@return closest integer for each number in a point for referencing
a particular pixel in an image
'''
x,y = [int(round(num)) for num in tup]
return (x,y)
def draw_rho_theta_pairs(target_im, pairs):
'''
@param opencv image
@param array of rho and theta pairs
Has the side-effect of drawing a line corresponding to a rho theta
pair on the image provided
'''
im_y_max, im_x_max, channels = np.shape(target_im)
for i in range(0, len(pairs), 1):
point = pairs[i]
rho = point[0]
theta = point[1] * np.pi / 180 # degrees to radians
# y = mx + b form
m = -np.cos(theta) / np.sin(theta)
b = rho / np.sin(theta)
# possible intersections on image edges
left = (0, b)
right = (im_x_max, im_x_max * m + b)
top = (-b / m, 0)
bottom = ((im_y_max - b) / m, im_y_max)
pts = [pt for pt in [left, right, top, bottom] if valid_point(pt, im_y_max, im_x_max)]
if len(pts) == 2:
cv2.line(target_im, round_tup(pts[0]), round_tup(pts[1]), (0,0,255), 1)
bw = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
edges = cv2.Canny(bw, threshold1 = 0, threshold2 = 50, apertureSize = 3)
rhos, thetas, H = hough_transform(edges)
rho_theta_pairs, x_y_pairs = top_n_rho_theta_pairs(H, 22, rhos, thetas)
im_w_lines = img.copy()
draw_rho_theta_pairs(im_w_lines, rho_theta_pairs)
# also going to draw circles in the accumulator matrix
for i in range(0, len(x_y_pairs), 1):
x, y = x_y_pairs[i]
cv2.circle(img = H, center = (x, y), radius = 12, color=(0,0,0), thickness = 1)
plt.subplot(141),plt.imshow(img)
plt.title('Original Image'), plt.xticks([]), plt.yticks([])
plt.subplot(142),plt.imshow(edges,cmap = 'gray')
plt.title('Image Edges'), plt.xticks([]), plt.yticks([])
plt.subplot(143),plt.imshow(H)
plt.title('Hough Transform Accumulator'), plt.xticks([]), plt.yticks([])
plt.subplot(144),plt.imshow(im_w_lines)
plt.title('Detected Lines'), plt.xticks([]), plt.yticks([])
plt.show()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment