% PSF2 - Generates point spread functions for use with deconvolution fns.
%
% This function can generate a variety function shapes based around the
% Butterworth filter. In plan view the filter can be elliptical and at
% any orientation. The 'squareness/roundness' of the shape can also be
% manipulated.
%
% Usage: h = psf2(sze, order, ang, lngth, width, sqrness)
%
% sze - two element array specifying size of filter [rows cols]
% order - an even integer specifying the order of the Butterworth filter.
% This controls the sharpness of the cutoff.
% ang - angle of rotation of the filter in radians.
% lngth - length of the filter in pixels along its major axis.
% width - width of the filter in pixels along its minor axis.
% sqrness - even integer specifying 'squareness' of the filter shape
% a value of 2 gives a circular filter (if lngth == width), higher
% values make the shape squarer.
%
% This function is almost identical to psf, it just has a different way of
% specifying the function shape whereby length and width are defined
% explicitly (rather than an average radius), this may be more convenient for
% some applications.
% Copyright (c) 1999-2003 Peter Kovesi
% School of Computer Science & Software Engineering
% The University of Western Australia
% http://www.csse.uwa.edu.au/
%
% Permission is hereby granted, free of charge, to any person obtaining a copy
% of this software and associated documentation files (the "Software"), to deal
% in the Software without restriction, subject to the following conditions:
%
% The above copyright notice and this permission notice shall be included in
% all copies or substantial portions of the Software.
%
% The Software is provided "as is", without warranty of any kind.
% June 1999
% May 2003 - Changed arguments so that psf is specified in terms of a length
% and width rather than an average radius.
function h = psf2(sze, order, ang, lngth, width, sqrness)
if mod(sqrness,2) ~=0
error('squareness parameter must be an even integer');
end
rows = sze(1);
cols = sze(2);
[x,y] = meshgrid([1:cols],[1:rows]);
% The following fiddles the origin to the correct position
% depending on whether we have and even or odd size
if mod(cols,2) == 0
x = x-cols/2-1;
else
x = x-(cols+1)/2;
end
if mod(rows,2) == 0
y = y-rows/2-1;
else
y = y-(rows+1)/2;
end
xp = x*cos(ang) - y*sin(ang); % Rotate axes by specified angle.
yp = x*sin(ang) + y*cos(ang);
rc = lngth/2; % Set cutoff radius to half the length
yp = yp*lngth/width; % Adjust y measure to give appropriate relative width.
radius = (xp.^sqrness + yp.^sqrness).^(1/sqrness); % Distort distance metric
% by squareness measure.
h = 1./(1+(radius./rc).^order); % Butterworth filter
h = h./(sum(sum(h))); % and normalise.