41

My question is very close to this question: How do I gaussian blur an image without using any in-built gaussian functions?

The answer to this question is very good, but it doesn't give an example of actually calculating a real Gaussian filter kernel. The answer gives an arbitrary kernel and shows how to apply the filter using that kernel but not how to calculate a real kernel itself. I am trying to implement a Gaussian blur in C++ or Matlab from scratch, so I need to know how to calculate the kernel from scratch.

I'd appreciate it if someone could calculate a real Gaussian filter kernel using any small example image matrix.

Community
  • 1
  • 1
gsgx
  • 12,020
  • 25
  • 98
  • 149
  • Have you read this: http://en.wikipedia.org/wiki/Gaussian_function? – Oliver Charlesworth Nov 20 '11 at 21:01
  • Or even this: http://en.wikipedia.org/wiki/Gaussian_blur – Bart Nov 20 '11 at 21:10
  • 1
    Yes, I've spent a lot of time trying to understand those. What I need is a step by step example. After I understand it, I'll probably add the example to the Gaussian Blur page. – gsgx Nov 20 '11 at 21:24
  • @gsingh2011: This is a nice thought, but probably pointless. Gaussian kernels are not created explicitly normally because they are separable if the correlation is 0. – thiton Nov 20 '11 at 21:33

8 Answers8

37

You can create a Gaussian kernel from scratch as noted in MATLAB documentation of fspecial. Please read the Gaussian kernel creation formula in the algorithms part in that page and follow the code below. The code is to create an m-by-n matrix with 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 that this can be done by the built-in fspecial as follows:

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

I think it is straightforward to implement this in any language you like.

EDIT: Let me also add the values of h1 and h2 for the given case, since you may be unfamiliar with meshgrid if you code 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
petrichor
  • 6,459
  • 4
  • 36
  • 48
  • 1
    I typed [h1,h2] = meshgrid(-(m-1)/2:(m-1)/2, -(n-1)/2:(n-1)/2) and got a h1 ranging from -2 to 2, not from -1.5 to 1.5. Same problem with h2. But my result is the same. Also, why did you use the values of mesh grid as the values in the formula? What does this represent if you were calculating this for an image? – gsgx Nov 20 '11 at 22:06
  • You're right! I changed `m` and `n` to 4 in order to see whether the code works and then copied the values for this case instead of giving them for value 5. I've fixed it, thanks. – petrichor Nov 20 '11 at 22:15
  • The values are computed on a grid where the one in the center is the origin, which is h1==0 and h2==0 in our case. All the other pairs represent the other coordinates when you look at the h1,h2 values element by element. During filtering, you can think that this grid will be put on a pixel of the image where the origin of the grid is fit exactly on the pixel. You can read Goz's answer in the link you gave in your question for the details. – petrichor Nov 20 '11 at 22:20
  • There is a small mistake. We should use `[h1, h2] = meshgrid(-(n-1)/2:(n-1)/2, -(m-1)/2:(m-1)/2);` So that the answer will the same with `fspecial` when m not equal n. – Panfeng Li Nov 02 '16 at 23:00
33

It's as simple as it sounds:

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;
rayryeng
  • 102,964
  • 22
  • 184
  • 193
thiton
  • 35,651
  • 4
  • 70
  • 100
  • 11
    This is flawed: you need to normalize the kernel too, or the image becomes darker depending on W and sigma. simply put: get the sum of kernel values and divide each kernel value by that sum. – Rookie May 24 '12 at 14:38
  • 3
    @Rookie - I've decided to modify this post and add the normalization. This is to allow for those who want a C / C++ solution to use this straight up. Good catch! – rayryeng Oct 17 '14 at 05:12
  • 1
    It seems not correct when m, n are even, comparing to the result of `fspecial`. – Panfeng Li Nov 02 '16 at 23:01
  • This solution is definitely wrong, because it just calculates the kernel for a larger range, but it does not 'spread' the kernel, i.e. it does not adjust the bandwidth (sigma). – Tobias Knauss Apr 22 '22 at 16:31
