% rotation matrix
%
% Image rotation, in 3d.
%
% A = matRot3(sz,rot,center,method)
%  rot is a 3 vector specifying the rotation
%      - n  = rot/norm(rot) is the (unit) axis of rotation
%      - al = norm(rot)     is the angle
%  center is a 3 vector                                       [default (sz+1)/2]
%  method is the interpolation method nearest, linear, cubic   [default nearest]
%
% matRot3 contains matRot2 as a special case
%   matRot2(sz,angle,center,method) = ...
%                               matRot3([sz,1],angle*[0,0,1],[center,1],method);
%
% (c) by Hannes Nickisch, MPI for Biological Cybernetics, 2011 February 17

function A = matRot3(sz, rot, center, method, complex)

  A.imsz = sz; sz = prod(A.imsz)*[1,1];
  if nargin<2, rot = zeros(3,1); end, A.rot = rot(:);
  if nargin<3, center = (A.imsz+1)/2; end         % default as Matlab's imrotate
  if numel(center)<3, 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  

  R = rotmat(A.rot);                                           % rotation matrix

  [i2,i1,i3] = meshgrid(1:A.imsz(2),1:A.imsz(1),1:A.imsz(3));
  i1=i1(:)'; i2=i2(:)'; i3=i3(:)';  A.g = [i1; i2; i3];            % target grid
  x = R'*[i1-center(1)    ; i2-center(2)    ; i3-center(3) ];      % non-integer
  A.x = [x(1,:)+center(1); x(2,:)+center(2); x(3,:)+center(3)]; % .. source grid

  A.type = 'mat';
  A.W = matResample(A.imsz,A.x,A.method,A.type);   % interpolation weight matrix
  A.di = 0;                                                    % derivative flag

  A = class(A,mfilename,mat(sz,complex));
