function ViewGrowthTool(arg,option)
    % ViewGrowthTool(arg,option)
    % view a bud mesh model
    %
    % arg, if 's' then force a stereo display
    %      otherwise it should automatically produce a 
    %      stereo display if it detects two monitors
    %
    %ViewGrowthTool
    %
    %J. Andrew Bangham, 2008

    if nargin<2
        option=0;
    end
    global FigureName % keep handle 
    if nargin<1 & ishandle(FigureName)
        delete(FigureName); % Delete existing 
        close all
    end
    isstereo=false;
    fromaamtoolbox_flag=false;
    Facets=[];
    if nargin>0
        if strcmp(lower(arg(1)),'s')
            isstereo=true;
            flag=true;
        elseif strcmp(lower(arg),'fromaamtoolbox')
            % assume it is to be in stereo of two monitors
            Facets=option;
            m=get(0,'MonitorPositions');
            if size(m,1)==2
                isstereo=true;
            end
            fromaamtoolbox_flag=true;
            flag=true;
        else
            isstereo=false;
            flag=false;
        end
    else
        flag=true;
    end
    if flag % then initialise nargin<1 | ~strcmp(lower(arg(1)),'s')
        if ~isempty(FigureName)
            if ishandle(FigureName )
                close(FigureName);
            end
        end
        IndexOfReferenceBud=1;%10;
        cwd=pwd;
        if ~exist('PointModels','dir')
            % then we are not within an appropriate project folder so move
            ok=false;
            while ~ok
                directory=uigetdir(pwd,'Please select a shape model project directory e.g. PRJ_*');
                [pathstr, name, ext, versn] = fileparts(directory);
                if strcmp(name(1:4),'PRJ_')
                    ok=true;
                end
            end
            cd(directory);
        end
        [Data_vertex,data_vertex_path]=uigetfile('*_pm.mat',...
            'Select example pointmodel from the set',...
            'PointModels\updated_pointmodels_wildtype\907awh9_x18p18y18p18z18p18cop_pm.mat');
        PathToFacets='Templates'; %fullfile(data_vertex_path,'FacetData');
        if ~exist(PathToFacets,'dir')
            mkdir(PathToFacets);
            disp('Please put filename_to_plastochron_mapping.txt file and ')
            disp(sprintf('template_triangles.xls file into %s',PathToFacets));
            uiwait(warndlg(sprintf(...
                'Please put filename_to_plastochron_mapping.txt file and template_triangles.xls file into %s',...
                PathToFacets)));
            return
        else
            d=dir(fullfile(PathToFacets,'*.xls'));
            if length(d)>0
                if isempty(Facets)
                    facets_file=fullfile(PathToFacets,d(1).name);
                else
                    facets_file=Facets.name;
                end
                [facets_file,PathToFacets]=uigetfile('*.xls',...
                    'Select facets_file.xls (facets on Sheet1)',...
                    facets_file);
                sheetname='Sheet1';
            else
                uiwait(warndlg(sprintf(...
                    'Please put template_triangles.xls file into %s',...
                    PathToFacets)));
                return;
            end
            d=dir(fullfile(PathToFacets,'*.txt'));
            if length(d)==0
                uiwait(warndlg(sprintf(...
                    'Please put filename_to_plastochron_mapping.txt file into %s',...
                    PathToFacets)));
                return;
            end
            
        end
        filename_to_plastochron_filename=fullfile(PathToFacets,'filename_to_plastochron_mapping.txt');
        %data_vertex_path='PRJ_wild_type\PointModels\910gfp163_x8p91y8p91z8p91pop_3D';
%         template_file='910gfp163_x8p91y8p91z8p91pop_3D.csv';
%         facets_file='template_triangles_jab.xml';
%         sheetname='Sheet1';
        framerate=20;
        AppData=initialise(IndexOfReferenceBud,data_vertex_path,...
            PathToFacets,facets_file,sheetname,framerate,isstereo,filename_to_plastochron_filename);% ********
        FigureName=AppData.FigureName;
        AppData.isstereo=isstereo;
        guidata(FigureName,AppData) % save the application data
    else
        AppData=guidata(FigureName); % get the application data
        if ~strcmp(arg,'Quit')
            switch arg % none of the cases must start with s
                case 'ChoosePetal'
                    AppData=ChoosePetal(AppData,option);     % ********
                case 'ChoosePart'
                    AppData=ChoosePart(AppData,option);
                case 'VertexNumbers'
                    AppData=VertexNumbers(AppData,option);
                case 'ExportObj'
                    AppData=ExportObj(AppData);
                case 'Filename'
                    AppData=Filename(AppData);
                case 'Align'
                    AppData=Align(AppData);
                case 'Clip'
                    AppData=Clip(AppData);
                case 'Grow'
                    AppData=Grow(AppData);
                case 'Rock'
                    AppData=Rock(AppData);
                case 'Rotate'
                    AppData=Rotate(AppData);
                case 'Zero'
                    AppData=Zero(AppData,option);
                case 'RotateXYZ'
                    AppData=RotateXYZ(AppData,option);
                case 'Disparity'
                    AppData=Disparity(AppData);
                case 'SaveAxes'
                    AppData=SaveAxes(AppData);
                case 'FlipUpDown'
                    AppData=FlipUpDown(AppData);
                case 'GrabRect'
                    AppData=GrabRect(AppData);
                case 'ViewAxes'
                    AppData=ViewAxes(AppData);
                case 'PlotAreaGrowth' 
                    AppData=PlotAreaGrowth(AppData);
                case 'MakeRockMovie'
                    AppData=MakeRockMovie(AppData);
                otherwise
                    error('unrecognised argument');
            end
            guidata(FigureName,AppData); % save the application data
        else
            if ishandle(AppData.RightPlotFigure)
                close(AppData.RightPlotFigure);
            end
            if ishandle(AppData.LeftPlotFigure)
                close(AppData.LeftPlotFigure);
            end
            if option==0
%             if ishandle(AppData.FigureName)
                close(AppData.FigureName);
            else
                eval(AppData.FigureCloseRequestFcn)
            end
        end
        %AppData=guidata(FigureName); % get the application data
    end
end

function AppData=ViewAxes(AppData)
    AppData=meshmaker_display(AppData);
end
    
function AppData=PlotAreaGrowth(AppData)
    set(AppData.handles.NormaliseScale,'value',0);
    set(AppData.handles.Grow,'Value',1);
    AppData=Grow(AppData);
    figure(5)
    clf
    subplot(2,2,3);
    plot(AppData.areas.p*10,AppData.areas.area);
    xlabel('hours');
    ylabel('mm^2');
    title('Smoothed growth areas');
    subplot(2,2,4);
    plot(AppData.areas.p*10,log(AppData.areas.area));
    xlabel('hours');
    ylabel('log(mm^2)');
    title('Smoothed growth areas');

    subplot(2,2,1);
    x=fliplr(AppData.areas.p_measured*10);
    y=fliplr(AppData.areas.area_measured);
    plot(x,y,'*r');
    xlabel('hours');
    ylabel('mm^2');
    title('Measured areas');
    subplot(2,2,2);
    plot(x,log(y),'*r');
    xlabel('hours');
    ylabel('log(mm^2)');
    hold on
    Y=polyfit(x,log(y),1);
    plot(x,(Y(2)+Y(1)*x),'-b');
    hold off
    R=Y(1)/log(2);
    D=1/R;
    title(sprintf('D=%3.1f hrs and  R=%3.2f',D,R));    
    subplot(2,2,1)
    hold on
    plot(x,exp(Y(2)+Y(1)*x),'b')
    hold off
end

function area=find_triangle_area(pts,verbose)
    if nargin<1
