Scribbling ….

Generic multivariate Gaussian kernel in any derivative order

Posted in Matlab by avan on May 27, 2010

Matlab’s image processing toolbox has fspecial function to create several 2D kernels, e.g., gaussian, laplacian, sobel, prewitt, etc. that can be used to filter an image, but I want more than that. I need a Gaussian kernel in any dimension (multivariate) and also in any derivative order. For instances, to create a triangular Hessian matrix, i.e.,

where g is an N-d Gaussian kernel, then I need up to 2nd order partial derivatives of the Gaussian kernel. So here’s my solution to Matlab.

Let’s define the Gaussian function as follows:

There are two important properties of a Gaussian function that helps me to create this program:

  1. A product of two Gaussian functions is another Gaussian function. So, I can create 2D Gaussian function from two 1D Gaussian functions as follows:
    Hence also for ND Gaussian function:
  2. A partial derivative of a Gaussian function can be generated by using Hermite polynomial:

    [see my previous post].
    Thus the hermite function that I’ve created before is a pre-requisite for this Gaussian kernel.

Now, all ingredients to create a multivariate Gaussian function are ready. Let’s assume we have created a grid. For instances, for 3D:

[xi,yi,zi] = meshgrid(-5:5);
X = [xi(:) yi(:) zi(:)];

then the algorithm to create multivariate Gaussian kernel in any derivative order can be written as follows:

% normalized 1D gaussian function
gfun = @(x,s) exp(-(x.*x)/(2*s.*s)) ./ (sqrt(2*pi)*s);

% create kernel
for i=1:dim
    c = sigma(i) * sqrt(2);
    gx = gfun(X(:,i),sigma(i));
    gx(gx<eps*max(X(:,i))) = 0;
    Hn = hermite(opt.der(i),X(:,i) ./ c);
    X(:,i) = Hn .* gx .* (-1/c).^opt.der(i);
end
g = prod(X,2);

Some examples

One dimensional Gaussian kernels in different derivative orders:

leg = {};
for i=0:4
    yi(:,i+1) = ndgauss(33,2*sqrt(2),'der',i);
    leg = [leg sprintf('g%d(x)',i)];
end
figure; plot(yi); legend(leg);

Result:

Two dimensional Gaussian kernels in different derivative orders:

%% partial derivatives 2D

figure('Color','w'); colormap(gray);

% zero derivative
zi = ndgauss([33 33],2*sqrt(2)*[1 1],'der',[0 0]);
subplot(4,4,1); imagesc(zi); axis image; title('g(x,y)');
set(gca,'XTick',[],'YTick',[]);

% dx
zi = ndgauss([33 33],2*sqrt(2)*[1 1],'der',[1 0]);
subplot(4,4,5); imagesc(zi); axis image; title('d/dx g(x,y)');
set(gca,'XTick',[],'YTick',[]);

% dy
zi = ndgauss([33 33],2*sqrt(2)*[1 1],'der',[0 1]);
subplot(4,4,6); imagesc(zi); axis image; title('d/dy g(x,y)');
set(gca,'XTick',[],'YTick',[]);

% dx2
zi = ndgauss([33 33],2*sqrt(2)*[1 1],'der',[2 0]);
subplot(4,4,9); imagesc(zi); axis image; title('d/dx2 g(x,y)');
set(gca,'XTick',[],'YTick',[]);

% dy2
zi = ndgauss([33 33],2*sqrt(2)*[1 1],'der',[0 2]);
subplot(4,4,10); imagesc(zi); axis image; title('d/dy2 g(x,y)');
set(gca,'XTick',[],'YTick',[]);

% dxdy
zi = ndgauss([33 33],2*sqrt(2)*[1 1],'der',[1 1]);
subplot(4,4,11); imagesc(zi); axis image; title('d/dxdy g(x,y)');
set(gca,'XTick',[],'YTick',[]);

% dx3
zi = ndgauss([33 33],2*sqrt(2)*[1 1],'der',[3 0]);
subplot(4,4,13); imagesc(zi); axis image; title('d/dx3 g(x,y)');
set(gca,'XTick',[],'YTick',[]);

% dy3
zi = ndgauss([33 33],2*sqrt(2)*[1 1],'der',[0 3]);
subplot(4,4,14); imagesc(zi); axis image; title('d/dy3 g(x,y)');
set(gca,'XTick',[],'YTick',[]);

% dx2dy
zi = ndgauss([33 33],2*sqrt(2)*[1 1],'der',[2 1]);
subplot(4,4,15); imagesc(zi); axis image; title('d/dx2dy g(x,y)');
set(gca,'XTick',[],'YTick',[]);

% dxdy2
zi = ndgauss([33 33],2*sqrt(2)*[1 1],'der',[1 2]);
subplot(4,4,16); imagesc(zi); axis image; title('d/dxdy2 g(x,y)');
set(gca,'XTick',[],'YTick',[]);

Result:

3D Gaussian kernel

[wi,xi,yi,zi] = ndgauss([33 33 33],2*sqrt(2)*[1 1 1],...
   'der',[0 1 0],'normalize',1);
figure('Renderer','zbuffer','DoubleBuffer','on','Color','w');
contourslice(yi,xi,zi,wi,-5:5,-5:5,-5:5);
camproj perspective;
set(gca,'Color','k');
axis vis3d;

Result:

The Laplacian of Gaussian kernel

can easily be created as follows:

Lg = ndgauss([15 15],[2 2],'der',[2 0]) + ...
     ndgauss([15 15],[2 2],'der',[0 2]);
figure; colormap(gray);
subplot(1,2,1); imagesc(Lg); axis image;
subplot(1,2,2); surf(Lg); axis vis3d;

resulting:

Source & links

  1. Inspiration of the algorithm: Front-End Vision and Multi-Scale Image Analysis book (see Chapter 3 & 4).
  2. Source code at Matlab Central’s File Exchange no. 27763.
Advertisements

3 Responses

Subscribe to comments with RSS.

  1. zunera said, on September 30, 2010 at 3:22 pm

    i want a answr of a question that
    to find the second order diff of 2D gussian function

  2. Essie said, on October 18, 2010 at 9:25 am

    really a great job!!!just wanna thank you because this is exactly what I need!!

  3. CaribeGirl said, on May 9, 2012 at 7:19 pm

    This is a great work!
    just one thing. I think you have the labels of d/dx g(x,y) and d/dy g(x,y) interchanged. d/dx g(x,y) produce vertical edges since you are checking the change in x (horizontal direction). While d/dy g(x,y) produce horizontal edges since you are checking the vertical change. Other than that this is Great!


Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Google+ photo

You are commenting using your Google+ account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s

%d bloggers like this: