허프 변환 - 원(1/2)

허프 변환을 이용한 원의 검출에 대해서 정리한다.

작업 순서는 다음과 같은데 설명은 기본 이론 요약 후 결과 이미지를 보며 설명한다.

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 값을 결정하는 방법 고민.