How do I transform the values of an accumulator [Hough Transformation] back to a line on a canvas?

Question

I am trying to detect the lines within an image using the Hough Transformation. Therefore I first create the accumulator like this:

from math import hypot, pi, cos, sin
from PIL import Image
import numpy as np
import cv2 as cv
import math

def hough(img):

    thetaAxisSize = 460 #Width of the hough space image
    rAxisSize = 360 #Height of the hough space image
    rAxisSize= int(rAxisSize/2)*2 #we make sure that this number is even

    img = im.load()
    w, h = im.size

    houghed_img = Image.new("L", (thetaAxisSize, rAxisSize), 0) #legt Bildgroesse fest
    pixel_houghed_img = houghed_img.load()

    max_radius = hypot(w, h)
    d_theta = pi / thetaAxisSize
    d_rho = max_radius / (rAxisSize/2) 


    #Accumulator
    for x in range(0, w):
        for y in range(0, h):

            treshold = 255
            col = img[x, y]
            if col >= treshold: #determines for each pixel at (x,y) if there is enough evidence of a straight line at that pixel.

                for vx in range(0, thetaAxisSize):
                    theta = d_theta * vx #angle between the x axis and the line connecting the origin with that closest point.
                    rho = x*cos(theta) + y*sin(theta) #distance from the origin to the closest point on the straight line
                    vy = rAxisSize/2 + int(rho/d_rho+0.5) #Berechne Y-Werte im hough space image
                    pixel_houghed_img[vx, vy] += 1 #voting

    return houghed_imgcode here

And then call the function like this:

im = Image.open("img3.pgm").convert("L")
houghed_img = hough(im)
houghed_img.save("ho.bmp")
houghed_img.show()

The result seems to be okay:

So here comes the problem. I know want to find the top 3 highest values in the hough space and transform it back to 3 lines. The highest values should be the strongest lines.

Therefore I am first looking for the highest values within the pixel array and take the X and Y values of the maxima I found. From my understading this X and Y values are my rho and theta. I finding the maxima like this:

def find_maxima(houghed_img):

    w, h = houghed_img.size
    max_radius = hypot(w, h)
    pixel_houghed_img = houghed_img.load()
    max1, max2, max3 = 0, 0, 0
    x1position, x2position, x3position = 0, 0, 0
    y1position, y2position, y3position = 0, 0, 0
    rho1, rho2, rho3 = 0, 0, 0
    theta1, theta2, theta3 = 0, 0, 0

    for x in range(1, w):
        for y in range(1, h):
            value = pixel_houghed_img[x, y]

            if(value > max1):

                max1 = value
                x1position = x
                y1position = y
                rho1 = x
                theta1 = y

            elif(value > max2):

                max2 = value
                x2position = x
                x3position = y
                rho2 = x
                theta2 = y

            elif(value > max3):

                max3 = value
                x3position = x
                y3position = y
                rho3 = x
                theta3 = y

    print('max', max1, max2, max3)
    print('rho', rho1, rho2, rho3)
    print('theta', theta1, theta2, theta3)

    # Results of the print:
    # ('max', 255, 255, 255)
    # ('rho', 1, 1, 1)
    # ('theta', 183, 184, 186)
    return rho1, theta1, rho2, theta2, rho3, theta3

And now I want to use this rho and theta values to draw the detected lines. I am doing this with the following code:

img_copy = np.ones(im.size)

rho1, theta1, rho2, theta2, rho3, theta3 = find_maxima(houghed_img)

a1 = math.cos(theta1)
b1 = math.sin(theta1)
x01 = a1 * rho1
y01 = b1 * rho1
pt11 = (int(x01 + 1000*(-b1)), int(y01 + 1000*(a1)))
pt21 = (int(x01 - 1000*(-b1)), int(y01 - 1000*(a1)))
cv.line(img_copy, pt11, pt21, (0,0,255), 3, cv.LINE_AA)

a2 = math.cos(theta2)
b2 = math.sin(theta2)
x02 = a2 * rho2
y02 = b2 * rho2
pt12 = (int(x02 + 1000*(-b2)), int(y02 + 1000*(a2)))
pt22 = (int(x02 - 1000*(-b2)), int(y02 - 1000*(a2)))
cv.line(img_copy, pt12, pt22, (0,0,255), 3, cv.LINE_AA)

a3 = math.cos(theta3)
b3 = math.sin(theta3)
x03 = a3 * rho3
y03 = b3 * rho3
pt13 = (int(x03 + 1000*(-b3)), int(y03 + 1000*(a3)))
pt23 = (int(x03 - 1000*(-b3)), int(y03 - 1000*(a3)))
cv.line(img_copy, pt13, pt23, (0,0,255), 3, cv.LINE_AA)

cv.imshow('lines', img_copy)
cv.waitKey(0)
cv.destroyAllWindows()

However, the result seems to be wrong:

So my assuption is that I either do something wrong when I declare the rho and theta values in the find_maxima() function, meaning that something is wrong with this:

   max1 = value
   x1position = x
   y1position = y
   rho1 = x
   theta1 = y

