图像处理算法均值滤波，高斯滤波，双边滤波的实现（处理灰度图像）

本文链接：https://blog.youkuaiyun.com/2301_79152899/article/details/139501513

均值滤波：

均值滤波的算法核心就是取一个滑动窗口，窗口是正方形的且保证边长为奇数（这样子就可以使得窗口中心必有一个元素

1	2	3
4	5	6
7	8	9

如图，我取了一个3X3的窗口，可知我现在要对中心元素5进行均值算法处理，顾名思义，我只需要将这九个元素全部相加再除以9，得到的值赋值给中心元素即可，新的值就是被均值处理过后的元素，以此类推。

这里有人可能会问如果是边界元素怎么办，这里我知道的到有两种方法，第一种是不管边界外的元素（本来也没有），只需正常将再窗口中的元素相加处理9就好，我使用的就是第一种。


	1	2
	4	5

例如我现在处理右上角的元素1，只需将1+2+4+5得到12，在除以9，得到值即为1的代替值

另外一种是对称延拓，在这里我简单说明一下就是将在窗口中的元素关于对称轴（边界）对称到边界外的元素点去，即在上图中，（我以数组下标为索引），（0，0）对称后为5，以此类推。

废话不多说上代码：

int kernelSize = 5;
    int kernelRadius = kernelSize / 2;
    U8 *filteredData = (U8 *)malloc(header.biSizeImage);
    memcpy(filteredData, data, header.biSizeImage); // 复制数据
    for (int y = 0; y < height; y++)
    {
        for (int x = 0; x < width; x++)
        {
            int sum = 0;
            int count = 0;
            for (int ky = -kernelRadius; ky <= kernelRadius; ky++)
            {
                for (int kx = -kernelRadius; kx <= kernelRadius; kx++)
                {
                    int nx = x + kx;
                    int ny = y + ky;
                    if (nx >= 0 && nx < width && ny >= 0 && ny < height)
                    {
                        sum += data[ny * width + nx];
                        count++;
                    }
                }
            }
            filteredData[y * width + x] = (U8)(sum / count);
        }
    }

高斯滤波：

高斯滤波的本质与均值滤波差不多，都是要求一个平均值来代替中心元素值，但是高斯滤波用到了高斯权矩阵与元素值相乘来使滤波平滑一点，高斯函数：，其中sigma可以决定滑动窗口的大小，x表示与窗口中心的距离。由高斯函数的出的矩阵还需进行归一化处理，也就是将已得出的矩阵中所有元素求和(设为sum)，再依次将每个元素除以和(sum)，得到的新矩阵即为高斯权矩阵。由于在进行高斯滤波时需要对X和Y两个方向进行滤波，如果在进行X方向的滤波同时进行Y方向的滤波将会非常耗时，故可以将滤波过程分解成两个一维方向的滤波，采取空间换取时间的方法提高运行效率。

// x方向一维运算
	for (int y = 0; y < height; y++)
	{
    	for (int x = 0; x < width; x++)
    	{
       		gaussSum = 0;
        	// 加权求和
        	for (int j = -kernelRadius; j <= kernelRadius; j++)
        	{
            	int idx = x + j;
            	if (idx >= 0 && idx < width)
                	gaussSum += gaussMatrix[j + kernelRadius] * data[y * width + idx];
        	}
        	// 记录X方向处理结果
        	filteredData[y * width + x] = (U8)gaussSum;
    	}
	}

	// y方向一维运算
	for (int x = 0; x < width; x++)
	{
    	for (int y = 0; y < height; y++)
    	{
        	gaussSum = 0;
        	// 加权求和
        	for (int j = -kernelRadius; j <= kernelRadius; j++)
        	{
            	int idy = y + j;
            	if (idy >= 0 && idy < height)
                	gaussSum += gaussMatrix[j + kernelRadius] * filteredData[idy * width + x];
        	}
        	// 更改原图像
        	data[y * width + x] = (U8)gaussSum;
    	}
	}
	//如果从第一个遍历到最后一个是height*width-1次，分成x与y两个方向则只需要width+height次

