深入剖析图像平滑与噪声滤波

噪声

在数字图像处理中，噪声是指在图像中引入的不希望的随机或无意义的信号。它是由于图像采集、传输、存储或处理过程中的各种因素引起的。

噪声会导致图像质量下降，使图像失真或降低细节的清晰度。它通常表现为图像中随机分布的亮度或颜色变化，类似于图像上的颗粒、斑点或像素偏移。

那么噪声如何避免呢？

在图像数字处理和图像分析的算法中，为了保持一定的稳健性，一般需要先对原始图像做一定的滤波处理—图像预处理

---

图像平滑

数字图像处理中，作用于去掉图像噪声的各种滤波方法统称为图像平滑，那么图像平滑的字面意思就是使一个像素到另一个像素的灰度变化是平滑的。

那么从空间域的观点看，图像平滑实际上就是去去掉突然变大或者变小的点，用一个合适的灰度值代替该值

我们来看一个例子

这样一串数字，他的曲线我们来画一下

空间域的图像滤波围绕邻域计算展开，所以我们图像平滑也将会围绕计算和邻域展开

---

均值滤波

均值滤波就是去当前像素点为中心取一个邻域，用该邻域的所有像素值的均值作为该像素滤波后的灰度值

就比如

对于每个像素，我们计算它自身以及左右相邻的两个像素的平均值，然后用这个平均值替代原始像素的值。

例如，对于原始数据中的第一个像素9，其邻域为：

3   3   3

将这些值相加并除以3，得到新的像素值：

(3 + 3 + 3) / 3 = 9 / 3 = 3

所以，第一个像素的新值为3。依此类推，对所有像素进行相同的操作。

最终，滤波后的数据为：

3   3   5   5   5   5   7   7   7   7   9   9   9

邻域与卷积运算

邻域是什么，很简单，相邻的域

那么在图像的均值滤波中，邻域肯定不能是前面一个数字和后面一个数字，邻域一定是二维的

就比如

可用一个小的图形来表示邻域，该图形称作模板（Template）。该图形中参与运算的像素其对应位置设为“1”，不参与运算的像素其对应位置设为“0”，有时也用不填值代表0

均值滤波的公式（3-3）用模板表示为

模板就是一个小的矩阵，像素邻域中的每个像素分别与该矩阵的每个元素对应相乘，这是一种卷积运算（Convolution），所以模板常被称为卷积核，均值滤波又被看成是一种卷积运算。

在图像处理中，常用的均值滤波模板的形状有方形和圆形两种，大小有3×3和5×5两种

均值滤波特点

咱们直接看图就行

总结一下就是：

• 滤波结果是灰度值大的像素在滤波后灰度值变小、灰度值小的像素在滤波后灰度值变大，图像总灰度值之和不变，即图像的均值（亮度）不变
• 但是，正是因为大的灰度值变小、小的灰度值变大，所以图像的标准差（对比度）变小，图像变模糊
• 而且在目标和背景的边界上，滤波前后灰度值的变化尤其大，所以边界模糊的尤其厉害
• 均值滤波的邻域越大则模糊越厉害

---

基于列积分的快速均值滤波

为什么一直没有上代码？

求均值的过程就是一个多次累加和一次除法的过程，如果邻域大小为101×101，难道均值滤波需要对每个像素求邻域均值都要进行1万多次（10201）加法吗？那么一幅1080p（1920×1080）灰度图像的均值滤波岂不是要进行200多亿次加法（21,152,793,600）？

因为计算量太大了，每个像素重复计算，计算量大的吓人

如何才能优化呢？

我们其实很容易就能看出，重复计算的内容很多，而且邻域越大重复计算的内容越多

我们先来看书上的代码怎么写

void RmwAvrFilterBySumCol( BYTE *pGryImg,  //原始灰度图像
	                       int width, int height,//图像的宽度和高度
	                       int M, int N, //滤波邻域：M列N行
	                       BYTE *pResImg //结果图像
                         )
{   //没有对边界上邻域不完整的像素进行处理，可以采用变窗口的策略
	BYTE *pAdd, *pDel, *pRes;
	int halfx, halfy;
	int x, y;
	int sum,c;
	int sumCol[4096]; //约定图像宽度不大于4096

	// step.1------------初始化--------------------------//
	M = M/2*2+1; //奇数化
	N= N/2*2+1; //奇数化
	halfx = M/2; //滤波器的半径x
	halfy = N/2; //滤波器的半径y
	c = (1<<23)/(M*N); //乘法因子
	memset(sumCol, 0, sizeof(int)*width);
	for (y = 0, pAdd = pGryImg; y<N; y++)
	{
		for (x = 0; x<width; x++) sumCol[x] += *(pAdd++);
	}
	// step.2------------滤波----------------------------//
	for (y = halfy, pRes = pResImg+y*width,pDel=pGryImg; y<height-halfy; y++)
	{
		//初值
		for (sum=0,x = 0; x<M; x++) sum += sumCol[x];
		//滤波
		pRes += halfx; //跳过左侧
		for (x = halfx; x<width-halfx; x++)
		{
			//求灰度均值
			//*(pRes++)=sum/(N*M);
			*(pRes++) = (sum*c)>>23; //用整数乘法和移位代替除法
			//换列,更新灰度和
			sum -= sumCol[x-halfx]; //减左边列
			sum += sumCol[x+halfx+1]; //加右边列
		}
		pRes += halfx;//跳过右侧
		//换行,更新sumCol
		for (x = 0; x<width; x++)
		{
			sumCol[x] -= *(pDel++); //减上一行
			sumCol[x] += *(pAdd++); //加下一行
		}
	}
	// step.3------------返回----------------------------//
	return;
}

这样子就避免了重复大量的运算，我们写成C语言试试

void RmwAvrFilterBySumCol(uint8_t *pGryImg,
                           int width, int height,
                           int M, int N,
                           uint8_t *pResImg) {
    uint8_t *pAdd, *pDel, *pRes;
    int halfx, halfy;
    int x, y;
    int sum, c;
    int sumCol[4096]; // 约定图像宽度不大于4096

    // step.1------------初始化--------------------------//
    M = M / 2 * 2 + 1; // 奇数化
    N = N / 2 * 2 + 1; // 奇数化
    halfx = M / 2; // 滤波器的半径x
    halfy = N / 2; // 滤波器的半径y
    c = (1 << 23) / (M * N); // 乘法因子
    memset(sumCol, 0, sizeof(int) * width);
    for (y = 0, pAdd = pGryImg; y < N; y++) {
        for (x = 0; x < width; x++) sumCol[x] += *(pAdd++);
    }
    // step.2------------滤波----------------------------//
    for (y = halfy, pRes = pResImg + y * width, pDel = pGryImg; y < height - halfy; y++) {
        // 初值
        for (sum = 0, x = 0; x < M; x++) sum += sumCol[x];
        // 滤波
        pRes += halfx; // 跳过左侧
        for (x = halfx; x < width - halfx; x++) {
            // 求灰度均值
            // *(pRes++)=sum/(N*M);
            *(pRes++) = (sum * c) >> 23; // 用整数乘法和移位代替除法
            // 换列,更新灰度和
            sum -= sumCol[x - halfx]; // 减左边列
            sum += sumCol[x + halfx + 1]; // 加右边列
        }
        pRes += halfx; // 跳过右侧
        // 换行,更新sumCol
        for (x = 0; x < width; x++) {
            sumCol[x] -= *(pDel++); // 减上一行
            sumCol[x] += *(pAdd++); // 加下一行
        }
    }
    // step.3------------返回----------------------------//
    return;
}

我们来找张图看一下实现效果

由于没有处理边界，所以边界都是白色的像素块

