画像処理

TOP / Articles

description

画像処理について。 ソースは24bitのビットマップ(及びその互換)、EclipseによるJava環境(org.eclipse.swt.graphicsとjava.utils.Arrays。具体的にはこんな感じ)を前提に記述していますが、ピクセル毎に入出力できる環境であれば任意に読み替えて問題ありません。 具体的には、画像から幅、高さ、指定場所の画素といったものが利用出来、なおかつそれから画像を生成できることが条件となります。 低機能でも良ければC++用の簡易クラス(.cppファイル.hファイル)とかあります。 ただし、簡易クラスの方では色データはRGB纏まった形で取得するのではなくて、GetR等のように一つ一つの色を取り出す必要がありますので、以下のソースコードそのままでは使用できない場合(や、面倒な前処理が不要となる場合)があります。 なお、ソースそれ自体は適当なものなので、変数宣言をループ外に出すなどの少しの心遣いでかなり高速化するかも知れません。 画像処理では数万回のループというのは普通に行われてしまうので、ループの改善は劇的な効果をもたらしたりします。 他にもpure Javaな画像処理だと、Appletへの移植も容易で面白いです。

なお、ソースコードの色付けにはprettifyを利用しています。 対応するブラウザでは多少カラフルに表示されます。

Trident(Internet Explorer等)でのcanvas対応にはExplorerCanvasを利用しています。 canvasは、WebKit(Safari等)、Gecko(Firefox等)、Presto(Opera等)ではネイティブ対応していますが、多少表示が異なることがあります。 確認環境は、Firefox 2.0.0.4、Opera 9.21、Safari 3.0.2ですが、小さな画面だと単位表示が壊れます。 なお、IEだとそもそも正常に単位表示がされない上、変なところにナビゲーションが現れますが、本来は対応さえしていないので、グラフが表示されるだけマトモだと思って諦めて下さい。

基本処理

画像の二値化、モノクロ化、反転などの技術。 これらは基本的には、一画素ごとに抽出して処理しているだけであり、原理としては単純である。 ただし処理の原理自体は単純であるとはいえ、例えば二値化の際の閾値の決定など、工夫を凝らすべきやり甲斐のある点は幾つか見られる。

基本部分のソースは以下の通り:

for (int y = 0; y < image.height; y++) {
	for (int x = 0; x < image.width; x++) {
		int color = image.getPixel(x, y);
		int new_color;  // 処理後の画素を格納
		
		// ここで処理
		
		image2.setPixel(x, y, new_color);
	}
}

また、一部の処理はJava Appletによる実装もしてみたので、アンチJavaアプレットでない方は具体的なイメージを掴むための参考にどうぞ。

反転

画素値を反転させ、例えば白(0xFFFFFF)を黒(0x000000)にするような処理。 具体的には、RGB各要素において、255から画素値を引く演算を行う。

なお、上記演算は各要素のXORを求めるのに等しい。 つまり、例えば赤要素だけを抽出した結果が二進数で00101101であったとすると、これを11010010にする処理に相当する。 これを利用すると、普通に処理するよりも高速な処理が出来る(といっても、レトロな環境か巨大な画像でないと実感できないかもしれない)。

ソースは以下の通り:

int red = 255 - (color & 0xff);
int green = 255 - ((color >> 8) & 0xff);
int blue = 255 - ((color >> 16) & 0xff);

new_color = (blue << 16) + (green << 8) + red;

或いは

new_color = color^0xffffff;  // XOR演算

或いは

new_color = ~color;  // XOR演算

モノクロ化

モノクロ化は、各ピクセルにおいてRGB各要素値を一定にする処理と言える。 具体的には、R121:G121:B121というような状態とすればよい。

モノクロ化に際しては、どの要素をどの程度重要視するかによって完成した画像が異なってくる。 一応、一般に良いとされているのはNTSCのY信号を利用するもので、Rが29.8912%、Gが58.6611%、Bが11.4478%という比率である。 ただし、ソースでは、全ての色を等価に扱う場合を考える。

ソースは以下の通り:

// 合成比率
double rRate = 0.34;
double gRate = 0.33;
double bRate = 0.33;

int red = color & 0xff;
int green = (color >> 8) & 0xff;
int blue = (color >> 16) & 0xff;

// 合成
int elem = (int)Math.round(red*rRate + green*gRate + b*bRate);

elem = (elem > 255) ? 255
     : elem;

new_color = (elem << 16) + (elem << 8) + elem;

モノクロ化に使う重みとしては、NTSCの他にもHDTVのY値を利用する場合がある。

二値化

モノクロ化の応用で、閾値以下の場合に黒、閾値以上の場合に白にする処理である。 閾値をどこに設定するかによって結果は大幅に変化する。

ソースは以下の通り:

// 合成比率
double rRate = 0.34;
double gRate = 0.33;
double bRate = 0.33;

// 閾値
int limen = 127;

int red = color & 0xff;
int green = (color >> 8) & 0xff;
int blue = (color >> 16) & 0xff;

int elem = (int)Math.round(red*rRate + green*gRate + b*bRate);

new_color = (elem > limen) ? 0xffffff
          : 0;

二値化は画像認識などで便利だが、実際に使用する場合にはHSV式を利用するなりして、重要となる情報を元に二値化することになる。

コントラスト変換

コントラスト変換は突き詰めるとかなり面倒な処理となるので、ここでは基礎のみを書く。

例えば、0~255までの値のうち、0~63までを切り捨て、残りを単純強調するコントラスト変換を考える。 このとき、処理式は

となる。

このときのソースコードは以下の通り:

for (int i = 0; i < 3; i++) {
	int cElem = (color >> (i*8)) & 0xff;
	
	if (cElem < 64) {
		cElem = 0;
	} else {
		cElem = (int)Math.round((255/191)*(cElem-64));
		cElem = (cElem > 255) ? 255
		         : cElem;
	}
	
	new_color += (cElem << (i*8));
}

実際に使用する際には、式をもっと複雑にするなり、曲線式にするなりすることになる。

ガンマ補正

コントラスト変換の一種で、ガンマ値を補正するもの。 ガンマ値の補正式はy=255×(x/255)1/γである。

ソースは以下の通り:

// ガンマ値の設定
double gamma = 2.2;

for (int i = 0; i < 3; i++) {
	int cElem = (color >> (i*8)) & 0xff;
	
	cElem = (int)Math.round(255*Math.pow((x/255), 1/gamma));
	cElem = (cElem > 255) ? 255
	         : cElem;
	
	new_color += (cElem << (i*8));
}

見ての通り、x/255として階調を0~1にした後に累乗しているため、ガンマに0或いはマイナス値を入れない限り、問題なく動く。

ちなみに、1/γは(現在のγ)/(目標γ)という意味合いである。 これを応用して、Windowsの一般的なγである2.2からMacintoshの一般的なγである1.8へ変換するということも可能である。

画像拡縮

画像を拡大したり、縮小したりする技術。 主流となっているものは、最近傍法バイリニア補間バイキュービック補完の三種類である。 その他、主として動画に使われることが多い有力な方式として、Lanczos法が存在する他、縮小専用の法式としては面積平均法平均画素法が存在する。 稀にB-SplineBesselMitchellなどが使われることもある。

以下に挙げる補間式にみられる特性として、全ての補間式においてその処理区間において定積分を行い、かつ負の領域となった部分を減法として処理した場合、全て同値の結果となる。 従って、補間式が充たさねばならない特性の一つは積分区間が一定であるということであることが予想される。 なおこの他、Lanczosの項で触れるように、sin(πD)/πDに近似するほど優秀な補間式となる。

canvasが有効であれば、各補間式の特性曲線が表示される。 canvasに対応したブラウザとしてはMozilla Firefoxなどが挙げられる。

