[python openCV] 이미지 처리 - 임계 처리 (2): otsu algorithm thresholding 논문 분석 및 구현

 

0. 들어가기 전

2020/09/02 - [코딩/opencv] - [python openCV] 이미지 처리 - 임계 처리 (1): inRange, threshold

[python openCV] 이미지 처리 - 임계 처리 (1): inRange, threshold

 

이전 포스트의 마지막에서 threshold의 여러 방법들 중에서 otsu 와 triangle에 대해 간단하게 알아보았다. 주어진 threshold를 사용하지 않고 자체적인 alogrithm을 따라 새로운 threshold를 찾는 방법이다. 이번 포스트에선 otsu algorithm을 자세히 다루고자 한다. 

 

목차는 다음과 같이 하여 작성하겠다.

 

• 논문 전체적인 설명
• 논문 내 수식 및 알고리즘
• 코드 구현

---

1. 논문

 

(1) 원본

링크: ieeexplore.ieee.org/document/4310076

A Threshold Selection Method from Gray-Level Histograms - IEEE Journals & Magazine

논문 제목:

N. Otsu. A threshold selection method from gray level histograms. IEEE Trans. systems. Man. and Cybernetics, 9:62–66, 1979.

 

 

(2) 설명

 

<abstract>

 

otsu algorithm의 장점 및 알고리즘에 대해 간단히 서술하고 있다.

 

요약하면 다른 임계처리 방식들과는 다르게 nonparametric 하고 unsupervised 한 algorithm을 사용한다는 것인데, 이미지의 픽셀 값에 대한 히스토그램만을 이용하므로 어떤 변수가 새로 필요하지 않고, 수식에 의해 자동으로 계산이 되므로 그러하다고 볼 수 있다. 경계에 해당하는 값을 기준으로 class0 과 class1으로 나눴을 때 분리되는 정도가 최대로 되게끔 하는 값이 최적화된 값이라고 설명하고 있다.

 

 

<introduction>

 

이러한 이미지 처리 방식이 필요한 이유 및 나오게 된 배경에 대해 설명하고 있다.

 

0: 검정 255: 흰색

이상적인 경우로 pixel intensity histogram이 위와 같은 모양으로 생겨서 background와 object의 구분이 쉬운 경우인데 (딱 가운데 극소값이 되는 지점) 실제로는 그렇지 않아서 여러 문제점들이 생긴다.

 

이를 해결하기 위해 valley sharpening, difference histogram method, gaussian distribuition으로 근사, 등 다른 방법들이 있지만, 계산 과정이 많을 뿐 아니라 error가 나올 확률도 높다고 지적하고 있다. 

 

이런 단점을 줄이기 위해 이미지의 픽셀 분포로만 해결하는 방법을 고안했다고 한다.

 

 

<formula>

 

경계값을 기준으로 나눴을 때 두 클래스들 각각의 분산과 공분산 등의 개념을 통해 풀어나가고 있다. 자세한건 뒤로 미루겠다. 

 

 

<discussion and remarks>

 

A. analysis of further important aspects

 

수식에 대한 시사점은 뒤로 미루고, 이 방법이 중요한 이유는 gray level에서 affine transformation을 취했을 때 이미 구한 최적화 값은 바뀌지 않는다는 것이다. 

 

B. Extension to Multithresholding

 

임계값을 여러개로 설정할 때의 경우에 대해 설명하고 있는데, 기존의 수식에서 구하려는 임계값의 수를 늘린 형태로 변형시킨 형태로 사용할 수 있다고 한다.

 

다만 클래스의 수가 늘어날수록 threshold의 선택 방법에 대한 신뢰도는 더 떨어지게 되는데, 기존의 식이 점점 더 복잡하게 바뀌기도 하고, 클래스 수가 늘어나면 클래스를 2개로 나눌 때로 고려하고 적용한 수식과 거리감이 생기기 때문이라고 설명할 수 있다. (쉽게 말해서 multi thresholding은 더 많은 클래스, 더 많은 변수를 의미하는데 이런 경우에 잘 맞지 않는다.)

 

