I = imread('P.tif');

BW = hEdge(I);

figure(1);subplot(1,2,2)
imshow(I)

[acc,T,R] = hough(BW, 'ThetaResolution', 1, 'RhoResolution', 1);

figure(1);subplot(1,2,1)
imshow(log(acc+1), [], 'XData', T, 'YData', R, 'InitialMagnification','fit');
xlabel('\theta'), ylabel('\rho'); axis on, axis normal;

P  = houghpeaks(acc, 18, 'Threshold', 0.2*max(acc(:)) );

[xs ys] = size(BW);

theta = T(P(:,2));   
s     = R(P(:,1));            

for i=1:size(P,1)
    figure(1);subplot(1,2,1)
    hold on
    plot(theta(i), s(i), 's', 'color', 'red');
    hold off
  
    figure(1);subplot(1,2,2);
    hold on
    ht_plot2( theta(i), s(i), xs, ys );
    hold off
end


figure(2); imshow(BW, [])