画像の拡縮の際には普通逆変換、つまり、拡縮後の画像の座標を基準に元画像の座標を参照し、拡大or縮小後の画像の画素を1ピクセルずつ求めていくという手法を取る。 これは、正変換の場合では元画像から拡大or縮小後の画像の画素を導いたとき、サイズの違いから飛び飛びの画素となってしまうことを防ぐためである。

最近傍法

画像拡縮の技術としては最も原始的なものの一つで、高速だがあまり美しいとはいえない。 方法は、拡大or縮小後の画素を拡縮率で除算し、その値を四捨五入して得られた座標の画素をそのまま利用するだけである。 最も近くにある色を代用しているに過ぎないので、拡大するとギザギザの目立つ汚い画像になるが、処理が速く、また二値化画像などでは二値化状態が保たれる(使用色数が変化しない)等の利点もある。

具体的には以下の通り:

for (int y = 0; y < image2.height; y++) {  // y座標系をスキャン
	for (int x = 0; x < image2.width; x++) {  // x座標系をスキャン
		// 四捨五入して最近傍を求める
		// scaleは拡大率
		int xp = (int)Math.round(x / scale);
		int yp = (int)Math.round(y / scale);
		
		if (xp < image.width && yp < image.height) {
			int color = image.getPixel(xp, yp);
			image2.setPixel(x, y, color);
		}
	}
}

最近傍法では四捨五入した際に存在しないピクセルを指し示してしまうことがある。 例えば、5x5(ピクセルとして取りうるのは0~4)の画像を2倍したとき、9x9地点の画素の元画像の最近傍は5x5地点となり、存在しない地点を示してしまう。 このとき、例えば4x4地点で代用するなど、何らかの処理を施す必要がある。

バイリニア補間法

画像拡大or縮小の技術としてはコストパフォーマンスに優れたものの一つで、4つの画素を合成して画素値を得る。 方法は、拡大or縮小後の画素を拡大率で除算して得た実数値から、その実数値の周辺4画素との距離を算出し重み付けを行う。 要するに、拡大したときに出来た隙間を、その周辺の画素のグラデーションで埋めてしまおうというものである。 最近傍法では少々汚いので、普通はこの方法で拡大を行うが……二値化画像を拡大すると中間色が出来てしまい二値化状態が崩れてしまう。

重み付け関数は、距離をDとすると、各座標毎に1-Dである。

具体的なソースは以下の通り:

for (int y = 0; y < image2.height; y++) {
	for (int x = 0; x < image2.width; x++) {
		int x0 = (int)(x/scale);
		int y0 = (int)(y/scale);
		
		if (x0 < (image.width-1) && y0 < (image.height-1)) {
			// 左上の画素との距離
			double a = x/scale-x0;
			double b = y/scale-y0;
			
			int newcolor = 0;
			int color[] = new int[4];
			
			color[0] = image.getPixel(x0, y0);
			color[1] = image.getPixel(x0+1, y0);
			color[2] = image.getPixel(x0, y0+1);
			color[3] = image.getPixel(x0+1, y0+1);
			
			// Color Element
			int c_element[] = new int[4];
			
			for (int k = 0; k < 3; k++) {
				for(int i = 0; i < 4; i++) {
					c_element[i] = (color[i] >> 8*k) & 0x00ff;
				}
				
				// Xの距離×Yの距離を各画素の重みとする
				int new_element = (int)Math.round((1-a)*(1-b)*c_element[0]
				            + a*(1-b)*c_element[1]
				            + (1-a)*b*c_element[2]
				            + a*b*c_element[3]);
				newcolor += (int)(new_element << 8*k);
			}
			
			image2.setPixel(x, y, newcolor);
		}
	}
}

バイリニア補間では周辺画素を多く利用する関係上、補間に必要な画素を得るのに必要な領域からはみ出してしまうという事態が発生する。 例えば、5x5の画像を2倍にしたとき、9x9地点の画素は4.5の周辺となってしまう。 このような場合には、例えば基点(上の例では左上)を変更したり、その他の補完法を利用したりして何とか凌ぐ必要がある。

バイキュービック補間法

一般的に用いられる画像拡大or縮小の技術としては最も美しいものの一つで、16画素を合成して画素値を得る。 方法は、基本的な考え方はバイリニアと同じだが、距離に応じて関数を使い分けて重み付けする。

使用する関数は、距離をDとすると、各座標毎に

であるが、一般的にはa=-1が用いられている。 aはより小になるにつれてより強いシャープ効果が得られる。

具体的なソースは以下の通り:

// 使い回す変数は最初に処理しておく
double parameter[] = new double[5];
parameter[0] = a + 3.0;
parameter[1] = a + 2.0;
parameter[2] = -a*4.0;
parameter[3] = a*8.0;
parameter[4] = a*5.0;

for (int y = 0; y < image2.height; y++) {
	for (int x = 0; x < image2.width; x++) {
		// 元画像における画素の対応する場所
		double x0 = x/sc;
		double y0 = y/sc;
		
		// 基点とする画素
		int xBase = (int)x0;
		int yBase = (int)y0;
		
		int color = 0;  // 補間する画素
		
		// バイキュービックの処理範囲
		if (xBase >= 1 && xBase < image.width - 2 && yBase >= 1 && yBase < image.height - 2) {
			double color_element[] = new double[3];  // RGB毎に変数を用意しておく
			Arrays.fill(color_element, 0.0);  // 0で初期化
			
			// 基点の周辺16画素を取得して処理
			for (int i = -1; i <= 2; i++) {
				for (int j = -1; j <= 2; j++) {
					// 実際に処理する画素を設定
					int xCurrent = xBase + i;
					int yCurrent = yBase + j;
					
					// 距離決定
					double distX = Math.abs(xCurrent - x0);
					double distY = Math.abs(yCurrent - y0);
					
					// 重み付け
					double weight = 0.0;  // 重み変数
					
					// まずはX座標の距離で重みを設定
					// 1以下、2以下のとき、それぞれ重み付け
					if (distX <= 1.0) {
						weight = 1.0 - parameter[0]*distX*distX + parameter[1]*distX*distX*distX;
					} else if (distX <= 2.0) {
						weight = parameter[2] + parameter[3]*distX - parameter[4]*distX*distX + a*distX*distX*distX;
					} else {
						continue;  // 何も処理しないので、次のループへ
					}
					
					// Y座標の距離で重みを設定
					if (distY <= 1.0) {
						weight *= 1.0 - parameter[0]*distY*distY + parameter[1]*distY*distY*distY;
					} else if (distY <= 2.0) {
						weight *= parameter[2] + parameter[3]*distY - parameter[4]*distY*distY + a*distY*distY*distY;
					} else {
						continue;
					}
					
					// 実際に画素を取得
					int color_process = image.getPixel(xCurrent, yCurrent);
					
					// 画素をRGB分割し、重みをかけて足し合わせる
					for (int k = 0; k < 3; k++) {
						color_element[k] += ((color_process >> k*8) & 0xff)*weight;
					}
				}
			}
			
			// 画素を結合する
			for (int i = 0; i < 3; i++) {
				color_element[i] = (color_element[i] > 255) ? 255
					   : (color_element[i] < 0)   ? 0
					   : color_element[i];
				color += (int)color_element[i] << i*8;
			}
		}
		
		image2.setPixel(x, y, color);
	}
}

バイキュービック補間では、バイリニアと同じく補間に必要な画素が確保できない領域を処理する必要に迫られることがある。 しかも、利用するが素数が多いだけより深刻である。

Lanczos(ランツォシュ)補間法

