最小値フィルタ

[usao]

●やりたいこと：
画像（２次元配列的に値が並んでいるデータ）上の各位置(x,y)に関して，
その近傍　：　上下左右r[pixel]分，すなわち，(x,y)を中心とした 縦横共に r*2+1 サイズの領域
の最小値を求める（→その結果を別の画像の(x,y)に出力する）

最もシンプル(?)に書くと↓のような感じです．

コード:

//※画像バッファは入力出力共にW*H画素分用意されているものとする
const int W = 画像幅:
const int H = 画像高さ;

//画素アクセス
inline unsigned char &SRC(x,y){	return 入力画像データバッファの(x,y)の位置;	}
inline unsigned char &DST(x,y){	return 出力画像データバッファの(x,y)の位置;	}

//シンプルな実装
void MinFilter( int r )
{
	r = (std::max)( abs( r ), 1 );	//※rは1以上の整数
	
	for( int y=0; y<H; y++ )
	{
		const int top = (std::max)( y-r, 0 );
		const int bottom = (std::min)( y+r, H-1 );
		
		for( int x=0; x<W; x++ )
		{
			const int left = (std::max)( x-r, 0 );
			const int right = (std::min)( x+r, W-1 );
			
			//(left,top)-(right,bottom)矩形内の最小値を見つける
			unsigned char MinVal = 255;
			for( int yy=top; yy<=bottom; yy++ )
			{
				for( int xx=left; xx<=right; xx++ )
				{
					if( SRC(xx,yy) < MinVal )
					{	MinVal = SRC(xx,yy);	}
				}
			}
			
			//結果格納
			DST(x,y) = MinVal;
		}
	}
}

●知りたいこと
効率の良いアルゴリズムが知りたいのです．
一般的な処理だと思うので，何か定石的な方法があるのではないか？と思うのですが
検索が下手なのか，そういったものを発見できないのです…
御存じの方おられましたら，御教授願いたくお願い申し上げます．

なお，コードや疑似コード的なものを示していただける場合，C,C++ な記述を用いていただけると助かります．

[usao]

このくらいのことなら自分でもすぐ思いつくのですが，
よりうまい方法が無いものか，と．

コード:

void MinFilterB( int r )
{
    r = (std::max)( abs( r ), 1 );  //※rは1以上の整数

    unsigned char *Buff = new unsigned char[W];
    
    for( int y=0; y<H; y++ )
    {
        const int top = (std::max)( y-r, 0 );
        const int bottom = (std::min)( y+r, H-1 );

        {//Buff[]に各X座標の幅1pixelな縦長領域の最小値を格納
            for( int i=0; i<W; i++ )Buff[i]=255;

            for( int yy=top; yy<=bottom; yy++ )
            {
                for( int x=0; x<W; x++ )
                {
                    if( SRC(x,yy) < Buff[x] ){    Buff[x]=SRC(x,yy);    }
                }
            }
        }
        
        //Buff[]を用いて(x,y)での最小を求む
        for( int x=0; x<W; x++ )
        {
            const int left = (std::max)( x-r, 0 );
            const int right = (std::min)( x+r, W-1 );
            
            //結果格納
            DST(x,y) = *( std::min_element(Buff+left, Buff+right+1) );
        }
    }

    delete[] Buff;
}

[みけCAT]

まず、次のような形式の入出力ができるように、提示されたサンプルを改造しました。
●入力
W : 入力されるデータの幅
H : 入力されるデータの高さ
r : 最小値を取る範囲(課題で指定されたもの)
data(i,j) : SRC(i,j)に代入するデータ

コード:

W H r
data(0,0) data(1,0) ... data(W-1,0)
data(0,1) data(1,1) ... data(W-1,1)
...
data(0,H-1) data(1,H-1) ... data(W-1,H-1)

●出力
data(i,j) : DST(i,j)に代入するデータ

コード:

