function VMSReport(varargin)
    % VMSReport(varargin)
    % Produce a printed report on results from running SilentGFtbox
    % before it finishes it invites you to rotate the figures
    % (rotate one using the sliders and the others follow suit)
    % it then invites you to print everything
    %
    % Arguments come in command,value pairs
    % 'Path', path to project default, pwd (default);
    % 'Project', name of project in the path
    % 'Experiment',  either 'All' (the default) where all the experiments in
    %            path/project/movies/result_directory/*
    %            are compiled into a printed report
    %                or {'experiment_name','experiment_name',...}
    % 'MORPHOGEN','name' default 'Polariser'
    % 'DRAWEDGES',integer, default 0 none, 1 outside edges only, 2 all edges
    % 'DRAWGRADIENTS',logical default true
    % 'THICK',logical default true
    % 'DRAWSECONDLAYER',logical default false (cell layer)
    % 'PRINT',logical default false
    % 'VISIBLE',logical default true if true do not display figures
    % 'ROTATION',[azimuth,elevation] default is as stored in project
    % 'CELLS',logical default true
    % 'STEREO',logical default false
    % usage
    % VMSreport('Project','HoleMotif4')
    %
    %J. Andrew Bangham 2008

%     VMSreport('Path','D:\ab\Matlab stuff\DArT_ModelstrunkAribidopsisLeaf\SciencePaper','project',...
%         'GPT_Lateinhpgrad_clust_20111204','Experiment','All','morphogen','KPAR',...
%         'SNAPFIG',false,'Rotation',[180,90],'PM_MORPHOGEN','ID_PERIM');
%     
    %VMSreport('Path','E:\DArT_Models\Aribidopsis\Leaf\SciencePaper','project',...
    %     'GPT_ArLeaf_ScienceResumission_20111205','Experiment','All','morphogen','KPAR',...
    %     'SNAPFIG',false,'Rotation',[180,90]);%,'PM_MORPHOGEN','ID_PERIM');
    % 'Layout', options are one of: 'TwoUp', 'FourUp', 'OneUp'
    %     close all
    pathname=pwd;
    projectname='HoleMotif4';
    layout='TWO_UP';
    experiment_list='ALL';
    printflag=false;
    morphogen='POLARISER';
    drawedges=1; %2;
    drawgradients=false;
    thick=true;
    drawsecondlayer=true;
    pngFlag=true;
    verbose=true;
    invisible=false;
    azimuth=135; %-45;
    elevation=30;
    viewcellsflag=true;
    stereo=false;
    flatten=false;
    NickAndFlatten=0;
    pm_morphogen='';
    snapdragonfigure=false;
    for i=1:2:length(varargin)
        command=upper(varargin{i});
        opt=varargin{i+1};
        switch command
            case 'PM_MORPHOGEN'
                pm_morphogen=opt;
            case 'SNAPFIG'
                snapdragonfigure=opt;
            case 'PATH'
                pathname=opt;
            case 'PROJECT'
                projectname=opt;
            case 'DRAWSECONDLAYER'
                drawsecondlayer=opt;
            case 'EXPERIMENT'
                experiment_list=opt;
            case 'VISIBLE'
                invisible=~opt;
            case 'PRINT'
                printflag= opt==true;
            case 'LAYOUT'
                layout=opt;
            case 'MORPHOGEN'
                morphogen=opt;
            case 'DRAWEDGES'
                drawedges=opt;
            case 'DRAWGRADIENTS'
                drawgradients=opt;
            case 'THICK'
                thick=opt;
            case 'ROTATION'
                azimuth=opt(1);%135; %-45;
                elevation=opt(2);%30;
            case 'CELLS'
                drawsecondlayer=opt;
            case 'STEREO'
                stereo=opt;
            case 'FLATTENTIME'
                if ischar(opt)
                    flatten=str2num(opt);
                else
                    flatten=opt;
                end
            case 'NICKANDFLATTEN'                
                if ischar(opt)
                    NickAndFlatten=str2num(opt);
                else
                    NickAndFlatten=opt;
                end
            otherwise
                error(sprintf('command %s is not recognised',command))
        end
    end
    if ~iscell(experiment_list)
        experiment_list={experiment_list};
    end
    if strcmpi(experiment_list{1},'ALL')
        experiment_list={};
        DD=dir(fullfile(pathname,projectname,'runs'));
        if isempty(DD)
            error(sprintf(...
                'Unable to locate the project directory: %s',...
                fullfile(pathname,projectname,'runs')))
        end
        for i=1:length(DD)
            if DD(i).isdir && ~strcmpi(DD(i).name,'.') && ~strcmpi(DD(i).name,'..')
                experiment_list{end+1}=DD(i).name;
            end
        end
    end
    allfigs=[];
    all_axes=[];
    for project_number=1:length(experiment_list)
        experiment=experiment_list{project_number};
        if project_number==1
            %             if ~invisible
            %                 figs=FlyleafFigs(pathname,projectname,experiment,pngFlag,printflag,invisible);
            %             end
            if printflag
                for i=1:length(figs)
                    figure(figs(i));
                    print(figs(i));
                end
            end
            if ~invisible
                %                 allfigs=[allfigs,figs];
                allfigs=[allfigs];
            end
        end
        % Set the current fig to the last one used for result plotting
        if isempty(allfigs)
            figsR = [];
        else
            figsR = allfigs(end);
        end
        [figsR,data]=ResultsFigs(pathname,projectname,experiment_list{project_number},...
            morphogen,drawedges,drawgradients,thick,pngFlag,printflag,invisible,...
            drawsecondlayer,azimuth,elevation,viewcellsflag,stereo,flatten,...
            snapdragonfigure, figsR,pm_morphogen);
        % Add the current fig to allfigs if it's not there already
        if ~isempty(figsR)
            if isempty(find(figsR(end))==allfigs)
                allfigs=[allfigs,figsR];
            end
        end
        %         all_axes=[all_axes, data.current_axes];
    end
    %     hlink=link_all_the_figures(all_axes);
    if printflag
        PrintFigs(allfigs);
        %     else
        %         ans=input('Axis locked together for rotations. Press Return to end or P to print ','s')
        %         if strcmp(upper(ans),'P') %printflag
        %             PrintFigs(allfigs);
        %         end
    end
    % Close the figures that have been used
    for i = 1:length(allfigs)
        close(allfigs(i)); 
    end
end

function PrintFigs(allfigs)
    disp('PrintResults');
    for i=1:length(allfigs)
        figure(allfigs(i));
        print(allfigs(i));
    end
end

function [figs,data]=ResultsFigs(pathname,projectname,experiment,... %3
        morphogen,drawedges,drawgradients,thick,pngFlag,printflag,...%6
        invisible,drawsecondlayer,azimuth,elevation,viewcellsflag,stereo,flatten,...%7
        snapdragonfigure, figs,pm_morphogen)
    if nargin<4
        morphogen='POLARISER';
    end
    if nargin<5
        drawedges=false;
    end
    if nargin<6
        drawgradients=true;
    end
    if nargin<7
        thick=true;
        pngFlag=false;
        printflag=false;
        invisible=true;
    end
    if nargin<16
        flatten=false;
    end
    if nargin<17
        snapdragonfigure=false;
    end
    %figs=[];
    data=[];
    commandname=fullfile(pathname,projectname,'runs',experiment,'CommandLine.txt');
    h=fopen(commandname,'r');
    if h==-1
        errorno=1;
    else
        errorno=0;
        commandstring=fgetl(h);
        if any(strfind(lower(commandstring),'flatten'))
            could_flatten=true;
        else
            could_flatten=false;
        end
        strind=strfind(lower(commandstring),'mutantnumber')
        mutantnumber_str='';
        if ~isempty(strind)
            s1=(commandstring(strind+length('mutantnumber')+2:end));
            ind2=strfind(s1,',');
            s2=s1(1:ind2-1);
            mutantnumber=str2num(s2);
            switch mutantnumber
                case 1 
                    mutantnumber_str='default';
                case 2
                    mutantnumber_str='no polariser from prox';
                case 3
                    mutantnumber_str='constitutive cenorg';
                case 4
                    mutantnumber_str='dich not influence distorg';
                case 5
                    mutantnumber_str='35S RAD';
                case 6
                    mutantnumber_str='backpetals';
                case 7
                    mutantnumber_str='backpetals 35S RAD';
                otherwise
                    mutantnumber_str='mutantnumber not recognised';
            end                        
            commandstring=[commandstring,' ',mutantnumber_str];
            % mutantnumber 1 no production of polariser from prox (still has
            % background)
            % mutantnumber 2 constitutive cenorg
            % mutantnumber 3, dich no longer influence distorg
            % mutantnumber 4, 35S RAD (constitutive rad)
            % mutantnumber 5, backpetals constitutive cyc
            % mutantnumber 6, constitutive cyc and rad
%         else
%             could_flatten=false;
        end
        
        fclose(h);
        
        strind=strfind(lower(commandstring),'mut_list')
        mutantname_str='Wild';
        if ~isempty(strind)
            s1=(commandstring(strind+length('mut_list')+3:end));
            ind2=strfind(s1,'}');
            s2=s1(1:ind2-1);
            ind3=strfind(s2,'''''');
            k=1;
            mutantname_str='';
           while ~isempty(ind3)
                mutantname_str=[mutantname_str,'_',s2(5:ind3(1)-1)];
                s2=s2(ind3(1)+1:end);
                ind3=strfind(s2,'''''');
           end
           mutantname_str=[mutantname_str,'_',s2(5:end-1)];
        end
        
        % loop around all .mat files in the experiment directory
        filenames=fullfile(pathname,projectname,'runs',experiment,'meshes','*.mat');
        DD=dir(filenames);
        %August 2011
        %figs=[];
        ss=get(0,'screensize');
        %         if invisible
        %             set(gcf,'visible','off');
        %         end
        position=[7   ss(4)-998-100   864   998];
        paperposition=[0.634518 1.17743 20.0 27.5];
        disp(sprintf('There are %d computed mesh files',length(DD)-1));
        for i=2:length(DD)
            disp(sprintf('%s',DD(i).name));
        end
        prefix='';
        for i=2:length(DD)
            printflag=false;
            file=DD(i).name;
            if true %~isempty(strfind(file,'057'))
                %             if isempty(azimuth)
                %name=sprintf('%s%d.png',projectname,i);
                name=sprintf('%s_%-2d%s_%s.png',projectname,i,lower(mutantname_str),lower(mutantnumber_str));
                prefix=experiment(end-3:end-1);
                printfilename=fullfile(pathname,projectname,'runs',experiment,'meshes',[prefix,name]);
                %             else
                %                 name=sprintf('%s%d_%d%d.png',projectname,i,round(azimuth),round(elevation));
                %                 printfilename=fullfile(pathname,projectname,'movies',experiment,'meshes',name);
                %             end

                if (exist(printfilename)==2)
                    disp(sprintf('Already created png %s',printfilename))
                else

                    m = loadmesh_anyfile( [], fullfile(pathname,projectname,'runs',experiment,'meshes',file));
                    if isfield(m,'nodes')
                        % If need to output data for shape model then this is where
                        if ~isempty(pm_morphogen)
                            %pm_morphogen
                            [id_perim_i,id_perim_p,id_perim_a,id_perim_l] = getMgenLevels( m, pm_morphogen );
                            pointmodel_filename=fullfile(pathname,projectname,'runs',experiment,'meshes',[file(1:end-4),'_pm.mat']);
                            x=m.nodes(find(id_perim_p),1);
                            y=m.nodes(find(id_perim_p),2);
                            xy=[x';y'];
                            save(pointmodel_filename,'-mat','xy');
                        end
                        %             if invisible
                        %                 m.plotdefaults.invisibleplot=true;
                        %             end
                        m.pictures=[];
                        flattened=false;
                        if could_flatten
                            ind=strfind(file,'_s000');
                            suffix=file(ind+5:end-4);
                            stagetime=stageTextToTime(suffix);
                            if flatten~=false
                                if stagetime>=flatten
                                    flattened=true;
                                end
                            end
                        end
                        if snapdragonfigure
                            morphogen = {
                                'ID_LIP'
                                'ID_PALATE'
                                'ID_CYC'
                                'ID_DICH'
                                'ID_PROXORG'
                                'ID_CENORG'
                                'ID_DISTORG'
                                'ID_DIV'};
                            morphogen = {
                                'ID_LIP'
                                'ID_DICH'
                                'ID_UPTUBE'
                                'ID_TUBE'};
                            m.mgenposcolors(:,11)=[
                                0.6
                                0.2
                                0];
                            m.mgenposcolors(:,10)=[
                                1
                                1
                                0];
                            m.mgenposcolors(:,12)=[
                                1
                                0
                                0];
                            m.mgenposcolors(:,21)=[
                                0.6824
                                0.4667
                                0];
                            m.mgenposcolors(:,35)=[
                                0
                                1
                                1];
                            m.mgenposcolors(:,41)=[
                                0
                                1
                                1];
                            m.mgenposcolors(:,40)=[
                                0.4941
                                0.4941
                                0.4941];
                            m.mgenposcolors(:,42)=[
                                0.4
                                0.4
                                0.4];
                            m = leaf_setbgcolor( m, [1, 1, 1] );
                        end
                        % New August 2011 - reuse figures.
                        if ~isempty(figs)
                            m.pictures = figs(end); 
                        end
                        if flattened
                            m.drawsecondlayer=viewcellsflag;
                            PaperType='A1';
                            %set(figs(end),'PaperType',PaperType);
                            %morphogen=m.plotdefaults.morphogen;
                            % temp force colours and morphogens

                            %                     m.mgenposcolors(:,[21,40,42,11])=[
                            %                         0.8706    1.0000         0 0.6
                            %                         0.4902         0    1.0000 0.2
                            %                         0    1.0000         0.5000  0];
                            %lip, uptube, palate dich
                            viewcellsflag=true; %false;
                            m = leaf_plot( m ,...
                                'clipbymgen',false,...
                                'doclip',false,...
                                'autozoom',true,...
                                'autocentre',true,...
                                'autoScale',true,...
                                'azimuth',0,...
                                'elevation',90,...
                                'morphogen',morphogen,...
                                'bgcolor', [1 1 1],...
                                'drawedges', drawedges,... % 0 or 2
                                'thick', thick,... % 0 or 2
                                'drawsecondlayer', viewcellsflag,... % 0 or 2
                                'drawgradients', drawgradients);%... % Boolean
                            pichandles = guidata( m.pictures(1) );
                            axes(pichandles.picture);
                            drawnow
                            azimuth=m.plotdefaults.ourViewParams.azimuth;
                            elevation=m.plotdefaults.ourViewParams.elevation;
                            %view(m.plotdefaults.ourViewParams.azimuth,m.plotdefaults.ourViewParams.elevation);
                            zoom(0.8);
                            drawnow
                        else
                            if isempty(azimuth)
                                %                            morphogen={...
                                %                                   'ID_DIV',...
                                %                                   'ID_LIP',...
                                %                                   'ID_PALATE',...
                                %                                   'ID_TUBE'};
                                %morphogen=m.plotdefaults.morphogen;
                                m = leaf_plot( m ,...
                                    'clipbymgen',false,...
                                    'doclip',false,...
                                    'autozoom',true,...
                                    'autocentre',true,...
                                    'autoScale',true,...
                                    'azimuth',43.8961,...
                                    'elevation',44.7321,... %33
                                    'morphogen',morphogen,...
                                    ...'morphogen',viewcellsflag,...
                                    'drawedges', drawedges,... % 0 or 2
                                    'thick', thick,... % 0 or 2
                                    'bgcolor', [1 1 1],...
                                    'drawsecondlayer', drawsecondlayer,... % 0 or 2
                                    'drawgradients', drawgradients);%... % Boolean
                                pichandles = guidata( m.pictures(1) );
                                axes(pichandles.picture);
                                PaperType='A4';
                            else
                                m = leaf_plot( m ,...
                                    'clipbymgen',false,...
                                    'doclip',false,...
                                    'autozoom',true,...
                                    'autocentre',true,...
                                    'autoScale',true,...
                                    'azimuth',azimuth,...
                                    'elevation',elevation,...
                                    'morphogen',morphogen,...
                                    'bgcolor', [1 1 1],...
                                    'drawedges', drawedges,... % 0 or 2
                                    'thick', thick,... % 0 or 2
                                    'drawsecondlayer', drawsecondlayer,... % 0 or 2
                                    'drawgradients', drawgradients);%... % Boolean
                                pichandles = guidata( m.pictures(1) );
                                axes(pichandles.picture);
                                azimuth=m.plotdefaults.ourViewParams.azimuth;
                                elevation=m.plotdefaults.ourViewParams.elevation;
                                view(m.plotdefaults.ourViewParams.azimuth,m.plotdefaults.ourViewParams.elevation);
                                zoom(0.8);
                                PaperType='A4';
                            end
                        end
                        figs(end+1)=get(get(pichandles.picture,'parent'),'parent');%(gcf;
                        set(figs(end),'Position',position);%,'visible','off');
                        set(figs(end),'PaperPosition',paperposition);
                        %set(figs(end),'PaperType',PaperType);
                        if invisible
                            set(gcf,'visible','off');
                        end
                        %             drawnow
                        %             shg
                        %                     data.current_axes(i)=m.pichandles.picture;%mainaxis(1);
                        %                     ch=pichandles.picture;
                        ch=m.pictures;
                        %                     ch=get(plotaxis(8),'Children');
                        %                     if ~ishandle(ch(1))
                        %                         keyboard
                        %                     end
                        %                     axes(ch(1));
                        %                     camlight('right')
                        %                 camlight('right')

                        ax=axis;
                        %axis(ax.*1.1);
                        titlestring=([experiment,' (t=',num2str(m.globalDynamicProps.currenttime),...
                            '  step',num2str(m.globalDynamicProps.currentIter),')']);
                        disp(titlestring);
                        cameratoolbar(gcf);
                        ha=axes('position',[0.03,0.8,0.9,0.08]);
                        dy=0.16;
                        y=0.9;
                        text(0.01,y,titlestring,'interpreter','none','fontsize',10);
                        y=y-dy;
                        text(0.01,y,file,'interpreter','none','fontsize',10);
                        y=y-dy;
                        cs={};
                        tempcs=commandstring;
                        while length(tempcs)>0
                            if length(tempcs)>100
                                cs{end+1}=tempcs(1:100);
                                tempcs=tempcs(101:end);
                            else
                                cs{end+1}=tempcs;
                                tempcs='';
                            end
                        end
                        for kk=1:length(cs)
                            text(0.01,y,cs{kk},'interpreter','none','fontsize',10);
                            y=y-dy;
                        end
                        axis off
                        if pngFlag
                            %                         set(gcf,'PaperPositionMode','auto')
                            %mainfig=get(get(ch,'parent'),'parent');
                            %could try ancestor
                            mainfig=ch;
                            figure(gcf);
                            set(mainfig,'PaperPositionMode','auto');
                            %                         print('-f',mainfig,'-dpng',printfilename
                            % Added 18/05/2011
                            s=sprintf('print -f%-1d -dpng ''%s''',mainfig,printfilename);

%                             if flatten
%                                 s=sprintf('print -r800 -f%-1d -dpng ''%s''',mainfig,printfilename);
%                             else
%                                 s=sprintf('print -r180 -f%-1d -dpng ''%s''',mainfig,printfilename);
%                             end
                            %s=sprintf('print -f%-1d -dpng %s',mainfig,printfilename);
                            eval(s);
                            printflag=true; % add print here
                            figurefilename=[printfilename(1:end-4),'.fig'];
                            %                         saveas(gcf,figurefilename,'fig');
                        end
                        if printflag
                            %                     print
                        end
%                         ind=find(figs==gcf);
%                         figs(ind)=[];
%                                             %close gcf
%                         close(ch);
                    else
                        disp('unable to read file');
                    end
                end
            end
        end
        if ~isempty(prefix)
            printcommandfilename=fullfile(pathname,projectname,'runs',experiment,'meshes',[prefix,'commandline.png']);
            % Comes here without figs if there are no meshes to print
            if isempty(figs)
                h = figure;
            end
            PrintCommandString(commandstring,position,paperposition,printcommandfilename);
            % If there were no meshes, the figure cannot be reused
            if isempty(figs)
                close(h);
            end
        end

        %         if printflag
        %         end
    end
end

function PrintCommandString(commandstring,position,paperposition,printcommandfilename)
    %     if exist(printcommandfilename)~=2
    % Get the current figure and hide the mesh if it exists
    fig=gcf; % figure
    picturepanel = findobj(gcf,'tag','picturepanel');
    if ~isempty(picturepanel)
        set(picturepanel,'visible','off');
    end
    set(fig,'Position',position,'paperposition',paperposition);
    ind=findstr(commandstring,',');
    str=commandstring;
    startind=1;
    endind=ind(1);
    cla
    axis off
    set(gca,'position',[0.01,0.01,0.99,0.99]);
    top=0.9;
    numberlines=20;
    linespacing=top/numberlines;
    toptemp=top;
    k=1;
    for i=2:2:length(ind)
        s1=str(startind:endind-1);
        startind=endind+1;
        endind=ind(i);
        text(0,toptemp,s1,'interpreter','none','fontsize',35);
        toptemp=(top-linespacing*k);
        k=k+1;
    end
%     print('-dpng',printcommandfilename);
    s=sprintf('print -f%-1d -dpng ''%s''',fig,printcommandfilename); % force correct figure to be printed
    eval(s);
    set(picturepanel,'visible','on');
    cla
    %close(fig);
    %     end
end

function   figs=FlyleafFigs(pathname,projectname,experiment,pngFlag,printflag,invisible)
    if ishandle(1)
        old_fig=gcf;
    else
        old_fig=[];
    end
    figind=1;
    filename=fullfile(pathname,projectname,'runs',experiment,'CommandLine.txt');
    h=fopen(filename,'r');
    if h==-1
        error('could not find CommandLine.txt in runs subdirectory');
    end
    commandline=fgetl(h);
    fclose(h);
    alltext={};
    filename2=fullfile(pathname,projectname,'runs',experiment,[projectname,'.txt']);
    if exist(filename2)==2
        h=fopen(filename2,'r');
        i=1;
        while ~feof(h)
            s=fgetl(h);
            alltext{end+1}=sprintf('%-10d%-100s',i,s);
            i=i+1;
        end
        fclose(h);
    end
    number_lines=70;
    lineheight=1/number_lines;
    figind=1;
    jj=1;
    textlen=length(alltext);
    linenumber=number_lines;
    figs(figind)=figure;
    %             if invisible
    %                 set(gcf,'Visible','off');
    %             end
    %     set(figs(figind),'Position',[7   523   864   998]);%,'visible','off');
    ss=get(0,'screensize');
    set(figs(end),'Position',[7   ss(4)-998-100   864   998]);%,'visible','off');
    set(figs(figind),'PaperPosition',[0.634518 1.17743 20.0 27.5]);
    axis([0,1,0,1],'off');
    set(gca,'Position',[0.0300    0.050    0.9750    0.93]);
    text(0,1,projectname,'FontSize',15,'FontWeight','bold','interpreter','none');
    linenumber=linenumber-3;
    text(0,linenumber*lineheight,commandline,'interpreter','none','FontWeight','bold');
    linenumber=linenumber-3;
    while jj<=textlen
        text(0,linenumber*lineheight,sprintf('%-s',alltext{jj}),'FontSize',10,...
            'FontWeight','normal','interpreter','none');
        jj=jj+1;
        linenumber=linenumber-1;
        if linenumber<=0
            linenumber=number_lines;
            %             printfilename=sprintf('%s%d.png',projectname,figind);
            %             if pngFlag
            %                 print('-dpng',printfilename);
            %             end
            figind=figind+1;
            figs(figind)=figure;
            if invisible
                set(gcf,'Visible','off');
            end
            %             set(figs(figind),'Position',[7   523   864   998]);%,'visible','off');
            ss=get(0,'screensize');
            set(figs(end),'Position',[7   ss(4)-998-100   864   998]);%,'visible','off');
            set(figs(figind),'PaperPosition',[0.634518 1.17743 20.3046 27.5]);
            axis([0,1,0,1],'off');
            set(gca,'Position',[0.0300    0.05    0.9750    0.93]);
            %                 set(gca,'Position',[0.1300    0.1100    0.7750    0.8150]);
        end
    end
    name=sprintf('%s%d.png',projectname,figind);
    printfilename=fullfile(pathname,projectname,'runs',experiment,name);
    %     if pngFlag
    %         print('-dpng',printfilename);
    %     end
    if ~isempty(old_fig)
        figure(old_fig);
    end
    close all
end

function hlink=link_all_the_figures(current_axes)
    list=[];
    for i=1:length(current_axes)
        if ishandle(current_axes(i))
            list(end+1)=current_axes(i);
        end
    end
    % set(gcf,'MenuBar','figure');
    list=list(:)';
    if length(current_axes)>1
        linkaxes(list);
        hlink=linkprop(list,{'CameraPosition','CameraUpVector'});
    else
        cameratoolbar(gcf);
    end
end

%------ what follows is stereo view and trackball control ------


function Initialise_stereoPairAxes( ax1, ax2, vergence, trackball, rotate_check_handle )
    %Entirely based on JRK stereoPairAxes
    %
    %Initialise_stereoPairAxes( ax1, ax2, vergence, trackball )
    % ax1, ax2 axes to be made a stereo pair
    % vergence, typically 1/12
    % trackball, default is false, else true and trackball will install
    %            i.e. left click and drag to rotate axes
    %                 right click to cancel rotate until next initialisation
    %
    %stereoPairAxes( ax1, ax2, vergence, trackball )
    %   Make the two axes objects into a stereo trackball pair, by storing in
    %   the guidata of each a reference to the other and the angle of
    %   convergence.  ax1 is assumed to be the left eye view and ax2 the right
    %   eye view, in which case vergence should be positive.
    %
    %   The trackball argument is optional and defaults to false.  If true,
    %   trackballButtonDownFcn is installed as the ButtonDownFcn in both axes
    %   objects and all their children.  Clicking and dragging on either
    %   window will rotate both views together.  If you only want trackball
    %   functionality part of the time, then you should not install it here, but
    %   call trackballButtonDownFcn from your own ButtonDownFcn when you want
    %   trackball behaviour.
    if nargin < 4
        trackball = false;
    end
    gd1 = guidata( ax1 );
    old_ButtonDownFcn=get(ax1,'ButtonDownFcn');
    child_old_ButtonDownFcn=get(get(ax1,'Children'),'ButtonDownFcn');
    gd1.stereodata = struct( 'otheraxes', ax2, 'vergence', vergence );
    gd1.old_ButtonDownFcn=old_ButtonDownFcn;
    gd1.child_old_ButtonDownFcn=child_old_ButtonDownFcn;
    gd1.rotate_check_handle=rotate_check_handle;
    guidata( ax1, gd1 );
    if trackball
        set( ax1, 'ButtonDownFcn', @trackballButtonDownFcn );
        set( get(ax1,'Children'), 'ButtonDownFcn', @trackballButtonDownFcn );
    end
    if ~isempty(ax2)
        gd2 = guidata( ax2 );
        old_ButtonDownFcn=get(ax2,'ButtonDownFcn');
        child_old_ButtonDownFcn=get(get(ax2,'Children'),'ButtonDownFcn');
        gd2.stereodata = struct( 'otheraxes', ax1, 'vergence', -vergence );
        gd2.old_ButtonDownFcn=old_ButtonDownFcn;
        gd2.child_old_ButtonDownFcn=child_old_ButtonDownFcn;
        gd2.rotate_check_handle=rotate_check_handle;
        guidata( ax2, gd2 );
        if trackball
            set( ax2, 'ButtonDownFcn', @trackballButtonDownFcn );
            set( get(ax2,'Children'), 'ButtonDownFcn', @trackballButtonDownFcn );
        end
    end
end

function trackballButtonDownFcn( hObject, eventData, cancel_trackball )
    %trackballButtonDownFcn( hObject, eventData, cancel_trackball)
    %   Installing this as the ButtonDownFcn for an axes object or its children
    %   will cause a click on the object to produce trackball-like behaviour
    %   when the mouse is dragged.
    %
    %   The trackball feature is cancelled either by clicking the right button
    %   or by making a call with hObject set to the left (or only) axis
    %   and cancel_trackball set to true.

    %   This routine temporarily stores all the information it needs to process the
    %   mouse movement, in the field 'trackballData' of the guidata of the figure.
    %   The field is deleted when the mouse button is released.
    %   As a side-effect, all of the camera mode properties of the axes are set
    %   to manual.
    %
    %   There is more than one way to map mouse movement to rotation, and the
    %   choice is a matter of design.  Three mappings are implemented, and the choice
    %   can be specified by setting the 'trackballMode' field of guidata.
    %   'global' (the default) maps trackball coordinates [x,y] to a rotation about
    %   an axis perpendicular to that vector in the (cameraright,cameraup) plane, by
    %   an amount proportional to the length of the vector.
    %   'local' maps incremental motion of the mouse to incremental rotation in
    %   the same way.
    %   'upright' maps trackball coordinates [x,y] to azimuth and elevation,
    %   maintaining the camera up vector equal to [0 0 1].  This is equivalent
    %   to the trackball mode available from the standard camera toolbar,
    %   except that the scaling may be different.
    %
    %   Global and Upright mode have the advantage that mouse positions are uniquely
    %   mapped to rotations: moving the mouse back to its starting point
    %   restores the axes to its original orientation.
    %   Local mode has the advantage that the current rotation axis is always
    %   perpendicular to the current mouse velocity.
    %   There is no mapping of trackball data to axes rotation that allows both
    %   of these.
    %
    %   The scaling factor is currently set so that movement by a distance
    %   equal to the smaller of the width and height of the axes object rotates
    %   it by 2 radians.  A different value can be set by setting the
    %   'trackballScale' field of the figure's guidata to the requried value.

    % Find the axes.
    if nargin<3
        cancel_trackball=false
    end
    while ~(isempty(hObject)) && ishandle(hObject) && ~strcmp(get(hObject,'Type'), 'axes')
        hObject = get( hObject, 'Parent' );
    end
    if ~ishandle(hObject), return; end

    % Find the figure.
    hFigure = get( hObject, 'Parent' );
    while ishandle(hFigure) && ~strcmp(get(hFigure,'Type'), 'figure')
        hFigure = get( hFigure, 'Parent' );
    end
    if ~ishandle(hFigure), return; end
    if cancel_trackball || strcmp(get( hFigure,'SelectionType'),'alt')
        % right button click to disable rotate until next time system is
        % initialised
        gd=guidata(hObject);
        ax1=gd.old_ButtonDownFcn;
        ax2=gd.child_old_ButtonDownFcn;
        if ~any(ishandle(ax2))
            ax2='';
        end
        set(gd.rotate_check_handle,'Value',0);
        set( hObject, 'ButtonDownFcn', ax1 );
        set( get(hObject,'Children'), 'ButtonDownFcn', ax2 );
        if isfield( gd, 'stereodata' )
            if ~isempty(gd.stereodata.otheraxes)
                gd2=guidata(gd.stereodata.otheraxes);
                ax1=gd2.old_ButtonDownFcn;
                ax2=gd2.child_old_ButtonDownFcn;
                if ~any(ishandle(ax2))
                    ax2='';
                end
                set(gd2.rotate_check_handle,'Value',0);
                set( gd.stereodata.otheraxes, 'ButtonDownFcn', ax1 );
                set( get(gd.stereodata.otheraxes,'Children'), 'ButtonDownFcn', ax2 );
            end
        end
    else
        % Get all the information we need.
        hitPointParent = get( hFigure, 'CurrentPoint' );
        oldUnits = get( hObject, 'Units' );
        set( hObject, 'Units', 'pixels' );
        axespos = get( hObject, 'Position' );
        set( hObject, 'Units', oldUnits );
        axessize = min(axespos([3 4]));

        cameraTarget = get( hObject, 'CameraTarget' );
        cameraPosition = get( hObject, 'CameraPosition' );
        [cameraLook, cameraUp, cameraRight] = cameraFrame( ...
            cameraPosition, ...
            cameraTarget, ...
            get( hObject, 'CameraUpVector' ) );

        % Store it in guidata.
        gd = guidata(hFigure);
        if isfield( gd, 'trackballScale' )
            trackballScale = gd.trackballScale;
        else
            trackballScale = 2;
        end
        if isfield( gd, 'trackballMode' )
            trackballMode = gd.trackballMode;
        else
            trackballMode = 'global';
        end
        gd.trackballData = struct( 'axes', hObject, ...
            'startpoint', hitPointParent, ...
            'currentpoint', hitPointParent, ...
            'axessize', axessize, ...
            'cameraLook', cameraLook, ...
            'cameraUp', cameraUp, ...
            'cameraRight', cameraRight, ...
            'cameraTarget', cameraTarget, ...
            'cameraPosition', cameraPosition, ...
            'trackballScale', trackballScale, ...
            'trackballMode', trackballMode, ...
            'oldWindowButtonMotionFcn', get( hFigure, 'WindowButtonMotionFcn' ), ...
            'oldWindowButtonUpFcn', get( hFigure, 'WindowButtonUpFcn' ) );
        guidata( hFigure, gd );
        set( hFigure, ...
            'WindowButtonMotionFcn', @trackballButtonMotionFcn, ...
            'WindowButtonUpFcn', @trackballButtonUpFcn );
        set( hObject, ...
            'CameraPositionMode', 'manual', ...
            'CameraTargetMode', 'manual', ...
            'CameraUpVector', cameraUp, ...
            'CameraUpVectorMode', 'manual', ...
            'CameraViewAngleMode', 'manual' );
        if isfield( gd, 'stereodata' )
            stereoTransfer( hObject, gd.stereodata.otheraxes, gd.stereodata.vergence );
        end
    end
end

function trackballButtonMotionFcn( hObject, eventData )
    gd = guidata( hObject );
    if isempty( gd ) || ~isfield( gd, 'trackballData' ), return; end
    currentpoint = get( hObject, 'CurrentPoint' );
    switch gd.trackballData.trackballMode
        case 'global'
            delta = (currentpoint - gd.trackballData.startpoint)/gd.trackballData.axessize;
        otherwise % case { 'local', 'upright' }
            delta = (currentpoint - gd.trackballData.currentpoint)/gd.trackballData.axessize;
    end
    [cameraPos,cameraUp,cameraRight] = ...
        trballView( gd.trackballData, [-delta(1), delta(2)], ...
        gd.trackballData.trackballScale );
    if strcmp( gd.trackballData.trackballMode, 'upright' )
        cameraUp = [0 0 1];
    end
    gd.trackballData.currentpoint = currentpoint;
    switch gd.trackballData.trackballMode
        case 'global'
            %nothing
        otherwise % case { 'local', 'upright' }
            gd.trackballData.cameraPosition = cameraPos;
            gd.trackballData.cameraUp = cameraUp;
            gd.trackballData.cameraRight = cameraRight;
    end
    set( gd.trackballData.axes, ...
        'CameraPosition', cameraPos, ...
        'CameraPositionMode', 'manual', ...
        ... % 'CameraTarget', ??, ...
        'CameraTargetMode', 'manual', ...
        'CameraUpVector', cameraUp, ...
        'CameraUpVectorMode', 'manual', ...
        ... % 'CameraViewAngle', ??, ...
        'CameraViewAngleMode', 'manual' );
    if isfield( gd, 'stereodata' )
        stereoTransfer( gd.stereodata.otheraxes, ...
            cameraPos, ...
            get( gd.trackballData.axes, 'CameraTarget' ), ...
            cameraUp, ...
            gd.stereodata.vergence );
    end
    guidata( hObject, gd );
end

function trackballButtonUpFcn( hObject, eventData )
    % Process the final mouse position.
    % Tests indicate that this is not necessary -- the final mouse position
    % always generates a ButtonMotion event.
    % trackballButtonMotionFcn( hObject, eventData )

    % Restore the original ButtonMotion and ButtonUp functions.
    % Delete the trackball data.
    gd = guidata( hObject );
    if ~isempty(gd) && isfield( gd, 'trackballData' )
        set( hObject, 'WindowButtonMotionFcn', gd.trackballData.oldWindowButtonMotionFcn );
        set( hObject, 'WindowButtonUpFcn', gd.trackballData.oldWindowButtonUpFcn );
        gd = rmfield( gd, 'trackballData' );
        guidata( hObject, gd );
    end
end

function [cameraPos,cameraUp,cameraRight] = trballView( trdata, trball, trballscale )
    %[cameraPos,cameraUp] = trballView( initview, trball, trballscale )
    %   initview is a structure containing
    %       cameraRight
    %       cameraUp
    %   trball is a pair of real numbers indicating the amount of trackball
    %   movement in two dimensions.  These represent a rotation about an axis
    %   in the view plane perpendicular to the trball vector.  The scaling is 1
    %   unit = trballscale radians. cameraPos and cameraUp are set to the new
    %   camera position and up-vector specified by the trackball value.

    if all(trball==0)
        cameraPos = trdata.cameraPosition;
        cameraUp = trdata.cameraUp;
        cameraRight = trdata.cameraRight;
    else
        rotAxis = trball(2)*trdata.cameraRight + trball(1)*trdata.cameraUp;
        rotAxis = rotAxis/norm(rotAxis);
        rotAmount = norm(trball) * trballscale;

        % Rotate cameraPosition about rotAxis through cameraTarget by
        % rotAmount.
        cameraPos = rotVec( trdata.cameraPosition, trdata.cameraTarget, rotAxis, rotAmount );

        % Rotate cameraPosition about rotAxis by rotAmount.
        cameraUp = rotVec( trdata.cameraUp, [0 0 0], rotAxis, rotAmount );

        % Rotate cameraRight about rotAxis by rotAmount.
        cameraRight = rotVec( trdata.cameraRight, [0 0 0], rotAxis, rotAmount );
    end
end

function [CameraPosition, CameraTarget, CameraUpVector] = stereoTransfer( varargin )
    %stereoTransfer( varargin )
    %   Set the viewpoint of one axes to that of another, rotated about the
    %   camera up vector.
    %
    %   stereoTransfer(AX) will transfer the view of the axes AX to the axes
    %       contained in AX's guidata element stereodata.otheraxes, rotated by
    %       stereodata.vergence (an angle in radians).
    %   stereoTransfer() is equivalent to stereoTransfer(gca).
    %   stereoTransfer( AX1, AX2, VERGENCE ) will transfer the view from AX1 to
    %       AX2, offset by the angle VERGENCE.
    %   stereoTransfer( AX, POSITION, TARGET, UP, VERGENCE ) will set the view
    %       of the axes AX to the result of offsetting the view by a rotation
    %       about UP by VERGENCE, from the view specified by the given camera
    %       POSITION, TARGET, and UP vectors.
    %
    %   In all cases, vergence must be an angle in radians.  If the view is
    %   being transferred from the left eye to the right eye, vergence should
    %   be positive.  If the distance between someone's eye is E and they are
    %   a distance D from the screen, vergence should be E/D radians.
    %
    %   The resulting view is returned in the output arguments.

    error(nargchk(0, 5, nargin, 'struct'));
    theaxes = [];
    switch nargin
        case { 0, 1 }
            if nargin==0
                theaxes = gca;
            else
                theaxes = varargin{1};
            end
            CameraPosition = get( theaxes,'CameraPosition' );
            CameraTarget = get( theaxes,'CameraTarget' );
            CameraUpVector = get( theaxes,'CameraUpVector' );
            gd = guidata( theaxes );
            if isempty(gd), return; end
            if ~isfield( gd, 'stereodata' ), return; end
            sd = gd.stereodata;
            otheraxes = sd.otheraxes;
            vergence = sd.vergence;
        case 3
            theaxes = varargin{1};
            CameraPosition = get( theaxes,'CameraPosition' );
            CameraTarget = get( theaxes,'CameraTarget' );
            CameraUpVector = get( theaxes,'CameraUpVector' );
            otheraxes = varargin{2};
            vergence = varargin{3};
        case 5
            otheraxes = varargin{1};
            CameraPosition = varargin{2};
            CameraTarget = varargin{3};
            CameraUpVector = varargin{4};
            vergence = varargin{5};
        otherwise
            error( 'stereoTransfer: wrong number of arguments, 3 or 5 expected.' );
    end
    CameraUpVector = makeperp( CameraTarget-CameraPosition, CameraUpVector );
    CameraUpVector = CameraUpVector/norm(CameraUpVector);

    % Rotate Position about Up by vergence
    CameraPosition = rotVec( CameraPosition, CameraTarget, CameraUpVector, vergence );
    set( otheraxes, ...
        'CameraPosition', CameraPosition, ...
        'CameraTarget', CameraTarget, ...
        'CameraUpVector', CameraUpVector, ...
        'CameraPositionMode', 'manual', ...
        'CameraTargetMode', 'manual', ...
        'CameraUpVectorMode', 'manual', ...
        'CameraViewAngleMode', 'manual' );
    if theaxes
        set( theaxes, 'CameraUpVector', CameraUpVector );
        set( otheraxes, 'CameraViewAngle', get( theaxes,'CameraViewAngle' ) );
    end
end

function stereoPairAxes( ax1, ax2, vergence, trackball )
    %stereoPairAxes( ax1, ax2, vergence, trackball )
    %   Make the two axes objects into a stereo trackball pair, by storing in
    %   the guidata of each a reference to the other and the angle of
    %   convergence.  ax1 is assumed to be the left eye view and ax2 the right
    %   eye view, in which case vergence should be positive.
    %
    %   The trackball argument is optional and defaults to false.  If true,
    %   trackballButtonDownFcn is installed as the ButtonDownFcn in both axes
    %   objects and all their children.  Clicking and dragging on either
    %   window will rotate both views together.  If you only want trackball
    %   functionality part of the time, then you should not install it here, but
    %   call trackballButtonDownFcn from your own ButtonDownFcn when you want
    %   trackball behaviour.

    if nargin < 4
        trackball = false;
    end
    gd1 = guidata( ax1 );
    ax1_old_ButtonDownFcn=get(ax1,'ButtonDownFcn');
    child_old_ButtonDownFcn=get(get(ax1,'Children'),'ButtonDownFcn');
    gd1.stereodata = struct( 'otheraxes', ax2, 'vergence', vergence );
    gd1.old_ButtonDownFcn=old_ButtonDownFcn;
    gd1.child_old_ButtonDownFcn=child_old_ButtonDownFcn;
    guidata( ax1, gd1 );
    gd2 = guidata( ax2 );
    old_ButtonDownFcn=get(ax2,'ButtonDownFcn');
    child_old_ButtonDownFcn=get(get(ax2,'Children'),'ButtonDownFcn');
    gd2.stereodata = struct( 'otheraxes', ax1, 'vergence', -vergence );
    gd2.old_ButtonDownFcn=old_ButtonDownFcn;
    gd2.child_old_ButtonDownFcn=child_old_ButtonDownFcn;
    guidata( ax2, gd2 );
    if trackball
        set( ax1, 'ButtonDownFcn', @trackballButtonDownFcn );
        set( get(ax1,'Children'), 'ButtonDownFcn', @trackballButtonDownFcn );
        set( ax2, 'ButtonDownFcn', @trackballButtonDownFcn );
        set( get(ax2,'Children'), 'ButtonDownFcn', @trackballButtonDownFcn );
    end
end

function [cameraPosition, cameraRight] = stereoOffset( cameraPosition, cameraTarget, cameraUp, cameraRight, offset )
    %[cameraPosition, cameraRight] = ...
    %   stereoOffset( cameraPosition, cameraTarget, cameraUp, cameraRight, offset )
    % Rotate cameraPosition and cameraRight by offset about cameraUp through cameraTarget.
    % offset is in degrees.  The vectors should be row vectors.

    m = matRot( cameraUp, offset*pi/180 )';
    cameraPosition = (cameraPosition - cameraTarget)*m + cameraTarget;
    cameraRight = cameraRight*m;
end

function v = rotVec( v, centre, rotAxis, rotAmount )
    %v = rotVec( v, centre, rotAxis, rotAmount )
    %   Rotate v about the axis rotAxis passing through centre by rotAmount in
    %   radians.
    %   v may be a 3-element row vector or an N*3 matrix.

    m = matRot( rotAxis, rotAmount );
    v = (v-centre)*m'+centre;
end

function m = matRot( rotAxis, rotAmount )
    %m = matRot( rotAxis, rotAmount )
    %   Calculate the rotation matrix that rotates about the given axis by the
    %   given amount in radians.  The axis must be a unit vector.
    %   Formula cribbed from Wikipedia.

    x = rotAxis(1);  y = rotAxis(2);  z = rotAxis(3);
    c = cos(rotAmount);  s = sin(rotAmount);  C = 1-c;
    xs = x*s; ys = y*s; zs = z*s;
    xC = x*C; yC = y*C; zC = z*C;
    xyC = x*yC; yzC = y*zC; zxC = z*xC;
    m = [ [ x*xC+c, xyC-zs, zxC+ys ]; ...
        [ xyC+zs, y*yC+c, yzC-xs ]; ...
        [ zxC-ys, yzC+xs, z*zC+c ] ];
end

function v2p = makeperp( v1, v2 )
    %v2p = makeperp( v1, v2 )
    %   Construct the vector v2p which is orthogonal to v1 and in the same plane
    %   as v2.

    v2p = v2 - v1*dot(v2,v1)/dot(v1,v1);
end

function [cameraLook, cameraUp, cameraRight] = cameraFrame( CameraPosition, CameraTarget, CameraUp )
    %[cameraLook, cameraUp, cameraRight] = cameraFrame( ...
    %       CameraPosition, CameraTarget, CameraUp )
    %   Construct the camera orthonormal frame of reference.  The arguments are
    %   fields of an axes object.
    %   If called as [cameraLook, cameraUp, cameraRight] = cameraFrame( axes ),
    %   the position, target, and up vectors are taken from the axes object.
    %   Calling it with no arguments is equivalent to cameraFrame( gca ).
    %   Note that cameraUp will not necessarily coincide with CameraUp, since
    %   cameraUp is constrained to be orthogonal to the look direction, while
    %   CameraUp is not.

    if nargin <= 1
        if nargin==0
            theaxes = gca;
        else
            theaxes = CameraPosition;
        end
        CameraPosition = get( theaxes,'CameraPosition' );
        CameraTarget = get( theaxes,'CameraTarget' );
        CameraUp = get( theaxes,'CameraUpVector' );
    end

    cameraLook = CameraTarget - CameraPosition;
    cameraLook = cameraLook/norm(cameraLook);

    cameraUp = makeperp( cameraLook, CameraUp );
    cameraUp = cameraUp/norm(cameraUp);

    cameraRight = cross( cameraLook, cameraUp );
end

