% Two and threedimensional (centered) Fourier transformation
% where the input can be subject to rigid body motion:
%    matrix vector multiplication
%
% (c) by Hannes Nickisch & Alexander Loktyushin, 
%                               MPI for Biological Cybernetics, 2011

function y = mvm(A,x,ctransp)

  if strcmp(A.dkind,'none')
    if ctransp                                                             % Transpose Slow CPU supplement
      tr_t = trans_diag(A.imsz,A.trans);                                   % Linear phase ramp multipliers for translations
      
      psz = fix(A.nuimsz/2);
      tr_F = matFFTN(2*A.nuimsz)*matZpad(A.nuimsz,psz,A.nuimsz-psz);       % Padd to double size + FFT
      
      nrot = A.ndims*(A.ndims-1)/2;                                        % number of rotational parameters
      dkR = cell(1,nrot);                                                  % derivatives of the rotation matrices
      [tr_k, dkR{:}] = grid_create(A.rot, A.imsz);                         % create grid
      
      y = 2*(conj(tr_t(:)).*x(:));
      y = resample_mex(tr_k(:),y(:),2*A.nuimsz,ctransp,0,2);               % Resample
      y = tr_F'*y;
    else 
      y = resample_gpu(A.trans,A.rot,x,A.nuimsz,0);                        % Fast GPU FFT motion
    end;
  elseif strcmp(A.dkind,'all')                                             % Compute image and derivatives in one shot
      if ctransp             
        error('No transposes of derivatives');
      end;      
      
      if A.ndims == 3
        y = zeros(prod(A.imsz),7);
        
        if strcmp(A.type, 'cubicgpumem')                                   % 30% slower, but memory efficient
          [y(:,1) y(:,2) y(:,3) y(:,4) y(:,5) y(:,6) y(:,7)] = resample_gpu_memory_efficient(A.trans,A.rot,x(:),A.nuimsz,1);
        elseif strcmp(A.type, 'cubicgpudouble')                            % double precision - two times slower
          [y(:,1) y(:,2) y(:,3) y(:,4) y(:,5) y(:,6) y(:,7)] = resample_gpu_double(A.trans,A.rot,x(:),A.nuimsz,1);
        else
          [y(:,1) y(:,2) y(:,3) y(:,4) y(:,5) y(:,6) y(:,7)] = resample_gpu(A.trans,A.rot,x(:),A.nuimsz,1);
        end
      else
        y = zeros(prod(A.imsz),4);
        if strcmp(A.type, 'cubicgpumem')                                   % 30% slower, but memory efficient
          [y(:,1) y(:,2) y(:,3) y(:,4)] = resample_gpu_memory_efficient(A.trans,A.rot,x(:),A.nuimsz,1);
        elseif strcmp(A.type, 'cubicgpudouble')                            % double precision - two times slower
            [y(:,1) y(:,2) y(:,3) y(:,4)] = resample_gpu_double(A.trans,A.rot,x(:),A.nuimsz,1);
        else
          [y(:,1) y(:,2) y(:,3) y(:,4)] = resample_gpu(A.trans,A.rot,x(:),A.nuimsz,1);
        end
      end;
      
  else
    error('No selective derivatives in fast mode');
  end
  
  
function [t,dt] = trans_diag(sz,tr)

  ndims = numel(sz); N = prod(sz);      % number of dimensions, number of pixels
  deriv = nargout>1;               % flag saying, whether we compute derivatives
  t = ones(sz); if deriv, dt = ones(N,ndims); end              % allocate memory
  for d=1:ndims
    rp = ((1:sz(d))-fix(1+sz(d)/2))/sz(d);    % ramp in direction of dimension d
    rp = -2i*pi*reshape(rp,[ones(1,d-1),sz(d),ones(1,ndims-d)]);
    erp = exp(bsxfun(@times,rp,tr(d,:,:)));           % tr is singleton in dim 1
    t = bsxfun(@times, t, erp);
    if deriv
      for dd=1:ndims
        dtdd = bsxfun(@times, reshape(dt(:,dd),sz), erp);
        if dd==d
          dtdd = bsxfun(@times, dtdd, rp);
        end
        dt(:,dd) = reshape(dtdd,N,1);
      end
    end
  end  