Otsu 방법을 사용해서 이미지 3개의 클래스로 분할하기 (matlab 소스코드 포함)

▶ Otsu로 이미지 3개의 클래스로 분할하기

 

http://blueskyvision.tistory.com/category/%EC%BB%B4%ED%93%A8%ED%84%B0%20%EB%B9%84%EC%A0%BC%20%26%20%EC%98%81%EC%83%81%EC%B2%98%EB%A6%AC에 Otsu 방법을 사용해서 이미지를 이진화하는 것을 정리했다. 이진화가 이미지를 두 개의 클래스로 분할하는 것이라면, 이번에는 세 개의 클래스로 분할하는 것에 대해 공부해보자. 두 개의 클래스로 이미지를 분할해줄 때 역치값이 하나 필요했다면, 세 개의 클래스로 분할하려면 역치값이 두 개가 필요하다. 만약에 4개의 클래스로 분할하려면 3개의 역치값이 필요할 것이다. 좀 더 나아가서, n개의 클래스로 분할하려면, n-1개의 역치값이 필요할 것이다. n개의 클래스로 가장 잘 분할할 수 있는 n-1개의 intensity값들을 찾는 것이 Otsu 방법의 목적이다. 

 

세 개의 클래스로 분할하는 것을 생각해보자. 간단히 말해서 between-class variance를 구해서 최대가 되게 하는 두 개의 역치값(k_1, k_2)을 구하면 된다. between-class variance는 클래스 간의 분리도를 의미한다. 

 

... (1: between-class 분산)

 

여기서 

 

...(2: 한 픽셀이 클래스 1에 속할 확률)

 

...(3: 한 픽셀이 클래스 2에 속할 확률)

 

...(4: 한 픽셀이 클래스 3에 속할 확률) 

 

...(5: 클래스 1에 속하는 픽셀들의 평균 intensity 값)

 

...(6: 클래스 2에 속하는 픽셀들의 평균 intensity 값)

 

...(7: 클래스 3에 속하는 픽셀들의 평균 intensity 값)

 

이다. 그리고 L-1은 최대 intensity값을 의미한다. 

결국 between-class 분산을 최대가 되게 하는 k_1, k_2를 찾아서, intensity값이 k_1이하인 픽셀들을 클래스 1로, k_1 초과 k_2이하인 픽셀들을 클래스 2로, 나머지는 클래스 3으로 분할해준다. 

 

 

▶ 매트랩으로 구현

 

매트랩에서 제공하는 함수(multithresh)가 있지만, 한번 직접 코드를 짜봤다. multithresh함수를 사용하는 것과 같은 결과를 산출한다.

 

img = imread('flowers.JPG');

imshow(img);

img = rgb2gray(img);

 

[M, N] = size(img);

hist = zeros(256, 1);

 

for i = 1:M

    for j = 1:N

            hist(img(i, j)+1) = hist(img(i, j)+1) + 1; % [0 255] => [1 256]

    end

end

 

hist = hist/(M*N);

 

figure,

plot(hist)

 

%% P1, P2, P3

 

P1 = double(zeros(256, 256));

P2 = double(zeros(256, 256));

P3 = double(zeros(256, 256));

 

for k1 = 2:254

    for k2 = k1+1:255

        P1(k1, k2) = sum(hist(1:k1));

        P2(k1, k2) = sum(hist(k1+1:k2));

        P3(k1, k2) = sum(hist(k2+1:256));

    end

end

 

%% m1, m2, m3, mG

 

m1 = double(zeros(256, 256));

m2 = double(zeros(256, 256));

m3 = double(zeros(256, 256));

intensity = [0:1:255]';

 

for k1 = 2:254

    for k2 = k1+1:255

        m1(k1, k2) = 1/P1(k1, k2)*sum(intensity(1:k1).*hist(1:k1));

        m2(k1, k2) = 1/P2(k1, k2)*sum(intensity(k1+1:k2).*hist(k1+1:k2));

        m3(k1, k2) = 1/P3(k1, k2)*sum(intensity(k2+1:256).*hist(k2+1:256));

    end

end

 

mG = sum(intensity(1:256).*hist(1:256));

 

%% between-class variance

 

squared_sigma_B = P1.*(m1-mG).^2 + P2.*(m2-mG).^2 + P3.*(m3-mG).^2;

temp = squared_sigma_B(:);

 

[sorted_temp, ind] = sort(temp, 'descend');

 

[maximum, index] = max(temp);

 

k1_opt = rem(index, 256);

if k1_opt == 0

    k1_opt = 256;

end

k2_opt = ceil(index/256);

 

k1_opt = k1_opt-1;

k2_opt = k2_opt-1;

 

T = [k1_opt, k2_opt];

seg_I = imquantize(img, T);

RGB = label2rgb(seg_I);

figure, imshow(RGB)

axis off

 

결과 이미지는 그림 1에서 확인할 수 있다. 이진화했을 때보다 좀 더 명확히 꽃을 배경에서부터 구분해냈다. 

 

그림 1. 이미지를 세 클래스로 분할한 결과. 순서대로 원본 이미지, 히스토그램, 분할 결과 이미지