% Figure to illustrate the Helmholtz coordinates for describing eye
% position
clear all
close all

fs=14; % fontsize to use for axis labels
fsa=18; % fontsize to use for panel labels
fst = 18; % fontsize to use for angle labels (V,H,T)
opticaxiscol = [1 1 0]*0.5; % colour to use for optic axis
opticaxislw = 3; % linewidth to use for optic axis
elcol = 'r'; % colour to use for elevation
azcol = 'b';
torcol = [0 0.5 0];
anglw = 2; % linewidth for marking gaze angles

% Head-centered coordinates of viewed object
objectP = [-6 7 10];

% Coordinates of eye:
eyes.radius = 1.75;
eyes.nsurface = 500;
eyes.eyeballcolor = [0.75 0.75 0.5];
eyes.eyeballtransparency = 0.98;
% Specify eye positions
eyes.eyepos.nodal = [0; 0 ; 0];
eyes.eyepos.azimuth = 45;
eyes.eyepos.elevation = 30;
eyes.eyepos.torsion = 30;
% Colors in which to draw coordinate systems
eyes.azcoordcol = 'k' ;
eyes.elcoordcol = 'k' ;
% Specify iris and pupil
eyes.irisradius = 0.8*eyes.radius;
eyes.iriscolor = [0.75 0.75 0.5]/2;
eyes.pupilcolor = [0 0 0 ];
eyes.pupilradius = 0.95*eyes.radius;
Rot = RotationMatrix(eyes.eyepos) ;
% Just the elevation component:
cosV = cosd(eyes.eyepos.elevation);
sinV = sind(eyes.eyepos.elevation);
RotV = [ 1 0 0 ; 0 cosV -sinV ; 0  sinV cosV ];
% Just the azimuth component:
cosH = cosd(eyes.eyepos.azimuth);
sinH = sind(eyes.eyepos.azimuth);
RotH = [cosH 0 sinH ; 0 1 0 ; -sinH 0 cosH];

% Planes to be drawn:
XZplanecol = [1 0.8 0.8]; % Colour for showing XZ plane in
ZYplanecol = [0.8 0.8 1]; % Colour for showing XZ plane in
gazeplanecol = [1 1 0.5]; % Colour for gaze plane
planetransp=1; % transparency of both planes.
XZplanetransp=0.8; % transparency of both planes.

% Specify the coordinate system to be painted on the retina
coordsystem.name = 'azlongellat';
coordsystem.nlines = 25;
coordsystem.spacing = 15;
coordsystem.sca  = 1.001; % amount by which to exaggerate eyeball radius so coord lines are visible
coordsystem.linewidth = 0.5;
coordsystem.linewidththick = 2; % used for plotting meridians

figure('pos',[ 479   161   686   529])

arrowheadlength = 0.2;  % increase to make arrows longer
arrowangle = 2; % increase to make arrows more pointed
% Draw interocular axis
Xaxis = [ [-1 1]*6;[0 0];[0 0] ] + [1 1 ; 0 0 ; 0 0];
plotvector(Xaxis,'k-','linew',2)
hold on
% Add on arrowhead
plot3(Xaxis(1,2) + [-1 0 -1]*arrowheadlength*arrowangle , Xaxis(2,2) * [1 1 1], [-1 0 1]*arrowheadlength,'col','k','linew',2)
text(Xaxis(1,2),0.,-0.2,'X','fonts',fs)
% Mark Y axis
Yaxis = [0 0;-1 3.5;0 0];
plotvector(Yaxis,'k-','linew',2)
% Add on arrowhead
plot3( Yaxis(1,2) * [1 1 1] , Yaxis(2,2) + [-1 0 -1]*arrowheadlength*arrowangle, [-1 0 1]*arrowheadlength,'col','k','linew',2)
text(0,Yaxis(2,2)-0.2,-0.2,'Y','fonts',fs)
% Mark Z axis
Zaxis = [0 0;0 0;-1 8];
plotvector(Zaxis,'k-','linew',2)
% Add on arrowhead
plot3( [-1 0 1]*arrowheadlength , [0 0 0], Zaxis(3,2) + [-1 0 -1]*arrowheadlength*arrowangle,'col','k','linew',2)
text(0,0.5,Zaxis(3,2),'Z','fonts',fs)