少々マイナーではあるが美しいと言われる補間法で、特に縮小に関してはLanczos-3が最も美しいと言われたりしている。 ただし計算コストが少々大きいため、速度が要求される場面では用いるのは避けた方がよい。 というか、離散解を生成して用いるなどして、その都度演算するのを避けた方がよい。 必要な画素数はLanczos-NのNによって異なり、(N×2)2の画素を処理する。

使用する関数は、距離をDとすると、各座標毎に(sin(πD)/πD)/(sin(πD/N)/(πD/N))である。 ちなみにsin(πD)/πDはsinc関数と呼ばれるもので、Dを無限遠点にまで応用すれば、実は理論的には最も完全な補間となり、バイキュービックの補間式はこれを近似したものなのであるが、無限遠点まで用いるのは現実的ではないため、普通は用いられない。

なお、注意しなければならない点として、重みを合計した値を最終結果から除算する処理を行う必要がある。 これを怠ると、斑目状の規則的な模様が現れることになる。

具体的なソースは以下の通り:

int n = 2;  // N値
int nx = n-1;

for (int y = 0; y < image2.height; y++) {
	for (int x = 0; x < image2.width; x++) {
		double x0 = x/sc;
		double y0 = y/sc;
		
		int xBase = (int)x0;
		int yBase = (int)y0;
		
		int color = 0;
		
		// ランツォシュの処理範囲
		if (xBase >= nx && xBase < image.width - n && yBase >= nx && yBase < image.height - n) {
			double color_element[] = new double[3];
			Arrays.fill(color_element, 0.0);
			
			double w_total = 0.0;
			
			// 周辺(a*2)^2画素を取得して処理
			for (int i = -nx; i <= n; i++) {
				for (int j = -nx; j <= n; j++) {
					int xCurrent = xBase + i;
					int yCurrent = yBase + j;
					
					// 距離決定
					double distX = Math.abs(xCurrent - x0);
					double distY = Math.abs(yCurrent - y0);
					
					// 重み付け
					double weight = 0.0;
					
					if (distX == 0.0) {
						weight = 1.0;
					} else if (distX < n) {
						double dPIx = Math.PI*distX;
						weight = (Math.sin(dPIx)*Math.sin(dPIx/n))/(dPIx*(dPIx/n));
					} else {
						continue;
					}
					
					if (distY == 0.0) {
						;
					} else if (distY < n) {
						double dPIy = Math.PI*distY;
						weight *= (Math.sin(dPIy)*Math.sin(dPIy/n))/(dPIy*(dPIy/n));
					} else {
						continue;
					}
					
					// 画素取得
					int color_process = image.getPixel(xCurrent, yCurrent);
					
					for (int k = 0; k < 3; k++) {
						color_element[k] += ((color_process >> k*8) & 0xff)*weight;
					}
					
					w_total += weight;
				}
			}
			
			for (int i = 0; i < 3; i++) {
				if (w_total != 0) color_element[i] /= w_total;
				color_element[i] = (color_element[i] > 255) ? 255
					   : (color_element[i] < 0)   ? 0
					   : color_element[i];
				color += (int)color_element[i] << i*8;
			}
		}
		
		image2.setPixel(x, y, color);
	}
}

変形処理

図形形状を変形させる処理。 ここで紹介する処理ではピクセル位置の移動が中心であり、実用上は図形の拡大縮小技術と組み合わせて使う。

アフィン変換

変形処理の中でも最も一般的なものの一つで、平行移動、拡大縮小、回転を一度に行うことが出来る。 それのみならず、反転や歪みまで表現できる。 基礎変形処理の殆ど全てを含んでいるので、画像編集とかやりたい場合にはぜひとも憶えておきたい処理である。

アフィン変換の一般式は以下の式で定義される。

x' = a*x + b*y + c;
y' = d*x + e*y + f;

しかし、この記述ではどのように拡縮回転平行移動を一度に行えるのかわかりにくい。 そこで、各々の処理がどのように行われているかが明確になるように、拡縮回転平行移動に絞ってアフィン変換式を記述してみると以下のようになる。

x' = A*(x*cos(B) - y*sin(B)) + C;
y' = D*(x*sin(B) + y*cos(B)) + E;

上式では、C,Eがそれぞれ平行移動を、A,Dがそれぞれ拡大縮小を、Bが回転量を表している。 この式の意味は明らかだし数式連ねるのも面倒なので詳しい説明は省くが、括弧内は任意の点を原点から任意の角度で回転させる一般式で、それをA,Dで拡縮し、最後にC,Eで平行移動するという流れになっている。 したがって、はじめに挙げた式において、

a = A*cos(B); b = -A*sin(B); c = C;
d = D*sin(B); e = D*cos(B); f = E;

と置いたとき、拡縮回転平行移動が一度に行える。

ところで、この式はx,yから新座標を求める式となっているから、このままでは変形後の画像に穴が空いてしまうので、実用上は逆変換式に変える必要がある。 この辺りは画像の拡大縮小と同じ。 逆変換の式は以下のようになる。

x = (e*x' - b*y' + b*f - c*e)/(a*e - b*d);
y = (-d*x' + a*y' + c*d - a*f)/(a*e - b*d);

見ての通り、数式的にはもっと纏められるのだけれども、画像処理に利用する場合、x'とy'の値が次々に変動してしまうため、纏める意味が無い。 それよりはa~fの値が固定値であることを考えて、可能な限りループ前に計算を終了させておく方が高速化に繋がる。 なお、ここで上に示したa = A*cos(B);等の変数を代入して更に式を纏めていくと、上式の分母が「コスモスコスモス咲いた咲いた」(或いは「コスプレコスプレしたいしたい」)化して、単なるA*Dになってしまうが、これは拡縮回転平行移動に絞って説明したための副作用で、アフィン変換一般に言えるものではないので注意して欲しい。

ソースコードは以下の通り:

int colour;  // 色保存用
int xx, yy;  // 元画像の座標

// 逆変換
for (int x = 0; x < image2.width; ++x) {
	for (int y = 0; y < image2.height; ++y) {
		// 元画像の座標を調べる
		// ちなみに、Math.roundは最近傍法を使っていることを意味するので、
		// ここを変えることで変換後の画像を綺麗にすることが出来る
		xx = (int)Math.round((e*x-b*y+b*f-c*e)/(a*e-b*d));
		yy = (int)Math.round((-d*x+a*y+c*d-a*f)/(a*e-b*d));
		
		// はみ出している場合は黒にする
		if (xx >= 0 && xx < image.width && yy >= 0 && yy < image.height) {
			colour = image.getPixel(xx, yy);
		} else {
			colour = 0;
		}
		
		image2.setPixel(x, y, colour);
	}
}

ただし上のソースコードでは逆変換後の画像から元画像の座標を調べているだけなので、変換後の画像が画面外にはみ出してしまったりする。 これを避けるためには、平行移動パラメータ(c,f)をとりあえず無視して画像の四隅のみをアフィン変換し、画像の描画範囲を調べた上で、はみ出さないように変換画像のサイズと平行移動パラメータを定めてやると良い。 また、多数のパラメータに惑わされがちだが、a*e-b*dの結果が0になるとエラーとなるので、予めチェックする。

フィルタ処理

フォトレタッチソフトなどで使われるフィルタ処理など。 ここでは、いわゆるカスタムフィルタと呼ばれるものや、その応用。 ぼかし処理、シャープネス、エッジ抽出、エンボス、ラプラシアンなどを基本とし、メディアンフィルタやソーベルフィルタなどの応用が存在する。

基本的な処理方法は、基点の周辺の画素を取得し、重み付け等を施して基点画素を設定するというものである。 従って、基礎部分を生成してしまえば係数の変更次第で多数の単純なフィルタを生み出すことが出来る。

基礎部分の処理は以下のような感じになる:

// フィルタの係数を作っておく
// ここではとりあえず原画像のまま無処理の場合を記述する
int weight[][] = {
		{0, 0, 0,},
		{0, 1, 0,},
		{0, 0, 0,},
	};
