matlab 均值滤波_数字图像处理基础 — 高斯滤波

高斯滤波，本文主要讲其如何通过 C语言实现。不太擅长写理论性质的文章，这里仅仅阐述自己怎么实现以及简单的优化过程。

通常我们对获取的图像进行进一步处理时，往往需要先进行一次降噪，而通常我们选择的是高斯滤波也称高斯模糊，其因为较好的降噪效果而被广泛使用，具体的原因大家可以做实验以及看相关的论文去理解。

如果你经常使用MATLAB或者Python，可能只需要调用一下库函数即可实现。以matlab为例，其高斯滤波的常见使用，如下：

% Generate Kernel[5, 5], sigma = 1.0.
gaussKernel = fspecial('gaussian', [5, 5], 1.0);
% Compute gaussian filter result.
gaussResult = imfilter(Image, gaussKernel, 'replicate');

真的是方便至极，但自己不实现一下，好像没什么意思。并且嵌入式设备没有相应的库函数让你使用，所有还是有必要对其进行理解和学习。以下为正文分割线。

理论公式与模板生成

高斯滤波作为一种线性滤波器，其本质是带权值的加权均值滤波，相对于均值滤波，高斯滤波除了中心位置为1，其他位置根据距离的远近，越远滤波系数权重越低。

对高斯滤波的公式进行理解学习，参考如下公式：

其中，

为模板中心坐标，即坐标系原点。有了该公式，我们就可以实现高斯模板的生成了，参考如下

C代码：

void kernelGenerate(float *gaussKernel, int k_size, float k_sigma) {
 
 float flag = 1 / (2 * PI * k_sigma * k_sigma);
 float coef;
 float sum  = 0;
 int center = (k_size - 1) / 2;

 for (int i = 0; i < k_size; i++) {
  float x2 = (float)pow((i - center), 2);
  for (int j = 0; j < k_size; j++) {
   float y2 = (float)pow((j - center), 2);
   coef = (-x2 - y2) * flag;
   gaussKernel[i * k_size + j] = flag * exp(coef);
   sum += gaussKernel[i * k_size + j];
  }
 }

 /* Normalize */
 if (sum) {
  for (int i = 0; i < k_size * k_size; i++) {
   gaussKernel[i] /= sum;
  }
 }
}

这个模板是生成的小数模板，要是生成整数模板呢？如下，将左上角置1。然后再进行取整操作。

void kernelGenerate(float *gaussKernel, int k_size, float k_sigma) {
 
 float flag = 1 / (2 * PI * k_sigma * k_sigma);
 float coef;
 int center = (k_size - 1) / 2;

 for (int i = 0; i < k_size; i++) {
  float x2 = (float)pow((i - center), 2);
  for (int j = 0; j < k_size; j++) {
   float y2 = (float)pow((j - center), 2);
   coef = (-x2 - y2) * flag;
   gaussKernel[i * k_size + j] = flag * exp(coef);
  }
 }

 /* Intergel Normalize */
 if (gaussKernel[0]) {
                float normal = 1.0 / gaussKernel[0];
  for (int i = 0; i < k_size * k_size; i++) {
   gaussKernel[i] *= normal;
  }
 }
}

有了模板之后，就需要进行滤波了。通过公式，可以发现滤波最关键的就是滤波的核的标准差，因此需要根据实际的需要对其进行相应的设置。

滤波处理

如果直接使用公式，进行滤波则可以写出如下代码：

void gaussian_filter(unsigned char *src, unsigned char *dst, int src_w, int src_h, int k_size, float k_sigma, int padd_type, int out_type) {

 int img_w, img_h;

 if (k_size == 1)
  dst = src;

 if (k_size % 2 == 0) {
  printf("The kernel size not odd and ksize > 1! Convert it to k_size + 1");
  k_size += 1;
 }
 int padd_size = (k_size - 1) / 2;

 img_w = src_w + padd_size * 2;
 img_h = src_h + padd_size * 2;

 unsigned char *pBlur = (unsigned char *)malloc(img_w * img_h);
 unsigned char *pDst  = dst;
 float *gaussKernel   = (float *)malloc(sizeof(float) * k_size * k_size);

 for (int i = 0; i < src_h; i++) {
  unsigned char *pTmp = pBlur + (padd_size + i)* img_w + padd_size;
  memcpy(pTmp, src + i * src_w, src_w);

  /* col padding */
  for (int j = 0; j < padd_size; j++) {
   if (padd_type == 0) { // zero padding.
    pTmp[j - padd_size] = 0;
    pTmp[src_w + j] = 0;
   } else {  // reflect padding.
    pTmp[j - padd_size] = pTmp[padd_size - j];
    pTmp[src_w + padd_size - j] = pTmp[src_w - padd_size + j];
   }
  }
 }

 /* rows padding. */
 if (padd_type == 0) {
  memset(pBlur, 0, padd_size * img_w);
  memset(pBlur + img_w * (img_h - 1 - padd_size), 0, padd_size * img_w);
 } else {
  for (int i = 0; i < padd_size; i++) {
   memcpy(pBlur + i * img_w, pBlur + (k_size - 1 - i) * img_w, img_w);
   memcpy(pBlur + (img_h - 1 - i) * img_w, pBlur + (img_h - k_size + i) * img_w, img_w);
  }
 }

 /* Generate filter kernel. */
 /* Eq: 1 / (2 * pi * Sigma^2) * exp( (x - m)^2 + (y - m)^2) */
 kernelGenerate(gaussKernel, k_size, k_sigma);

 /* Gaussian filter. */
 for (int row = padd_size; row < img_h - padd_size; row++) {
  for (int col = padd_size; col < img_w - padd_size; col++) {
   float sum = 0;
   for (int i = -padd_size; i < padd_size; i++) {
    for (int j = -padd_size; j < padd_size; j++) {
     sum += pBlur[(row + i) * img_w + col + j] * gaussKernel[(i + padd_size) * k_size + j + padd_size];
    }
   }
   pDst[(row - padd_size) * src_w + col - padd_size] = sum;
  }
 }

 free(gaussKernel);
 free(pBlur);
}