画像のエッジ抽出用フィルタを実装してみる その2 Cannyエッジ検出編

こんにちは。そらです。
画像処理の勉強は変わらず続けています。

前回の記事で紹介したCannyエッジ検出器の実装をしていきたいと思います。空間フィルタリングのときに比べて処理は増えますが、閾値となるパラメータ2つを上手く設定してあげることで空間フィルタリングの時よりもエッジ検出が綺麗にすることができるようになります。




Cannyエッジ検出器の実装のための勉強

処理の流れについて

以下の流れで処理をしていきます。

  1. ノイズ低減
  2. 微分処理(ソーベルフィルタ)
  3. 勾配の最大位置の検出
  4. 細線化処理
  5. 閾値処理(ヒステリシスを使った閾値処理)

1.と2.は過去にやった空間フィルタリングを使用します。簡単に説明をしていきます。
ノイズ処理では、ガウシアンフィルタをかけて画像の平滑化処理を行います。ノイズがあった場合、次の微分処理でノイズ周りにエッジが強く検出されてしまうことが予想されるためです。
次に、ソーベルフィルタで縦方向、横方向でフィルタをかけます。これは、x方向、y方向それぞれに微分をしたことになるので、画像をf(x,y)となる関数で考えると偏微分をしたことになります。
問題は3.~です。

勾配の最大位置の検出と細線化処理について

勾配の最大位置の検出について考えていきます。ここでは、エッジの勾配とエッジの方向(角度情報)を求めます。それぞれ、以下の式を使用します。
続いて、細線化処理を行います。名前の通り線を補足する処理です。具体的には、先程の求めた角度情報を使用しながら、現在注目している画素のまわりのピクセルと比較をして、最大値であれば値を保持し、最大値でなければ0にするという処理になります。
角度情報が0のときは左右と比較、90度の時は上下と比較、45度、135度のときは、斜め右または、斜め左と比較になります。

例 注目画素を中心、角度情報を90度のとき1を考えます。

上下の画素と比較して下の画素のほうが値が大きかったため0にするという処理を行いました。
実際に実装をするときは、全ての画素を0にした状態で処理をするため、最大値のときに値を代入するとして実装をしたほうがよさそうです。

閾値処理(ヒステリシスを使った閾値処理)

閾値処理といえば、デジタル信号で2.0V以上は1,それ以下は0のように閾値が一つあり、それよりも大きいかどうかで判断することが多いですが、ヒステリシスを使った閾値処理では閾値を2つ使用します。値が大きい方をHT、小さい方をLTと書きます。それぞれ、high/low thresholdの略です。
エッジ検出では線を検出したいが、エッジが弱いものは検出をしたくないということが考えられます。この処理では、線を検出すると考えてもらえるとわかりやすいとおもいます。
具体的な処理について説明をします。ある画素に注目したときを考えます。

  1. 画素の値をHTと比較します。HTよりも大きい場合は画素を255で確定して、ここから先の処理は行いません。
  2. 画素の値をLTと比較します。LTよりも小さい場合は画素の値を0として、ここから先の処理は行いません。
  3. 画素の8近傍の画素を確認します。8近傍の画素の中にHTよりも大きい値が存在すれば値を255として確定します。

実装編

前処理の実装

前処理の実装。以前紹介をしたガウシアンフィルタでは画像の端の処理は行いませんでしたが、画像の端よりも外のピクセルは0として実装をし直したため、ガウシアンフィルタなどの前処理の実装を以下に置いておきます。

cv::Mat convert_gray_scale(cv::Mat img)
{
	cv::Mat gray = cv::Mat::zeros(cv::Size(img.cols, img.rows), CV_8UC1);

	for (int y = 0; y < img.rows; y++) {
		for (int x = 0; x < img.cols; x++) {
			// GrayScale = 0.2126R + 0.7152G + 0.0722B
			gray.at<uchar>(y, x) = 0.2126f* (float)img.at<cv::Vec3b>(y, x)[2] + 0.7152f* (float)img.at<cv::Vec3b>(y, x)[1] + 0.0722f* (float)img.at<cv::Vec3b>(y, x)[0];
		}
	}
	return gray;
}