int offset = 0;
int div = 1;

// 処理範囲に注意
for (int y = 1; y < image.height-1; y++) {
	for (int x = 1; x < image.width-1; x++) {
		int red, green, blue;
		red = green = blue = 0;
		
		// 処理画素を取り出し、処理する
		for (int a = 0; a < 3; a++) {
			for (int b = 0; b < 3; b++) {
				// 処理画素を取り出す
				int tempColor = image.getPixel(x+a-1, y+b-1);
				int tempRed = tempColor & 0xff;
				int tempGreen = (tempColor >> 8) & 0xff;
				int tempBlue = (tempColor >> 16) & 0xff;
				
				// 重み付け
				red += tempRed * weight[a][b];
				green += tempGreen * weight[a][b];
				blue += tempBlue * weight[a][b];
			}
		}
		// 除算する
		red /= div;
		green /= div;
		blue /= div;
		
		// オフセット値を足す
		red += offset;
		green += offset;
		blue += offset;
		
		// 数値チェック
		red = (red > 255) ? 255
		    : (red < 0)   ? 0
			: red;
		green = (green > 255) ? 255
		      : (green < 0)   ? 0
			  : green;
		blue = (blue > 255) ? 255
		     : (blue < 0)   ? 0
			 : blue;
		
		// 画素を生成
		int color = (blue << 16) + (green << 8) + red;
		
		image2.setPixel(x, y, color);
	}
}

このとき、重み付け係数は以下のように表現される(weight変数と同一である)。

0 0 0
0 1 0
0 0 0

一般には、これを関数化して用いると良いだろう。

ぼかし処理

基本的な考え方は、基点画素と周辺画素とを足し合わせ、超過した重みに従って除算するというだけである。

以下に重み付け係数の例を示す:

0 1 0
1 1 1
0 1 0

除数:÷5

ラプラシアンフィルタ

いわゆるエッジ抽出フィルタの一種。 基点画素に大きな重み付けを行い、周辺画素を引くことで、基点が周辺と比較して異常な画素であった場合にだけ強い画素値が残るようにしたものである。

以下に重み付け係数の例を示す:

0 -1 0
-1 4 -1
0 -1 0

シャープネス

ラプラシアンフィルタなどのエッジ抽出フィルタの効果を応用し、エッジを際だたせる。 基点画素に大きな重み付けを行い、周辺画素を引くことで、基点が周辺と比較して異常な画素であった場合にはより強い画素値が残るようにしたものである。

以下に重み付け係数の例を示す:

0 -1 0
-1 5 -1
0 -1 0

エンボスフィルタ

エッジ抽出を応用して、斜め方向のエッジだけを抽出したものに画素の取りうる値の中央値(0~255なら128)を足したもの。 重み付けを行う画素の距離を変更することで、エンボスの深度を変えることが出来る。

以下に重み付け係数の例を示す:

-1 0 0
0 1 0
0 0 0

オフセット:+128

メディアンフィルタ

ノイズ削除に有用なフィルタで、周辺画素をソートした結果の中央値をそのまま画素値とする。 重み付けのような単調な処理ではないため、上記のフィルタのように係数を変更するだけでは実現できない。

言葉だけでは理解しづらいので、例を示す:

例えば、画素が

1 3 7
2 9 4
5 8 6

という配置となっていたとする。 これをソートすると、1, 2, 3, 4, 5, 6, 7, 8, 9となり、5が中央値である。 従って、5を用いる。

具体的なソースは以下の通り:

for (int y = 1; y < (image.height-1); y++) {
	for (int x = 1; x < (image.width-1); x++) {
		int[] csort = new int[9];  // 色値(ソート用)
		int[] ysort = new int[9];  // 輝度値(ソート用)
		
		for (int a = 0; a < 3; a++) {
			for (int b = 0; b < 3; b++) {
				int temp_color = image.getPixel(x+a-1, y+b-1);
				int red = temp_color & 0xff;
				int green = (temp_color >> 8) & 0xff;
				int blue = (temp_color >> 16) & 0xff;
				int Y = (int)(0.299*red + 0.587*green + 0.114*blue);
				int i = a+b*3;
				
				csort[i] = temp_color;
				ysort[i] = Y;
			}
		}
		
		// 中央値を求めるまで輝度値を使ってバブルソート
		for (int i = 0; i < 5; i++) {
			for (int j = i + 1; j < 9; j++) {
				if (ysort[i] < ysort[j]) {
					int tmpy = ysort[i];
					ysort[i] = ysort[j];
					ysort[j] = tmpy;
					
					int tmpc = csort[i];
					csort[i] = csort[j];
					csort[j] = tmpc;
				}
			}
		}
		
		int color = csort[4];
		image2.setPixel(x, y, color);
	}
}

なお、上記の場合、ソートに必要なループは8+7+6+5+4=30回であり、コストがかかりすぎている感がある。 そこで、メディアンフィルタを高速化するアルゴリズムが幾つか提案されている。

代表的な高速化アルゴリズムとしては、浜村氏、入江氏の高速化アルゴリズムなどがある。

ソーベルフィルタ

ソーベルフィルタは輪郭線抽出フィルタの一種。 ラプラシアンと比較して、ノイズを抑え、しかも強いエッジを得ることが出来る。 その代わりカスタムフィルタのような一段階の重み付けで実現できるものではなく、軸毎に重み付けを行う必要がある。

以下に重み付け係数を示す:

X軸:

-1 0 1
-2 0 2
-1 0 1

Y軸:

-1 -2 -1
0 0 0
1 2 1

具体的なソースは以下の通り:

// 重み
int weightX[][] = {
	{-1, 0, 1,},
	{-2, 0, 2,},
	{-1, 0, 1,},
};
int weightY[][] = {
	{-1, -2, -1,},
	{ 0,  0,  0,},
	{ 1,  2,  1,},
};

for (int y = 1; y < (image.height-1); y++) {
	for (int x = 1; x < (image.width-1); x++) {
		int rX = 0;
		int gX = 0;
		int bX = 0;
		
		int rY = 0;
		int gY = 0;
		int bY = 0;
		
		// 処理画素を取り出し、処理する
		for (int a = 0; a < 3; a++) {
			for (int b = 0; b < 3; b++) {
				// 処理画素を取り出す
				int tempColor = image.getPixel(x+a-1, y+b-1);
				int tempRed = tempColor & 0xff;
				int tempGreen = (tempColor >> 8) & 0xff;
				int tempBlue = (tempColor >> 16) & 0xff;
				
				// 重み付け
				rX += tempRed * weightX[a][b];
				gX += tempGreen * weightX[a][b];
				bX += tempBlue * weightX[a][b];
				
				rY += tempRed * weightY[a][b];
				gY += tempGreen * weightY[a][b];
				bY += tempBlue * weightY[a][b];
			}
		}
		
		// 累計値を算出
		int red = Math.abs(rX) + Math.abs(rY);
		int green = Math.abs(gX) + Math.abs(gY);
		int blue = Math.abs(bX) + Math.abs(bY);
		
		// 0~255の範囲に揃える
		red = red > 255 ? 255
			: red < 0   ? 0
			: red;
		
		green = green > 255 ? 255
			: green < 0   ? 0
			: green;
		
		blue = blue > 255 ? 255
			: blue < 0   ? 0
			: blue;
		
		int color = (blue << 16) + (green << 8) + red;
		image2.setPixel(x, y, color);
	}
}

なお、重み付けを斜めにも行うとされている資料も存在した。 斜めの重み付けを行った場合、更に強いエッジが得られるだろう。

斜めの重み付けは以下のようなものである:

0 1 2
-1 0 1
-2 -1 0
2 1 0
1 0 -1
0 -1 -2

