円や曲線のhough変換

下記の構文によって直線の検出を行い、交点を見つけ出すことができるのですが、円や曲線においては検出ができないので改良方法について意見をいただけると幸いです。

clear all　 %ワークスペースから全てのオブジェクトの消去

format long E %小数点以下表示

tic %計測始め(処理にかかる時間を計り始める)

% I=imread('0117_5_1.png'); %画像データの読み込み

% figure(1);imshow(I)

% BW=im2bw(I ,0.65); %画像を二値化

I=imread('mc2.jpg'); %画像データの読み込み

I2=rgb2gray(I);

BW=edge(I2,'Canny'); %画像のエッジ化

figure(1);imshow(I);

hold on

[M,N]=size(BW);

i=1;

j=1;

MMM=[]; %空の配列作成

MMM2=[];

for n=1:N

for m=1:M

if BW(m, n)==1; MMM(i, 1)=n; MMM(i, 2)=m;

i=i+1;

end

end

for m=(round(M/2)):M

if BW(m, n)==1; MMM2(j, 1)=n; MMM2(j, 2)=m;

j=j+1;

end

end

end

Data_N = size(MMM2); %データの行列の大きさ

X = MMM(:,1); %データの1列目をXとする

Y = MMM(:,2); %データの2列目をYとする

X2 = MMM2(:,1); %データの1列目をXとする

Y2 = MMM2(:,2); %データの2列目をYとする

figure(2)

subplot(1,2,1);imshow(BW);

subplot(1,2,2);plot(X2,Y2,'.')

set(gca,'YDir','reverse')

xlim([0 N])

ylim([0 M])

daspect([1 1 1])

% x = DATA(:,1); %データの1列目をXとする

% y = DATA(:,2); %データの2列目をYとする

% X=x*0.01;

% Y=y*0.01;

d_pi = 0.01; %傾きπの分割間隔

Pai = 0:d_pi:pi; %傾きπ(0～3.14)

Pai_N = size(Pai); %傾きπの行列の大きさ

i = Data_N(1,1); %i=1～6302(DataNumber)

j = Pai_N(2); %j=1～315(Pai)

hough_1 = zeros(i , j); %各ボックスフィルタの初期化

for j = 1:Pai_N(2)

for i = 1:Data_N(1,1)

rou = X2(i)*cos(Pai(j)) +Y2(i)*sin(Pai(j)); %ρの計算

rou2 = round(rou); %ρの値を四捨五入(整数値化)

rou3 = rou2-rou;

rou4 = rou2 + 301;

if rou3>0

hough_1(rou4,j) = hough_1(rou4,j) + (1-abs(rou3));

hough_1(rou4-1,j) = hough_1(rou4-1,j) + (abs(rou3));

else

hough_1(rou4,j) = hough_1(rou4,j) + (1-abs(rou3));

hough_1(rou4+1,j) = hough_1(rou4+1,j) + abs(rou3);

end

% hough_1(rou2, j) = hough_1(rou2, j) + 1; %Hough空間において、曲線交わるたびにその座標に重みを足していく(重み付け)

end

end

% グラフのhyouzi

figure(3) %二つのヒストグラムの比較

%ax1=subplot(2,1,1); %同時に二つのグラフを表示(上側に表示)

plot(Pai,max(hough_1)) %各θの最大の重さを示すヒストグラム

title('Hough Transform') %タイトル「Max Histogram」

xlabel('⊿θ(rad)','FontSize',20) %x軸のラベル

ylabel('Frequency','FontSize',20) %y軸のラベル

set(gca,'FontSize',20); %フォントの大きさ

xlim([0 3.14]) %x軸の目盛の範囲の指定

%y軸の目盛の範囲の指定

% ax2=subplot(2,1,2); %同時に二つのグラフを表示(下側に表示)

% plot(ax2,Pai,max(fm2_3)) %ボックスフィルタをかけたあとのヒストグラム

% title(ax2,'After') %タイトル「m推定」

% xlabel('⊿θ(rad)','FontSize',20) %x軸のラベル

% ylabel('Frequency','FontSize',20) %y軸のラベル

% set(gca,'FontSize',20); %フォントの大きさ

% xlim([0 3.14]) %x軸の目盛の範囲の指定

figure(4)

findpeaks(max(hough_1),Pai,'NPeaks',10,'SortStr','descend');%,'MinPeakDistance',0.4); %Class20000以上かつ0.5rad間隔という条件ででピークを検出

pks=findpeaks(max(hough_1),Pai,'NPeaks',10,'SortStr','descend');%,'MinPeakDistance',0.4); %ピーク値を表示(コマンドウィンドウに'pks'と入力)

xlabel('⊿θ(rad)','FontSize',20) %x軸のラベル

ylabel('Class','FontSize',20) %y軸のラベル

set(gca,'FontSize',20); %フォントの大きさ

hold on

plot([0.17 0.17],[-100000 400000],'r--',[1.4 1.4],[-100000 400000],'r--',[1.74 1.74],[-100000 400000],'m--',[2.97 2.97],[-100000 400000],'m--')

figure(5)

plot(X,Y,'.')

hold on

set(gca,'YDir','reverse')

xlim([0 N])

ylim([0 M])

%axis tight

xlabel('Width (pixel)','FontSize',20) %x軸のラベル

