Truncated gaussian kernel implementation Matlab,right?

Question

I have the defination of Truncated gaussian kernel as:

So I confuse which is correct implementation of truncated gaussian kernel. Let see two case and let me know, thank you so much

Case 1:

G_truncated=fspecial('gaussian',round(2*sigma)*2 + 1,sigma); %  kernel

Case 2:

G=fspecial('gaussian',round(2*sigma)*2 + 1,sigma); %  normal distribution kernel
B = ones(round(2*sigma)*2 + 1,round(2*sigma)*2 + 1);
G_truncated=G.*B;
G_truncated = G_truncated/sum(G_truncated(:)); %normalized for sum=1

Answer 1

The answer of rayryeng is very useful for me. I only extend the gaussian kernel to ball kernel. The ball kernel is defined :

So based on answer of rayryeng. We can do it by

sigma=2;
rho=sigma;
tau=sigma;
%// Find grid of points
[X,Y] = meshgrid(-rho : rho, -rho : rho)

dists = (X.^2 + Y.^2); %// Find distances from the centre (Euclidean distance squared)
ballVal=dists;
ballVal(dists>sigma)=0;
ballVal(dists<=sigma)=1;

%// Now normalize
ballMask = ballVal / (sum(ballVal(:)));

Let me know, if it has any error or problem. Thank you

Answer 2

To add on to the previous post, there is a question of how to implement the kernel. You could use fspecial, truncate the kernel so that anything outside of the radius is zero, then renormalize it, but I'm assuming you'll want to do this from first principles.... so let's figure that out then. First, you need to generate a spatial map of distances from the centre of the mask. In conjunction, you use this to figure out what the Gaussian values (un-normalized) would be. You filter out those values in the un-normalized mask based on the spatial map of distances, then normalize that. As such, given your standard deviation tau, and your radius rho, you can do this:

%// Find grid of points
[X,Y] = meshgrid(-rho : rho, -rho : rho)

dists = (X.^2 + Y.^2); %// Find distances from the centre (Euclidean distance squared)
gaussVal = exp(-dists / (2*tau*tau)); %// Find unnormalized Gaussian values

%// Filter out those locations that are outside radius and set to 0
gaussVal(dists > rho^2) = 0;

%// Now normalize
gaussMask = gaussVal / (sum(gaussVal(:)));

---

Here is an example with using rho = 2 and tau = 2 with the outputs at each stage:

Stage #1 - Find grid co-ordinates

>> X

X =

-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

>> Y

Y = 

-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

Step #2 - Find distances from centre and unnormalized Gaussian values

>> dists

dists =

 8     5     4     5     8
 5     2     1     2     5
 4     1     0     1     4
 5     2     1     2     5
 8     5     4     5     8

>> gaussVal

gaussVal =

0.3679    0.5353    0.6065    0.5353    0.3679
0.5353    0.7788    0.8825    0.7788    0.5353
0.6065    0.8825    1.0000    0.8825    0.6065
0.5353    0.7788    0.8825    0.7788    0.5353
0.3679    0.5353    0.6065    0.5353    0.3679

Step #3 - Filter out locations that don't belong within the radius and set to 0

>> gaussVal =

     0         0    0.6065         0         0
     0    0.7788    0.8825    0.7788         0
0.6065    0.8825    1.0000    0.8825    0.6065
     0    0.7788    0.8825    0.7788         0
     0         0    0.6065         0         0

Step #4 - Normalize so sum is equal to 1

>> gaussMask =

     0         0    0.0602         0         0
     0    0.0773    0.0876    0.0773         0
0.0602    0.0876    0.0993    0.0876    0.0602
     0    0.0773    0.0876    0.0773         0
     0         0    0.0602         0         0

---

To verify that the mask sums to 1, just do sum(gaussMask(:)) and you'll see it's equal to 1... more or less :)

Answer 3

Your definition of truncated gaussian kernel is different than how MATLAB truncates filter kernels, though it generally won't matter in practice for sizable d.

fspecial already returns truncated AND normalized filter, so the second case is redundant, because it generates exactly the same result as case 1.

From MATLAB help:

H = fspecial('gaussian',HSIZE,SIGMA) returns a rotationally
symmetric Gaussian lowpass filter  of size HSIZE with standard
deviation SIGMA (positive). HSIZE can be a vector specifying the
number of rows and columns in H or a scalar, in which case H is a
square matrix.
The default HSIZE is [3 3], the default SIGMA is 0.5.

You can use fspecial('gaussian',1,sigma) to generate a 1x1 filter and see that it is indeed normalized.

To generate a filter kernel that fits your definition, you need to make B in your second case a matrix that has ones in a circular area. A less strict (but nonetheless redundant in practice) solution is to use fspecial('disk',size) to truncate your gaussian kernel. Don't forget to normalize it in either case.