% rotation matrix
%
% Image rotation, in 2d.
%
% A = matRot2(sz,angle,center,method)
%  angle  from [0,2*pi)
%  center is a 2 vector                                       [default (sz+1)/2]
%  method is the interpolation method nearest,linear,cubic,shear   [def. linear]
%
%  shear uses the method by Paeth that does lossless rotation by three
%  consecutive shearing operations
%
% (c) by Hannes Nickisch, MPI for Biological Cybernetics, 2012 August 15

function A = matRot2(sz, angle, center, method, complex)

  A.imsz = sz; sz = prod(A.imsz)*[1,1];
  if nargin<2, angle = 0; end
  A.angle = angle;                                       % min angle is required
  if nargin<3, center = (A.imsz+1)/2; end         % default as Matlab's imrotate
  if numel(center)<2, center = (A.imsz+1)/2; end, A.center  = center;
  if nargin<4, method = 'linear'; end, A.method = method;
  if nargin<5, complex = 1; end
  if numel(complex)==1, complex = complex*[1,1]; end

  A.type = 'mat';
  A.w1=[]; A.d1=[]; A.w2=[]; A.d2=[]; A.g=[]; A.x=[]; A.W=[];  % make vars known
  if strcmp(A.method,'shear')
    ramp = @(n) ((1:n)-fix(1+n/2))/n;
    s = reshape(ramp(A.imsz(1))*A.imsz(1),[A.imsz(1),1]);
    r = reshape(ramp(A.imsz(2)),[1,A.imsz(2)]); w = -2i*pi*s*r;
    A.w1 = exp(tan(angle/2)*w);
    A.d1 = 0.5/cos(angle/2)^2*w .* A.w1;
    s = reshape(ramp(A.imsz(2))*A.imsz(2),[1,A.imsz(2)]);
    r = reshape(ramp(A.imsz(1)),[A.imsz(1),1]); w = -2i*pi*r*s;
    A.w2 = exp(-sin(angle)*w); A.d2 = -cos(angle)*w .* A.w2;
  else
    R = rotmat(angle(1));                                      % rotation matrix
    [i2,i1]=meshgrid(1:A.imsz(2),1:A.imsz(1)); i1=i1(:)'; i2=i2(:)'; % targ grid
    x   = R'*[i1-center(1)    ; i2-center(2)    ];        % non-integer src grid
    A.x = [x(1,:)+center(1); x(2,:)+center(2)]; A.g = [i1; i2];
    A.W = matResample(A.imsz,A.x,A.method,A.type); % interpolation weight matrix
  end
  
  A.di = 0;                                                    % derivative flag
  A = class(A,mfilename,mat(sz,complex));