data(0,0) data(1,0) ... data(W-1,0)
data(0,1) data(1,1) ... data(W-1,1)
...
data(0,H-1) data(1,H-1) ... data(W-1,H-1)

●サンプル入力

コード:

3 3 1
1 2 3
4 5 6
7 8 9

●サンプル出力

コード:

1 1 2
1 1 2
4 4 5

改造したコードがこちらになります。

► スポイラーを表示

コード:

#include <cstdio>
#include <algorithm>
#include <cstdlib>

//※画像バッファは入力出力共にW*H画素分用意されているものとする
//const int W = 画像幅:
//const int H = 画像高さ;
int W,H;

unsigned char src[4096][4096];
unsigned char dest[4096][4096];

//画素アクセス
//inline unsigned char &SRC(x,y){ return 入力画像データバッファの(x,y)の位置;  }
//inline unsigned char &DST(x,y){ return 出力画像データバッファの(x,y)の位置;  }
#define SRC(x,y) src[(y)][(x)]
#define DST(x,y) dest[(y)][(x)]
 
//シンプルな実装
void MinFilter( int r )
{
    r = (std::max)( abs( r ), 1 );  //※rは1以上の整数
    
    for( int y=0; y<H; y++ )
    {
        const int top = (std::max)( y-r, 0 );
        const int bottom = (std::min)( y+r, H-1 );
        
        for( int x=0; x<W; x++ )
        {
            const int left = (std::max)( x-r, 0 );
            const int right = (std::min)( x+r, W-1 );
            
            //(left,top)-(right,bottom)矩形内の最小値を見つける
            unsigned char MinVal = 255;
            for( int yy=top; yy<=bottom; yy++ )
            {
                for( int xx=left; xx<=right; xx++ )
                {
                    if( SRC(xx,yy) < MinVal )
                    {   MinVal = SRC(xx,yy);    }
                }
            }
            
            //結果格納
            DST(x,y) = MinVal;
        }
    }
}

int main(void) {
	int r;
	scanf("%d%d%d",&W,&H,&r);
	for(int i=0;i<H;i++) {
		for(int j=0;j<W;j++) {
			int temp;
			scanf("%d",&temp);
			src[i][j]=(unsigned char)(temp&0xFF);
		}
	}
	MinFilter(r);
	for(int i=0;i<H;i++) {
		for(int j=0;j<W;j++) {
			printf("%d%c",dest[i][j],j+1<W?' ':'\n');
		}
	}
	return 0;
}

そして、肝心の考え方としては、セグメント木を使います。
まず、各行ごとに、横は指定された範囲、縦はその行だけという範囲の最小値を求め、一旦DSTに代入しておきます。
次に、各列ごとに、横はその列だけ、縦は指定された範囲という範囲の最小値をDSTから求め、結果をDSTに代入します。
これでDSTに求める結果が代入されます。

コード:

#include <cstdio>

#define W_MAX 4096
#define H_MAX 4096

int W,H;

unsigned char src[4096][4096];
unsigned char dest[4096][4096];

#define SRC(x,y) (src[(y)][(x)])
#define DST(x,y) (dest[(y)][(x)])

/**
 * 指定した範囲の最小値を求めるセグメント木を構築する
 * @param segtree セグメント木データへのポインタ
 * @param segtreeSize セグメント木のサイズ(2の非負整数乗でないといけない)
 */
void MakeSegmentTree(unsigned char* segtree,int segtreeSize) {
	// セグメント木を構築する
	for(int i=segtreeSize-2;i>=0;i--) {
		segtree[i]=segtree[i*2+1];
		if(segtree[i*2+2]<segtree[i])segtree[i]=segtree[i*2+2];
	}
}

/**
 * セグメント木から指定した範囲の最小値を求める
 * @param segtree セグメント木のデータ
 * @param segtreeSize セグメント木のサイズ(2の非負整数乗でないといけない)
 * @param start 範囲の開始点(含まれる)
 * @param end 範囲の終了点(含まれない)
 * @return 最小値
 */
