噪声
在数字图像处理中，噪声是指在图像中引入的不希望的随机或无意义的信号。它是由于图像采集、传输、存储或处理过程中的各种因素引起的。
噪声会导致图像质量下降，使图像失真或降低细节的清晰度。它通常表现为图像中随机分布的亮度或颜色变化，类似于图像上的颗粒、斑点或像素偏移。
那么噪声如何避免呢？
在图像数字处理和图像分析的算法中，为了保持一定的稳健性，一般需要先对原始图像做一定的滤波处理—图像预处理
图像平滑
数字图像处理中，作用于去掉图像噪声的各种滤波方法统称为图像平滑，那么图像平滑的字面意思就是使一个像素到另一个像素的灰度变化是平滑的。
那么从空间域的观点看，图像平滑实际上就是去去掉突然变大或者变小的点，用一个合适的灰度值代替该值
我们来看一个例子
这样一串数字，他的曲线我们来画一下
空间域的图像滤波围绕邻域计算展开，所以我们图像平滑也将会围绕计算和邻域展开
均值滤波
均值滤波就是去当前像素点为中心取一个邻域，用该邻域的所有像素值的均值作为该像素滤波后的灰度值
就比如
对于每个像素，我们计算它自身以及左右相邻的两个像素的平均值，然后用这个平均值替代原始像素的值。
例如，对于原始数据中的第一个像素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------------初始化---------------------