• 本算法求像素的灰度均值时仅需要2个加法、2个减法、1个乘法、1个移位，共6个基本的整数运算，与邻域的大小n×m无关。此算法的巧妙之处在于采用了一个称之为“列积分”的数组sumCol。
• 本算法并没有使用除法运算sum/(N∗M)，而是使用了(sum∗c)>>23，这涉及到了一个编程技巧，即“整数除法或者浮点乘法除法变为整数乘法和移位”。

这个编程技巧涉及到了一种利用位运算来近似代替除法运算的方法，特别是在需要高效率的情况下，这种方法非常有用。让我详细解释一下：

假设我们有一个整数除法表达式 sum / divisor，其中 sum 是被除数，divisor 是除数。在计算机中，整数除法通常比其他基本运算（比如加法、减法、乘法）要慢，特别是对于一些处理器来说。

而将整数除法替换为整数乘法和移位操作的技巧可以显著提高计算效率，尤其是在一些嵌入式系统或者对计算性能要求很高的场景下。

具体来说，假设我们想计算 sum / 2^k，其中 k 是一个正整数。这个除法运算可以通过将 sum 乘以 2^k 的倒数来实现，即 sum * (1 / 2^k)。在计算机中，乘法运算通常比除法运算更快。

而 1 / 2^k 可以表示为右移 k 位，即 1 >> k。所以整个表达式 sum / 2^k 可以近似地表示为 sum * (1 >> k)。

对于给定的 k，我们可以用位移运算 sum * (1 >> k) 来代替除法运算，从而提高计算效率。

具体情况中，使用了 (sum * c) >> 23 的形式来代替除法运算。这里的 c 是一个与除数相关的常数。这样做的好处是，右移操作比除法操作更快速，因此可以提高整体的计算效率。

基于积分图的快速均值滤波

灰度图像积分图中任意一个像素s(x,y)的值是从灰度图像的左上角(0,0)与当前位置(x,y)所围成的矩形区域内的像素灰度值之和

在得到积分图后，均值滤波就变得非常简单

比如，下图中灰色区域的灰度值之和为：s( 3,3)-s(0,3)-s(3,0)+s(0,0)

那么积分图怎么来呢？

基于列积分的积分图实现

假设在y_0行上，已经知道了每列之和sumCol[x,y_0]，则积分图**s(x,y_0 )**的值为：

也就是说：

其实也很好理解

我们通过列积分来算积分图来减少计算量

void RmwDoSumGryImg( BYTE *pGryImg,  //原始灰度图像
	                 int width, //图像的宽度 
	                 int height,//图像的高度
	                 int *pSumImg //计算得到的积分图
                   )
{
	BYTE *pGry;
	int *pRes;
	int x, y;
	int sumCol[4096]; //约定图像宽度不大于4096

	memset(sumCol, 0, sizeof(int)*width);
	for (y = 0, pGry = pGryImg, pRes=pSumImg; y<height; y++)
	{
		//最左侧像素的特别处理
		sumCol[0] += *(pGry++);
		*(pRes++) = sumCol[0];
		//正常处理
		for (x = 1; x<width; x++)
		{
			sumCol[x] += *(pGry++); //更新列积分
			*(pRes++) = *(pRes-1)+sumCol[x];
		}
	}
	return;
}

其实书上也给了彩色图片的积分图计算，我们在这里也放一下，不过还是以灰度值为主

void RmwDoSumRGBImg( BYTE *pRGBImg,  //原始灰度图像
	                 int width, //图像的宽度 
	                 int height,//图像的高度
	                 int *pSumImg //计算得到的积分图
                   )
{
	BYTE *pRGB;
	int *pRes;
	int x, y;
	int sumCol[4096*3]; //约定图像宽度不大于4096

	memset(sumCol, 0, sizeof(int)*width*3);
	for (y = 0, pRGB = pRGBImg, pRes=pSumImg; y<height; y++)
	{
		//最左侧像素的特别处理
		sumCol[0] += *(pRGB++);//blue
		sumCol[1] += *(pRGB++);
		sumCol[2] += *(pRGB++);
		*(pRes++) = sumCol[0]; //blue
		*(pRes++) = sumCol[1];
		*(pRes++) = sumCol[2];
		//正常处理
		for (x = 3; x<width*3; x++)
		{
			//更新列积分
			sumCol[x] += *(pRGB++);
			//更新积分图
			*(pRes++) = *(pRes-3)+sumCol[x];
		}
	}
	return;
}

我们把灰度值的积分图函数改为C语言

void RmwDoSumGryImg(uint8_t *pGryImg, // 原始灰度图像
                     int width,       // 图像的宽度 
                     int height,      // 图像的高度
                     int *pSumImg     // 计算得到的积分图
                    )
{
    uint8_t *pGry;
    int *pRes;
    int x, y;
    int sumCol[4096]; // 约定图像宽度不大于4096

    memset(sumCol, 0, sizeof(int) * width);
    for (y = 0, pGry = pGryImg, pRes = pSumImg; y < height; y++)
    {
        // 最左侧像素的特别处理
        sumCol[0] += *(pGry++);
        *(pRes++) = sumCol[0];
        // 正常处理
        for (x = 1; x < width; x++)
        {
            sumCol[x] += *(pGry++);       // 更新列积分
            *(pRes++) = *(pRes - 1) + sumCol[x];
        }
    }
    return;
}

其实还有一种实现方式来实现

基于SSE的积分图实现

图像数据是一种非常特别的数据，灰度值一般是8bits的，颜色分量也是8bits的，但是现在计算机的数据总线宽一般都是64位，这就是说在访问一个像素时，数据总线上传输的64个bit中只有8个bit是有效的，另外的都白白浪费了。为此，ARM为处理图像等多媒体数据设计了NEON指令集；Intel为处理图像等多媒体数据设计了MMX和SSE指令集。MMX、SSE、AVX都是利用CPU内部的寄存器进行计算，MMX寄存器的宽度是64位，SSE寄存器的宽度为128位，AVX寄存器的宽度为256位。

采用MMX或者SSE实现C/C++程序优化有两种方式：一种是嵌入式汇编的方式，需要将汇编代码嵌入到C/C++语句中，但这样的程序可读性很差；另外一种是内建函数的方式，可以像其他函数一样直接调用，程序可读性较好。在C/C++编程中，通常使用第二种方式，这两种方式的执行效率是相等的。

​ 内建函数是按照约定语法规则的函数，如果各家编译器支持该语法规则，则必须为使用者提供其函数，这些函数包含在编译器的运行库中，程序员不必单独书写代码，只需要调用这些函数即可，它们的实现由编译器厂商完成，比如在VC++程序设计中，只需包含<nmmintrin.h>即可。

void RmwDoSumGryImg_SSE(uint8_t *pGryImg,  //原始灰度图像
                         int width,        //图像的宽度，必须是4的倍数
                         int height,       //图像的高度
                         int *pSumImg      //计算得到的积分图
                        )
{
    int sumCol[4096]; //约定图像宽度不大于4096
    __m128i *pSumSSE, A;
    uint8_t *pGry;
    int *pRes;
    int x, y;

    memset(sumCol, 0, sizeof(int) * width);
    for (y = 0, pGry = pGryImg, pRes = pSumImg; y < height; y++)
    {
        // 0:需要特别处理
        sumCol[0] += *(pGry++);
        *(pRes++) = sumCol[0];
        // 1
        sumCol[1] += *(pGry++);
        *(pRes++) = *(pRes - 1) + sumCol[1];
        // 2
        sumCol[2] += *(pGry++);
        *(pRes++) = *(pRes - 1) + sumCol[2];
        // 3
        sumCol[3] += *(pGry++);
        *(pRes++) = *(pRes - 1) + sumCol[3];
        // [4...width-1]
        for (x = 4, pSumSSE = (__m128i *)(sumCol + 4); x < width; x += 4, pGry += 4)
        {
            // 把变量的低32位(有4个8位整数组成)转换成32位的整数
            A = _mm_cvtepu8_epi32(_mm_loadl_epi64((__m128i *)pGry));
            // 4个32位的整数相加
            *(pSumSSE++) = _mm_add_epi32(*pSumSSE, A);
            // 递推
            *(pRes++) = *(pRes - 1) + sumCol[x + 0];
            *(pRes++) = *(pRes - 1) + sumCol[x + 1];
            *(pRes++) = *(pRes - 1) + sumCol[x + 2];
            *(pRes++) = *(pRes - 1) + sumCol[x + 3];
        }
    }
    return;
}

