%% Copyright (C) 2005 Peter Kovesi
%%
%% This program is free software; you can redistribute it and/or modify
%% it under the terms of the GNU General Public License as published by
%% the Free Software Foundation; either version 2, or (at your option)
%% any later version.
%%
%% This program is distributed in the hope that it will be useful, but
%% WITHOUT ANY WARRANTY; without even the implied warranty of
%% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
%% General Public License for more details.
%%
%% You should have received a copy of the GNU General Public License
%% along with this program. If not, write to the Free Software
%% Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
%% 02110-1301, USA.

%% FSPECIAL - Create spatial filters for image processing
%%
%% Usage:  f = fspecial(filtertype, optionalparams)
%%
%% filtertype can be
%%
%%    'average'   - Rectangular averaging filter
%%    'disk'      - Circular averaging filter.
%%    'gaussian'  - Gaussian filter.
%%    'log'       - Laplacian of Gaussian.
%%    'laplacian' - Crude 3x3 spproximation to laplacian
%%    'unsharp'   - Crude sharpening filter
%%    'motion'    - Motion blur filter of width 1 pixel
%%    'sobel'     - Horizontal Sobel edge filter
%%    'prewitt'   - Horizontal Prewitt edge filter
%%
%% The parameters that need to be specified depend on the filtertype
%%
%% Examples of use and associated default values:
%%
%%  f = fspecial('average',sze)           sze can be a 1 or 2 vector
%%                                        default is [3 3].
%%  f = fspecial('disk',radius)           default radius = 5
%%  f = fspecial('gaussian',sze, sigma)   default sigma is 0.5
%%  f = fspecial('laplacian',sze, sigma)  default sze is [5 5]
%%                                        default sigma is 0.5
%%  f = fspecial('log');
%%  f = fspecial('motion', length, angle) default length is 9
%%                                        default angle is 0 (degrees)
%%  f = fspecial('sobel');
%%  f = fspecial('prewitt');
%%
%% Where sze is specified as a single value the filter will be square.

%% Author:   Peter Kovesi <pk@csse.uwa.edu.au>
%% Keywords: image processing, spatial filters
%% Created:  August 2005


function f = my_fspecial(varargin)
    
    [filtertype, sze, sigma, radius, len, angle] = checkargs(varargin(:));    
    
    rows = sze(1); cols = sze(2);
    r2 = (rows-1)/2; c2 = (cols-1)/2;

    if strcmpi(filtertype, 'average')
	f = ones(sze)/(rows*cols);
    
    elseif strcmpi(filtertype, 'disk')
	[x,y] = meshgrid(-c2:c2, -r2:r2);
	rad = sqrt(x.^2 + y.^2);	
	f = rad <= radius;
	f = f/sum(f(:));	
	
    elseif strcmpi(filtertype, 'gaussian')
	[x,y] = meshgrid(-c2:c2, -r2:r2);
	radsqrd = x.^2 + y.^2;
	f = exp(-radsqrd/(2*sigma^2));
	f = f/sum(f(:));
	
    elseif strcmpi(filtertype, 'log') 
	[x,y] = meshgrid(-c2:c2, -r2:r2);
	radsqrd = x.^2 + y.^2;	
	f = -1/(pi*sigma^4)*(1-radsqrd/(2*sigma^2))...
	     .*exp(-radsqrd/(2*sigma^2));	   
	f = f-mean(f(:));    % Ensure 0 DC

    elseif strcmpi(filtertype, 'laplacian')	    	
	f = [1  1  1 
	     1 -8  1 
	     1  1  1]; 
    
    elseif strcmpi(filtertype, 'unsharp')	    		
	f = -fspecial('log') + [0 0 0 
		                0 1 0 
		                0 0 0];
	
    elseif strcmpi(filtertype, 'sobel')
	f = [1 2 1; 0 0 0; -1 -2 -1];
    
    elseif strcmpi(filtertype, 'prewitt')
	f = [1 1 1; 0 0 0; -1 -1 -1];	
	
    elseif strcmpi(filtertype, 'motion')	
	% First generate a horizontal line across the middle
	f = zeros(sze);
	f(floor(len/2)+1,1:len) = 1;   

	% Then rotate to specified angle
	f = imrotate(f,angle,'bilinear','loose');
	f = f/sum(f(:));
	
    else
	error('Unrecognized filter type');
	
    end	
	
%%----------------------------------------------------------------	
%%
%% Sort out input arguments, setting defaults and checking for errors.
	
function [filtertype, sze, sigma, radius, len, angle] = checkargs(arg)
    
% Set defaults
    
    sze = [3 3];
    sigma = 0.5;
    radius = 5;
    len = 9;
    angle = 0;

    narg = length(arg);    

    filtertype = arg{1};
    if ~ischar(filtertype)
	error('filtertype must be specified as a string');
    end

    if strcmpi(filtertype, 'log')    
      if narg == 1
	  sze = [5 5];
      end
    end
    
    if strcmpi(filtertype, 'average')   | ...
       strcmpi(filtertype, 'gaussian')  | ...
       strcmpi(filtertype, 'log')
	if narg >= 2
	    sze = arg{2};
	    if isscalar(sze)
		sze = [sze sze];
	    end
	end
    end

    if strcmpi(filtertype, 'gaussian')  | ...
       strcmpi(filtertype, 'log')    
	if narg >= 3
	    sigma = arg{3};
	end
    end
    
    if strcmpi(filtertype, 'disk')  
	if narg >= 2
	    radius = arg{2};
	end
	sze = [2*radius+1 2*radius+1];
    end
    
    if strcmpi(filtertype, 'motion')      
	if narg >= 2
	    len = arg{2};
	end

	if narg >= 3
	    angle = arg{3};
	end
	
	% Ensure size is odd so that there is a middle point 
	% about which to rotate the filter by angle.
	if mod(len,2)         
	    sze = [len len];
	else
	    sze = [len+1 len+1];
	end
    end
    
    if ~isscalar(len) | len < 1
	error('length must be a scalar >= 1');
    end
    
    if ~isscalar(angle)
	error('angle must be a scalar');
    end

    if ~isscalar(radius) | radius < 1
	error('radius must be a scalar >= 1');
    end
    
    if ~isscalar(sigma) | sigma < 0.5
	error('sigma must be a scalar >= 0.5');
    end    
    
    if length(sze) > 2
	error('filter size must be a scalar or 2-vector');
    end
    
    if any(fix(sze) ~= sze)
	error('filter size must be integer');
    end	