【问题标题】:Implementing Gaussian Blur - How to calculate convolution matrix (kernel)实现高斯模糊 - 如何计算卷积矩阵(内核)
【发布时间】:2012-01-02 12:41:37
【问题描述】:

我的问题非常接近这个问题:How do I gaussian blur an image without using any in-built gaussian functions?

这个问题的答案很好,但是没有给出实际计算真正的高斯滤波器内核的例子。答案给出了一个任意内核,并展示了如何使用该内核应用过滤器,而不是如何计算真正的内核本身。我正在尝试从头开始在 C++ 或 Matlab 中实现高斯模糊,所以我需要知道如何从头开始计算内核。

如果有人可以使用任何小的示例图像矩阵计算出真正的高斯滤波器内核,我将不胜感激。

【问题讨论】:

  • 你读过这个吗:en.wikipedia.org/wiki/Gaussian_function
  • 是的,我花了很多时间试图理解这些。我需要的是一步一步的例子。理解之后,我大概会把例子加到 Gaussian Blur 页面上。
  • @gsingh2011:这是个好主意,但可能毫无意义。高斯核通常不会显式创建,因为如果相关性为 0,它们是可分离的。

标签: c++ matlab filter blur gaussian


【解决方案1】:

您可以按照fspecial 的 MATLAB 文档中的说明从头开始创建高斯内核。请阅读该页面算法部分中的高斯核创建公式,并按照以下代码进行操作。代码是创建一个 sigma = 1 的 m×n 矩阵。

m = 5; n = 5;
sigma = 1;
[h1, h2] = meshgrid(-(m-1)/2:(m-1)/2, -(n-1)/2:(n-1)/2);
hg = exp(- (h1.^2+h2.^2) / (2*sigma^2));
h = hg ./ sum(hg(:));

h =

    0.0030    0.0133    0.0219    0.0133    0.0030
    0.0133    0.0596    0.0983    0.0596    0.0133
    0.0219    0.0983    0.1621    0.0983    0.0219
    0.0133    0.0596    0.0983    0.0596    0.0133
    0.0030    0.0133    0.0219    0.0133    0.0030

请注意,这可以通过内置的fspecial 来完成,如下所示:

fspecial('gaussian', [m n], sigma)
ans =

    0.0030    0.0133    0.0219    0.0133    0.0030
    0.0133    0.0596    0.0983    0.0596    0.0133
    0.0219    0.0983    0.1621    0.0983    0.0219
    0.0133    0.0596    0.0983    0.0596    0.0133
    0.0030    0.0133    0.0219    0.0133    0.0030

我认为用任何你喜欢的语言来实现它都很简单。

编辑:让我也为给定的情况添加 h1h2 的值,因为如果您使用 C++ 编写代码,您可能不熟悉 meshgrid

h1 =

    -2    -1     0     1     2
    -2    -1     0     1     2
    -2    -1     0     1     2
    -2    -1     0     1     2
    -2    -1     0     1     2

h2 =

    -2    -2    -2    -2    -2
    -1    -1    -1    -1    -1
     0     0     0     0     0
     1     1     1     1     1
     2     2     2     2     2

【讨论】:

  • 我输入了 [h1,h2] = meshgrid(-(m-1)/2:(m-1)/2, -(n-1)/2:(n-1)/ 2) 得到的 h1 范围是 -2 到 2,而不是 -1.5 到 1.5。 h2也有同样的问题。但我的结果是一样的。另外,你为什么使用网格网格的值作为公式中的值?如果您为图像计算此值,这代表什么?
  • 你是对的!我将 mn 更改为 4,以查看代码是否有效,然后复制此案例的值,而不是为值 5 提供它们。我已修复它,谢谢。
  • 这些值是在一个网格上计算的,其中中心的一个是原点,在我们的例子中是 h1==0 和 h2==0。当您逐个元素查看 h1,h2 值时,所有其他对表示其他坐标。在过滤期间,您可以认为此网格将放置在图像的一个像素上,其中网格的原点恰好适合该像素。您可以在问题中提供的链接中阅读 Goz 的答案以了解详细信息。
  • 有个小错误。我们应该使用[h1, h2] = meshgrid(-(n-1)/2:(n-1)/2, -(m-1)/2:(m-1)/2); 这样当m 不等于n 时,答案与fspecial 相同。
