%  PSF - 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 = psf(sze, order, ang, eccen, rc, 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
%   eccen - ratio of eccentricity of the filter shape (major/minor axis ratio)
%   rc    - mean radius of the filter in pixels
%   sqrness - even integer specifying `squareness' of the filter shape
%             a value of 2 gives a circular filter (if eccen = 1), higher
%             values make the shape squarer.       

% Copyright (c) 1999 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

function h = psf(sze, order, ang, eccen, rc, sqrness)

    if mod(sqrness,2) ~=0
	error('squareness parameter must be an even integer');
    end
    
    rows = sze(1);
    cols = sze(2);
    
    x = ones(rows,1) * [1:cols]  - (fix(cols/2)+1);
    y = [1:rows]' * ones(1,cols) - (fix(rows/2)+1);
    
    xp = x*cos(ang) - y*sin(ang);   % Rotate axes by specified angle.
    yp = x*sin(ang) + y*cos(ang);
    
    x = sqrt(eccen)*xp;             % Distort x and y according to eccentricity.
    y = yp/sqrt(eccen);
    
    radius = (x.^sqrness + y.^sqrness).^(1/sqrness);  % Distort distance metric
						      % by squareness measure.
    h = 1./(1+(radius./rc).^order);  % Butterworth filter 
    h = h./(sum(sum(h)));            % and normalise.