经实际测试，在某款CPU上，后者比前者的速度约提高了2.8倍**，这是因为后者使用了SSE内建函数，**能够同时更新4个sumCol的值。

---

那么滤波的函数也是呼之欲出了：

基于积分图的快速均值滤波

void RmwAvrFilterBySumImg(int *pSumImg, // 计算得到的积分图
                           int width, int height, // 图像的宽度和高度
                           int M, int N, // 滤波邻域：M列N行
                           uint8_t *pResImg // 结果图像
                          )
{
    // 没有对边界上邻域不完整的像素进行处理，可以采用变窗口的策略
    int *pY1, *pY2;
    uint8_t *pRes;
    int halfx, halfy;
    int y, x1, x2;
    int sum, c;

    // step.1------------初始化--------------------------//
    M = M / 2 * 2 + 1; // 奇数化
    N = N / 2 * 2 + 1; // 奇数化
    halfx = M / 2;      // 滤波器的半径x
    halfy = N / 2;      // 滤波器的半径y
    c = (1 << 23) / (M * N); // 乘法因子
    // step.2------------滤波----------------------------//
    for (y = halfy + 1, pRes = pResImg + y * width, pY1 = pSumImg, pY2 = pSumImg + N * width;
         y < height - halfy;
         y++, pY1 += width, pY2 += width)
    {
        pRes += halfx + 1; // 跳过左侧
        for (x1 = 0, x2 = M; x2 < width; x1++, x2++) // 可以简化如此，但不太容易读
        {
            sum = *(pY2 + x2) - *(pY2 + x1) - *(pY1 + x2) + *(pY1 + x1);
            *(pRes++) = (uint8_t)((sum * c) >> 23); // 用整数乘法和移位代替除法
        }
        pRes += halfx; // 跳过右侧
    }
    // step.3------------返回----------------------------//
    return;
}

可以看到也是可以实现的

这里要注意一个问题，就是申请积分图内存时要记得类型转换，避免指针越界访问出现的问题

int* pSumImg = (int*)malloc(grayScaleWidth * grayScaleHeight * sizeof(int));//申请积分图的内存

中值滤波

中值滤波的由来

为了解决图像变模糊的问题

• 如何解决这个问题呢？首先来分析一下产生模糊的原因是什么。模糊的原因就是没有区分像素的类别，而将不同类别像素的灰度值进行了平均，比如没有区分是牛还是老鼠，而将牛的体重和老鼠的体重进行了累加，将其均值作为牛（老鼠）的均值滤波后的体重。那么，怎么区分是牛还是老鼠呢？难道在均值滤波以前还得先做个分类？
• 体重在不同类别之间求平均是不对的（会生成不存在的物种），因此可以先在邻域内进行判断，如果牛的数量多就取牛的体重，如果老鼠的数量多就取老鼠的体重，是从里面挑出一个体重，而不是求均值去生成一个新的体重，这样就把分类的问题转换成了谁多的问题。那么在邻域内是牛多还是老鼠多呢？显然，将它们的体重进行从小到大排序，则谁多就谁出现在排序的最中间位置。

将n(n为奇数)个数据按其值d_i进行从大到小或者从小到大排列后得到一个有序序列d_0 d_1…d_(n-1)，则d_⌊n/2⌋ 称为中值。例如：有序序列10,11,12,13,14,15,16,17,18的n=9，有⌊9/2⌋=4，则中值为d_4，即14。

也就是说：

中值滤波就是以当前像素为中心取一个邻域，用该区域的所有像素灰度值的中值作为该像素滤波后的灰度值

所以最后的图像就是

可以看到确实变得很平滑

中值滤波的特点

中值滤波在实际应用时，其邻域大小的选择要保证该邻域内最多有2类目标。至于受到噪声干扰的像素，其灰度值在邻域内要么最大、要么最小，一般不会是中值，所以会被滤除。当邻域不大时，可以认为邻域内最多有2类目标。

• 中值一定是邻域内某个像素的灰度值而不是某几个像素灰度的生成值（相比于它，均值就像是个伪造的值），所以中值滤波不会使图像变模糊，这是它的优点；
• 但是，中值滤波需要排序，而排序的复杂度远比相加求和要大得多，所以中值滤波的速度要比均值滤波慢很多，这是它的缺点。（其实，在图像处理中，排除编程技巧等的因素后，几乎可以说算法越复杂，则执行速度越慢，则处理效果越好）。

中值滤波的快速实现

利用直方图数据本身的有序特性，使得求图像的中值非常简单，大大减少了比较次数。利用直方图求中值最多进行256次比较，因此利用直方图来快速得到中值。

void GetMedianGry(int *histogram, int N, int *medGry)
{
	int g;
	int num;

	// step.1-------------求灰度中值------------------------//
	num = 0;
	for (g = 0; g<256; g++)
	{
		num += histogram[g];
		if (2*num>N) break;  //num>N/2
	}
	*medGry = g;
	// step.2-------------结束------------------------------//
	return;
}

好的，那么接下来就可以实现中值滤波了

double RmwMedianFilter(uint8_t *pGryImg, int width, int height,int M, int N,uint8_t *pResImg) {
    uint8_t *pCur, *pRes;
    int halfx, halfy, x, y, i, j, y1, y2;
    int histogram[256];
    int wSize, j1, j2;
    int num, med, v;
    int dbgCmpTimes = 0; // 搜索中值所需比较次数的调试

    M = M / 2 * 2 + 1; // 奇数化
    N = N / 2 * 2 + 1; // 奇数化
    halfx = M / 2;      // x半径
    halfy = N / 2;      // y半径
    wSize = (halfx * 2 + 1) * (halfy * 2 + 1); // 邻域内像素总个数

    for (y = halfy, pRes = pResImg + y * width; y < height - halfy; y++) {
        // step.1----初始化直方图
        y1 = y - halfy;
        y2 = y + halfy;
        memset(histogram, 0, sizeof(int) * 256);

        for (i = y1, pCur = pGryImg + i * width; i <= y2; i++, pCur += width) {
            for (j = 0; j < halfx * 2 + 1; j++) {
                histogram[*(pCur + j)]++;
            }
        }

        // step.2-----初始化中值
        num = 0; // 记录着灰度值从0到中值的个数
        for (i = 0; i < 256; i++) {
            num += histogram[i];
            if (num * 2 > wSize) {
                med = i;
                break;
            }
        }

        // 滤波
        pRes += halfx; // 没有处理图像左边界侧的像素
        for (x = halfx; x < width - halfx; x++) {
            // 赋值
            *(pRes++) = med;

            // step.3-----直方图递推: 减去当前邻域最左边的一列,添加邻域右侧的一个新列
            j1 = x - halfx;     // 最左边列
            j2 = x + halfx + 1; // 右边的新列

            for (i = y1, pCur = pGryImg + i * width; i <= y2; i++, pCur += width) {
                // 减去最左边列
                v = *(pCur + j1);
                histogram[v]--;  // 更新直方图
                if (v <= med) num--; // 更新num

                // 添加右边的新列
                v = *(pCur + j2);
                histogram[v]++; // 更新直方图
                if (v <= med) num++; // 更新num
            }

            // step.4-----更新中值
            if (num * 2 < wSize) { // 到上次中值med的个数不够了,则med要变大
                for (med = med + 1; med < 256; med++) {
                    dbgCmpTimes += 2; // 总的比较次数,调试用
                    num += histogram[med];
                    if (num * 2 > wSize) break;
                }
                dbgCmpTimes += 1; // 总的比较次数,调试用
            } else { // 到上次中值med的个数多了,则med要变小
                while ((num - histogram[med]) * 2 > wSize) { // 若减去后,仍变小
                    dbgCmpTimes++; // 总的比较次数,调试用
                    num -= histogram[med];
                    med--;
                }
                dbgCmpTimes += 2; // 总的比较次数,调试用
            }
        }
        pRes += halfx; // 没有处理图像右边界侧的像素
    }
    // 返回搜索中值需要的平均比较次数
    return dbgCmpTimes * 1.0 / ((width - halfx * 2) * (height - halfy * 2));
}