【解决方案2】:

听起来很简单:

double sigma = 1;
int W = 5;
double kernel[W][W];
double mean = W/2;
double sum = 0.0; // For accumulating the kernel values
for (int x = 0; x < W; ++x) 
    for (int y = 0; y < W; ++y) {
        kernel[x][y] = exp( -0.5 * (pow((x-mean)/sigma, 2.0) + pow((y-mean)/sigma,2.0)) )
                         / (2 * M_PI * sigma * sigma);

        // Accumulate the kernel values
        sum += kernel[x][y];
    }

// Normalize the kernel
for (int x = 0; x < W; ++x) 
    for (int y = 0; y < W; ++y)
        kernel[x][y] /= sum;

【讨论】:

  • 这是有缺陷的:您也需要对内核进行归一化,否则图像会根据 W 和 sigma 变得更暗。简单地说:获取内核值的总和,然后将每个内核值除以该总和。
  • @Rookie - 我决定修改这篇文章并添加规范化。这是为了让那些想要 C/C++ 解决方案的人直接使用它。好收获!
  • fspecial的结果相比,当m,n为偶数时似乎不正确。
【解决方案3】:

要实现gaussian blur,您只需获取gaussian function 并为内核中的每个元素计算一个值。

通常您希望将最大权重分配给内核中的中心元素,并为内核边界处的元素分配接近零的值。 这意味着内核应该有一个奇数高度(分别是宽度)以确保实际上有一个中心元素。

要计算实际的内核元素,您可以将高斯钟形缩放到内核网格(选择任意范围,例如 sigma = 1 和任意范围,例如 -2*sigma ... 2*sigma)并将其标准化,s.t.元素之和为一。 为了实现这一点,如果您想支持任意内核大小,您可能需要调整 sigma 以适应所需的内核大小。

这是一个 C++ 示例:

#include <cmath>
#include <vector>
#include <iostream>
#include <iomanip>

double gaussian( double x, double mu, double sigma ) {
    const double a = ( x - mu ) / sigma;
    return std::exp( -0.5 * a * a );
}

typedef std::vector<double> kernel_row;
typedef std::vector<kernel_row> kernel_type;

kernel_type produce2dGaussianKernel (int kernelRadius) {
  double sigma = kernelRadius/2.;
  kernel_type kernel2d(2*kernelRadius+1, kernel_row(2*kernelRadius+1));
  double sum = 0;
  // compute values
  for (int row = 0; row < kernel2d.size(); row++)
    for (int col = 0; col < kernel2d[row].size(); col++) {
      double x = gaussian(row, kernelRadius, sigma)
               * gaussian(col, kernelRadius, sigma);
      kernel2d[row][col] = x;
      sum += x;
    }
  // normalize
  for (int row = 0; row < kernel2d.size(); row++)
    for (int col = 0; col < kernel2d[row].size(); col++)
      kernel2d[row][col] /= sum;
  return kernel2d;
}

int main() {
  kernel_type kernel2d = produce2dGaussianKernel(3);
  std::cout << std::setprecision(5) << std::fixed;
  for (int row = 0; row < kernel2d.size(); row++) {
    for (int col = 0; col < kernel2d[row].size(); col++)
      std::cout << kernel2d[row][col] << ' ';
    std::cout << '\n';
  }
}

输出是:

$ g++ test.cc && ./a.out
0.00134 0.00408 0.00794 0.00992 0.00794 0.00408 0.00134 
0.00408 0.01238 0.02412 0.03012 0.02412 0.01238 0.00408 
0.00794 0.02412 0.04698 0.05867 0.04698 0.02412 0.00794 
0.00992 0.03012 0.05867 0.07327 0.05867 0.03012 0.00992 
0.00794 0.02412 0.04698 0.05867 0.04698 0.02412 0.00794 
0.00408 0.01238 0.02412 0.03012 0.02412 0.01238 0.00408 
0.00134 0.00408 0.00794 0.00992 0.00794 0.00408 0.00134 

