function [Data,latest_outputs_from_model]=ReadGFtboxMonitor(projectpath,projectname,template,pointmodelsdirectory,verbose)
    if nargin<1
        projectpath='C:\clusterstuff';
    end
    if nargin<2
        projectname='GPT_ArLeaf_ProxNamesRT_110321A';
    end
    if nargin<3
        template=struct;
    end
    if nargin<4
        pointmodelsdirectory='PointModels\1015';
    end
    if nargin<5
        verbose=true;
    end
    snapshotsname=fullfile(projectpath,projectname,'snapshots');
    D=dir(fullfile(snapshotsname,'monitor-*'));
    [d,ind]=sort({D.date});
    % find the default monitor directory to process
    latest_outputs_from_model=D(ind(end)).name;
    % and set up a list box to select another one if needed - cancel and 
    % the default will be used
    MD=dir(fullfile(projectpath,projectname,'snapshots','monitor*'));
    mondirs=fliplr({MD.name});
    [selection,ok]=listdlg('PromptString','Select a monitor file directory:',...
                'SelectionMode','single',...
                'ListString',mondirs,'listsize',[400,200],'Name','Select directory');
    if ok
        latest_outputs_from_model=mondirs{selection};
    else
        % keep default latest_outputs_from_model
    end

    D2=dir(fullfile(snapshotsname,latest_outputs_from_model,'Data*.mat'));
    filenames={D2.name};
    
    range=ceil(linspace(1,length(filenames),100));
    Data=ReadData(fullfile(snapshotsname,latest_outputs_from_model),filenames(range));
    %Data=ReadData(fullfile(snapshotsname,latest_outputs_from_model),filenames(end-2:end),interval);
    fighandle=figure(2);
    pts=get(template,'pts');
    Data=PlotData(Data,fighandle,pts);
end
function SmoothedData=InterpolateVertices(ModelVerts,Dists,Region)
    % ModelVerts, vertices from GFtbox
    % Dists, distances along the vertices
    % Region, morphogen name designating the region
    SmoothedData.DistsWanted=[]; % distances between vertices
    SmoothedData.VertsWanted=[]; % interpolated vertices 
    SmoothedData.LmksWanted=[];  % associated landmark indexes?
    