我们来看效果

极值滤波

如果事先能够知道噪声像素的灰度值是比周围像素的灰度值大，那么就可以使用邻域内的最小值来替代当前像素的灰度值，称为最小值滤波（Minimum Filter）。反之，称为最大值滤波（Maximum Filter）。

这个和中值滤波差不多，大家可以自己修改代码试试

---

高斯滤波

什么是高斯滤波？

对3个像素的灰度值求平均值即式

其对应的模板如下图所示。

下图虽然是5邻域均值滤波的模板，但是该模板中最左侧和最右侧的值为“0”，相当于没有起作用，它实际上也是3邻域的均值滤波。

可以这样认为，因为该邻域内的最左侧像素、最右侧像素到中心像素的距离远了，所以把它们对G(x)的贡献度设为0，令它们对均值没有贡献。

简单地将权重设置成要么“1”要么“0”，没有柔和地体现**“近处贡献度大、远处贡献度小”**的特性。那么，用什么函数能够更好地体现贡献度的变换呢？均值为0的高斯函数是一个常用的权重函数，如图3-32所示，其表达式如下

使用高斯函数加权均值滤波时，邻域的半径一般定为3σ，半径超过3σ时，对滤波结果的影响不大

习惯上把使用高斯函数作为权重函数的均值滤波，称为高斯滤波（Gaussian Average Filter）或者高斯平滑。同其他均值滤波一样，高斯滤波也会导致结果图像变模糊，所以有时也称为高斯模糊

高斯滤波的特点是：当σ越大时，高斯函数也越平缓；当σ越大时，邻域的尺寸也越大；邻域尺寸越大，就有更多的像素参加了滤波，图像模糊得也就越厉害

为了快速计算，二维的高斯滤波可以被分解为两个一维的高斯滤波，比如先对图像按行进行一维的高斯滤波，再对得到的结果图像按列进行一维的高斯滤波。另外，为了提高Cache的命中率，从而提高速度，还可以把行滤波得到的结果图像进行90度转置，把按列进行的一维高斯滤波变为按行的一维高斯滤波，最后再把得到的最终结果图像转置回来

二值图像滤波与数学形态学滤波

基于均值滤波的二值图像极值与中值滤波

二值图像：

• 二值图像只有2种灰度值，比如分别是“0”或者“255”。
• 二值图像的灰度最小值的和灰度最大值也不用求取；最小值为0，最大值为255。
• 二值图像可以进行均值滤波，但是不能在结果图像中给像素赋值为非“0”且非“255”的灰度值。
• 二值图像可以进行中值滤波，但是肯定不用进行排序。在一个邻域内，若“0”的个数多，中值就是0，此时邻域内灰度均值小于128；反之，中值就是255，此时邻域内灰度均值大于等于128。

所以，我们就能得出以下的概念

（1）当邻域内灰度均值大于等于255时，赋值255；否则，赋值0；这就是最小值滤波。

（2）当邻域内灰度均值大于0时，赋值255；否则，赋值0；这就是最大值滤波。

（3）当邻域内灰度均值大于或者等于128时，赋值255；否则，赋值0；这就是中值滤波。

上代码

void RmwBinImgFilter(uint8_t *pBinImg,  //原始二值图像
                     int width, int height,//图像的宽度和高度
                     int M, int N, //滤波邻域：M列N行
                     double threshold, //灰度阈值,大于等于该值时结果赋255
                     uint8_t *pResImg //结果图像
                    )
{
    // 没有对边界上邻域不完整的像素进行处理，可以采用变窗口的策略
    uint8_t *pAdd, *pDel, *pRes;
    int halfx, halfy;
    int x, y, sum, sumThreshold;
    int sumCol[4096]; //约定图像宽度不大于4096

    // step.1------------初始化--------------------------//
    M = M / 2 * 2 + 1; //奇数化
    N = N / 2 * 2 + 1; //奇数化
    halfx = M / 2; //滤波器的x半径
    halfy = N / 2; //滤波器的y半径
    sumThreshold = max(1, (int)(threshold * M * N)); //转换成邻域内灰度值之和的阈值
    memset(sumCol, 0, sizeof(int) * width);
    for (y = 0, pAdd = pBinImg; y < N; y++)
    {
        for (x = 0; x < width; x++)
            sumCol[x] += *(pAdd++);
    }
    // step.2------------滤波----------------------------//
    for (y = halfy, pRes = pResImg + y * width, pDel = pBinImg; y < height - halfy; y++)
    {
        //初值
        for (sum = 0, x = 0; x < M; x++)
            sum += sumCol[x];
        //滤波
        pRes += halfx; //跳过左侧
        for (x = halfx; x < width - halfx; x++)
        {
            //求灰度均值
            /*if (sum>=sumThreshold)
            {
                *(pRes++) = 255;
            }
            else  *(pRes++) = 0;*/
            *(pRes++) = (sum >= sumThreshold) * 255; //请理解这个表达式的含义
            //换列,更新灰度和
            sum -= sumCol[x - halfx];     //减左边列
            sum += sumCol[x + halfx + 1]; //加右边列
        }
        pRes += halfx; //跳过右侧
        //换行,更新sumCol
        for (x = 0; x < width; x++)
        {
            sumCol[x] -= *(pDel++); //减上一行
            sumCol[x] += *(pAdd++); //加下一行
        }
    }
    // step.3------------返回----------------------------//
    return;
}

可以看到效果也还是很好的

二值图像的数学形态学滤波

对于二值图像的滤波，还有一种称之为“数学形态学”（Mathematical Morphology）的滤波方法，简称形态学滤波。下面不具体讲述数学形态学的内涵，只是讲一下形态学滤波的基本运算。

​ 二值图像一般都是原始图像经过某种处理后得到的，它的像素具有明确的含义，比如黑色像素 (以灰度值“0”表示)代表背景，白色像素(以灰度值“1”表示)代表目标，而目标往往具有某种几何形状，即具有某种形态。形态学滤波就是以目标形状为出发点来进行滤波，并发展出了膨胀、腐蚀、开运算、闭运算等方法。

膨胀运算：

和膨胀对应的还有腐蚀运算

形态学滤波有一套较为晦涩的符号体系，不再细讲。可以认为形态学滤波就是使用布尔代数和集合运算的二值图像滤波，形态学滤波能做到的，常规的滤波方法也能做到。在本质上，空间域滤波就是邻域的选定和邻域内计算函数的选定，形态学滤波和常规的滤波方法是一致的。