双边滤波：

双边滤波的算法需要用到两个矩阵，一个是高斯权矩阵，一个是像素差值矩阵，将两个矩阵对应元素相乘再进行归一化即可得到双边滤波之后的像素元素，公式如下：

sigmaS矩阵即为高斯权矩阵(又称空间域核)，sigmar即为像素差值矩阵(又称值域核)，像素差值矩阵就是当前遍历元素与窗口中心像素的差值所组成的。(注意：权矩阵并不是唯一的，因为高斯权矩阵可以更改sigma的值从而改变运算结果)。

空域权重 $\omega _{d}$ 和值域权重 $\omega _{r}$ 的意义

空域权重 $\omega _{d}$ 衡量的是两点之间的距离，距离越远权重越低；

值域权重 $\omega _{r}$ 衡量的是两点之间的像素值相似程度，越相似权重越大。

先计算出高斯权矩阵备用，待会儿计算时只需要查找即可

float gaussian(float x, float sigma)
{
    return exp(-(x * x) / (2 * sigma * sigma)) / (sqrt(2 * Pi) * sigma);
}

void createGaussianKernel(float *kernel, int windowSize, float sigma)
{
    int halfSize = windowSize / 2;
    float sum = 0.0;
    for (int i = -halfSize; i <= halfSize; i++)
    {
        for (int j = -halfSize; j <= halfSize; j++)
        {
            float value = gaussian(sqrt(i * i + j * j), sigma);
            kernel[(i + halfSize) * windowSize + (j + halfSize)] = value;
            sum += value;
        }
    }

    // 归一化
    for (int i = 0; i < windowSize * windowSize; i++)
    {
        kernel[i] /= sum;
    }
}

在编写双边滤波时我使用了上面提到过的对称延拓

// 对称延拓
	for (int i = 0; i < height + 2 * halfWindowSize; i++)
	{
    	for (int j = 0; j < width + 2 * halfWindowSize; j++)
    	{
        	int offsetX = i - halfWindowSize;
        	int offsetY = j - halfWindowSize;

        	// 计算对称位置的坐标
        	int symX = (offsetX < 0) ? -offsetX : ((offsetX >= height) ? (2 * height - offsetX - 1) : offsetX);
        	int symY = (offsetY < 0) ? -offsetY : ((offsetY >= width) ? (2 * width - offsetY - 1) : offsetY);

        	extendedData[i * (width + 2 * halfWindowSize) + j] = data[symX * width + symY];
    	}
	}

最后进行双边滤波

// 双边滤波参数
    int windowSize = 5;
    int halfWindowSize = windowSize/2; 
    float sigmaSpatial = 2.0; // 空域高斯标准差
    float sigmaRange = 30.0; // 值域高斯标准差
// 双边滤波处理
    for (int x = halfWindowSize; x < height + halfWindowSize; x++)
    {
        for (int y = halfWindowSize; y < width + halfWindowSize; y++)
        {
            float sumW = 0.0;
            float sumGrayW = 0.0;

            for (int k = -halfWindowSize; k <= halfWindowSize; k++)
            {
                for (int l = -halfWindowSize; l <= halfWindowSize; l++)
                {
                    int neighborX = x + k;
                    int neighborY = y + l;
                    float intensityDifference = extendedData[x * (width + 2 * halfWindowSize) + y] - extendedData[neighborX * (width + 2 * halfWindowSize) + neighborY];
                    float rangeWeight = exp(-intensityDifference * intensityDifference / (2 * sigmaRange * sigmaRange));
                    float spatialWeight = gaussianKernel[(k + halfWindowSize) * windowSize + (l + halfWindowSize)];
                    float weight = spatialWeight * rangeWeight;
                    sumGrayW += extendedData[neighborX * (width + 2 * halfWindowSize) + neighborY] * weight;
					sumW += weight;
                }
            }

            filteredData[(x - halfWindowSize) * width + (y - halfWindowSize)] = (U8)(sumGrayW / sumW);
        }
    }

代码中的sigmaSpatial可理解为空域权重，sigmaRange可理解为值域权重。