Implementing Hough transform line detection in 2D coordinate system

Question

I want to implement line detection in a simple coordinate system. I roughly followed an article about how to implement the Hough Transform, but the results I get are quite far off from what I want.

Given a 3 x 3 matrix like this:

X X X
X X X
- - -

I want to detect the line starting at 0,0 going to 2,0. I represent the coordinate system as a simple array of tuples, the first item in the tuple is x, the second is y, the third is the type of the point (canvas or line).

I thought it would be relatively easy to detect the line using Hough, because the edge detection is basically just a binary decision: either the tuple is of type line, or not.

I implemented the following program in Rust:

use std::f32;

extern crate nalgebra as na;
use na::DMatrix;

#[derive(Debug, PartialEq, Clone)]
enum Representation {
   Canvas,
   Line,
}

fn main () {
    let image_width = 3;
    let image_height = 3;

    let grid = vec![
        (0, 0, Representation::Line), (1, 0, Representation::Line), (2, 0, Representation::Line),
        (0, 1, Representation::Canvas), (1, 1, Representation::Canvas), (2, 1, Representation::Canvas),
        (0, 2, Representation::Canvas), (1, 2, Representation::Canvas), (2, 2, Representation::Canvas),
    ];

    //let tmp:f32 = (image_width as f32 * image_width as f32) + (image_height as f32 * image_height as f32);
    let max_line_length = 3;
    let mut accumulator = DMatrix::from_element(180, max_line_length as usize, 0);

    for y in 0..image_height {
        for x in 0..image_width {
            let coords_index = (y * image_width) + x;
            let coords = grid.get(coords_index as usize).unwrap();

            // check if coords is an edge
            if coords.2 == Representation::Line {
                for angle in 0..180 {
                    let r = (x as f32) * (angle as f32).cos() + (y as f32) * (angle as f32).sin();
                    let r_scaled = scale_between(r, 0.0, 2.0, -2.0, 2.0).round() as u32;

                    accumulator[(angle as usize, r_scaled as usize)] += 1;
                }
            }
        }
    }

    let threshold = 3;

    // z = angle
    for z in 0..180 {
        for r in 0..3 {
            let val = accumulator[(z as usize, r as usize)];

            if val < threshold {
                continue;
            }

            let px = (r as f32) * (z as f32).cos();
            let py = (r as f32) * (z as f32).sin();

            let p1_px = px + (max_line_length as f32) * (z as f32).cos();
            let p1_py = py + (max_line_length as f32) * (z as f32).sin();

            let p2_px = px - (max_line_length as f32) * (z as f32).cos();
            let p2_py = px - (max_line_length as f32) * (z as f32).cos();

            println!("Found lines from {}/{} to {}/{}", p1_px.ceil(), p1_py.ceil(), p2_px.ceil(), p2_py.ceil());
        }
    }
}

fn scale_between(unscaled_num: f32, min_allowed: f32, max_allowed: f32, min: f32, max: f32) -> f32 {
    (max_allowed - min_allowed) * (unscaled_num - min) / (max - min) + min_allowed
}

The result is something like:

Found lines from -1/4 to 1/1
Found lines from 2/4 to 0/0
Found lines from 2/-3 to 0/0
Found lines from -1/4 to 1/1
Found lines from 1/-3 to 0/0
Found lines from 0/4 to 1/1
...

Which is actually quite a lot, given that I only want to detect a single line. My implementation clearly is wrong, but I don't know where to look, my maths-fu is not high enough to debug further.

I think the first part, the actual Hough Transform, seems kind of correct, because the linked article says:

for each image point p 
{
  if (p is part of an edge)
  {
    for each possible angle
    {
     r = x * cos(angle) + y * sin(angle);
     houghMatrix[angle][r]++;
    }
  }
}

I'm stuck at mapping and filtering, which is according to the article:

1. Each point in Hough space is given by angle a and distance r. Using these values, one single point p(x,y) of the line can be calculated by px = r * cos(angle) py = r * sin(angle).

2. The maximum length of a line is restricted by sqrt(imagewidth2 + imageheight2).

3. The point p, the angle a of the line and the maximum line length 'maxLength' can be used to calculate two other points of the line. The maximum length here ensures that both points to be calculated are lying outside of the actual image, resulting in the fact that if a line is drawn between these two points, the line goes from image border to image border in any case and is never cropped somewhere inside the image.

4. These two points p1 and p2 are calculated by: p1_x = px + maxLength * cos(angle); p1_y = py + maxLength * sin(angle); p2_x = px - maxLength * cos(angle); p2_y = py - maxLength * sin(angle);

5. ...

EDIT

Updated version that takes the image size into account, as suggested by @RaymoAisla

use std::f32;

extern crate nalgebra as na;
use na::DMatrix;

