Faisal Qureshi
Professor
Faculty of Science
Ontario Tech University
Oshawa ON Canada
http://vclab.science.ontariotechu.ca
© Faisal Qureshi
This work is licensed under a Creative Commons Attribution-NonCommercial 4.0 International License.
import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
import matplotlib.lines as mlines
import matplotlib.cm as cm
from sklearn.cluster import KMeans
P.V.C. Hough, “Machine Analysis of Bubble Chamber Pictures,” Proc. Int. Conf. High Energy Accelerators and Instrumentation, 1959
def sample_circle(n=1, center=np.zeros(2), radius=1):
"""
Returns an nx2 array
"""
angles = np.random.rand(n)*2*np.pi
pts = np.empty([2,n])
for i in range(n):
a = angles[i]
r = np.array([np.cos(a), np.sin(a), -np.sin(a), np.cos(a)]).reshape(2,2)
pts[:,i] = center + np.dot(r, np.array([radius,0]).reshape(2,1)).reshape(2)
return pts.T
# x = sample_circle(10000, center=np.array([.5,.5]), radius=.3)
# plt.figure(figsize=(10,10))
# plt.plot(x[:,0], x[:,1], '.')
# plt.xlim(-2,2)
# plt.ylim(-2,2)
def draw_line(pt1, pt2, ax, **kwargs):
"""
Draws a line passing through two points.
"""
xmin, xmax = ax.get_xbound()
ymin, ymax = ax.get_ybound()
scale = np.sqrt((xmax - xmin)**2 + (ymax - ymin)**2)
d = pt2 - pt1
if np.isclose(np.linalg.norm(d), 0):
print(f'draw_line warning: pt1={pt1}, pt2={pt2}')
s = pt1 - d * scale
e = pt2 + d * scale
l = mlines.Line2D([s[0], e[0]], [s[1], e[1]], **kwargs)
ax.add_line(l)
def draw_line_mc(theta, ax, **kwargs):
"""
Draws a line with slope m and y-intercept c
"""
m = theta[0]
c = theta[1]
xmin, xmax = ax.get_xbound()
ymin = m * xmin + c
ymax = m * xmax + c
l = mlines.Line2D([xmin, xmax], [ymin,ymax], **kwargs)
ax.add_line(l)
def draw_line_polar(coord, ax, **kwargs):
"""
Draws a line expressed in polar coordinates. Here coord[0] is theta and coord[1] is rho. Theta is in radians.
"""
theta = coord[0]
rho = coord[1]
pt = np.array([
rho * np.cos(theta),
rho * np.sin(theta)
])
n = pt / np.linalg.norm(pt)
c = np.dot(n, pt)
x2 = pt[0]+1
y2 = (c - n[0]*x2) / n[1]
return draw_line(pt, np.array([x2, y2]), ax, **kwargs)
# plt.figure(figsize=(10,10))
# draw_line(np.array([1,2]), np.array([1,2]), plt.gca())
# plt.xlim(-2,2)
# plt.ylim(-2,2)
def compute_m_and_c(pt1, pt2):
A = np.array([pt1[0], 1, pt2[0], 1]).reshape(2,2)
y = np.array([pt1[0], pt2[1]]).reshape(2,1)
return np.dot(np.linalg.inv(A), y)
Lets generate some toy date to play with
def toy_data1():
true_l1 = [1, 3]
x1 = np.linspace(1,5,5)
y1 = true_l1[0] * x1 + true_l1[1]
true_l2 = [2, -5]
x2 = np.linspace(5,7,3)
y2 = true_l2[0] * x2 + true_l2[1]
data = np.hstack([x1, x2, y1, y2]).reshape(2,-1).T
return data, (0, 12, 0, 12)
data, extents = toy_data1()
plt.figure(figsize=(7,7))
plt.title('Data')
plt.xlabel('x')
plt.ylabel('y')
plt.scatter(data[:,0], data[:,1], c='red')
plt.xlim(extents[:2])
plt.ylim(extents[-2:]);
Lets pick a point and draw the set of lines that pass through them
id = 3
pt = data[id,:]
np.random.seed(0)
num_lines = 10
tmp = sample_circle(n=num_lines, center=pt, radius=1)
cmap = cm.get_cmap('Spectral')
mc = np.empty((num_lines, 2))
plt.figure(figsize=(7,7))
plt.title(f'Lines through pt=({pt[0]},{pt[1]})')
plt.xlabel('x')
plt.ylabel('y')
plt.scatter(data[:,0], data[:,1], c='blue', alpha=0.1)
plt.xlim(extents[:2])
plt.ylim(extents[-2:])
for i in range(num_lines):
draw_line(pt, tmp[i,:], plt.gca(), color=cmap(i/num_lines))
mc[i,:] = compute_m_and_c(pt, tmp[i,:]).reshape(2)
plt.scatter(data[id,0], data[id,1], c='red');
We can plot $m$ and $c$ for every line that passes through a point $(x,y)$ in the Cartesian space in Hough space.
smallest = np.min(mc.flatten())
largest = np.max(mc.flatten())
plt.figure(figsize=(15,7))
plt.subplot(121)
plt.xlabel('x')
plt.ylabel('y')
plt.scatter(data[:,0], data[:,1], c='blue', alpha=0.1)
plt.xlim(extents[:2])
plt.ylim(extents[-2:])
for i in range(num_lines):
draw_line(pt, tmp[i,:], plt.gca(), color=cmap(i/num_lines))
plt.scatter(data[id,0], data[id,1], c='red')
plt.title(f'Lines through pt=({pt[0]},{pt[1]})')
plt.subplot(122)
plt.xlim(smallest, largest)
plt.ylim(smallest, largest)
draw_line(mc[0,:], mc[1,:], plt.gca())
plt.scatter(mc[:,0], mc[:,1], c=range(num_lines), cmap=cmap)
plt.xlabel('m')
plt.ylabel('c')
plt.title('Hough space');
The set of lines through a poing $(x,y)$ is $y = m x + c$. Any such line can be represented as a point $(m,c)$ in the Hough space. Say lines $\{(m_1,c_1), (m_2,c_2), \cdots\}$ all pass through point $(x,y)$, then these points form a line in the Hough space. This can be seen if we re-write $y = m x + c$ as $c = - x m + y$. This suggests that a point $(x,y)$ in Cartesian space is a line with slope $-x$ and y-intercept $y$ in the Hough space.
n = data.shape[0]
mc = np.empty((n,2))
for i in range(n):
mc[i,0] = - data[i,0]
mc[i,1] = data[i, 1]
smallest = np.min(mc.flatten())
largest = np.max(mc.flatten())
plt.figure(figsize=(15,15))
plt.subplot(221)
plt.xlabel('x')
plt.ylabel('y')
plt.title('Data points in Cartesian coordinates')
plt.scatter(data[-3:,0], data[-3:,1], c='blue', alpha=0.5)
plt.scatter(data[:5,0], data[:5,1], c='red', alpha=0.9)
plt.xlim([0,12])
plt.ylim([0,12])
plt.subplot(222)
plt.title('Corresponding lines in Hough space')
plt.xlim(-8, 8)
plt.ylim(-8, 8)
plt.grid(which='both')
for i in range(n-3):
draw_line_mc(mc[i,:], plt.gca(), c='gray')
for i in range(5,8):
draw_line_mc(mc[i,:], plt.gca(), c='gray')
plt.xlabel('m')
plt.ylabel('c');
plt.subplot(223)
plt.title('Lines for data points shown in red (Zoomed in)')
plt.grid(which='both')
plt.xlim(0, 2)
plt.ylim(2, 4)
plt.xlabel('m')
plt.ylabel('c');
for i in range(n-3):
draw_line_mc(mc[i,:], plt.gca(), c='red')
for i in range(n-3,n):
draw_line_mc(mc[i,:], plt.gca(), c='gray', alpha=0.2)
plt.scatter([1], [3], label='(1,3)', c='black', alpha=1)
plt.legend()
plt.subplot(224)
plt.title('Lines for data points shown in blue (Zoomed in)')
plt.grid(which='both')
plt.xlim(1, 3)
plt.ylim(-4, -7)
plt.xlabel('m')
plt.ylabel('c');
for i in range(n-3):
draw_line_mc(mc[i,:], plt.gca(), c='gray', alpha=0.2)
for i in range(n-3,n):
draw_line_mc(mc[i,:], plt.gca(), c='blue')
plt.scatter([2], [-5], label='(2,-5)', c='black')
plt.legend();
It is as if each 2D point in Cartesian space is voting for multiple possible values for $m$ and $c$ in Hough space. We can identify the line that passes through the Cartesian points by picking values for $m$ and $c$ that have the most votes. In the above example
fig = plt.figure(figsize=(7,7))
ax = fig.add_subplot(111)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_title('Data points in Cartesian coordinates')
ax.scatter(data[-3:,0], data[-3:,1], c='blue', alpha=0.5)
ax.scatter(data[:5,0], data[:5,1], c='red', alpha=0.9)
ax.set_xlim(extents[:2])
ax.set_ylim(extents[-2:])
# Draw first line
draw_line_mc([1,3], ax, label='m=1, c=3', c='red')
# Draw second line
draw_line_mc([2, -5], ax, label='m=2, c=-5', c='blue')
ax.legend();
In more realistic settings, data is corrupted by noise. This makes counting votes in Hough space challenging.
def toy_data2():
np.random.seed(0)
true_l1 = [0, 3]
x1 = np.linspace(-5,5,40)
y1 = true_l1[0] * x1 + true_l1[1] + np.random.randn(40)*.2
true_l2 = [0, -10]
x2 = np.linspace(-10,10,50)
y2 = true_l2[0] * x2 + true_l2[1] + np.random.randn(50)*.3
true_l3 = [2, -14]
x3 = np.linspace(-0,10,60)
y3 = true_l3[0] * x3 + true_l3[1] + np.random.randn(60)*.4
data = np.hstack([x1, x2, x3, y1, y2, y3]).reshape(2,-1).T
return data, (-15, 15, -15, 15)
data, extents = toy_data2()
plt.figure(figsize=(15,15))
plt.subplot(221)
plt.xlabel('x')
plt.ylabel('y')
plt.title('Noisy data')
plt.scatter(data[:,0], data[:,1], c='blue', alpha=0.5)
plt.xlim(extents[:2])
plt.ylim(extents[-2:]);
n = data.shape[0]
mc = np.empty((n,2))
for i in range(n):
mc[i,0] = - data[i,0]
mc[i,1] = data[i, 1]
cmap = cm.get_cmap('Spectral')
plt.figure(figsize=(15,7))
plt.subplot(121, aspect=True)
plt.xlabel('x')
plt.ylabel('y')
plt.title('Data points in Cartesian coordinates')
plt.scatter(data[:,0], data[:,1], c='blue', alpha=0.5)
plt.xlim([-15, 15])
plt.ylim([-15, 15])
plt.subplot(122, aspect=True)
plt.title('Corresponding lines in Hough space')
plt.xlim(-30, 30)
plt.ylim(-30, 30)
plt.grid(which='both')
for i in range(n):
draw_line_mc(mc[i, :], plt.gca(), c=cmap(i/n), alpha=0.3)
plt.xlabel('m')
plt.ylabel('c');
Lets rasterize these lines in Hough space on a c-m grid. The idea is as follows. Each cell in this c-m grid starts at 0. Each line that passes through a c-m grid cell increments its count by one. Recall that given a point $(x,y)$ in the Cartesian space, we can draw the corresponding line in the Hough space (i.e., c-m space) using
$$ c = -mx + y. $$or
$$ m = \frac{y - c}{x} $$# def define_grid(i, j):
# i_values = np.linspace(i[0], i[1], i[2])
# j_values = np.linspace(j[0], j[1], j[2])
# ii, jj = np.meshgrid(i_values, j_values, indexing='ij')
# ii = ii.flatten()
# jj = jj.flatten()
# di = np.abs(i_values[1] - i_values[0])
# dj = np.abs(j_values[1] - j_values[0])
# return {
# 'extents': np.vstack([i,j]),
# 'i': i_values,
# 'j': j_values,
# 'ii': ii,
# 'jj': jj
# }
m_min, m_max, m_num = -30, 30, 1024
c_min, c_max, c_num = -30, 30, 1024
m = np.linspace(m_min, m_max, m_num)
c = np.linspace(c_min, c_max, c_num)
print(f'Shape of m = {m.shape}')
print(f'Shape of c = {c.shape}')
mm, cc = np.meshgrid(m, c, indexing='ij')
print(f'Shape of mm = {mm.shape}')
print(f'Shape of cc = {cc.shape}')
mm = mm.flatten()
cc = cc.flatten()
dc = c[1] - c[0]
print(f'dc = {dc}')
dm = m[1] - m[0]
print(f'dm = {dm}')
We are using a particularly inefficient scheme for doing it. I suggest you refer to your "computer graphics" notes to find more efficient schemes for rasterizing lines: 1) Midpoint Algorithm (Pitteway, 1967), (van Aken and Novak, 1985) and 2) Bresenham Line Algorithm (Bresenham, 1965).
Our approach is simple: for each cell $i$ in the c-m grid, we evaluate if the line satisfies that cell.
#counts = np.zeros((m_num, c_num))
# mcgrid = np.vstack([mm, cc])
# print(f'Shape of mc = {mc.shape}')
# for i in range(n):
# x = data[i,0]
# y = data[i,1]
# if x > 1 or x < -1:
# estimated_m = (y - mcgrid[1,:]) / x
# diffs = np.abs(estimated_m - mcgrid[0,:])
# line = diffs < .5*dm
# else:
# estimated_c = y - x * mcgrid[0,:]
# diffs = np.abs(estimated_c - mcgrid[1,:])
# line = diffs < .5*dc
# line = line.astype('int').reshape(m_num, c_num)
# counts = counts + line
def raster_hough_line(data, counts, mm, cc, dm, dc):
m_num, c_num = counts.shape
n = data.shape[0]
for i in range(n):
x = data[i,0]
y = data[i,1]
if x > 1 or x < -1:
estimated_m = (y - cc) / x
diffs = np.abs(estimated_m - mm)
votes = diffs < .5*dm
else:
estimated_c = y - x * mm
diffs = np.abs(estimated_c - cc)
votes = diffs < .5*dc
votes = votes.astype('int').reshape(m_num, c_num)
counts = counts + votes
return counts
counts = np.zeros((m_num, c_num))
counts = raster_hough_line(data, counts, mm, cc, dm, dc)
plt.figure(figsize=(10,10))
plt.title('Hough space - counting votes')
plt.imshow(counts.T, origin='lower', extent=[m_min, m_max, c_min, c_max])
plt.xlabel('m')
plt.ylabel('c');