여기서 한번더 장점에 대해 언급하는데, M - thresholding에 대해 otsu는 M-1, parametric method는 3M-1개가 필요하다고 한다. 그만큼 다른 방법들에 비해 필요한 수식과 과정이 적다는 것을 알 수 있다. 

 

C. Experimental Results

 

여러 사진들에 대해 처리한 결과를 보여준다.

 

이미지 임계 처리에 대한 방식이 세포 관찰같이 한정된 빛과 관측 환경에서 성능을 높이고자해서 나왔다 정도를 알 수 있다. 

 

<conclusion>

abstract와 내용이 비슷하므로 skip

 

 

---

2. 수식 및 알고리즘

 

1) 기본 개념에 대해 설명.

 

gray level = [ k, ... , L ] 로 나타낼 수 있는데 k는 이미지 중에서 가장 어두운 값, L은 가장 밝은 값이 되겠다. 복잡해보이지만 [0, 1, ... , 255] 라고 생각하면 된다. 앞으로 gray level은 [0,1, ... , 255]라 하겠다. 

 

N 은 이미지 픽셀의 전체 수라고 하자. N = 이미지의 가로 길이 x 세로길이 로 생각할 수 있다.

 

이때 gray level의 어떤 값을 i에 해당하는 픽셀 수를 $$n_i$$ 라 한다면 N에 대하여는 다음이 성립한다.

 

$$ N=n_0 + n_1 + ... + n_{255} $$

 

 

여기서!

 

양변을 N으로 나누면

 

$$ 1 = \frac{n_{0}}{N} + \frac{n_{1}}{N} + . . . + \frac{n_{255}}{N} $$

 

로 나타낼 수 있다. 우변 항들의 합이 1이므로 우변의 각 항을 전체 이미지 중에서 어떤 픽셀 값이 나올 확률로 생각할 수 있다. 즉,

 

$$ 1 = p_{0} + p_{1} + ... + p_{255}$$

 

로 나타낼 수 있다. 

 

 

2) 클래스에 대한 식

 

어떤 값 k ( threshold가 될 값) 을 기준으로 클래스를 2개로 나눈다고 하자. C0, C1이라 하고 각각 background class 또는 object class라 할 수 있다. (또는 그 반대)

 

C0 = [0, 1, ... k]

C1 = [k+1, k+2, ... , L]

 

여기서 1)에서 알게된 점을 이용하여 probability of class occurence를 다음과 같이 정의할 수 있다.

 

 

$$w_{0} = Pr(C_{0}) = \sum_{i=0}^k P_{i} = w(k) $$

$$w_{1} = Pr(C_{1}) = \sum_{i=k+1}^{255(=L)} P_{i} = 1-w(k)$$

 

 

class mean value에 대해서도 정의할 수 있다.

 

$$\mu_{0} = \sum_{i=1}^{k} i \cdot Pr(i|C_{0}) =
\sum_{i=1}^{k} i \cdot p_{i} / w_{0} = \mu(k) / w(k)$$

 

$$\mu_{1} = \sum_{i=k+1}^{L} i \cdot Pr(i|C_{1}) =
\sum_{i=k+1}^{L} i \cdot p_{i} / w_{1} = \frac{\mu_{T} - \mu(k)}{1-w(k)}$$

 

여기서 

 

$$w(k) = \sum_{i=1}^{k} p_i  \space \space \space\space ... (1)$$

$$\mu(k) = \sum_{i=1}^{k} i \cdot p_i  \space \space \space\space ... (2)$$

$$\mu_{T} = \mu(L) = \sum_{i=1}^{L} i \cdot p_i  \space \space \space\space... (3)$$

 

(1)는 0차 cumulative moment, (2)는 1차 cumulative moment, (3)는 이미지 전체의 평균이라고 생각하면 된다.

 

여기서 class mean value와 class occurence에 대해 식을 정리하면 다음 두 식이 k 값에 상관없이 무조건 성립한다. 

 

$$ w_0\mu_0 + w_1\mu_1 = \mu_T$$

$$ w_0 + w_1 = 1 $$

 

 

3) class의 분산에 관한 식

 

각 클래스의 분산은 다음과 같이 정리할 수 있다. 

 

