关于MATLAB中Hough变换提取直线

，本人正在利用Hough变换进行直线提取，现在已经对图像进行了分割和数学形态学的处理，可是在Hough变换提取直线是怎么都不对，求一个正确的Hough 变换的程序!急啊！先谢谢各位了！附上经数学形态学处理后的照片。目的是把道路的两条边界线提取出来。

• 
  
  1. clear all
     
  2. close all
     
  3. clc
     
  4. I=imread('2220.jpg');
     
  5. BW=im2bw(I);
     
  6. BW=edge(BW,'canny');
     
  7. % I=imdilate(I,strel('disk',8));
     
  8. % I=imerode(I,strel('disk',8));
     
  9. % BW=bwmorph(I,'thin',20);
     
  10. %%
      
  11. 
      
  12. [H,T,R] = hough(BW);
      
  13. imshow(H,[],'XData',T,'YData',R,...
      
  14. 'InitialMagnification','fit');
      
  15. xlabel('\theta'), ylabel('\rho');
      
  16. axis on, axis normal, hold on;
      
  17. P = houghpeaks(H,10,'threshold',ceil(0.3*max(H(:))));
      
  18. x = T(P(:,2)); y = R(P(:,1));
      
  19. plot(x,y,'s','color','white');
      
  20. lines = houghlines(BW,T,R,P,'FillGap',5,'MinLength',7);
      
  21. figure,imshow(BW), hold on
      
  22. max_len = 0;
      
  23. %%
      
  24. 
      
  25. for k = 1:length(lines)
      
  26. xy = [lines(k).point1; lines(k).point2];
      
  27. plot(xy(:,1),xy(:,2),'LineWidth',2,'Color','green');
      
  28. 
      
  29. % Plot beginnings and ends of lines
      
  30. plot(xy(1,1),xy(1,2),'x','LineWidth',2,'Color','yellow');
      
  31. plot(xy(2,1),xy(2,2),'x','LineWidth',2,'Color','red');
      
  32. 
      
  33. % Determine the endpoints of the longest line segment
      
  34. len = norm(lines(k).point1 - lines(k).point2);
      
  35. Len(k)=len
      
  36. if ( len > max_len)
      
  37. max_len = len
      
  38. xy_long = xy
      
  39. end
      
  40. end
      
  41. 
      
  42. % highlight the longest line segment
      
  43. plot(xy_long(:,1),xy_long(:,2),'LineWidth',2,'Color','blue');
      
  44. 
      
  45. %%
      
  46. [L1 Index1]=max(Len(:))
      
  47. Len(Index1)=0
      
  48. [L2 Index2]=max(Len(:))
      
  49. 
      
  50. %%
      
  51. %
      
  52. x1=[lines(Index1).point1(1) lines(Index1).point2(1)]
      
  53. y1=[lines(Index1).point1(2) lines(Index1).point2(2)]
      
  54. x2=[lines(Index2).point1(1) lines(Index2).point2(1)]
      
  55. y2=[lines(Index2).point1(2) lines(Index2).point2(2)]
      
  56. 
      
  57. %%
      
  58. K1=(lines(Index1).point1(2)-lines(Index1).point2(2))/(lines(Index1).point1(1)-lines(Index1).point2(1))
      
  59. K2=(lines(Index2).point1(2)-lines(Index2).point2(2))/(lines(Index2).point1(1)-lines(Index2).point2(1))
      
  60. %%
      
  61. hold on
      
  62. [m,n] = size(BW); % 尺寸
      
  63. BW1=zeros(m,n);
      
  64. b1=y1(1)-K1*x1(1)
      
  65. b2=y2(1)-K2*x2(1)
      
  66. for x=1:n
      
  67. for y=1:m
      
  68. if y==round(K1*x+b1)|y==round(K2*x+b2)
      
  69. BW1(y,x)=1;
      
  70. 
      
  71. end
      
  72. 
      
  73. end
      
  74. end
      
  75. for x=1:n
      
  76. for y=1:m
      
  77. 
      
  78. if ceil(K1*x+b1)==ceil(K2*x+b2)
      
  79. y1=round(K1*x+b1)
      
  80. BW1(1:y1-1,:)=0;
      
  81. end
      
  82. end
      
  83. end
      
  84. figure,imshow(BW1)