その他、最終的な値決定を上記のように絶対値で行わず、自乗を加算して平方根を取ることで行うとされている資料も存在する。 従って、同じソーベルフィルタと記されている場合においても、その出力結果は実装によって大きく異なっていることがありうるだろう。

ロバーツフィルタ

斜め差分の抽出に優れた輪郭線抽出フィルタで、エッヂが強調されず、結果ソーベルフィルタよりも細線が抽出される傾向がある。

以下に重み付け係数を示す:

X軸:

0 0 0
0 1 0
0 0 -1

Y軸:

0 0 0
0 0 1
0 -1 0

プリューウィットフィルタ

ソーベルフィルタの強調しない版といった感じ。

以下に重み付け係数を示す:

X軸:

-1 0 1
-1 0 1
-1 0 1

Y軸:

-1 -1 -1
0 0 0
1 1 1

領域処理

ペイントソフトの塗り潰しや領域選択などにあたる処理。

領域拡張法

リージョングローイング等とも呼ばれる。 領域選択の手法の中でも汎用性が高く、使いやすい手法である。

基本的な考え方は、次の通り:

  1. 何処か一点を注目点に指定する。
  2. 注目点を調べ、領域条件に合格したらマーカを付ける。
  3. 注目点の周辺の点を新たな注目点とし、2-3を繰り返す。
  4. 注目点が無くなったときマーカの付いているところを領域とする。

つまり、次々と隣を見て領域を地道に広げていけば、そのうち領域全体を選択できるんじゃないか、っていう手法である。 特定の条件を充たす領域(矩形とか)の場合にはもっと効率的な選択方法があるのだけれど、領域拡張法では領域がどんな形状をしていようと、領域が繋がってさえいれば使えるという自由度がある。 しかも周辺画素を調べて再帰するだけという単純なアルゴリズム。

で、さっそく実装……といきたいところだが、実は落とし穴があったりする。 領域拡張法を素直に再帰処理で実装すると、領域が大きいときにスタックオーバーフローを起こして死ぬ。 コールスタックはそんなに大きくないんですね。 そこで、Stack(C++ではstd::vector)等を利用して再帰処理を置き換える必要がある。

具体的には以下のようなソースコードとなる:

// 領域の開始点を決める。ここでは(0, 0)としている
int sx = 0;
int sy = 0;

// 領域の開始点の色を領域条件とする
int colour = image.getPixel(sx, sy);

// 開始点をあらかじめStackに収めておく
Pos pos = new Pos(sx, sy);
Stack<Pos> stack = new Stack<Pos>();
stack.push(pos);

// リージョングローイング処理
// Stackが空になるまで繰り返す
while (!stack.empty()) {
	// 注目点を取り出す
	Pos p = stack.pop();
	
	// 注目点にまだマーカが付いていない場合
	if (image2.getPixel(p.x, p.y) != 0x0000ff) {
		
		// マーカを付ける
		image2.setPixel(p.x, p.y, 0x0000ff);
		
		// 右隣を見て、領域条件に合えばStackに収める
		if (p.x+1 < image.width && image.getPixel(p.x+1, p.y) == colour) {
			stack.push(new Pos(p.x+1, p.y));
		}
		
		// 左隣を見て、領域条件に合えばStackに収める
		if (p.x-1 >= 0 && image.getPixel(p.x-1, p.y) == colour) {
			stack.push(new Pos(p.x-1, p.y));
		}
		
		// 下を見て、領域条件に合えばStackに収める
		if (p.y+1 < image.height && image.getPixel(p.x, p.y+1) == colour) {
			stack.push(new Pos(p.x, p.y+1));
		}
		
		// 上を見て、領域条件に合えばStackに収める
		if (p.y-1 >= 0 && image.getPixel(p.x, p.y-1) == colour) {
			stack.push(new Pos(p.x, p.y-1));
		}
	}
}

条件を複雑にしたい場合には「== colour」となっている部分を任意の条件に書き換えると良い。 また、ここではimage2に領域部分のマーカを付けているが、これをマスクとして利用して実際の領域抽出を行うことになる。

また、ソースコード中にあるPosなるクラスは自作の小クラスで、x,yの二つの値を収めている。 こちらのソースは次のような感じ:

class Pos {
	public int x, y;
	public Pos(int x, int y)
	{
		this.x = x; this.y = y;
	}
}

波形解析

画像を波形と捉えて各種の処理を行う手法。 イメージし難いようなら、まず画像を輝度情報のみのグレースケールにして、その輝度情報を(XYZ軸で表現される空間のZ軸の)高さに見立てて画像を三次元グラフにしたと考えるとわかりやすいかもしれない。 こういった技術は主に解析や圧縮に利用される。

波形解析処理には、大きく「フーリエ変換」系の変換と「ウェーブレット変換」系の変換とがある。 フーリエ変換系の変換では、離散フーリエ変換と離散コサイン変換が有名。 ウェーブレット変換系の変換では、ガボールウェーブレット変換とドベシィウェーブレットシリーズの変換が有名。 この他にも一般調和解析なんかがあるけれども、あまり使われていない。

離散フーリエ変換

フーリエ変換の離散バージョン。 主な使い道は周波数解析とコンボリュージョンフィルタの高速演算。 (畳み込み演算の逆処理を応用した)劣化画像の復元処理なんかも行える。 信号処理的な意味合いにおける画像処理の基礎にして最も奥深きものの一つ。

フーリエ変換とは、周期性のある任意の波形を無数のsin波とcos波の合成で表現してしまう技術。 つまり、もとの波形をfとするとき、f=a0*cos(0)+a1*cos(x)+a2*cos(2*x)+...+b1*sin(x)+b2*sin(2*x)+...で表現してしまう技術。 まあ、幾つでもsinとcos使って良いなら、漠然と、どんな波形でも作れそうな気がしないでもないって気分になることと思う。

具体的には、sin(N1x),cos(N2x)(N1=1,2,...,∞, N2=1,2,...∞)とするとき、全てのN1,N2においてsin(N1x)とcos(N2x)が直交し、更にsin同士、cos同士であってもN1(或いはN2)さえ異なれば直交するというような性質があるので、これを利用し、解析波(sin(Nx)或いはcos(Nx))を生成しては内積を取り、また別の解析波を生成しては内積を取り……とすることで解析する。 言葉の説明じゃあ意味不明だけれども。 ここで、周期の等しいcosとsinとの関係は単位円上における実数と虚数との関係に等しいから、オイラーの公式を利用してe^iθの形に纏めると、教科書的なフーリエ変換式となる。

ところで、普通のフーリエ変換式は一次元だけれども、画像は二次元なので、二次元フーリエ変換しなければならない。 そこで二次元フーリエ変換式を探すと、Wikipediaに載っているのが見つかる……のだが、これをそのままプログラムにすると悲惨なことになる。 というのもX^2*Y^2の四重ループ(512x512の画像なら687億ループくらい)になるので、ちょっとカップラーメン作るかとかいう次元を通り越して、軽くメインフレーム時代の実機デバッグを体感出来てしまう。 これじゃあ余りにも遅すぎるので、2の累乗のサイズに限定してFFTなる高速化アルゴリズム(素朴に組んでも512x512の画像で142万ループくらいで、48000万倍くらい高速)が編み出されているのだけれども、アルゴリズム的にちょっと面倒な上ネット上に優秀なソースが幾らでも転がっているので、ぶっちゃけ作る気が起きない(フーリエ変換自体は改造しても楽しいものでもないし)。 なので、ちょっとだけ高速な二次元DFTを作ってみる。

