% Resampling matrix: explicit sparse matrix
%
% Signal resampling in 2d and 3d.
%
% method needs to contain {lin, cub, epa<l>, biw<l>, trc<l>, *} to match
%  {linear, cubic, Epanechnikov(l), biweight(l), tricubic(l), nearest}
%
% (c) by Hannes Nickisch, MPI for Biological Cybernetics, 2011 February 24

function W = resample_mat(sz,x,method,deriv)

  ndims = numel(sz); N = prod(sz); M = size(x,2);                   % dimensions

  if nargin<3, method = 'lin'; end                          % set default values
  if nargin<4, deriv  = 0; end
  if     strfind(lower(method),'lin')                     % set window functions
    l = 1; k = @(x,d) lin(x,d);
  elseif strfind(lower(method),'cub')
    l = 2; k = @(x,d) cub(x,d);
  elseif strfind(lower(method),'epa')
    if numel(method)==3, l = 2; else l = str2double(method(4)); end
    if numel(l)==0; l = 2; end, k = @(x,d) epa(x,d,l);
  elseif strfind(lower(method),'biw')
    if numel(method)==3, l = 2; else l = str2double(method(4)); end
    if numel(l)==0; l = 2; end, k = @(x,d) biw(x,d,l);
  elseif strfind(lower(method),'trc')
    if numel(method)==3, l = 2; else l = str2double(method(4)); end
    if numel(l)==0; l = 2; end, k = @(x,d) trc(x,d,l);
  else
    l = 1; k = @(x,d) nea(x,d);
  end

  W = sparse(M,N);                                 % interpolation weight matrix
  for n1 = 1-l:l,   j1 = floor(x(1,:))+n1;    % linear interpolation from source
    for n2 = 1-l:l, j2 = floor(x(2,:))+n2;   % .. neighbours; iterate over those
      if ndims==2
        ok = 0<j1 & j1<=sz(1)  &  0<j2 & j2<=sz(2);     %   indices inside image
        j = j1(ok) + (j2(ok)-1)*sz(1);
        w =    k(x(1,ok)-j1(ok),deriv==1);
        w = w.*k(x(2,ok)-j2(ok),deriv==2);
        W = W + sparse(find(ok),j,w,M,N);                % add upp contributions
      else
        for n3=1-l:l, j3 = floor(x(3,:))+n3;
          ok = 0<j1 & j1<=sz(1) & 0<j2 & j2<=sz(2) ...
                                & 0<j3 & j3<=sz(3);     %   indices inside image
          j = j1(ok) + (j2(ok)-1)*sz(1) + (j3(ok)-1)*prod(sz(1:2));
          w =    k(x(1,ok)-j1(ok),deriv==1);
          w = w.*k(x(2,ok)-j2(ok),deriv==2);
          w = w.*k(x(3,ok)-j3(ok),deriv==3);
          W = W + sparse(find(ok),j,w,M,N);              % add upp contributions
        end
      end
    end
  end


% Robert G. Keys, Cubic Convolution Interpolation for Digital Image Processing,
% IEEE ASSP, 29:6, December 1981, p. 1153-1160.
function y = cub(x,d)
  if ~d
    y = zeros(size(x)); x = abs(x);
    q = (x<=1);          % Coefficients:  1.5, -2.5,  0.0, 1.0
    y(q) =            (( 1.5 * x(q) - 2.5) .* x(q)      ) .* x(q) + 1.0;
    q = (1<x & x<=2);    % Coefficients: -0.5,  2.5, -4.0, 2.0
    y(q) =            ((-0.5 * x(q) + 2.5) .* x(q) - 4.0) .* x(q) + 2.0;
  else
    y = sign(x); x = abs(x);
    q = (x<=1);          % Coefficients:  1.5, -2.5,  0.0, 1.0
    y(q) = y(q) .*  ( 4.5 * x(q) - 5.0) .* x(q);
    q = (1<x & x<=2);    % Coefficients: -0.5,  2.5, -4.0, 2.0
    y(q) = y(q) .* ((-1.5 * x(q) + 5.0) .* x(q) - 4.0);
  end
  y(x>2) = 0;

function y = lin(x,d)                                           % linear kernel
  if ~d
    y = max(1-abs(x),0);
  else
    y = -sign(x); y(abs(x)>1) = 0;
  end

function y = nea(x,d)                           % nearest neighbor (box) kernel
  y = zeros(size(x));
  if ~d
    y(abs(x)<=1/2) = 1;
  end

function y = epa(x,d,l)       % Epanechnikov kernel, variance=l/5, support -l..l
  if ~d
    y = ( 3/4/l)  *(1-x.*x/l^2);
  else
    y = (-3/2/l^3)*      x     ;
  end

function y = biw(x,d,l)           % biweight kernel, variance=l/7, support -l..l
  y = (1-x.*x/l^2);
  if ~d
    y = ( 15/16/l  ) * y.*y;
  else
    y = (-15/ 4/l^3) * x.*y;
  end

function y = trc(x,d,l)            % tricube kernel, variance=l/9, support -l..l
  y = (1-x.*x/l^2);
  if ~d
    y = (  35/32/l  ) * y.*y.*y;
  else
    y = (-105/16/l^3) * x.*y.*y;
  end