cv::Mat gaussian_filter(cv::Mat img, double sigma, int kernel_size)
{
	cv::Mat kernel = cv::Mat::zeros(cv::Size(kernel_size, kernel_size), CV_32F);

	float gause_sum = 0.0f;

	for (int x = -kernel_size / 2; x <= kernel_size / 2; x++) {
		for (int y = -kernel_size / 2; y <= kernel_size / 2; y++) {
			kernel.at<float>(y + kernel_size / 2, x + kernel_size / 2) = (float)exp(-((x * x + y * y) / (2 * sigma * sigma))) / (2 * PI * sigma * sigma);
			gause_sum += kernel.at<float>(y + kernel_size / 2, x + kernel_size / 2);
		}
	}

	kernel /= gause_sum;

	//std::cout << kernel << std::endl;

	double value = 0.0;

	cv::Mat out = cv::Mat::zeros(cv::Size(img.cols, img.rows), CV_8UC3);

	for (int y = 0; y < img.rows; y ++) {
		for (int x = 0; x < img.cols; x ++) {
			for (int c = 0; c < img.channels(); c++) {
				value = 0.0;
				for (int dx = -kernel_size / 2; dx <= kernel_size / 2; dx++) {
					for (int dy = -kernel_size / 2; dy <= kernel_size / 2; dy++) { if (x + dx >= 0 && y + dy >= 0 && x + dx < img.cols && y + dy < img.rows)
							value += (double)img.at<cv::Vec3b>(y + dy, x + dx)[c] * kernel.at<float>(dy + kernel_size / 2, dx + kernel_size / 2);
					}
				}
				out.at<cv::Vec3b>(y, x)[c] = (uchar)value;
			}
		}
	}

	return out;
}

ケニーフィルタの実装

実装に使用した関数は、エッジ強度を求める関数、角度を求める関数、細線化処理、ヒステリシス処理、ケニーフィルタのすべての処理を実行して出力画像を返す関数です。

#define PI 3.14159265359

cv::Mat get_edge_strength(cv::Mat fx, cv::Mat fy)
{
	int width = fx.cols;
	int height = fx.rows;

	cv::Mat out = cv::Mat::zeros(cv::Size(width, height), CV_8UC1);

	double _fx, _fy;

	for (int y = 0; y < height; y++){
		for (int x = 0; x < width; x++) {
			_fx = (double)fx.at<uchar>(y, x);
			_fy = (double)fy.at<uchar>(y, x);

			out.at<uchar>(y, x) = (uchar)fmin(fmax(sqrt(_fx * _fx + _fy * _fy), 0), 255);
		}
	}

	return out;
}

cv::Mat get_angle(cv::Mat fx, cv::Mat fy)
{
	int width = fx.cols;
	int height = fx.rows;
	cv::Mat out = cv::Mat::zeros(cv::Size(width, height), CV_8UC1);

	double _fx, _fy;
	double angle;

	for (int y = 0; y < height; y++) {
		for (int x = 0; x < width; x++) {
			_fx = fmax((double)fx.at(y, x), 0.000001);
			_fy = (double)fy.at(y, x);

			angle = atan2(_fy, _fx);
			angle = angle / PI * 180.0;

			if (angle < -22.5) { angle = 180.0 + angle; } else if (angle >= 157.5) {
				angle = angle - 180.0;
			}

			if (angle <= 22.5) {
				out.at<uchar>(y, x) = 0;
			}
			else if (angle <= 67.5) {
				out.at<uchar>(y, x) = 45;
			}
			else if (angle <= 112.5) {
				out.at<uchar>(y, x) = 90;
			}
			else {
				out.at<uchar>(y, x) = 135;
			}
		}
	}

	return out;
}

