function I=EulerAngles(A)
% The code calculates Euler angles needed to rotate a coordinate system into
% the coordinate system described by A. The X-axis of the A coordinate
% system is defined in the first column of A, the Y-axis in the second
% column, and the Z-axis in the third column.

% Perform normalisation. (In principal not needed)
A=A./repmat(sqrt(sum(A.*A)),3,1);

% Then check that data is delivered in a right-handed coordinate system
% (i.e. (v1×v2)⋅v3>0 must hold)
if sum(abs(A(:,3) - cross(A(:,1),A(:,2))))>.01
  % Not a right handed coordinate system
  disp('Can not create Euler Angles. Require a right handed coordinate system.')
  I=[];
  return;
end

beta_angle=acos(A(3,3));
if abs(sin(beta_angle))>(1e-8)
  % The variables with the suffix _rounding_protection
  % are to secure possible rounding error that might bring the
  % rotation out of the range [-1,1].
  % Note: You might note that the variable naming doesn't align with the
  % naming if you consider the rotation matrix to go from the average
  % direction into the current given ROI matrix. However the naming of the
  % angles makes sense if you consider the transformation from the given
  % ROI into the average Direction matrix because it is just the inverse
  % operation and thus the inverse matrix. (c.f. Note on Euler angles and
  % how to calculculate them.docx). So instead of using A(1,3) in the next
  % line we actually use A(3,1).
  gamma_angle_rounding_protection = max(-1, min(1, A(3,1)/sin(beta_angle) )); % this is sin(gamma)
  gamma_angle_sin = asin(gamma_angle_rounding_protection);
  gamma_angle_rounding_protection=max(-1, min(1, A(3,2)/sin(beta_angle) ));
  gamma_angle_cos=acos(gamma_angle_rounding_protection);
  [gamma_angle_cos,gamma_angle_sin]=CheckQuadrants(gamma_angle_cos,gamma_angle_sin);

  % Here we do the same as above. First a rounding protection is applied
  % and then the inverse operations are applied to obtain alpha. Note again
  % that it is the INVERSE here. So instead of A(3,1) we select A(1,3)
  % etc..
  temp=max(-1,min(1,A(1,3)/sin(beta_angle)));
  alpha_angle_sin=asin(temp);
  temp=max(-1,min(1,-A(2,3)/sin(beta_angle)));
  alpha_angle_cos=acos(temp);
  [alpha_angle_cos,alpha_angle_sin]=CheckQuadrants(alpha_angle_cos,alpha_angle_sin);
else
  % This is the special case where the equation for the B matrix in Note on Euler angles and
  % how to calculculate them.docx because much simpler :-)
  gamma_angle_cos=0;
  gamma_angle_sin=0;
  alpha_angle_sin=asin(A(2,1));
  alpha_angle_cos=acos(A(1,1));
  [alpha_angle_cos,alpha_angle_sin]=CheckQuadrants(alpha_angle_cos,alpha_angle_sin);
end
I=[alpha_angle_cos,beta_angle,gamma_angle_cos];

%Validation in principel not needed

if abs(gamma_angle_sin-gamma_angle_cos)>0.1*pi/180 || abs(alpha_angle_sin-alpha_angle_cos)>0.1*pi/180
  I=[]; %Error has occured
  warning('An error occured in the Euler angle calculation. Possibly due to machine precision')
else
  %Calc full matrix rotation to transform one set of coordinates to the
  %values in the new cooridnate system
  RotAroundZ=[cos(I(1)),sin(I(1)),0;-sin(I(1)),cos(I(1)),0;0,0,1];
  RotAroundNewX=[1,0,0;0,cos(I(2)),sin(I(2));0,-sin(I(2)),cos(I(2))];
  RotAroundNewZ=[cos(I(3)),sin(I(3)),0;-sin(I(3)),cos(I(3)),0;0,0,1];
  TotTransform=RotAroundNewZ*RotAroundNewX*RotAroundZ;
  TestMat=TotTransform-A';
  if sum(abs(TestMat(:)))>.01
    I=[]; %Error has occured
    warning('An error occured in the Euler angle calculation. Possibly due to machine precision')
  end
end
end
function [cosval_angle_new,sinval_angle_new]=CheckQuadrants(cosval_angle,sinval_angle)
cosval_angle_new=cosval_angle;
sinval_angle_new=sinval_angle;
minval=1e-5*pi/180;
cosarray=[cosval_angle,-cosval_angle,cosval_angle,-cosval_angle];
sinarray=[sinval_angle,sinval_angle,pi-sinval_angle,pi-sinval_angle];
index=abs((mod(cosarray,2*pi)-mod(sinarray,2*pi)))<minval;
if index(1)
  %all ok
elseif index(2)
  cosval_angle_new=-cosval_angle;
elseif index(3)
  sinval_angle_new=pi-sinval_angle;
elseif index(4)
  cosval_angle_new=-cosval_angle;
  sinval_angle_new=pi-sinval_angle;
end
cosval_angle_new=mod(cosval_angle_new,2*pi);
sinval_angle_new=mod(sinval_angle_new,2*pi);
if cosval_angle_new>pi
  cosval_angle_new=cosval_angle_new-2*pi;
end
if sinval_angle_new>pi
  sinval_angle_new=sinval_angle_new-2*pi;
end
end
