matlab - Help understanding Hough transform

I'm trying to get a Hough transform to work in MATLAB, but I'm having problems. I have a really bad way of detecting peaks that needs to be fixed, but before that I need to be able to reverse the hough transform to create the lines again properly. This is the type of stuff I'm getting right now:

looks like its rotated by 90 degrees, but I'm not sure why. I'm not sure if it's my Hough space that's wrong, or if it's the way I de-Hough and draw the lines. Also could someone help improve my peak detection? The images used in the code are here

%% load a sample image; convert to grayscale; convert to binary

%create 'x' image (works well)
a = eye(255);
b = flipud(eye(255));
x = a + b;
x(128,128) = 1;


%image = rgb2gray(imread('up.png')) < 255;
%image = rgb2gray(imread('hexagon.png')) < 255;
%image = rgb2gray(imread('traingle.png')) < 255;
%%% these work
%image = x;
%image = a;
image = b;    

%% set up variables for hough transform

theta_sample_frequency = 0.01;                                             
[x, y] = size(image);
rho_limit = norm([x y]);                                                
rho = (-rho_limit:1:rho_limit);
theta = (0:theta_sample_frequency:pi);
num_thetas = numel(theta);
num_rhos = numel(rho);
hough_space = zeros(num_rhos, num_thetas);

%% perform hough transform

for xi = 1:x
    for yj = 1:y
        if image(xi, yj) == 1 
            for theta_index = 1:num_thetas
                th = theta(theta_index);
                r  = xi * cos(th) + yj * sin(th);
                rho_index = round(r + num_rhos/2);                      
                hough_space(rho_index, theta_index) = ...
                     hough_space(rho_index, theta_index) + 1;
            end

        end
    end
end  

%% show hough transform
subplot(1,2,1);
imagesc(theta, rho, hough_space);
title('Hough Transform');
xlabel('Theta (radians)');
ylabel('Rho (pixels)');

colormap('gray');

%% detect peaks in hough transform
r = [];
c = [];
[max_in_col, row_number] = max(hough_space);
[rows, cols] = size(image);
difference = 25;
thresh = max(max(hough_space)) - difference;
for i = 1:size(max_in_col, 2)

   if max_in_col(i) > thresh
       c(end + 1) = i;
       r(end + 1) = row_number(i);
   end
end

%% plot all the detected peaks on hough transform image
hold on;
plot(theta(c), rho(r),'rx');
hold off;



%% plot the detected line superimposed on the original image
subplot(1,2,2)
imagesc(image);
colormap(gray);
hold on;

for i = 1:size(c,2)
    th = theta(c(i));

    rh = rho(r(i));
    m = -(cos(th)/sin(th));
    b = rh/sin(th);
    x = 1:cols;
    plot(x, m*x+b);
    hold on;
end

Linked: How to do De-Houghing of a Hough transform'ed Image?

Answer

Firstly, Matlab has a Built in Hough Transform: no need to reinvent the wheel.

[H,T,R] = hough(BW,'RhoResolution',0.5,'Theta',-90:0.5:89.5);

Although your image doesn't necessarily require edge detection you could improve the processing time and effectiveness of the algorithm my using it. Your triangle has fat regions of white and black. Ideally, the triangle would be 1px thick marking the edges of the triangle. Use Canny Edge Detection

BW = edge(Image,'canny');

the result is \theta is in the range $-90 < \theta < 90$ degrees while your plot is $0 < \theta < 180$ (or $0 < \theta < \pi$) so subtract $90$ degrees ($\pi/2$) .

There is the potential for you to choose the wrong peak because there are neighboring peaks that may be larger in the accumulator matrix. While there are many algorithms here is one that I've seen used in Hough Transforms in the past:

1) Define a region shape (typically its square) 
2) Define an accumulator threshold  
3) Select one pixel in the accumulator matrix
4) If it is a peak (i.e., larger than neighboring values above a threshold)
       a) look at the points within the region shape.
       b) If there is a larger value
              this is not a peak
          Else
              this is a peak
 5) Move to next pixel in accumulator matrix.

Look into HoughLines for displaying the hough transform lines, results:

http://www.mathworks.com/help/toolbox/images/ref/houghlines.html

Effects of Using Canny Edge Detector

Edge Detection can potentially turn each side of the triangle into two lines.

The goal of canny edge detection is to produce maximally thin/narrow edges by using nonmaximal supression

Canny Edge Detection in a Nutshell (Source: Digital Image Processing, Gonazalez)

1) Smooth input Image using a Gaussian Filter
2) Compute the Gradient magnitude and angle (Sobel, Perwitt or robert cross filters)
3) Apply Nonmaxima suppression (this is where the thinning happens) 

   a) Figure out which direction the edge is
   b) If the edge's magnitude is smaller than one of its two neighbors in the direction of the edge
          set the edge point to zero
      Else
          leave it alone
4) Use double thresholding and connectivity analysis to detect and link edges