function Y = linear_filter(X,F,n,b)

% linear_filter - linear filtering
% --------------------------------
% 
% Y = linear_filter(X,F,n,b)
%   = linear_filter(X,F,Z,b)
%
% Input:
% ------
%  X - input image or handle to parent figure
%  F - linear filter
%  n - number of iterations (def: 1)
%  Z - computation mask image (def: [])
%  b - boundary behavior (def: -1)
%
% Output:
% -------
%  Y - linear filtered image

%--------------------------------
% Author: Harold Figueroa
%--------------------------------
% $Date: 2004-12-21 14:49:07 -0500 (Tue, 21 Dec 2004) $
% $Revision: 326 $
%--------------------------------

% TODO: extend mex to compute with multiple filters

% TODO: extend to handle filter decomposition at the mex level

%----------------------------------
% HANDLE INPUT
%----------------------------------

%--
% consider image or handle input
%--

% NOTE: this is very old code and has been put into a separate function

flag = 0;

if ((length(X) == 1) && ishandle(X))
	
    h = X;
	
    hi = get_image_handles(h);
	
    if (length(hi) == 1)
        X = get(hi,'cdata');
        flag = 1;
    else
        warning('Figure has more than one image. No computation performed.');
        return;
	end  
	
end

%--
% set boundary behavior
%--

% NOTE: the boundary behavior code is documented in 'image_pad'

if (nargin < 4)
	b = -1;
end

%--
% iteration or mask
%--

if ((nargin < 3) | isempty(n))
	
	%--
	% set default behavior of single unmasked iteration
	%--
	
	n = 1;
	Z = [];
	
else
	
	%--
	% check whether input may be computation mask rather than iterations
	%--
	
	Z = [];
	
	[r,c,s] = size(X);
	
	if (all(size(n) == [r,c]))
		Z = n;
		n = 1;
	end
	
end

%---------------------------------------
% MULTIPLE PLANE IMAGES
%---------------------------------------

if (ndims(X) > 2)
	
	[r,c,s] = size(X);
	
	for k = 1:s
		if (isempty(Z))
			Y(:,:,k) = linear_filter(X(:,:,k),F,n,b);
		else
			Y(:,:,k) = linear_filter(X(:,:,k),F,Z,b);
		end
	end	
	
%---------------------------------------
% SCALAR IMAGES
%---------------------------------------

else 
	
	%--
	% iterate operator
	%--

	if (isempty(Z))

		%--
		% simple filter computation
		%--

		if (isreal(F))

			pq = (size(F) - 1)./2;

			for j = 1:n
				X = image_pad(X,pq,b);
				Y = linear_filter_(X,F);
				if (n > 1)
					X = Y;
				end
			end

		%--
		% decomposed filter computation
		%--

		else

			for j = 1:n
				Y = linear_filter_decomp(X,F,b);
				if (n > 1)
					X = Y;
				end
			end

		end

	%--
	% mask operator
	%--

	else

		% NOTE: the decomposed filter is not used on masked computation

		if (~issreal(F))
			F = F.H;
		end

		pq = (size(F) - 1)./2;

		X = image_pad(X,pq,b);
		Z = image_pad(Z,pq,0);
		Y = linear_filter_(X,F,uint8(Z));

	end

end

%---------------------------------------
% DISPLAY OUTPUT
%---------------------------------------

% NOTE: this was part of the old image processing framework

% NOTE: this framework permits the processing of figure displaye images

if (flag && view_output)

	switch(view_output)

		%--
		% update display of image from figure
		%--
		
		case (1)
			figure(h);
			set(hi,'cdata',Y);
			set(gca,'clim',fast_min_max(Y));

		%--
		% produce a new image display
		%--
		
		otherwise
			fig;
			image_view(Y);

	end

end


%-----------------------------------------------
% LINEAR_FILTER_DECOMP
%-----------------------------------------------

function Y = linear_filter_decomp(X,F,b)

% linear_filter_decomp - apply decomposed linear filter
% -----------------------------------------------------
%
% Y = linear_filter_decomp(X,F,b)
%
% Input:
% ------
%  X - input image
%  F - decomposed filter structure
%  tol - rank computation tolerance
%  b - boundary behavior

%--------------------------------
% Author: Harold Figueroa
%--------------------------------
% $Date: 2004-12-21 14:49:07 -0500 (Tue, 21 Dec 2004) $
% $Revision: 326 $
%--------------------------------

% TODO: consider whether the approximation needs to be renormalized

%--
% pad input image (only need to do it once)
%--

pq = (size(F.H) - 1)./2;

X = image_pad(X,pq,b);

%--
% apply first component
%--

fx = F.X(:,1)';
fy = F.S(1) * F.Y(:,1);	

Y = linear_filter_(linear_filter_(X,fx),fy);

%--
% apply remaining components if needed
%--

% NOTE: the rank is computed using the tolerance parameter

for k = 2:F.rank	

	fx = F.X(:,k)';
	fy = F.S(k) * F.Y(:,k);	

	Y = Y + linear_filter_(linear_filter_(X,fx),fy);	

end
