实现高斯模糊 - 如何计算卷积矩阵(内核)

问题描述 投票:32回答:7

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

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

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

c++ matlab filter blur gaussian
7个回答
35
投票

您可以从头开始创建高斯内核,如fspecial的MATLAB文档中所述。请阅读该页面算法部分中的高斯内核创建公式,并遵循以下代码。代码是用sigma = 1创建m-by-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

28
投票

这听起来很简单:

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;

17
投票

要实现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内核。更容易实现并且计算效率更高是使用两个正交的1d内核。由于这种类型的线性卷积(线性可分离性)的相关性,这是可能的。您可能还想查看相应维基百科文章的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']

0
投票

使用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

0
投票
 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)

0
投票
// 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);

    }
}

0
投票

好的,迟到的答案,但在...的情况下

使用@moooeeeep答案,但有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])

只是少了几行。

© www.soinside.com 2019 - 2024. All rights reserved.