$$ \sigma_{0}^2 = \sum_{i=1}^{k}(i-\mu_{0})^2\cdot Pr(i|C_0) = \sum_{i=1}^{k}(i-\mu_{0})^2\cdot p_i / w_0 $$

$$ \sigma_{1}^2 = \sum_{i=k+1}^{L}(i-\mu_{0})^2\cdot Pr(i|C_1) = \sum_{i=k+1}^{L}(i-\mu_{1})^2\cdot p_i / w_1$$

 

여기까지는 통계나 확률에 대해 좀 공부하신 분들이라면 이해하는데 무리가 없을 것이라 생각한다. 지금까지는 discrete random variable에 대한 평균, 분산 구하기 과정이었다. 

 

자 이제 경계값을 기준으로 해서 두 클래스로 나누고 각 클래스의 평균과 분산을 구하는 식까지 모두 살펴보았다. 이제 무엇이 필요한가? 

 

어떤 기준으로 두 클래스가 잘 구분되었다고 판단할 수 있는지 그 식이 있어야 한다. 원래 의도는 background와 object에 해당하는 픽셀들을 구분하는 것이므로.

 

이 논문에서는 다음과 같은 개념을 설명한다. (여기서부터 중요)

 

$$ \sigma_{w}^2 = w_0\sigma_0^2 + w_1\sigma_1^2  \space \space \space\space ... (4) $$ 

$$ \sigma_{B}^2 = w_0(\mu_0 - \mu_T)^2 + w_1(\mu_1 - \mu_T)^2=w_0w_1(\mu_0 - \mu_1)^2 \space \space \space\space ... (5) $$

$$ \sigma_T^2 = \sum_{i=1}^L(i-\mu_T)^2\cdot p_i \space \space \space\space ... (6) $$

 

(4)번이 의미하는 within-class variance를 뜻한다. (5)는 between-class variance를 뜻한다. (6)은 전체의 분산을 뜻한다. 

(6)번 식이야 원래 알던 공식이니까 넘어가고 대관절 (4)랑 (5)란 식은 어디서 나왔느냐...

 

+ 이 부분에 대해서는 논문의 reference 논문을 참고하여 읽어보다가, LDA (linear discriminant analysis)에 대한 개념을 공부하여 부족하나마 이해하게 되었다. 원본과 추가 자료에 대한 정보를 접은 글로 남기겠다.

 

원본(otsu paper의 reference):

K, Fukunage, Introduction to Statisticul Pattern Recogniition. New York:
Academic, 1972, pp. 260-267.

 

위의 페이지 말고도 pdf기준 ch10.2, p464정도부터 discriminant analysis에 대한 개념이 나온다.

 

 

이해하기 도움되는 글 링크:

 

ratsgo.github.io/machine%20learning/2017/03/21/LDA

 

다시 본론으로 와서, (4)와 (5)의 개념은 class seperability (클래스끼리 구분되는 정도) 을 측정하기 위한 값이다. 접은 글에 있는 책 464쪽을 참고해보면 다음과 같은 설명이 나온다.

 

지금부터 나올 수식에서 sigma의 i=1부터 L까지는 클래스의 수를 의미한다. 위의 gray level과 혼동주의.

 

within-class scatter matrix: 각 클래스의 평균을 기준으로 샘플들이 얼마나 흩어져 있는지 정도를 나타낸다. 

 

$$ S_w = \sum_{i=1}^L P_i E[(X-M_i)(X-M_i)^T]|w_i| $$

 

 

between_class scatter matrix: mixture mean에 대하여 샘플들이 얼마나 흩어져 있는지 정도를 나타낸다. 

 

$$ S_b = \sum_{i=1}^L P_i(M_i-M_0)(M_i-M_0)^T $$

 

여기서 M_0가 mixture distribuition의 expected vector를 나타낸다는데 수식으론 다음처럼 나타난다.

 

$$ M_0 = E(X) = \sum_{i=1}^L P_i M_i $$

 