为了简化,您不需要使用 2d 内核。更容易实现且计算效率更高的是使用两个正交的一维内核。由于这种类型的线性卷积(线性可分离性)的关联性,这是可能的。 您可能还想查看相应维基百科文章的this section


在 Python 中也是如此(希望有人会觉得它有用):

from math import exp

def gaussian(x, mu, sigma):
  return exp( -(((x-mu)/(sigma))**2)/2.0 )

#kernel_height, kernel_width = 7, 7
kernel_radius = 3 # for an 7x7 filter
sigma = kernel_radius/2. # for [-2*sigma, 2*sigma]

# compute the actual kernel elements
hkernel = [gaussian(x, kernel_radius, sigma) for x in range(2*kernel_radius+1)]
vkernel = [x for x in hkernel]
kernel2d = [[xh*xv for xh in hkernel] for xv in vkernel]

# normalize the kernel elements
kernelsum = sum([sum(row) for row in kernel2d])
kernel2d = [[x/kernelsum for x in row] for row in kernel2d]

for line in kernel2d:
  print ["%.3f" % x for x in line]

产生内核:

['0.001', '0.004', '0.008', '0.010', '0.008', '0.004', '0.001']
['0.004', '0.012', '0.024', '0.030', '0.024', '0.012', '0.004']
['0.008', '0.024', '0.047', '0.059', '0.047', '0.024', '0.008']
['0.010', '0.030', '0.059', '0.073', '0.059', '0.030', '0.010']
['0.008', '0.024', '0.047', '0.059', '0.047', '0.024', '0.008']
['0.004', '0.012', '0.024', '0.030', '0.024', '0.012', '0.004']
['0.001', '0.004', '0.008', '0.010', '0.008', '0.004', '0.001']