条件滤波

超限平滑

所谓超限平滑，就是均值滤波得到的值u(x,y)和原值g(x,y)相比，超过了一定的程度C才使用新值，否则保持原值

此种情况使用超限平滑的合理性在于，在实际拍摄得到的图像中，目标和背景之间至少存在着1个像素的过渡区，即目标和背景的边界像素灰度值是由目标的灰度逐渐变化到背景的灰度的。边界上的像素的灰度值g(x,y)与其滤波均值u(x,y)是接近的，而噪声的灰度值则与均值有较大的差异。

K个邻点平均法

所谓K个邻点平均法就是不使用邻域的全部像素，而是只使用其中的K个像素求均值；假设邻域中有A个像素，则K<A。

在A个像素中，到底选取那K个呢？有一个基本原则就是选取与当前像素的灰度值最接近的那K个像素。

此种情况使用K个邻点平均法的合理性在于，因为噪声的灰度值跟正常像素的灰度值不接近，所以噪声的灰度值能被其周围（即邻域）的灰度值修改掉；而目标像素或者背景像素的邻域内，肯定有多个同类像素，设为K个，用这K个同类像素求均值，肯定不会产生模糊。

在邻域内选择与当前像素的灰度值最接近的K个灰度值，肯定需要比较操作，因此这个滤波器可以看成是既有中值滤波的排序操作，又有均值滤波的求和操作。

多邻域枚举法均值滤波

在邻域内选择K个像素是不太容易的。不妨认为，在一个邻域内，属于同一类别K个像素的分布形式一共有N种，如果把每一种都作为一个模板，就是N个模板。这样每个模板都得到一个均值，总共得到了N个均值，从而从这N个均值中选择一个最佳模板的均值作为滤波结果。

如何选择最佳模板呢？

若模板内包含中心像素的K个像素都来自相同类，则该模板肯定是最佳的。

如何评价模板内的像素来自相同类呢？

模板内像素灰度值的均方差越小，则说明该模板越具有灰度一致性，即越有可能不包含边缘和不包含多类目标，因此该模板就越佳。

上述做法是同时采用N种邻域，因此可将其称为多邻域枚举法均值滤波。述。

计算量一定巨大

源码

IDP.h(函数声明)

#pragma once

#include <stdio.h>
#include <stdint.h>
#include <stdlib.h>
#include <math.h>
#include <nmmintrin.h>
 
uint8_t* readGrayScaleBMP(const char* filename, int* width, int* height);//读取8位灰度图片
void saveGrayScaleBMP(const char* filename, const uint8_t* imageData, int width, int height);// 将8位灰度图像数据保存为BMP文件
uint8_t* readColorBMP(const char* filename, int* width, int* height);//读取24位彩色图像的BMP文件
void saveColorBMP(const char* filename, const uint8_t* imageData, int width, int height);//将24位彩色图像数据保存为BMP文件
void LinearStretchDemo(uint8_t* pGryImg, int width, int height, double k, double b);//灰度线性拉伸
void GetHistogram(uint8_t* pImg, int width, int height, int* histogram);//统计图像灰度值
void GetBrightContrast(int* histogram, double* bright, double* contrast);//亮度和对比度
void RmwHistogramEqualize(uint8_t* pGryImg, int width, int height);//直方图均衡化
void RmwLogTransform(uint8_t* pGryImg, int width, int height);//对数变换
void RmwAvrFilterBySumCol(uint8_t* pGryImg, int width, int height, int M, int N, uint8_t* pResImg);//基于列积分的快速均值滤波
void RmwDoSumGryImg(uint8_t* pGryImg, int width, int height, int* pSumImg);//基于列积分的积分图实现
void RmwDoSumGryImg_SSE(uint8_t* pGryImg, int width, int height, int* pSumImg);//基于SSE的积分图实现
void RmwAvrFilterBySumImg(int* pSumImg, int width, int height, int M, int N, uint8_t* pResImg);//基于积分图的快速均值滤波  
void GetMedianGry(int* histogram, int N, int* medGry);//求灰度值中值
double RmwMedianFilter(uint8_t* pGryImg, int width, int height, int M, int N, uint8_t* pResImg);//中值滤波
void RmwBinImgFilter(uint8_t* pBinImg, int width, int height, int M, int N, double threshold, uint8_t* pResImg);//二值滤波

IDP.c(函数主要实现)

#define _CRT_SECURE_NO_WARNINGS 1

#include "IDP.h"

//读取8位灰度图片
//filename：字符数组的指针，用于指定要保存的图像文件的名称或路径。
//imageData：无符号 8 位整型数据的指针，代表要保存的图像数据。
//width：图像的宽度。
//height：图像的高度。
uint8_t* readGrayScaleBMP(const char* filename, int* width, int* height) 
{
    FILE* file = fopen(filename, "rb");
    if (!file) {
        fprintf(stderr, "Error opening file %s\n", filename);
        return NULL;
    }

    // 读取BMP文件头部信息
    uint8_t bmpHeader[54];
    fread(bmpHeader, 1, 54, file);

    // 从文件头部提取图像宽度和高度信息
    *width = *(int*)&bmpHeader[18];
    *height = *(int*)&bmpHeader[22];

    // 分配存储图像数据的内存
    uint8_t* imageData = (uint8_t*)malloc(*width * *height);
    if (!imageData) {
        fprintf(stderr, "内存分配失败\n");
        fclose(file);
        return NULL;
    }

    // 计算调色板的大小
    int paletteSize = *(int*)&bmpHeader[46];
    if (paletteSize == 0)
        paletteSize = 256;

    // 读取调色板数据
    uint8_t palette[1024];
    fread(palette, 1, paletteSize * 4, file);

    // 读取图像数据
    fseek(file, *(int*)&bmpHeader[10], SEEK_SET);
    fread(imageData, 1, *width * *height, file);

    fclose(file);

    return imageData;
}

// 将8位灰度图像数据保存为BMP文件
//filename：字符数组的指针，用于指定要保存的图像文件的名称或路径。
//imageData：无符号 8 位整型数据的指针，代表要保存的图像数据。
//width：图像的宽度。
//height：图像的高度。
void saveGrayScaleBMP(const char* filename, const uint8_t* imageData, int width, int height) 
{
    FILE* file = fopen(filename, "wb");
    if (!file) {
        fprintf(stderr, "Error creating file %s\n", filename);
        return;
    }

    // BMP文件头部信息
    uint8_t bmpHeader[54] = {
        0x42, 0x4D,             // 文件类型标识 "BM"
        0x36, 0x00, 0x0C, 0x00, // 文件大小（以字节为单位，此处假设图像数据大小不超过4GB）
        0x00, 0x00,             // 保留字段
        0x00, 0x00,             // 保留字段
        0x36, 0x00, 0x00, 0x00, // 位图数据偏移（以字节为单位）
        0x28, 0x00, 0x00, 0x00, // 位图信息头大小（40字节）
        0x00, 0x00, 0x00, 0x00, // 图像宽度
        0x00, 0x00, 0x00, 0x00, // 图像高度
        0x01, 0x00,             // 目标设备的级别（此处为1，不压缩）
        0x08, 0x00,             // 每个像素的位数（8位）
        0x00, 0x00, 0x00, 0x00, // 压缩类型（此处为不压缩）
        0x00, 0x00, 0x00, 0x00, // 图像数据大小（以字节为单位，此处为0，表示不压缩）
        0x00, 0x00, 0x00, 0x00, // 水平分辨率（像素/米，此处为0，表示未知）
        0x00, 0x00, 0x00, 0x00, // 垂直分辨率（像素/米，此处为0，表示未知）
        0x00, 0x00, 0x00, 0x00, // 使用的颜色索引数（0表示使用所有调色板项）
        0x00, 0x00, 0x00, 0x00  // 重要的颜色索引数（0表示所有颜色都重要）
    };

    // 更新BMP文件头部信息中的宽度和高度
    *(int*)&bmpHeader[18] = width;
    *(int*)&bmpHeader[22] = height;

    // 写入BMP文件头部信息
    fwrite(bmpHeader, 1, 54, file);

    // 写入调色板数据
    for (int i = 0; i < 256; i++) {
        fputc(i, file);  // 蓝色分量
        fputc(i, file);  // 绿色分量
        fputc(i, file);  // 红色分量
        fputc(0, file);  // 保留字节
    }

    // 写入图像数据
    fwrite(imageData, 1, width * height, file);

    fclose(file);
}