왜 갑자기 matrix인가... 싶을 수 있겠지만 gray level 0부터 255까지의 값과 두 개의 클래스만 따지기 때문에 matrix도 아니고 단순 값이 된다. 그래서 matrix라 생각하지 말고 단순 값을 생각하고 식을 정리하면 (4)와 (5)의 식이 나옴을 알 수가 있다. 식들의 의미와 나오게 된 배경은 이러하다. 

 

(5)번 식이 왜갑자기 저렇게 바뀌냐에 대한 풀이는 이렇다. 직접 식을 풀어보고 싶은 사람은 참고

 

$$\sigma_B^2 = w_0(\mu_0-\mu_T)^2 + w_1(\mu_1 - \mu_T)^2 = w_0(\mu_0^2 - 2\mu_0\mu_T + \mu_T^2) + w_1(\mu_1^2 - 2\mu_1\mu_T + \mu_T^2) = \mu_T^2(w_0 +w_1) + w_0\mu_0^2 + w_1\mu_1^2 - 2\mu_T(w_0\mu_0 + w_1\mu_1)$$

 

여기서 3) 바로 전에 언급한 식들을 이용하면

 

$$ w_0\mu_0^2 + w_1\mu_1^2 - \mu_T^2(=(w_0\mu_0 + w_1\mu_1)^2) = \mu_2^2w_0(1-w_0) + \mu_1^2w_1(1-w_1) - 2\mu_0\mu_1w_0w_1 = w_0w_1(\mu_0^2 + \mu_1^2 -2\mu_0\mu_1) = w_0w_1(\mu_0-\mu_1)^2 $$

 

이란 식이 나오게 된다. 

 

이 (4)~(6) 식들을 이용하여 다음의 값들을 설정한다.

 

$$ \lambda = \sigma_B^2/\sigma_w^2 $$

$$ K = \sigma_T^2/\sigma_w^2 $$

$$ \eta = \sigma_B^2 / \sigma_T^2 $$ 

 

이 값들이 최대가 되도록 하는 threshold 값이 목표다. (두 번째 식 K는 경계값 k와 다르다.)

 

여기서 논문에서는 eta를 기준으로 잡는 것이 가장 계산하기 쉬우므로 이를 기준으로 하는데, 분모는 고정된 값이므로 분자즉 between-class variance가 최대가 되는 값을 찾고자 한다. 분자에 대한 식은 최종적으로 다음과 같이 정리된다.

 

$$ \sigma_B^2 = \frac{[\mu_Tw(k) - \mu(k)]^2}{w(k)[1-w(k)]} \space \space \space\space ... (*) $$ 

 

필요한 변수는 w(k)와 u(k). 즉 경계값 k로 나눴을 때 첫번째 클래의 발생 확률과 평균이다. 지금까지 거쳐온 과정은 복잡하지만 실제로 구해야할 값은 2개로 줄어든다.

 

 

 

---

3. 코드 구현

(1) 계산 간소화

 

이 과정은 사실 필요없을 수 있는데, '코드에 맞게 (*)에 해당하는 식을 간단하게 나타낼 수 없을까?' 라는 궁금증으로 시도해보았다. 어찌됐든 여기서 픽셀에 대한 확률이란 해당 픽셀값의 개수와 전체 이미지 픽셀 수의 비율로 구분하는 것이니까.

 

따라서 (*) 식에서

 

$$ \mu(k) = \sum_{i=1}^Ki\cdot P_i = \frac{\sum_{i=1}^K i \cdot n_i}{N} $$ 

$$ \mu(T) = \sum_{i=1}^Li\cdot P_i = \frac{\sum_{i=1}^L i \cdot n_i}{N} $$

$$ w(k) = \frac{\sum_{i=1}^K n_i}{N} $$ 

 

위에서 언급했지만 n_i는 픽셀값을 0부터 255라고 했을 때 어떤 값 i의 빈도수에 대한 값이다. 히스토그램을 함수라 생각하면 되겠다. N은 전체 픽셀 수.

 

으로 정하고 식을 대입해보면

 

$$ \frac{[(\sum_{i=1}^Li n_i) \cdot (\sum_{i=1}^K n_i) - N( \sum_{i=1}^Ki n_i)]^2}{\sum_{i=1}^K n_i\cdot(T-\sum_{i=1}^K n_i) \cdot N^2} $$ 

 

