%CONVOLVE   General convolution filter
%
% SYNOPSIS:
%  image_out = convolve(image_in,filter)
%
% DEFAULTS:
%  filter = ones(3)/9 (which only works with a 2D image!)
%
% NOTE:
%  The convolution is computed through the Fourier Domain.

% (C) Copyright 1999-2000               Pattern Recognition Group
%     All rights reserved               Faculty of Applied Physics
%                                       Delft University of Technology
%                                       Lorentzweg 1
%                                       2628 CJ Delft
%                                       The Netherlands
%
% Lucas van Vliet, May 2000.
% 24 October 2001 (CL): convolve(a,b) == convolve(b,a)
%                       Made dimensionality-independent.
% Jan 2006, output can be complex now (BR)
% July 2009, for real inputs the output is now forced to be real (BR) 

function out = convolve(varargin)

d = struct('menu','Filters',...
           'display','General convolution',...
           'inparams',struct('name',       {'in',         'flt'},...
                             'description',{'Input image','Filter image'},...
                             'type',       {'image',      'image'},...
                             'dim_check',  {0,            0},...
                             'range_check',{[],           []},...
                             'required',   {1,            0},...
                             'default',    {'a',          'ones(3)/9'}...
                            ),...
           'outparams',struct('name',{'out'},...
                              'description',{'Output image'},...
                              'type',{'image'}...
                              )...
          );
if nargin == 1
   s = varargin{1};
   if ischar(s) & strcmp(s,'DIP_GetParamList')
      out = d;
      return
   end
end
try
   [in,flt] = getparams(d,varargin{:});
catch
   if ~isempty(paramerror)
      error(paramerror)
   else
      error(firsterr)
   end
end

if ndims(in) ~= ndims(flt)
   error('Dimensionalities do not match')
end
insz = size(in);
fltsz = size(flt);
if insz ~= fltsz
   sz = max(insz,fltsz);
   if fltsz ~= sz
      left = floor(sz/2)-floor(fltsz/2);
      right = left+fltsz-1;
      indx = cell(length(left));
      for ii=1:length(left)
         indx{ii} = left(ii):right(ii);
      end
      flt = subsasgn(newim(sz,datatype(flt)),substruct('()',indx),flt);
   end
   if insz ~= sz
      left = floor(sz/2)-floor(insz/2);
      right = left+insz-1;
      indx = cell(length(left));
      for ii=1:length(left)
         indx{ii} = left(ii):right(ii);
      end
      in = subsasgn(newim(sz,datatype(in)),substruct('()',indx),in);
   end
end

out = dip_convolveft(in,flt,'spatial','spatial','spatial');
if isreal(in) & isreal(flt)
   out = real(out); %return real, through away the small round-off erros
end

%??? Crop output image to fit with input image?
