画像補間で習う高速化入門

プログラム C++

タイトルが毎回大げさになるのはそれだけ内容に自信がな(ry

画像(bitmap)は、格子点のみ定義されている。なので中途半端な座標の場合は何かしらの補間が必要である。
補間の種類はいろいろあるのだがとりあえず有名な３つを取り上げる。

  • NearestNeighbor
  • Linear
  • Cubic

簡単に言えば上から順に０次近似、１次近似、３次近似である。あくまで代表的なだけでこれ以外にもいくらでもあるし場合によっては自作することもあるが詳しくは例によってwikipediaにでも聞いてください。
これらの補間にはアフィン変換と合わせて使うことが多いのだが今回は補間を単体で扱ってみる。つまり、ある画像に関して中途半端な座標の値を得るクラスを作ってみる。

template< class I >
class Interpolation2D
{
public:
	enum Type{
		NN,LINEAR,CUBIC
	};

	Interpolation2D( I*** aaaImage, int iSizeX, int iSizeY, int iSizeZ, Type type, I out = 0 ):
		m_aaaImage(aaaImage),m_iSizeX(iSizeX),m_iSizeY(iSizeY),m_iSizeZ(iSizeZ),m_type(type),m_Out(out){}

	I operator()( double dX, double dY, double dZ )const{
		switch( m_type ){
			case NN:
				return GetNN( dX, dY, dZ );
			case LINEAR:
				return GetLinear( dX, dY, dZ );
			case CUBIC:
				return GetLinear( dX, dY, dZ );
			default:
				throw;
		}
	};
protected:
	I GetNN( double dX, double dY, double dZ )const{
		const int iX = static_cast< int >( dX + 0.5 );
		const int iY = static_cast< int >( dY + 0.5 );
		const int iZ = static_cast< int >( dZ + 0.5 );
		if( 0 <= iX && iX < m_iSizeX &&
			0 <= iY && iY < m_iSizeY &&
			0 <= iZ && iZ < m_iSizeZ ){
			return m_aaaImage[iZ][iY][iX];
		}else{
			return m_Out;
		}
	}
	I GetLinear( double dX, double dY, double dZ )const{
		const int iX = static_cast< int >( dX );
		const int iY = static_cast< int >( dY );
		const int iZ = static_cast< int >( dZ );
		const double delX = dX - iX;
		const double delY = dY - iY;
		const double delZ = dZ - iZ;
		if( 0 <= iX && iX+1 < m_iSizeX &&
			0 <= iY && iY+1 < m_iSizeY &&
			0 <= iZ && iZ+1 < m_iSizeZ ){
			return static_cast< I >(
				((m_aaaImage[ iZ ][ iY ][iX]*(1.0f-delX) + m_aaaImage[ iZ ][ iY ][iX+1]*delX)*(1.0f-delY) +
				 (m_aaaImage[ iZ ][iY+1][iX]*(1.0f-delX) + m_aaaImage[ iZ ][iY+1][iX+1]*delX)*delY )*(1.0f-delZ) +
				((m_aaaImage[iZ+1][ iY ][iX]*(1.0f-delX) + m_aaaImage[iZ+1][ iY ][iX+1]*delX)*(1.0f-delY) +
				 (m_aaaImage[iZ+1][iY+1][iX]*(1.0f-delX) + m_aaaImage[iZ+1][iY+1][iX+1]*delX)*delY )*delZ);
		}else{
			return m_Out;
		}
	}
	I GetCubic( double dX, double dY, double dZ )const{
		//長くなるので省略orz
	}
private:
	I*** m_aaaImage;
	int m_iSizeX, m_iSizeY, m_iSizeZ;
	Type m_type;
	I m_Out;
};

オブジェクト指向の初心者みたいなクラスになったがこの設計だと欠点がある。まず、switchによる比較が毎回あり、速度に影響がでる。もう一つが致命的で、補間はなにもこの3種類だけではない。世に出回っているだけでも多数あり、場合によっては自作することもあるだろう。そういったものに対応できない。
というわけで拡張可能なようにもうちょっとOOPっぽく作ってみた。

template< class I >
class Interpolation2D
{
public:
	Interpolation2D( I*** aaaImage, int iSizeX, int iSizeY, int iSizeZ, I out = 0 ):
		m_aaaImage(aaaImage),m_iSizeX(iSizeX),m_iSizeY(iSizeY),m_iSizeZ(iSizeZ),m_Out(out){}

	I operator()( double dX, double dY, double dZ )const{
		return Get( dX, dY, dZ );
	};
protected:
	virtual I Get( double dX, double dY, double dZ ) = 0;

	I*** m_aaaImage;
	int m_iSizeX, m_iSizeY, m_iSizeZ;
	I m_Out;
};

template< class I >
class Interpolation2D_NN:public Interpolation2D< I >
{
protected:
	virtual I Get( double dX, double dY, double dZ ){
		const int iX = static_cast< int >( dX + 0.5 );
		const int iY = static_cast< int >( dY + 0.5 );
		const int iZ = static_cast< int >( dZ + 0.5 );
		if( 0 <= iX && iX < m_iSizeX &&
			0 <= iY && iY < m_iSizeY &&
			0 <= iZ && iZ < m_iSizeZ ){
			return m_aaaImage[iZ][iY][iX];
		}else{
			return m_Out;
		}
	}
};

template< class I >
class Interpolation2D_Linear:public Interpolation2D< I >
{
protected:
	virtual I Get( double dX, double dY, double dZ ){
		//やっぱり長くなるので(ry
	}
};

典型的なTemplateMethodPatternですね。必要な補間は継承して作成すればよい。
蛇足だが、これを使うクラスとしてtransformationみたいなクラスを作って、継承先にAffineTransformationとかその他非線形な画像変形クラスを作れば典型的なBrigePatternに！
でもまあこれで拡張性は解決しているんですが、今度は中間に仮想関数を挟んでしまったためにパフォーマンスに影響が。んじゃあどうすればいいかというと、関数オブジェクトを使うんです。

template< class I, class O >
class Interpolation2D
{
public:
	Interpolation2D( I*** aaaImage, int iSizeX, int iSizeY, int iSizeZ ):
		m_aaaImage(aaaImage),m_iSizeX(iSizeX),m_iSizeY(iSizeY),m_iSizeZ(iSizeZ){}

	I operator()( double dX, double dY, double dZ )const{
		return m_Obj( m_aaaImage, m_iSizeX, m_iSizeY, m_iSizeZ, dX, dY, dZ );
	};
private:
	I*** m_aaaImage;
	int m_iSizeX, m_iSizeY, m_iSizeZ;
	const O m_Obj;
};

template< class I, I out >
struct Interpolation2D_NN{
	inline I operator( I*** aaaImage, int iSizeX, int iSizeY, int iSizeZ,
		double dX, double dY, double dZ )const{
		const int iX = static_cast< int >( dX + 0.5 );
		const int iY = static_cast< int >( dY + 0.5 );
		const int iZ = static_cast< int >( dZ + 0.5 );
		if( 0 <= iX && iX < iSizeX &&
			0 <= iY && iY < iSizeY &&
			0 <= iZ && iZ < iSizeZ ){
			return aaaImage[iZ][iY][iX];
		}else{
			return out;
		}
	}
};

template< class I, I out >
struct Interpolation2D_NN{
	inline I operator( I*** aaaImage, int iSizeX, int iSizeY, int iSizeZ,
		double dX, double dY, double dZ )const{
		//例によって省略
	}
};

こうすればインライン展開されて関数呼び出しがなくなり高速になります。かつ、関数オブジェクトを増やせばどんな補間も使えます。
「最初から関数だけで直接叩けばいいんじゃね？」と思われる方もいるかも知れません。

template< class I, I out >
inline I Interpolation2D_NN( I*** aaaImage, int iSizeX, int iSizeY, int iSizeZ,
	double dX, double dY, double dZ )const{
	const int iX = static_cast< int >( dX + 0.5 );
	const int iY = static_cast< int >( dY + 0.5 );
	const int iZ = static_cast< int >( dZ + 0.5 );
	if( 0 <= iX && iX < iSizeX &&
		0 <= iY && iY < iSizeY &&
		0 <= iZ && iZ < iSizeZ ){
		return aaaImage[iZ][iY][iX];
	}else{
		return out;
	}
};

確かに補間単体のみならそれでもいいかもしれません。しかし最初に少し触れたとおり、補間は通常はアフィン変換とセットで使われます。それどころかアフィン変換に限らず一般的な変形に関してもやはり補間は必要となるのでそうなってくると威力が発揮されます。アフィン変換クラスのテンプレート引数として補間の関数オブジェクトを渡せば関数呼び出しなしでつくることが出来ます。

というわけで、関数オブジェクトの入門みたいな内容になりました（＾＾；
で、今回まだちょっと欠点があります。画像範囲のチェックの為に画像サイズを受け取っていますが使う側としては画像の内部であることが保証されている場合、範囲チェックの分だけ重くなり、さらに引数が多い分だけ遅くなります。その分まで減らせる実装を次回紹介します・・・がいつになるやらorz