24

To implement the gaussian blur you simply take the gaussian function and compute one value for each of the elements in your kernel.

Usually you want to assign the maximum weight to the central element in your kernel and values close to zero for the elements at the kernel borders. This implies that the kernel should have an odd height (resp. width) to ensure that there actually is a central element.

To compute the actual kernel elements you may scale the gaussian bell to the kernel grid (choose an arbitrary e.g. sigma = 1 and an arbitrary range e.g. -2*sigma ... 2*sigma) and normalize it, s.t. the elements sum to one. To achieve this, if you want to support arbitrary kernel sizes, you might want to adapt the sigma to the required kernel size.

Here's a C++ example:

#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';
  }
}

The output is:

$ 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 

As a simplification you don't need to use a 2d-kernel. Easier to implement and also more efficient to compute is to use two orthogonal 1d-kernels. This is possible due to the associativity of this type of a linear convolution (linear separability). You may also want to see this section of the corresponding wikipedia article.


Here's the same in Python (in the hope someone might find it useful):

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]

produces the 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']
Soonts
  • 20,079
  • 9
  • 57
  • 130
moooeeeep
  • 31,622
  • 22
  • 98
  • 187
  • I think it's not sufficient to use just one sample value of the kernel function (let it be gaussian or another kernel). The values in the center would be too small, and the values on the borders would be too large. You should sample in a much smaller grid (e.g. 0.001) in the range -3...+3, split the range into the desired number of kernel elements, and integrate all values in each range. This is also suggested by https://homepages.inf.ed.ac.uk/rbf/HIPR2/gsmooth.htm. I have added a C# example below: https://stackoverflow.com/a/71968356/2505186 – Tobias Knauss Apr 22 '22 at 12:09
  • @TobiasKnauss Your suggested process results in the convolution of a Gaussian with a small square window, and therefore is not a Gaussian any more. Point sampling is the correct procedure. – Cris Luengo Apr 22 '22 at 13:46
2

OK, a late answer but in case of...

Using the @moooeeeep answer, but with 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])

Just a bit less of lines.

Metal3d
  • 2,905
  • 1
  • 23
  • 29
  • I think you can save even more lines if you use scipy or np.convolve as suggested here: https://stackoverflow.com/q/29920114/1025391 – moooeeeep Apr 27 '21 at 06:57
1

Gaussian blur in python using PIL image library. For more info read this: 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
Vlad
  • 11
  • 2
1
// 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
 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)
Panfeng Li
  • 3,321
  • 3
  • 26
  • 34
0

Here's a calculation in C#, which does not take single samples from the gaussian (or another kernel) function, but it calculates a large number of samples in a small grid and integrates the samples in the desired number of sections.
The calculation is for 1D, but it may easily be extended to 2D.

This calculation uses some other functions, which I did not add here, but I have added the function signatures so that you will know what they do.

This calculation produces the following discrete values for the limits +/- 3 (sum areaSum of integral is 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

This calculation produces the following discrete values for the limits +/- 2 (sum areaSum of integral is 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

Code:

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)
{ ... }
Tobias Knauss
  • 3,361
  • 1
  • 21
  • 45
  • Repeating my earlier comment here: This process results in the convolution of a Gaussian with a small square window, and therefore is not a Gaussian any more. Point sampling is the correct procedure. – Cris Luengo Apr 22 '22 at 13:47
  • @CrisLuengo: https://homepages.inf.ed.ac.uk/rbf/HIPR2/gsmooth.htm suggests an integration, which makes more sense to me than taking samples. – Tobias Knauss Apr 22 '22 at 14:33
  • Wouldn't be the first misconception on that site. All of the signal processing theory (Fourier analysis etc) that supports convolution is based on point sampling of signals. Integrating a signal over the area of a pixel is equivalent to point-sampling the signal after convolving it with a pixel-sized uniform kernel. It is clear to me that this will make the result less correct, not more correct. – Cris Luengo Apr 22 '22 at 14:41