허프 변환을 이용한 원의 검출에 대해서 정리한다.
작업 순서는 다음과 같은데 설명은 기본 이론 요약 후 결과 이미지를 보며 설명한다.
1) 이론 이해.
2) 이미지 로딩.
3) 이진화.
4) 엣지 검출.
5) HT 공간 변환.
6) accumulator local maxima 찾기.
7) ( x, y )좌표에 원 표시.
8) 문제점.
1) 이론 이해.
일단 원을 수식으로 표현하면 다음과 같다.
원의 중심은 ( a, b ), 반지름은 r 인 원.
( x, y ) 평면을 ( a, b, r ) 공간으로 변환시키기 위해서 원의 식을 다른 방법으로 표현한다.
θ는 0 ≤ θ ≤ 2π .
( x, y ) 평면 상의 후보 픽셀을 ( a, b, r ) 공간으로 변환하는 작업을 아래와 같이 수행한다.
for r from r_min to r_max
for x from x_min to x_max
for y form y_min to y_max
accumulator( a, b, r ) += 1
end
end
end
( a, b, r ) 공간의 값들 중에 지역 최대 값을 찾아 다시 ( x, y )평면으로 변환하면 원을 구할 수 있다.
2) 이미지 로딩
일단 중간 크기 동전을 찾도록 반지름을 45으로 고정해서 계산한다. 이 후에 반지름 고정 없이
모든 원을 찾을 수 있도록 업그레이드 시킨다.
filename = 'coins3.jpg';
image = imread( filename );
3) 이진화 검출
컬러 이미지를 grayscale로 변환하고 Otsu's 방법을 이용해서 threshold 값을 찾고 이진화 후 'open' 모폴로지 연산을 적용하여 작은 조각들을 제거한다.
gray_image = rgb2gray( image );
level = graythresh( uint8(gray_image) );
binary_image = im2bw( uint8(gray_image), level );
binary_image = bwmorph( binary_image, 'open' );
4) 엣지 검출
sobel 필터를 사용해서 엣지 검출.
edge_image = edge( binary_image, 'sobel' );
5) HT 공간 변환.
윗 수식에 맞게 공간 변환 작업 수행.
수식 그대로 적용하면 싸인 함수 때문에 마이너스 인덱스를 갖을 수 있어서 에러가 발생한다.
최대 반지름 만큼 수평 이동시켜 인덱스를 구하고 나중에 값에서 최대 반지름 값을 빼준다.
accumulator = zeros( height+2*r_idx, width+2*r_idx );
for y_idx=1 : height
for x_idx=1 : width
if( edge_image(y_idx,x_idx) > 0 )
for degValue=5 : 5 : 360
a_idx = round( x_idx-r_idx*cos(pi*degValue/180) );
b_idx = round( y_idx-r_idx*sin(pi*degValue/180) );
accumulator( b_idx+r_idx, a_idx+r_idx ) = accumulator( b_idx+r_idx, a_idx+r_idx )+1;
end
end
end
end
6) accumulator local maxima 찾기.
아래 이미지의 가운데 상 하 원의 중심점 찾기.
thresholdValue = 5;
smoothAccum = imfilter( accumulator, fspecial('disk',5), 'replicate' );
[ height2 width2 depth2 ] = size( smoothAccum );
for y=1 : height2
for x=1 : width2
if( smoothAccum(y,x) < thresholdValue )
smoothAccum(y,x) = 0;
end
end
end
posImg = imregionalmax( smoothAccum );
posIdx = zeros(5,2);
posCnt = 0;
for y=1 : height2
for x=1 : width2
if( posImg(y,x) > 0 )
posCnt = posCnt+1;
posIdx(posCnt,1) = x-r_idx;
posIdx(posCnt,2) = y-r_idx;
end
end
end
7) (x,y)좌표에 원 표시.
원 표시.
NOP = 120;
THETA = linspace(0,2*pi,NOP);
RHO = ones(1,NOP)*r_idx;
[ origX, origY ] = pol2cart( THETA, RHO );
for circ=1 : posCnt
X = origX+posIdx(circ,1);
Y = origY+posIdx(circ,2);
plot(X,Y,'g', 'LineWidth', 2);
end
8) 문제점
① 이진화의 정밀도에 따라서 이후 작업에 영향
- 노이즈의 비율에 따라 끼치는 영향의 정량적인 분석 필요.
② 원주의 가려지는 비율과 노이즈 또는 threshold 값 결정에 주는 영향 정도 파악.
③ accumulator 공간에서 local maxima를 찾기위한 threshold 값을 결정하는 방법 고민.
