허프 변환 - 원(2/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 값을 결정하는 방법 고민.