中值滤波Median filtering

背景简述

在概率论中，中值将概率分布的高半部分与低半部分分开，对一个随机变量x 而言，x < M 的概率为0.5。对于有限实数集，其排序后中间的数值即为它的中值。集合的大小通常取奇数，以保证中值的唯一性。中值滤波是一种减小边缘模糊的非线性平滑的方法。

基本理论

中值滤波(Median filtering) 的思想是用邻域中亮度的中值取代图像当前的点。邻域中亮度的中值不受个别噪声毛刺的影响。因此，中值平滑相当好地消除了冲击噪声。更进一步，由于中值滤波并不明显地模糊边缘，因此可以迭代使用。显然，在每个像素位置上都要对一个矩形内部的所有像素进行排序，这样的开销很大。

因此，一个更有效方法是注意到当窗口沿着行移一列时，窗口内容的变化只是舍去了最左边的列取代为一个新的右侧列，对于M x N 的列的中值窗口，mn - 2*m个像素没有变化，并不需要重新排序。算法如下所示：

算法：高效的中值滤波

1. 置 th = nm/2;如果m和n都为奇数，则对th取整，这样我们总是可以避免不必要的浮点数运算。
2. 将窗口移至一个新行的开始，对其内容排序。建立窗口像素的直方图H，确定其中值m ，记下亮度小于或等于m 的像素数目Nm.
3. 对于左列亮度是Pg的每一个像素p,做： H[Pg] = H[Pg]-1.进一步，如果Pg< m,置Nm = Nm - 1.
4. 将窗口右移一列，对于右列亮度是Pg的每一个像素p,做： H[Pg] = H[Pg]-1.如果Pg< m,置Nm = Nm + 1.
5. 如果Nm = t 则跳转到至步骤8.
6. 如果Nm > t 则跳转到至步骤7.重复m  = m + 1Nm = Nm + H[m] 直到Nm大于等于t，则跳转到至步骤8.
7.(此时有Nm > t)重复m  = m － 1Nm = Nm － H[m] 直到Nm 小于等于 t.
8. 如果窗口的右侧列不是图像的右边界，跳转至步骤3.
9  如果窗口的低行不是图像的下边界，跳转至步骤2.

上述，我们叙述了中值滤波的算法的基本步骤思想。现给出参考代码，仅作参考。

参考代码

median_filter.m文件

function d=median_filter(x,n)
[height, width]=size(x);
x1=double(x);
x2=x1;
for i=1:height-n+1  for j=1:height-n+1  c=x1(i:i+(n-1),j:j+(n-1));   e=c(1,:);       for u=2:n  e=[e,c(u,:)];     end  mm=median(e);       x2(i+(n-1)/2,j+(n-1)/2)=mm; end
end
d=uint8(x2);

C版Median Filter

void MedianFilter(unsigned char *srcdata, unsigned char *resdata,int nHeight, int nWidth, int nR)
{ int i,j,m,n,r; unsigned tmp; unsigned char* mdata = new unsigned char[(2*nR+1)*(2*nR+1)];for (i=0;i<nRows;i++)for (j=0;j<nCols;j++){if((i<nR) || (i>nHeight-nR-1) ||(j<nR) || (j>nWidth-nR-1))resdata[i*nWidth+j] = 0;else { for(m=-nR;m<=nR;m++)for(n=-nR;n<=nR;n++) mdata[(m+nR)*(2*nR+1)+n+nR]=srcdata[(i+m)*nWidth+(j+n)]; for(m=0;m<(2*nR+1)*(2*nR+1)-2;m++){ r = 1; for(n=m+1;n<(2*nR+1)*(2*nR+1);n++){ if (mdata[n]<mdata[n+1]){tmp = mdata[n];mdata[n]= mdata[n+1];mdata[n+1] = tmp;r=0; } }if (r) break; } resdata[i*nWidth+j] = mdata[nR*(2*nR+1)+nR]; } }
}

