Implementing the Hough Transform from Scratch

How to implement the Hough Transform from scratch and some practical tips.

In this blog post, I want to teach you how to implement a powerful line detection tool: the Hough Transform.

All the code used for writing this tutorial can be found at https://github.com/AlbertoFormaggio1/Hough-Transform-From-Scratch.

For any question, don’t hesitate also to reach out at alb.formaggio@gmail.com.

The Hough Transform is a widely applied algorithm in computer vision for feature extraction. In theory, it can detect any kind of shape, e.g. lines, circles, ellipses, etc. In this tutorial, however, we will stick with the detection of lines only. People in general refrain from using this approach due to its complexity (which we will see later on).

Algorithm

There exist many tutorials over the internet explaining very well how the Hough transform works, you can check this one for example.

First and foremost, we need to say that we will work on the result produced by an edge detector since they extract the edges from the image and thus the possible candidates for being lines.

The main idea is the following: the bundle of lines passing through a point (xi,yi) in the image space is defined as yi = a * xi + b.
This bundle corresponds to exactly one line when you’re moving to the parameter space where b = -xi * a + yi. This is because (a,b) are the independent and dependent variables, respectively, and (xi, yi) are known from the image space.

We can therefore extend this to the fact that one line passing through 2 points in the image space y = a’ * x + b’, can be represented as 2 lines that intersect in one point (a’, b’) in the parameter space.
If we have n points, we will have n lines intersecting in the parameter space at the same point.

The procedure thus will be the following:

1. Transform the image in the parameter space
2. Look for the points where a high number of lines intersect: they correspond to a large number of points aligned over the same line.

Easy, right? Well, this is not the final process we will adopt. In fact, there is a big flaw: what would happen if the line to be detected is vertical? The coefficient of such a line is infinite, making the line impossible to be represented.

Normalized representation:

Instead of considering the image space as (x,y), we can transform it into polar coordinates (rho, theta). This way, there is no problem in the detection of any inclination of the lines.

Another important thing to point out is that we need to discretize the space since because of noise and other external factors it’s almost impossible to have perfect lines and we also need a computable way to find these lines: we can’t use infinite precision! We therefore divide the parameter space in accumulation cells: every line that passes through that cell will give a “vote” to that cell. Finally, the positions with many votes correspond to the position where multiple lines intersect in the parameter space and thus many points are aligned in the image space.

Pseudocode

1. Compute the edge detection to get the edge points (e.g. Canny).
2. The parameter space is quantized in cells, and a counter is associated with each cell. Theta is in the range [-pi/2, pi/2] and rho in the range [-D, D] where D is the size of the diagonal of the image.
3. For each edge point (x,y):
   a. Let theta assume all the values in the quantized range (-pi/2, pi/2) and compute rho = x * cos(theta) + y * sin(theta).
   b. For each cell crossed increment by 1 the counter in position (theta, rho)
4. The counter for each cell contains the number of pixels collinear on that line.
   Pick all the cells that have a number of counts higher than some user-defined threshold.

The computation time for this algorithm is O ( edge pixels * num of quantization cells).

Implementing the Hough Transform

I talked more than I wanted, so let’s start writing the code.

Setting up the workspace

Before starting, we need to import all the needed libraries and import the image we want to use.

import numpy as np
import cv2 as cv
from matplotlib import pyplot as plt

image_path = "path_to_your_img.jpg"
image = cv.imread(image_path)

cv.namedWindow("image")
cv.imshow("image", image)
cv.waitKey(0)
cv.destroyAllWindows()

Input image from which we want to extract the lines of the road lane the car is traveling on.

Apply the Canny Edge Detector

Among all the existing edge detectors, Canny is the one that is most widely used and that leads to the best performance, so we will use it for this exercise.
Note: even though the OpenCV implementation of Canny works also with color images, it is known that edge detection achieves better results when processing images with one single channel.

grey_image = cv.cvtColor(image, cv.COLOR_BGR2GRAY)

# The threshold values should be adapted to the image you have 
# in order to achieve meaningful edges extraction
tl = 634
th = 854

cv.Canny(grey_image, self.get_tl(), self.get_th())

win_name = "canny"
cv.namedWindow(win_name)
cv.imshow(win_name, canny_edges)
cv.waitKey(0)
cv.destroyAllWindows()

The result of the Canny Edge detector: in our case, we are left with mainly the lines in which we are interested (i.e. the lines of the lane)

Hough Implementation

We finally get to the interesting part of this post: we will create a function getting to input the extracted edges, the count threshold, and the min and max theta that will be used for the cell quantization. In the end, it will return the lines found defined by (rho, theta).
Remember, (rho, theta) are in the parameter space and they correspond to a line in the image space.

rho = 9
theta = 0.261
threshold = 101