【讨论】:

    【解决方案4】:

    使用 PIL 图像库在 python 中进行高斯模糊。欲了解更多信息,请阅读:http://blog.ivank.net/fastest-gaussian-blur.html

    from PIL import Image
    import math
    
    # img = Image.open('input.jpg').convert('L')
    # r = radiuss
    def gauss_blur(img, r):
        imgData = list(img.getdata())
    
        bluredImg = Image.new(img.mode, img.size)
        bluredImgData = list(bluredImg.getdata())
    
        rs = int(math.ceil(r * 2.57))
    
        for i in range(0, img.height):
            for j in range(0, img.width):
                val = 0
                wsum = 0
                for iy in range(i - rs, i + rs + 1):
                    for ix in range(j - rs, j + rs + 1):
                        x = min(img.width - 1, max(0, ix))
                        y = min(img.height - 1, max(0, iy))
                        dsq = (ix - j) * (ix - j) + (iy - i) * (iy - i)
                        weight = math.exp(-dsq / (2 * r * r)) / (math.pi * 2 * r * r)
                        val += imgData[y * img.width + x] * weight
                        wsum += weight 
                bluredImgData[i * img.width + j] = round(val / wsum)
    
        bluredImg.putdata(bluredImgData)
        return bluredImg
    

    【讨论】:

      【解决方案5】:
      // my_test.cpp : Defines the entry point for the console application.
      //
      
      #include "stdafx.h"
      
      #include <cmath>
      #include <vector>
      #include <iostream>
      #include <iomanip>
      #include <string>
      
      //https://stackoverflow.com/questions/8204645/implementing-gaussian-blur-how-to-calculate-convolution-matrix-kernel
      //https://docs.opencv.org/2.4/modules/imgproc/doc/filtering.html#getgaussiankernel
      //http://dev.theomader.com/gaussian-kernel-calculator/
      
      double gaussian(double x, double mu, double sigma) {
          const double a = (x - mu) / sigma;
          return std::exp(-0.5 * a * a);
      }
      
      typedef std::vector<double> kernel_row;
      typedef std::vector<kernel_row> kernel_type;
      
      kernel_type produce2dGaussianKernel(int kernelRadius, double sigma) {
          kernel_type kernel2d(2 * kernelRadius + 1, kernel_row(2 * kernelRadius + 1));
          double sum = 0;
          // compute values
          for (int row = 0; row < kernel2d.size(); row++)
              for (int col = 0; col < kernel2d[row].size(); col++) {
                  double x = gaussian(row, kernelRadius, sigma)
                      * gaussian(col, kernelRadius, sigma);
                  kernel2d[row][col] = x;
                  sum += x;
              }
          // normalize
          for (int row = 0; row < kernel2d.size(); row++)
              for (int col = 0; col < kernel2d[row].size(); col++)
                  kernel2d[row][col] /= sum;
          return kernel2d;
      }
      
      char* gMatChar[10] = {
          "          ",
          "         ",
          "        ",
          "       ",
          "      ",
          "     ",
          "    ",
          "   ",
          "  ",
          " "
      };
      
      static int countSpace(float aValue)
      {
          int count = 0;
          int value = (int)aValue;
          while (value > 9)
          {
              count++;
              value /= 10;
          }
          return count;
      }
      
      int main() {
          while (1)
          {
              char str1[80]; // window size
              char str2[80]; // sigma
              char str3[80]; // coefficient
              int space;
      
              int i, ch;
              printf("\n-----------------------------------------------------------------------------\n");
              printf("Start generate Gaussian matrix\n");
              printf("-----------------------------------------------------------------------------\n");
              // input window size
              printf("\nPlease enter window size (from 3 to 10) It should be odd (ksize/mod 2 = 1 ) and positive: Exit enter q \n");
              for (i = 0; (i < 80) && ((ch = getchar()) != EOF)
                  && (ch != '\n'); i++)
              {
                  str1[i] = (char)ch;
              }
      
              // Terminate string with a null character
              str1[i] = '\0';
              if (str1[0] == 'q')
              {
                  break;
              }
              int input1 = atoi(str1);
              int window_size = input1 / 2;
              printf("Input window_size was: %d\n", input1);
      
              // input sigma
              printf("Please enter sigma. Use default press Enter . Exit enter q \n");
              str2[0] = '0';
              for (i = 0; (i < 80) && ((ch = getchar()) != EOF)
                  && (ch != '\n'); i++)
              {
                  str2[i] = (char)ch;
              }
      
              // Terminate string with a null character
              str2[i] = '\0';
              if (str2[0] == 'q')
              {
                  break;
              }
              float input2 = atof(str2);
              float sigma;
              if (input2 == 0)
              {
                  // Open-CV sigma � Gaussian standard deviation. If it is non-positive, it is computed from ksize as sigma = 0.3*((ksize-1)*0.5 - 1) + 0.8 .
                  sigma = 0.3*((input1 - 1)*0.5 - 1) + 0.8;
              }
              else
              {
                  sigma = input2;
              }
              printf("Input sigma was: %f\n", sigma);
      
              // input Coefficient K
              printf("Please enter Coefficient K. Use default press Enter . Exit enter q \n");
              str3[0] = '0';
              for (i = 0; (i < 80) && ((ch = getchar()) != EOF)
                  && (ch != '\n'); i++)
              {
                  str3[i] = (char)ch;
              }
      
              // Terminate string with a null character
              str3[i] = '\0';
              if (str3[0] == 'q')
              {
                  break;
              }
              int input3 = atoi(str3);
              int cK;
              if (input3 == 0)
              {
                  cK = 1;
              }
              else
              {
                  cK = input3;
              }
              float sum_f = 0;
              float temp_f;
              int sum = 0;
              int temp;
              printf("Input Coefficient K was: %d\n", cK);
      
              printf("\nwindow size=%d | Sigma = %f Coefficient K = %d\n\n\n", input1, sigma, cK);
      
              kernel_type kernel2d = produce2dGaussianKernel(window_size, sigma);
              std::cout << std::setprecision(input1) << std::fixed;
              for (int row = 0; row < kernel2d.size(); row++) {
                  for (int col = 0; col < kernel2d[row].size(); col++)
                  {
                      temp_f = cK* kernel2d[row][col];
                      sum_f += temp_f;
                      space = countSpace(temp_f);
                      std::cout << gMatChar[space] << temp_f << ' ';
                  }
                  std::cout << '\n';
              }
              printf("\n Sum array = %f | delta = %f", sum_f, sum_f - cK);
      
              // rounding
              printf("\nRecommend use round(): window size=%d | Sigma = %f Coefficient K = %d\n\n\n", input1, sigma, cK);
              sum = 0;
              std::cout << std::setprecision(0) << std::fixed;
              for (int row = 0; row < kernel2d.size(); row++) {
                  for (int col = 0; col < kernel2d[row].size(); col++)
                  {
                      temp = round(cK* kernel2d[row][col]);
                      sum += temp;
                      space = countSpace((float)temp);
                      std::cout << gMatChar[space] << temp << ' ';
                  }
                  std::cout << '\n';
              }
              printf("\n Sum array = %d | delta = %d", sum, sum - cK);
      
              // recommented
              sum_f = 0;
              int cK_d = 1 / kernel2d[0][0];
              cK_d = cK_d / 2 * 2;
              printf("\nRecommend: window size=%d | Sigma = %f Coefficient K = %d\n\n\n", input1, sigma, cK_d);
              std::cout << std::setprecision(input1) << std::fixed;
              for (int row = 0; row < kernel2d.size(); row++) {
                  for (int col = 0; col < kernel2d[row].size(); col++)
                  {
                      temp_f = cK_d* kernel2d[row][col];
                      sum_f += temp_f;
                      space = countSpace(temp_f);
                      std::cout << gMatChar[space] << temp_f << ' ';
                  }
                  std::cout << '\n';
              }
              printf("\n Sum array = %f | delta = %f", sum_f, sum_f - cK_d);
      
              // rounding
              printf("\nRecommend use round(): window size=%d | Sigma = %f Coefficient K = %d\n\n\n", input1, sigma, cK_d);
              sum = 0;
              std::cout << std::setprecision(0) << std::fixed;
              for (int row = 0; row < kernel2d.size(); row++) {
                  for (int col = 0; col < kernel2d[row].size(); col++)
                  {
                      temp = round(cK_d* kernel2d[row][col]);
                      sum += temp;
                      space = countSpace((float)temp);
                      std::cout << gMatChar[space] << temp << ' ';
                  }
                  std::cout << '\n';
              }
              printf("\n Sum array = %d | delta = %d", sum, sum - cK_d);
      
          }
      }
      

      【讨论】:

      • 在您的解决方案中为您的代码添加一些上下文可能对其他人有所帮助。
      【解决方案6】:

      好的,迟到的答案,但万一……

      使用@moooeeeeep 答案,但使用 numpy;

      import numpy as np
      radius = 3
      sigma = radius/2.
      
      k = np.arange(2*radius +1)
      row = np.exp( -(((k - radius)/(sigma))**2)/2.)
      col = row.transpose()
      out = np.outer(row, col)
      out = out/np.sum(out)
      for line in out:
          print(["%.3f" % x for x in line])
      

      行数少一点。

      【讨论】:

      【解决方案7】:
       function kernel = gauss_kernel(m, n, sigma)
       % Generating Gauss Kernel
      
       x = -(m-1)/2 : (m-1)/2;
       y = -(n-1)/2 : (n-1)/2;
      
       for i = 1:m
           for j = 1:n
               xx(i,j) = x(i);
               yy(i,j) = y(j);
           end
       end
      
       kernel = exp(-(xx.*xx + yy.*yy)/(2*sigma*sigma));
      
       % Normalize the kernel
       kernel  = kernel/sum(kernel(:));
      
       % Corresponding function in MATLAB
       % fspecial('gaussian', [m n], sigma)
      

      【讨论】:

      • 在您的代码中添加一些 cmets,对其他人会有所帮助。
      猜你喜欢
      • 2011-03-22
      • 1970-01-01
      • 1970-01-01
      • 2021-12-15
      • 2013-10-25
      • 2013-09-15
      • 1970-01-01
      • 2014-05-25
      • 1970-01-01
      相关资源
      最近更新 更多