Implementazione della sfocatura gaussiana – Come calcolare la matrice di convoluzione (kernel)

La mia domanda è molto vicina a questa domanda: come posso sfocare un’immagine gaussiana senza utilizzare funzioni gaussiane integrate?

La risposta a questa domanda è molto buona, ma non fornisce un esempio di calcolo effettivo di un vero kernel di filtro gaussiano. La risposta fornisce un kernel arbitrario e mostra come applicare il filtro usando quel kernel ma non come calcolare un vero kernel. Sto provando a implementare una sfocatura gaussiana in C ++ o Matlab da zero, quindi ho bisogno di sapere come calcolare il kernel da zero.

Sarei grato se qualcuno potesse calcolare un vero kernel di filtro gaussiano usando una piccola matrice immagine di esempio.

È ansible creare un kernel gaussiano da zero come indicato nella documentazione MATLAB di fspecial . Si prega di leggere la formula di creazione del kernel gaussiana nella parte degli algoritmi in quella pagina e seguire il codice sottostante. Il codice è creare una matrice 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 

Osservare che questo può essere fatto dalla fspecial come segue:

 fspecial('gaussian', [mn], 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 

Penso che sia semplice implementarlo in qualsiasi lingua tu voglia.

EDIT: Consentitemi di aggiungere anche i valori di h1 e h2 per il caso dato, poiché potreste non avere familiarità con meshgrid se si meshgrid codice in 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 

È semplice come sembra:

 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; 

Per implementare la sfocatura gaussiana, prendi semplicemente la funzione gaussiana e calcola un valore per ciascuno degli elementi nel tuo kernel.

Normalmente si desidera assegnare il peso massimo all’elemento centrale nel kernel e valori prossimi allo zero per gli elementi ai bordi del kernel. Ciò implica che il kernel dovrebbe avere un’altezza dispari (resp. Larghezza) per garantire che ci sia effettivamente un elemento centrale.

Per calcolare gli effettivi elementi del kernel puoi ridimensionare la campana gaussiana alla griglia del kernel (scegli un arbitrario eg sigma = 1 e un intervallo arbitrario es -2*sigma ... 2*sigma ) e normalizzalo, st gli elementi sumno a uno . Per ottenere ciò, se si desidera supportare le dimensioni arbitrarie del kernel, si potrebbe voler adattare il sigma alla dimensione del kernel richiesta.

Ecco un esempio in C ++:

 #include  #include  #include  #include  double gaussian( double x, double mu, double sigma ) { const double a = ( x - mu ) / sigma; return std::exp( -0.5 * a * a ); } typedef std::vector kernel_row; typedef std::vector 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'; } } 

L'output è:

 $ 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 

Come semplificazione non è necessario utilizzare un kernel 2d. Più semplice da implementare e anche più efficiente da calcolare è utilizzare due kernel 1d ortogonali. Ciò è ansible a causa dell'associatività di questo tipo di convoluzione lineare (separabilità lineare). Potresti anche voler vedere questa sezione dell'articolo wikipedia corrispondente.


Ecco lo stesso in Python (nella speranza che qualcuno possa trovarlo utile):

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

Sfocatura gaussiana in python usando la libreria di immagini PIL. Per maggiori informazioni leggi questo: 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 
  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', [mn], sigma) 
 // my_test.cpp : Defines the entry point for the console application. // #include "stdafx.h" #include  #include  #include  #include  #include  //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 kernel_row; typedef std::vector 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); } }