// 读取24位彩色图像的BMP文件
//filename：字符数组的指针，用于指定要读取的 BMP 格式图像文件的名称或路径。
//width：整型变量的指针，用于存储读取的图像的宽度。
//height：整型变量的指针，用于存储读取的图像的高度。
uint8_t* readColorBMP(const char* filename, int* width, int* height) 
{
    FILE* file = fopen(filename, "rb");
    if (!file) {
        fprintf(stderr, "Error opening file %s\n", filename);
        return NULL;
    }

    // 读取BMP文件头部信息
    uint8_t bmpHeader[54];
    fread(bmpHeader, 1, 54, file);

    // 从文件头部提取图像宽度和高度信息
    *width = *(int*)&bmpHeader[18];
    *height = *(int*)&bmpHeader[22];

    // 分配存储图像数据的内存
    uint8_t* imageData = (uint8_t*)malloc(*width * *height * 3);
    if (!imageData) {
        fprintf(stderr, "Memory allocation failed\n");
        fclose(file);
        return NULL;
    }

    // 读取图像数据
    fseek(file, *(int*)&bmpHeader[10], SEEK_SET);
    fread(imageData, 1, *width * *height * 3, file);

    fclose(file);

    return imageData;
}

//将24位彩色图像数据保存为BMP文件
//filename：字符数组的指针，用于指定要保存的图像文件的名称或路径。
//imageData：无符号 8 位整型数据的指针，代表要保存的图像数据。
//width：图像的宽度。
//height：图像的高度。
void saveColorBMP(const char* filename, const uint8_t* imageData, int width, int height) 
{
    FILE* file = fopen(filename, "wb");
    if (!file) {
        fprintf(stderr, "Error creating file %s\n", filename);
        return;
    }

    // BMP文件头部信息
    uint8_t bmpHeader[54] = {
        0x42, 0x4D,             // 文件类型标识 "BM"
        0x00, 0x00, 0x00, 0x00, // 文件大小（占位，稍后计算）
        0x00, 0x00,             // 保留字段
        0x00, 0x00,             // 保留字段
        0x36, 0x00, 0x00, 0x00, // 位图数据偏移（以字节为单位）
        0x28, 0x00, 0x00, 0x00, // 位图信息头大小（40字节）
        0x00, 0x00, 0x00, 0x00, // 图像宽度
        0x00, 0x00, 0x00, 0x00, // 图像高度
        0x01, 0x00,             // 目标设备的级别（此处为1，不压缩）
        0x18, 0x00,             // 每个像素的位数（24位）
        0x00, 0x00, 0x00, 0x00, // 压缩类型（此处为不压缩）
        0x00, 0x00, 0x00, 0x00, // 图像数据大小（占位，稍后计算）
        0x00, 0x00, 0x00, 0x00, // 水平分辨率（像素/米，此处为0，表示未知）
        0x00, 0x00, 0x00, 0x00, // 垂直分辨率（像素/米，此处为0，表示未知）
        0x00, 0x00, 0x00, 0x00, // 使用的颜色索引数（0表示使用所有调色板项）
        0x00, 0x00, 0x00, 0x00  // 重要的颜色索引数（0表示所有颜色都重要）
    };

    // 更新BMP文件头部信息中的宽度和高度
    *(int*)&bmpHeader[18] = width;
    *(int*)&bmpHeader[22] = height;

    // 计算图像数据大小
    uint32_t imageDataSize = width * height * 3 + 54; // 加上文件头部大小
    bmpHeader[2] = (uint8_t)(imageDataSize & 0xFF);
    bmpHeader[3] = (uint8_t)((imageDataSize >> 8) & 0xFF);
    bmpHeader[4] = (uint8_t)((imageDataSize >> 16) & 0xFF);
    bmpHeader[5] = (uint8_t)((imageDataSize >> 24) & 0xFF);

    // 写入BMP文件头部信息
    fwrite(bmpHeader, 1, 54, file);

    // 写入图像数据
    fwrite(imageData, width * height * 3, 1, file);

    fclose(file);
}

//灰度线性拉伸
//pGryImg：灰度图像数据的指针。
//width：图像的宽度。
//height：图像的高度。
//k：线性拉伸的斜率。它控制着拉伸的速率或程度。当(k) 大于 1 时，图像的对比度增加；当(k) 小于 1 时，对比度降低。
//b：线性拉伸的偏移。它控制着拉伸后灰度值的起始位置。当(b) 大于 0 时，图像的整体亮度增加；当(b) 小于 0 时，整体亮度减小。
void LinearStretchDemo(uint8_t* pGryImg, int width, int height, double k, double b)
{
    uint8_t* pCur, * pEnd;
    int LUT[256];    //因为只有[0,255]共256个灰度值

    //step1. 生成查找表
    for (int g = 0; g < 256; g++)
    {
        LUT[g] = max(0, min(255, k * g + b));
    }

    //step2. 进行变换
    for (pCur = pGryImg, pEnd = pGryImg + width * height; pCur < pEnd; pCur++)
    {
        *pCur = LUT[*pCur];
    }
    //step3. 结束
    return;
}

//统计图像灰度值
//pImg：灰度图像数据的指针。
//width：图像的宽度。
//height：图像的高度。
//* histogram：数组首元素地址，需要一个能储存256个变量的整型数组
void GetHistogram(uint8_t* pImg, int width, int height, int* histogram)
{
    uint8_t* pCur;
    uint8_t* pEnd = pImg + width * height;

    // 初始化直方图数组
    memset(histogram, 0, sizeof(int) * 256);

    // 直方图统计
    for (pCur = pImg; pCur < pEnd;)
    {
        histogram[*pCur]++;
        pCur++;
    }

    // 函数结束
    return;
}

//亮度和对比度
//储存histogram灰度直方图的指针
//接收亮度的变量地址
//接收对比度的变量地址
void GetBrightContrast(int* histogram, double* bright, double* contrast)
{
    int g;
    double sum, num; //书上说图像很亮时，int有可能会溢出，所以我这里直接用double
    double fsum;

    //step.1 求亮度
    for (sum = num = 0, g = 0; g < 256; g++)
    {
        sum += histogram[g] * g;
        num += histogram[g];
    }
    *bright = sum * 1.0 / num;

    //step.2 求对比度
    for (fsum = 0.0, g = 0; g < 256; g++)
    {
        fsum += histogram[g] * (g - *bright) * (g - *bright);
    }
    *contrast = sqrt(fsum / (num - 1)); //即Std Dev

    //step.3 结束
    return;
}

