Radon Transform Line Detection

Question

I'm trying to detect lines in a grayscale image. For that purpose, I'm using Radon transform in MATLAB. An example of my m-file is like below. I can detect multiple lines using this code. I also draw lines using shift and rotation properties for lines. However, I didn't understand how to get the start and end points of the detecting lines after getting rho and theta values.

It is easy for Hough transform since there is a function called houghlines() that returns the list of the lines for the given peaks. Is there any function that i can use for Radon transform similar to this function?

    % Radon transform line detection algorithm
    clear all; close all;

    % Determine the path of the input image
    str_inputimg = '3_lines.png' ;

    % Read input image 
    I = imread(str_inputimg) ;

    % If the input image is RGB or indexed color, convert it to grayscale 
    img_colortype = getfield(imfinfo(str_inputimg), 'ColorType') ;
    switch img_colortype
       case 'truecolor'
         I = rgb2gray(I) ;
       case 'indexedcolor'
         I = ind2gray(I) ;
    end

    figure;
    subplot(2,2,1) ;
    imshow(I) ;
    title('Original Image') ;

    % Convert image to black white 
    %BW = edge(I,'Sobel');
    BW=im2bw(I,0.25) ;
    subplot(2,2,2) ;
    imshow(BW); 
    title('BW Image') ;

    % Radon transform 
    % Angle projections  
    theta = [0:179]' ;
    [R, rho] = radon(BW, theta) ;
    subplot(2,2,3) ;
    imshow(R, [], 'XData', theta, 'YData', rho, 'InitialMagnification', 'fit');
    xlabel('\theta'), ylabel('\rho');
    axis on, axis normal, hold on;

    % Detect the peaks of transform output  
    % Threshold value for peak detection
    threshold_val = ceil(0.3*max(R(:))) ;
    % Maximum nof peaks to identify in the image
    max_nofpeaks = 5 ;
    max_indexes = find(R(:)>threshold_val) ;
    max_values = R(max_indexes) ;
    [sorted_max, temp_indexes] = sort(max_values, 'descend') ;
    sorted_indexes = max_indexes(temp_indexes) ;

    % Get the first highest peaks for the sorted array
    if (length(sorted_max) <= max_nofpeaks)
        peak_values = sorted_max(1:end) ; 
        peak_indexes = sorted_indexes(1:end) ;
    else
        peak_values = sorted_max(1:max_nofpeaks) ;
        peak_indexes = sorted_indexes(1:max_nofpeaks) ;
    end
    [y, x]  = ind2sub(size(R), peak_indexes ) ;
    peaks = [rho(y) theta(x)] ;
    plot(peaks(:,2), peaks(:,1), 's', 'color','white');
    title('Radon Transform & Peaks') ;

    % Detected lines on the image
    subplot(2,2,4), imshow(I), title('Detected lines'), hold on

    x_center = floor(size(I, 2)/2) ;
    y_center = floor(size(I, 1)/2) ;
    for p=1:length(peaks)

        x_1 = [-x_center, x_center] ;
        y_1 = [0, 0] ;

        % Shift at first
        x_1_shifted = x_1 ;
        y_1_shifted = [y_1(1)-peaks(p,1), y_1(2)-peaks(p,1)] ;

        % Rotate 
        peaks(p,2) = 90 - peaks(p,2) ;
        t=peaks(p,2)*pi/180;
        rotation_mat = [ cos(t) -sin(t) ; sin(t) cos(t) ] ;
        x_y_rotated = rotation_mat*[x_1_shifted; y_1_shifted] ;
        x_rotated = x_y_rotated(1,:) ;
        y_rotated = x_y_rotated(2,:) ;
        plot( x_rotated+x_center, y_rotated+y_center, 'b', 'linewidth', 2 );
   end
   hold off;

Answer 1

There's a suggestion at math.SE which might help. Then there's a rather complicated-looking research paper "Sharp endpoint estimates for the X-ray transform and the Radon transform in finite fields", which appears just to show certain bounds on estimation accuracy.

From skimming other papers, it appears that it's a nontrivial problem. I suspect it may be simpler (if less accurate) to use some adaptation of a Sobel-operation to identify high gradient points along the discovered line, and claim those as endpoints.