function I=anaview(varargin)
% ANAVIEW   Generate anaglyph from 3-D plot
%  I=ANAVIEW generates the anaglyph version of the 3-D plot in the
%  current axis. 
%
%  The left-right view disparity is computed assuming that the
%  distance from the camera position to the camera target is of one
%  meter, that the two eyes have distance of 0.1 meters and that the
%  screen is at a distance of 1 meter from the observer. 
%
%  Target and image distance can be changed (see options
%  below). Generally speaking:
%
%  + Increasing/decreasing the target distance will move the object
%    away/closer to the observer.
%
%  + The target distance is represented relative to the screen
%    distance. So increasing/decreasing the screen distance while
%    keeping the object distance and the true distance to the screen
%    constant, will make the object to appear closer/farther.
%
%  + Object points that are farther/closer than the screen will
%    appear behind/above the screen.
%
%  All these operations do not change the amount of perspective
%  distortion as the object is scaled up/down and the focal length
%  stays the same. The focal length can also be changed (see options
%  below). A smaller focal length will cause the object to appear
%  deeper.
%
%  The following option-value pairs are supported:
%  'TargetDistance' [1]
%    Distance of camera target. The closer, the more pronounced the
%    disparity.
%
%  'ImageDistance' [1]
%    Distance of monitor surface. The closer, the deeper the
%    objects are pushed.
%
%  'ScaleFocal' [1]
%    Multiply camera focal length by this factor, while moving
%    closer/farther to the object to keep the target of the same
%    size. This step is applied before setting the target and
%    screen distances.
%
%  'Background' [0.1]
%    Intensity of the background tone (green channel), in the 0-1
%    range.

% AUTORIGHTS
% Copyright (C) 2006 Andrea Vedaldi
%       
% This file is part of Anaview.
% 
% Anaview is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2, or (at your option)
% any later version.
% 
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
% 
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software Foundation,
% Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.

h  = gca ;

% PARAMETERS
%   Lengths are measured in meters [m], angles in radians.
%   ned = nominal eye distance (distance between eyes)
%   nsd = nominal screen distance (distance to the screen)
%   ntd = nominal target distance (distance to the 3-D target)
%   nsh = nominal screen height
%   a   = aperture (vertical, full angle)
%   sf  = rescale focal factor
%   bt  = background tone

ned = 0.05 ;
nsd = 1.0 ;
ntd = 1.0 ;
sf  = 1.0 ;
bt  = 0.1 ;

for k=1:2:length(varargin)
	switch lower(varargin{k})
    case 'targetdistance'
      ntd = varargin{k+1} ;
      
    case 'imagedistance'
      nsd = varargin{k+1} ;
    
    case 'scalefocal'
      sf = varargin{k+1} ;
      
    case 'background'
      bt = varargin{k+1} ;
      
    otherwise
      error (['Unknown option ''', varargin{k} '''.']) ;
	end
end

% --------------------------------------------------------------------
%                                                            Copy axes
% --------------------------------------------------------------------

% get current figure
sfig = gcf;

% create a new figure with same backgorund color and size,
fig = figure ;
set(fig,'Color',get(sfig,'Color'),...
				'Position',get(sfig,'Position') ) ;

% copy the current axis onto the new figure
h   = copyobj(h,fig) ;
axes(h) ;

% make sure we use perspective projection
set(h,'Projection','perspective') ; 

% --------------------------------------------------------------------
%                                                      Adjust aperture
% --------------------------------------------------------------------

% fix ratios and view angle
da  = get(h,'DataAspectRatio') ;
pba = get(h,'PlotBoxAspectRatio') ;
a   = get(h,'CameraViewAngle') * pi/180 ;

% maps axes units to 3-D world units
A = diag(da) ;

% camera position and target
p = A * campos' ;
t = A * camtarget' ;

% adjust
f = 1 / tan(a / 2) ;
f = f * sf ;
p = t + (p - t) * sf ;
a = 2 * atan(1 / f) ;

% set
campos(    inv(A) * p ) ;
camtarget( inv(A) * t ) ;
set(h,'DataAspectRatio',   da) ;
set(h,'PlotBoxAspectRatio',pba) ;
set(h,'CameraViewAngle',   a * 180/pi) ;

% expand axis to fill the whole figure
set(h,'Position',[0 0 1 1]) ;

% --------------------------------------------------------------------
%                                                       Compute optics
% --------------------------------------------------------------------

% compute equivalent eye distance
p   = A * campos' ;
t   = A * camtarget' ;
d   = norm(t - p) ;
eed = ned * d/ntd ;

% this is also size (in pixels) of the axes
pos = get(fig, 'Position') ;

% compute projected eye distance in pixels
nsh = 2 * tan(a/2) * nsd ;
psh = pos(4) ; % screen height in pixel
ped = round( ned * psh/nsh ) ;

% make room
set(fig, 'Position', pos + [0 0 ped 0]) ;

% --------------------------------------------------------------------
%                                                   Generate two views
% --------------------------------------------------------------------

% get current camera tf matrix and distance from target
[R,T] = camget(h) ;

% compute displaced cameras
Tl = T - [ eed/2 ; 0 ; 0] ;
Tr = T + [ eed/2 ; 0 ; 0] ;

% left image
camset(h,R,Tl) ;
drawnow ;
Il = getframe(fig) ;
Il = rgb2gray(Il.cdata) ;
Il = Il(:,ped+1:end,:) ;

% right image
camset(h,R,Tr) ; 
drawnow ;
Ir = getframe(fig) ;
Ir = rgb2gray(Ir.cdata) ;
Ir = Ir(:,1:end-ped,:) ;

% we do not this anymore
close(fig) ;

% --------------------------------------------------------------------
%                                                            Mix views
% --------------------------------------------------------------------

I = cat(3,Il,uint8(255*bt*ones(size(Il))),Ir) ;

% --------------------------------------------------------------------
function [R,T] = camget(h)
% --------------------------------------------------------------------

A = diag(get(h,'DataAspectRatio')') ;

up     = A * camup(h)' ;
target = A * camtarget(h)' ;
pos    = A * campos(h)' ;

z = target - pos ;
x = cross(up,z) ;
y = cross(z,x) ;

x = x(:) / norm(x) ;
y = y(:) / norm(y) ;
z = z(:) / norm(z) ;

T = pos(:) ;
R = [x y z] ;

R = R' ;
T = - R * T ;

% --------------------------------------------------------------------
function camset(h,R,T)
% --------------------------------------------------------------------

A = diag(get(h,'DataAspectRatio')') ;

R = R' ;
T = - R * T ;

pos = T ;

x = R(:,1) ;
y = R(:,2) ;
z = R(:,3) ;

target = pos + z ;
up     = y ;

campos(   h, inv(A) * pos    ) ;
camtarget(h, inv(A) * target ) ;
camup(    h, inv(A) * up     ) ;
