我的问题非常接近这个问题:我如何在不使用任何内置高斯函数的情况下高斯模糊图像?
这个问题的答案非常好，但是它没有给出一个实际计算真实高斯滤波器核的例子。答案给出了一个任意的内核，并展示了如何使用该内核应用过滤器，但没有说明如何计算真正的内核本身。我正试图从头开始在c++或Matlab中实现高斯模糊，所以我需要知道如何从头开始计算内核。
如果有人能计算一个真正的高斯滤波器核使用任何小的例子图像矩阵。
您可以按照fspecial的MATLAB文档中所述从头创建高斯内核。请阅读该页算法部分的高斯核创建公式，并遵循下面的代码。代码是创建一个m × n矩阵，sigma = 1。
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
我认为在任何你喜欢的语言中实现它都是很简单的。
编辑:让我也为给定的情况添加h1和h2的值，因为如果您使用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
听起来很简单:
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;
要实现高斯模糊，只需取高斯函数并为内核中的每个元素计算一个值。
通常你想给内核的中心元素分配最大的权重，内核边界的元素的值接近于零。这意味着内核的高度应该是奇数。宽度)，以确保确实有一个中心元素。
为了计算实际的内核元素，你可以将高斯钟形缩放到内核网格(选择一个任意的例如sigma = 1和一个任意的范围，例如-2*sigma ... 2*sigma)并将其规范化，s.t.元素之和为1。为了实现这一点，如果您想支持任意内核大小，您可能需要调整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内核。更容易实现和更有效的计算是使用两个正交的一维核。这是可能的，由于这种类型的线性卷积的结合性(线性可分性)。您可能还想查看相应维基百科文章的这一部分。
在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']
好的，迟了一点答复，但是万一…
使用@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])
行数少一点。
使用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
// 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 matrixn");
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] = ' ';
if (str1[0] == 'q')
{
break;
}
int input1 = atoi(str1);
int window_size = input1 / 2;
printf("Input window_size was: %dn", 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] = ' ';
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: %fn", 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] = ' ';
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: %dn", cK);
printf("nwindow size=%d | Sigma = %f Coefficient K = %dnnn", 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 = %dnnn", 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 = %dnnn", 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 = %dnnn", 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);
}
}
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)
这是c#中的一个计算，它不从高斯(或另一个核)函数中获取单个样本，而是在一个小网格中计算大量样本，并将样本集成到所需的部分数量中。
计算是一维的，但可以很容易地扩展到二维。
这个计算使用了一些其他的函数，我没有在这里添加，但是我已经添加了函数签名，这样你就会知道它们是做什么的。
对于+/- 3的极限(积分的和areaSum为0.997300)，该计算产生以下离散值:
kernel size: normalized kernel values, rounded to 6 decimals
3: 0.157731, 0.684538, 0.157731
5: 0.034674, 0.238968, 0.452716, 0.238968, 0.034674
7: 0.014752, 0.083434, 0.235482, 0.332663, 0.235482, 0.083434, 0.014752
对于+/- 2的极限，该计算产生以下离散值(积分的和areaSum为0.954500):
kernel size: normalized kernel values, rounded to 6 decimals
3: 0.240694, 0.518612, 0.240694
5: 0.096720, 0.240449, 0.325661, 0.240449, 0.096720
7: 0.056379, 0.124798, 0.201012, 0.235624, 0.201012, 0.124798, 0.056379
using System.Linq;
private static void Main ()
{
int positionCount = 1024; // Number of samples in the range 0..1.
double positionStepSize = 1.0 / positionCount;
double limit = 3; // The calculation range of the kernel. +/- 3 is sufficient for gaussian.
int sectionCount = 3; // The number of elements in the kernel. May be 1, 3, 5, 7, ... (n*2+1)
// calculate the x positions for each kernel value to calculate.
double[] positions = CreateSeries (-limit, positionStepSize, (int)(limit * 2 * positionCount + 1));
// calculate the gaussian function value for each position
double[] values = positions.Select (pos => Gaussian (pos)).ToArray ();
// split the values into equal-sized sections and calculate the integral of each section.
double[] areas = IntegrateInSections (values, positionStepSize, sectionCount);
double areaSum = areas.Sum ();
// normalize to 1
double[] areas1 = areas.Select (a => a / areaSum).ToArray ();
double area1Sum = areas1.Sum (); // just to check it's 1 now
}
///-------------------------------------------------------------------
/// <summary>
/// Create a series of <paramref name="i_count"/> numbers, starting at <paramref name="i_start"/> and increasing by <paramref name="i_stepSize"/>.
/// </summary>
/// <param name="i_start">The start value of the series.</param>
/// <param name="i_stepSize">The step size between values in the series.</param>
/// <param name="i_count">The number of elements in the series.</param>
///-------------------------------------------------------------------
public static double[] CreateSeries (double i_start,
double i_stepSize,
int i_count)
{ ... }
private static readonly double s_gaussian_Divisor = Math.Sqrt (Math.PI * 2.0);
/// ------------------------------------------------------------------
/// <summary>
/// Calculate the value for the given position in a Gaussian kernel.
/// </summary>
/// <param name="i_position"> The position in the kernel for which the value will be calculated. </param>
/// <param name="i_bandwidth"> The width factor of the kernel. </param>
/// <returns> The value for the given position in a Gaussian kernel. </returns>
/// ------------------------------------------------------------------
public static double Gaussian (double i_position,
double i_bandwidth = 1)
{
double position = i_position / i_bandwidth;
return Math.Pow (Math.E, -0.5 * position * position) / s_gaussian_Divisor / i_bandwidth;
}
/// ------------------------------------------------------------------
/// <summary>
/// Calculate the integrals in the given number of sections of all given values with the given distance between the values.
/// </summary>
/// <param name="i_values"> The values for which the integral will be calculated. </param>
/// <param name="i_distance"> The distance between the values. </param>
/// <param name="i_sectionCount"> The number of sections in the integration. </param>
/// ------------------------------------------------------------------
public static double[] IntegrateInSections (IReadOnlyCollection<double> i_values,
double i_distance,
int i_sectionCount)
{ ... }