function y = filter2(b,x,shape)
%FILTER2 Two-dimensional digital filter.
%   Y = FILTER2(B,X) filters the data in X with the 2-D FIR
%   filter in the matrix B.  The result, Y, is computed 
%   using 2-D correlation and is the same size as X. 
%   
%   Y = FILTER2(B,X,'shape') returns Y computed via 2-D
%   correlation with size specified by 'shape':
%   'same'  - (default) returns the central part of the 
%              correlation that is the same size as X.
%   'valid' - returns only those parts of the correlation
%             that are computed without the zero-padded
%             edges, size(Y) < size(X).
%   'full'  - returns the full 2-D correlation, 
%             size(Y) > size(X).
%   
%   FILTER2 uses CONV2 to do most of the work.  2-D correlation
%   is related to 2-D convolution by a 180 degree rotation of the
%   filter matrix.
%   
%   B and X must be single or double arrays. Both real and complex types are supported.
%   
%   Example:
%    
%       X = gpuArray.rand(3, 4);
%       B = gpuArray.rand(5, 8);
%       Y = filter2(B,X)
%   
%   See also FILTER2, GPUARRAY, GPUARRAY/FILTER, GPUARRAY/CONV2.
%   


%   Copyright 2010-2013 The MathWorks, Inc.
% The logic in this file copies the builtin filter2

narginchk(2,3);
if nargin<3, shape = 'same'; end

if (~isfloat(b)),  b = double(b); end
if (~isfloat(x)),  x = double(x); end

if ~strcmp(shape, 'full') && ...
   ~strcmp(shape, 'same') && ...
   ~strcmp(shape, 'valid')         
    error(message('MATLAB:filter2:InvalidParam'));
end

[mx,nx] = size(x);
stencil = rot90(b,2);
[ms,ns] = size(stencil);

% 1-D stencil?
if (ms == 1)
  y = conv2(1,stencil,x,shape);
elseif (ns == 1)
  y = conv2(stencil,1,x,shape);
else
  if (ms*ns > mx*nx)
    % The filter is bigger than the input.  This is a non-typical
    % case, and it may be counterproductive to check the
    % separability of the stencil.
    y = conv2(x,stencil,shape);
  else
    separable = false;
    
    if all(isfinite(subsref(stencil, substruct('()',{':'}))))
      % Check rank (separability) of stencil: use SVD on CPU.
      [u,s,v] = svd(gather(stencil));
      s = diag(s);
      tol = length(stencil) * eps(max(s));
      rank = sum(s > tol);   
      separable = (rank ==1);
    end
    if separable
      % Separable stencil.
      % Bring SVD products back to GPU.
      u = gpuArray(u);
      s = gpuArray(s);
      v = gpuArray(v);
      idu = substruct('()', {':',1});
      id1 = substruct('()', {1});
      hcol = subsref(u, idu) * sqrt(subsref(s, id1));
      hrow = conj(subsref(v, idu)) * sqrt(subsref(s, id1));
      if (all(all((round(stencil) == stencil))) && all(all((round(x) == x))))
        % Output should be integer
        y = round(conv2(hcol, hrow, x, shape));
      else
        y = conv2(hcol, hrow, x, shape);
      end
    else
      % Non-separable stencil
      y = conv2(x,stencil,shape);
    end
  end
end

end