OR that I am doing something wrong when translating the rho and theta value back to a line.

I would be very thankful if someone can help me with that!

Edit1: As request please finde the original Image where I want to finde the lines from below:

Edit2: Thanks to the input of @Alessandro Jacopson and @Cris Luegno I was able to make some changes that definitely give me some hope!

In my def hough(img): I was setting the threshold to 255, which means that I only voted for white pixels, which is wrong since I want to look at the black pixels, since these pixels will indicate lines and not the white background of my image. So the calculation of the accumlator in def hough(img): looks like this now:

#Accumulator
    for x in range(0, w):
        for y in range(0, h):

            treshold = 0
            col = img[x, y]
            if col <= treshold: #determines for each pixel at (x,y) if there is enough evidence of a straight line at that pixel.

                for vx in range(0, thetaAxisSize):
                    theta = d_theta * vx #angle between the x axis and the line connecting the origin with that closest point.
                    rho = x*cos(theta) + y*sin(theta) #distance from the origin to the closest point on the straight line
                    vy = rAxisSize/2 + int(rho/d_rho+0.5) #Berechne Y-Werte im hough space image
                    pixel_houghed_img[vx, vy] += 1 #voting

    return houghed_img

This leads to the following Accumulator and the following rho and thea values, when using the find_maxima() function:

# Results of the prints: (now top 8 instead of top 3)
# ('max', 155, 144, 142, 119, 119, 104, 103, 98)
# ('rho', 120, 264, 157, 121, 119, 198, 197, 197)
# ('theta', 416, 31, 458, 414, 417, 288, 291, 292)

The Lines that I can draw from this values look like this:

So this results are much more better but something seems to be still wrong. I have a strong suspicion that still something is wrong here:

for x in range(1, w):
    for y in range(1, h):
        value = pixel_houghed_img[x, y]

        if(value > max1):

            max1 = value
            x1position = x
            y1position = y
            rho1 = value
            theta1 = x

Here I am setting rho and theta equals [0...w] respectively [0...h]. I think that this is wrong since in the hough space values of X and why Y are not 0, 1,2,3... since we are in a another space. So I assume, that I have to multiply X and Y with something to bring them back in hough space. But this is just an assumption, maybe you guys can think of something else?

Again thank you very much to Alessandro and Cris for helping me out here!

Edit3: Working Code, thanks to @Cris Luengo

from math import hypot, pi, cos, sin
from PIL import Image
import numpy as np
import cv2 as cv
import math

def hough(img):

    img = im.load()
    w, h = im.size

    thetaAxisSize = w #Width of the hough space image
    rAxisSize = h #Height of the hough space image
    rAxisSize= int(rAxisSize/2)*2 #we make sure that this number is even

    houghed_img = Image.new("L", (thetaAxisSize, rAxisSize), 0) #legt Bildgroesse fest
    pixel_houghed_img = houghed_img.load()

    max_radius = hypot(w, h)
    d_theta = pi / thetaAxisSize
    d_rho = max_radius / (rAxisSize/2) 

    #Accumulator
    for x in range(0, w):
        for y in range(0, h):

            treshold = 0
            col = img[x, y]
            if col <= treshold: #determines for each pixel at (x,y) if there is enough evidence of a straight line at that pixel.

                for vx in range(0, thetaAxisSize):
                    theta = d_theta * vx #angle between the x axis and the line connecting the origin with that closest point.
                    rho = x*cos(theta) + y*sin(theta) #distance from the origin to the closest point on the straight line
                    vy = rAxisSize/2 + int(rho/d_rho+0.5) #Berechne Y-Werte im hough space image
                    pixel_houghed_img[vx, vy] += 1 #voting

    return houghed_img, rAxisSize, d_rho, d_theta

def find_maxima(houghed_img, rAxisSize, d_rho, d_theta):

    w, h = houghed_img.size
    pixel_houghed_img = houghed_img.load()
    maxNumbers = 9
    ignoreRadius = 10
    maxima = [0] * maxNumbers
    rhos = [0] * maxNumbers
    thetas = [0] * maxNumbers

    for u in range(0, maxNumbers):

        print('u:', u)
        value = 0 
        xposition = 0
        yposition = 0

        #find maxima in the image
        for x in range(0, w):
            for y in range(0, h):

                if(pixel_houghed_img[x,y] > value):

                    value = pixel_houghed_img[x, y]
                    xposition = x
                    yposition = y

        #Save Maxima, rhos and thetas
        maxima[u] = value
        rhos[u] = (yposition - rAxisSize/2) * d_rho
        thetas[u] = xposition * d_theta

        pixel_houghed_img[xposition, yposition] = 0

        #Delete the values around the found maxima
        radius = ignoreRadius

        for vx2 in range (-radius, radius): #checks the values around the center
            for vy2 in range (-radius, radius): #checks the values around the center
                x2 = xposition + vx2 #sets the spectated position on the shifted value 
                y2 = yposition + vy2

                if not(x2 < 0 or x2 >= w):
                    if not(y2 < 0 or y2 >= h):

                        pixel_houghed_img[x2, y2] = 0
                        print(pixel_houghed_img[x2, y2])

    print('max', maxima)
    print('rho', rhos)
    print('theta', thetas)

    return maxima, rhos, thetas