unsigned char GetMinFromSegmentTree(const unsigned char* segtree,
		int segtreeSize,int start,int end,int now=0,int kstart=0,int kend=-1) {
	if(kend<0)kend=segtreeSize;
	// 範囲外なので無視する
	if(end<=kstart || kend<=start)return 255;
	// 範囲全てが含まれるのでそのまま返す
	if(start<=kstart && kend<=end)return segtree[now];
	// 下の階層に丸投げする
	unsigned char candidate1,candidate2;
	candidate1=GetMinFromSegmentTree(segtree,segtreeSize,start,end,
			now*2+1,kstart,(kstart+kend)/2);
	candidate2=GetMinFromSegmentTree(segtree,segtreeSize,start,end,
			now*2+2,(kstart+kend)/2,kend);
	return candidate1<candidate2?candidate1:candidate2;
}

void MinFilter(int r) {
	if(r<1)r=1;
	unsigned char* segtree;
	// セグメント木のサイズを決定する
	int segtreeSize;
	segtreeSize=(W>=H?W:H);
	for(;;) {
		int nextss=segtreeSize-(segtreeSize & (-segtreeSize));
		if(nextss==0)break;
		segtreeSize=nextss;
	}
	if(segtreeSize!=r)segtreeSize*=2;
	segtree=new unsigned char[segtreeSize*2];
	unsigned char* segtreeData=segtree+segtreeSize-1;

	// 横方向の最小値を求める
	for(int y=0;y<H;y++) {
		// セグメント木にデータを設定する
		for(int x=0;x<W;x++)segtreeData[x]=SRC(x,y);
		for(int x=W;x<segtreeSize;x++)segtreeData[x]=255;
		MakeSegmentTree(segtree,segtreeSize);
		// セグメント木から最小値を求める
		for(int x=0;x<W;x++) {
			int start=x-r;
			int end=x+r+1;
			if(start<0)start=0;
			if(end>W)end=W;
			DST(x,y)=GetMinFromSegmentTree(segtree,segtreeSize,start,end);
		}
	}
	// 横方向の最小値の最小値を、縦方向で求める
	for(int x=0;x<W;x++) {
		// セグメント木にデータを設定する
		for(int y=0;y<H;y++)segtreeData[y]=DST(x,y);
		for(int y=H;y<segtreeSize;y++)segtreeData[y]=255;
		MakeSegmentTree(segtree,segtreeSize);
		// セグメント木から最小値を求める
		for(int y=0;y<H;y++) {
			int start=y-r;
			int end=y+r+1;
			if(start<0)start=0;
			if(end>H)end=H;
			DST(x,y)=GetMinFromSegmentTree(segtree,segtreeSize,start,end);
		}
	}

	delete[] segtree;
}

int main(void) {
	int r;
	scanf("%d%d%d",&W,&H,&r);
	for(int i=0;i<H;i++) {
		for(int j=0;j<W;j++) {
			int temp;
			scanf("%d",&temp);
			src[i][j]=(unsigned char)(temp&0xFF);
		}
	}
	MinFilter(r);
	for(int i=0;i<H;i++) {
		for(int j=0;j<W;j++) {
			printf("%d%c",dest[i][j],j+1<W?' ':'\n');
		}
	}
	return 0;
}

添付のテストケースでテストしたところ、最初のプログラムで答えが出せなかったcase16.txtを除き、出力が一致しました。

[みけCAT]

「区間 スライド最小値 O(N)」でググッたところ、以下のサイトが見つかりました。
http://d.hatena.ne.jp/kyuridenamida/tou ... 1325870418
正方形領域の最小値を求めるコードを含むので、これが応用できるかもしれません。

[usao]

ありがとうございます．

