【画像処理100本ノックに挑戦】Q.4. 大津の二値化

使用したライブラリ

【画像処理100本ノック】独自の画像入出力クラスを作る

Q.4. 大津の二値化

大津の二値化を実装せよ。 大津の二値化とは判別分析法と呼ばれ、二値化における分離の閾値を自動決定する手法である。

分離度という指標を最大化すればよいのですが、分離度というものがいまいちピンときません。分離度を最大化することは

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;
}