im = Image.open("img5.pgm").convert("L")
houghed_img, rAxisSize, d_rho, d_theta = hough(im)
houghed_img.save("houghspace.bmp")
houghed_img.show()

img_copy = np.ones(im.size)

maxima, rhos, thetas = find_maxima(houghed_img, rAxisSize, d_rho, d_theta)

for t in range(0, len(maxima)):
    a = math.cos(thetas[t])
    b = math.sin(thetas[t])
    x = a * rhos[t]
    y = b * rhos[t]
    pt1 = (int(x + 1000*(-b)), int(y + 1000*(a)))
    pt2 = (int(x - 1000*(-b)), int(y - 1000*(a)))
    cv.line(img_copy, pt1, pt2, (0,0,255), 3, cv.LINE_AA)

cv.imshow('lines', img_copy)
cv.waitKey(0)
cv.destroyAllWindows()

Original Image:

Accumulator:

Successful Line Detection:

Your 3 maxima having a value of 255 is fishy. Is you accumulator maybe of type uint8? Are you saturating the accumulator bins? — Cris Luengo
Hey Cris, I knew that I could count on you again! Thank you very much! You made me think with your input and I think that I found the mistake in my code regarding the 255 values. However, something is still wrong... Please see my updated question. — Leonard Michalas

Alessandro Jacopson Alessandro Jacopson · Accepted Answer · 2018-06-24T17:30:58

First of all, following How to create a Minimal, Complete, and Verifiable example you should post or give a link to your image img3.pgm, if possible.

Then, you wrote that:

# Results of the print:
# ('max', 255, 255, 255)
# ('rho', 1, 1, 1)
# ('theta', 183, 184, 186)

so rho is the same for the three lines and theta is not so different varying between 183 and 186; so the three lines are almost equal each other and this fact does not depend on the method you use to get the line equation and draw it.

According to the tutorial Hough Line Transform it seems to me that your method for finding two points on a line is correct. That's is what the tutorial is suggesting and it seems to me equivalent to your code:

lines = cv2.HoughLines(edges,1,np.pi/180,200)
for rho,theta in lines[0]:
    a = np.cos(theta)
    b = np.sin(theta)
    x0 = a*rho
    y0 = b*rho
    x1 = int(x0 + 1000*(-b))
    y1 = int(y0 + 1000*(a))
    x2 = int(x0 - 1000*(-b))
    y2 = int(y0 - 1000*(a))

    cv2.line(img,(x1,y1),(x2,y2),(0,0,255),2)

I suspect the peak finding algorithm may not be correct. Your peak finding algorithm finds the location of the largest peak and then the two locations very close to that maximum.

For the sake of simplicity see what happens in just one dimension, a peak finding algorithm is expected to find three peak locations at x=-1, x=0 and x=1 and the peak values should be close to .25, .5 and 1.

import numpy as np
import matplotlib.pyplot as plt


x = np.linspace(-2, 2, 1000)
y = np.exp(-(x-1)**2/0.01)+.5*np.exp(-(x)**2/0.01)+.25*np.exp(-(x+1)**2/0.01)

max1, max2, max3 = 0, 0, 0
m1 = np.zeros(1000)
m2 = np.zeros(1000)
m3 = np.zeros(1000)
x1position, x2position, x3position = 0, 0, 0
for i in range(0,1000):
    value = y[i]

    if(value > max1):

        max1 = value
        x1position = x[i]

    elif(value > max2):

        max2 = value
        x2position = x[i]

    elif(value > max3):

        max3 = value
        x3position = x[i]

    m1[i] = max1
    m2[i] = max2
    m3[i] = max3



print('xposition',x1position, x2position, x3position )
print('max', max1, max2, max3)

plt.figure()
plt.subplot(4,1,1)
plt.plot(x, y)
plt.ylabel('$y$')
plt.subplot(4,1,2)
plt.plot(x, m1)
plt.ylabel('$max_1$')
plt.subplot(4,1,3)
plt.plot(x, m2)
plt.ylabel('$max_2$')
plt.subplot(4,1,4)
plt.plot(x, m3)
plt.xlabel('$x$')
plt.ylabel('$max_3$')
plt.show()

the output is

('xposition', 0.99899899899899891, 1.0030030030030028, 1.0070070070070072)
('max', 0.99989980471948192, 0.99909860379824966, 0.99510221871862647)

and it is not what expected.

Here you have a visual trace of the program:

To detect multiple peaks in a 2D field you should have a look for example at this Peak detection in a 2D array

How do I transform the values of an accumulator [Hough Transformation] back to a line on a canvas?

2 Answers