Laplacian Filter Study

Construct a 256 x 256 image of the function gaus(r) = exp(-pi*r^2) over the domain -1.5 < (x,y) < 1.5. Then construct an image of the Laplacian of gaus(r) over the same domain. Find the rms difference between this image and those constructed by applying a Laplacian convolution filter to the original image. Plot the rms difference as a function of alpha. Find the value of alpha that minimizes this difference.

Image of Gaussian function

global  m1 m2 z lz range steps;

range = 3;
steps = 256;
x = linspace(-range/2,range/2,steps);
y = -x;
x = ones(steps,1)*x;
y = y' * ones(1,steps);

r2 = x.*x + y.*y;
z = exp(-pi*r2);
view = imnorm(z);
imshow(view);
imwrite(view,'gaus.jpg');

Image of Laplacian of Gaussian

The image constructed is the negative of the Laplacian.

lz = -4*pi*(pi*r2-1).*exp(-pi*r2);
view = imnorm(lz);
imshow(view);
imwrite(view,'lgaus.jpg');

Plot RMS Error vs. Alpha

m1 = [0 1 0; 1 -4 1; 0 1 0];
m2 = [1 0 1; 0 -4 0; 1 0 1];

clear rms_error
lmin = -.5; lmax = 0.5; lsteps = 31;

for i = 0:(lsteps-1)
   a = lmin+ (lmax-lmin)*i/(lsteps-1);
   rms_error(i+1) = func(a);
end

%Plot the Error Graph
plot(linspace(lmin,lmax,lsteps),rms_error);
title('RMS Error vs. Alpha');
xlabel('Alpha');
ylabel('RMS Error');

Determine the Best Alpha to Minimize RMS Error

%Determine the Best Alpha based upon RMS Error
amin = fminsearch(@func,0)
%Show the Best Laplacian Matrix
lap = ((1-amin).*m1+amin.*m2)/(1+amin)

The results are

amin =

   -0.2174

lap =

   -0.2779    1.5557   -0.2779
    1.5557   -5.1114    1.5557
   -0.2779    1.5557   -0.2779

func.m

function error = func(v)
global  m1 m2 z lz range steps
a = v(1);
lap = ((1-a)*m1+a*m2)/(1+a);
zlap = -filter2(lap,z,'valid')/(range/(steps-1))^2;
diff = zlap-lz;
error = sqrt(mean(diff(:).^2));

imnorm.m

function out = imnorm(inp)
% IMNORM
%	out = imnorm(inp)
%   normalizes input image between (0..1) 
%     for unipolar image divide by the maximum value
%   bipolar image mapped such that zero -> 0.5 and image scales between (0..1)

maxv = max(inp(:));
minv = min(inp(:));
if (minv<0.0)
	vn = max(maxv,-minv);
	out = (inp/vn+1)/2;
else
	out = inp/maxv;
end

Maintained by John Loomis, last updated 14 March 2000