cv::Mat non_maximum_suppression(cv::Mat angle, cv::Mat edge)
{
	int width = edge.cols;
	int height = edge.rows;

	int dx1, dx2, dy1, dy2;
	int now_angle;

	cv::Mat out = cv::Mat::zeros(cv::Size(width, height), CV_8UC1);

	for (int y = 0; y < height; y++) {
		for (int x = 0; x < width; x++) {
			now_angle = angle.at<uchar>(y, x);

			if (now_angle == 0) {
				dx1 = -1;
				dy1 = 0;
				dx2 = 1;
				dy2 = 0;
			}
			else if (now_angle == 45) {
				dx1 = -1;
				dy1 = 1;
				dx2 = 1;
				dy2 = -1;
			}
			else if (now_angle == 90) {
				dx1 = 0;
				dy1 = -1;
				dx2 = 0;
				dy2 = 1;
			}
			else {
				dx1 = -1;
				dy1 = -1;
				dx2 = 1;
				dy2 = 1;
			}

			if (x == 0) {
				dx1 = fmax(dx1, 0);
				dx2 = fmax(dx2, 0);
			}
			if (x == (width - 1)) {
				dx1 = fmin(dx1, 0);
				dx2 = fmin(dx2, 0);
			}
			if (y == 0) {
				dy1 = fmax(dy1, 0);
				dy2 = fmax(dy2, 0);
			}
			if (y == (height - 1)) {
				dy1 = fmin(dy1, 0);
				dy2 = fmin(dy2, 0);
			}

			if (fmax(fmax(edge.at<uchar>(y, x), edge.at<uchar>(y + dy1, x + dx1)), edge.at<uchar>(y + dy2, x + dx2)) == edge.at<uchar>(y, x)) {
				out.at<uchar>(y, x) = edge.at<uchar>(y, x);
			}
		}
	}

	return out;
}

cv::Mat histerisis(cv::Mat edge, int HT, int LT)
{
	int width = edge.cols;
	int height = edge.rows;

	cv::Mat out = cv::Mat::zeros(cv::Size(width, height), CV_8UC1);

	int now_pixel;

	for (int y = 0; y < height; y++) {
		for (int x = 0; x < width; x++) {
			now_pixel = edge.at<uchar>(y, x);

			if (now_pixel > HT) {
				out.at<uchar>(y, x) = 255;
			}
			else if(now_pixel > LT) {
				for (int dy = -1; dy <= 1; dy++) {
					for (int dx = -1; dx <= 1; dx++) {
						if (x + dx < 0 || x + dx > width - 1 || y + dx < 0 || y + dy > height - 1)
							continue;
						else {
							if (edge.at<uchar>(x + dx, y + dy) > HT)
								out.at<uchar>(y, x) = 255;
						}	
					}
				}
			}
		}
	}

	return out;
}

cv::Mat canny_defect(cv::Mat img, int HT, int LT)
{
	// gaussian fileter
	cv::Mat gaussian = gaussian_filter(img, 1.2, 5);
	// gray scale
	cv::Mat gray = convert_gray_scale(gaussian);
	// sobel filter (vertical)
	cv::Mat fy = sobel_filter(gray, 3, false);
	// sobel filter (horizontal)
	cv::Mat fx = sobel_filter(gray, 3, true);
	// get edge strength
	cv::Mat edge = get_edge_strength(fx, fy);
	// get angle
	cv::Mat angle = get_angle(fx, fy);
	// non_maximum_suppression
	cv::Mat nms_edge = non_maximum_suppression(angle, edge);
	// th check
	cv::Mat out = histerisis(nms_edge, HT, LT);
	
	return out;
}

使用するときは、canny_defectにカラー画像とHTとLTの設定をしてあげればOKです。

実行結果

ちゃんとエッジ検出ができました。

さいごに

ケニーのエッジ検出を実装してみました。処理は以外と単純ですが今までの画像処理のプログラムと比べるといろいろとやったため、途中でバグを作ってしまうとどこでおかしくなったのかを探すのがとても大変でした。

OpenCVっていう便利なライブラリはありますが、一度自分の手で実装をしてみると実装力と引数に与えるべきパラメータの決め方が以前よりもわかるようになること間違いなしかなと思っています。どう思いますか?

参考文献

ディジタル画像処理[改訂新版]

  1. パラメータは適当なので実際とは違うと思います。