def hough_lines(edges: np.ndarray, threshold: float, min_theta: float, max_theta: float) -> np.ndarray:
    # Initialize the counter matrix in polar coordinates
    diagonal = np.sqrt(image.shape[0]**2 + image.shape[1]**2)

    # Compute the values for the thetas and the rhos
    theta_angles = np.arange(min_theta, max_theta, theta)
    rho_values = np.arange(-diagonal, diagonal, rho)
    # Compute the dimension of the accumulator matrix
    num_thetas = len(theta_angles)
    num_rhos = len(rho_values)
    accumulator = np.zeros([num_rhos, num_thetas])
    print('Accumulator shape (rhos x thetas):' + str(accumulator.shape))

Here we just created the accumulator matrix: the angles are in a range [min_theta, max_theta] by increasing in each cell the angle value by theta, a hyperparameter. The same thing is done for rho between [-diagonal, diagonal].

    # Pre-compute sin and cos
    sins = np.sin(theta_angles)
    coss = np.cos(theta_angles)
    
    # Consider edges only
    xs, ys = np.where(edges > 0)
    
    for x,y in zip(xs,ys):
        for t in range(num_thetas):
            # compute the rhos for the given point for each theta
            current_rho = x * coss[t] + y * sins[t]
            # for each rho, compute the closest rho among the rho_values below it
            # the index corresponding to that rho is the one we will increase
            rho_pos = np.where(current_rho > rho_values)[0][-1]
            #rho_pos = np.argmin(np.abs(current_rho - rho_values))
            accumulator[rho_pos, t] += 1

We pre-compute the sines and cosines since the angles will be always the same since the angles theta are fixed. This way we can create our loop more efficiently. Then, we extract the indices of the positions in the image where there is an edge.

What we do then, is simply identifying the accumulation cells that are intersected by each of the points in the parameter space and increasing the corresponding counter.

    # Take the polar coordinates most matched
    final_rho_index, final_theta_index = np.where(accumulator > threshold)
    final_rho = rho_values[final_rho_index]    
    final_theta = theta_angles[final_theta_index]
    
    polar_coordinates = np.vstack([final_rho, final_theta]).T

At the end, we extract the cells where the counters are above the threshold and the corresponding rho and theta so that in the end we will be able to reconstruct the lines in the image space.

And that’s it! We only need to plot our result.

Plotting

For plotting, we define these 2 functions.

#Funtion to add lines to an image
def draw_lines(img: np.ndarray, lines: np.ndarray, color: tp.List[int] = [0, 0, 255], thickness: int = 1) -> tp.Tuple[np.ndarray]:
    new_image = np.copy(img)
    empty_image = np.zeros(img.shape[:2])

    if len(lines.shape) == 1:
        lines = lines[None, ...]

    # Draw found lines
    for rho, theta in lines:
        x0 = polar2cartesian(rho, theta)
        direction = np.array([x0[1], -x0[0]])
        pt1 = np.round(x0 + 1000*direction).astype(int)
        pt2 = np.round(x0 - 1000*direction).astype(int)
        empty_image = cv.line(img=empty_image,pt1=pt1, pt2=pt2, color=255, thickness=thickness)

    # Keep lower part of each line until intersection
    mask_lines = empty_image != 0
    min_diff = np.inf
    max_line = 0
    for i in range(mask_lines.shape[0]):
        line = mask_lines[i]
        indices = np.argwhere(line)
        if indices[-1] - indices[0] < min_diff:
            min_diff = indices[-1] - indices[0]
            max_line = i

    mask_boundaries = np.zeros_like(empty_image)
    mask_boundaries[max_line:] = 1
    mask = (mask_lines * mask_boundaries).astype(bool)

    new_image[mask] = np.array(color)
    
    return new_image, mask

# Function to perform the conversion between polar and cartesian coordinates
def polar2cartesian(radius: np.ndarray, angle: np.ndarray, cv2_setup: bool = True) -> np.ndarray
    return radius * np.array([np.sin(angle), np.cos(angle)])

And we use them for plotting the results:

lines = hough_lines(best_canny_res, threshold, -np.pi/2, np.pi/2)

# Show the image with the lines found
lines_img, mask = draw_lines(image, lines)

win_name = "hough"
cv.namedWindow(win_name)
cv.imshow(win_name, lines_img)
cv.waitKey(0)
cv.destroyAllWindows()

Detected lines of the road lanes are marked in red

Conclusions

In this tutorial, we demystified the Hough Transform for line detection, breaking down its complexities into manageable steps. The implementation in Python, using OpenCV, highlighted the algorithm’s practical application.

Understanding the mechanics of the Hough Transform, from edge detection to parameter space transformation, contributes to a deeper grasp of computer vision concepts. While OpenCV provides a convenient solution, this hands-on approach fosters a more profound appreciation for the algorithm’s versatility.

Happy coding!