% ZY plane:
ZYplane = [ [0 0 0 0 0 ] ; [-1 1 1 -1 -1]*3 ; [0 0 1 1 0]*6];
patch(ZYplane(1,:),ZYplane(2,:),ZYplane(3,:),ZYplanecol,'edgecol','none','facealpha',planetransp,'edgealpha',planetransp);
% This is a rectangle in the XZ plane:
XZlimit = eyes.radius*2.5;
XZplane = [[-1 1 1 -1 -1]*XZlimit ; [0 0 0 0 0] ;  [0 0 1 1 0]*6];
patch(XZplane(1,:),XZplane(2,:),XZplane(3,:),XZplanecol,'edgecol','none','facealpha',XZplanetransp,'edgealpha',XZplanetransp);
% Rotate it through the elevation angle, so as to get the gaze plane:
gazeplane = RotV*XZplane;
patch(gazeplane(1,:),gazeplane(2,:),gazeplane(3,:),gazeplanecol,'edgecol','none','facealpha',planetransp,'edgealpha',planetransp);

% Draw the eyes
DrawEyes(eyes);
% Mark on coordinate system:
eyecoordhandles = DrawEyeCoords(eyes,coordsystem);

% Plot optic axis:
fovea = Rot * [0 ; 0 ; -eyes.radius];
opticaxis = Rot * [0 ; 0 ; 6];
plotvector([fovea opticaxis],'col',opticaxiscol,'linew',opticaxislw)

% Show elevation angle:
th=[0:0.05:1]*eyes.eyepos.elevation;
r = 1.5;
plot3(ones(length(th))*XZlimit, -r*sind(th),r*cosd(th),'col',elcol,'linew',anglw)
plot3([1 1 1]*XZlimit,[0 0 -r*sind(eyes.eyepos.elevation)],[r 0 r*cosd(eyes.eyepos.elevation)],'col',elcol,'linew',anglw)
text(XZlimit,-0.3,1.3,'V','fonts',fst,'col',elcol,'fontw','bold');

% Show torsion angle
% Show intersection between eyeball and gaze plane, between 0 and gaze azimuth:
th=[eyes.eyepos.azimuth:0.05:90];
nrm = sqrt(tand(eyes.eyepos.elevation).^2 + tand(th).^2 + 1 );
X = eyes.radius * tand(th)./ nrm;
Y = -eyes.radius * tand(eyes.eyepos.elevation)./ nrm;
Z = eyes.radius./ nrm;
plot3(X,Y,Z,'col',torcol,'linew',anglw)
% Eyeball crosses X axis at X=eyes.radius,Y=0,Z=0. Find what lat & long on
% the rotated eyeball this corresponds to
Reyeball = inv(Rot) * [1; 0 ; 0];
azlat = asind(Reyeball(1));
ellat = asind(Reyeball(2));
ellong = asind(Reyeball(2)./cosd(azlat));
azlong = asind(Reyeball(1)./cosd(ellat))
% Plot line of elevation longitude from 0 to ellong, at constant azlong:
th=[0:-0.01:ellong];
nrm = sqrt(1+tand(th).^2+tand(azlong).^2);
curve = eyes.radius.*[tand(azlong)./nrm;tand(th)./nrm;1./nrm];
plotvector(Rot*curve,'-','col',torcol,'linew',anglw)
% Plot retinal horizontal meridian 
th=[0:0.05:azlong];
% before eyeball is rotated:
X = sind(th);Z = cosd(th);
curve = 1.001*coordsystem.sca*eyes.radius.*[X;zeros(size(th));Z];
plotvector(Rot*curve,'-','col',torcol,'linew',anglw)
text(2.,0.7,0.2,'T','fonts',fst,'col',torcol,'fontw','bold');


% Show gaze azimuth:
th=[-eyes.eyepos.elevation:0.05:90];
curve = [ zeros(size(th)) ; eyes.radius*sind(th) ; eyes.radius*cosd(th) ];
plotvector(curve,'-','col',azcol)
% Rotate through azimuthal angle:
curve = RotH * curve;
plotvector(coordsystem.sca*curve,'-','col',azcol,'linew',anglw)
text(1.2,0.2,2.2,'H','fonts',fst,'col',azcol,'fontw','bold');
% Show intersection between eyeball and XZ plane, between 0 and gaze
% azimuth:
th=[0:0.01:eyes.eyepos.azimuth];
nrm = sqrt( tand(th).^2 + 1 )/coordsystem.sca;
X = eyes.radius * tand(th)./ nrm;
Y = zeros(size(th));
Z = eyes.radius./ nrm;
plot3(X,Y,Z,'col',azcol,'linew',anglw)



camup([0 1 0])
camtarget([ 1.6704    0.0105    2.5561])
campos([10 5 5])
camva(40)
lighting phong
%light('pos',Zaxis(:,2)/2+Yaxis(:,2)*2)
camlight left
camlight headlight
%    set(gca,'Projection','persp') % NB apparently some Matlab bug --
%    persepctive rendering gets the intersections wrong, such that the
%    horizontal plane no longer intersects the horizontal meridians. Leaves
%    as orthogonal.
axis equal off


set(gcf,'color','w')
print -dtiff -r700 -zbuffer Fig_HelmholtzEyeCoords