end
function Data=PlotData(Data,fighandle,template_pts)
    if nargin<2
        fighandle=figure(2);
    end
    for i=1:length(Data)
        modelled_vertices=zeros(size(template_pts,2),3);
        D=Data(i).MonData;
        Regions=D.Results.REGN;
        XYZ=D.Results.XYZ;
        DISTS=D.Results.D;
        clf(fighandle)
        hold on
        % first process the midlines
        indk=[];
        for k=1:length(Regions)
            if regexp(Regions{k},'MID')
                indk(end+1)=k;
            end
        end
        modelled_pts=zeros(size(template_pts));
        Petiole=[];
        Lamina=[];
        % so now go through the vertices along the midlines
        for j=1:length(indk)
            ind=indk(j);
            region=Regions{ind};
            verts=XYZ{ind};
            dists=cumsum(DISTS{ind}')';
            flat_verts=zeros(size(verts));
            flat_verts(:,2)=dists;
            ThisDISTS=[0;sqrt(sum((diff(verts).^2)')')];
            landmarks=lookupmodel(region); % get details of the landmarks of interest
            lmks=[landmarks.endpoints(1),landmarks.arc,landmarks.endpoints(2)]; % list of landmarks
            dists_wanted=linspace(dists((1)),dists((end)),length(lmks)); % evenly space secondary landmarks
            if size(verts,1)>1
                plot3(verts(:,1),verts(:,2),verts(:,3),'.--b');
                plot3(verts(1,1),verts(1,2),verts(1,3)+0.01,'+g','markersize',10);
                plot3(verts(end,1),verts(end,2),verts(end,3)+0.01,'+g','markersize',10);
                linverts(i,j)=interp1(dists,verts,'linear','pp'); % interpolate a smooth line through middle of thick surface
                Dists_wanted{j}=dists_wanted; % these are the interpolated distances between vertices
                MidlineVerts_wanted{j}= (ppval(linverts(i,j),Dists_wanted{j})); % and interpolated vertex positions around curve
                MidlineLmks_wanted{j}=lmks; % the associated landmark indexes
                % also need interpolated thickness
                
                MidlineThickness_wanted{j}=lmks; % interpolated thicknesses at each point
                if isempty(regexp(region,'PET'))
                    Lamina =j;
                else
                    Petiole =j;
                end
                VertsW=MidlineVerts_wanted{j};
                plot3(VertsW(:,1),VertsW(:,2),VertsW(:,3),'Om');
                text(VertsW(1,1),VertsW(1,2),VertsW(1,3),'1'); %sprintf('%d %d %d',VertsW(1)));
            else
                plot3(verts(:,1),verts(:,2),verts(:,3),'*r');
            end
        end
        % now that we have the two midline arcs we can process the margins
        for j=1:length(Regions)
            
            if isempty(find(j==indk))
                region=Regions{j};
                verts=XYZ{j};
                dists=cumsum(DISTS{j}')';
                ThisDISTS=[0;sqrt(sum((diff(verts).^2)')')];
                landmarks=lookupmodel(region);
                lmks=[landmarks.endpoints(1),landmarks.arc,landmarks.endpoints(2)];
                margin_dists_wanted=linspace(dists((1)),dists((end)),length(lmks));
                if size(verts,1)>1
                    plot3(verts(:,1),verts(:,2),verts(:,3),'.--b');
                    plot3(verts(1,1),verts(1,2),verts(1,3)+0.01,'+g','markersize',10);
                    plot3(verts(end,1),verts(end,2),verts(end,3)+0.01,'+g','markersize',10);
                    linverts(i,j)=interp1(dists,verts,'linear','pp');
                    ind1=regexp(region,'PET');
                    if isempty(ind1) % it must be lamina so compare Verts with MidlineVerts_wanted
                        NumberVerts=length(MidlineVerts_wanted{Lamina}); % from first pass above
                        MidlineLmks=MidlineLmks_wanted{Lamina}; % from first pass above
                        MidLineWanted=MidlineVerts_wanted{Lamina};
                    else
                        NumberVerts=length(MidlineVerts_wanted{Petiole}); % from second pass above
                        MidlineLmks=MidlineLmks_wanted{Petiole}; % from second pass above
                        MidLineWanted=MidlineVerts_wanted{Petiole};
                    end
                    if NumberVerts~=length(margin_dists_wanted)
                        error('number of landmarks in midline not equal to margin')
                    end
                    margin_verts_wanted= (ppval(linverts(i,j),margin_dists_wanted));
                    plot3(margin_verts_wanted(:,1),margin_verts_wanted(:,2),margin_verts_wanted(:,3),'Om');
                    plot3(MidLineWanted(:,1),MidLineWanted(:,2),MidLineWanted(:,3),'py');
                    for kk=1:length(MidLineWanted)
                        plot3([MidLineWanted(kk,1), margin_verts_wanted(kk,1)],[MidLineWanted(kk,2), margin_verts_wanted(kk,2)],...
                            [MidLineWanted(kk,3), margin_verts_wanted(kk,3)],'-g');
                        modelled_vertices(lmks(kk),:)=margin_verts_wanted(kk,:);
                        modelled_vertices(MidlineLmks(kk),:)=MidLineWanted(kk,:);
                    end
                end
            else
                plot3(verts(:,1),verts(:,2),verts(:,3),'*r');
            end
        end
        Data(i).ModelledVertices=modelled_vertices;
        drawnow
        disp(sprintf('Done %d out of %d',i,length(Data)));
    end
end

function landmarks=lookupmodel(region)
    switch upper(region)
        case 'V_RIMLAMINA1'
            landmarks.endpoints=[24, 14];[4, 40];
            landmarks.arc=23:-1:15; %5:13;
        case 'V_RIMPETIOL1'
            landmarks.endpoints=[27 24];%[1, 4];
            landmarks.arc=26:-1:25; %2:3;
        case 'V_RIMLAMINA2'
            landmarks.endpoints=[4, 14];%[17, 40];
            landmarks.arc=5:13; %18:26;
        case 'V_RIMPETIOL2'
            landmarks.endpoints=[1, 4]; %[14, 17];
            landmarks.arc=2:3;
        case 'V_MIDLINE'
            landmarks.endpoints=[31, 14]; %[30, 40];
            landmarks.arc=32:40; %31:39;
        case 'V_PETMIDLINE'
            landmarks.endpoints=[28, 31];%[27, 30];
            landmarks.arc=29:30;%[28,29];
        otherwise
            error('region label mismatch');
    end
end

function Data=ReadData(pathname,filenames)
    for i=1:length(filenames)
        filename=filenames{i};
        stage_t = stageTextToTime( filename(5:end-4));
        dat=load(fullfile(pathname,filename));
        Data(i).MonData=dat.MonData;
        Data(i).stage_t=stage_t;
        Data(i).stage_suffix=filename(5:end-4);
        disp(sprintf('%s t=%f',filename,Data(i).stage_t));
    end
end