ちょっとだけ高速なDFTの考え方は、二次元DFTを一次元DFTを二回の処理に置き換えるというもの。 つまり、はじめ全てのYにおいてX軸方向に一次元DFTを行う。 その後、全てのXにおいてY軸方向に一次元DFTを行う。 これで二次元DFT処理を行うことが出来る。 こうすることで三重ループ2回の処理X^2*Y+Y^2*X(512x512の画像なら2億6843万ループくらいで、250倍くらい高速だがFFTと比べると190倍くらい時間が掛かる)の計算量に出来るので、FFTと比べると超遅いけれども、普通に二次元で組むよりは相当高速に計算できる。 ……二次元DFTのまま高速化出来そうな気もするんだけれど……、ちょっと思いつかないので、これで行こうと思う。

……などと書き連ねてもウザいだけなのでソースコード:

// このルーチンは順変換(画像→フーリエ)、逆変換両方に使い回せる
// realには実数(順変換では元画像の値)、imageには虚数(はじめは0)が入っている
// invは順変換ではfalse、逆変換ではtrue
// 変換後の値はreal、imageを書き換えて保存される
double M2 = (inv) ? (2*Math.PI) : (-2*Math.PI);
double wnorm = (inv) ? (1.0/wid) : 1.0;
double hnorm = (inv) ? (1.0/hgt) : 1.0;
double rw = 1.0/wid;
double rh = 1.0/hgt;

// X軸方向のみDFT
double tre[][] = new double[wid][hgt];
double tim[][] = new double[wid][hgt];
{
	// テーブル生成
	// Y軸毎に使い回せる部分は予め作っておく
	double wrtable[][] = new double[wid][wid];
	double witable[][] = new double[wid][wid];
	
	for (int i = 0; i < wid; ++i) {
		for (int j = 0; j < wid; ++j) {
			wrtable[i][j] = Math.cos(M2*i*j*rw);
			witable[i][j] = Math.sin(M2*i*j*rw);
		}
	}
	
	// 変換
	// y軸毎に変換する
	for (int y = 0; y < hgt; ++y) {
		
		// 一次元DFT
		for (int x = 0; x < wid; ++x) {
			tre[x][y] = tim[x][y] = 0.0;
			for (int xx = 0; xx < wid; ++xx) {
				tre[x][y] += real[xx][y]*wrtable[x][xx] - image[xx][y]*witable[x][xx];
				tim[x][y] += real[xx][y]*witable[x][xx] + image[xx][y]*wrtable[x][xx];
			}
			// 逆変換の場合には1/widをかける
			tre[x][y] *= wnorm;
			tim[x][y] *= wnorm;
		}
	}
}

// Y軸方向のみDFT
// X軸と同じ要領で行う
{
	// テーブル生成
	double hrtable[][] = new double[hgt][hgt];
	double hitable[][] = new double[hgt][hgt];
	
	for (int i = 0; i < hgt; ++i) {
		for (int j = 0; j < hgt; ++j) {
			hrtable[i][j] = Math.cos(M2*i*j*rh);
			hitable[i][j] = Math.sin(M2*i*j*rh);
		}
	}
	
	// 変換
	for (int x = 0; x < wid; ++x) {
		for (int y = 0; y < hgt; ++y) {
			real[x][y] = image[x][y] = 0.0;
			for (int yy = 0; yy < hgt; ++yy) {
				real[x][y] += tre[x][yy]*hrtable[y][yy] - tim[x][yy]*hitable[y][yy];
				image[x][y] += tre[x][yy]*hitable[y][yy] + tim[x][yy]*hrtable[y][yy];
			}
			real[x][y] *= hnorm;
			image[x][y] *= hnorm;
		}
	}
}

この結果をスペクトル表示にするには、次のようなソースコードを用いる:

double spec[][] = new double[wid][hgt];
double max = 0.0;
double min = Double.MAX_VALUE;

// パワースペクトル(実数の二乗+虚数の二乗を生成し、対数)を生成
// 対数を取っているのは、そのままだと直交成分(よく見かけるフーリエ変換画像の真ん中の明るい部分)が強すぎて
// 殆どの画像は真ん中一点だけ白、その他真っ黒の画像になってしまうので、緩和させる意味合いがある
for (int x = 0; x < wid; ++x) {
	for (int y = 0; y < hgt; ++y) {
		spec[x][y] = Math.log(Math.pow(re[x][y], 2.0) + Math.pow(im[x][y], 2.0))/Math.log(10.0);
		if (spec[x][y] > max) max = spec[x][y];
		if (spec[x][y] < min) min = spec[x][y];
	}
}

int w2 = wid/2;
int h2 = hgt/2;

// 画像データに書き込み
// 直交成分は最初スペクトルの四隅にあるので、一般的な変換画像を見たければスペクトルをずらす必要がある
for (int x = 0; x < wid; ++x) {
	for (int y = 0; y < hgt; ++y) {
		int col = (int)Math.round((spec[x][y] - min)/(max - min)*255);
		if (col > 255) col = 255;
		else if (col < 0) col = 0;
		
		// 直交成分が真ん中に来るようにずらして書き込み
		int x2 = x + w2;
		int y2 = y + h2;
		x2 = (x2 >= wid) ? x2 - wid : x2;
		y2 = (y2 >= hgt) ? y2 - hgt : y2;
		
		image2.setPixel(x2, y2, (col << 16) + (col << 8) + col);
	}
}

なお、ここではモノクロ画像を前提に作ってしまっているが、カラー画像の場合には色毎に処理すれば良いだろう。 或いは各色を取り出さず(RGB全てが渾然となった)intのまま変換するのもありだが、応用する際の利便性が落ちると思われるのでやめておいた方が良いと思う。

ただし、まだ十分遅いので、小手先のテクでもう少しだけ高速化とかしてみようと思う。 そのためには、フーリエ変換の特性の一つで、実数列をフーリエ変換すると(数列長をN-1とし、n≠0のとき、C(n)とC(N-n)とが)複素共役になるというものを利用する。 これを利用すれば、ループ回数が(X/2)*X*Y+Y^2*X(512x512の画像のとき、2億133万ループくらい)となり、何もしない場合と比べて1.3倍くらい高速になる。 しかしそのままでは逆変換できなくなったり、いろいろと問題も多い。 何故か私が試したらスペクトル表示さえ出来なかったけれども、逆変換したら戻ったので成功しているらしい。

ソースコードとしては先に記した変換コードのX軸変換部分を以下のコードに書き換える:

// (共役のための)折り返し地点
int w2 = (int)((wid-1)*0.5 + 0.51);

// X軸方向のみDFT
double tre[][] = new double[wid][hgt];
double tim[][] = new double[wid][hgt];
{
	// テーブル生成
	double wrtable[][] = new double[w2][wid];
	double witable[][] = new double[w2][wid];
	
	for (int i = 0; i < w2; ++i) {
		for (int j = 0; j < wid; ++j) {
			wrtable[i][j] = Math.cos(M2*i*j*rw);
			witable[i][j] = Math.sin(M2*i*j*rw);
		}
	}
	
	// 変換
	for (int y = 0; y < hgt; ++y) {
		// 0のときは共役が無いので外に出す
		tre[0][y] = tim[0][y] = 0.0;
		for (int xx = 0; xx < wid; ++xx) {
			tre[0][y] += real[xx][y]*wrtable[0][xx] - image[xx][y]*witable[0][xx];
			tim[0][y] += real[xx][y]*witable[0][xx] + image[xx][y]*wrtable[0][xx];
		}
		
		// 共役のある部分
		for (int x = 1; x < w2; ++x) {
			tre[x][y] = tim[x][y] = 0.0;
			for (int xx = 0; xx < wid; ++xx) {
				tre[x][y] += real[xx][y]*wrtable[x][xx] - image[xx][y]*witable[x][xx];
				tim[x][y] += real[xx][y]*witable[x][xx] + image[xx][y]*wrtable[x][xx];
			}
			
			// 共役を埋める
			tre[wid - x][y] = tre[x][y];
			tim[wid - x][y] = -tim[x][y];
		}
	}
}

