%ROTATION3D   Rotate a 3D image
%
%  There are 2 different methods avaiable to do the rotation
%  by DIRECT mapping or by 9 SHEARS. Both have different advantages.
%
%  Rotation in 3D by 3 Euler angles
%  R= R_{3''}(\gamma) R_{2'}(\beta) R}_3(\alpha)
%   Note: There is another definition of the
%   Euler angle. The second rotation is then about the intermediate 1-axis.
%   This convention is traditionally used in the mechanics of the rigid body.
%   The angles are called \psi, \theta, \phi and their relation with the
%   other angles is
%   \phi=\gamma-\pi/2 mod 2\pi, \theta=\beta,
%   \psi=\alpha+\pi/2 mod 2\pi.
%
% LITERATURE: D. Hearn and M.P. Baker, Computer Graphics, Prentice Hall
%             J. Foley and A. van Dam, Computer Graphics, Addison-Wesly
%             F. Scheck, Mechanics, Springer 
%
% SYNOPSIS:
%  image_out = rotation3d(image_in, method, RotationParameters, CenterOfRotation, interploation, bgval)
%
%  method
%      Computation method: 'direct', '9 shears'
%  RotationParameters:
%      DIRECT: 1x3 array containing the Euler angles [alpha beta gamma]
%              3x3 array containing the rotation matrix
%              4x4 array containing the transformation matrix
%                  in homogenous coordinates
%      9 SHEARS: 1x3 array containing the Euler angles [alpha beta gamma]
%  CenterOfRotation:
%      DIRECT: 1x3 array containing the rotation center
%              if CenterOfRotation==[] then CenterOfRotation=center of the image
%      9 SHEARS: no choice (center of the image)
%  interpolation
%      DIRECT: no choice (linear)
%      9 SHEARS: 'default', 'bspline', '4-cubic', '3-cubic', 'linear', 'zoh'.
%  bgval
%      DIRECT: no meaning
%      9 SHEARS: Value used to fill up the background. String containing
%      one of the following values: 'zero', 'min', 'max'.
%
%  All angles in rad.
%
% SEE ALSO: rotation, rotate

% (C) Copyright 1999-2001               Pattern Recognition Group
%     All rights reserved               Faculty of Applied Physics
%                                       Delft University of Technology
%                                       Lorentzweg 1
%                                       2628 CJ Delft
%                                       The Netherlands
%
% Bernd Rieger, Feb 2001
% 18-19 September 2007: Modified the "RotationParameters" input parameter, GETPARAMS
%                       now tests the sizes. "CenterOfRotation" now needs to be empty
%                       to use the default value.

function out = rotation3d(varargin)

d = struct('menu','Manipulation',...
           'display','3D general rotation ',...
           'inparams',struct('name',       {'in',         'method','rotpara',                     'center',         'inter',               'bgval'},...
                             'description',{'Input image','Method','Euler Angles/Rotation Matrix','Rotation Center','Interpolation Method','Background Value'},...
                             'type',       {'image',      'option','array',                       'array',          'option',              'option'},...
                             'dim_check',  {0,            0,       {[1,3],[3,3],[4,4]},           {[],1},           0,                     0},...
                             'range_check',{[],{'direct','9 shears'},'R','R+',{'bspline','4-cubic','3-cubic','linear','zoh'},{'zero','min','max'}},...
                             'required',   {1,            1,       1,                             0,                0,                     0},...
                             'default',    {'a',          'direct',[0 0 0],                       [],               'bspline',             'zero'}...
                              ),...
           'outparams',struct('name',{'out'},...
                              'description',{'Output image'},...
                              'type',{'image'}...
                              )...
           );
if nargin == 1
   s = varargin{1};
   if ischar(s) & strcmp(s,'DIP_GetParamList')
      out = d;
      return
   end
end
if nargin >= 5 & isequal(varargin{5},'bilinear')
   varargin{5} = 'linear';
end
try
   [in,method, rotpara, center,inter, bgval ] = getparams(d,varargin{:});
catch
   if ~isempty(paramerror)
      error(paramerror)
   else
      error(firsterr)
   end
end

if ndims(in)~=3
   error('Input image not 3D.');
end
switch method
   case 'direct'
      %default for rotation center is round(size(in)/2)
      if isempty(center)
         center=round(size(in)/2);
      end
      if prod(size(rotpara))==3 %3 euler angles
         R=zeros(4,4);T=zeros(4,4);
         R(4,4)=1;T(4,4)=1;
         T(1,1)=1;T(2,2)=1;T(3,3)=1;
         %Translate
         T(1,4)=-center(1);
         T(2,4)=-center(2);
         T(3,4)=-center(3);
         %Rotate, erster index zeile, zweiter spalte
         %Rotation Axes: 3, 2', 3''
         alpha=rotpara(1);
         beta=rotpara(2);
         gamma=rotpara(3);
         R(1,1)=cos(gamma)*cos(alpha)-cos(beta)*sin(alpha)*sin(gamma);
         R(2,1)=-sin(gamma)*cos(alpha)-cos(beta)*sin(alpha)*cos(gamma);
         R(3,1)=sin(beta)*sin(alpha);
         R(1,2)=cos(gamma)*sin(alpha)+cos(beta)*cos(alpha)*sin(gamma);
         R(2,2)=-sin(gamma)*sin(alpha)+cos(beta)*cos(alpha)*cos(gamma);
         R(3,2)=-sin(beta)*cos(alpha);
         R(1,3)=sin(beta)*sin(gamma);
         R(2,3)=sin(beta)*cos(gamma);
         R(3,3)=cos(beta);
         %combine step
         B=R*T;
         %Translate back
         T(1,4)=center(1);
         T(2,4)=center(2);
         T(3,4)=center(3);
         %combine step to overall Matrix and go
         R=T*B;
      elseif isequal(size(rotpara),[3,3]) %only rotation matrix
         R = rotpara;
         R(4,4)=1;
         T=zeros(4,4);T(4,4)=1;
         T(1,1)=1;T(2,2)=1;T(3,3)=1;
         %Translate
         T(1,4)=-center(1);
         T(2,4)=-center(2);
         T(3,4)=-center(3);
         %combine step
         B=R*T;
         %Translate back
         T(1,4)=center(1);
         T(2,4)=center(2);
         T(3,4)=center(3);
         %combine step to overall Matrix and go
         R=T*B;
      elseif isequal(size(rotpara),[4,4]) %homogenous coordinates
         R = rotpara;
      else
         error('Rotation parameters invalid.');
      end
      %R
      %two low level routines for different data types too save memory
      if strcmp(datatype(in),'dfloat')
         %disp('double routine');
         out=dip_image(rot_euler_low(double(in),R));
      else
         %disp('sfloat routine');
         out=dip_image(rot_euler_low_sfloat(single(in),R));
      end
   case '9 shears'
      out=dip_rotation3d(in,rotpara(1),rotpara(2),rotpara(3), inter, bgval);
   otherwise
      error('Unkown rotation method.');
end