//pGryImg：灰度图像数据的指针。
//width：图像的宽度。
//height：图像的高度。
void RmwHistogramEqualize(uint8_t* pGryImg, int width, int height)
{
    uint8_t* pCur, * pEnd = pGryImg + width * height; // 指针变量，指向当前像素和图像末尾
    int histogram[256], LUT[256], A, g; // 直方图数组、查找表数组、累积直方图、灰度级

    // step.1-------------求直方图--------------------------//
    memset(histogram, 0, sizeof(int) * 256); // 初始化直方图数组为0
    for (pCur = pGryImg; pCur < pEnd;)
        histogram[*(pCur++)]++; // 统计每个灰度级出现的频率

    // step.2-------------求LUT[g]-------------------------//
    A = histogram[0]; // 初始化累积直方图的值为第一个灰度级的频率
    LUT[0] = 255 * A / (width * height); // 计算第一个灰度级对应的均衡化后的灰度值
    for (g = 1; g < 256; g++) {
        A += histogram[g]; // 更新累积直方图的值
        LUT[g] = 255 * A / (width * height); // 计算当前灰度级对应的均衡化后的灰度值
    }

    // step.3-------------查表------------------------------//
    for (pCur = pGryImg; pCur < pEnd;)
        *(pCur++) = LUT[*pCur]; // 使用查找表对每个像素进行灰度映射

    // step.4-------------结束------------------------------//
    return;
}

//对数变换
//pGryImg：灰度图像数据的指针。
//width：图像的宽度。
//height：图像的高度。
void RmwLogTransform(uint8_t* pGryImg, int width, int height)
{
    uint8_t* pCur, * pEnd = pGryImg + width * height; // 指向灰度图像数据的当前指针和结束指针
    int histogram[256], LUT[256], gmax, g; // 声明直方图数组、查找表数组、最大灰度值、当前灰度值
    double c; // 声明常数c

    // step.1-------------求直方图--------------------------//
    memset(histogram, 0, sizeof(int) * 256); // 初始化直方图数组为0
    for (pCur = pGryImg; pCur < pEnd;)
        histogram[*(pCur++)]++; // 遍历图像数据，统计每个灰度级的像素数量

    // step.2-------------最大值---------------------------//
    for (gmax = 255; gmax >= 0; gmax++)
        if (histogram[gmax]) break; // 从最大灰度级开始向低灰度级搜索，找到第一个非零灰度级，即最大灰度值

    // step.3-------------求LUT[g]-------------------------//
    c = 255.0 / log(1 + gmax); // 计算常数c
    for (g = 0; g < 256; g++)
    {
        LUT[g] = (int)(c * log(1 + g)); // 根据对数变换公式计算查找表中每个灰度级的映射值
    }

    // step.4-------------查表------------------------------//
    for (pCur = pGryImg; pCur < pEnd;)
        *(pCur++) = LUT[*pCur]; // 使用查找表将图像数据进行对数变换

    // step.5-------------结束------------------------------//
    return; // 函数结束
}

//基于列积分的快速均值滤波
//原始灰度图像
//图像的宽度和高度
//滤波邻域：M列N行
//结果图像
void RmwAvrFilterBySumCol(uint8_t* pGryImg,int width, int height,int M, int N,uint8_t* pResImg) 
{
    uint8_t* pAdd, * pDel, * pRes;
    int halfx, halfy;
    int x, y;
    int sum, c;
    int sumCol[4096]; // 约定图像宽度不大于4096

    // step.1------------初始化--------------------------//
    M = M / 2 * 2 + 1; // 奇数化
    N = N / 2 * 2 + 1; // 奇数化
    halfx = M / 2; // 滤波器的半径x
    halfy = N / 2; // 滤波器的半径y
    c = (1 << 23) / (M * N); // 乘法因子
    memset(sumCol, 0, sizeof(int) * width);
    for (y = 0, pAdd = pGryImg; y < N; y++) {
        for (x = 0; x < width; x++) sumCol[x] += *(pAdd++);
    }
    // step.2------------滤波----------------------------//
    for (y = halfy, pRes = pResImg + y * width, pDel = pGryImg; y < height - halfy; y++) {
        // 初值
        for (sum = 0, x = 0; x < M; x++) sum += sumCol[x];
        // 滤波
        pRes += halfx; // 跳过左侧
        for (x = halfx; x < width - halfx; x++) {
            // 求灰度均值
            // *(pRes++)=sum/(N*M);
            *(pRes++) = (sum * c) >> 23; // 用整数乘法和移位代替除法
            // 换列,更新灰度和
            sum -= sumCol[x - halfx]; // 减左边列
            sum += sumCol[x + halfx + 1]; // 加右边列
        }
        pRes += halfx; // 跳过右侧
        // 换行,更新sumCol
        for (x = 0; x < width; x++) {
            sumCol[x] -= *(pDel++); // 减上一行
            sumCol[x] += *(pAdd++); // 加下一行
        }
    }
    // step.3------------返回----------------------------//
    return;
}

//基于列积分的积分图实现
//pGryImg, // 原始灰度图像
//width,       // 图像的宽度 
//height,      // 图像的高度
//pSumImg     // 计算得到的积分图
void RmwDoSumGryImg(uint8_t* pGryImg,int width,int height, int* pSumImg)
{
    uint8_t* pGry;
    int* pRes;
    int x, y;
    int sumCol[4096]; // 约定图像宽度不大于4096

    memset(sumCol, 0, sizeof(int) * width);
    for (y = 0, pGry = pGryImg, pRes = pSumImg; y < height; y++)
    {
        // 最左侧像素的特别处理
        sumCol[0] += *(pGry++);
        *(pRes++) = sumCol[0];
        // 正常处理
        for (x = 1; x < width; x++)
        {
            sumCol[x] += *(pGry++);       // 更新列积分
            int temp = *(pRes - 1);
            *(pRes++) = temp + sumCol[x];
        }
    }
    return;
}

//基于SSE的积分图实现
//pGryImg原始灰度图像
//width图像的宽度，必须是4的倍数
//height图像的高度
//pSumImg计算得到的积分图
void RmwDoSumGryImg_SSE(uint8_t* pGryImg,int width,int height,int* pSumImg)
{
    int sumCol[4096]; //约定图像宽度不大于4096
    __m128i* pSumSSE, A;
    uint8_t* pGry;
    int* pRes;
    int x, y;

    memset(sumCol, 0, sizeof(int) * width);
    for (y = 0, pGry = pGryImg, pRes = pSumImg; y < height; y++)
    {
        // 0:需要特别处理
        sumCol[0] += *(pGry++);
        *(pRes++) = sumCol[0];
        // 1
        sumCol[1] += *(pGry++);
        *(pRes++) = *(pRes - 1) + sumCol[1];
        // 2
        sumCol[2] += *(pGry++);
        *(pRes++) = *(pRes - 1) + sumCol[2];
        // 3
        sumCol[3] += *(pGry++);
        *(pRes++) = *(pRes - 1) + sumCol[3];
        // [4...width-1]
        for (x = 4, pSumSSE = (__m128i*)(sumCol + 4); x < width; x += 4, pGry += 4)
        {
            // 把变量的低32位(有4个8位整数组成)转换成32位的整数
            A = _mm_cvtepu8_epi32(_mm_loadl_epi64((__m128i*)pGry));
            // 4个32位的整数相加
            *(pSumSSE++) = _mm_add_epi32(*pSumSSE, A);
            // 递推
            *(pRes++) = *(pRes - 1) + sumCol[x + 0];
            *(pRes++) = *(pRes - 1) + sumCol[x + 1];
            *(pRes++) = *(pRes - 1) + sumCol[x + 2];
            *(pRes++) = *(pRes - 1) + sumCol[x + 3];
        }
    }
    return;
}