%         pts = [0 0;0 1;0.5 0.5];
        pts = [0 0 0;0 1 0.5;0.5 0.5 0];
    end
    if nargin<2
        verbose=false;
    end
    [q,r,e] = qr((pts(2:3,:) - repmat(pts(1,:),2,1))');
    area=abs(prod(diag(r)))/2;
    if verbose
        disp(['Area=',num2str(area)]);
        if size(pts,2)==2
            tri=patch(pts(:,1),pts(:,2),'r');
        elseif size(pts,2)==3
            tri=patch(pts(:,1),pts(:,2),pts(:,3),'r');
        end
    end
end

function AppData=GrabRect(AppData)
    set(gco,'BackgroundColor',[1,0.7,0.7],'string','RectSet');
    f=getframe(AppData.handles.axes1);
    figure(5)
    imshow(f.cdata,'border','tight');
    title('drag a rectangle over area needed for movie');
    disp('drag a rectangle over area needed for movie');
    beep
    AppData.rect = round(getrect(gca));
    f.cdata=f.cdata(AppData.rect(2):AppData.rect(2)+AppData.rect(4),...
        AppData.rect(1):AppData.rect(1)+AppData.rect(3),:);
    imshow(f.cdata);
    drawnow
    pause(1)
    close 5
end

function AppData=SaveAxes(AppData)
    
end

function AppData=FlipUpDown(AppData)
    i=get(AppData.handles.Filename,'Value');
    AppData=meshmaker_display(AppData,AppData.file_indexes(i));
    AppData=ChoosePetal(AppData,0);
end

function AppData=Disparity(AppData)
    disparity=get(AppData.handles.Disparity,'Value');
    plotaxes(1)=AppData.handles.axes1;
    plotaxes(2)=AppData.handles.axes2;
    stereoTransfer( plotaxes(2), ...
        get( plotaxes(1),'CameraPosition'), ...
        get( plotaxes(1),'CameraTarget'), ...
        get( plotaxes(1),'CameraUpVector'), ...
        disparity );
end

function AppData=RotateXYZ(AppData,option)
    Ramp=1.0;
    if option==0
        v=get(AppData.handles.RotX,'Value');
        amp=Ramp*sign(v);
        camorbit(AppData.handles.axes1,amp,0,'data',[0 0 1]);
        if ~isempty(AppData.handles.axes2)
            camorbit(AppData.handles.axes2,amp,0,'data',[0 0 1]);
        end
        set(AppData.handles.RotX,'Value',0)
    elseif option==1
        v=get(AppData.handles.RotY,'Value');
        amp=Ramp*sign(v);
        camorbit(AppData.handles.axes1,0,amp,'data',[0 0 1]);
        if ~isempty(AppData.handles.axes2)
            camorbit(AppData.handles.axes2,0,amp,'data',[0 0 1]);
        end
        set(AppData.handles.RotY,'Value',0)
    else
        v=get(AppData.handles.RotZ,'Value');
        amp=Ramp*sign(v);
        camorbit(AppData.handles.axes1,0,amp,'data',[0 0 1]);
        if ~isempty(AppData.handles.axes2)
            camorbit(AppData.handles.axes2,0,amp,'data',[0 0 1]);
        end
        set(AppData.handles.RotZ,'Value',0)
    end
    drawnow
end


function AppData=Rotate(AppData)
    last_point=[];
    leftaxis=AppData.handles.axes1;
    rightaxis=AppData.handles.axes2;
    disp('Rotate')
    if get(gco,'Value')
        vergence=get(AppData.handles.Disparity,'Value'); %1/12;
        trackball=true;
        %set mouse button callbacks ... WindowButtonDownFcn which should
        Initialise_stereoPairAxes( leftaxis, rightaxis, vergence, trackball, gco );
    else
        %clear mouse button callbacks
        cancel_trackball=true;
        eventData=[];
        trackballButtonDownFcn( leftaxis, eventData,  cancel_trackball);
    end
end

function AppData=Zero(AppData,option) % 38 40
%     Disparity=get(AppData.handles.Disparity,'Value');
    plotaxes(1)=AppData.handles.axes1;
    if ~isempty(AppData.handles.axes2)
        plotaxes(2)=AppData.handles.axes2;
    end
    AppData.axis=[-AppData.MaxLim/3 AppData.MaxLim/3 -AppData.MaxLim/3 AppData.MaxLim/3 -(+1000) (AppData.MaxLim-0)];
    if option==0
        set(AppData.handles.axes1,'View',[90 0]);
    else
        set(AppData.handles.axes1,'View',[38 40]);
    end
    AppData.axis=[-AppData.MaxLim/3 AppData.MaxLim/3 -AppData.MaxLim/3 AppData.MaxLim/3 -(+1000) (AppData.MaxLim-0)];
    for ii=1:length(plotaxes)
        set(plotaxes(ii),'XGrid','off','YGrid','off','ZGrid','off','View',[90 0],...
                    'XLim',AppData.axis(1:2),'YLim',AppData.axis(1:2),'ZLim',AppData.axis(5:6))
        if ii==1
            disparity=0; %get(AppData.handles.Disparity,'Value');
        elseif ii==2
            disparity=get(AppData.handles.Disparity,'Value');
        end
        stereoTransfer( plotaxes(ii), ...
            get( plotaxes(1),'CameraPosition'), ...
            get( plotaxes(1),'CameraTarget'), ...
            get( plotaxes(1),'CameraUpVector'), ...
            disparity );
    end
    
%     if option==0
%         set(AppData.handles.axes1,'View',[90 0]);
%         if ~isempty(AppData.handles.axes2)
%             set(AppData.handles.axes2,'View',[90 0+Disparity]);
%         end
%     else
%         set(AppData.handles.axes1,'View',[38 40]);
%         if ~isempty(AppData.handles.axes2)
%             set(AppData.handles.axes2,'View',[38 40+Disparity]);
%         end
%     end
end

function AppData=Rock(AppData)
    if get(AppData.handles.Rock,'Value')
        AppData.RockDone=false;
        movie_i=1;
        tic;
        i=0;
        j=0;
        di=0.5; %0.07;
        dj=di*0.95;
        Ramp=10; %0.1;
        AlphaDorsal=1;
        AlphaLateral=1;
        axis(AppData.handles.axes1,'square');
        axis(AppData.handles.axes1,'vis3d');
        if ~isempty(AppData.handles.axes2)
%             axis(AppData.handles.axes2,'square');
            axis(AppData.handles.axes2,'vis3d');
        end
        drawnow;
%         ax=axis(AppData.handles.axes1);
        theaxes=AppData.handles.axes1;
        %[az,el,roll] = getview(theaxes);
        set( theaxes, 'CameraViewAngleMode', 'manual' );
%         axis(AppData.handles.axes1,'square');
%         hp=get(theaxes,'parent');
%         [X,Y,Z]=sphere;
%         hold on
%         hs=mesh(theaxes,X,Y,Z);
        ax=320;
        el=15;
        i=ax;
        %setview(theaxes,az,el,roll);
        while AppData.RockDone==false
            t=toc;
            amp=Ramp*get(AppData.handles.RockAmp,'Value');
            si=i+ax; %rem(i,2)-1;%sin(i)*amp;
            sj=el; %cos(j)*amp/100;
            %camorbit(AppData.handles.axes1,si,sj,'data',[0 0 1]);
            %setview(theaxes,si,el,roll);
            view(theaxes,si,sj);
            set( theaxes, 'CameraViewAngleMode', 'manual' );
            set( theaxes, 'DataAspectRatio', [1 1 1] );
%             axis(AppData.handles.axes1,'square');
%             if ~isempty(AppData.handles.axes2)
%                 camorbit(AppData.handles.axes2,si,sj,'data',[0 0 1]);
%             end
            i=i+di;
            j=j+dj;
            %axis(AppData.handles.axes1,ax);
            drawnow;
            if get(AppData.handles.MakeRockMovie,'value')==1 && ~AppData.MakingRockMovie
                % start a movie
                AppData.MakingRockMovie=true;
                f=getframe(get(AppData.handles.axes1,'parent'));
%                 f=getAframe(AppData);
                AppData.MovieSize=size(f.cdata);
                AppData.mov(1)=f;%getframe(AppData.handles.axes1);
            elseif get(AppData.handles.MakeRockMovie,'value')==0 && AppData.MakingRockMovie
                % stop making movie
                AppData.MakingRockMovie=false;
                % and write it out
                set(AppData.handles.Movie,'Value',0);
                moviename=['rotate-',datestr(now),'.avi'];
                ind=findstr(moviename,' ');
                moviename(ind)='T';
                ind=findstr(moviename,':');
                moviename(ind)='-';
                moviename=fullfile(AppData.DataFacetsPath,moviename);
                disp(sprintf('\nPlease wait as writing %s will take some time',moviename));
                movie2avi(AppData.mov(2:end-1),moviename,'fps',AppData.framerate,'Quality',90,'compression','None');
                disp(sprintf('Finished writing %s\n',moviename));
            elseif get(AppData.handles.MakeRockMovie,'value')==1 && AppData.MakingRockMovie
                % continue making movie
%                 f=getAframe(AppData);
                %f=getframe(gcf);%AppData.handles.axes1);
                f=getframe(get(AppData.handles.axes1,'parent'));
                size(f)
                f.cdata=imresize(f.cdata,AppData.MovieSize(1:2));
                %size(f.cdata);
                AppData.mov(end+1)=f;%getframe(AppData.handles.axes1);
            end
            while t<0.2 %0.002
                t=toc;
            end
            AppData.RockDone=~get(AppData.handles.Rock,'Value');
            tic;
        end
    else
        AppData.RockDone=true;
        if get(AppData.handles.MakeRockMovie,'value')==0 && AppData.MakingRockMovie
            % stop making movie
            AppData.MakingRockMovie=false;
            % and write it out
            set(AppData.handles.Movie,'Value',0);
            moviename=['gyrate-',datestr(now),'.avi'];
            ind=findstr(moviename,' ');
            moviename(ind)='T';
            ind=findstr(moviename,':');
            moviename(ind)='-';
            moviename=fullfile(AppData.DataFacetsPath,moviename);
            disp(sprintf('\nPlease wait as writing %s will take some time',moviename));
            movie2avi(AppData.mov(2:end-1),moviename,'fps',AppData.framerate,'Quality',90,'compression','None');
            disp(sprintf('Finished writing %s\n',moviename));
        end
    end
end

function AppData=MakeRockMovie(AppData)
    get(AppData.handles.Rock,'Value')
        AppData.RockDone=false;
        movie_i=1;
        tic;
        i=0;
        j=0;
        di=2; %0.03; %0.07;
        dj=di*0.95;
        Ramp=10; %0.1;
        AlphaDorsal=1;
        AlphaLateral=1;
        while AppData.RockDone==false
            t=toc;
            amp=Ramp*get(AppData.handles.RockAmp,'Value');
            si=sin(i)*amp;
            sj=cos(j)*amp/100;
            camorbit(AppData.handles.axes1,si,sj,'data',[0 0 1]);
            if ~isempty(AppData.handles.axes2)
                camorbit(AppData.handles.axes2,si,sj,'data',[0 0 1]);
            end
            i=i+di;
            j=j+dj;
            drawnow;
            if get(AppData.handles.MakeRockMovie,'value')==1 && ~AppData.MakingRockMovie
                % start a movie
                AppData.MakingRockMovie=true;
                f=getAframe(AppData);
                AppData.MovieSize=size(f.cdata);
                AppData.mov(1)=f;%getframe(AppData.handles.axes1);
            elseif get(AppData.handles.MakeRockMovie,'value')==0 && AppData.MakingRockMovie
                % stop making movie
                AppData.MakingRockMovie=false;
                % and write it out
                set(AppData.handles.Movie,'Value',0);
                moviename=['gyrate-',datestr(now),'.avi'];
                ind=findstr(moviename,' ');
                moviename(ind)='T';
                ind=findstr(moviename,':');
                moviename(ind)='-';
                moviename=fullfile(AppData.DataFacetsPath,moviename);
                disp(sprintf('\nPlease wait as writing %s will take some time',moviename));
                movie2avi(AppData.mov,moviename,'fps',AppData.framerate,'Quality',90,'compression','None');
                disp(sprintf('Finished writing %s\n',moviename));
            elseif get(AppData.handles.MakeRockMovie,'value')==1 && AppData.MakingRockMovie
                % continue making movie
                f=getAframe(AppData);
                f.cdata=imresize(f.cdata,AppData.MovieSize(1:2));
                AppData.mov(end+1)=f;%getframe(AppData.handles.axes1);
            end