OpenCV版Median Filter

#coding=utf-8
import cv2
import numpy as np    def salt(img, n):    for k in range(n):    i = int(np.random.random() * img.shape[1]);    j = int(np.random.random() * img.shape[0]);    if img.ndim == 2:     img[j,i] = 255    elif img.ndim == 3:     img[j,i,0]= 255    img[j,i,1]= 255    img[j,i,2]= 255    return img   img = cv2.imread("test.jpg", 0)
result = salt(img, 100)
median = cv2.medianBlur(result, 5)  cv2.imshow("Median", median)
cv2.imwrite("median.jpg",median)cv2.waitKey(0)

由于中值滤波不会处理最大和最小值，所以就不会受到噪声的影响。相反，如果直接采用blur进行均值滤波，则不会区分这些噪声点，滤波后的图像会受到噪声的影响。中值滤波器在处理边缘也有优势。但中值滤波器会清除掉某些区域的纹理（如背景中的树）。

CUDA版Median Filter

#define I_Width 640
#define I_Height 480
#define BLOCK_X  16
#define BLOCK_Y  16// Exchange trick
#define s2(a,b)            { tmp = a; a = min(a,b); b = max(tmp,b); }
#define mn3(a,b,c)         s2(a,b); s2(a,c);
#define mx3(a,b,c)         s2(b,c); s2(a,c);#define mnmx3(a,b,c)       mx3(a,b,c); s2(a,b); // 3 exchanges
#define mnmx4(a,b,c,d)     s2(a,b);s2(c,d);s2(a,c);s2(b,d);  // 4 exchanges
#define mnmx5(a,b,c,d,e)   s2(a,b);s2(c,d);mn3(a,c,e);mx3(b,d,e); // 6 exchanges
#define mnmx6(a,b,c,d,e,f) s2(a,d);s2(b,e);s2(c,f);mn3(a,b,c);mx3(d,e,f);//7 exchanges#define SMEM(x,y)  smem[(x)+1][(y)+1]
#define IN(x,y)    d_in[(y)*nx + (x)]__global__ void Kernel_MedianFilter(int nx, int ny, char *d_out, char *d_in)
{int tx = threadIdx.x, ty = threadIdx.y;// guards: is at boundary?bool is_x_top = (tx == 0), is_x_bot = (tx == BLOCK_X-1);bool is_y_top = (ty == 0), is_y_bot = (ty == BLOCK_Y-1);__shared__ char smem[BLOCK_X+2][BLOCK_Y+2];// clear out shared memory (zero padding)if (is_x_top)           SMEM(tx-1, ty  ) = 0;else if (is_x_bot)      SMEM(tx+1, ty  ) = 0;if (is_y_top) {         SMEM(tx  , ty-1) = 0;if (is_x_top)       SMEM(tx-1, ty-1) = 0;else if (is_x_bot)  SMEM(tx+1, ty-1) = 0;} else if (is_y_bot) {  SMEM(tx  , ty+1) = 0;if (is_x_top)  SMEM(tx-1, ty+1) = 0;else if (is_x_bot)  SMEM(tx+1, ty+1) = 0;}// guards: is at boundary and still more image?int x = blockIdx.x * blockDim.x + tx;int y = blockIdx.y * blockDim.y + ty;is_x_top &= (x > 0); is_x_bot &= (x < nx - 1);is_y_top &= (y > 0); is_y_bot &= (y < ny - 1);// each thread pulls from imageSMEM(tx  , ty  ) = IN(x  , y  ); // selfif (is_x_top)           SMEM(tx-1, ty  ) = IN(x-1, y  );else if (is_x_bot)      SMEM(tx+1, ty  ) = IN(x+1, y  );if (is_y_top) {         SMEM(tx  , ty-1) = IN(x  , y-1);if (is_x_top)       SMEM(tx-1, ty-1) = IN(x-1, y-1);else if (is_x_bot)  SMEM(tx+1, ty-1) = IN(x+1, y-1);} else if (is_y_bot){  SMEM(tx  , ty+1) = IN(x  , y+1);if (is_x_top)       SMEM(tx-1, ty+1) = IN(x-1, y+1);else if (is_x_bot)  SMEM(tx+1, ty+1) = IN(x+1, y+1);}// pull top six from shared memorychar v[6] = {SMEM(tx-1, ty-1), SMEM(tx  , ty-1), SMEM(tx+1, ty-1),SMEM(tx-1, ty  ), SMEM(tx  , ty  ), SMEM(tx+1, ty  ) };// with each pass, remove min and max values and add new valuechar tmp;mnmx6(v[0], v[1], v[2], v[3], v[4], v[5]);v[5] = SMEM(tx-1, ty+1); // add new contestantmnmx5(v[1], v[2], v[3], v[4], v[5]);v[5] = SMEM(tx  , ty+1);mnmx4(v[2], v[3], v[4], v[5]);v[5] = SMEM(tx+1, ty+1);mnmx3(v[3], v[4], v[5]);//pick the middle oned_out[y*nx + x] = v[4];
}

