c++ - fspecial - Implementación de desenfoque gaussiano: cómo calcular la matriz de convolución(kernel)

gaussian filter matlab

Question

Mi pregunta es muy similar a esta pregunta: ¿cómo difumino Gauss una imagen sin utilizar ninguna función gaussiana incorporada?

La respuesta a esta pregunta es muy buena, pero no da un ejemplo de calcular realmente un kernel de filtro gaussiano real. La respuesta proporciona un kernel arbitrario y muestra cómo aplicar el filtro utilizando ese núcleo, pero no cómo calcular un kernel real. Estoy tratando de implementar un desenfoque gaussiano en C ++ o Matlab desde cero, así que necesito saber cómo calcular el kernel desde cero.

Agradecería que alguien pudiera calcular un núcleo de filtro Gaussiano real usando cualquier matriz de imagen de ejemplo pequeña.

Answer 1

Desenfoque gaussiano en python utilizando la biblioteca de imágenes PIL. Para obtener más información, lea esto: 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

Answer 2

Es tan simple como suena:

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;

Answer 3

Para implementar el en.wikipedia.org/wiki/Gaussian_blur , simplemente toma la en.wikipedia.org/wiki/Gaussian_function y calcula un valor para cada uno de los elementos en tu kernel.

Por lo general, desea asignar el peso máximo al elemento central en su kernel y valores cercanos a cero para los elementos en los bordes del kernel. Esto implica que el núcleo debe tener una altura impar (o ancho) para garantizar que realmente exista un elemento central.

Para calcular los elementos reales del núcleo, puede escalar la campana gaussiana a la cuadrícula del kernel (elija un arbitrario, por ejemplo, sigma = 1 y un rango arbitrario, por ejemplo, -2*sigma ... 2*sigma ) y normalícelo, st los elementos suman uno . Para lograr esto, si desea admitir tamaños de kernel arbitrarios, es posible que desee adaptar sigma al tamaño de kernel requerido.

Aquí hay un ejemplo de 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'';
  }
}

El resultado es:

$ 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 

Como simplificación, no necesita usar un kernel 2d. Más fácil de implementar y también más eficiente de calcular es usar dos núcleos 1d ortogonales. Esto es posible debido a la asociatividad de este tipo de convolución lineal (separabilidad lineal). También es posible que desee ver esta sección del artículo de wikipedia correspondiente.

Esto es lo mismo en Python (con la esperanza de que alguien lo encuentre útil):

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]

produce el kernel:

[''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'']

Answer 4

Puede crear un núcleo gaussiano desde cero, como se indica en la documentación de MATLAB de fspecial . Lea la fórmula de creación de kernel de Gauss en la parte de algoritmos de esa página y siga el código que se muestra a continuación. El código es para crear una matriz m-by-n con 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

Observe que esto se puede hacer con el fspecial siguiente manera:

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

Creo que es sencillo implementar esto en el idioma que desee.

EDITAR: Permítame también agregar los valores de h1 y h2 para el caso dado, ya que puede que no esté familiarizado con meshgrid si codifica en 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

Answer 5

 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)

Answer 6

// 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://.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);

    }
}