% Compute a sparse interpolation matrix D of size (nt,ns) from a source grid s
% of size ns to a target array t of size nt.
% The interpolation matrix D has rows summing up to 1.
%
% Copyright (c) by Hannes Nickisch 2014-11-24.

function D = intermat(s,t,d)
  if nargin<3, d = 3; end                                    % set default value
  switch d
    case 0, k = @(x) abs(x)<0.5;      kI = -1:0;
    case 1, k = @(x) max(1-abs(x),0); kI = -1:0;
    case 3, k = @(x) kcub(x);         kI = -2:1;
    otherwise, error('Unknown interpolation kernel.')
  end

  nt = numel(t); ns = numel(s); ds = s(2)-s(1);
  j = floor((t-s(1))/ds)+1;                % index of closest smaller grid point
  w = (t-s(1))/ds-j+1;   % relative distance to closest smaller grid point [0,1]
  D = 0; for i=kI, D = D + sparse(1:nt,j-i,k(w+i),nt,ns); end

% Robert G. Keys, Cubic Convolution Interpolation for Digital Image Processing,
% IEEE ASSP, 29:6, December 1981, p. 1153-1160.
function y = kcub(x)
  y = zeros(size(x)); x = abs(x);
  q = x<=1;          % Coefficients:  1.5, -2.5,  0, 1
  y(q) =            (( 1.5 * x(q) - 2.5) .* x(q)    ) .* x(q) + 1;
  q = 1<x & x<=2;    % Coefficients: -0.5,  2.5, -4, 2
  y(q) =            ((-0.5 * x(q) + 2.5) .* x(q) - 4) .* x(q) + 2;