测试输出

矩形邻域中值滤波的主要不足是图像中细线和显著解点分遭到时损坏，如果使用其他形状的邻域是可以避免。

补充

(2013年10月30日)

在论文作者的原始代码中，可以看用到SSE2的相关函数进行处理即：HistgramAdd( ) 和HistgramSub( )函数.这里_mm_add_epi16是一组Instructions指令，编译器在编译的时候会将其编译为SSE对应的汇编代码，而由于这些指令能实现指令级别并行，比如上述_mm_add_epi16可以在同一个指令周期对8个16位数据同时进行处理，并且HistgramAdd这些函数在程序里会大量使用到，因此程序的速度能大幅提高。在此贴出作者原始代码(供大家学习)：

/** ctmf.c - Constant-time median filtering* Copyright (C) 2006  Simon Perreault** Reference: S. Perreault and P. HÃ©bert, "Median Filtering in Constant Time",* IEEE Transactions on Image Processing, September 2007.** This program has been obtained from http://nomis80.org/ctmf.html. No patent* covers this program, although it is subject to the following license:**  This program is free software: you can redistribute it and/or modify*  it under the terms of the GNU General Public License as published by*  the Free Software Foundation, either version 3 of the License, or*  (at your option) any later version.**  This program is distributed in the hope that it will be useful,*  but WITHOUT ANY WARRANTY; without even the implied warranty of*  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the*  GNU General Public License for more details.**  You should have received a copy of the GNU General Public License*  along with this program.  If not, see <http://www.gnu.org/licenses/>.** Contact:*  Laboratoire de vision et systÃ¨mes numÃ©riques*  Pavillon Adrien-Pouliot*  UniversitÃ© Laval*  Sainte-Foy, QuÃ©bec, Canada*  G1K 7P4**  perreaul@gel.ulaval.ca*//* Standard C includes */
#include <assert.h>
#include <math.h>
#include <stdlib.h>
#include <string.h>/* Type declarations */
#ifdef _MSC_VER
#include <basetsd.h>
typedef UINT8 uint8_t;
typedef UINT16 uint16_t;
typedef UINT32 uint32_t;
#pragma warning( disable: 4799 )
#else
#include <stdint.h>
#endif/* Intrinsic declarations */
#if defined(__SSE2__) || defined(__MMX__)
#if defined(__SSE2__)
#include <emmintrin.h>
#elif defined(__MMX__)
#include <mmintrin.h>
#endif
#if defined(__GNUC__)
#include <mm_malloc.h>
#elif defined(_MSC_VER)
#include <malloc.h>
#endif
#elif defined(__ALTIVEC__)
#include <altivec.h>
#endif/* Compiler peculiarities */
#if defined(__GNUC__)
#include <stdint.h>
#define inline __inline__
#define align(x) __attribute__ ((aligned (x)))
#elif defined(_MSC_VER)
#define inline __inline
#define align(x) __declspec(align(x))
#else
#define inline
#define align(x)
#endif#ifndef MIN
#define MIN(a,b) ((a) > (b) ? (b) : (a))
#endif#ifndef MAX
#define MAX(a,b) ((a) < (b) ? (b) : (a))
#endif/*** This structure represents a two-tier histogram. The first tier (known as the* "coarse" level) is 4 bit wide and the second tier (known as the "fine" level)* is 8 bit wide. Pixels inserted in the fine level also get inserted into the* coarse bucket designated by the 4 MSBs of the fine bucket value.** The structure is aligned on 16 bytes, which is a prerequisite for SIMD* instructions. Each bucket is 16 bit wide, which means that extra care must be* taken to prevent overflow.*/
typedef struct align(16)
{uint16_t coarse[16];uint16_t fine[16][16];
} Histogram;/*** HOP is short for Histogram OPeration. This macro makes an operation \a op on* histogram \a h for pixel value \a x. It takes care of handling both levels.*/
#define HOP(h,x,op) \h.coarse[x>>4] op; \*((uint16_t*) h.fine + x) op;#define COP(c,j,x,op) \h_coarse[ 16*(n*c+j) + (x>>4) ] op; \h_fine[ 16 * (n*(16*c+(x>>4)) + j) + (x & 0xF) ] op;/*** Adds histograms \a x and \a y and stores the result in \a y. Makes use of* SSE2, MMX or Altivec, if available.*/
#if defined(__SSE2__)
static inline void histogram_add( const uint16_t x[16], uint16_t y[16] )
{*(__m128i*) &y[0] = _mm_add_epi16( *(__m128i*) &y[0], *(__m128i*) &x[0] );*(__m128i*) &y[8] = _mm_add_epi16( *(__m128i*) &y[8], *(__m128i*) &x[8] );
}
#elif defined(__MMX__)
static inline void histogram_add( const uint16_t x[16], uint16_t y[16] )
{*(__m64*) &y[0]  = _mm_add_pi16( *(__m64*) &y[0],  *(__m64*) &x[0]  );*(__m64*) &y[4]  = _mm_add_pi16( *(__m64*) &y[4],  *(__m64*) &x[4]  );*(__m64*) &y[8]  = _mm_add_pi16( *(__m64*) &y[8],  *(__m64*) &x[8]  );*(__m64*) &y[12] = _mm_add_pi16( *(__m64*) &y[12], *(__m64*) &x[12] );
}
#elif defined(__ALTIVEC__)
static inline void histogram_add( const uint16_t x[16], uint16_t y[16] )
{*(vector unsigned short*) &y[0] = vec_add( *(vector unsigned short*) &y[0],*(vector unsigned short*) &x[0] );*(vector unsigned short*) &y[8] = vec_add( *(vector unsigned short*) &y[8],*(vector unsigned short*) &x[8] );
}
#else
static inline void histogram_add( const uint16_t x[16], uint16_t y[16] )
{int i;for ( i = 0; i < 16; ++i ) {y[i] += x[i];}
}
#endif/*** Subtracts histogram \a x from \a y and stores the result in \a y. Makes use* of SSE2, MMX or Altivec, if available.*/
#if defined(__SSE2__)
static inline void histogram_sub( const uint16_t x[16], uint16_t y[16] )
{*(__m128i*) &y[0] = _mm_sub_epi16( *(__m128i*) &y[0], *(__m128i*) &x[0] );*(__m128i*) &y[8] = _mm_sub_epi16( *(__m128i*) &y[8], *(__m128i*) &x[8] );
}
#elif defined(__MMX__)
static inline void histogram_sub( const uint16_t x[16], uint16_t y[16] )
{*(__m64*) &y[0]  = _mm_sub_pi16( *(__m64*) &y[0],  *(__m64*) &x[0]  );*(__m64*) &y[4]  = _mm_sub_pi16( *(__m64*) &y[4],  *(__m64*) &x[4]  );*(__m64*) &y[8]  = _mm_sub_pi16( *(__m64*) &y[8],  *(__m64*) &x[8]  );*(__m64*) &y[12] = _mm_sub_pi16( *(__m64*) &y[12], *(__m64*) &x[12] );
}
#elif defined(__ALTIVEC__)
static inline void histogram_sub( const uint16_t x[16], uint16_t y[16] )
{*(vector unsigned short*) &y[0] = vec_sub( *(vector unsigned short*) &y[0],*(vector unsigned short*) &x[0] );*(vector unsigned short*) &y[8] = vec_sub( *(vector unsigned short*) &y[8],*(vector unsigned short*) &x[8] );
}
#else
static inline void histogram_sub( const uint16_t x[16], uint16_t y[16] )
{int i;for ( i = 0; i < 16; ++i ) {y[i] -= x[i];}
}
#endifstatic inline void histogram_muladd( const uint16_t a, const uint16_t x[16],uint16_t y[16] )
{int i;for ( i = 0; i < 16; ++i ) {y[i] += a * x[i];}
}static void ctmf_helper(const unsigned char* const src, unsigned char* const dst,const int width, const int height,const int src_step, const int dst_step,const int r, const int cn,const int pad_left, const int pad_right)
{const int m = height, n = width;int i, j, k, c;const unsigned char *p, *q;Histogram H[4];uint16_t *h_coarse, *h_fine, luc[4][16];assert( src );assert( dst );assert( r >= 0 );assert( width >= 2*r+1 );assert( height >= 2*r+1 );assert( src_step != 0 );assert( dst_step != 0 );/* SSE2 and MMX need aligned memory, provided by _mm_malloc(). */
#if defined(__SSE2__) || defined(__MMX__)h_coarse = (uint16_t*) _mm_malloc(  1 * 16 * n * cn * sizeof(uint16_t), 16 );h_fine   = (uint16_t*) _mm_malloc( 16 * 16 * n * cn * sizeof(uint16_t), 16 );memset( h_coarse, 0,  1 * 16 * n * cn * sizeof(uint16_t) );memset( h_fine,   0, 16 * 16 * n * cn * sizeof(uint16_t) );
#elseh_coarse = (uint16_t*) calloc(  1 * 16 * n * cn, sizeof(uint16_t) );h_fine   = (uint16_t*) calloc( 16 * 16 * n * cn, sizeof(uint16_t) );
#endif/* First row initialization */for ( j = 0; j < n; ++j ) {for ( c = 0; c < cn; ++c ) {COP( c, j, src[cn*j+c], += r+1 );}}for ( i = 0; i < r; ++i ) {for ( j = 0; j < n; ++j ) {for ( c = 0; c < cn; ++c ) {COP( c, j, src[src_step*i+cn*j+c], ++ );}}}for ( i = 0; i < m; ++i ) {/* Update column histograms for entire row. */p = src + src_step * MAX( 0, i-r-1 );q = p + cn * n;for ( j = 0; p != q; ++j ) {for ( c = 0; c < cn; ++c, ++p ) {COP( c, j, *p, -- );}}p = src + src_step * MIN( m-1, i+r );q = p + cn * n;for ( j = 0; p != q; ++j ) {for ( c = 0; c < cn; ++c, ++p ) {COP( c, j, *p, ++ );}}/* First column initialization */memset( H, 0, cn*sizeof(H[0]) );memset( luc, 0, cn*sizeof(luc[0]) );if ( pad_left ) {for ( c = 0; c < cn; ++c ) {histogram_muladd( r, &h_coarse[16*n*c], H[c].coarse );}}for ( j = 0; j < (pad_left ? r : 2*r); ++j ) {for ( c = 0; c < cn; ++c ) {histogram_add( &h_coarse[16*(n*c+j)], H[c].coarse );}}for ( c = 0; c < cn; ++c ) {for ( k = 0; k < 16; ++k ) {histogram_muladd( 2*r+1, &h_fine[16*n*(16*c+k)], &H[c].fine[k][0] );}}for ( j = pad_left ? 0 : r; j < (pad_right ? n : n-r); ++j ) {for ( c = 0; c < cn; ++c ) {const uint16_t t = 2*r*r + 2*r;uint16_t sum = 0, *segment;int b;histogram_add( &h_coarse[16*(n*c + MIN(j+r,n-1))], H[c].coarse );/* Find median at coarse level */for ( k = 0; k < 16 ; ++k ) {sum += H[c].coarse[k];if ( sum > t ) {sum -= H[c].coarse[k];break;}}assert( k < 16 );/* Update corresponding histogram segment */if ( luc[c][k] <= j-r ) {memset( &H[c].fine[k], 0, 16 * sizeof(uint16_t) );for ( luc[c][k] = j-r; luc[c][k] < MIN(j+r+1,n); ++luc[c][k] ) {histogram_add( &h_fine[16*(n*(16*c+k)+luc[c][k])], H[c].fine[k] );}if ( luc[c][k] < j+r+1 ) {histogram_muladd( j+r+1 - n, &h_fine[16*(n*(16*c+k)+(n-1))], &H[c].fine[k][0] );luc[c][k] = j+r+1;}}else {for ( ; luc[c][k] < j+r+1; ++luc[c][k] ) {histogram_sub( &h_fine[16*(n*(16*c+k)+MAX(luc[c][k]-2*r-1,0))], H[c].fine[k] );histogram_add( &h_fine[16*(n*(16*c+k)+MIN(luc[c][k],n-1))], H[c].fine[k] );}}histogram_sub( &h_coarse[16*(n*c+MAX(j-r,0))], H[c].coarse );/* Find median in segment */segment = H[c].fine[k];for ( b = 0; b < 16 ; ++b ) {sum += segment[b];if ( sum > t ) {dst[dst_step*i+cn*j+c] = 16*k + b;break;}}assert( b < 16 );}}}#if defined(__SSE2__) || defined(__MMX__)_mm_empty();_mm_free(h_coarse);_mm_free(h_fine);
#elsefree(h_coarse);free(h_fine);
#endif
}/*** \brief Constant-time median filtering** This function does a median filtering of an 8-bit image. The source image is* processed as if it was padded with zeros. The median kernel is square with* odd dimensions. Images of arbitrary size may be processed.** To process multi-channel images, you must call this function multiple times,* changing the source and destination adresses and steps such that each channel* is processed as an independent single-channel image.** Processing images of arbitrary bit depth is not supported.** The computing time is O(1) per pixel, independent of the radius of the* filter. The algorithm's initialization is O(r*width), but it is negligible.* Memory usage is simple: it will be as big as the cache size, or smaller if* the image is small. For efficiency, the histograms' bins are 16-bit wide.* This may become too small and lead to overflow as \a r increases.** \param src           Source image data.* \param dst           Destination image data. Must be preallocated.* \param width         Image width, in pixels.* \param height        Image height, in pixels.* \param src_step      Distance between adjacent pixels on the same column in*                      the source image, in bytes.* \param dst_step      Distance between adjacent pixels on the same column in*                      the destination image, in bytes.* \param r             Median filter radius. The kernel will be a 2*r+1 by*                      2*r+1 square.* \param cn            Number of channels. For example, a grayscale image would*                      have cn=1 while an RGB image would have cn=3.* \param memsize       Maximum amount of memory to use, in bytes. Set this to*                      the size of the L2 cache, then vary it slightly and*                      measure the processing time to find the optimal value.*                      For example, a 512 kB L2 cache would have*                      memsize=512*1024 initially.*/
void ctmf(const unsigned char* const src, unsigned char* const dst,const int width, const int height,const int src_step, const int dst_step,const int r, const int cn, const long unsigned int memsize)
{/** Processing the image in vertical stripes is an optimization made* necessary by the limited size of the CPU cache. Each histogram is 544* bytes big and therefore I can fit a limited number of them in the cache.* That number may sometimes be smaller than the image width, which would be* the number of histograms I would need without stripes.** I need to keep histograms in the cache so that they are available* quickly when processing a new row. Each row needs access to the previous* row's histograms. If there are too many histograms to fit in the cache,* thrashing to RAM happens.** To solve this problem, I figure out the maximum number of histograms* that can fit in cache. From this is determined the number of stripes in* an image. The formulas below make the stripes all the same size and use* as few stripes as possible.** Note that each stripe causes an overlap on the neighboring stripes, as* when mowing the lawn. That overlap is proportional to r. When the overlap* is a significant size in comparison with the stripe size, then we are not* O(1) anymore, but O(r). In fact, we have been O(r) all along, but the* initialization term was neglected, as it has been (and rightly so) in B.* Weiss, "Fast Median and Bilateral Filtering", SIGGRAPH, 2006. Processing* by stripes only makes that initialization term bigger.** Also, note that the leftmost and rightmost stripes don't need overlap.* A flag is passed to ctmf_helper() so that it treats these cases as if the* image was zero-padded.*/int stripes = (int) ceil( (double)(width - 2*r)/(memsize /sizeof(Histogram)-2*r));int stripe_size = (int) ceil( (double)( width + stripes*2*r - 2*r )/ stripes );int i;for ( i = 0; i < width; i += stripe_size - 2*r ) {int stripe = stripe_size;/* Make sure that the filter kernel fits into one stripe. */if (i + stripe_size-2*r>= width||width-(i + stripe_size - 2*r)<2*r+1 ) {stripe = width - i;}ctmf_helper( src + cn*i, dst + cn*i, stripe, height, src_step,dst_step,r, cn,i == 0, stripe == width - i );if ( stripe == width - i ) {break;}}
}