ylabel('Height (pixel)','FontSize',20) %y軸のラベル

set(gca,'FontSize',20) %フォントの大きさ

solx = zeros();

solxr = zeros();

srou = zeros();

stheta = zeros();

j = 1;

A = zeros();

B = zeros();

for i=1:numel(pks)

solx(i)=find(max(hough_1)==pks(:,i)); %ピーク値に対応するPaiの座標抽出

solxr(i)=find(hough_1(:,solx(i))==pks(:,i)); %ピーク値に対応するrou2の座標抽出

%if solx(i) ~= 140:175

if 18<solx(i) && solx(i)<140 || 175<solx(i) && solx(i)< 298

%srou(j)=(solxr(i)-4000)/500; %得られたρの値(谷部)+離散化から元に戻す(+8して*500をしたため)

srou(j)=(solxr(i)-301);

stheta(j)=Pai(solx(i)); %得られたθの値(谷部)

XX = 0:1:X(end);

A(j)=-(cos(stheta(j))/sin(stheta(j)));

B(j)=srou(j)/sin(stheta(j));

Z = A(j)*XX + B(j); %直線

j=j+1;

figure(5)

plot(XX,Z,'LineWidth',2)

end

end

c_x=zeros();

c_y=zeros();

i=1;

for p=1:(j-1)

for q=1:(j-1)

if A(p)<0 && A(q)>0 || A(p)>0 && A(q)<0

c_x(i) = (B(p)-B(q))/(A(q)-A(p));

c_y(i) = (A(q)*B(p)-A(p)*B(q))/(A(q)-A(p));

figure(5)

plot(c_x(i),c_y(i),'ro')

i=i+1;

end

end

end

t=10;

p = 200/t;

q = 300/t;

koten = zeros(p , q);

fm1_1 = zeros(p , q);

fm1_2 = zeros(p , q);

fm1_3 = zeros(p , q);

fm2_1 = zeros(p , q);

fm2_2 = zeros(p , q); %各ボックスフィルタの初期化

fm2_3 = zeros(p , q);

for i = 1:numel(c_x)

if 0<c_x(1,i) && c_x(1,i)<200 && 0<c_y(1,i) && c_y(1,i)<200

koten( fix(c_y(1,i)/t)+1 , fix(c_x(1,i)/t)+1 ) = koten( fix(c_y(1,i)/t)+1 , fix(c_x(1,i)/t)+1 ) + 1;

end

end

m = 1;

%一方向目(θ方向)

for j=m+1:q-m %ボックスフィルタ1回目

for i=m+1:p-m

fm1_1(i,j)=fm1_1(i,j)+koten(i,j); %初期位置の指定

for v=1:m

fm1_1(i,j)=fm1_1(i,j)+koten(i,j+v)+koten(i,j-v); %両隣の配列を参照していく

end

end

end

for j=m+1:q-m %ボックスフィルタ2回目

for i=m+1:p-m

fm1_2(i,j)=fm1_2(i,j)+fm1_1(i,j);

for v=1:m

fm1_2(i,j)=fm1_2(i,j)+fm1_1(i,j+v)+fm1_1(i,j-v);

end

end

end

for j=m+1:q-m %ボックスフィルタ3回目

for i=m+1:p-m

fm1_3(i,j)=fm1_3(i,j)+fm1_2(i,j);

for v=1:m

fm1_3(i,j)=fm1_3(i,j)+fm1_2(i,j+v)+fm1_2(i,j-v);

end

end

end

%2方向目(ρ方向)

for j=m+1:q-m %ボックスフィルタ1回目

for i=m+1:p-m

fm2_1(i,j)=fm2_1(i,j)+fm1_3(i,j);

for v=1:m

fm2_1(i,j)=fm2_1(i,j)+fm1_3(i+v,j)+fm1_3(i-v,j);

end

end

end

for j=m+1:q-m %ボックスフィルタ2回目

for i=m+1:p-m

fm2_2(i,j)=fm2_2(i,j)+fm2_1(i,j);

for v=1:m

fm2_2(i,j)=fm2_2(i,j)+fm2_1(i+v,j)+fm2_1(i-v,j);

end

end

end

for j=m+1:q-m %ボックスフィルタ3回目

for i=m+1:p-m

fm2_3(i,j)=fm2_3(i,j)+fm2_2(i,j);

for v=1:m

fm2_3(i,j)=fm2_3(i,j)+fm2_2(i+v,j)+fm2_2(i-v,j);

end

end

end

[M,N] = max(fm2_3(:)); %Mは最大要素、Nは最大要素を含むfm2_3(:)のインデックス

[gyo,retsu] = ind2sub(size(fm2_3),N); %最大要素に対応する行と列のインデックスを抽出

figure(1)

r = t/2;

xc = t*retsu-(t/2);

yc = t*gyo-(t/2);

theta = linspace(0,2*pi);

x = r*cos(theta) + xc;

y = r*sin(theta) + yc;

plot(x,y,'m','LineWidth',2)

axis equal

figure(5)

r = t/2;

xc = t*retsu-(t/2);

yc = t*gyo-(t/2);

theta = linspace(0,2*pi);

x = r*cos(theta) + xc;

y = r*sin(theta) + yc;

plot(x,y,'m','LineWidth',2)

%axis equal

toc %計測終わり