中值滤波器减少噪音以参数维纳滤波器减少噪音

之前的bug：

当灰度为255或者0时，出现灰度溢出的bug，导致灰度黑白颠倒，现已修复，并重新将函数改为无返回值类型，原有的带有返回图像的函数不规范，容易忘记释放空间。新的函数再最后面

经测试，我的程序计算速度比OpenCV耗时多多了，，，~~(>_<)~~ ，，真心比不起OpenCV，因为这货用intel 的IPP库了，谁让人家是亲儿子呢，，，，

---

（1）中值滤波

由于中值滤波是一个非线性滤波器，中值的求取需要排序，传统的中值滤波算法中无非是在排序上做文章，争取做到最优的快速排序，但是本文所实现的方法不同，该方法充分利用了直方图。
从编程的观点看，直方图是一种很有效的数据结构，所占内存空间很少，又能反映出图像中的灰度分布和目标特性等等，且直方图本身就是有序的。基于直方图可以很容易、很高效地得到图像中亮度、对比度、最大亮度、最小亮度及亮度中值。
本文章旨在实现高效的中值滤波（对于这个大众化的图像增强方法没有深究，也不知该方法到底源自哪一片文章，鄙人知识有幸找到了一份老师给的PDF文档，上面介绍了该方法，就抽时间实现了一下，效果和OpenCV里实现的效果是一样的，至于速度，只会快，不会慢）

（2）原理

如上图所示，在进行横向滑动窗口滤波时，窗口中的像素仅仅是丢掉了左侧一列，增加了右侧一列数据，如果丢掉中间重叠的这一部分数据，到下个窗口再重新寻址和读取数据，无疑是计算的沉重负担，所以，该算法的核心思想就是充分利用重叠部分，使用直方图来计算中值，不需要排序算法，快，且高效。

（3）算法流程

假设窗口尺寸为 N × M

1. 设置门限th=N∗M2+1
2. 将窗口移动到一个新行的开始，建立窗口像素的直方图，通过直方图确定中值 med，记下亮度小于或等于med 的像素数目到 mNum 中
3. 对于最左列的每个像素，设其灰度为 gl ，执行：（即去掉左侧）
   H[gl]=H[gl]−1
   
   若 gl≤med，则 mNum=mNum−1
4. 将窗口右移一列，对于最右列的每个像素，设其灰度为gr，执行：（即添上右侧）
   H[gr]=H[gr]+1
   
   若 gr≤med，则 mNum=mNum+1
5. 若 mNum>th ,则转 6；否则( mNum≤th )：
   while(mNum<th)
   
   {med=med+1;mNum=mNum+H[med];}
   
   转 7
6. 重复 mNum=mNum−H[med]
   
   med=med−1
   直到 mNum≤th
7. 如果窗口的右侧不是图像的右边界，则转 3
8. 若果窗口的底行不是图象的下边界，则转 2

（4）实现

/**
 * @brief Fast median filter
 * @param src
 * @param width
 * @param height
 * @param mask_width
 * @param mask_height
 * @return address of data after median blur
 */
unsigned char* mxnMedianBlur(unsigned char* src, int width, int height, int mask_width, int mask_height)
{
    int hist[256];
    int i, j, k, nMin;
    int width1, height1;
    int th = mask_width * mask_height / 2 + 1;
    unsigned char *src1, *dst, *pCurOrg, *pCurObj, *pl, *pr;
    unsigned char med;
    src1 = boundary1(src, width, height, mask_width, mask_height, &width1, &height1);
//  cout << "new width = " << width1 << endl;
//  cout << "new height = " << height1 << endl;
//  write8bitBmp(src1, width1, height1, "extended.bmp");
    if((dst = (unsigned char*)malloc(width * height * sizeof(unsigned char))) == NULL)
    {
        printf("No space for destination data distributed !\n");
        return NULL;
    }
    memset(dst, 0, width * height * sizeof(unsigned char));
    for(i = 0; i < height; i++)
    {
        nMin = 0;
        pCurObj = dst + i * width;
        pCurOrg = src1 + (i + 1) * width1;
        getHist_local(pCurOrg, hist, width1, mask_width, mask_height);
        med = getMedianGry(hist, mask_width * mask_height);
        for(j = 0; j <= med; j++)
            nMin += hist[j];
        (*pCurObj++) = med;
        for(k = 0; k < width - 1; k++, pCurOrg++, pCurObj++)
        {
            for(j = 0, pl = pCurOrg, pr = pCurOrg + mask_width; j < mask_height; j++, pl += width1, pr += width1)
            {
                if(*pl <= med ) nMin -= 1;
                if(*pr <= med ) nMin += 1;
                hist[*pl] -= 1;
                hist[*pr] += 1;
            }
//          pCurObj++;
            if(nMin > th)
            {
                while(nMin > th)
                {
                    nMin -= hist[med];
                    med--;
                }
                *pCurObj = med;
                continue;
            }
            else
            {
                while(nMin < th)
                {
                    med++;
                    nMin += hist[med];
                }
                *pCurObj = med;
                continue;
            }
        }
    }
    free(src1);
    return dst;
}