(2013年11月2日)

代码来自网络。个人感觉不错C++版中值滤波，共大家分享。呵呵......

/**
** method to remove noise from the corrupted image by median value
* @param corrupted input grayscale binary array with corrupted info
* @param smooth output data for smooth result, the memory need to be
*                           allocated outside of the function
* @param width width of the input grayscale image
* @param height height of the input grayscale image
*/
void medianFilter (unsigned char* corrupted,unsigned char* smooth,int width,int height)
{memcpy ( smooth, corrupted, width*height*sizeof(unsigned char) );for (int j=1;j<height-1;j++){for (int i=1;i<width-1;i++){int k = 0;unsigned char window[9];for (int jj = j - 1; jj < j + 2; ++jj)for (int ii = i - 1; ii < i + 2; ++ii)window[k++] = corrupted[jj * width + ii];//   Order elements (only half of them)for (int m = 0; m < 5; ++m){int min = m;for (int n = m + 1; n < 9; ++n)if (window[n] < window[min])min = n;//   Put found minimum element in its placeunsigned char temp = window[m];window[m] = window[min];window[min] = temp;}smooth[ j*width+i ] = window[4];}}
}

参考文献

[1] S. Perreault and P. Hebert, "Median Filtering in Constant Time",IEEE Transactions on Image Processing, September 2007.

[2] Svoboda T.nKybic J., and Hlavac V. "Image Processing Analysis and Machine Vision". Thomson Engineering 2008.

[3] Aaftab Munshi , Benedict R. Gaster , Timothy G. Mattson , James Fung and Dan Ginsburg “ OpenCL Programming Guide”.

关于Image Engineering & Computer Vision的更多讨论与交流，敬请关注本博和新浪微博songzi_tea.