まあでも、実際に利用したい場合にはこんなDFTじゃなくてFFTを何処かから取ってくるのが良い。 本当に速度に天と地との差があるので……。

さて、フーリエ変換を行うことは出来るようになった。 でも、パワースペクトル見るだけじゃ正直つまらない。 これをどう使うか……。

ハイパスフィルタ/ローパスフィルタ

フーリエ変換の結果は周波数帯域毎の強度を示す係数だが、これは概ね、低周波が画像の概要、高周波が画像の詳細(エッジ)に相当する情報を有している。 そこで低周波のみ、或いは高周波のみ取り出すことで、ぼけ画像やエッジ画像を得ることが出来る。 これを、低周波を取り出す場合をローパスフィルタ、高周波を取り出す場合をハイパスフィルタという。

注意点として、周波数情報は、スペクトル表示の時には中心に低周波が来るように移動させているが、本来は四隅に低周波があることを考慮して構築する必要がある。

ソースコードは次のとおり:

// フーリエ変換しておく
// それぞれ、DFTで書いたところのreal, image, invに相当
dft(real, image, false);

// ローパスフィルタの半径
double dist = 15.0;

// ローパスフィルタを収める配列
double lp[][] = new double[wid][hgt];

// 中心座標を求める
int hw = wid/2;
int hh = hgt/2;

// ローパスフィルタ生成
for (int x = 0; x < wid; ++x) {
	for (int y = 0; y < hgt; ++y) {
		// 中心座標からの距離を調べ、半径に収まっていれば1を指定する
		lp[x][y] = (Math.sqrt(Math.pow(x-hw, 2) + Math.pow(y-hh, 2)) <= dist) ? 1 : 0;
	}
}

// 周波数位置にあわせたフィルタ
double lpf[][] = new double[wid][hgt];

// フィルタを周波数位置にあわせる
for (int x = 0; x < wid; ++x) {
	for (int y = 0; y < hgt; ++y) {
		int x2 = x - hw;
		int y2 = y - hh;
		x2 = (x2 < 0) ? x2 + wid : x2;
		y2 = (y2 < 0) ? y2 + hgt : y2;
		
		lpf[x2][y2] = lp[x][y];
	}
}

// フィルタとフーリエ係数をかける
// ここで(1-lpf[x][y])とかければハイパスフィルタになる
for (int x = 0; x < wid; ++x) {
	for (int y = 0; y < hgt; ++y) {
		real[x][y] *= lpf[x][y];
		image[x][y] *= lpf[x][y];
	}
}

// 逆変換する
dft(real, image, true);

ハイブリッドイメージ的なもの

ハイブリッドイメージとはMITで開発された騙し絵作成手法の一つで、異なる画像の低周波と高周波を合成することで、近くで見た場合と遠くで見た場合の印象を大幅に変える技術である。 原典には当たっていないので研究者たちが調査した最適手法は知らないが、おおよその作り方はけっこう簡単で、近くで見えるようにしたい画像にハイパスフィルタをかけ、遠くで見えるようにしたい画像にローパスフィルタをかけ、両者を合成すればよい。 ただしこの手法では若干ノイジーなので、ガウシアンフィルタを代わりに用いるなどしてみる。

特許とかの関係がどうなっているのか知らないので注意する必要があるかも知れない。

ソースコードは次の通り:

// 二枚の画像をフーリエ変換しておく
// 1の方が低周波(下地)、2が高周波(エッジ)になる
dft(real1, image1, false);
dft(real2, image2, false);

double g[][] = new double[wid][hgt];

// 中心座標を調べる
int hw = wid/2;
int hh = hgt/2;

double s = 8.0; // ガウシアンフィルタの係数σ

// ガウシアンフィルタを作る
for (int x = 0; x < wid; ++x) {
	for (int y = 0; y < hgt; ++y) {
		g[x][y] = Math.exp(-(Math.pow(x-hw, 2)+Math.pow(y-hh, 2))/(2*Math.pow(s, 2)));
	}
}

// 周波数位置にあわせたフィルタ
double gf[][] = new double[wid][hgt];

// フィルタを周波数位置にあわせる
for (int x = 0; x < wid; ++x) {
	for (int y = 0; y < hgt; ++y) {
		int x2 = x - hw;
		int y2 = y - hh;
		x2 = (x2 < 0) ? x2 + wid : x2;
		y2 = (y2 < 0) ? y2 + hgt : y2;
		
		gf[x2][y2] = g[x][y];
	}
}

// ハイブリッドイメージ
double real[][] = new double[wid][hgt];
double image[][] = new double[wid][hgt];

// ハイブリッドイメージ生成
for (int x = 0; x < wid; ++x) {
	for (int y = 0; y < hgt; ++y) {
		real[x][y] = real1[x][y]*gf[x][y] + real2[x][y]*(1-gf[x][y]);
		image[x][y] = image1[x][y]*gf[x][y] + image2[x][y]*(1-gf[x][y]);
	}
}

dft(real, image, true);

作ってみるとわかるが、意外にσの設定が難しい。 ここでは固定しているが、実際には利用する画像のサイズや元々のボケ具合(ボケ画像は殆どの成分が低周波なので)によって設定を変える必要があるかも知れない。 それから、利用する画像も慎重に選び、上手く重ね合わさるような画像としなければならず、意外に面倒だったりする。

フーリエ変換を使ったフィルタ処理

まずはフィルタ処理の項に記したようなフィルタ処理を、フーリエ変換を使って実装してみることにする。 先に記したフィルタ処理は、実は「畳み込み(コンボリューション)」というものなのだけれど、これに関係して「畳み込みの定理」なるものがあって、「画像のフィルタリング結果は、画像のフーリエ変換結果とフィルタのフーリエ変換の積に等しい」という。 こう書くと何のことだかイマイチよくわからないが、要するに、今まで一画素処理するのに周辺画素まで使って何度も計算して処理していたのを、たった一度の計算で済ませてしまえるってこと……って、もっとわからないかも。

で、こんなことして何が嬉しいのかというと、計算を高速化出来ること。 まあ、3x3くらいならあまり気にかからないかもだけれど、例えば、ちょっと大袈裟に書くけれども、99x99のぼかしフィルタとかを作りたい場合、普通に処理すると一画素あたり9800回くらい計算しなくちゃならないので、あり得ないくらい時間が掛かってしまう。 そこをフーリエ変換を使って処理すると、一画素につき一度だけの計算になるので、超速くなる。

実装の際に気を付けなければならない点は、フィルタの位置だろう。 フィルタのサイズは画像と同じにする必要があるのだが、そうすると(3x3サイズのフィルタ等だと)サイズが全く違うので何処にフィルタ係数を置けば良いのか迷うと思う。 これは、左上が中心になるようにし、マイナスの座標となった分を逆端にずらし込んで配置するとうまくいく。 この結果、フィルタ係数は四隅に配置されることになる。 何故はみ出した分を逆端に配置するかというと、フーリエ変換は基本的に定常波解析なので、同じ波(画像)が繰り返し繰り返し現れ無限に続くことを想定している。 このため、画像の左端と右端は(フーリエ変換の理論上)繋がっているし、上端と下端もまた繋がっていることになっているので、はみ出した分が逆側に現れることには何の矛盾もない。 ただし、この性質が災いし、端の画素の処理に注意しなければならなくなる。

あと、DFTの項にも書いたけれども、実際に処理するときには何処かからFFT拾ってきて下さい。 DFTだと、常識的なサイズのフィルタだと、むしろ異様に遅くなったりしますので。

で、ソースコードは以下:

