Otsu 방법을 사용해서 이미지 이진화하기 (matlab 소스코드 포함)

▶ OTSU란?

Otsu 방법은 1979년에 일본의 Otsu라는 사람이 개발한 이미지 분할 방법이다. 우선 이 알고리즘의 원리를 간단히 정리하자면, 이미지의 히스토그램을 이용해서 이미지를 두 개의 클래스로 가장 잘 분할할 수 있는 intensity값이 얼마인지를 찾는 것이다. 그 값을 Threshold 값으로해서 클래스 1과 클래스 2로 분류한다. 그렇다면 어떻게 최적의 Threshold 값을 찾을 수 있을까? 이것이 Otsu 방법의 핵심이다. 

그럼 차근차근 한 단계 한 단계 살펴보자. 먼저 이미지의 히스토그램을 구한다. 

{0, 1, 2, ..., L-1}과 같이 L 개의 intensity 레벨을 갖고 있는 M x N 크기의 이미지를 생각해보자. intensity값 i를 갖고 있는 픽셀의 갯수를 n_i라고 하자 (아래 첨자를 쓰면 너무 번거로워서 n_i로 표시). 그러면 이미지의 픽셀의 갯수, 즉 MN은 n_0 + n_1 + n_2 + ... + n_{L-1}와 같을 것이다. 그리고 전체 픽셀의 갯수로 각 n_i를 나눠주면 normalized 히스토그램을 얻을 수 있다. 

p_i = n_i/(MN) ... (1: intensity값 i를 가지는 픽셀이 이미지에 포함되어 있는 비율)

자연스럽게 p_1 + p_2 + p_3 + ... + p_{L-1} = 1이 된다. 여기서 p_i는 intensity값 i를 가지는 픽셀이 이미지에 포함되어 있는 비율로 볼 수 있다. 만약 p_i = 0.1이라면 전체 이미지의 십분의 일에 해당하는 픽셀들의 intensity 값이 i라는 것이다.  

좀 전에 말했듯이 이미지를 두 개의 클래스(C_1, C_2)로 분할하기 위해서는 적절한 역치값 T(k) = k (여기서, 0 < k < L-1을 만족한다)을 찾아야 한다. 어떤 값이 적절한지 아직은 알 수 없지만, 일단 찾게 되면 [0, k] 내의 intensity 값들을 가진 픽셀들은 클래스 1(C_1)으로 분류될 것이고, [k+1, L-1] 내의 intensity 값들은 가진 픽셀들은 클래스 2(C_2)로 분류될 것이다. 이제 아래로 지루한 식들이 전개되는데, 결국 목적은 between-class variance를 구하는 공식에 필요한 요소들이기 때문임을 기억하고 차근차근 따라가면 된다. 그렇게 어려운 전개는 없다. 일단 between-class variance는 두 개의 클래스 사이의 분리도라고 생각하면 좀 더 이해에 좋을 것 같다. 이 분리도를 크게 만드는 역치값을 찾는 것이 우리의 궁극적 목표다. 

한 픽셀이 클래스 1에 속하게 될 확률은 

... (2: 한 픽셀이 클래스 1에 속하게 될 확률, 다른 말로 하면 클래스 1이 발생할 확률)

이다. 동시에 한 픽셀이 클래스 2에 속하게 될 확률 P_2(k) 는 1 - P_1(k)가 된다. 모든 픽셀은 클래스 1 아니면 클래스 2에 속하므로 두 확률의 합은 1이다.다음으로 클래스 1에 속하는 픽셀들의 평균 intensity 값은 아래와 같이 나타낼 수 있다.   

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

