Реализация размытия по Гауссу - как вычислить матрицу свертки (ядро)

Мой вопрос очень близок к этому вопросу: Как размыть изображение по Гауссу без использования встроенных функций Гаусса?

Ответ на этот вопрос очень хорош, но он не дает примера реального вычисления ядра гауссовского фильтра. Ответ дает произвольное ядро ​​и показывает, как применить фильтр, используя это ядро, но не как вычислить само реальное ядро. Я пытаюсь реализовать размытие по Гауссу в C ++ или Matlab с нуля, поэтому мне нужно знать, как рассчитать ядро ​​с нуля.

Я был бы признателен, если бы кто-нибудь мог вычислить реальное ядро ​​фильтра Гаусса, используя любую небольшую матрицу изображения примера.


person gsingh2011    schedule 20.11.2011    source источник
comment
Вы читали это: en.wikipedia.org/wiki/Gaussian_function?   -  person Oliver Charlesworth    schedule 21.11.2011
comment
Или даже это: en.wikipedia.org/wiki/Gaussian_blur   -  person Bart    schedule 21.11.2011
comment
Да, я потратил много времени, пытаясь понять это. Мне нужен пошаговый пример. После того, как я это пойму, я, вероятно, добавлю пример на страницу Gaussian Blur.   -  person gsingh2011    schedule 21.11.2011
comment
@ gsingh2011: Хорошая мысль, но, вероятно, бессмысленная. Ядра Гаусса обычно не создаются явно, потому что они отделимы, если корреляция равна 0.   -  person thiton    schedule 21.11.2011


Ответы (7)


Вы можете создать гауссовское ядро ​​с нуля, как указано в документации MATLAB fspecial . Прочтите формулу создания ядра Гаусса в части алгоритмов на этой странице и следуйте приведенному ниже коду. Код должен создать матрицу размером 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 для данного случая, поскольку вы, возможно, не знакомы с _ 7_, если вы пишете на C ++.

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
person petrichor    schedule 20.11.2011
comment
Я набрал [h1, h2] = meshgrid (- (m-1) / 2: (m-1) / 2, - (n-1) / 2: (n-1) / 2) и получил h1 в диапазоне от От -2 до 2, а не от -1,5 до 1,5. Та же проблема с h2. Но мой результат такой же. Кроме того, почему вы использовали значения сетки в качестве значений в формуле? Что это означает, если вы рассчитываете это для изображения? - person gsingh2011; 21.11.2011
comment
Ты прав! Я изменил m и n на 4, чтобы посмотреть, работает ли код, а затем скопировал значения для этого случая вместо того, чтобы давать им значение 5. Я исправил это, спасибо. - person petrichor; 21.11.2011
comment
Значения вычисляются в сетке, где точка в центре является началом координат, в нашем случае это h1 == 0 и h2 == 0. Все остальные пары представляют другие координаты, когда вы смотрите на значения h1, h2 поэлементно. Во время фильтрации вы можете подумать, что эта сетка будет помещена на пиксель изображения, где начало координат сетки точно соответствует пикселю. Вы можете прочитать ответ Гоза по ссылке, которую вы указали в своем вопросе, чтобы узнать подробности. - person petrichor; 21.11.2011
comment
Есть небольшая ошибка. Мы должны использовать [h1, h2] = meshgrid(-(n-1)/2:(n-1)/2, -(m-1)/2:(m-1)/2);, чтобы ответ был таким же с fspecial, когда m не равно n. - person Panfeng Li; 03.11.2016

Это так просто, как кажется:

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;
person thiton    schedule 20.11.2011
comment
Это недоработано: нужно также нормализовать ядро, иначе изображение станет темнее в зависимости от W и сигмы. Проще говоря: получить сумму значений ядра и разделить каждое значение ядра на эту сумму. - person Rookie; 24.05.2012
comment
@Rookie - Я решил изменить этот пост и добавить нормализацию. Это сделано для того, чтобы те, кто хочет, чтобы решение C / C ++ могли использовать это прямо. Хороший улов! - person rayryeng; 17.10.2014
comment
Это кажется неправильным, когда m, n четные, по сравнению с результатом fspecial. - person Panfeng Li; 03.11.2016

Чтобы реализовать размытие по Гауссу, вы просто берете gaussian function и вычислить одно значение для каждого элемента в вашем ядре.

Обычно вы хотите назначить максимальный вес центральному элементу в ядре и значения, близкие к нулю, для элементов на границах ядра. Это означает, что ядро ​​должно иметь нечетную высоту (соответственно ширину), чтобы гарантировать наличие центрального элемента.

Чтобы вычислить фактические элементы ядра, вы можете масштабировать гауссовский колокол до сетки ядра (выберите произвольный, например, sigma = 1, и произвольный диапазон, например, -2*sigma ... 2*sigma) и нормализовать его, s.t. сумма элементов равна одному. Для этого, если вы хотите поддерживать произвольные размеры ядра, вы можете адаптировать сигму к требуемому размеру ядра.

Вот пример 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-ядра. Это возможно благодаря ассоциативности этого типа линейной свертки (линейной отделимости). Вы также можете просмотреть этот раздел соответствующей статьи в Википедии.


Вот то же самое в 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']
person moooeeeep    schedule 20.11.2011

Размытие по Гауссу в Python с использованием библиотеки изображений PIL. Для получения дополнительной информации прочтите это: 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
person Vlad    schedule 13.12.2015

Хорошо, поздний ответ, но в случае ...

Используя ответ @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])

Чуть меньше строк.

person Metal3d    schedule 22.03.2019
comment
Я думаю, вы можете сэкономить еще больше строк, если используете scipy или np.convolve, как предлагается здесь: stackoverflow.com/q/29920114/1025391 < / а> - person moooeeeep; 27.04.2021

person    schedule
comment
Другим может быть полезно добавить некоторый контекст к вашему коду в вашем решении. - person Grant Miller; 02.06.2018

person    schedule
comment
Добавьте комментарии к своему коду, это будет полезно другим людям. - person HDJEMAI; 03.11.2016