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