여기서 P(i|C_1)은 intensity값 i가 클래스 1에서 온 것이라는 조건 하에 intensity값이 i일 확률이다. 첫번째 행에서 두번째 행으로의 전개는 베이즈 공식에 의해 이뤄진다. 두번째 행에서 P(C_1|i)는 intensity 값으로 i가 주어졌을 때 클래스 1인 확률을 의미하는데, 우리는 클래스 1에서 나온 i의 값만 다루고 있기 때문에 1이다. P(i)는 픽셀의 intensity값이 i일 확률이므로 p_i와 동일하다. 그리고 P(C_1)은 클래스 1이 나올 확률이므로, P_1(k)와 같다. 그래서 두번째 행에서 세번째 행으로 전개가 가능하다. 마찬가지로 한 픽셀이 클래스 2에 속하는 픽셀들의 평균 intensity값도 아래와 같이 표현할 수 있다.

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

intensity 레벨 k까지의 평균 intensity는 

... (5: intensity 레벨 k까지의 평균 intensity) 

이 되고, 전체 이미지의 평균 intensity는 

... (6: 전체 이미지의 평균 intensity)

이다. 지금까지 P_1, P_2, m_1, m_2를 구했다. 다시 정리하면 P_1은 한 픽셀이 클래스 1에 속하게 될 확률이고, P_2는 클래스 2에 속하게 될 확률이다. 그리고 m_1, m_2는 클래스 1, 2에 각각 속하는 픽셀들의 평균 intensity 값이다. 확률변수의 기댓값 공식에 따라 (P_1)(m_1) + (P_2)(m_2) = m_G이 성립된다. 이번에는 전체 이미지의 분산 intensity값을 구해보자.

... (7: 전체 이미지의 분산 intensity)

그러면 여기서 잠시 도대체 어떻게 최적의 역치값을 찾는 것인가? 이를 위해서 Otsu는 between-class variance라는 개념을 도입한다. 

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

between-class 분산은 각 이 식의 두번째 행을 보면 m_1과 m_2의 차가 클수록 between-class 분산이 커짐을 알 수 있다(클래스 1에 속하는 픽셀들의 평균 intensity 값과 클래스 2에 속하는 픽셀들의 평균 intensity 값의 차가 클수록 between-class 분산이 커진다). 즉, between-class 분산이 커질수록 클래스간의 분리도가 커지는 것이다. 따라서 이 between-class 분산을 최대화하는 k를 찾는 것이 바로 우리의 목적이다! 

공식 (8)을 k를 집어넣어서 다시 쓰면, 

... (9) between-class 분산

이 된다. 최적의 k(k*로 표기)를 찾는 것은 간단하다. k에 intensity 범위인 0부터 L-1의 값을 하나씩 공식(9)에 넣어서 최대가 되게 하는 k가 무엇인지 찾으면 된다. 이렇게 구한 k*을 기준으로 이미지를 분할하는 것이 바로 Otsu 방법의 원리다. 

▶ OTSU 코드 구현

연습삼아 직접 매트랩 코드를 짜봤다. 매트랩에서 제공하는 otsu함수(otsuthresh)와 결과를 비교했을 때 동일한 결과가 산출되는 것으로 봐서 코딩상의 오류는 없는 것으로 보인다. 

clc, clear, close all

img = imread('flowers.JPG');

img = rgb2gray(img);

imshow(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 = zeros(256, 1);

for k = 1:256

    P1(k) = sum(hist(1:k));    

end

P2 = 1 - P1;

m = zeros(256, 1);

intensity = [0:1:255]';

for k = 1:256

    m(k) = sum(intensity(1:k).*hist(1:k));

end

m_G = m(256);

squared_sigma_B = ((m_G*P1 - m).^2)./(P1.*P2);

[maximum, k_opt] = max(squared_sigma_B);

BW = imbinarize(img, (k_opt-1)/255);

figure, imshow(BW)

구현된 결과는 그림 1에 있다. 꽃과 배경을 분리하고 싶었는데, 부분적으로 성공했다. 세 개의 클래스로 분할하면 꽃만 따로 분할하는 것이 가능하다 (다음 포스팅 참고). 

 

그림 1. otsu 이진화 결과. 순서대로 원본 이미지, 히스토그램, 분할된 결과 이미지

<참고 자료>

[1] 곤잘레스, Digital Image Processing, 3판, p. 764-769.