//基于积分图的快速均值滤波
//pSumImg计算得到的积分图
//width,height,图像的宽度和高度
//M, N,滤波邻域：M列N行
//pResImg 结果图像
void RmwAvrFilterBySumImg(int* pSumImg,int width, int height,int M, int N,uint8_t* pResImg)
{
    // 没有对边界上邻域不完整的像素进行处理，可以采用变窗口的策略
    int* pY1, * pY2;
    uint8_t* pRes;
    int halfx, halfy;
    int y, x1, x2;
    int sum, c;

    // step.1------------初始化--------------------------//
    M = M / 2 * 2 + 1; // 奇数化
    N = N / 2 * 2 + 1; // 奇数化
    halfx = M / 2;      // 滤波器的半径x
    halfy = N / 2;      // 滤波器的半径y
    c = (1 << 23) / (M * N); // 乘法因子
    // step.2------------滤波----------------------------//
    for (y = halfy + 1, pRes = pResImg + y * width, pY1 = pSumImg, pY2 = pSumImg + N * width;
        y < height - halfy;
        y++, pY1 += width, pY2 += width)
    {
        pRes += halfx + 1; // 跳过左侧
        for (x1 = 0, x2 = M; x2 < width; x1++, x2++) // 可以简化如此，但不太容易读
        {
            sum = *(pY2 + x2) - *(pY2 + x1) - *(pY1 + x2) + *(pY1 + x1);
            *(pRes++) = (uint8_t)((sum * c) >> 23); // 用整数乘法和移位代替除法
        }
        pRes += halfx; // 跳过右侧
    }
    // step.3------------返回----------------------------//
    return;
}

void GetMedianGry(int* histogram, int N, int* medGry)
{
    int g;
    int num;

    // step.1-------------求灰度中值------------------------//
    num = 0;
    for (g = 0; g < 256; g++)
    {
        num += histogram[g];
        if (2 * num > N) break;  //num>N/2
    }
    *medGry = g;
    // step.2-------------结束------------------------------//
    return;
}

//中值滤波
//pGryImg：指向待处理灰度图像数据的指针。
//width、height：表示图像的宽度和高度。
//M、N：分别表示中值滤波器的水平和垂直邻域大小（以像素为单位）。
//pResImg：指向存储结果图像数据的指针。
double RmwMedianFilter(uint8_t* pGryImg, int width, int height, int M, int N, uint8_t* pResImg) 
{
    uint8_t* pCur, * pRes;
    int halfx, halfy, x, y, i, j, y1, y2;
    int histogram[256];
    int wSize, j1, j2;
    int num, med, v;
    int dbgCmpTimes = 0; // 搜索中值所需比较次数的调试

    M = M / 2 * 2 + 1; // 奇数化
    N = N / 2 * 2 + 1; // 奇数化
    halfx = M / 2;      // x半径
    halfy = N / 2;      // y半径
    wSize = (halfx * 2 + 1) * (halfy * 2 + 1); // 邻域内像素总个数

    for (y = halfy, pRes = pResImg + y * width; y < height - halfy; y++) {
        // step.1----初始化直方图
        y1 = y - halfy;
        y2 = y + halfy;
        memset(histogram, 0, sizeof(int) * 256);

        for (i = y1, pCur = pGryImg + i * width; i <= y2; i++, pCur += width) {
            for (j = 0; j < halfx * 2 + 1; j++) {
                histogram[*(pCur + j)]++;
            }
        }

        // step.2-----初始化中值
        num = 0; // 记录着灰度值从0到中值的个数
        for (i = 0; i < 256; i++) {
            num += histogram[i];
            if (num * 2 > wSize) {
                med = i;
                break;
            }
        }

        // 滤波
        pRes += halfx; // 没有处理图像左边界侧的像素
        for (x = halfx; x < width - halfx; x++) {
            // 赋值
            *(pRes++) = med;

            // step.3-----直方图递推: 减去当前邻域最左边的一列,添加邻域右侧的一个新列
            j1 = x - halfx;     // 最左边列
            j2 = x + halfx + 1; // 右边的新列

            for (i = y1, pCur = pGryImg + i * width; i <= y2; i++, pCur += width) {
                // 减去最左边列
                v = *(pCur + j1);
                histogram[v]--;  // 更新直方图
                if (v <= med) num--; // 更新num

                // 添加右边的新列
                v = *(pCur + j2);
                histogram[v]++; // 更新直方图
                if (v <= med) num++; // 更新num
            }

            // step.4-----更新中值
            if (num * 2 < wSize) { // 到上次中值med的个数不够了,则med要变大
                for (med = med + 1; med < 256; med++) {
                    dbgCmpTimes += 2; // 总的比较次数,调试用
                    num += histogram[med];
                    if (num * 2 > wSize) break;
                }
                dbgCmpTimes += 1; // 总的比较次数,调试用
            }
            else { // 到上次中值med的个数多了,则med要变小
                while ((num - histogram[med]) * 2 > wSize) { // 若减去后,仍变小
                    dbgCmpTimes++; // 总的比较次数,调试用
                    num -= histogram[med];
                    med--;
                }
                dbgCmpTimes += 2; // 总的比较次数,调试用
            }
        }
        pRes += halfx; // 没有处理图像右边界侧的像素
    }
    // 返回搜索中值需要的平均比较次数
    return dbgCmpTimes * 1.0 / ((width - halfx * 2) * (height - halfy * 2));
}

//二值滤波
//pBinImg,  原始二值图像
// width, height,图像的宽度和高度
// M, N, 滤波邻域：M列N行
// threshold, 灰度阈值,大于等于该值时结果赋255
// pResImg 结果图像
void RmwBinImgFilter(uint8_t* pBinImg,int width, int height,int M, int N,double threshold,uint8_t* pResImg )
{
    // 没有对边界上邻域不完整的像素进行处理，可以采用变窗口的策略
    uint8_t* pAdd, * pDel, * pRes;
    int halfx, halfy;
    int x, y, sum, sumThreshold;
    int sumCol[4096]; //约定图像宽度不大于4096

    // step.1------------初始化--------------------------//
    M = M / 2 * 2 + 1; //奇数化
    N = N / 2 * 2 + 1; //奇数化
    halfx = M / 2; //滤波器的x半径
    halfy = N / 2; //滤波器的y半径
    sumThreshold = max(1, (int)(threshold * M * N)); //转换成邻域内灰度值之和的阈值
    memset(sumCol, 0, sizeof(int) * width);
    for (y = 0, pAdd = pBinImg; y < N; y++)
    {
        for (x = 0; x < width; x++)
            sumCol[x] += *(pAdd++);
    }
    // step.2------------滤波----------------------------//
    for (y = halfy, pRes = pResImg + y * width, pDel = pBinImg; y < height - halfy; y++)
    {
        //初值
        for (sum = 0, x = 0; x < M; x++)
            sum += sumCol[x];
        //滤波
        pRes += halfx; //跳过左侧
        for (x = halfx; x < width - halfx; x++)
        {
            //求灰度均值
            /*if (sum>=sumThreshold)
            {
                *(pRes++) = 255;
            }
            else  *(pRes++) = 0;*/
            *(pRes++) = (sum >= sumThreshold) * 255; //请理解这个表达式的含义
            //换列,更新灰度和
            sum -= sumCol[x - halfx];     //减左边列
            sum += sumCol[x + halfx + 1]; //加右边列
        }
        pRes += halfx; //跳过右侧
        //换行,更新sumCol
        for (x = 0; x < width; x++)
        {
            sumCol[x] -= *(pDel++); //减上一行
            sumCol[x] += *(pAdd++); //加下一行
        }
    }
    // step.3------------返回----------------------------//
    return;
}

感谢您的阅读！！！