ilyakava · March 1, 2020 07:08
diff --git a/hough_by_hand.py b/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 https://stackoverflow.com/a/15074748/2256243

 # https://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()
	# 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 https://stackoverflow.com/a/15074748/2256243

	# https://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(-qrho_res, qrho_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 = colIdxnp.cos(theta[thIdx]np.pi/180.0) + \
	rowIdxnp.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()