% Figure to illustrate the meaning of the distances R and S.
% NB I have found the pupils don't come out OK on some Matlab installations
% - no idea why; Matlab bug I assume.
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

% Interocular distance:
interoculardistance = 8;

% Coordinates of fixation point
fix = [3 9 12]';

% Head-centered coordinates of viewed object
objectP = [2 7 6]';
objectcol = [0.8 0 0];

% Deduce eye coords necessary to fixate desired point
V = -atand(fix(2)/fix(3));
nrm = norm(fix(2:3));
HL = atand((fix(1)-interoculardistance/2)/nrm);
HR = atand((fix(1)+interoculardistance/2)/nrm);

% Coordinates of eye:
eyes.radius = 1.3;
eyes.nsurface = 600;
eyes.eyeballcolor = [0.75 0.75 0.5];
eyes.eyeballtransparency = 0.8;
% Specify eye positions
eyes.eyepos.nodal = [0.5; 0 ; 0]*interoculardistance;
eyes.eyepos.azimuth = HL;
eyes.eyepos.elevation =V;
eyes.eyepos.torsion = 0;
% 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;
eyes.eyepos(2) = eyes.eyepos;
eyes.eyepos(2).azimuth = HR;
eyes.eyepos(2).nodal = [-0.5; 0 ; 0]*interoculardistance;

% Cyclopean eye:
cyclopean.radius = eyes.radius;
cyclopean.eyeballcolor = [0.7 0.7 1];
cyclopean.eyeballtransparency = 0.8;
cyclopean.eyepos.nodal = [0; 0 ; 0];
cyclopean.eyepos.azimuth = (HL+HR)/2;
cyclopean.eyepos.elevation =V;
cyclopean.eyepos.torsion = 0;


% 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.

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 = [-8 8;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 8;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)


% Draw the eyes
DrawEyes(eyes);
DrawEyes(cyclopean);

% Plot optic axes:
plotvector(fix,'o','markerfacecol',opticaxiscol,'col',opticaxiscol)
plotvector([eyes.eyepos(1).nodal fix],'col',opticaxiscol,'linew',opticaxislw)
plotvector([eyes.eyepos(2).nodal fix],'col',opticaxiscol,'linew',opticaxislw)
% Plot line from origin to fixation point
plotvector([[0 0 0]' fix],'col','b','linew',opticaxislw,'markersize',12)
text(fix(1),fix(2),fix(3)+0.5,'F','fonts',fs)

text(0.2,-0.2,0,'O','fonts',fs,'col','k')

% plot object
plotvector(objectP,'o','markerfacecol',objectcol,'col',objectcol)
plotvector([[0 0 0]' objectP],'-','col',objectcol)
text(objectP(1),objectP(2),objectP(3)+0.5,'P','fonts',fs,'col',objectcol)
% Project object P onto cyclopean gaze:
% lambda fix = norm(objectP) cos(theta) * fix/norm(fix)
% norm(fix) norm(objectP) cos(theta) = dot(fix,objectP)
lambda = dot(fix,objectP) / norm(fix)^2;
projction = lambda * fix;
% plot line from origin to projection of P onto gaze axis, shifted slightly
% so as to overwrite gaze axis:
plotvector([[1e-3 0 0]' projction+[1e-3 0 0]'],'-','col',objectcol,'linew',opticaxislw)
plotvector([objectP projction],'-','col',objectcol)
% Mark right angle:
tmp = 0.03;
A = projction*(1-tmp);
B = projction*(1-tmp) + norm(projction*tmp)*(objectP-projction)/norm(objectP-projction);
C = projction + norm(projction*tmp)*(objectP-projction)/norm(objectP-projction);
plotvector([A B C],'-','col',objectcol)

% Mark R and S
text(objectP(1)/2,objectP(2)/2+0.5,objectP(3)/2,'R','fonts',fs,'col',objectcol,'fontname','times','fontangle','italic')
text(projction(1)/2,projction(2)/2-0.1,projction(3)/2+0.5,'S','fonts',fs,'col',objectcol,'fontname','times','fontangle','italic')


camup([0 1 0])
camtarget([0.6005    3.5846    1.9623])
campos([20 6 17])
camva(20)
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')
% Nb wouldn't print correctly; had to use screen capture and/or cut and
% paste the good bits from the following two files:
print -dtiff -r1000 -opengl Fig_DistancesRSb
