function y = convolve2(x, m, shape, tol)
%CONVOLVE2 Two dimensional convolution.
% Y = CONVOLVE2(X, M) performs the 2-D convolution of matrices X and
% M. If [mx,nx] = size(X) and [mm,nm] = size(M), then size(Y) =
% [mx+mm-1,nx+nm-1]. Values near the boundaries of the output array are
% calculated as if X was surrounded by a border of zero values.
%
% Y = CONVOLVE2(X, M, SHAPE) where SHAPE is a string returns a
% subsection of the 2-D convolution with size specified by SHAPE:
%
% 'full' - (default) returns the full 2-D convolution,
% 'same' - returns the central part of the convolution
% that is the same size as X (using zero padding),
% 'valid' - returns only those parts of the convolution
% that are computed without the zero-padded
% edges, size(Y) = [mx-mm+1,nx-nm+1] when
% size(X) > size(M),
% 'wrap' - as for 'same' except that instead of using
% zero-padding the input X is taken to wrap round as
% on a toroid.
% 'reflect' - as for 'same' except that instead of using
% zero-padding the input X is taken to be reflected
% at its boundaries.
%
% CONVOLVE2 is fastest when mx > mm and nx > nm - i.e. the first
% argument is the input and the second is the mask.
%
% If the rank of the mask M is low, CONVOLVE2 will decompose it into a
% sum of outer product masks, each of which is applied efficiently as
% convolution with a row vector and a column vector, by calling CONV2.
% The function will often be faster than CONV2 or FILTER2 (in some
% cases much faster) and will produce the same results as CONV2 to
% within a small tolerance.
%
% Y = CONVOLVE2(... , TOL) where TOL is a number in the range 0.0 to
% 1.0 computes the convolution using a reduced-rank approximation to
% M, provided this will speed up the computation. TOL limits the
% relative sum-squared error in the effective mask; that is, if the
% effective mask is E, the error is controlled such that
%
% sum(sum( (M-E) .* (M-E) ))
% -------------------------- <= TOL
% sum(sum( M .* M ))
%
% See also CONV2, FILTER2.
% Copyright David Young, Feb 2002, revised Jan 2005, Jan 2009, Apr 2011
% Deal with optional arguments
error(nargchk(2,4,nargin));
if nargin < 3
shape = 'full'; % shape default as for CONV2
tol = 0;
elseif nargin < 4
if isnumeric(shape)
tol = shape;
shape = 'full';
else
tol = 0;
end
end;
% Set up to do the wrap & reflect operations, not handled by conv2
if ismember(shape, {'wrap' 'reflect'})
x = extendarr(x, m, shape);
shape = 'valid';
end
% do the convolution itself
y = doconv(x, m, shape, tol);
end
%-----------------------------------------------------------------------
function y = doconv(x, m, shape, tol)
% Carry out convolution
[mx, nx] = size(x);
[mm, nm] = size(m);
% If the mask is bigger than the input, or it is 1-D already,
% just let CONV2 handle it.
if mm > mx || nm > nx || mm == 1 || nm == 1
y = conv2(x, m, shape);
else
% Get svd of mask
if mm < nm; m = m'; end % svd(..,0) wants m > n
[u,s,v] = svd(m, 0);
s = diag(s);
rank = trank(m, s, tol);
if rank*(mm+nm) < mm*nm % take advantage of low rank
if mm < nm; t = u; u = v; v = t; end % reverse earlier transpose
vp = v';
% For some reason, CONV2(H,C,X) is very slow, so use the normal call
y = conv2(conv2(x, u(:,1)*s(1), shape), vp(1,:), shape);
for r = 2:rank
y = y + conv2(conv2(x, u(:,r)*s(r), shape), vp(r,:), shape);
end
else
if mm < nm; m = m'; end % reverse earlier transpose
y = conv2(x, m, shape);
end
end
end
%-----------------------------------------------------------------------
function r = trank(m, s, tol)
% Approximate rank function - returns rank of matrix that fits given
% matrix to within given relative rms error. Expects original matrix
% and vector of singular values.
if tol < 0 || tol > 1
error('Tolerance must be in range 0 to 1');
end
if tol == 0 % return estimate of actual rank
tol = length(m) * max(s) * eps;
r = sum(s > tol);
else
ss = s .* s;
t = (1 - tol) * sum(ss);
r = 0;
sm = 0;
while sm < t
r = r + 1;
sm = sm + ss(r);
end
end
end
%-----------------------------------------------------------------------
function y = extendarr(x, m, shape)
% Extend x so as to wrap around on both axes, sufficient to allow a
% "valid" convolution with m to return the cyclical convolution.
% We assume mask origin near centre of mask for compatibility with
% "same" option.
[mx, nx] = size(x);
[mm, nm] = size(m);
mo = floor((1+mm)/2); no = floor((1+nm)/2); % reflected mask origin
ml = mo-1; nl = no-1; % mask left/above origin
mr = mm-mo; nr = nm-no; % mask right/below origin
if strcmp(shape, 'wrap')
y = exindex(x, 1-ml:mx+mr, 1-nl:nx+nr, 'circular');
else
y = exindex(x, 1-ml:mx+mr, 1-nl:nx+nr, 'symmetric');
end
end