Hough transform

Faisal Qureshi
Professor
Faculty of Science
Ontario Tech University
Oshawa ON Canada
http://vclab.science.ontariotechu.ca

Copyright information

© Faisal Qureshi

License

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

Lesson plan

• Hough transform
• Fitting lines to data
• Polar representation of a line
• Counting votes
• Applications

Hough transform

P.V.C. Hough, “Machine Analysis of Bubble Chamber Pictures,” Proc. Int. Conf. High Energy Accelerators and Instrumentation, 1959

Helper routines

Sampling points on a circle

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)

Line drawing

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)

Hough transform

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)

Toy data

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:]);

The set of lines passing through a point

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');

Plotting $m$ and $c$ in Hough space

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');

Relationship between Cartesian space and 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();

Counting votes in Hough space

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

• values $m = 1$ and $c=3$ received $5$ votes;
• values $m = 2$ and $c=-5$ received $2$ votes; and
• all other values for $m$ and $c$ received $0$ votes.

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

Hough transform in Practice

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)

Toy data 2

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:]);

Transforming Cartesian points into 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]

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');

Rastering c-m lines in Hough space

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} $$

Step 1: set up the c-m grid

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

Shape of m = (1024,)
Shape of c = (1024,)
Shape of mm = (1024, 1024)
Shape of cc = (1024, 1024)

dc = c[1] - c[0]
print(f'dc = {dc}')

dm = m[1] - m[0]
print(f'dm = {dm}')

dc = 0.058651026392961825
dm = 0.058651026392961825

Step 2: rasterize lines in the c-m grid.

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');