hough_by_hand.py

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