% Resampling matrix
%
% Signal resampling in 2d and 3d.
%
% A = matResample(sz,x,g,method,type)
%  sz     size of the image
%  x      new coordinates [d,m]
%  method is the interpolation method: nearest, linear, cubic   [default linear]
%  type   is the implementation type: mex, mat                     [default mat]
%
% The type='mat' allows in addition for method='epa<l>', 'biw<l>', 'trc<l>',
% which corresponds to Epanechnikov, biweight and tricubic kernel of width 2*l.
%
% (c) by Hannes Nickisch & Alexander Loktyushin, 
%                               MPI for Biological Cybernetics, 2011 February 24

% cd code
% mex resample_mex.c CFLAGS="\$CFLAGS -fopenmp" LDFLAGS="\$LDFLAGS -fopenmp"
% cd ..

function  A = matResample(sz,x,method,type,complex)

  A.imsz = sz;
  A.ndims = numel(sz);                            % Are we doing 2d or 3d stuff?
  if A.ndims<2 || A.ndims>3, error('We support only 2d and 3d data.'), end
  A.x = x;
  sz = [size(x,2),prod(A.imsz)];  
  if nargin<2, error('We need at least two input arguments.'), end
  if nargin<3, A.method = 'linear'; else A.method = method; end    % set default
  if nargin<4, A.type   = 'mat';    else A.type   = type;   end    % set default

  A.di = 0; A.W = []; A.dW = cell(A.ndims,1);           % which derivative, init
  if any(strfind(lower(A.type),'mex'))
    setenv('OMP_NUM_THREADS', '1');                         % multi core support
    % determine #cores automatically :
    % setenv('OMP_NUM_THREADS', num2str(feature('numCores')));  
  else
    A.W = resample_mat(A.imsz,A.x,A.method);          % precompute sparse matrix
  end

  if nargin<5, complex = 1; end
  if numel(complex)==1, complex = complex*[1,1]; end
  A = class(A,mfilename,mat(sz,complex));