%             while t<0.2 %0.002
%                 t=toc;
%             end
            drawnow
            AppData.RockDone=~get(AppData.handles.Rock,'Value');
            tic;
        end
%     else
%         AppData.RockDone=true;
%         if get(AppData.handles.MakeRockMovie,'value')==0 && AppData.MakingRockMovie
%             % stop making movie
%             AppData.MakingRockMovie=false;
%             % and write it out
%             set(AppData.handles.Movie,'Value',0);
%             moviename=['gyrate-',datestr(now),'.avi'];
%             ind=findstr(moviename,' ');
%             moviename(ind)='T';
%             ind=findstr(moviename,':');
%             moviename(ind)='-';
%             moviename=fullfile(AppData.DataFacetsPath,moviename);
%             disp(sprintf('\nPlease wait as writing %s will take some time',moviename));
%             movie2avi(AppData.mov,moviename,'fps',AppData.framerate,'Quality',90,'compression','None');
%             disp(sprintf('Finished writing %s\n',moviename));
%         end
    
end

function AppData=Grow(AppData)
    if get(AppData.handles.Grow,'Value')
        disp('Grow')
        s=get(AppData.handles.Filename,'String');
        kk=1;
        for i=1:length(s)
            p=str2num(s{i});
            ind=i; 
            [AppData,vert,index,area]=meshmaker_display(AppData,AppData.file_indexes(ind));
            AppData.areas.p_measured(kk)=p;
            AppData.areas.area_measured(kk)=sum(area);
            kk=kk+1;
            drawnow % vert is 120x3
            tempvert=vert'; %reshape(vert,AppData.rows,AppData.cols);
            I.vert(:,i)=tempvert(:);
            I.plast(i)=p;
        end
        AppData.I.vert=fliplr(I.vert);
        AppData.I.plast=fliplr(I.plast);
        h=fspecial('gaussian',7,1);
        hh=h(4,:);
        hh=hh/sum(hh);
        for i=1:size(AppData.I.vert,1)
            xx=AppData.I.vert(i,:);
            smoothed1(i,:)=convn(xx,hh,'same');
        end
        smoothed1(:,end-2:end)=AppData.I.vert(:,end-2:end);
        % drop last point
        smoothed=smoothed1;
        % interpolate
        i=get(AppData.handles.FrameRate,'Value');
        s=get(AppData.handles.FrameRate,'String');
        frstr=s{i};
        framerate=str2num(frstr(1:3));
        increment=0.1;
        if framerate==40
            framerate=20;
            increment=0.2;
        elseif framerate==80
            framerate=20;
            increment=0.4;
        end
        movie_i=1;
        if get(AppData.handles.Movie,'Value')
            [f,AppData]=getAframe(AppData);
            AppData.mov(movie_i)=f;%getframe(AppData.handles.axes1);
            movie_i=movie_i+1;
        end
        start_plast=AppData.I.plast(1);
        flag=true;
        tic;
        ka=1;
        while flag
            flag=false;
            for p=start_plast:increment:AppData.I.plast(end-1)
                vert=interp1(AppData.I.plast,smoothed',p,'linear')';
%                vert=interp1(AppData.I.plast,AppData.I.vert',p)';
                vert=reshape(vert,AppData.rows,AppData.cols);
%                 vert(3,:)=-vert(3,:);
%                 [AppData]=meshmaker_display(AppData,vert);
                [AppData,vert,index,area]=meshmaker_display(AppData,vert);
                AppData.areas.p(ka)=p;
                AppData.areas.area(ka)=sum(area);
                ka=ka+1;
                if ~isfield(AppData,'ph')%==AppData.I.plast(1)
                    AppData.ph=text(0,-8000,0,[num2str(round(p*10)),' hrs'],...
                        'FontSize',20,'Parent',AppData.handles.axes1,'color','w');
                else
                    set(AppData.ph,'String',[num2str(round(p*10)),' hrs'],'color','w');
                end
                t=toc;
                while t<(1/framerate)
                    t=toc;
                end
                tic
                drawnow;
                if get(AppData.handles.Movie,'Value')
%                     f=getframe(AppData.handles.axes1);
                    f=getAframe(AppData);
%                     if isfield(AppData,'rect')
%                         f.cdata=f.cdata(AppData.rect(2):AppData.rect(2)+AppData.rect(4),AppData.rect(1):AppData.rect(1)+AppData.rect(3),:);
%                     end
                    AppData.mov(movie_i)=f;%getframe(AppData.handles.axes1);
                    movie_i=movie_i+1;
                end
                if ~get(AppData.handles.Grow,'Value') % pause during growth
                    start_plast=p;
                    flag=true;
                    set(AppData.FigureName,'pointer','watch');
                    cs=get(AppData.handles.Grow,'Callback');
                    set(AppData.handles.Grow,'Callback','')
                    while ~get(AppData.handles.Grow,'Value')
                        pause(.5)
                        %disp('w')
                    end
                    set(AppData.FigureName,'pointer','arrow');
                    set(AppData.handles.Grow,'Callback',cs)
                    tic
                    break
                end
            end
        end
    end
    drawnow;
    set(AppData.handles.Grow,'Value',0);
    if get(AppData.handles.Movie,'Value')
        tic;
        i=0;
        j=0;
        di=0.07;
        dj=di*0.95;
        Ramp=0.1;
        for k=1:100% AppData.RockDone==false
            t=toc;
            amp=Ramp*get(AppData.handles.RockAmp,'Value');
            si=sin(i)*amp;
            sj=cos(j)*amp;
            camorbit(AppData.handles.axes1,si,sj,'data',[0 0 1]);
            if ~isempty(AppData.handles.axes2)
                camorbit(AppData.handles.axes2,si,sj,'data',[0 0 1]);
            end
            i=i+di;
            j=j+dj;
            drawnow;
            while t<0.01
                t=toc;
            end
            [f,AppData]=getAframe(AppData);
%             f=getframe(AppData.handles.axes1);
            if isfield(AppData,'rect')
                f.cdata=f.cdata(AppData.rect(2):AppData.rect(2)+AppData.rect(4),AppData.rect(1):AppData.rect(1)+AppData.rect(3),:);
            end
            AppData.mov(movie_i)=f;
            movie_i=movie_i+1;
            tic;
        end        
                
        set(AppData.handles.Movie,'Value',0);
        moviename=['bud-',datestr(now),'.avi'];
        ind=findstr(moviename,' ');
        moviename(ind)='T';
        ind=findstr(moviename,':');
        moviename(ind)='-';
        moviename=fullfile(AppData.DataFacetsPath,moviename);
        disp(sprintf('\nPlease wait as writing %s will take some time',moviename));
        movie2avi(AppData.mov,moviename,'fps',AppData.framerate,'Quality',90,'compression','None');
        disp(sprintf('Finished writing %s\n',moviename));
    end
end

function [ff,AppData]=getAframe(AppData)
    f=getframe(AppData.handles.axes1);
    if isfield(AppData,'rect')
        f.cdata=f.cdata(AppData.rect(2):AppData.rect(2)+AppData.rect(4),AppData.rect(1):AppData.rect(1)+AppData.rect(3),:);
    end
    if AppData.isstereo
        f2=getframe(AppData.handles.axes2);
        if isfield(AppData,'rect')
            f2.cdata=f2.cdata(AppData.rect(2):AppData.rect(2)+AppData.rect(4),AppData.rect(1):AppData.rect(1)+AppData.rect(3),:);
        end
        % create a composite cdata3 embedded in a 2048 x 768 image
        [r,c,p]=size(f.cdata);
        aspectratio=c/r;
        if false %aspectratio > 1024/768
            % limit width
            newc=1024;
            newr=newc/aspectratio;
            widthflag=false;
        else
            newr=768;
            newc=floor(newr*aspectratio);
            widthflag=true;
        end
        cdata1rs=imresize(f.cdata,[newr,newc]);
        cdata2rs=imresize(f2.cdata,[newr,newc]);
        cdata3rs=[cdata1rs,cdata2rs];
        [rr,cc,pp]=size(cdata3rs);
        combined=uint8(zeros(768,2048,3));
        cmid1=512;
        left1=cmid1/2;
        cmid2=cmid1+1024;
        left2=left1+1024;
        rmid=768/2;
        combined(1:rr,left1:(left1+newc-1),1:3)=cdata1rs;
        combined(1:rr,left2:(left2+newc-1),1:3)=cdata2rs;
        
%         if ~isfield(AppData,'cdatafig')
            AppData.cdatafig=figure('position',[48         293        1827         785]);
%         elseif ~ishandle(AppData.cdatafig)
%             AppData.cdatafig=figure('position',[48         293        1827         785]);
%         end
        iptsetpref('ImshowBorder','tight');
        iptsetpref('ImshowInitialMagnification','fit');
        imshow(combined)
        shg
        drawnow
        ff=getframe(AppData.cdatafig);
%         pause(1)
        close(AppData.cdatafig)
    else
        ff=f;
    end
end

function AppData=Clip(AppData)
    if get(AppData.handles.Clip,'Value')==get(AppData.handles.Clip,'Max')
        set(AppData.Mesh.cliplines_h,'Visible','on')
    else
        set(AppData.Mesh.cliplines_h,'Visible','off')
    end
end

function vert_out=Align(AppData,vert_in,template)
    % start by getting offset and rotation roughly right
    % by aligning scale, translation and rotation to template
    opts.scaling = 1;
    opts.rotation = 1;
    opts.translation = 1;
    % align the data to template using all vertices
    % first reshape
    V=vert_in';
    T=template';
    disp('Initial procrustes with scaling');
    [vert_temp1, template, XX_aligned, scale_ratio] = pcalib_GPA_jab('data', V(:), 'opts', opts,...
        'dimension', 3,'template',T(:));
    % next centre on the base and align accurately to y axis
    vert_temp2= reshape(vert_temp1, 3, length(vert_temp1)/3);
    if isfield(AppData.Mesh,'base_mark_i')
        vert=vert_temp2';
        vxy=vert(AppData.Mesh.base_mark_i(:),:);
        vert_mean_xy=mean(vxy);
        vert=vert-repmat(vert_mean_xy,size(vert,1),1);
        % decide which way using vertices 7 (bottom) and 89 (top)
        if vert(7,3)>=vert(89,3)
            vert(:,3)=-vert(:,3);
        end
        minvert3=min(vert(:,3));
        vert(:,3)=vert(:,3)-minvert3;
    else
        vert=vert_temp2';
    end
    % centred on middle of base
%    vert_origin=vert(AppData.Mesh.between_i(5),:);
    % perform procrustes on rotation only
%     if ind==1
%         vert=rotate_through_180(3,vert);
%         vert=rotate_through_90(2,vert);
%         vert=rotate_through_90(3,vert);
%     end
%     vert_out=vert;
    if ~get(AppData.handles.NormaliseScale,'Value')
        vert_out=vert/scale_ratio;
    else
        vert_out=vert;
    end

    
    %     opts.scaling = 0;
%     opts.rotation = 1;
%     opts.translation = 0;
%     % zero base to z axis
%     axis_to_zero=3;
%     X1=vxy';
%     X1(axis_to_zero,:)=0; % zero z dimension
%     X2=vxy';
%     X=[X1(:),X2(:)];
%     XX1=vert'; % dummy to replace template
%     XX2=vert';
%     XX=[XX1(:),XX2(:)];
%     disp('align base to z axis');
%     [X_aligned, template, XX_aligned] = pcalib_GPA_jab('data', X, 'opts', opts,...
%         'dimension', 3,'extra_data',XX,'axis_to_zero',axis_to_zero);
% 
%     axis_to_zero=2;
%     vz=vert([AppData.Mesh.between_i(:),AppData.Mesh.centre_ventral_i(:)],:);
%     X1=vz';
%     X1(axis_to_zero,:)=0; % zero z dimension
%     X2=vz';
%     X=[X1(:),X2(:)];
%     XX=XX_aligned; %[XX1(:),XX2(:)];
%     disp('align dorso ventral axis to x');
%     [X_aligned, template, XX_aligned] = pcalib_GPA_jab('data', X, 'opts', opts,...
%         'dimension', 3,'extra_data',XX,'axis_to_zero',axis_to_zero);
% 
%     ZZ2=reshape(XX_aligned(:,2),3,length(XX_aligned(:,1))/3);
% 
%     % now rescale the vertices back to their proper size
%     if ~get(AppData.handles.NormaliseScale,'Value')
%         vert_out=ZZ2'/scale_ratio;
%     else
%         vert_out=ZZ2';
%     end
    %AppData.all_vert{ind}=ZZ2';
end


function plot3axis(plot3axis,Z1,colour)
    plot3([Z1(1,1) Z1(1,3)],[Z1(2,1) Z1(2,3)],[Z1(3,1) Z1(3,3)],['D--',colour],'parent',plot3axis);
    plot3([Z1(1,1) Z1(1,2)],[Z1(2,1) Z1(2,2)],[Z1(3,1) Z1(3,2)],['O-',colour],'parent',plot3axis);
    plot3([Z1(1,1) ],[Z1(2,1) ],[Z1(3,1) ],['*',colour],'parent',plot3axis);
end

function AppData=ChoosePart(AppData,option)
    AppData=ChoosePetal(AppData,option)
end

function AppData=VertexNumbers(AppData,option)
    state=get(AppData.handles.VertexNumbers,'Value')==get(AppData.handles.VertexNumbers,'Max');
    if state
        set(AppData.Mesh.numbers.h,'Visible','on');
    else
        set(AppData.Mesh.numbers.h,'Visible','off');
    end
end

function AppData=ExportObj(AppData)
    %read in the facets file - don't need to do this anymore
    %[fname,path] = uigetfile('*.xls','Facet File',AppData.DataFacetsPath);
    %[num, txt]=xlsread([path,fname],1);
    for i = 1:length(AppData.all_vert)
        %open file for writing
        objFileName = [AppData.pathname,AppData.files{i},'.obj'];
        writeFile = fopen(objFileName,'w');
        %write out all the verticies
        for j = 1:length(AppData.all_vert{i})
            tmpString = num2str(AppData.all_vert{i}(j,1:3));
            fprintf(writeFile,['v ',tmpString,'\n']);
        end
        %write out all the faces 
        for j = 1:length(AppData.fac)
            tmpString = num2str(AppData.fac(j,1:3));
            fprintf(writeFile,['f ',tmpString,'\n']);
        end

%         %face = num+1;
%         for j = 1:length(face)
%             tmpString = num2str(face(j,1:3));
%             fprintf(writeFile,['f ',tmpString,'\n']);
%         end
        
        %close the file
        fclose(writeFile);
    end
end

function AppData=Filename(AppData)
    i=get(AppData.handles.Filename,'Value');
    AppData=meshmaker_display(AppData,AppData.file_indexes(i));
    AppData=ChoosePetal(AppData,0);
end

function HideBud(AppData)
    f=fieldnames(AppData.Mesh);
    for i=1:length(f)
        set(AppData.Mesh.(f{i}).h,'Visible','off');
    end
end

function ShowBud(AppData)
    f=fieldnames(AppData.Mesh);
    for i=1:length(f)
        set(AppData.Mesh.(f{i}).h,'Visible','on');
    end
end

function ShowBudParts(AppData)
    b=fieldnames(AppData.handles);
    for i=1:length(f)
        for j=1:length(b)
            switch b{i}
                case 'V'
                case 'VTube'
                    %                 if get(AppData.handles.(b{i},'Value')==1
                    %                 end
            end
        end
        set(AppData.Mesh.(f{i}).h,'Visible','on');

    end
end

function AppData=ChoosePetal(AppData,option)
    if strcmp(AppData.Mesh.OrganType,'SnapdragonFlower')
        i=get(AppData.handles.Filename,'Value');
        ii=AppData.file_indexes(i);
        b=fieldnames(AppData.handles);
        for i=1:length(b)
            if strcmp(get(AppData.handles.(b{i}),'Type'),'uicontrol')
                names=get(AppData.handles.(b{i}),'String');
                if iscell(names)
                    name=names{1};
                else
                    name=names;%{ii};
                end
                switch upper(name)
                    case {upper('Dorsal Right'),upper('Lateral Right')...
                            ,upper('Ventral Left'),upper('Ventral Right')...
                            ,upper('Dorsal Left'),upper('Lateral Left')}
                        if (get(AppData.handles.(b{i}),'Value'))==(get(AppData.handles.(b{i}),'Max'))
                            set(AppData.Mesh.(b{i}).hn,'Visible','on')
                            %Clip(AppData)
                        else
                            set(AppData.Mesh.(b{i}).hn,'Visible','off')
                            %Clip(AppData)
                        end
                    otherwise
                end
            end
        end
    end
    % HideBud(AppData);
    % ShowBudParts(AppData);
end

function AppData=initialise(index,data_vertex_path,DataFacetsPath,...
        facets_file,sheetname,framerate,isstereo,filename_to_plastochron_filename)
    global AppData
    clear AppData
    AppData.FigureName=openfig('ViewGrowthTool','reuse'); % created using 'guide TestInterface'
    set(AppData.FigureName,'Units','pixels','name','ViewGrowthTool')
    m=get(0,'screensize');
    p=floor(get(AppData.FigureName,'position'));
    set(AppData.FigureName,'Position',[1   m(1,4)-p(4)-70   p(3)   p(4)]);
%     set(AppData.FigureName,'Position',[0.4000   86.6923   47.4000   32.0769])
%    cameratoolbar(AppData.FigureName);
    FigureName=AppData.FigureName;
    handles=guihandles(AppData.FigureName); % access the gui components
    %set(handles.StartStop,'callback','SecretMessage(''StartStop'');');
    set(handles.VR, 'callback',[mfilename, '(''ChoosePetal'',1);'],'Value',1);
    set(handles.VL, 'callback',[mfilename, '(''ChoosePetal'',6);'],'Value',1);
    set(handles.LL,'callback',[mfilename, '(''ChoosePetal'',2);'],'Value',1);
    set(handles.LR,'callback',[mfilename, '(''ChoosePetal'',3);'],'Value',1);
    set(handles.DL,'callback',[mfilename, '(''ChoosePetal'',4);'],'Value',1);
    set(handles.DR,'callback',[mfilename, '(''ChoosePetal'',5);'],'Value',1);

    set(handles.VertexNumbers,'callback',[mfilename, '(''VertexNumbers'',0);'],'Value',0);
    
    set(handles.ExportObj,'callback',[mfilename, '(''ExportObj'',0);'],'Enable','on');

    set(handles.Filename,'callback',[mfilename, '(''Filename'',1);']);
    set(handles.Grow,'callback',[mfilename, '(''Grow'',1);'],'Enable','on');
    set(handles.Rock,'callback',[mfilename, '(''Rock'',1);']);
%     set(handles.MakeRockMovie,'callback',[mfilename, '(''MakeRockMovie'',1);']);
%     set(handles.Movie,'callback',[mfilename, '(''Movie'',1);'],'Enable','on');
    set(handles.Rotate,'callback',[mfilename, '(''Rotate'',1);'],'Enable','on');

    set(handles.RotX,'callback',[mfilename, '(''RotateXYZ'',1);'],'Enable','on');
%     set(handles.RotY,'callback',[mfilename, '(''RotateXYZ'',1);'],'Enable','on');
%     set(handles.RotZ,'callback',[mfilename, '(''RotateXYZ'',2);'],'Enable','off','Visible','off');

    set(handles.Quit,'callback',[mfilename, '(''Quit'',0);'],'Enable','on');
%     set(handles.Align,'callback',[mfilename, '(''Align'',0);'],'Value',0,'Visible','on');
    set(handles.NormaliseScale,'Value',0,'Visible','on');
    set(handles.Zero,'callback',[mfilename, '(''Zero'',0);'],'Enable','on');
    set(handles.SaveAs,'callback',[mfilename, '(''SaveAs'',0);'],'Enable','on');
    set(handles.FlipUpDown,'callback',[mfilename, '(''FlipUpDown'',0);'],'Enable','on');
    set(handles.Oblique,'callback',[mfilename, '(''Zero'',1);'],'Enable','on');
    set(handles.GrabRect,'callback',[mfilename, '(''GrabRect'',0);'],'Enable','on');
    set(handles.ViewAxes,'callback',[mfilename, '(''ViewAxes'',0);'],'Enable','on');
    set(handles.PlotAreaGrowth,'callback',[mfilename, '(''PlotAreaGrowth'',0);'],'Enable','on');
    AppData.FigureCloseRequestFcn=get(AppData.FigureName,'CloseRequestFcn');
    set(AppData.FigureName,'CloseRequestFcn',[mfilename, '(''Quit'',1);']);

    if isstereo
        set(handles.Disparity,'callback',[mfilename, '(''Disparity'',1);'],'Enable','on','Value',0.02,'Max',0.15);
        set(handles.Disparity2,'Enable','on');
    else
        set(handles.Disparity,'Enable','off');
        set(handles.Disparity2,'Enable','off');
    end
    pr=get(handles.RockAmp,'Position');
    pr(4)=1.15385;
    set(handles.RockAmp,'Position',pr,'Value',0.25);
    set(handles.Clip,'callback',[mfilename, '(''Clip'',0);'],'Value',0,'Visible','on');
    AppData.handles=handles; % keep the handles
    gui_position=get(handles.figure1,'Position');
    screensize=get(0,'screensize');
    AppData.LeftPlotFigure=figure('NumberTitle','off','Name','ViewGrowthTool left',...
        'WindowStyle','normal','DockControls','off');
    set(AppData.LeftPlotFigure,'Position',[250,screensize(4)-780,700,700],'color','k')
    AppData.handles.axes1=axes; %AppData.LeftPlotFigure;
    if isstereo
        m=get(0,'MonitorPositions');
        if size(m,1)==1
            AppData.RightPlotFigure=figure('NumberTitle','off','Name','ViewGrowthTool right',...
                'WindowStyle','normal','DockControls','off','color','k');
            set(AppData.RightPlotFigure,'Position',[960,screensize(4)-780,700,700]);
        else
            leftpos=get(AppData.LeftPlotFigure,'Position');
            rightpos=leftpos;
            rightpos(1)=rightpos(1)-1+m(2,1);
            AppData.RightPlotFigure=figure('NumberTitle','off','Name','ViewGrowthTool right',...
                'WindowStyle','normal','DockControls','off','color','k');
            set(AppData.RightPlotFigure,'position',rightpos);
        end
        AppData.handles.axes2=axes;
    else
        AppData.RightPlotFigure=[];
        AppData.handles.axes2=[];
    end
    drawnow
    
    AppData.data_vertex_path=data_vertex_path;
    AppData.DataFacetsPath=DataFacetsPath;
    AppData.framerate=framerate;
    AppData=get_vertices(AppData,data_vertex_path,filename_to_plastochron_filename);
    [AppData.fac,AppData.labels]=get_facets(fullfile(DataFacetsPath,facets_file),sheetname);
    AppData=SpecifyPetalParts(AppData);
    AppData.IndexOfReferenceBud=AppData.file_indexes(ceil(1*size(AppData.all_vert,2)/4));%index;
    AppData=meshmaker_display(AppData,index,AppData.IndexOfReferenceBud);
    AppData=SpecifyPetalParts(AppData);
    AppData=Filename(AppData);
    if strcmp(AppData.Mesh.OrganType,'none')
        set(handles.VR, 'enable','off');
        set(handles.VL, 'enable','off');
        set(handles.LL,'enable','off');
        set(handles.LR,'enable','off');
        set(handles.DL,'enable','off');
        set(handles.DR,'enable','off');
    end
    AppData.MakingRockMovie=false;
%     AppData=Align(AppData,0);
%     AppData=meshmaker_display(AppData);
end

function result=isstereo
    result=false;
end

function value=cost_of_translation(params)
    global AppData
    disp(params)
    vert=AppData.vert-repmat(params,size(AppData.vert,1),1);
    value=sum(sum(AppData.vert.^2));
end

function [fac, labels]=get_facets(filename1,sheet)
    if nargin<1
        filename1='template_triangles_jab.xml';
        sheet='Sheet1';
    end
    [num, txt]=xlsread(filename1,sheet);
    if size(num,1)~=size(txt,1)
        s1=(sprintf('the first two rows of'));
        s2=filename1;
        s3=(sprintf('Should look like \nPetal Side Part 0 0 0'));
        s4=(sprintf('Petal Side Part 0 0 0'));
        s5=('Then something like');
        s6=('D L L 77 5 105');
        s={s1;s2;s3;s4;s5;s6};
        uiwait(warndlg(s));
        error('Error');
    end
    labels=txt(3:end,:);
    fac=num(3:end,1:end)+1;   
end


function AppData=get_vertices(AppData,data_vertex_path,filename_to_plastochron_filename)
    %     if nargin<1
%     Scan_stage=load('filename_to_plastochron_mapping.txt','ascii');
    Scan_stage=load(filename_to_plastochron_filename,'ascii');
    % the template is
%     filename2=template_file;%'910gfp163_x8p91y8p91z8p91pop_3D.csv';
%     vert{1}=load(filename2);
%     stage(1)=27;
%     stage_str=num2str(stage(1));
%     files(1)={filename2};
%     set(AppData.handles.Filename,'String',{stage_str});
    pathname=data_vertex_path; %'PRJ_wild_type\PointModels\910gfp163_x8p91y8p91z8p91pop_3D';
    AppData.pathname=pathname;
    d=dir(fullfile(pathname,'*.mat'));
%paul in here (28/08/08)      
    tmpIdx = 1;
%paul out here (28/08/08)
    for i=1:length(d)
        filename2=fullfile(pathname,d(i).name);
        vert1=load(filename2);
        [rr,cc]=size(vert1.pts);
        cols=rr/3;
        if rem(rr,3)
            error('vertex data is wrong');
        end
        if i==1
            AppData.rows=3;
            AppData.cols=cols;
        end
        if AppData.cols~=cols
            error('number of vertices has changed')
        end
        v=reshape(vert1.pts,[3,cols]);
        vert{tmpIdx}=v'/1000; % imported units are microns, convert to mm
%         vert{i}=v';
        % decide on the plastochron stage number
        %the following is old logic
        %begin old code
%         c=d(i).name(5);
%         if findstr('0123456789',c)
%             num=str2num(d(i).name(1:5));
%         else
%             c=d(i).name(4);
%             if findstr('0123456789',c)
%                 num=str2num(d(i).name(1:4));
%             else
%                 num=str2num(d(i).name(1:3));
%             end
%         end
        %end old code
        %trying to write some new logic
        idxs = isstrprop(d(i).name,'digit');
        idxs = ~idxs;
        idx = min(find(idxs == 1));
        num = str2num(d(i).name(1:idx-1));
        %end new logic
                
        ind=find(Scan_stage(:,1)==num);
        
%Paul in here (28/08/08)       
%         if isempty(ind)
%             disp(sprintf('Error trying to match filename %s',d(i).name));
%             disp(sprintf('in particular the leading number %d',num));
%             disp(sprintf('with an entry in scan_stage file %s',filename_to_plastochron_filename));
%             error('Please make the filenames agree');
%         end

      
        if(~isempty(ind))
%         stage(i+1)=Scan_stage(ind,2);
            stage(tmpIdx)=Scan_stage(ind,2);
            stage_str=num2str(Scan_stage(ind,2));
            s=get(AppData.handles.Filename,'String');
            if strcmp(s,'Filename')
                set(AppData.handles.Filename,'String',{stage_str});
            else
                s=[s;stage_str];
                set(AppData.handles.Filename,'String',s);
            end
%         files(i+1)={d(i).name};
            files(tmpIdx)={d(i).name};
            vert2{tmpIdx}=v'/1000;
            tmpIdx = tmpIdx + 1;
        end
%Paul out here (28/08/08)  
    end
    %     end
    %paul in here 28/08/08
    AppData.all_vert=vert2;
    %paul out 28/08/08
    AppData.all_stages=stage;
    AppData.files=files;
    [f,index]=sort(stage,'descend');
    AppData.file_indexes=index;
    s=get(AppData.handles.Filename,'String');
    ss=cell(size(s,1),1);
    for i=1:length(s)
        ss(i)=s(AppData.file_indexes(i));
    end
    set(AppData.handles.Filename,'String',ss);

    % where i=get(AppData.handles.Filename,'Value');
    % to recover the index to the relevant file,
    % AppData.all_vert(AppData.file_indexes(i))
    % to recover the filename, AppData.files(AppData.file_indexes(i))
end

function [M,area]=PlotPartPetal(fac,vert,colour,plotaxes,M,Alpha)
    if nargin<6
        Alpha=0.5;
    end
    %find those faces with no vertices with negative y
    indx=1:size(fac,1);
    if isempty(M.hn)
        if length(plotaxes)==1 % plot mono
            M.hn=trisurf(fac(indx,:),vert(:,1),vert(:,2),vert(:,3),...
                'FaceColor',colour,'Parent',plotaxes,'EdgeColor',colour/1.3);
        else % plot stereo
            M.hn(1)=trisurf(fac(indx,:),vert(:,1),vert(:,2),vert(:,3),...
                'FaceColor',colour,'Parent',plotaxes(1),'EdgeColor',colour/1.3);
            M.hn(2)=trisurf(fac(indx,:),vert(:,1),vert(:,2),vert(:,3),...
                'FaceColor',colour,'Parent',plotaxes(2),'EdgeColor',colour/1.3);
        end
    else
        if length(plotaxes)==1
            set(M.hn,'Vertices',vert)
        else
            set(M.hn(1),'Vertices',vert)
            set(M.hn(2),'Vertices',vert)
        end
    end
%     set(M.hn(1),'FaceVertexAlphaData',Alpha);
%     set(M.hn(2),'FaceVertexAlphaData',Alpha);
    % calculate area of the whole part by looping through the triangles
    area=0;
    for i=1:length(fac)
        area=area+find_triangle_area(vert(fac(i,:),:));
    end
end

function AppData=SpecifyPetalParts(AppData)
    a=unique(AppData.labels);
    astr=upper(strcat(a{:}));
    AppData.Mesh.OrganType='Unknown';
    if strcmp(astr,'DLRTV') % use this a signature for Amelia's snapdragon
        fprintf(1,'From the pattern of letters in the first three columns of the\n');
        fprintf(1,'This appears to be a Snapdragon bud.\n');
        fprintf(1,'If it is not a Snapdragon bud, please change one of the letters globally to something else.\n')
        AppData.Mesh.OrganType='SnapdragonFlower';
        AppData.Mesh.base_marks=[  12   56    7    57     5    58    35    59    33    60  ... % one side left
            27   51    26   52    21    53    19    54    14    55];  % y is negative
        AppData.Mesh.between_dorsal_marks=[8 9 10 11 12]; % top to bottom
        AppData.Mesh.topright_dorsal_marks=[8,75,74,73,13];
        AppData.Mesh.centre_ventral_all=[27,93 92 91 38 90 89 88 28]; % top to bottom
        AppData.Mesh.centre_ventral=[27,93 92 91 38]; % bottom to top
        AppData.Mesh.Ventral_left=[33,32,31,30,36,29,60,61,62,63];
        AppData.Mesh.Ventral_right=[26,51,27,93,92,91,37 90 89 88 28,64,65,66,22,23,24,25];       
%Paul in here 28/08/08
%     else
%         filename=fullfile(pwd,'Templates','landmarks_on_base.xls');
%         if exist(filename)==2
%             x=xlsread(filename);
%             AppData.Mesh.base_marks=x';
%         else
%             fprintf(1,'In ''Templates'' please create a file called ''landmarks_on_base.xls''\n');
%             fprintf(1,'it should contain a single column, top left, of the landmarks that\n');
%             fprintf(1,'that form the base of the organ. ');
%             fprintf(1,'They are used to orient the organ.\n');
%             error('Error reading file ''landmarks_on_base.xls'' ');
%         end
%Paul out here 28/08/08
    end
end

function [AppData,vert,index,area]=meshmaker_display(AppData,index,IndexOfReferenceBud)
    % function [AppData,vert,index]=meshmaker_display(AppData,index,IndexOfReferenceBud)
    %
    % index, access the vertices in AllData.all_vert
    % IndexOfReferenceBud, normally the largest (first)
    %
    % vert is one row per vertex
    if nargin<2
        ind=get(AppData.handles.Filename,'Value');
        index=AppData.file_indexes(ind);
        vert=AppData.all_vert{index};
        template_index=index; %AppData.IndexOfReferenceBud;
        vert(:,3)=-vert(:,3);
    end
    if nargin<3
        if prod(size(index))==1
            vert=AppData.all_vert{index};
            template_index=AppData.IndexOfReferenceBud;
            template_vert=AppData.all_vert{template_index};
            vert=Align(AppData,vert,template_vert);
        else
            vert=index';
            index=1;
        end
%         vert(:,3)=-vert(:,3);
    else
        vert=AppData.all_vert{index};%{1};%{index};
%         template_vert=AppData.all_vert{IndexOfReferenceBud};
    end
    if get(AppData.handles.FlipUpDown,'Value')~=0
        vert(:,3)=-vert(:,3);        
    end
%     minvert3=min(vert(:,3));
%     vert(:,3)=vert(:,3)-minvert3;
%     AppData=Align(AppData,0)
%     vert(:,1)=-vert(:,1);
    set(AppData.FigureName,'Pointer','watch');
    set(gca,'visible','off');
    labels=AppData.labels;

    fac=AppData.fac;
    if isfield(AppData,'handles')
        plotaxes(1)=AppData.handles.axes1;
        if ~isempty(AppData.handles.axes2)
            plotaxes(2)=AppData.handles.axes2;
        end
    else
        plotaxes=gca;
    end
    %     cla(plotaxes)
    % where i=get(AppData.handles.Filename,'Value');
    % to recover the index to the relevant file,
    % AppData.all_vert(AppData.file_indexes(i))
    % to recover the filename, AppData.files(AppData.file_indexes(i))

    set(AppData.handles.Panel,'Title',sprintf('%s  (Item %d)',AppData.files{1},1));

    for i=1:length(plotaxes)
        hold(plotaxes(i),'on');
    end
%     colour=dorsal_lobe_colour;
    area=[];
    if strcmp(AppData.Mesh.OrganType,'SnapdragonFlower')
        lateral_lobe_colour=[1,0.5,0];
        lateral_tube_colour=[0.7,0.4,0];
        dorsal_lobe_colour=[1,0.2,0.0];
        dorsal_tube_colour=[0.8,0.1,0.0];
        ventral_lobe_colour=[1,1,0.0];
        ventral_tube_colour=[0.7,0.66,0.0];
        petals=[labels{:,1}]';
        sides=[labels{:,2}]';
        parts=[labels{:,3}]';
        if ~isfield(AppData.Mesh,'DLLobe')
            cla(plotaxes)
            AppData.Mesh.DLLobe.hn=[];
            AppData.Mesh.DLTube.hn=[];
            AppData.Mesh.DL.hn=[];
            AppData.Mesh.DRLobe.hn=[];
            AppData.Mesh.DRTube.hn=[];
            AppData.Mesh.DR.hn=[];
            AppData.Mesh.LLLobe.hn=[];
            AppData.Mesh.LLTube.hn=[];
            AppData.Mesh.LL.hn=[];
            AppData.Mesh.LRLobe.hn=[];
            AppData.Mesh.LRTube.hn=[];
            AppData.Mesh.LR.hn=[];
            AppData.Mesh.VLLobe.hn=[];
            AppData.Mesh.VLTube.hn=[];
            AppData.Mesh.VRLobe.hn=[];
            AppData.Mesh.VRTube.hn=[];
            AppData.Mesh.V.hn=[];
            AppData.Mesh.numbers.h=[];
        end
        ii=find((petals=='D' & sides=='L' & parts=='L'));
        [AppData.Mesh.DLLobe,area(end+1)]=PlotPartPetal(fac(ii,:),vert,dorsal_lobe_colour,plotaxes,AppData.Mesh.DLLobe);
        ii=find((petals=='D' & sides=='L' & parts=='T'));
        [AppData.Mesh.DLTube,area(end+1)]=PlotPartPetal(fac(ii,:),vert,dorsal_tube_colour,plotaxes,AppData.Mesh.DLTube);
        AppData.Mesh.DL.hn=[AppData.Mesh.DLLobe.hn,AppData.Mesh.DLTube.hn];

        ii=find((petals=='D' & sides=='R' & parts=='L'));
        [AppData.Mesh.DRLobe,area(end+1)]=PlotPartPetal(fac(ii,:),vert,dorsal_lobe_colour,plotaxes,AppData.Mesh.DRLobe);
        ii=find((petals=='D' & sides=='R' & parts=='T'));
        [AppData.Mesh.DRTube,area(end+1)]=PlotPartPetal(fac(ii,:),vert,dorsal_tube_colour,plotaxes,AppData.Mesh.DRTube);
        AppData.Mesh.DR.hn=[AppData.Mesh.DRLobe.hn,AppData.Mesh.DRTube.hn];

        ii=find((petals=='L' & sides=='L' & parts=='L'));
        [AppData.Mesh.LLLobe,area(end+1)]=PlotPartPetal(fac(ii,:),vert,lateral_lobe_colour,plotaxes,AppData.Mesh.LLLobe);
        ii=find((petals=='L' & sides=='L' & parts=='T'));
        [AppData.Mesh.LLTube,area(end+1)]=PlotPartPetal(fac(ii,:),vert,lateral_tube_colour,plotaxes,AppData.Mesh.LLTube);
        AppData.Mesh.LL.hn=[AppData.Mesh.LLLobe.hn,AppData.Mesh.LLTube.hn];

        ii=find((petals=='L' & sides=='R' & parts=='L'));
        [AppData.Mesh.LRLobe,area(end+1)]=PlotPartPetal(fac(ii,:),vert,lateral_lobe_colour,plotaxes,AppData.Mesh.LRLobe);
        ii=find((petals=='L' & sides=='R' & parts=='T'));
        [AppData.Mesh.LRTube,area(end+1)]=PlotPartPetal(fac(ii,:),vert,lateral_tube_colour,plotaxes,AppData.Mesh.LRTube);
        AppData.Mesh.LR.hn=[AppData.Mesh.LRLobe.hn,AppData.Mesh.LRTube.hn];

        ii=find((petals=='V'  & sides=='L' & parts=='L'));
        [AppData.Mesh.VLLobe,area(end+1)]=PlotPartPetal(fac(ii,:),vert,ventral_lobe_colour,plotaxes,AppData.Mesh.VLLobe);
        ii=find((petals=='V'  & sides=='L' & parts=='T'));
        [AppData.Mesh.VLTube,area(end+1)]=PlotPartPetal(fac(ii,:),vert,ventral_tube_colour,plotaxes,AppData.Mesh.VLTube);
        AppData.Mesh.VL.hn=[AppData.Mesh.VLLobe.hn,AppData.Mesh.VLTube.hn];

        ii=find((petals=='V'  & sides=='R' & parts=='L'));
        [AppData.Mesh.VRLobe,area(end+1)]=PlotPartPetal(fac(ii,:),vert,ventral_lobe_colour,plotaxes,AppData.Mesh.VRLobe);
        ii=find((petals=='V'  & sides=='R' & parts=='T'));
        [AppData.Mesh.VRTube,area(end+1)]=PlotPartPetal(fac(ii,:),vert,ventral_tube_colour,plotaxes,AppData.Mesh.VRTube);
        AppData.Mesh.VR.hn=[AppData.Mesh.VRLobe.hn,AppData.Mesh.VRTube.hn];

        AppData.Mesh.All.hn=[AppData.Mesh.VL.hn,AppData.Mesh.VR.hn,AppData.Mesh.LR.hn,AppData.Mesh.LL.hn,AppData.Mesh.DR.hn,AppData.Mesh.DL.hn];
        AppData.Mesh.All.h=[AppData.Mesh.All.hn];
        graycolour=[0.8,0.8,0.8];
        if strcmp(AppData.Mesh.OrganType,'SnapdragonFlower')
            if isfield(AppData.Mesh,'cliplines_h')
                for i=1:length(plotaxes)
                    set(AppData.Mesh.cliplines_h(i,1),'XData',vert(AppData.Mesh.centre_ventral_all,1),'YData',vert(AppData.Mesh.centre_ventral_all,2),...
                        'ZData',vert(AppData.Mesh.centre_ventral_all,3));
                    set(AppData.Mesh.cliplines_h(i,2),'XData',vert(AppData.Mesh.between_dorsal_marks,1),'YData',vert(AppData.Mesh.between_dorsal_marks,2),...
                        'ZData',vert(AppData.Mesh.between_dorsal_marks,3));
                end
            else
                for i=1:length(plotaxes)
                    hh1=plot3(vert(AppData.Mesh.centre_ventral_all,1),vert(AppData.Mesh.centre_ventral_all,2),...
                        vert(AppData.Mesh.centre_ventral_all,3),'-','Color',graycolour,'Parent',plotaxes(i),'Linewidth',3.5);
                    hh3=plot3(vert(AppData.Mesh.between_dorsal_marks,1),vert(AppData.Mesh.between_dorsal_marks,2),...
                        vert(AppData.Mesh.between_dorsal_marks,3),'-','Color',graycolour,'Parent',plotaxes(i),'Linewidth',2.5);
                    hh=[hh1,hh3];
                    AppData.Mesh.cliplines_h(i,1:length(hh))=hh;
                end
            end
        end
    else
        colour=[0.5,1,0.5];
        if ~isfield(AppData.Mesh,'Organ')
            cla(plotaxes)
            AppData.Mesh.Organ.hn=[];
        end
        [AppData.Mesh.Organ,area(end+1)]=PlotPartPetal(fac,vert,colour,plotaxes,AppData.Mesh.Organ);
    end
    if get(AppData.handles.Black,'value')
        if isfield(AppData.Mesh,'cliplines_h')
            set(AppData.Mesh.cliplines_h,'Color',[1,1,1]);
        end
        set(AppData.LeftPlotFigure,'color','k');
        if ~isempty(AppData.RightPlotFigure)
            set(AppData.RightPlotFigure,'color','k');
            Clip(AppData);
        end
    else
        if isfield(AppData.Mesh,'cliplines_h')
            set(AppData.Mesh.cliplines_h,'Color',[0.4,0.2,0.2]);
            Clip(AppData);
        end
        set(gcf,'color','w');
        set(AppData.LeftPlotFigure,'color','w');
        if ~isempty(AppData.RightPlotFigure)
            set(AppData.RightPlotFigure,'color','w');
        end
    end
    % shading interp
    % lighting phong
    for i=1:length(plotaxes)
        xlabel(plotaxes(i),'x')
        ylabel(plotaxes(i),'y right is positive')
        zlabel(plotaxes(i),'z')
    end
    graycolour=[0.4,0.2,0.2];
    for_C_code=0; % make this 1 if vertex numbers refer to C code
    AppData.MaxLim=25000/1000; % convert microns to mm
    if isfield(AppData.Mesh,'numbers') && ~isempty(AppData.Mesh.numbers.h)
        for i=1:length(vert)
            set(AppData.Mesh.numbers.h(i),'Position',vert(i,:));
        end        
    else
        if strcmp(AppData.Mesh.OrganType,'SnapdragonFlower')
            AppData.Mesh.numbers.h=[];
            AppData.Mesh.base_mark_h=[];
            AppData.Mesh.between_h=[];
            AppData.Mesh.Ventral_left_h=[];
            AppData.Mesh.Ventral_right_h=[];
            AppData.Mesh.centre_ventral_h=[];
            AppData.Mesh.base_mark_i=[];
            AppData.Mesh.between_i=[];
            AppData.Mesh.Ventral_left_i=[];
            AppData.Mesh.Ventral_right_i=[];
            AppData.Mesh.centre_ventral_i=[];
            for_C_code=0; % or 1 if you want to check vertex numbers with those in C code
            for ii=1:length(plotaxes)
                h=[];
                for i=1:length(vert)
                    h(end+1)=text(vert(i,1),vert(i,2),vert(i,3),...
                        num2str(i-for_C_code),'Parent',plotaxes(ii),'color','w');
                    if ~isempty(find(i==AppData.Mesh.base_marks))
                        AppData.Mesh.base_mark_h(end+1)=h(end); %AppData.Mesh.numbers.h(end);
                        AppData.Mesh.base_mark_i(end+1)=i;
                    end
                    if ~isempty(find(i==AppData.Mesh.centre_ventral))
                        AppData.Mesh.centre_ventral_h(end+1)=h(end); %AppData.Mesh.numbers.h(end);
                        AppData.Mesh.centre_ventral_i(end+1)=i;
                    end
                    if ~isempty(find(i==AppData.Mesh.between_dorsal_marks))
                        AppData.Mesh.between_h(end+1)=h(end); %AppData.Mesh.numbers.h(end);
                        AppData.Mesh.between_i(end+1)=i;
                    end
                    if ~isempty(find(i==AppData.Mesh.Ventral_right))
                        AppData.Mesh.Ventral_left_h(end+1)=h(end); %AppData.Mesh.numbers.h(end);
                        AppData.Mesh.Ventral_left_i(end+1)=i;
                    end
                    if ~isempty(find(i==AppData.Mesh.Ventral_left))
                        AppData.Mesh.Ventral_left_h(end+1)=h(end); %AppData.Mesh.numbers.h(end);
                        AppData.Mesh.Ventral_left_i(end+1)=i;
                    end
                end
                AppData.Mesh.numbers.h(ii,1:length(h))=h;
            end
        else
            AppData.Mesh.numbers.h=[];
            for ii=1:length(plotaxes)
                h=[];
                for i=1:length(vert)
                    h(end+1)=text(vert(i,1),vert(i,2),vert(i,3),...
                        num2str(i-for_C_code),'Parent',plotaxes(ii),'color','w');
                    if isfield(AppData.Mesh,'base_marks')
                        if i==1
                            AppData.Mesh.base_mark_h=[];
                            AppData.Mesh.base_mark_i=[];
                        end
                        if ~isempty(find(i==AppData.Mesh.base_marks))
                            AppData.Mesh.base_mark_h(end+1)=h(end); %AppData.Mesh.numbers.h(end);
                            AppData.Mesh.base_mark_i(end+1)=i;
                        end
                    end
                end
                AppData.Mesh.numbers.h(ii,1:length(h))=h;
            end
        end
    end
    for ii=1:length(plotaxes)
        set(AppData.Mesh.numbers.h(ii,:),'Visible','off');
        hold(plotaxes(ii),'off');
        axis(plotaxes(ii),'vis3d')
        axis(plotaxes(ii),'equal')
        AppData.axis=[-AppData.MaxLim/3, AppData.MaxLim/3,...
            -AppData.MaxLim/3, AppData.MaxLim/3,...
            -(+0) (AppData.MaxLim-0)];
        %             AppData.axis=[-AppData.MaxLim/3 AppData.MaxLim/3 -AppData.MaxLim/3 AppData.MaxLim/3 -(+1000) (AppData.MaxLim-0)];
%         set(plotaxes(ii),'XGrid','off','YGrid','off','ZGrid','off','View',[90 0],...
        set(plotaxes(ii),'XGrid','off','YGrid','off','ZGrid','off',...
            'XLim',AppData.axis(1:2),'YLim',AppData.axis(1:2),'ZLim',AppData.axis(5:6))

        if ii==1
            disparity=0; %get(AppData.handles.Disparity,'Value');
        elseif ii==2
            disparity=get(AppData.handles.Disparity,'Value');
        end
        stereoTransfer( plotaxes(ii), ...
            get( plotaxes(1),'CameraPosition'), ...
            get( plotaxes(1),'CameraTarget'), ...
            get( plotaxes(1),'CameraUpVector'), ...
            disparity );
    end
    if get(AppData.handles.ViewAxes,'Value')
        axis(plotaxes,'on');
    else
        axis(plotaxes,'off');
    end
    set(AppData.FigureName,'Pointer','Arrow');
end

function [X, template, XX, scale_ratio] = pcalib_GPA_jab(varargin)
    % function [X] = pcalib_GPA(varargin)
    %
    % Generalized procrustes analysis
    %
    % Arguments presented in pairs,
    % string,      value
    % 'data',X     one column per mesh
    % 'template',  T  one column per vertex and one row per coordinate
    % 'opts', Opt  3 element vector of [0 1] flags
    %                [scaling,translation,rotation]
    % 'dimension',scalar  the number of dimension in the mesh
    % 'extra_data',one column per mesh. Procrustes will be applied to these
    %                using the same parameters as for X
    %
    % X,  one column per mesh of input X aligned to first column
    % template,      the template used (should be the first column)
    % XX, one column per mesh of input XX aligned according to X
    %
    % figure(1); clf;
    % subplot(1,2,1);
    % template = [0 1 1 0; 0 0 1 1];
    % template = template(:);
    % X = repmat(template, 1, 100);
    % X = X + rand(size(X))/10;
    % D = pcalib_randDeformMatrix2D(X);
    % pcalib_plotshapematrix(D, 'g', gca);
    % opts.scaling = 0;
    % opts.rotation = 1;
    % opts.translation = 1;
    % [X, template] = pcalib_GPA('data',D, 'template', template, 'opts', opts);
    % subplot(1,2,2);
    % pcalib_plotshapematrix(X, 'r', gca);
    %
    % Dr. A. I. Hanna (2006) J. Andrew Bangham (2008)
    error(nargchk(0,inf,nargin));
    template = [];
    template_flag = 0;
    X = [];
    opts.scaling = 0;
    opts.translation = 0;
    opts.rotation = 1;
    maxiter = 100;
    dim = 2;
    extra_flag=false;
    axis_to_zero=3;
    axis_to_zero_flag=false;
    if mod(length(varargin),2) ~= 0
        % input args have not com in pairs, woe is me
        error('Arguments to LISTDLG must come param/value in pairs.')
    end
    for i=1:2:length(varargin)
        switch lower(varargin{i})
            case 'data'
                X = varargin{i+1};
            case 'template'
                template = varargin{i+1};
                template_flag = 1;
                maxiter = 0;
            case 'opts'
                opts = varargin{i+1};
            case 'dimension'
                dim = varargin{i+1};
            case 'extra_data'
                XX = varargin{i+1};
                extra_flag=true;
            case 'axis_to_zero'
                axis_to_zero = varargin{i+1};
                axis_to_zero_flag=true;
            otherwise
                error(['Unknown parameter name passed to pcalib_GPA.  Name was ' varargin{i}])
        end
    end
    if isempty(X)
        error('Data parameter is required.')
    end
    if extra_flag
        if size(X,2)~=size(XX,2)
            error('extra_data must have the same number of columns as data')
        end
    end
    if size(template,2)==1
        template = reshape(template, dim, length(template)/dim);
    end
    if size(template,1)>3
        template = template';
    end
    % if the template is empty, then use the first shape in the data set
    if isempty(template)
        template = reshape(X(:,1), dim, length(X(:,1))/dim);
        if axis_to_zero_flag
            template(axis_to_zero,:)=0;
        end
    end
    % centroid is the mean and scale is the standard deviation
    [template, temp_centroid, temp_scale] = pcalib_centershape(template(:), dim);
    scaling = opts.scaling;
    rotation = opts.rotation;
    translation = opts.translation;
    %temp_scale  = pcalib_shapesizemetric(template, dim);
    dt = realmax;
    thr = 10e-10;
    iter = 0;
    while (dt>thr) && (iter <= maxiter)
        for i=1:size(X,2) % loop through the individual meshes
            x = X(:,i); % extract the mesh
            [x, mu] = pcalib_centershape(x, dim); % and convert to one vertex per column
            if extra_flag
                xx=XX(:,i);
                [xx, mu] = pcalib_centershape(xx, dim, mu);
            end
            % opts(1) is the scaling flag
            if scaling==1
                [x,scale_ratio] = pcalib_scaleshape(x, temp_scale);
                if extra_flag
                    xx = xx*scale_ratio;
                end
                %[x] = pcalib_scaleshape(x, 1);
            end
            % opts(2) is the rotation flag. Here R refers to the rotation
            % matrix [cos(theta) -sin(theta); sin(theta) cos(theta)] to align
            % the points.
            if rotation==1
                [x, R] = pcalib_rotateshape(x, template);
                if extra_flag
                    xx = R*xx;
                end
            end
            % opts(3) is the translation flag
            if translation==1
                %[x] = pcalib_transshape(x, zeros(dim,1));
                [x] = pcalib_transshape(x, temp_centroid);
                if extra_flag
                    xx= xx + temp_centroid*ones(1, size(xx, 2));
                end
            else
                [x] = pcalib_transshape(x, mu);
                if extra_flag
                    xx= xx + mu*ones(1, size(xx, 2));
                end
            end
            X(:,i) = x(:);
            if extra_flag
                XX(:,i)=xx(:);
            else
                XX=x;
            end
        end
        if template_flag==0
            if axis_to_zero_flag
                new_template=x;
                new_template(axis_to_zero,:)=0;
            else
                [new_template, temp_centroid] = pcalib_estimate_template(X, dim);
            end
            dt = sum(sqrt(sum((template - new_template).^2,1)));
            disp(sprintf('%5d: dt=%12.10f',iter,dt));
            template = new_template;
        end
        iter = iter + 1;
    end
    [new_template, temp_centroid] = pcalib_estimate_template(X, dim);
    dt = sum(sqrt(sum((template - new_template).^2,1)));
    fprintf('pcalib_GPA: Converged to %5.5f in %d iterations\n', dt, iter);
    template = mean(X,2);
    template = reshape(template, dim, length(template)/dim);
end

function [template, temp_centroid] = pcalib_estimate_template(X, dim)
    % function [template, temp_centroid] = estimate_template(X, dim)
    %
    % Dr. A. I. Hanna (2007);
    template = mean(X,2);
    [template, temp_centroid] = pcalib_centershape(template, dim);
end

function [x] = pcalib_transshape(x, temp_centroid)
    % function [x] = pcalib_transshape(x, temp_centroid)
    %
    % Dr. A. I. Hanna (2007)
    x = x + temp_centroid*ones(1, size(x, 2));
end

function [x, R] = pcalib_rotateshape(x, template)
    % function [x, R] = pcalib_rotateshape(x, template)
    %
    % Dr. A. I. Hanna (2007)
    A = x*template';
    [U, S, V] = svd(A);
    R = V*U';
    x = R*x;
end

function [x,scale_ratio] = pcalib_scaleshape(x, temp_scale)
    % function [xscale_ratio] = pcalib_scaleshape(x, temp_scale)
    %
    % metric here is sqrt{ sum[ (x-xbar)^2 + (y-ybar)^2 ] }
    % since we have already centered the shape, (i.e. xbar = ybar = 0),
    % we can scale using the ratio of standard deviations.
    % i.e. between the template and x
    %
    % Dr. A. I. Hanna (2007)
    Sx = pcalib_shapesizemetric(x, size(x,1));
    scale_ratio = temp_scale/Sx;
    x = x.*scale_ratio;
end

function [S] = pcalib_shapesizemetric(x, dim)
    % function [S] = pcalib_shapesizemetric(x, dim)
    %
    % this is S = sqrt( sum (x(j) - mu_x)^2 + (y(j) - mu_y)^2 )
    %
    % Dr. A. I. Hanna (2007);
    [x, mu] = pcalib_centershape(x(:), dim);
    S = sqrt(sum(x(:).^2));
end

function [x, mu, sd] = pcalib_centershape(x, dim, mu)
    % function [x, mu, sd] = pcalib_centershape(x, dim)
    %
    % subtract the mean from shape in column x
    %
    % x, one column representing vertex coordinates
    % dim, number of dimensions
    % mu, optional by default it finds the mean of the data
    %
    % x, one column representing vertices with mean subtracted
    % mu, the mean
    % sd, the standard deviation
    %
    % Dr. A. I. Hanna (2007)

    x = reshape(x, dim, length(x)/dim);
    if nargin<3
        mu = mean(x, 2); % find mean [x,y,z]'
    end
    x = x - mu*ones(1, size(x, 2));
    sd = sqrt(sum(x(:).^2));
end

function R = rot3dmat(xrot, yrot, zrot)
    %generate the rotation matrix required for xrot, yrot, zrot
    Rx  = [1 0 0; 0 cos(xrot) -sin(xrot); 0 sin(xrot) cos(xrot)];
    Ry  = [cos(yrot) 0 sin(yrot); 0 1 0;-sin(yrot) 0 cos(yrot)];
    Rz  = [cos(zrot) -sin(zrot) 0; sin(zrot) cos(zrot) 0; 0 0 1];
    R = Rx*Ry*Rz;
end

function pcalib_plotshapematrix(X, c, axis_h)
    % function pcalib_plotshapematrix(X, c, axis_h)
    %
    % A simple script to plot the shapes given in the matrix X.
    %
    % Inputs:
    %  X - a matrix where each column is of the form x = [x_1, y_1, x_2, y_2, ..., x_N, y_N]';
    %  c - the color to plot the shape (default: c = [0 0 1])
    %  axis_h - the handle of the axis to plot the shape (default: axis_h = gca)
    %
    % Dr. A. I. Hanna (2006)
    if nargin<2
        c = [0 0 1];
    end
    if nargin<3
        axis_h = gca;
    end
    hold(axis_h, 'on');
    for i=1:size(X,2)
        x = X(:,i);
        pcalib_plotshape(x, c, axis_h);
    end
    hold(axis_h, 'off');
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,'CameraUp' );
    end

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

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

    cameraRight = cross( cameraLook, cameraUp );
end