// フィルタ(ここではラプラシアンフィルタ)
double conv[][] = {{0,  1, 0},
                   {1, -4, 0},
                   {0,  1, 0}};

// フィルタのフーリエ変換結果を収める配列
double conv_real[][] = new double[wid][hgt];
double conv_image[][] = new double[wid][hgt];

// 初期化
for (int x = 0; x < wid; ++x) {
	for (int y = 0; y < hgt; ++y) {
		conv_real[x][y] = 0.0;
		conv_image[x][y] = 0.0;
	}
}

// フィルタ係数の半分(切り捨て)サイズを調べる
int half_cw = (int)(conv.length*0.5);
int half_ch = (int)(conv[0].length*0.5);

// 一時変数
int nx, ny;

// フィルタを画像と同じサイズに拡張
// フィルタ係数を収める
for (int x = 0; x < conv.length; ++x) {
	for (int y = 0; y < conv[0].length; ++y) {
		nx = wid - half_cw + x;
		ny = hgt - half_ch + y;
		
		if (nx >= wid) nx -= wid;
		if (ny >= hgt) ny -= hgt;
		
		conv_real[nx][ny] = conv[x][y];
	}
}

// 拡張したフィルタをフーリエ変換する
// それぞれ、DFTで書いたところのreal, image, invに相当
dft(conv_real, conv_image, false);

// 一次変数
double r, i;

// 画像とフィルタ両方のフーリエ変換結果をかける
// realとimageには既に画像のフーリエ変換結果が収められているものとする
for (int x = 0; x < wid; ++x) {
	for (int y = 0; y < hgt; ++y) {
		r = real[x][y]*conv_real[x][y] - image[x][y]*conv_image[x][y];
		i = real[x][y]*conv_image[x][y] + image[x][y]*conv_real[x][y];
		
		real[x][y] = r;
		image[x][y] = i;
	}
}

なお、フィルタは掛けるだけではなくて、フィルタ自体のフーリエ変換結果をスペクトル表示してみるのも面白い。 各フィルタがどのような周波数をどの程度通すのか、一目でわかる。

デコンボリューション

前項でフーリエ変換を利用したフィルタ処理を扱ったが、フーリエ変換を利用するとフィルタの畳み込みは単純な掛け算になるのだった。 ということは、フィルタ処理した画像を割ってやれば、元の画像に戻ると考えられる。 まあ実際にはゼロ除算だとかの問題があって、完全に戻るとは限らない。 更に悪いことには、丸め誤差だの何だのでフィルタ演算の結果は不完全な保存をされているため、ゼロ除算でなくても問題が生じる。 しかもフーリエ空間で一部に影響すると、逆変換したときに全体に影響が広がってしまうため、下手をするとよりノイジーな画像になってしまう。

これを応用すると、アレな画像のぼかしを解除するという男の永遠のロマンピンぼけ画像の復元なんかを行うことが出来る。 しかしながら、適当な画像に対して応用しようとすると、劣化関数(フィルタ係数)が未知だったり、アレな画像は普通JPEGなので様々な環境要因によって画像は劣化しているので、ノイズの影響を考慮しなくてはならないなど、そう簡単に男のロマンは叶わない。 まあ、期待に夢膨らませているうちが華なので、どうせ、叶ってしまえば結果に失望するのだろうけれども。

やり方は、基本的には単に割るだけ:

// フィルタ指定(ぼかし)
double conv[][] = {{0,   0.2, 0},
                   {0.2, 0.2, 0.2},
                   {0,   0.2, 0}};

double conv_real[][] = new double[wid][hgt];
double conv_image[][] = new double[wid][hgt];

int half_cw = (int)(conv.length*0.5);
int half_ch = (int)(conv[0].length*0.5);

for (int x = 0; x < conv.length; ++x) {
	for (int y = 0; y < conv[0].length; ++y) {
		int nx = wid - half_cw + x;
		int ny = hgt - half_ch + y;
		if (nx >= wid) nx -= wid;
		if (ny >= hgt) ny -= hgt;
		
		conv_real[nx][ny] = c[x][y];
	}
}

dft(conv_real, conv_image, false);

double r, i;

// ゼロ除算除けと誤差回避のために
// 予め逆数を作っておくことにする
double rdiv;

for (int x = 0; x < wid; ++x) {
	for (int y = 0; y < hgt; ++y) {
		rdiv = Math.pow(combre[x][y], 2.0) + Math.pow(combim[x][y], 2.0);
		
		// 誤差を考慮して少しだけ大きめに範囲指定
		if (rdiv >= 0.001) rdiv = 1.0 / rdiv;
		else rdiv = 0.0;
		
		// 除算
		r = (real[x][y]*conv_real[x][y] + image[x][y]*conv_image[x][y])*rdiv;
		i = (image[x][y]*conv_real[x][y] - real[x][y]*conv_image[x][y])*rdiv;
		
		real[x][y] = r;
		image[x][y] = i;
	}
}

実際にやってみるとわかるが、完全ではないものの結構それらしく復元できる。 ぼけ画像を先鋭化によって擬似的に復元するのとは違う、自然な復元が出来ているのがわかるだろう。 また、先鋭化とぼかしを両方同時に行うような微妙なフィルタに対しても、綺麗に逆変換できるのが確認できる。

離散コサイン変換

離散フーリエ変換の特殊バージョンで、コサインのみを使って実数のみを変換する技術。 フーリエ変換よりも高速で、変換結果がより直交成分に集中しやすいため、圧縮技術として利用されている。 ただし普通は幾つものブロックに分割後ブロック毎に変換するのが一般的で、たとえばJPEGでは8x8のブロックで分割している。

ソースコードはDFTに近いので比べてみると面白い:

// 高速化のため計算しておく
double M2;

// 逆変換時に使用する係数等
double norm1 = (inv) ? 0.5 : 0;
double norm2;
int startx = (inv) ? 1 : 0;

// X軸方向のみDCT
double tre[][] = new double[wid][hgt];
{
	M2 = Math.PI/wid;
	norm2 = (inv) ? 2.0/wid : 1;
	
	// テーブル生成
	double wrtable[][] = new double[wid][wid];
	
	for (int i = 0; i < wid; ++i) {
		for (int j = startx; j < wid; ++j) {
			wrtable[i][j] = (inv) ? Math.cos(M2*(i + 0.5)*j) : Math.cos(M2*(j + 0.5)*i);
		}
	}
	
	// 変換
	for (int y = 0; y < hgt; ++y) {
		for (int x = 0; x < wid; ++x) {
			tre[x][y] = 0.0;
			for (int xx = startx; xx < wid; ++xx) {
				tre[x][y] += real[xx][y]*wrtable[x][xx];
			}
			tre[x][y] += real[0][y]*norm1;
			tre[x][y] *= norm2;
		}
	}
}

// Y軸方向のみDCT
{
	M2 = Math.PI/hgt;
	norm2 = (inv) ? 2.0/hgt : 1;
	
	// テーブル生成
	double hrtable[][] = new double[hgt][hgt];
	
	for (int i = 0; i < hgt; ++i) {
		for (int j = startx; j < hgt; ++j) {
			hrtable[i][j] = (inv) ? Math.cos(M2*(i + 0.5)*j) : Math.cos(M2*(j + 0.5)*i);
		}
	}
	
	// 変換
	for (int x = 0; x < wid; ++x) {
		for (int y = 0; y < hgt; ++y) {
			real[x][y] = 0.0;
			for (int yy = startx; yy < hgt; ++yy) {
				real[x][y] += tre[x][yy]*hrtable[y][yy];
			}
			real[x][y] += tre[x][0]*norm1;
			real[x][y] *= norm2;
		}
	}
}

また、スペクトル表示ではフーリエ変換とは異なり、ズラさないのが普通である。