Friday, May 25. 2007
Dilated + rotated 2D Gaussian
I’ve finally improved my script that computes a dilated and rotated 2D Gaussian with respect to low memory consumption. It applies a dilation and rotation matrix to the 2D domain. I published it under the GNU GPLv2, so if you find any improvements, you are forced to share them. 
I still have the impression that it is rather slow, but I currently don’t know how to do it better. Of course, it has to be evaluated on
function g=nsgauss(p,q,vdil,hdil,rot)
% Computes a non-separable (dilated + rotated) 2D Gaussian
% Usage: g = nsgauss(p, q, vdil, hdil, rot);
% Input: p,q .... size of g
% vdil ... vertical dilation factor (before rotation)
% hdil ... horizontal dilation factor (before rotation)
% rot .... rotation angle, e.g. pi/4
% Example:
% norm(nsgauss(p,q,1,1,0) - gaussnk(p)’*gaussnk(q)) == eps
%
% Version 0.2-20070525
% by Stephan Paukner {stephan+math at paukner dot cc}
% Licensed under the GNU General Public License v2
% $Id: nsgauss.m,v 1.2 2007/05/25 10:14:52 ps Exp ps $
D=[1/vdil 0; 0 1/hdil]; %dilation matrix
R=[cos(rot) -sin(rot); sin(rot) cos(rot)]’; %’%rotation matrix
sp=sqrt(p); sq=sqrt(q);
g=zeros(1,p*q);
for jp=-3:3
for jq=-3:3
[x y]=meshgrid( (0:p-1)/sp + jp*sp , (0:q-1)/sq + jq*sq );
v=D*R*[x(:)’; y(:)’];
g=g+exp(-pi*(v(1,:).^2 + v(2,:).^2));
end
end
g=reshape(g,q,p)’; %’
g=g/norm(g,’fro’);
% Computes a non-separable (dilated + rotated) 2D Gaussian
% Usage: g = nsgauss(p, q, vdil, hdil, rot);
% Input: p,q .... size of g
% vdil ... vertical dilation factor (before rotation)
% hdil ... horizontal dilation factor (before rotation)
% rot .... rotation angle, e.g. pi/4
% Example:
% norm(nsgauss(p,q,1,1,0) - gaussnk(p)’*gaussnk(q)) == eps
%
% Version 0.2-20070525
% by Stephan Paukner {stephan+math at paukner dot cc}
% Licensed under the GNU General Public License v2
% $Id: nsgauss.m,v 1.2 2007/05/25 10:14:52 ps Exp ps $
D=[1/vdil 0; 0 1/hdil]; %dilation matrix
R=[cos(rot) -sin(rot); sin(rot) cos(rot)]’; %’%rotation matrix
sp=sqrt(p); sq=sqrt(q);
g=zeros(1,p*q);
for jp=-3:3
for jq=-3:3
[x y]=meshgrid( (0:p-1)/sp + jp*sp , (0:q-1)/sq + jq*sq );
v=D*R*[x(:)’; y(:)’];
g=g+exp(-pi*(v(1,:).^2 + v(2,:).^2));
end
end
g=reshape(g,q,p)’; %’
g=g/norm(g,’fro’);
Here’s a comparison of computing time and accuracy:
> tic; g1=gaussnk(600)’ * gaussnk(800); toc
Elapsed time is 0.100037 seconds.
> tic; g2=nsgauss(600,800,1,1,0); toc
Elapsed time is 19.163170 seconds.
> compnorm(g1,g2);
quotient of norms: norm(x)/norm(y) = 1
difference of normalized versions = 1.344e-16
And some plots using different parameters:
 |  |
 |  |
 |  |
Categories
Follow me
