사용한 라이브러리

【화상 처리 100개 노크】 독자적인 화상 입출력 클래스를 만든다

Q.4. 오쓰의 2치화

오쓰의 이진화를 구현하라. 오쓰의 이진화는 판별 분석법이라 불리며, 이진화에 있어서의 분리의 역치를 자동 결정하는 수법이다.

분리도라고 하는 지표를 최대화하면 됩니다만, 분리도라고 하는 것이 지금 핀 때 없습니다. 분리도를 극대화하는 것은

w_0\sigma_0^2 + w_1\sigma_1^2

최대화하는 것과 같은 것 같습니다. 그러나 $w$는 각 클래스의 비율입니다. 음, 알 수있는 모르는 ....
좀 더 살펴보면 원본 이미지와 이진 이미지의 제곱 오차를 줄이는 임계값을 찾는 것과 동일하다는 설명이 있습니다. 식에 의한 상세한 설명은 발견되지 않았습니다만, 이것이 제일 납득할 수 있을 것 같은 생각이 듭니다.
다음과 같은 이진화를 생각합니다.

g_i =   \left\{
    \begin{array}{ll}
      \mu_0 & i<t\\
      \mu_1 & i\geq t
    \end{array}
  \right.

그러나 $\mu$는 각 클래스의 평균값입니다. 원본 이미지의 픽셀 값을 ${f_i}$로 설정하면 제곱 오차는

\begin{eqnarray}
\sum_i(f_i-g_i)^2 &=& \sum_{\rm class0}(f_i-\mu_0)^2 + \sum_{\rm class1}(f_i-\mu_1)^2 \\
&=& n_0\sigma_0^2+n_1\sigma_1^2\\
&\propto& w_0\sigma_0^2 + w_1\sigma_1^2
\end{eqnarray}

가 오즈의 방법의 식과 같은 식이 나옵니다. 다만 $n$는 각각의 클래스에 속하는 화소수입니다.

코드는 다음과 같습니다.

int main()
{
    PPM ppm("imori.pnm");
    int width = ppm.Get_width();
    int height = ppm.Get_height();
    PPM ppm2(width, height);

    std::vector<int> x;
    for (int j = 0; j < height; j++)
        for (int i = 0; i < width; i++)
        {
            int r = ppm(i, j, 'r');
            int g = ppm(i, j, 'g');
            int b = ppm(i, j, 'b');
            int y = (std::round)(0.2126 * r + 0.7152 * g + 0.0722 * b);
            x.push_back(y);
        }

    double maxeval = 0;
    double T;
    for (int t = 1; t < 255; t++)
    {
        int n[2] = {0};
        double mu[2] = { 0 };
        for (int i = 0; i < width*height; i++)
        {
            if (x[i] < t)
            {
                n[0]++;
                mu[0] += x[i];
            }
            else
            {
                n[1]++;
                mu[1] += x[i];
            }
        }
        mu[0] /= (double)n[0];
        mu[1] /= (double)n[1];
        double w0 = n[0] / (double)(width * height);
        double w1 = n[1] / (double)(width * height);
        double eval = w0 * w1 * (mu[0] - mu[1]) * (mu[0] - mu[1]);
        if (eval > maxeval)
        {
            maxeval = eval;
            T = t;
        }

    }

    for(int j=0; j<height; j++)
        for (int i = 0; i < width; i++)
        {
            int r = ppm(i, j, 'r');
            int g = ppm(i, j, 'g');
            int b = ppm(i, j, 'b');
            int y = (std::round)(0.2126 * r + 0.7152 * g + 0.0722 * b);
            y = (y < T ? 0 : 255);
            ppm2(i, j, 'r') = y;
            ppm2(i, j, 'g') = y;
            ppm2(i, j, 'b') = y;
        }

    ppm2.Flush("out.ppm");
    return 0;
}