学步园

参考：桑卡，《图像处理分析与机器学习》

中值滤波每次都需要对kernel内的像素进行排序，效率较低。kernel每次向右移动时，最左一列消失，同时增加最右一列，所以可以对kernel内的像素建立一个直方图，计算median value以及小于等于median value的像素个数nm。kernel每次更新时，重新计算median value和nm。算法如下：

#include "opencv2/imgproc/imgproc.hpp"
#include "opencv2/highgui/highgui.hpp"
#include <iostream>
#include <cmath>

using namespace cv;
using namespace std;

void MyEfficientMedianFiltering(Mat& src, Mat& dst, int ksize);

int main( int argc, char** argv )
{
  // load image
  char* imageName = "images/Lenna_noise.jpg";

  Mat image;
  image = imread(imageName,1);
  if(!image.data)
  {
	  cout << "No image data" << endl;
	  getchar();
	  return -1;
  }

  // show image
  namedWindow("image", CV_WINDOW_AUTOSIZE);
  imshow("image", image);
  cvtColor( image, image, CV_RGB2GRAY );
  
  Mat dst;
  MyEfficientMedianFiltering(image, dst, 1);
 
  namedWindow("out_image", CV_WINDOW_AUTOSIZE);
  imshow("out_image", dst);
  waitKey(0);

  return 0;
}


void MyEfficientMedianFiltering(Mat& src, Mat& dst, int ksize)
{
	CV_Assert(src.data);
	CV_Assert(src.depth() == CV_8U);

	// create dst image
	dst.create(src.size(), src.type());
	
	for( int i = ksize; i < dst.rows-ksize; ++i )
	{
		uchar* p_dst = dst.ptr<uchar>(i);

		int median = 0;// median value
		int nm = 0;    // number of elements whose value is less than median

		int hist[256] = {0};
		for ( int j = ksize; j < dst.cols-ksize; ++j )
		{
			if( j == ksize )
			{
				//construct histogram
				for( int tmp_i = -ksize; tmp_i <= ksize; tmp_i++ )
				{
					uchar* p_src = src.ptr<uchar>(i+tmp_i);
					for( int tmp_j = -ksize; tmp_j <= ksize; tmp_j++ )
						hist[p_src[j+tmp_j]] += 1;	
				}

				// calculate median and nm
				nm = 0;
				int element_nums = (2*ksize+1)*(2*ksize+1);
				int gray_level = 0;
				for( gray_level = 0; gray_level <= 255; gray_level++ )
				{
					if( hist[gray_level] > 0 )
						nm += hist[gray_level];
					if( nm >= element_nums/2 )
						break;
				}
				median = gray_level;

				p_dst[j] = (uchar)median;
				continue;
			}

			// update histogram
			for( int tmp_i = -ksize; tmp_i <= ksize; tmp_i++ )
			{
				uchar* p_src = src.ptr<uchar>(i+tmp_i);
				hist[p_src[j-ksize-1]] -= 1;
				if( p_src[j-ksize-1] <= median )
					nm --;
				hist[p_src[j+ksize]] += 1;
				if( p_src[j+ksize] <= median )
					nm ++;					
			}

			// calculate new median and nm
			int element_nums = (2*ksize+1)*(2*ksize+1);
			while( nm > element_nums/2 )
			{	
				nm -= hist[median];
				median --;
			}
			while( nm < element_nums/2 )
			{
				median ++;
				nm += hist[median];
			}

			p_dst[j] = (uchar)median;

		}
	}
}