위와 같이 나타낼 수 있다. 이제 histogram을 iteration하면 되므로 계산이 한층 간단해진다. 

 

(2) 코드

 

import cv2
import time
import numpy as np
import matplotlib.pyplot as plt


road = cv2.imread('./road.jpg')
gray = cv2.cvtColor(road, cv2.COLOR_BGR2GRAY)
height, width = gray.shape


'''
make histogram of pixel intensity
variable 'hist' is presented as list of which shape is (256, 1).
each element of hist shows how many pixel there are.
ex) hist[70] = [618.]
'''
hist = cv2.calcHist([gray], [0], None, [256], [0,256])


class Otsu:
    def __init__(self, histogram):
        """
        hist : image의 pixel histogram
        where: threshold 값.
        stack: otsu algorithm을 쓸 때 계산 과정을 줄이기 위함인데 없어도 됨.
        all: 전체 픽셀 수. N_T에 해당
        bin: 픽셀 범위. 0부터 255
        total: 픽셀 값 * 픽셀 값의 빈도수를 전부 합한 것.
        k: total을 all로 나눈 값으로 image의 평균 픽셀 값을 의미한다.
        """
        self.hist = histogram.flatten()
        self.where = 0
        self.stack = 0
        self.all = height * width
        self.bin = np.arange(0, 256)
        self.total = int(np.dot(self.bin, self.hist))
        self.k = int(np.rint(self.total / self.all))

    def show_initial(self):
        plt.plot(self.hist)
        plt.scatter(self.k, self.hist[self.k], s=50, c='g')
        plt.annotate('initial k: '+str(self.k), xy=(self.k, self.hist[self.k]), xytext=(self.k, 0),
                     arrowprops={'color': 'green'})
        plt.show()

    def parameter_set(self, k):
        """
        set first class (0~K) parameters
        """
        num0 = np.sum(hist[:k+1])
        moment0 = np.dot(self.bin[:k+1], self.hist[:k+1])

        return num0, moment0

    def process(self):
        w_ini, u_ini = self.parameter_set(self.k)

        w_l = w_ini
        w_r = w_ini
        u_l = u_ini
        u_r = u_ini
        i = 0
        btw_r, btw_l = 0, 0
        large = 0
        while self.k - i > 0 and self.k + i < 255 or self.stack < 10:

            if i == 0:
                pass
            else:
                w_l += -self.hist[self.k-i]
                w_r += self.hist[self.k+i]
                u_l += -self.hist[self.k-i]*(self.k-i)
                u_r += self.hist[self.k+i]*(self.k+i)

            bottom_l = w_l * (self.all-w_l)
            upper_l = np.power((self.total*w_l - self.all * u_l), 2)
            btw_l = upper_l / bottom_l
            bottom_r = w_r * (self.all-w_r)
            upper_r = np.power((self.total * w_r - self.all * u_r), 2)
            btw_r = upper_r / bottom_r

            larger = max(btw_l, btw_r)
            if large < larger:
                self.stack = 0
                large = larger
                if btw_l == larger:
                    self.where = self.k - i
                elif btw_r == larger:
                    self.where = self.k + i
                else:
                    pass
            else:
                self.stack += 1
                pass
            i += 1

    def result(self):
        self.process()
        print(self.where)
        return self.where

    def show_result(self):
        plt.plot(self.hist)
        plt.scatter(self.where, self.hist[self.where], s=50, c='r')
        plt.annotate('optimal k: '+str(self.where), xy=(self.where, self.hist[self.where]), xytext=(self.where, 0),
                     arrowprops={'color': 'red'})
        plt.show()


start = time.perf_counter()
ot = Otsu(histogram=hist)
ot.result()
finish = time.perf_counter()
ot.show_result()
print('elapsed: ', finish - start)

## comparison
print('compare with opencv code')
start = time.perf_counter()
ret, binary = cv2.threshold(gray, 128, 255, cv2.THRESH_OTSU)
finish = time.perf_counter()
print(ret)
print('cv2 code: ', finish - start)

 

결과는 다음과 같다.