其中boundary1(src, width, height, mask_width, mask_height, &width1, &height1);是我自己实现的根据滑动窗口尺寸对原图像的边界进行扩展的函数，扩展方式采用复制最临近的像素的方法进行，此处就不贴出来了。
getHist_local(pCurOrg, hist, width1, mask_width, mask_height)是实现的局部直方图函数，pCurOrg为要统计直方图的区域的左上角像素地址，width1是原图像的行步长。

/**
 * @brief Get median grayscale
 * @param histogram
 * @param imgSize
 * @return medianGray
 */
unsigned char getMedianGry(int *histogram, int imgSize)
{
    int sum;
    unsigned char medGry;
    for(sum = 0, medGry = 0; medGry < 256; medGry++)
    {
        sum += histogram[medGry];
        if(sum * 2 > imgSize)break;
    }
    return medGry;
}

该函数实现通过直方图寻找中值的功能，在上一个程序中有调用。

（5）结果

原图像（女神图）

中值滤波后的图像（只抽取并处理了灰度图）

处理结果与photoshop与OpenCV以及Matlab中基本上都是一致的。

改进和修复bug的版本

/**
 * @brief Fast median filter
 * @param src
 * @param width
 * @param height
 * @param mask_width
 * @param mask_height
 */
void medianBlur(unsigned char* src, int width, int height, int mask_width, int mask_height)
{
    int hist[256];
    int i, j, k, nMin;
    int width1 = width + mask_width - 1;
    int height1 = height + mask_height - 1;
    int th = mask_width * mask_height / 2 + 1;
    unsigned char *pCurOrg, *pCurObj, *pl, *pr;
    unsigned char med;
    unsigned char *src1 = Malloc(unsigned char, width1 * height1);
    //获得带边界扩展的图像副本
    boundaryExtention(src, src1, width, height, mask_width, mask_height, width1, height1);
    //清空原图像数据
    memset(src, 0, width * height * sizeof(unsigned char));
    for(i = 0; i < height; i++)
    {
        nMin = 0;
        pCurObj = src + i * width;
        pCurOrg = src1 + i * width1;
        getHistLocal(pCurOrg, hist, width1, mask_width, mask_height);
        med = getMedianGry(hist, mask_width * mask_height);
        for(j = 0; j <= med; j++)
            nMin += hist[j];
        (*pCurObj++) = med;
        for(k = 0; k < width - 1; k++, pCurOrg++, pCurObj++)
        {
            for(j = 0, pl = pCurOrg, pr = pCurOrg + mask_width; j < mask_height; j++, pl += width1, pr += width1)
            {
                if(*pl <= med ) nMin -= 1;
                if(*pr <= med ) nMin += 1;
                hist[*pl] -= 1;
                hist[*pr] += 1;
            }
            while(nMin > th)
            {
                if(med==0)
                    break;
                else
                {
                    nMin -= hist[med];
                    med--;
                }
            }
            while(nMin < th)
            {
                if(med==255)
                    break;
                else
                {
                    med++;
                    nMin += hist[med];
                }
            }
            *pCurObj = med;
        }
    }
    free(src1);
    return;
}