'circle detection'에 해당되는 글 2건

  1. 2016.05.11 허프 변환 - 원(2/2)
  2. 2016.05.03 허프 변환 - 원(1/2)

허프 변환을 이용한 원의 검출 중에 반지름을 고정시키지 않은 방법에 대해서 정리한다.


기본 이론은 앞에서 정리했고 보완된 작업 위주로 정리한다. 


1) 반지름 범위 반영.

2) 스무딩.

3) 로컬 극대값 계산.

4) 근접 원들 합치기.

5) 결과.

6) 문제점.


1) 반지름 범위 반영

찾고자 하는 원의 크기를 범위로 지정할 수 있는데, 범위가 넓어질 수록 메모리와 계산량이 증가한다. 또한 반지름의 크기를 어느 정도 세밀하게 찾을 것인가를 지정할 수 있는데 이 역시 세밀하게 할 수록 부하가 심해질 것이다. 일단 예제는 반지름이 30~60 픽셀 원을 찾도록 하고 간격은 5 픽셀 정도로 했다. accumulator 공간이 반지름 변수의 추가로 3차원이 되었다.

rmin_idx = 30;

rmax_idx = 60;

r_step = 5;

accumulator = zeros( height+2*rmax_idx, width+2*rmax_idx, rmax_idx );


for r_idx=rmin_idx : r_step : rmax_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, r_idx ) = accumulator( b_idx+r_idx,                                                                                            a_idx+r_idx, r_idx )+1;

                end

            end

        end

    end

end

z축의 관점에서 보면 6군데에서 값이 크게 보인다. 

보이는 각도를 변경시키면 더 명확하게 확인할 수 있다. 반지름이 30, 45, 60정도에 주로 분포한다. 

각 반지름에서의 누적된 분포를 보면 다음과 같다. 







2) 스무딩

로컬 최대값을 찾기 위해서 스무딩처리를 한 후 임계값 이하는 0으로 처리한다. 

thresholdValue = 5.0;

for fig=0 : (rmax_idx-rmin_idx)/r_step

    figure(fig+1);

    

    smoothAccum(:,:,rmin_idx+fig*r_step) = imfilter( accumulator(:,:,rmin_idx+fig*r_step),                                                                       fspecial('disk',5), 'replicate' );

    

    for row=1:yd

        for col=1:xd

            if( smoothAccum(row,col,rmin_idx+fig*r_step) < thresholdValue )

      smoothAccum(row,col,rmin_idx+fig*r_step) = 0.0;

            end

        end

    end


    mesh( smoothAccum(:,:,rmin_idx+fig*r_step) );

    view(-15,70);

    title( strcat( 'radius : ', num2str(rmin_idx+fig*r_step) ), 'color', 'black' );

end


3) 로컬 극대값 계산

로컬 극대값을 찾는 함수를 사용해 극대값 위치를 찾는데 편평한( 모든 값이 일정한 ) 평면의  경우 모든 점을 로컬 극대값으로 계산하기 때문에 편평하지 않은 이미지에 대해서만 계산하도록 한다.

posCnt = 0;

posIdx = zeros(10,3); % ( x, y, radius )

for fig=0 : (rmax_idx-rmin_idx)/r_step

    posImg = imregionalmax( smoothAccum(:,:,rmin_idx+fig*r_step) );

    if( max(max(posImg)) == min(min(posImg)) )

        disp( 'no data' );        

    else

        for y=1 : yd

            for x=1 : xd

                if( posImg(y,x) > 0 )

                    posCnt = posCnt+1;

                    posIdx(posCnt,1) = x-(rmin_idx+fig*r_step);

                    posIdx(posCnt,2) = y-(rmin_idx+fig*r_step);

                    posIdx(posCnt,3) = rmin_idx+fig*r_step;

                end

            end

        end

    end    

end


4) 근접 원들 합치기.

같은 원을 다른 원으로 찾은 것 들을 하나로 합치는데 중심간의 거리가 최소 반지름보다 작은 것들은 하나라고 간주한다.( 중심이나 반지름이 약간 다른 원을 합치면서 정확한 원을 약간 벗어나게 찾을 수 있다. 이 부분은 개선의 여지가 있다. )

[ posCnt col ] = size( posIdx );  

cirIdx = posIdx;

cirIdx(:,4) = zeros();

cirCnt = 0;


for outer=1 : posCnt

if( cirIdx(outer,4) == 0 )

cirCnt = cirCnt+1;

cirIdx(outer,4) = cirCnt;

for inner=outer+1:posCnt

eucDist = sqrt( (cirIdx(outer,1)-cirIdx(inner,1))^2+(cirIdx(outer,2)-cirIdx(inner,2))^2 );

if( eucDist < rmin_idx )

cirIdx(inner,4) = cirIdx(outer,4);

end       

end

end

end


mergedMax = max( cirIdx(:,4) );

mergeIdx = 1;


mergeCnt = zeros( mergedMax, 1 );

xtotal = zeros( mergedMax, 1 );

ytotal = zeros( mergedMax, 1 );

rtotal = zeros( mergedMax, 1 );


for idx=1 : posCnt

xtotal( cirIdx(idx,4) ) = xtotal( cirIdx(idx,4) ) + cirIdx(idx,1);

ytotal( cirIdx(idx,4) ) = ytotal( cirIdx(idx,4) ) + cirIdx(idx,2);

rtotal( cirIdx(idx,4) ) = rtotal( cirIdx(idx,4) ) + cirIdx(idx,3);

mergeCnt( cirIdx(idx,4) ) = mergeCnt( cirIdx(idx,4) )+1;

end

for circle=1 : mergedMax

circleIndex(circle,1) = xtotal(circle)/mergeCnt(circle);

circleIndex(circle,2) = ytotal(circle)/mergeCnt(circle);

circleIndex(circle,3) = rtotal(circle)/mergeCnt(circle);

end



5) 결과



6) 문제점

① 이진화의 정밀도에 따라서 이후 작업에 영향 

    - 노이즈의 비율에 따라 끼치는 영향의 정량적인 분석 필요.

② 원주의 가려지는 비율과 노이즈 또는 threshold 값 결정에 주는 영향 정도 파악.

③ accumulator 공간에서 local maxima를 찾기위한 threshold 값을 결정하는 방법 고민.


'프로그래밍 > 영상처리' 카테고리의 다른 글

ubuntu octave 설치  (0) 2016.07.22
허프 변환 - 원(1/2)  (0) 2016.05.03
허프 변환 - 직선  (0) 2016.04.28
모달리티(Modality)  (2) 2016.04.13
Posted by 게 르 니 카

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


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


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



'프로그래밍 > 영상처리' 카테고리의 다른 글

ubuntu octave 설치  (0) 2016.07.22
허프 변환 - 원(2/2)  (0) 2016.05.11
허프 변환 - 직선  (0) 2016.04.28
모달리티(Modality)  (2) 2016.04.13
Posted by 게 르 니 카
이전버튼 1 이전버튼

블로그 이미지
게 르 니 카

공지사항

Yesterday
Today
Total

달력

 « |  » 2024.5
1 2 3 4
5 6 7 8 9 10 11
12 13 14 15 16 17 18
19 20 21 22 23 24 25
26 27 28 29 30 31

최근에 올라온 글

최근에 달린 댓글

최근에 받은 트랙백

글 보관함