縦と横で分けて２回走査するのは他のフィルタでもやる手ですね（なのになぜ自分は気づかなかったのか！？）

ちょっと今現在時間とれないのですが
縦横に分離するのを基本方針として，必要時間が
　セグメント木　＜　その他てきとーな方法
となるのかどうかなど，時間が取れるときに試してみようかと思います．


＃思いついているてきとーな方法：
　ある(x,y)を中心とした領域内の最小値の座標が(x',y')のとき，(x+1,y)を中心とした領域内に(x',y')が含まれているなら
　２領域の共通領域内の要素比較を再度やる必要はないから比較回数を減らせるかな　みたいな．
　（ワーストケースでは全く働かないだろうけど）

[みけCAT]

http://techtipshoge.blogspot.jp/2012/03/blog-post.html
を参考に実装してみました。
セグメント木のプログラムより速い可能性が高いと思います。
(test16.txtでテストした結果：セグメント木→約690ms　今回のプログラム→約590ms)

コード:

#include <cstdio>

#define W_MAX 4096
#define H_MAX 4096

int W,H;

unsigned char src[4096][4096];
unsigned char dest[4096][4096];

#define SRC(x,y) (src[(y)][(x)])
#define DST(x,y) (dest[(y)][(x)])

// バッファに新しい値を挿入する
void insertValue(unsigned char* buffer,int bufferStart,int& bufferEnd,
		unsigned char value) {
	while(bufferStart<bufferEnd && buffer[bufferEnd-1]>value)bufferEnd--;
	buffer[bufferEnd++]=value;
}

// バッファの先頭のデータが消えるデータと一致すれば、バッファの先頭のデータを消す
void deleteValue(unsigned char* buffer,int& bufferStart,int bufferEnd,
		unsigned char value) {
	if(bufferStart<bufferEnd && buffer[bufferStart]==value)bufferStart++;
}

void MinFilter(int r) {
	if(r<1)r=1;

	int bufferSize=(W>=H?W:H);
	unsigned char* buffer=new unsigned char[bufferSize];
	unsigned char* backup=new unsigned char[H];
	int bufferStart,bufferEnd;

	// 横方向の最小値を求める
	for(int y=0;y<H;y++) {
		bufferStart=bufferEnd=0;
		for(int x=0;x<r;x++)insertValue(buffer,bufferStart,bufferEnd,SRC(x,y));
		for(int x=0;x<W;x++) {
			if(x-r-1>=0)deleteValue(buffer,bufferStart,bufferEnd,SRC(x-r-1,y));
			if(x+r<W)insertValue(buffer,bufferStart,bufferEnd,SRC(x+r,y));
			DST(x,y)=buffer[bufferStart];
		}
	}
	// 横方向の最小値の最小値を、縦方向で求める
	for(int x=0;x<W;x++) {
		bufferStart=bufferEnd=0;
		for(int y=0;y<r;y++)insertValue(buffer,bufferStart,bufferEnd,DST(x,y));
		for(int y=0;y<H;y++) {
			backup[y]=DST(x,y);
			if(y-r-1>=0)deleteValue(buffer,bufferStart,bufferEnd,backup[y-r-1]);
			if(y+r<H)insertValue(buffer,bufferStart,bufferEnd,DST(x,y+r));
			DST(x,y)=buffer[bufferStart];
		}
	}

	delete[] buffer;
	delete[] backup;
}

int main(void) {
	int r;
	scanf("%d%d%d",&W,&H,&r);
	for(int i=0;i<H;i++) {
		for(int j=0;j<W;j++) {
			int temp;
			scanf("%d",&temp);
			src[i][j]=(unsigned char)(temp&0xFF);
		}
	}
	MinFilter(r);
	for(int i=0;i<H;i++) {
		for(int j=0;j<W;j++) {
			printf("%d%c",dest[i][j],j+1<W?' ':'\n');
		}
	}
	return 0;
}

[usao]

