% Two and threedimensional (centered) Fourier transformation
% where the input can be subject to rigid body motion
%
% A = matGPUFFTmotion(sz,trans,rot,complex,type)
%
% sz: size of the array; either 2d or 3d
%
% 2d i.e. numel(sz)=2
%  trans: [2,1] or [2,sz(2)], the translation in pixels,           default [0;0]
%  rot:   [1,1] or [1,sz(2)], the rotation angle in radian,            default 0
%
% 3d i.e. numel(sz)=3
%  trans: [3,1] or [3,sz(2),sz(3)], the translation in pixels,   default [0;0;0]
%  rot:   [3,1] or [3,sz(2),sz(3)], the rotation axis, whose norm equals the 
%                                    rotation angle in radian,   default [0;0;0]
% Note that by convention, translation is done after rotation.
% Note further that the matrix relies on the fact that the pixels/voxels are
% cubes i.e. we have the same resolution along every dimension.
%
% complex: [1,2] complex number treatment A:IN^n->OUT^m,           default [1,1]
%                 first number specifies output field OUT^m and the second the
%                 input field IN^m of the matrix A
%                 1 -  real(z) + 1i*imag(z)    IC:   native complex numbers
%                 2 - [real(z); imag(z)]       IR^2: stacked real numbers
%                 3    real(z)                 IR:   real part only
%
%
% (c) Alexander Loktyushin,     MPI for Biological Cybernetics, 2011

function A = matGPUFFTmotion(sz, trans, rot, complex, type, pre_op)

  %% general info
  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.'), end
  A.imsz = sz; sz = prod(A.imsz)*[1,1];
  if nargin<4, complex = 1; end, if numel(complex)==0, complex = 1; end
  if numel(complex)==1, complex = complex*[1,1]; end
  center = ceil((A.imsz-1)/2)+1; A.center = center;          % center in k-space
  if nargin<5, type = 'linear'; end, A.type = type;   % set (default) nufft type
  if nargin<6, pre_op = []; end;
  
  %% derivatives
  A.dkind = 'none';                    % kind of derivative, either rot or trans
  A.di = 0;               % which derivative, 0 means all of the respective kind

  %% rotation
  nrot = A.ndims*(A.ndims-1)/2;                % number of rotational parameters
  if nargin<3,        rot = zeros(nrot,1); end
  if numel(rot)<nrot, rot = zeros(nrot,1); end, A.rot = rot(:);
    if numel(A.rot)>nrot                                 % bring into proper shape
    A.rot = reshape(A.rot,[nrot,prod(A.imsz(2:end))]);
  end
  
  nusz = A.imsz;
  if numel(nusz)==3 && nusz(3)==1, nusz = nusz(1:2); k = k(:,1:2); end % pseudo 3d image
  A.nuimsz = nusz;
  A.nundims = numel(nusz);                            % Are we doing 2d or 3d stuff?

  %% translation
  if nargin<2,             trans = zeros(A.ndims,1); end
  if numel(trans)<A.ndims, trans = zeros(A.ndims,1); end, A.trans = trans(:);
  if numel(A.trans)>A.ndims                            % bring into proper shape
    A.trans = reshape(A.trans,[A.ndims,A.imsz(2:end)]);
  end
  
  %keyboard
  
    
  %% construct matrix class, inherit from mat
  A = class(A,mfilename,mat(sz,complex));