fn main () {
    let image_width = 3;
    let image_height = 3;

    let mut grid = DMatrix::from_element(image_width as usize, image_height as usize, 0);
    grid[(0, 0)] = 1;
    grid[(1, 0)] = 1;
    grid[(2, 0)] = 1;

    let accu_width = 7;
    let accu_height = 3;
    let max_line_length = 3;

    let mut accumulator = DMatrix::from_element(accu_width as usize, accu_height as usize, 0);


    for y in 0..image_height {
        for x in 0..image_width {
            let coords = (x, y);
            let is_edge = grid[coords] == 1;

            if !is_edge {
                continue;
            }

            for i in 0..7 {
                let angle = i * 30;

                let r = (x as f32) * (angle as f32).cos() + (y as f32) * (angle as f32).sin();
                let r_scaled = scale_between(r, 0.0, 2.0, -2.0, 2.0).round() as u32;

                accumulator[(i as usize, r_scaled as usize)] += 1;

                println!("angle: {}, r: {}, r_scaled: {}", angle, r, r_scaled);
            }
        }
    }

    let threshold = 3;

    // z = angle index
    for z in 0..7 {
        for r in 0..3 {
            let val = accumulator[(z as usize, r as usize)];

            if val < threshold {
                continue;
            }

            let px = (r as f32) * (z as f32).cos();
            let py = (r as f32) * (z as f32).sin();

            let p1_px = px + (max_line_length as f32) * (z as f32).cos();
            let p1_py = py + (max_line_length as f32) * (z as f32).sin();

            let p2_px = px - (max_line_length as f32) * (z as f32).cos();
            let p2_py = px - (max_line_length as f32) * (z as f32).cos();

            println!("Found lines from {}/{} to {}/{} - val: {}", p1_px.ceil(), p1_py.ceil(), p2_px.ceil(), p2_py.ceil(), val);
        }
    }
}

fn scale_between(unscaled_num: f32, min_allowed: f32, max_allowed: f32, min: f32, max: f32) -> f32 {
    (max_allowed - min_allowed) * (unscaled_num - min) / (max - min) + min_allowed
}

The reported output is now:

angle: 0, r: 0, r_scaled: 1
angle: 30, r: 0, r_scaled: 1
angle: 60, r: 0, r_scaled: 1
angle: 90, r: 0, r_scaled: 1
angle: 120, r: 0, r_scaled: 1
angle: 150, r: 0, r_scaled: 1
angle: 180, r: 0, r_scaled: 1
...
Found lines from 3/4 to -1/-1
Found lines from -3/1 to 2/2 

I plotted the lines on a coordinate system, the lines are very far off from the line that I would expect. I wonder if the conversion back to points is still off.

Answer 1

Your angles are in degrees rather than radians!

Rust, like all other programming languages, uses radians for its trigonometry functions. Running

let ang_d = 30.0;
let ang_r = ang_d * 3.1415926 / 180.0;
println!("sin(30) {} sin(30*pi/180) {}", (ang_d as f32).sin(), (ang_r as f32).sin());

gives the results

sin(30) -0.9880316 sin(30*pi/180) 0.5

You need to convert all your angles to radians before calling cos and sin.

In the first loop I've got

let angle = (i as f32) * 30.0 * 3.1415926 / 180.0;
let r = (x as f32) * (angle as f32).cos() + (y as f32) * (angle as f32).sin();

and in the second where you calculate the points on the lines

let ang = (z as f32) * 30.0 * 3.1415926 / 180.0;
let px = (r as f32) * (ang as f32).cos();
let py = (r as f32) * (ang as f32).sin();
let p1_px = px + (max_line_length as f32) * (ang as f32).cos();          
let p1_py = py + (max_line_length as f32) * (ang as f32).sin();
let p2_px = px - (max_line_length as f32) * (ang as f32).cos();
let p2_py = px - (max_line_length as f32) * (ang as f32).cos();

My Rust is rusty (actually non-existent) so there is nicer ways of doing the conversion and there should be a constant with the exact value of pi somewhere.

Answer 2

Hough Transform principle is to search all the lines passing through each considered point, and counting these lines occurrences thanks to the accumulator.

However, we cannot determine all these lines because their number is infinite. Moreover the image is discretized, so calculating all lines does not make sense.

And the problem comes from this discretization. The angle discretization needs to be relevant with the image size. Here, calculating the radius for 180 angles is overcalculation, because the image only make 9 pixels, and the possible angles for any line in this image are restricted to a dozen value.

So here, for the first point (0,0), for each angle, the associated radius is r = 0

For the second (1,0), the associated radius is r = cos(angle)

For the third (2,0), the associated radius is r = 2 cos(angle)

With the rounding, numerous values will have an associated radius of 0 for the same angle, and it cause over-detection. Discretization causes a spreading of Hough Accumulator.

So the radius and the angle discretization needs to be calculated depending on the image size. Here, a 30° step, so a 7*3 accumulator will be sufficient to detect a line.