自分の思いつきコードと　No.6 のコードの処理速度を test16.txt で比較してみました．
（Ｎ回連続して処理させると１，２回目あたり(?)は処理時間がすこしかかるようなので，
　処理はループで15回連続して行い，6回目～15回目の10回の試行の平均を取った）

No.6コード　：　4.93035[ms]
私のコード　：　5.14685[ms]

…くやしいです！
が，4096*4096とられていたバッファサイズを 512*512 に変更して見ると以下のように逆転．

No.6コード　：　3.66212[ms]
私のコード　：　2.70538[ms]

なんだろう？　バッファへのアクセス順の違いとかが影響しているのかな？


本件は定石的な方法があるかどうか的な内容だったのですが，その後も探してみたもののそういった明確な情報は無いようでした．
が，頂いた情報およびアイデアを参考に最初のシンプルな２つに比べるとかなりの高速化は達成できましたし，
とりあえず本件を閉じる意味で解決を付けさせていただきます．
ありがとうございました．

↓私のコード

コード:

void MinFilterC_( int r )
{
	if( r<1 )r=1;
	
	//幅(r*2+1)×高(1) の領域の最小値をDSTにいれる
	for( int y=0; y<H; y++ )
	{
		int valid_count = 0;
		unsigned char min_val = 255;

		for( int x=0; x<W; x++ )
		{
			const int right = (std::min)( x+r, W-1 );

			if( valid_count > 0 )
			{
				const unsigned char val = SRC(right,y);
				if( min_val >= val )
				{
					min_val = val;
					valid_count = right + r - x;
				}
				else
				{	--valid_count;	}
			}
			else
			{
				min_val = 255;
				int min_xx=right;
				const int left = (std::max)( x-r, 0 );

				for( int xx=left; xx<=right; xx++ )
				{
					const unsigned char val = SRC(xx,y);
					if( min_val >= val )
					{	
						min_val = val;
						min_xx = xx;
					}
				}
				valid_count = min_xx + r - x;
			}

			DST(x,y) = min_val;
		}
	}
		
	//最終結果をつくる
	unsigned char *Buff = new unsigned char[H];	//１列分バッファ

	for( int x=0; x<W; x++ )
	{
		//列バッファにDST(x,*)をコピー
		for( int y=0; y<H; y++ )
		{	Buff[y] = DST(x, y);	}

		int valid_count = 0;
		unsigned char min_val = 255;

		for( int y=0; y<H; y++ )
		{
			const int bottom = (std::min)( y+r, H-1 );

			if( valid_count > 0 )
			{
				if( min_val >= Buff[bottom] )
				{
					min_val = Buff[bottom];
					valid_count = bottom + r - y;
				}
				else
				{	--valid_count;	}
			}
			else
			{
				min_val = 255;
				int min_yy = bottom;

				const int top = (std::max)( y-r, 0 );
				for( int yy=top; yy<=bottom; yy++ )
				{
					if( min_val >= Buff[yy] )
					{	
						min_val = Buff[yy];
						min_yy = yy;
					}
				}
				valid_count = min_yy + r - y;
			}

			DST(x,y) = min_val;
		}
	}

	delete[] Buff;
}

[みけCAT]

本件は定石的な方法があるかどうか的な内容だったのですが，その後も探してみたもののそういった明確な情報は無いようでした．

プログラミングコンテストチャレンジブック（蟻本）の300～301ページに、今回のNo.6で使用した方法が載っています。
(端の処理や、蟻本ではデックに要素のインデックスを入れているという点が異なりますが、基本は同じです)
(No.6では、この方法を横と縦の2回使用しています)
即ち、これが定石的な方法であると考えることができると思います。

[usao]

>即ち、これが定石的な方法
と言えるほど有名な本なのでしょうか．
しかし，なんか楽しそうな本ですね．
アルゴリズムの解説とかが詳しく述べられていたりするのでしょうか．