
function [] = PointCloudProcessor(arg,option)
% function [] = PointCloudProcessor(arg,option)
%
% Usage:
%            >> PointCloudProcessor('init');
%
% This function displays point clouds in 3D and then fits a polynomial
% function to the cloud. Additionally, you can wrap a 4D surface around the
% cloud which is visualised by level surfaces.
%
%
% Dr. A. I. Hanna (2005)
global AAMShapemodel_figure_handle
global AAMModelGenerator_handle
global AAMModelBuilder_handle
global AAMSetPicker_handle
global UpdateModel_handle
global PMTemplateBuilder_handle
global PDMWalk_handle
global ShapeSpaceTool_handle
global PreProcess_gui_handle
global axesFig_handle
global axesFig_handleA

if nargin<2
    option='';
end
if strcmp(class(arg), 'AAM')
    pdm = get(arg, 'pdm');
    if isempty(pdm)
        uiwait(warndlg('No PDM found','PDM Error'));
        return;
    end
    % Open the GUI for display
    if exist('AAMShapemodel_figure_handle')
        if ishandle(AAMShapemodel_figure_handle)
            set(AAMShapemodel_figure_handle,'visible','off');
        end
    end
    fig = openfig(mfilename, 'reuse');

    %set(fig, 'Pos', [17 2.615 167.2 26.1]);
    Data.handles = guihandles(fig);
    Data.handles.AAM = arg;
    Data.handles.groups = [];
    Data.handles.groupfilename = 'groups.mat';
    Data.option=option;
    Data.handles.figMain = fig;
    ShapeSpaceTool_handle = fig;
    set(Data.handles.modelNametxt, 'String', [getCurrentAAMProjectName, filesep, get(Data.handles.AAM, 'modelDirec')]);
    % Set up the shape axis
    Data.handles.axesFig =  figure;
    Data.handles.shape_axes_handle = gca;
    set(Data.handles.axesFig,'Name','Shape Model', 'Pos', [889   301   775   658])
    axesFig_handle=Data.handles.axesFig;
    set(Data.handles.axesFig, 'closereq', 'closereq');
    %set(gca,'Position',[0.1 0.1 0.8 0.8])
    axis on;

    Data.handles.UIContextMenu=uicontextmenu('Parent',Data.handles.axesFig);
    uimenu(Data.handles.UIContextMenu,'Label','Original image','Callback',[mfilename, '(''ShowOriginalImage'');']);
   % uimenu(Data.handles.UIContextMenu,'Label', 'Use this baseline shape + shape and appearance of the following image (click on target)','Callback', [mfilename, '(''ShowFullModelImage'',0);']);
   % uimenu(Data.handles.UIContextMenu,'Label', 'Use this baseline appearance + shape and appearance of the following image (click on target)','Callback', [mfilename, '(''ShowFullModelImage'',1);']);
   % uimenu(Data.handles.UIContextMenu,'Label', 'Use this shape and the following appearance (click on target)','Callback', [mfilename, '(''ShowFullModelImage'',2);']);

    % Set up the appearance axis
    Data.handles.axesFigA =  figure;
    set(Data.handles.axesFigA, 'closereq', 'closereq');
    Data.handles.app_axes_handle = gca;
    set(Data.handles.axesFigA, 'Visible', 'off');
    set(Data.handles.axesFigA,'Name','Appearance Model', 'Position',[808   314   771   678])
    axesFig_handleA=Data.handles.axesFigA;
    axis on
    Data.handles.UIContextMenuA=uicontextmenu('Parent',Data.handles.axesFigA);
    uimenu(Data.handles.UIContextMenuA,'Label','Original image','Callback',[mfilename, '(''ShowOriginalImageA'');']);
    %uimenu(Data.handles.UIContextMenuA,'Label', 'Use this baseline shape + shape and appearance of the following image (click on target)','Callback', [mfilename, '(''ShowFullModelImageA'',0);']);
    %uimenu(Data.handles.UIContextMenuA,'Label', 'Use this baseline appearance + shape and appearance of the following image (click on target)','Callback', [mfilename, '(''ShowFullModelImageA'',1);']);
    %uimenu(Data.handles.UIContextMenuA,'Label', 'Use this shape and the following appearance (click on target)','Callback', [mfilename, '(''ShowFullModelImageA'',2);']);

    if isinteger(max(get(0,'Children'))) % get(0),'Children' seems to return inconsistent results
        Data.FigNo=max(get(0,'Children'));
    else
        Data.FigNo=3;
    end
    Data.FigNoStart=Data.FigNo;
    Data.PC = [];
    Data.RegCurve = {};
    Data.handles.groups = [];
    Data.SurfaceData.X = [];
    Data.SurfaceData.Y = [];
    Data.SurfaceData.Z = [];
    Data.SurfaceData.V = [];

    Data.Pcomponents=1:3;
    Data.PcomponentsA=1:3;


    Data.FitOpt = 'Data';
    Data.handles.AAM_renderer = 'opengl'

    %Data.handles = LoadGroupInfo(Data.handles); %Data.Pcomponents,option);
    Data.handles.PC = [];
Data.handles.bAll = [];
Data.handles.bVar = [];
Data.handles.bAllA = [];
Data.handles.bVarA = [];
    if isfield(Data.handles, 'bAll')
        %         if isempty(Data.handles.bAll) | (size(Data.handles.bAll{1}, 2)<3)
        %             Data.handles.groups = [];
        %             guidata(fig, Data);
        %             uiwait(msgbox(sprintf('%s\n%s', 'WARNING::There are less than 3 principle components.', 'Use PDM walk to view result.'),'Insufficient Data','modal'));
        %             PointCloudProcessor('Quit');
        %             return;
        %         end
    else
        Data.handles.bAll = [];
        Data.handles.bAllA = [];
    end
    Data.PC = Data.handles.PC;
    Data.RegCurve = cell(1, length(Data.PC));
    Data.SurfaceData = cell(1, length(Data.PC));
    for i=1:length(Data.PC)
        Data.SurfaceData{i}.X = [];
        Data.SurfaceData{i}.Y = [];
        Data.SurfaceData{i}.Z = [];
        Data.SurfaceData{i}.V = [];
    end
    %Data.handles =PlotData(Data.handles, Data.PC, Data.RegCurve, Data.SurfaceData,Data.Pcomponents,Data.option,Data.PcomponentsA);
    guidata(fig, Data);
    set(Data.handles.figMain, 'DoubleBuffer', 'on');%, 'Toolbar', 'figure');
else
    Data = guidata(findobj('Tag', 'ControlGui'));
    %    Data = guidata(gcbo);
    switch arg
%         case 'Use'
%             Data=Use(Data,option);
%             Data=CloseExtras(Data);
%         case 'CloseExtras'
%             Data=CloseExtras(Data);
        %case 'ShowPlot'
        %    Data=ShowPlot(Data);
        case 'LoadGroups'
            Data=LoadGroups(Data);
        case 'SaveGroups'
            Data=SaveGroups(Data);
%         case 'ShowModel'
%             Data=ShowModel(Data);
        case 'RandomColorGroups'
            Data=RandomColorGroups(Data);
            %Data.handles = LoadGroupInfo(Data.handles);

            Data.handles = PlotData(Data.handles, Data.PC, Data.RegCurve, Data.SurfaceData,Data.Pcomponents,Data.option, Data.Pcomponents);
        case 'ResetAxis'
            ResetAxis(Data.handles);
        case 'LDAUsePC'
            Data.handles=LDAUsePC(Data.handles);
        case 'CloseAppWindow'
            set(Data.handles.hide_app_window_chk, 'Value', 1);
            Data.handles = HideAppWindow(Data.handles);

        case 'CloseShapeWindow'
            set(Data.handles.hide_shape_window_chk, 'Value', 1);
            Data.handles = HideShapeWindow(Data.handles);
        case 'ShowOriginalImage'
            Data=ShowOriginalImage(Data);
        case 'ShowPCVariance'
            Data=ShowPCVariance(Data);
        case 'ShowFullModelImage'
            Data=ShowFullModelImage(Data);
        case 'DefineBaseAndTargetGrps'
            Data=DefineBaseAndTargetGrps(Data);
        case 'ShiftToCommonMean'
            Data=ShiftToCommonMean(Data);
            Data.handles = PlotData(Data.handles, Data.PC, Data.RegCurve, Data.SurfaceData, Data.Pcomponents,Data.option, Data.PcomponentsA);
        case 'PlotData'
            Data.handles = PlotData(Data.handles, Data.PC, Data.RegCurve, Data.SurfaceData, Data.Pcomponents,Data.option, Data.PcomponentsA);
        case 'Pcomponents'
            Data = Pcomponents(Data,option);
            Data.handles = PlotData(Data.handles, Data.PC, Data.RegCurve, Data.SurfaceData, Data.Pcomponents,Data.option, Data.PcomponentsA);
        case 'PrintEPS'
            PrintEPS(Data.handles);

        case 'ImportGroup'
            Data.handles = ImportGroup(Data.handles);
            %Data.handles = LoadGroupInfo(Data.handles); %Data.Pcomponents,option);
        case 'ProjectData'
            Data.handles = ProjectData(Data.handles);
        case 'AddGroup'
            Data.handles = AddGroup(Data.handles);
        case 'DeleteGroup'
            Data.handles = DeleteGroup(Data.handles);
            %Data.handles = LoadGroupInfo(Data.handles); %Data.Pcomponents,option);
        case 'UpdateGroupInfo'
            Data.handles = UpdateGroupInfo(Data.handles);
            %Data.handles = PlotData(Data.handles, Data.PC, Data.RegCurve, Data.SurfaceData, Data.Pcomponents,Data.option, Data.PcomponentsA);
            %Data.handles = LoadGroupInfo(Data.handles); %Data.Pcomponents,option);
%         case 'ShowStats'
%             Data = ShowStats(Data);
%             set(Data.handles.CloseExtras,'Enable','on');
        case 'ClearGroups'
            Data.handles.groups = [];
            Data.handles = UpdateGroupInfo(Data.handles);

        case 'ShowMeans'
            Data = ShowMeans(Data);
        case 'Rotate'
            Data = Rotate(Data);
        case 'Zoom'
            Data = Zoom(Data);
        case 'Pan'
            Data = Pan(Data);
        case 'Select'
            Data = Select(Data,option);
            Data=CloseExtras(Data);
        case 'HideShapeWindow'
            Data.handles = HideShapeWindow(Data.handles);
        case 'HideAppWindow'
            Data.handles = HideAppWindow(Data.handles);
        case 'FitData'
            [Data.RegCurve, Data.SurfaceData] = FitData(Data.handles, Data.handles.bAll, Data.Pcomponents, Data.FitOpt);
            Data.handles = PlotData(Data.handles, Data.PC, Data.RegCurve, Data.SurfaceData,Data.Pcomponents,Data.option, Data.Pcomponents);
        case 'ChangeColor'
            [Data.handles] = ChangeColor(Data.handles);
            %Data.handles = LoadGroupInfo(Data.handles);
            Data.handles = PlotData(Data.handles, Data.PC, Data.RegCurve, Data.SurfaceData,Data.Pcomponents,Data.option, Data.Pcomponents);
        case 'ShowNumbers'
            Data.handles = PlotData(Data.handles, Data.PC, Data.RegCurve, Data.SurfaceData,Data.Pcomponents,Data.option, Data.Pcomponents);
        case 'UpdateGroups'
            Data.handles = UpdateGroups(Data.handles);
           % Data.handles = PlotData(Data.handles, Data.PC, Data.RegCurve, Data.SurfaceData, Data.Pcomponents,Data.option, Data.PcomponentsA);
        case 'ChangeIsoLevel'
            Data.handles = ChangeIsoLevel(Data.handles, Data.SurfaceData);
        case 'ToggleEllipsoid'
            Data.handles = ToggleEllipsoid(Data.handles);
        case 'OrderGroups'
            Data.handles = OrderGroups(Data.handles);
        case 'CalcProjVar'
            Data.handles = CalcProjVar(Data.handles);
        case 'Quit'
%            Data=CloseExtras(Data);
            groups = Data.handles.groups;
            %if ~isempty(groups)
            % end
            if exist('AAMShapemodel_figure_handle')
                if ishandle(AAMShapemodel_figure_handle)
                    set(AAMShapemodel_figure_handle,'visible','on');
                end
            end
            [Data.handles] = Quit(Data.handles);
        otherwise
            error('Unknown Option in PointCloudProcessor');
    end
    if isfield(Data.handles, 'figMain')
        guidata(Data.handles.figMain, Data);
    end
end



function handles = HideShapeWindow(handles)
if get(handles.hide_shape_window_chk, 'value') == 1
    set(handles.axesFig, 'Visible', 'off');

else
    set(handles.axesFig, 'Visible', 'on');

end
return
function handles = HideAppWindow(handles)
if get(handles.hide_app_window_chk, 'value') == 1
    set(handles.axesFigA, 'Visible', 'off');

else
    set(handles.axesFigA, 'Visible', 'on');

end
return



% function  Data=Use(Data,option)
% if option<2 % set appearance options
%     if option==0
%         set(Data.handles.UsePC0,'Value',1);
%         set(Data.handles.UsePC1,'Value',0);
%     else
%         set(Data.handles.UsePC0,'Value',0);
%         set(Data.handles.UsePC1,'Value',1);
%     end
% else
%     if option==2
%         set(Data.handles.UsePC2,'Value',1);
%         set(Data.handles.UsePC3,'Value',0);
%     else
%         set(Data.handles.UsePC2,'Value',0);
%         set(Data.handles.UsePC3,'Value',1);
%     end
% end

% function Data=CloseExtras(Data)
% if isempty(Data)
%     return;
% end
% for i=Data.FigNoStart+1:Data.FigNo
%     if ishandle(i)
%         close(i);
%     end
% end
% Data.FigNo=Data.FigNoStart;
% if isfield(Data.handles,'orig_model') % use as proxy for all three
%     if ishandle(Data.handles.orig_model)
%         delete(Data.handles.orig_model);
%         delete(Data.handles.orig_image);
%         delete(Data.handles.orig_plot);
%         rmfield(Data.handles,'orig_model','orig_image','orig_plot');
%     end
% end
% set(Data.handles.CloseExtras,'Enable','off');
%%%%%%
%
%
%%%%%%
function  Data=SaveGroups(Data)
modelDir = get(Data.handles.AAM, 'modelDirec');
[filename, pathname] = uiputfile([pwd, filesep, '*.mat'], 'Pick a group file.');
if isequal(filename,0) || isequal(pathname,0)
    disp('User pressed cancel')
else
    groups = Data.handles.groups;
    save([pathname, filesep, filename], 'groups');
    disp(['User selected ', fullfile(pathname, filename)])
end
return;
%%%%%%
%
%
%%%%%%
function  Data=LoadGroups(Data)
modelDir = get(Data.handles.AAM, 'modelDirec');
[filename, pathname] = uigetfile([pwd, filesep, '*.mat'], 'Pick a group file.');
if isequal(filename,0) || isequal(pathname,0)
    disp('User pressed cancel')
else
    groups = load([pathname, filesep, filename]);
    groups = groups.groups;
    Data.handles.groups = groups;
    Data.handles.groupfilename = [pathname, filesep, filename];
    disp(['User selected ', fullfile(pathname, filename)])
    [Data.handles] = UpdateGroupInfo(Data.handles);
end
return
%%%%%%
%
%
%%%%%%
% function  Data=ShowPlot(Data)
% if isfield(Data.handles,'orig_plot')
%     if get(Data.handles.ShowPlot,'Value')
%         set(Data.handles.orig_plot,'Visible','on');
%     else
%         set(Data.handles.orig_plot,'Visible','off');
%     end
% end
%%%%%
%
%
%%%%%
function  handles=LDAUsePC(handles)
direc = [pwd, filesep, get(handles.AAM, 'modelDirec')];
indices = get(handles.AAM, 'pointindices');

template = get(handles.AAM, 'PointModelTemplate');
template_name = get(template, 'name');
[templatepath, template_name, ext, vers] = fileparts(template_name);
pmdirectory = ['PointModels', filesep, template_name];

if ~exist(pmdirectory, 'dir')
    fprintf('There are no point models for this template');
    return;
end
pdmfile = [direc, filesep, 'mod_pdm.mat'];
if ~exist(pdmfile)
    return;
end
pdm = load(pdmfile);
pdm = pdm.pdm;
[details, detailstr] = getDetails(direc);
projectLDAviaPCA(handles.groups, pmdirectory, pdm, indices, details);
return;
%%%%%
%
%
%%%%%
function  Data=RandomColorGroups(Data)
if isempty(Data.handles.groups)
   uiwait(msgbox('Please load some groups','No Group Info.','modal'));
    return;
end
groups = Data.handles.groups;
for i=1:length(groups)
    color = rand(1,3);
    groups(i).color = color;
    Data.handles.PC{i}.color = color;
end
%setgroups(Data.handles.AAM, groups);
Data.handles.groups = groups;
Data.handles = UpdateGroupInfo(Data.handles);
return
%%%%%
%
%
%%%%%
% function  Data=ShowModel(Data)
% if isfield(Data.handles,'orig_model')
%     if get(Data.handles.ShowModel,'Value')
%         if ishandle(Data.handles.orig_model)
%             set(Data.handles.orig_model,'Visible','on');
%             set(Data.handles.orig_plot,'Color','r');
%         end
%     else
%         if ishandle(Data.handles.orig_model)
%             set(Data.handles.orig_model,'Visible','off');
%             set(Data.handles.orig_plot,'Color','k');
%         end
%     end
% end
%%%%%
%
%
%%%%%
function Data=Plot3D2(Data)
if get(Data.handles.Plot3D2,'Value')
    set(Data.handles.Plot3D2,'Value',1);
    set(Data.handles.Plot3D2,'String','Click to Plot 2D');
else
    set(Data.handles.Plot3D2,'Value',0);
    set(Data.handles.Plot3D2,'String','Click to Plot 3D');
end
%%%%%
%
%
%%%%%
function [Data,positions,pts1,fignos,model_img,sourcestr,img]=ShowImageData(Data,ExternalXm,positions,fignos,picture_img,sourcestr2,Externalimg,flag)
dat=get(gco,'UserData');
mp=get(0,'MonitorPositions');
wid=mp(1,3);
ht=mp(1,4);
dht=ht/8;
figure(Data.handles.axesFig)
set(gcf,'Pointer','watch');
sourcestr=sprintf('%d/%d %s',dat.region,dat.num,dat.filename);
disp(sprintf('Working on image %s',sourcestr));
if nargin<2
    %Data.FigNo=Data.FigNo+1;
    fignos(1:5)=Data.FigNo+1:Data.FigNo+5;
    Num=Data.FigNo-Data.FigNoStart;
    positions{1}=[(wid/20)+((Num-1)/2)*wid/10,ht-1*dht*1.5,wid/10,dht]; % original image
    positions{2}=[(wid/20)+((Num-1)/2)*wid/10,ht-2*dht*1.5,wid/10,dht]; % shape and appearance representation of original
    positions{3}=[(wid/20)+((Num-1)/2)*wid/10,ht-3*dht*1.5,wid/10,dht]; % source baseline shape + shape and appearance of target
    positions{4}=[(wid/20)+((Num-1)/2)*wid/10,ht-4*dht*1.5,wid/10,dht]; % source baseline appearance + shape and appearance of target
    positions{5}=[(wid/20)+((Num-1)/2)*wid/10,ht-5*dht*1.5,wid/10,dht]; % source shape + appearance of target
    sourcestr2='';
else
    sourcestr2=sprintf('(applied to %s)',sourcestr2);
    sourcestr=[sourcestr,sourcestr2];
    disp(sprintf('Working on images %s',sourcestr2));
end
if nargout<2
    figure(fignos(1));
    Data.FigNo=Data.FigNo+1;
    set(fignos(1),'Name',sourcestr,'Position', positions{1});
    imfilename=dat.filename(1:findstr(dat.filename,'_pm.mat')-1);
    file=fullfile('Cropped',[imfilename,'.jpg']);
    if exist(file)==2
        picture_img=imread(file);
    else
        picture_img=[];
    end
    if ~isfield(Data.handles,'orig_model')
        Data.handles.orig_model=[];
    else
        if ~ishandle(Data.handles.orig_model)
            Data.handles.orig_model=[];
        end
    end
    if ~isempty(picture_img)
        Data.handles.orig_model(end+1)=imshow(uint8(picture_img));
        iptsetpref('ImshowBorder','tight');
        hold on
        text(5,0.08*size(picture_img,1),sourcestr,'color','r','Interpreter','none');
        hold off
    end
    if ~isfield(Data.handles,'orig_plot')
        Data.handles.orig_plot=[];
    else
        if ~ishandle(Data.handles.orig_plot)
            Data.handles.orig_plot=[];
        end
    end
    if dat.num>0
        pointindices=get(Data.handles.AAM,'pointindices');
        file=fullfile('PointModels',dat.templatename,dat.filename);
        temp=load(file);
        hold on
        xp=temp.pts(1:2:end);
        yp=temp.pts(2:2:end);
        Data.handles.orig_plot(end+1)=plot(xp(pointindices),yp(pointindices),'.b');
        hold off
    end
%    set(Data.handles.CloseExtras,'Enable','on');
end
% AAM=Data.handles.AAM;
% md=get(AAM,'modelDirec');
% pmt=get(AAM,'PointModelTemplate');
% loops=get(pmt,'loops');
% edges=loops{1};
% temp=load(fullfile(md,'mod_sfam.mat'),'-mat');
% aam=temp.sfam;
% pdm=aam.pdm;
% if dat.num<0 % overall mean
%     img = double(aam.Am); %+aam.P * bWA;
%     pts1=pdm.Xm;
%     if nargin>2
%         pts2=ExternalXm; %+ pdm.P * bW;
%         img2=Externalimg(:); %+ aam.P * bWA;
%     else
%         img2=img;
%     end
% else
%     %if get(Data.handles.UsePC3,'Value')==1 %& nargin>2
%     %    bW=zeros(size(dat.bAll(:)));
%     %    bW(dat.Pcomponents)=dat.bAll(dat.Pcomponents);
%     %else
%         bW=dat.bAll(:);
%     %end
%     pts1=pdm.Xm+ pdm.P * bW;
%     if get(Data.handles.UsePC0,'Value')==1 %& nargin>2
%         bWA=dat.bAllA(:);
%     else
%         bWA=zeros(size(dat.bAllA(:)));
%         bWA(dat.PcomponentsA)=dat.bAllA(dat.PcomponentsA);
%     end
%     img = double(aam.Am)+ aam.P * bWA;
%     if nargin>2
%         pts2=ExternalXm+ pdm.P * bW;
%         img2=Externalimg(:)+ aam.P * bWA;
%     else
%         img2=img;
%     end
% end
% img = reshape(img, [aam.siz(2), aam.siz(1), 3]);
% img2 = reshape(img2, [aam.siz(2), aam.siz(1), 3]);
% pts(:,1)=pts1(1:2:end);
% pts(:,2)=pts1(2:2:end);
% ptsz(:,1)=pts(:,1)-min(pts(:,1)); % offset to zero on axis
% ptsz(:,2)=pts(:,2)-min(pts(:,2));
% model_img=triwarp(img,aam.ompts,ptsz); % current shape and appearance models
% if nargin>2
%     img = reshape(img, [aam.siz(2), aam.siz(1), 3]);
%     pts(:,1)=pts2(1:2:end);
%     pts(:,2)=pts2(2:2:end);
%     ptsz2(:,1)=pts(:,1)-min(pts(:,1)); % offset to zero on axis
%     ptsz2(:,2)=pts(:,2)-min(pts(:,2));
%     model_img2=triwarp(img,aam.ompts,ptsz2); % new baseline + current shape and appearance models
%     model_img3=triwarp(img2,aam.ompts,ptsz); % new baseline appearance + current shape and appearance models
%     pts3=ExternalXm;
%     pts(:,1)=pts3(1:2:end);
%     pts(:,2)=pts3(2:2:end);
%     ptsz3(:,1)=pts(:,1)-min(pts(:,1)); % offset to zero on axis
%     ptsz3(:,2)=pts(:,2)-min(pts(:,2));
%     model_img4=triwarp(img,aam.ompts,ptsz3); % new shape + current appearance models
% end
% if nargout<2
%     figure(fignos(2));
%     Data.FigNo=Data.FigNo+1;
%     set(fignos(2),'Name',sourcestr,'Position', positions{2});
%     clf
%     if ~isempty(model_img)
%         Data.handles.orig_model(end+1)=imshow(uint8(model_img));
%     end
%     iptsetpref('ImshowBorder','tight');
%     hold on
%     text(5,0.08*size(model_img,1),sourcestr,'color','r','Interpreter','none');
%     inds=get(AAM,'pointindices'); % only use points within current set
%     for i=1:length(edges)
%         if ~isempty(intersect(inds,edges(i,:)))
%             ii=find(inds==edges(i,1));
%             jj=find(inds==edges(i,2));
%             Data.handles.orig_plot(end+1)=plot([ptsz(ii,1),ptsz(jj,1)],[ptsz(ii,2),ptsz(jj,2)],'-r');
%         end
%     end
%     hold off
%     if nargin>2
%         if flag==0
%             figure(fignos(3));
%             Data.FigNo=Data.FigNo+1;
%             set(fignos(3),'Name',sourcestr,'Position', positions{3});
%             clf
%             if ~isempty(model_img2)
%                 Data.handles.orig_model(end+1)=imshow(uint8(model_img2));
%             end
%             iptsetpref('ImshowBorder','tight');
%             hold on
%             text(5,0.08*size(model_img2,1),sourcestr,'color','r','Interpreter','none','Interpreter','none');
%             inds=get(AAM,'pointindices'); % only use points within current set
%             for i=1:length(edges)
%                 if ~isempty(intersect(inds,edges(i,:)))
%                     ii=find(inds==edges(i,1));
%                     jj=find(inds==edges(i,2));
%                     Data.handles.orig_plot(end+1)=plot([ptsz2(ii,1),ptsz2(jj,1)],[ptsz2(ii,2),ptsz2(jj,2)],'-r');
%                 end
%             end
%             hold off
%         end
%         if flag==1
%             figure(fignos(4));
%             Data.FigNo=Data.FigNo+1;
%             set(fignos(4),'Name',sourcestr,'Position', positions{4});
%             clf
%             if ~isempty(model_img3)
%                 Data.handles.orig_model(end+1)=imshow(uint8(model_img3));
%             end
%             iptsetpref('ImshowBorder','tight');
%             hold on
%             text(5,0.08*size(model_img3,1),sourcestr,'color','r','Interpreter','none','Interpreter','none');
%             hold on
%             for i=1:length(edges)
%                 if ~isempty(intersect(inds,edges(i,:)))
%                     ii=find(inds==edges(i,1));
%                     jj=find(inds==edges(i,2));
%                     Data.handles.orig_plot(end+1)=plot([ptsz(ii,1),ptsz(jj,1)],[ptsz(ii,2),ptsz(jj,2)],'-r');
%                 end
%             end
%             hold off
%         end
%         if flag==2
%             figure(fignos(5));
%             Data.FigNo=Data.FigNo+1;
%             set(fignos(5),'Name',sourcestr,'Position', positions{5});
%             clf
%             if ~isempty(model_img4)
%                 Data.handles.orig_model(end+1)=imshow(uint8(model_img4));
%             end
%             iptsetpref('ImshowBorder','tight');
%             hold on
%             text(5,0.08*size(model_img3,1),sourcestr,'color','r','Interpreter','none','Interpreter','none');
%             hold on
%             for i=1:length(edges)
%                 if ~isempty(intersect(inds,edges(i,:)))
%                     ii=find(inds==edges(i,1));
%                     jj=find(inds==edges(i,2));
%                     Data.handles.orig_plot(end+1)=plot([ptsz3(ii,1),ptsz3(jj,1)],[ptsz3(ii,2),ptsz3(jj,2)],'-r');
%                 end
%             end
%             hold off
%         end
%     end
% end
figure(Data.handles.axesFig)
set(gcf,'Pointer','arrow');

%%%%%%
%
%
%%%%%%
function Data=DefineBaseAndTargetGrps(Data)
fprintf('Select groups to define reference mean.\n');
bAll = Data.handles.bAll;
grps = Data.handles.groups;
str = {grps.name};
[s,v] = listdlg('Name', 'GroupSelector', 'PromptString','Select Reference Groups:', 'SelectionMode','multiple','ListString',str, 'ListSize', [300 160]);
if isempty(s)
    return;
end
groups.base = s;
fprintf('Select groups to define target mean.\n');
[s,v] = listdlg('Name', 'GroupSelector', 'PromptString','Select Target Groups:', 'SelectionMode','multiple','ListString',str, 'ListSize', [300 160]);
if isempty(s)
    return;
end
groups.target = s;
fprintf('Select groups associated with target mean.\n');
[s,v] = listdlg('Name', 'GroupSelector', 'PromptString','Select Associated Target Groups:', 'SelectionMode','multiple','ListString',str, 'ListSize', [300 160]);
if isempty(s)
    return;
end
groups.associated = s;
shiftgrpsfilename = getShiftedGroupsFilename(Data.handles.AAM);
save(shiftgrpsfilename, 'groups');
return;
%%%%
%
%
%%%%
function shiftgrpsfilename = getShiftedGroupsFilename(aam)
modelDir = [pwd, filesep, get(aam, 'modelDirec')];
shiftgrpsfilename = [modelDir, filesep, 'shifted_grps.mat'];
%%%%%%
%
%
%%%%%%
function Data=ShiftToCommonMean(Data)
if get(Data.handles.shifted_chk, 'Value')==0
    Data.handles = UpdateGroups(Data.handles);
    return;
end
shiftgrpsfilename = getShiftedGroupsFilename(Data.handles.AAM);
if ~exist(shiftgrpsfilename)
    return;
end
shiftgrps = load(shiftgrpsfilename);
shiftgrps = shiftgrps.groups;
bAll = Data.handles.bAll;
shiftmean = calcGroupMeans(bAll, shiftgrps.target);
basemean = calcGroupMeans(bAll, shiftgrps.base);
diff_vec = shiftmean - basemean;
bAll = adjustGroups(bAll, diff_vec, shiftgrps.associated);
Data.handles.bAll = bAll;
return;

%%%%%%
%
%
%%%%%%
function bAll = adjustGroups(bAll, diff_vec, shiftmean)
for i=shiftmean
    bAll{i} = bAll{i} - ones(size(bAll{i},1),1)*diff_vec;
end
return;

%%%%%%
%
%
%%%%%%
function mu = calcGroupMeans(bAll, grps)
X = [];
for i=grps
    X = [X; bAll{i}];
end
mu = mean(X,1);
return;
%%%%%%
%
%
%%%%%%
function Data=ShowPCVariance(Data)
direc = [pwd, filesep, get(Data.handles.AAM, 'modelDirec')];
pdmfile = [direc, filesep, 'mod_pdm.mat'];
pdm = load(pdmfile);
pdm = pdm.pdm;
b = pdm.pca.b;
b = 100*(b./sum(b));
x = 1:length(b);
figure;
stem(x, b, 'MarkerFaceCOlor', 'y');
xlabel('Principle Component');
ylabel('Percentage of Variance Covered (%)');
axis tight;
grid on;
return;
%%%%%%
%
%
%%%%%%
function Data=ShowOriginalImage(Data)
Option=0;
Data=ShowImageData(Data);
%%%%%%
%
%
%%%%%%
function Data=ShowFullModelImage(Data,option)
[Data,positions,model,fignos,model_img,sourcestr,img]=ShowImageData(Data);
figure(Data.handles.axesFig)
obj1=gco;
if option==0
    disp('Use this baseline shape + shape and appearance of the following image (click on target)');
    s1='B';
    s2='SA';
elseif option==1
    disp('Use this baseline appearance + shape and appearance of the following image');
    s1='A';
    s2='SA';
else
    disp('Use this shape and the following appearance (click on target)');
    s1='S';
    s2='A';
end
hold on
if strcmp(get(obj1,'type'),'text')
    pos=get(obj1,'Position');
    if get(Data.handles.Plot3D2,'Value') % 3D
        x1=pos(1);
        y1=pos(2);
        z1=pos(3);
        ph=text(x1+5,y1+5,z1+5,s1, 'Interpreter', 'None');
    else
        x1=pos(1);
        y1=pos(2);
        ph=text(x1+5,y1+5,s1, 'Interpreter', 'None');
    end
else
    if get(Data.handles.Plot3D2,'Value') % 3D
        x1=get(obj1,'Xdata');
        y1=get(obj1,'Ydata');
        z1=get(obj1,'Zdata');
        ph=text(x1,y1,z1,s1, 'Interpreter', 'None');
    else
        x1=get(obj1,'Xdata');
        y1=get(obj1,'Ydata');
        ph=text(x1,y1,s1, 'Interpreter', 'None');
    end
end
dat=get(obj1,'UserData');
set(ph,'UserData',dat,'UIContextMenu',Data.handles.UIContextMenu);
disp('Wait for cross-hairs then right-click on the target point');
[a,b]=ginput(1);
obj2=gco;
disp('Head of ''arrow'' is black')
if strcmp(get(obj1,'type'),'text')
    pos=get(obj1,'Position');
    if get(Data.handles.Plot3D2,'Value') % 3D
        x2=pos(1);
        y2=pos(2);
        z2=pos(3);
        ph=text(x2,y2,z2,s2, 'Interpreter', 'None');
        plot3([x1,x2],[y1,y2],[z1,z2],'-r');
        x3=x2-((x2-x1)/7);
        y3=y2-((y2-y1)/7);
        z3=z2-((z2-z1)/7); % head of 'arrow' is black
        plot3([x3,x2],[y3,y2],[z3,z2],'-k');
        %plot3(x3,y3,z3,'Ok')
    else % 2D
        x2=pos(1);
        y2=pos(2);
        ph=text(x1+5,y1+5,s1, 'Interpreter', 'None');
        plot([x1,x2],[y1,y2],'-r');
        x3=x2-((x2-x1)/7);
        y3=y2-((y2-y1)/7); % head of 'arrow' is black
        plot([x3,x2],[y3,y2],'-k');
    end
else
    if get(Data.handles.Plot3D2,'Value') % 3D
        x2=get(obj2,'Xdata');
        y2=get(obj2,'Ydata');
        z2=get(obj2,'Zdata');
        ph=text(x2,y2,z2,s2, 'Interpreter', 'None');
        plot3([x1,x2],[y1,y2],[z1,z2],'-r');
        x3=x2-((x2-x1)/7);
        y3=y2-((y2-y1)/7);
        z3=z2-((z2-z1)/7); % head of 'arrow' is black
        plot3([x3,x2],[y3,y2],[z3,z2],'-k');
        %plot3(x3,y3,z3,'Ok')
    else % 2D
        x2=get(obj2,'Xdata');
        y2=get(obj2,'Ydata');
        ph=text(x2,y2,s2, 'Interpreter', 'None');
        plot([x1,x2],[y1,y2],'-r');
        x3=x2-((x2-x1)/7);
        y3=y2-((y2-y1)/7); % head of 'arrow' is black
        plot([x3,x2],[y3,y2],'-k');
        %x4=x3-dat.meanXYZ(1)
    end
end
dat=get(obj2,'UserData');
set(ph,'UserData',dat,'UIContextMenu',Data.handles.UIContextMenu);
hold off
Data=ShowImageData(Data,model,positions,fignos,model_img,sourcestr,img,option);

function Data=ShowFullModelImageA(Data,option)
[Data,positions,model,fignos,model_img,sourcestr,img]=ShowImageData(Data);
figure(Data.handles.axesFigA)
if option==0
    disp('Use this baseline shape + shape and appearance of the following image (click on target)');
elseif option==1
    disp('Use this baseline appearance + shape and appearance of the following image');
else
    disp('Use this shape and the following appearance (click on target)');
end
obj1=gco;
if option==0
    disp('Use this baseline shape + shape and appearance of the following image (click on target)');
    s1='B';
    s2='SA';
elseif option==1
    disp('Use this baseline appearance + shape and appearance of the following image');
    s1='A';
    s2='SA';
else
    disp('Use this shape and the following appearance (click on target)');
    s1='S';
    s2='A';
end
hold on
if get(Data.handles.Plot3D2,'Value') % 3D
    x1=get(obj1,'Xdata');
    y1=get(obj1,'Ydata');
    z1=get(obj1,'Zdata');
    text(x1,y1,z1,s1, 'Interpreter', 'None');
else
    x1=get(obj1,'Xdata');
    y1=get(obj1,'Ydata');
    text(x1,y1,s1, 'Interpreter', 'None');
end
disp('Wait for cross-hairs then right-click on the target point');
[a,b]=ginput(1);
obj2=gco;
disp('Head of ''arrow'' is black')
if get(Data.handles.Plot3D2,'Value') % 3D
    x2=get(obj2,'Xdata');
    y2=get(obj2,'Ydata');
    z2=get(obj2,'Zdata');
    text(x2,y2,z2,s2, 'Interpreter', 'None');
    plot3([x1,x2],[y1,y2],[z1,z2],'-r');
    x3=x2-((x2-x1)/7);
    y3=y2-((y2-y1)/7);
    z3=z2-((z2-z1)/7); % head of 'arrow' is black
    plot3([x3,x2],[y3,y2],[z3,z2],'-k');
    %plot3(x3,y3,z3,'Ok')
else % 2D
    x2=get(obj2,'Xdata');
    y2=get(obj2,'Ydata');
    text(x2,y2,s2, 'Interpreter', 'None');
    plot([x1,x2],[y1,y2],'-r');
    x3=x2-((x2-x1)/7);
    y3=y2-((y2-y1)/7); % head of 'arrow' is black
    plot([x3,x2],[y3,y2],'-k');
end
hold off
Data=ShowImageData(Data,model,positions,fignos,model_img,sourcestr,img,option);

%%%%%%
%
% Rotate
%
%%%%%%
function [Data] = Rotate(Data)
%set(Data.handles.Rotate,'Value',0);
set(Data.handles.Zoom,'Value',0);
set(Data.handles.Pan,'Value',0);
set(Data.handles.SelectShape,'Value',0);
set(Data.handles.SelectAppearance,'Value',0);
rotate3d off
figure(Data.handles.axesFig);
rotate3d off
zoom off
pan off
figure(Data.handles.axesFigA);
rotate3d off
zoom off
pan off
set(Data.handles.Rotate,'Value',1);
figure(Data.handles.axesFig);
rotate3d on
figure(Data.handles.axesFigA);
rotate3d on
%%%%%%
%
% Zoom
%
%%%%%%
function [Data] = Zoom(Data)
set(Data.handles.Rotate,'Value',0);
%set(Data.handles.Zoom,'Value',0);
set(Data.handles.Pan,'Value',0);
set(Data.handles.SelectShape,'Value',0);
set(Data.handles.SelectAppearance,'Value',0);
figure(Data.handles.axesFig);
rotate3d off
zoom off
pan off
figure(Data.handles.axesFigA);
rotate3d off
zoom off
pan off
figure(Data.handles.axesFig);
zoom on
figure(Data.handles.axesFigA);
zoom on
%%%%%%
%
% Pan
%
%%%%%%
function [Data] = Pan(Data)
set(Data.handles.Rotate,'Value',0);
set(Data.handles.Zoom,'Value',0);
%set(Data.handles.Pan,'Value',0);
set(Data.handles.SelectShape,'Value',0);
set(Data.handles.SelectAppearance,'Value',0);
figure(Data.handles.axesFig);
rotate3d off
zoom off
pan off
figure(Data.handles.axesFigA);
rotate3d off
zoom off
pan off
figure(Data.handles.axesFig);
pan on
figure(Data.handles.axesFigA);
pan on

%%%%%%
%
% Pcomponents
%
%%%%%%
function [Data] = Pcomponents(Data,option)
if isempty(Data.handles.bAll)
    s = sprintf('There are no data to plot. Try clicking ''%s''.', get(Data.handles.load_btn, 'String'));
    uiwait(msgbox(s,'No Groups','modal'));
    fprintf('Data.handles.bAll is empty, function Pcomponents\n');
    return;
end
M = size(Data.handles.bAll{1},2);
if M==0
   s = sprintf('WARNING: There are no deviations to plot, aborting.', get(Data.handles.load_btn, 'String'));
   uiwait(msgbox(s,'ERROR','modal'));
   return;
end
for i=1:M; str{i} = num2str(i); end;
[temp_Pcomponents,v] = listdlg('Name', 'PC Selection', 'PromptString','Select PC''s:', 'SelectionMode','multiple','ListString',str, 'ListSize', [300 160]);
if isempty(temp_Pcomponents)
    return;
end
if length(temp_Pcomponents)>3
    temp_Pcomponents = temp_Pcomponents(1:3);
end
switch length(temp_Pcomponents)
    case 1
        set(Data.handles.oneDradio, 'Value', 1);
    case 2
        set(Data.handles.twoDradio, 'Value', 1);
    case 3
        set(Data.handles.threeDradio, 'Value', 1);
end
if option==0
    Data.Pcomponents=temp_Pcomponents;
    set(Data.handles.shapeaxistxt, 'String', num2str(temp_Pcomponents));
else
    Data.PcomponentsA=temp_Pcomponents;
end

return
%%%%%%
%
% ChangeColor
%
%%%%%%
function [handles] = ChangeColor(handles)
groups = handles.groups;
v = get(handles.groupList, 'value');
if isempty(v)
    return;
end
color = uisetcolor;
group = groups(v);
group.color = color;
groups(v) = group;
handles.groups = groups;
handles.PC{v}.color = color;
%save([get(handles.AAM, 'modelDirec'), filesep, 'groups.mat'],'groups' );
handles = UpdateGroupInfo(handles);
return;
%%%%%%
%
% DeleteGroup
%
%%%%%%
function [handles] = DeleteGroup(handles)
groups = handles.groups;
v = get(handles.groupList, 'value');
if isempty(v)
    return;
end
groups(v) = [];
handles.groups = groups;
%save([get(handles.AAM, 'modelDirec'), filesep, 'groups.mat'],'groups' );
handles = UpdateGroupInfo(handles);
%%%%%
%
% UpdateGroupInfo
%
%%%%%
function [handles] = UpdateGroupInfo(handles)
groups = handles.groups;
if isempty(groups)
    set(handles.groupMemberList, 'String', '');
    set(handles.groupList, 'String', '');
    return;
else
    indx = get(handles.groupList, 'value');
    if isempty(indx)
        indx = 1;
    end
    indx = min(indx, length(groups));
    group = groups(indx);
    if ~isempty(group.elements)
        for i=1:length(group.elements)
            tempstring{i}=sprintf('%3d) %s',i,group.elements{i});
        end
        set(handles.groupMemberList, 'String', tempstring);
        set(handles.groupMemberList, 'foregroundcolor', group.color);

        set(handles.groupList, 'value', indx);
        set(handles.groupList, 'String', {groups.name});

    else
        set(handles.groupMemberList, 'String', 'No members in this group');
        set(handles.groupList, 'value', indx);
        set(handles.groupList, 'String', {groups.name});
    end
end
%%%%%
%
% ProjectData
%
%%%%%
function [handles] = ProjectData(handles)
direc = [pwd, filesep, get(handles.AAM, 'modelDirec')];
indices = get(handles.AAM, 'pointindices');

template = get(handles.AAM, 'PointModelTemplate');
template_name = get(template, 'name');
[templatepath, template_name, ext, vers] = fileparts(template_name);
if ~exist(['PointModels', filesep, template_name], 'dir')
    fprintf('There are no point models for this template');
    return;
end
pmdirectory = ['PointModels', filesep, template_name];
pdmfile = [direc, filesep, 'mod_pdm.mat'];
if ~exist(pdmfile)
    return;
end
pdm = load(pdmfile);
pdm = pdm.pdm;
[details, detailstr] = getDetails(direc);
max_var_ratio(handles.AAM, direc, pdm, indices, pmdirectory, details, handles.groups);
return;
%%%%%
%
%
%%%%%
function [handles] = ImportGroup(handles)
[filename, pathname] = uigetfile('*.mat', 'Pick Group File');
if isequal(filename,0) || isequal(pathname,0)
    disp('User pressed cancel')
    return;
end
groups = load([pathname, filesep, filename]);
if ~isfield(groups, 'groups')
    uiwait(msgbox('Not a valid groups file.','Error','modal'));
    return;
end
groups = groups.groups;
handles.groups = [handles.groups, groups];
handles = UpdateGroupInfo(handles);
return;
%%%%%
%
%
%%%%%
function [handles] = AddGroup(handles)
template = get(handles.AAM, 'PointModelTemplate');
template_name = get(template, 'name');
[templatepath, template_name, ext, vers] = fileparts(template_name);
if ~exist(['PointModels', filesep, template_name], 'dir')
    fprintf('There are no point models for this template');
    return;
end
files = dir(['PointModels', filesep, template_name,filesep, '*_pm.mat']);
files = {files.name};

% At this point (for now) load the structure that Nicolas provides and use
% that.

[s, v] = CustomListDlg(files);
if isempty(v)
    return;
end
prompt={'Name:'};
name='Group Name';
numlines=1;
defaultanswer={'New group'};
answer=inputdlg(prompt,name,numlines,defaultanswer);
if isempty(answer)
    return;
end
groups = handles.groups;
groups(end+1).name = answer{1};
groups(end).elements = s;
%colours='rgbcmyk';
colours=[...
    1 0 0;...
    0 1 0;...
    0 0 1;...
    1 0 1;...
    0 1 1;...
    1 1 0;...
    0.5 0.5 0.5;...
    0.5 0.0 0.5;...
    0.0 0.5 0.5;...
    1 0.5 0.5;...
    1 0.0 0.5;...
    0.0 1 0.5;...
    0.5 1.0 0.5;...
    0.5 0.0 1.0;...
    0.0 0.5 1.0;...
    ];
cc=rem(length(groups)-1,size(colours,1))+1;
groups(end).color = colours(cc,:);%[1 0 0];
{groups.name}
set(handles.groupList, 'value', 1);
set(handles.groupList, 'String', {groups.name});
handles.groups = groups;
handles = UpdateGroupInfo(handles);
%save([get(handles.AAM, 'modelDirec'), filesep, 'groups.mat'],'groups' );

%%%%%
%
%
%%%%%
function [handles] = UpdateGroups(handles)
handles = LoadGroupInfo(handles); %Data.Pcomponents,option);

%%%%%
%
%
%%%%%
function [groups] = getgroups(aam)
% Now we must load the group information
groups = [];
modelDir = get(aam, 'modelDirec');
[filename, pathname] = uigetfile([modelDir, filesep, '*.mat'], 'Pick a group file.');
if isequal(filename,0) || isequal(pathname,0)
    disp('User pressed cancel')
else
    %groups = Data.handles.groups;
    groups = load([pathname, filesep, filename]);
    groups = groups.groups;
    %Data.handles.groups = groups;
    disp(['User selected ', fullfile(pathname, filename)])
end
return
% model_dir = get(aam, 'modelDirec');
% if exist([model_dir, filesep, 'groups.mat'])
%     groups = load([model_dir, filesep, 'groups.mat']);
%     groups = groups.groups;
% end
% return;
%%%%%
%
%
%%%%%
% function setgroups(aam, groups)
% % Now we must save the group information
% if ~isa(aam, 'AAM') | isempty(groups)
%     return;
% end
% model_dir = get(aam, 'modelDirec');
% if exist([model_dir])
%     save([model_dir, filesep, 'groups.mat'], 'groups');
% end
% return;
%%%%%
%
%
%%%%%
function [details, detailstr] = getDetails(modelDirec, modelname, projectname)
details =[];
detailstr = '';
if nargin>1
    detailstr=sprintf('%s  %s  ',projectname,modelname);

end

if exist(fullfile(modelDirec,'details.mat'))==2
    temp2=load(fullfile(modelDirec,'details.mat'),'-mat');
    details=temp2.details;
    detailstr=sprintf('%s %d pixels, normalised: ',detailstr,details.pixelnumber);
    if details.scaling==1
        detailstr=sprintf('%s Scaling',detailstr);
    end
    if details.rotation==1
        detailstr=sprintf('%s Rotation',detailstr);
    end
    if details.translation==1
        detailstr=sprintf('%s Translation',detailstr);
    end
    if isfield(details,'FractionPCs')
        detailstr=sprintf('%s (%2.1f%% Var.)',detailstr,100*details.FractionPCs);
    end
end
return;
%%%%%
%
%
%%%%%
function [pdm] = loadPDMFromDisk(pdmdirec)
pdmfilename = [pdmdirec, filesep, 'mod_pdm.mat'];
pdm = [];
if exist(pdmfilename)
    pdm = load(pdmfilename);
    pdm = pdm.pdm;
end
return;
%%%%%
%
%
%%%%%
function [sfam] = loadSFAMFromDisk(sfamdirec)
sfamfilename = [sfamdirec, filesep, 'mod_sfam.mat'];
sfam = [];
if exist(sfamfilename)
    sfam = load(sfamfilename);
    sfam = sfam.sfam;
end
return;
%%%%%
%
%
%%%%%
function [handles] = LoadGroupInfo(handles,Pcomponents,option,PcomponentsA)
handles.PC = [];
ldaX = [];
handles.bAll = [];
handles.bVar = [];
handles.bAllA = [];
handles.bVarA = [];
set(gcf,'Pointer','watch');
if nargin<2
    Pcomponents=1:3;
end
if nargin<3
    option='pdm';
end
if nargin<4
    PcomponentsA=1:3;
end
aam = handles.AAM;
PMT = get(aam, 'PointModelTemplate');
templatename = get(PMT, 'name');
[path, templatename, ext, vers] = fileparts(templatename);
if strcmp(option,'pdm')
    pdm = get(aam, 'pdm');
else
    pdm = aam; %get(aam, 'pdm');
end
if isempty(pdm)
    s = sprintf('The pdm is empty');
    uiwait(msgbox(s,'Error','modal'));
    fprintf('Data.handles.bAll is empty, function Pcomponents\n');
    set(gcf,'Pointer','arrow');
    return;
end

if length(Pcomponents)>length(pdm.b)
    Pcomponents = 1:length(pdm.b);
end


Xm = pdm.Xm;
P = pdm.P;
PC = [];
aam = handles.AAM;
modelDirec=get(aam,'modelDirec');
PMT = get(aam, 'PointModelTemplate');
loops=get(PMT,'loops');
edges=loops{1};
templatename = get(PMT, 'name');
[path, templatename, ext, vers] = fileparts(templatename);
% Now we must load the group information
groups = handles.groups;
%groups = getgroups(handles.AAM);
if isempty(groups)
    handles.groups = [];
    handles.PC = [];
    handles = UpdateGroupInfo(handles);
    %s = sprintf('There are no groups to load.\nTry clicking ''Open''');
    %uiwait(msgbox(s,'!!','modal'));
    set(gcf,'Pointer','arrow');
    s = sprintf('There are no groups open. Try clickin ''Open''.');
    uiwait(msgbox(s,'No Groups','modal'));
    fprintf('Data.handles.bAll is empty, function Pcomponents\n');
    return;
end
%---
projectname=get(PMT,'name');
modelDirec=get(aam,'modelDirec'); %temp(1:findstr(temp,'mod_pdm.mat')-2);
modelname=modelDirec(length('StatisticalModels\')+1:end);

[details, detailstr] = getDetails(modelDirec, modelname, projectname);
set(handles.Details_panel,'Title',detailstr);

Data.details=details;

% Load the PDM first
pdm = loadPDMFromDisk([pwd, filesep, modelDirec]);
if isempty(pdm)
    fprintf('Could not load the PDM\n');
    dbstop in PointCloudProcessor.m at 1189;
end

% Load the SFAM first
sfam = loadSFAMFromDisk([pwd, filesep, modelDirec]);
% if isempty(sfam)
%     fprintf('Could not load the SFAM\n');
%     dbstop in PointCloudProcessor.m at 1208;
% end

handles.PC = {};
waitbarh = waitbar(0, 'Loading Groups please wait...');
set(waitbarh, 'renderer', handles.AAM_renderer);
point_set_directory = get(handles.AAM, 'modelDirec');
[point_set_directory, image_set_name, ext, vers] = fileparts(point_set_directory);

numsamples = 0;
for g=1:length(groups)
    numsamples = numsamples + length(groups(g).elements);
end


for g=1:length(groups)
    PC = [];
    bAll=[];
    bAllStdDev = [];
    bVar=[];
    PCA = [];
    bAllA=[];
    bVarA=[];
    LDAPC = [];
    group = groups(g);
    imagenames = group.elements;
    imagefilename={};
    thumbnail={};
    for i=1:length(imagenames)
        timetogo=(g/length(groups));
        name = imagenames{i};%s = regexprep('str', 'expr', 'repstr')
        % replace _ by - because it seems to be impossible to turn of the
        % 'interpreter'
        %waitbar(timetogo, waitbarh, sprintf('Loading %s',regexprep(name,'_','-')));
        waitbar(timetogo, waitbarh, sprintf('Loading %s', group.name));
        set(findall(waitbarh,'type','text'),'Interpreter','none')

        [path, name, ext, vers] = fileparts(name);
        pmfilename = ['PointModels', filesep,templatename, filesep, name,'.mat'];
        % ptsa are the aligned points (very bad naming)
        if exist([pwd, filesep, pmfilename])
            [pts, ptsa, b, ba] = getAlignedPoints(pmfilename, point_set_directory, Data.details, pdm, get(aam, 'pointindices'));

            ldab = pdm.pca.P'*(ptsa-pdm.Xm);
            LDAPC(end+1,:) = ldab(1:min(length(ldab), numsamples));

            if ~isempty(pts)
                bAllStdDev(end+1,:)=b;
                bAll(end+1,:)=ba;
                bVar(end+1,:)=pdm.b;
                PC(end+1,:) = b(Pcomponents); %(1:3)';  % can we take more than first 3
                %disp(b(Pcomponents)')
                imname=[name(1:findstr(name,'_pm')-1),'.jpg'];
                imagefilename{i}=fullfile(modelDirec,imname);
                if ~exist(imagefilename{i})
                    thumbnail{i} = zeros(20,20);
                else
                    thumbnail{i}=imread(imagefilename{i});
                end
                norm_name=fullfile(modelDirec,['Norm_',imname]);
                if isempty(sfam)
                    baA = [];
                    bA = [];
                    bAllA = [];
                    bVarA = [];
                    PCA = [];
                else
                    if ~exist(norm_name)
                        meansize_img = zeros(length(sfam.Am),1);
                    else
                        meansize_img=double(imread(norm_name));
                    end
                    baA = sfam.P'*(meansize_img(:)-sfam.Am);
                    bA = baA./sqrt(sfam.b);
                    bAllA(end+1,:)=baA;
                    bVarA(end+1,:)=sfam.b;
                    PCA(end+1,:) = bA(PcomponentsA); %(1:3)';  % can we take more than first 3
                end

            end
            handles.PC{g}.ptsA = PCA;
            handles.PC{g}.pts = PC;
            ldaX{g} = LDAPC;
            handles.PC{g}.color = group.color;
            handles.PC{g}.name = group.name;
            handles.bAll{g}=bAll;
            handles.bAllStdDev{g}=bAllStdDev;
            handles.bVar{g}=bVar;
            handles.bAllA{g}=bAllA;
            handles.bVarA{g}=bVarA;
            handles.imname{g}=imname;
            handles.thumbnail{g}=thumbnail;
        end

    end
end

if get(handles.lda_chk, 'Value')
    handles.bAll = projectUsingLDA(handles.bAll, ldaX, getDimension(handles));

end

close(waitbarh);
handles.groups = groups;
handles.pdm = pdm;
handles = UpdateGroupInfo(handles);
set(gcf,'Pointer','arrow');
%%%%%%
%
% projectUsingLDA
%
%%%%%%
function PC = projectUsingLDA(PC, X, N)
[Xm, P, b, pcaDat] = lda(X, N);
for i=1:length(X)
    PC{i} = (P*X{i}')';
end
return;
%%%%%%
%
% getAlignedPoints
%
%%%%%%
function [pts, aligned_pts, b, ba] = getAlignedPoints(pmfilename, point_set_directory, details, pdm, indices)
pts = [];
aligned_pts = [];
b = [];
ba = [];
opts = [details.scaling, details.rotation, details.translation];
if exist(pmfilename)
    temp =load(pmfilename);
    pts_temp = temp.pts;
    pts_temp = reshape(pts_temp, 2, length(pts_temp)/2);
    if ~exist([point_set_directory, filesep, 'context.mat'])
        pts = pts_temp(:, indices); % take just the subset of points
    else
        pts = pts_temp(:);
    end
    pts = pts(:); % next align the points using the appropriate details
    aligned_pts= pmalign(pts,opts,pdm.Xm);
    %ba = P'*(pts-Xm);
    ba = pdm.pca.P'*(aligned_pts-pdm.Xm);
    b = ba./sqrt(pdm.pca.b);
end
return;
%%%%%%
%
% Quit
%
%%%%%%
function [handles] = Quit(handles)
delete(handles.figMain);
if isfield(handles, 'axesFig')
    % if exist('handles.axesFig')
    delete(handles.axesFig); %axesFig_handleA
    % end
end
if isfield(handles, 'axesFigA')
    %if exist('handles.axesFigA')

    delete(handles.axesFigA); %axesFig_handleA
    %end
end
handles = rmfield(handles, 'figMain');
return;

function [handles] = TogglePan(handles)
if get(handles.togglePan,'Value')
    axes(handles.mainAxes)
    rotate3d off;
    axes(handles.mainAxesA)
    rotate3d off;
    zoom off;
    set(handles.toggleZoom, 'Value', 0);
    pan on;
else
    axes(handles.mainAxes)
    rotate3d on;
    axes(handles.mainAxesA)
    rotate3d on;
    pan off;
end
return;
%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Toggle Zoom
%
%%%%%%%%%%%%%%%%%%%%%%%%%%
function [handles] = ToggleZoom(handles)
if get(handles.toggleZoom,'Value')
    axes(handles.mainAxes)
    rotate3d off;
    axes(handles.mainAxesA)
    rotate3d off;
    zoom on;
    set(handles.togglePan, 'Value', 0);
else
    axes(handles.mainAxes)
    rotate3d on;
    axes(handles.mainAxesA)
    rotate3d on;
    zoom off;
end
return;
%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% PlotData
%
%%%%%%%%%%%%%%%%%%%%%%%%%%
function [handles] = PlotData(handles, PCS, RegCurveS, SDS, Pcomponents,option,PcomponentsA)
handles.ellipsoid_handle = [];
if get(handles.rawdataradio, 'value')
    bShape = handles.bAll;
else
    bShape = handles.bAllStdDev;
end
if ~isempty(handles.bAll)
    handles =PlotAxisData(handles.axesFig, handles.shape_axes_handle, bShape, handles.UIContextMenu, handles, PCS, RegCurveS, SDS, Pcomponents, option);
end
if ~isempty(handles.bAllA)
    handles =PlotAxisData(handles.axesFigA, handles.app_axes_handle, handles.bAllA, handles.UIContextMenuA, handles, PCS, RegCurveS, SDS, PcomponentsA, option);
end
return;
%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% ResetAxis
%
%%%%%%%%%%%%%%%%%%%%%%%%%%
function ResetAxis(handles)
axis(handles.shape_axes_handle, 'image', 'tight');
axis(handles.app_axes_handle, 'image', 'tight');
%%%%%%
%
%
%%%%%%
function handles = ToggleEllipsoid(handles)
showstr = 'off';
if get(handles.showEllipsoid, 'Value')
    showstr = 'on';
end
for i=1:length(handles.ellipsoid_handle)
   set(handles.ellipsoid_handle(i), 'Visible', showstr); 
end
%%%%%%
%
%
%%%%%%
function handles = OrderGroups(handles)
groups = handles.groups;
if isempty(groups)
    return;
end
[s, order, v] = SelectOrderListDlg('ListString', {groups.name});
for i=1:length(order)
    pos = find(order==i);
    q(i) = pos(1);
end
order = q;
if isempty(s)
    return;
end
modelDirec = get(handles.AAM, 'modelDirec');
save([pwd, filesep, modelDirec, filesep, 'usergrouporder.mat'], 'order');
return;
%%%%%%
%
%
%%%%%%
function handles = ChangeIsoLevel(handles, SurfaceData)

if ishandle(handles.surface_handle)
    delete(handles.surface_handle);
end
PCS = handles.PC;
for region  = 1:length(SurfaceData)
    SD = SurfaceData{region};
    COLOR = PCS{region}.color;
    if (get(handles.showSurface,'value')==1)
        if (size(SD.V,1) > 0) & (length(size(SD.V)) == 3)
            scale = get(handles.radiusSlider, 'value');
            set(handles.radiusTxt, 'String', num2str(scale));
            p = patch(isosurface(SD.X, SD.Y, SD.Z, SD.V, scale), 'Parent', handles.shape_axes_handle);
            isonormals(SD.X, SD.Y, SD.Z, SD.V, p);
            set(p, 'FaceAlpha', .5, 'FaceColor', COLOR, 'EdgeColor', 'none');
            handles.surface_handle = p;
            lighting(handles.shape_axes_handle, 'phong');

        end
    end
end
return;
%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% showSurface
%
%%%%%%%%%%%%%%%%%%%%%%%%%%
function p = showSurface(handles, dat, SD, axis_handle)
COLOR = dat.colour;
if (size(SD.V,1) > 0) & (length(size(SD.V)) == 3)
    scale = get(handles.radiusSlider, 'value');
    set(handles.radiusTxt, 'String', num2str(scale));
    p = patch(isosurface(SD.X, SD.Y, SD.Z, SD.V, scale), 'Parent', axis_handle);
    isonormals(SD.X, SD.Y, SD.Z, SD.V, p);
    set(p, 'FaceAlpha', .5, 'FaceColor', COLOR, 'EdgeColor', 'none');
end
return;
%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% showEllipsoid
%
%%%%%%%%%%%%%%%%%%%%%%%%%%
function h = showEllipsoid(handles, PC, COLOR, axis_handle)
h = [];
if get(handles.oneDradio, 'Value')
    %uiwait(msgbox('Cannot show ellipse for less than 2D.','Not enough PCs','modal'));
    return;
end
if get(handles.twoDradio, 'Value')
    [D, mu] = fit_ellipse(PC, 2);
    h = patch(D.x, D.y, COLOR, 'Parent', axis_handle);
    set(h, 'FaceAlpha', .5, 'FaceColor', COLOR, 'EdgeColor', COLOR, 'Visible', 'off');

    %h = plot(axis_handle, D.x, D.y, '-', 'Color', COLOR);
end
if get(handles.threeDradio, 'Value')
    [D, mu] = fit_ellipse(PC);
    h = surf(axis_handle, D.x, D.y, D.z);
    set(h, 'FaceAlpha', .5, 'FaceColor', COLOR, 'EdgeColor', COLOR, 'Visible', 'off');
end
axis(axis_handle, 'normal');
%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Plot Particles
%
%%%%%%%%%%%%%%%%%%%%%%%%%%
function [handles, dat] = PlotParticles(handles, dat, PC,axis_handle)
COLOR = dat.colour;
region = dat.region;
elements = handles.groups(region).elements;
symbol='o';
if get(handles.threeDradio,'Value')
    ph = plot3dParticles(dat, PC, COLOR, symbol, elements,axis_handle);
    if get(handles.ShowNumbers,'Value')
        plot3dText(PC, COLOR, elements, axis_handle);
    end
end
if get(handles.twoDradio,'Value')
    ph = plot2dParticles(dat, PC, COLOR, symbol, elements,axis_handle);
    if get(handles.ShowNumbers,'Value')
        plot2dText(PC, COLOR, elements, axis_handle);
    end
end
if get(handles.oneDradio,'Value')
    ph = plot1dParticles(dat, PC, COLOR, symbol, elements,axis_handle);
    if get(handles.ShowNumbers,'Value')
        plot2dText(PC, COLOR, elements, axis_handle);
    end
end
%set(ph,'UserData',dat,'UIContextMenu',dat.UIContextMenu);

%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% plot3dText
%
%%%%%%%%%%%%%%%%%%%%%%%%%%
function  plot3dText(PC, COLOR, elements, axis_handle)
for i=1:size(PC,1) % PC images
    dat.filename= elements{i};
    ph=text(PC(i,1), PC(i,2), PC(i,3),sprintf('%3d',i),'color',COLOR, 'BackgroundColor', 'w', 'EdgeColor', COLOR, 'Parent', axis_handle, 'Interpreter', 'None');
    set(axis_handle,'Clipping','off');
    set(ph,'UserData',dat,'UIContextMenu',UIContextMenu);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% plot3dParticles
%
%%%%%%%%%%%%%%%%%%%%%%%%%%
function  ph = plot3dParticles(dat, PC, COLOR, symbol, elements, axis_handle)
for i=1:size(PC,1)
    dat.filename=elements{i};
    dat.num=i;
    dat.PC=PC(i,:);
    %dat.bAll=dat.bAll(i,:);
    %dat.PCA=dat.PCA(i,:);
    %dat.bAllA=dat.bAllA(i,:);
    ph=plot3(axis_handle, PC(i,1), PC(i,2), PC(i,3), symbol, 'MarkerFaceColor', COLOR, 'MarkerEdgeColor', COLOR, 'MarkerSize', 5);
    set(axis_handle,'Clipping','off');
    axis(axis_handle, 'normal', 'tight');
    set(ph,'UserData',dat,'UIContextMenu',dat.UIContextMenu);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% plot2dText
%
%%%%%%%%%%%%%%%%%%%%%%%%%%
function  plot2dText(PC, COLOR, elements, axis_handle)
for i=1:size(PC,1) % PC images
    dat.filename= elements{i};
    text(PC(i,1), PC(i,2), sprintf('%3d',i),'color',COLOR, 'BackgroundColor', 'w', 'EdgeColor', COLOR, 'Parent', axis_handle, 'Interpreter', 'None');
    set(axis_handle,'Clipping','off');
end
%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% plot2dParticles
%
%%%%%%%%%%%%%%%%%%%%%%%%%%
function  ph = plot2dParticles(dat, PC, COLOR, symbol, elements, axis_handle)
bAll = dat.bAll;
for i=1:size(PC,1)
    dat.filename=elements{i};
    dat.num=i;
    dat.PC=PC(i,:);
    dat.bAll=bAll(i,:);
    %    dat.PCA=dat.PCA(i,:);
    %    dat.bAllA=dat.bAllA(i,:);
    ph=plot(axis_handle, PC(i,1), PC(i,2), symbol, 'MarkerFaceColor', COLOR, 'MarkerEdgeColor', COLOR, 'MarkerSize', 5);
    set(axis_handle,'Clipping','off');
    axis(axis_handle, 'normal', 'tight');
    set(ph,'UserData',dat,'UIContextMenu',dat.UIContextMenu);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% plot1dParticles
%
%%%%%%%%%%%%%%%%%%%%%%%%%%
function  ph = plot1dParticles(dat, PC, COLOR, symbol, elements, axis_handle)
bAll = dat.bAll;
sigma = 2*std(PC(:,1));
sigma = max(abs(PC(:,1)));

barheight = 0.1;
mu = mean(PC(:,1));
centerPC = PC(:,1)-mu;
minsigma = abs(min(centerPC));
maxsigma = max(centerPC);

ph=plot(axis_handle, [mu-minsigma mu], [PC(1,2) PC(1,2)], 'Color', COLOR, 'MarkerFaceColor', COLOR, 'MarkerEdgeColor', COLOR, 'MarkerSize', 5, 'LineWidth', 3);
ph=plot(axis_handle, [mu mu+maxsigma], [PC(1,2) PC(1,2)], 'Color', COLOR, 'MarkerFaceColor', COLOR, 'MarkerEdgeColor', COLOR, 'MarkerSize', 5, 'LineWidth', 3);
ph=plot(axis_handle, [mu-minsigma mu-minsigma], [PC(1,2)-barheight PC(1,2)+barheight], 'Color', COLOR, 'MarkerFaceColor', COLOR, 'MarkerEdgeColor', COLOR, 'MarkerSize', 5, 'LineWidth', 3);
ph=plot(axis_handle, [mu+maxsigma mu+maxsigma], [PC(1,2)-barheight PC(1,2)+barheight], 'Color', COLOR, 'MarkerFaceColor', COLOR, 'MarkerEdgeColor', COLOR, 'MarkerSize', 5, 'LineWidth', 3);


function d = getDimension(handles)
d = 1;
if get(handles.oneDradio, 'Value')
    d = 1;
end
if get(handles.twoDradio, 'Value')
    d = 2;
end
if get(handles.threeDradio, 'Value')
    d = 3;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% PlotAxisData
%
%%%%%%%%%%%%%%%%%%%%%%%%%%
function [handles] = PlotAxisData(fig_handle, axis_handle, PCAb, UIContextMenu, handles, PCS, RegCurveS, SDS, Pcomponents,option)
title(axis_handle, 'Right click on a point to show the original image');
if isempty(PCAb{1})
    return;
end
if nargin<5
    Pcomponents=1:3;
end
if nargin<6
    option='pdm';
end

PCS= handles.PC;
if isempty(PCS)
    return;
end
numPCs = size(PCAb{1},2);
if max(Pcomponents)>numPCs
    uiwait(msgbox('Thera are less than 3 Principle Components.','Not enough PCs','modal'));
    if get(handles.threeDradio, 'Value')
        set(handles.twoDradio, 'Value', 1);
    end
    Pcomponents = Pcomponents(1:numPCs);
end
%set(handles.plotDataBtn,'String',sprintf('Plot %d %d %d',Pcomponents(1),Pcomponents(2),Pcomponents(3)));

pmt=get(handles.AAM,'PointModelTemplate');
template=get(pmt,'name');
[templatepath, templatename, ext, vers] = fileparts(template);

RegCurveS = cell(1, length(PCS));
if isempty(SDS)
SDS = cell(1, length(PCS));
for i=1:length(PCS)
    SDS{i}.X = [];
    SDS{i}.Y = [];
    SDS{i}.Z = [];
    SDS{i}.V = [];
end
end
set(fig_handle, 'RendererMode', 'auto');
cla(axis_handle);
hold(axis_handle, 'on');
handles.meanPC = [];
userorderfilename = [get(handles.AAM, 'modelDirec'), filesep, 'usergrouporder.mat'];
order = [];
if get(handles.usergrouporderchk, 'Value') & exist(userorderfilename)
    order = load(userorderfilename);
    if isfield(order, 'order')
        order = order.order;
    end
end
if length(order)<length(PCS)
    order = orderData(PCAb, Pcomponents);
end
wh = waitbar(0, 'Plotting data, please wait...');
for region  = 1:length(PCS)
    waitbar(region/length(PCS), wh);
    clear PC
    clear PCA
    COLOR = PCS{region}.color;
    RegCurve = RegCurveS{region};
    SD = SDS{region};
    dat.region=region;
    dat.groupname=PCS{region}.name;
    dat.colour=COLOR;
    dat.templatename=templatename;
    bAll=PCAb{region};

    dat.Pcomponents=Pcomponents;
    if size(bAll,1)>0
        if get(handles.oneDradio, 'Value')
            variance = handles.pdm.pca.b;
            for i=1:size(bAll,2)
               bAll(:,i) = bAll(:,i)/ sqrt(variance(1));
            end
        end
        PC=bAll(:,Pcomponents);

        mr = mean(bAll,1);
        mr = mr(Pcomponents);
        if size(PC,2)==1
            set(handles.oneDradio, 'Value', 1)
        end
        if get(handles.oneDradio, 'Value')

            %PC = [PC(:,1), region*ones(size(PC,1),1)];
            
            PC = [PC(:,1), order(region)*ones(size(PC,1),1)];
            mr(2) = order(region);
            set(handles.showEllipsoid, 'Value', 0);
        end
        mall=mean(bAll,1);
        dat.meanXYZ=mr;
        if(size(PC,1)>0) & (get(handles.showParticles,'value')==1)
            dat.bAll = bAll;
            dat.UIContextMenu = UIContextMenu;
            dat.PCAb = PCAb;
            [handles, dat] = PlotParticles(handles, dat, PC, axis_handle);
        end
        if (size(RegCurve,1)>0) & (get(handles.showFit,'value')==1)
            plot3(axis_handle, RegCurve(:,1), RegCurve(:,2), RegCurve(:,3), 'r', 'LineWidth', 2);
        end
        if (size(PC,1)>1) & (get(handles.showSurface,'value')==1)
            handles.surface_handle(region) = showSurface(handles, dat, SD, axis_handle);
        end
        %if (size(PC,1)>1) & (get(handles.showEllipsoid,'value')==1)
            ellip_handle = showEllipsoid(handles, PC, COLOR, axis_handle);
            if ~isempty(ellip_handle)
            	handles.ellipsoid_handle(end+1) = ellip_handle;
            end
        %end
        if (size(PC,1)>1) & (get(handles.showConvHull,'value')==1)
            fprintf('Region %d , Size: %d\n', region, size(PC,1));
            C = convhulln(PC);
            for i = 1:size(C,1)
                j = C(i,[1 2 3 1]);
                patch(PC(j,1),PC(j,2),PC(j,3),COLOR,'FaceAlpha',0.6, 'Parent', axis_handle);
            end
        end
        handles.meanPC(region,:)=mr;
        dat.templatename=templatename;
        dat.num=0;
        dat.PC=mr;
        dat.filename=sprintf('Group_%d',region);
        dat.bAll=mall;
        %figure(handles.axesFig);
        if get(handles.showmeanschk, 'Value');
            plotGroupMeans(handles, mr, axis_handle, COLOR);
        end

        if get(handles.showlabelschk, 'Value')
            plotGroupLabels(mr, PCS, COLOR, axis_handle, region, getDimension(handles), order);
        end
        if region==1
            bAllPC=PC;
        else
            bAllPC=[bAllPC;PC];
        end
    end
end
if exist('bAllPC')
    %figure(handles.axesFig);
    meanbAllPC=mean(bAllPC);
    ph = [];
    if get(handles.threeDradio,'Value')
        ph=plot3(axis_handle, meanbAllPC(1), meanbAllPC(2), meanbAllPC(3),'X','MarkerFaceColor',COLOR,'MarkerEdgeColor','k','MarkerSize',12);
    end
    if get(handles.twoDradio, 'Value')
        ph=plot(axis_handle, meanbAllPC(1), meanbAllPC(2),'X','MarkerFaceColor',COLOR,'MarkerEdgeColor','k','MarkerSize',12);
    end
    dat.templatename=templatename;
    dat.num=-1;
    dat.PC=zeros(size(meanbAllPC,2));
    dat.bAll=zeros(size(meanbAllPC,2));
    dat.region=0;
    dat.groupname='Overall Mean';
    dat.filename='Overall_Mean';
    dat.Pcomponents=Pcomponents;
    set(ph,'UserData',dat,'UIContextMenu',UIContextMenu);
    %hold(axis_handle, 'off');
end


set(axis_handle, 'Color', [1 1 1], 'XColor', [0 0 0], 'YColor', [0 0 0], 'ZColor', [0 0 0]);
if get(handles.oneDradio, 'Value')
    view(axis_handle,2);
    axis(axis_handle, 'normal');
    axis(axis_handle, 'tight');
    xlim = get(axis_handle, 'XLim');
    axis(axis_handle, [xlim(1) xlim(2) 0 length(PCS)+1]);
    axis(axis_handle, 'normal');
else
    axis(axis_handle, 'normal', 'equal', 'tight');
end
handles = ToggleEllipsoid(handles);
grid(axis_handle, 'on');
drawnow; hold(axis_handle, 'off');
show_axis_labels(axis_handle, Pcomponents, PCS, handles, order);
%set(axis_handle,'DataAspectRatio',[1 1 1]);
lighting(axis_handle, 'phong');
cameratoolbar(fig_handle);
close(wh);
%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% orderData
%
%%%%%%%%%%%%%%%%%%%%%%%%%%
function order = orderData(PCAb, Pcomponents)

m = [];
for region  = 1:length(PCAb)
    bAll=PCAb{region};
    if size(bAll,1)>0
        PC=bAll(:,Pcomponents);
    end
    m(region) = mean(PC(:,1));
end
[val, ind] = sort(m);
[val, order] = sort(ind);
%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% plotGroupMeans
%
%%%%%%%%%%%%%%%%%%%%%%%%%%
function plotGroupMeans(handles, mr, axis_handle, COLOR)
if get(handles.threeDradio,'Value')
    ph=plot3(axis_handle, mr(1), mr(2), mr(3),'O','MarkerFaceColor',COLOR,'MarkerEdgeColor','k','MarkerSize',10);
end
if get(handles.twoDradio, 'Value')
    ph=plot(axis_handle, mr(1), mr(2),'O','MarkerFaceColor',COLOR,'MarkerEdgeColor','k','MarkerSize',10);
end
if get(handles.oneDradio, 'Value')
    ph=plot(axis_handle, mr(1), mr(2)+.2,'v','MarkerFaceColor',COLOR,'MarkerEdgeColor','k','MarkerSize',10);
end
%set(ph,'UserData',dat,'UIContextMenu',UIContextMenu);

%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% plotGroupLabels
%
%%%%%%%%%%%%%%%%%%%%%%%%%%
function plotGroupLabels(mr, PCS, COLOR, axis_handle, region, dimension, order)
switch dimension
    case 1
        text(mr(1), order(region), PCS{region}.name,'Color','k', 'Parent', axis_handle, 'BackgroundColor', 'w', 'EdgeColor', COLOR, 'FontSize', 8, 'Interpreter', 'None');
    case 2
        text(mr(1), mr(2)+0.3, PCS{region}.name,'Color','k', 'Parent', axis_handle, 'BackgroundColor', 'w', 'EdgeColor', COLOR, 'FontSize', 8, 'Interpreter', 'None');
    case 3
        text(mr(1), mr(2)+0.3, mr(3)+0.3, PCS{region}.name,'Color','k', 'Parent', axis_handle, 'BackgroundColor', 'w', 'EdgeColor', COLOR, 'FontSize', 8, 'Interpreter', 'None');

end


%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% show_axis_labels
%
%%%%%%%%%%%%%%%%%%%%%%%%%%
function show_axis_labels(axis_handle, Pcomponents, PCS, handles, order)
str = 'PC';
if get(handles.lda_chk, 'Value')
    str = 'LDA';
end
if get(handles.oneDradio, 'Value')
    xlabel(axis_handle, sprintf('%s %d',str, Pcomponents(1)));
    ylabel(axis_handle, sprintf('Groups'));
    ytick = 1:length(PCS);
    yticklabel = {};
    for i=1:length(PCS)
        yticklabel{end+1} = PCS{find(order==i)}.name;
    end
    set(axis_handle, 'Ytick', ytick, 'YTickLabel', yticklabel);
end
if get(handles.twoDradio, 'Value') | get(handles.threeDradio, 'Value')
    xlabel(axis_handle, sprintf('%s %d',str, Pcomponents(1)));
    ylabel(axis_handle, sprintf('%s %d',str, Pcomponents(2)));
    ylim = get(axis_handle, 'YLim');
    yticklabel = linspace(ylim(1), ylim(2), 5);
    ytick = yticklabel;
    set(axis_handle, 'Ytick', ytick, 'YTickLabel', yticklabel);
end
if get(handles.threeDradio, 'Value')
    zlabel(axis_handle, sprintf('%s %d',str, Pcomponents(3)));
end

% %%%%%%%%%%%%%%%%%%%%%%%%%%
% %
% % ShowStats
% %
% %%%%%%%%%%%%%%%%%%%%%%%%%%
% function Data = ShowStats(Data)
% Data.FigNo=Data.FigNo+1;
% figure(Data.FigNo);
% clf
% Data.FigNo=Data.FigNo+1;
% c='mkrgbcy';
% first=Data.FigNo;
% Data.FigNo=Data.FigNo+1;
% for i  = 1:length(Data.handles.bAll)
%     figure(Data.FigNo)
%     Data.FigNo=Data.FigNo+1;
%     bAll=Data.handles.bAll{i};%one cell per group,
%     bVar=Data.handles.bVar{i}; %one row per image, one column per component
%     sd=bAll./sqrt(bVar);
%     pc_means=mean(sd);
%     pc_std=std(sd);
%     pc_ste=pc_std/sqrt(size(sd,1));
%     errorbar(Data.FigNo, 1:length(pc_means),pc_means,pc_ste,'color',c(1+rem(i,length(c)-1)));
%     axis([0 length(pc_means) -4 4]);
%     title(Data.handles.groups(i).name);
%     xlabel('Shape Component');
%     figure(first);
%     hold on
%     errorbar(first, 1:length(pc_means),pc_means,pc_ste,'color',c(1+rem(i,length(c)-1)));
%     axis([0 length(pc_means) -4 4]);
%     hold off
% end
% first=Data.FigNo;
% Data.FigNo=Data.FigNo+1;
% for i  = 1:length(Data.handles.bAllA)
%     figure(Data.FigNo)
%     Data.FigNo=Data.FigNo+1;
%     bAll=Data.handles.bAllA{i};%one cell per group,
%     bVar=Data.handles.bVarA{i}; %one row per image, one column per component
%     sd=bAll./sqrt(bVar);
%     pc_means=mean(sd);
%     pc_std=std(sd);
%     pc_ste=pc_std/sqrt(size(sd,1));
%     errorbar(1:length(pc_means),pc_means,pc_ste,'color',c(1+rem(i,length(c)-1)));
%     axis([0 length(pc_means) -4 4]);
%     title(Data.handles.groups(i).name);
%     xlabel('Appearance Component');
%     figure(first);
%     hold on
%     errorbar(1:length(pc_means),pc_means,pc_ste,'color',c(1+rem(i,length(c)-1)));
%     axis([0 length(pc_means) -4 4]);
%     hold off
% end
%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% PrintEPS
%
%%%%%%%%%%%%%%%%%%%%%%%%%%
% function [] = PrintEPS(handles);
% [filename, pathname] = uiputfile('*.eps', 'Choose filename');
% if isequal(filename,0) | isequal(pathname,0)
%     disp('User cancelled printing.')
% else
%     h1=figure;
%     set(h1, 'Visible', 'off', 'Units',get(handles.figMain, 'Units'),  'Position', get(handles.figMain, 'Position'), 'PaperPositionMode','auto');
%     h2=copyobj(handles.mainAxes,h1); % handles.axes1 is the object to print
%     set(h1, 'PaperType', 'A4', 'PaperPosition', [0 0 10 10]);
%     fn = [pathname, filesep, filename];
%     print(h1, '-depsc', fn);
%     close(h1);
% end
%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Load Data
%
%%%%%%%%%%%%%%%%%%%%%%%%%%
% function [PC] = LoadData()
% [filename, pathname] = uigetfile('*.pcdat', 'Please open a valid Point Cloud File');
% if isequal(filename, 0) | isequal(pathname, 0)
%     disp('User pressed cancel to load point cloud');
%     PC = [];
% else
%     PC = load([pathname, filesep, filename],'-mat');
%     PC = PC.Data;
%     if strcmp(class(PC), 'double')
%         PC{1}.pts =PC;
%         PC{1}.color = [1 0 0];
%     end
% end
%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FitData
%
%%%%%%%%%%%%%%%%%%%%%%%%%%
function [REGCURVES, SDS] = FitData(handles, PCS, Pcomp, fitOpt)
for region = 1:length(PCS)
    PC = PCS{region};
    if ~isempty(PC)
        PC = PC(:, Pcomp);
        if ~isempty(PC)
            [X,Y,Z] = poly3dfit(PC(:,1), PC(:,2), PC(:,3), [get(handles.xpolyorder, 'value'),get(handles.ypolyorder, 'value')], 1000);
            RegCurve = [X, Y, Z];
            dx = str2num(get(handles.dx_txt, 'String'));
            dy = str2num(get(handles.dy_txt, 'String'));
            dz = str2num(get(handles.dz_txt, 'String'));
            xmin = min(PC(:,1))-1; xmax = max(PC(:,1))+1;
            ymin =min(PC(:,2))-1; ymax = max(PC(:,2))+1;
            zmin = min(PC(:,3))-1; zmax = max(PC(:,3))+1;
            [X, Y, Z] = meshgrid([xmin:dx:xmax],[ymin:dy:ymax],[zmin:dz:zmax] );
            V = zeros(size(X));
            waitbar_h = waitbar(0, 'Estimating surface...');
            indx =1;
    N = 2;

    for i=1:size(X,1)
        for j=1:size(X,2)
            for k = 1:size(X,3)
                waitbar(indx/(size(X,1)*size(X,2)*size(X,3)), waitbar_h); indx = indx+1;
                if get(handles.fitToData, 'value')
                    d = sqrt((PC(:,1) - X(i,j,k)).^2 + (PC(:,2) - Y(i,j,k)).^2 + (PC(:,3) - Z(i,j,k)).^2 );
					[d, ind] = sort(d);
                    V(i, j, k) = 1/((sum(d(1:N))/N));
                 else
                     [v, p] = min(sqrt((RegCurve(:,1) - X(i,j,k)).^2 + (RegCurve(:,2) - Y(i,j,k)).^2 + (RegCurve(:,3) - Z(i,j,k)).^2 ));
                     V(i, j, k) = sum(sqrt((RegCurve(p,1) - X(i,j,k)).^2 + (RegCurve(p,2) - Y(i,j,k)).^2 + (RegCurve(p,3) - Z(i,j,k)).^2 )) ;
                 end
            end
        end
    end
    close(waitbar_h);
    V = smooth3(V, 'gaussian',11, 1);
    V = smooth3(V, 'gaussian',11, 1);
    V = smooth3(V, 'box', 3);
    REGCURVES{region} = RegCurve;
    SDS{region}.X = X;
    SDS{region}.Y = Y;
    SDS{region}.Z = Z;
    SDS{region}.V = V;
        end
    else
        SDS{region}.X = [];
        SDS{region}.Y = [];
        SDS{region}.Z = [];
        SDS{region}.V = [];
    end
end
% function [data] = fit_ellipse(PC);
%
% Input: PC is a Nx(2 or 3) matrix of the raw 3D scatter data
% Output: data is a structure containing either
%             data.x and data.y for 2D input
%             and
%             data.x, data.y and data.z for 3D input
% Dr. A. I. Hanna (2005)
function [data, mean_PC] = fit_ellipse(PC, dim);
% Check to see that there is more than one data point.
if nargin<2
    dim = 3;
end
if dim==2
    PC = PC(:,1:2);
end
if(size(PC,1) < 2)
    return;
end
mean_PC = mean(PC);
PC = PC - ones(size(PC,1),1)*mean(PC);
[V, D] = eig(PC'*PC);
ProjC = (V'*PC')';
smax = max(abs(ProjC));
if dim == 2

    theta = linspace(0, 2*pi, 100);
    x = cos(theta);
    y = sin(theta);
    C = [x(:)*smax(1) y(:)*smax(2)];
    C = (V*C')';
    x = C(:,1); y = C(:,2);;
    data.x = x + mean_PC(1);
    data.y = y + mean_PC(2);
else
    % Create a sphere
    [x,y,z] = sphere(20);
    [mm,nn] = size(x);
    C = [x(:)*smax(1) y(:)*smax(2) z(:)*smax(3)];
    C = (V*C')';
    x = C(:,1); y = C(:,2); z = C(:,3);
    x = reshape(x,mm,nn);
    y = reshape(y,mm,nn);
    z = reshape(z,mm,nn);
    data.x = x + mean_PC(1);
    data.y = y + mean_PC(2);
    data.z = z + mean_PC(3);
end
return;
% function [X, Y, Z] = poly3dfit(x, y, z, O, N)
%
% Input Params:
%            x - X Coordinates
%            y - Y Coordinates
%            z - Z Coordinates
%            O - a 1x2 vector conatining the order of the polynomials used
%                for the 'y' and 'z' axis respectively, (default O = [2,2]).
%            N - the number of data points used in the interpolation,
%                (default N = 100)
%
% Output Params :
%           X, Y, Z - the N coordinates of the fitted curve.
%
% Example:
%         x = -10:.03:10;%
%         y = exp(x);
%         z = abs(x);
%         y = y + .3*randn(size(y));
%         z = z + 1*randn(size(z));
%         plot3(x, y, z, 'o'); grid on; hold on;
%         [X,Y,Z] = poly3dfit(x, y, z, [5,10], 100);
%         plot3(X, Y, Z, 'r', 'LineWidth', 2);
%
%  Dr. A. I. Hanna 2005.
function [X, Y, Z] = poly3dfit(x, y, z, O, N)
if (nargin < 3)
    fprintf('Must have at least 3 input vectors, x, y, z!!!');
    X = []; Y = []; Z = [];
    return;
end
x = x(:); y = y(:); z = z(:);
if (nargin == 3)
    O = [2, 2];
    N = 100;
end
if (nargin == 4)
    O = O(:);
    if(length(O) >2)
        fprintf('O only needs 2 components!!!');
        X = []; Y = []; Z = [];
    end
    N = 100;
end
% Caclulate the best fit against the z-data
[pz, s] = polyfit(x, z, O(1));
fz = polyval(pz, x);
% Caclulate the best fit against the y-data
[py, s] = polyfit(x, y, O(2));
fy = polyval(py, x);
X = linspace(min(x), max(x), N);
Y = polyval(py, X);
Z = polyval(pz, X);
X = X(:); Y = Y(:); Z = Z(:);
return;
%%%%%%
%
%
%%%%%%
function max_var_ratio(aam, modelDirec, pdm, indices, pmdirectory, details, groups)
% function [P, b] = max_var_ratio(modelDirec, P, b)
%
% Dr. A. I. Hanna (2006)
X = [];
%groups_file = [modelDirec, filesep, 'groups.mat'];
if ~isempty(groups)
    %groups = load(groups_file);
    %groups = groups.groups;
    str = {groups.name};
    [s,v] = listdlg('PromptString','Select a group to project:', 'SelectionMode','multiple','ListString',str)
    if v == 0
        return;
    end
    indx = 1;
    for g = 1:length(s)
        elements = groups(s(g)).elements;
        for ii=1:length(elements)
            name=  elements{ii};
            name =name(1:end-7);
            name = [pwd, filesep, pmdirectory, filesep, name, '_pm.mat'];
            if exist(name)
                [pts, ptsa, b, ba] = getAlignedPoints(name, pmdirectory, details, pdm, indices);
                X(:,indx) = ba(:)';
                indx = indx + 1;
            end

        end
    end
end
filenames = get(aam, 'elements');
indx = 1;
modelX = [];
for e =1:length(filenames)
    element = filenames{e};
    %pts = get(element, 'alignedPts');
    name = get(element, 'imagefilename');
    [pathname, name, ext, vers] = fileparts(name);
    name = [pwd, filesep, pmdirectory, filesep, name, '_pm.mat'];
    if exist(name)
        [pts, ptsa, b, ba] = getAlignedPoints(name, pmdirectory, details, pdm, indices);
        pts = ptsa;
        modelX(:, indx) = pdm.pca.P'*(pts(:) - pdm.Xm);
        indx = indx+1;
    else
        fprintf('%s doesn''t exist.\n', name);
    end
end
if isempty(modelX)
    fprintf('WARNING: Cannot project this data onto this model\n');
    return;
end
for a = 1:size(X,1)
    varX = var(X(a,:));
    if varX == 0
        r(a) = intmax;
    else
        r(a) = pdm.pca.b(a)/varX;
    end
end
[sr, ind] = sort(r, 'descend');
% Just use the first 10.
ind = ind(1:min(24, length(ind)));
N = 3;
M = 4;
figure;
indx = 1;
for PC_to_draw = ind
    leo = X(PC_to_draw, :);
    rem = modelX(PC_to_draw, :);
    range = [min(min(leo), min(rem)), max(max(leo), max(rem))];
    [gxl, gyl] = gaussian_1d(mean(leo), std(leo), range, 100);
    [gxr, gyr] = gaussian_1d(mean(rem), std(rem), range, 100);
    if indx>(M*N)
        indx = 1;
        figure;
    end
    subplot(M, N, indx); hold on; indx = indx + 1;
    plot(leo, zeros(length(leo),1)+0.01, 'o', 'MarkerFaceColor', 'r', 'MarkerEdgeColor', 'r');
    plot(gxl, gyl, 'r-');
    plot(rem, zeros(length(rem),1), 'ro', 'MarkerFaceColor', 'b', 'MarkerEdgeColor', 'b');
    plot(gxr, gyr, 'b-');
    title(sprintf('PC: %d, ratio: %f\n', PC_to_draw, r(PC_to_draw)));
    grid on;
    axis tight;
end
return;
%%%%%%
%
%
%%%%%%
function projectLDAviaPCA(groups, pmdirectory, pdm, indices, details)
% function [P, b] = max_var_ratio(modelDirec, P, b)
%
% Dr. A. I. Hanna (2006)
G = {};
P = pdm.pca.P;

for g=1:length(groups)
    group = groups(g);
    elements = group.elements;
    X = [];
    for e = 1:length(elements)
        ptsfilename = [pwd, filesep, pmdirectory, filesep, elements{e}];
        [pts, ptsa, b, ba] = getAlignedPoints(ptsfilename, pmdirectory, details, pdm, indices);
        b = P'*(ptsa-pdm.Xm);
        X = [X, b];
    end
    G{end+1} = X;
end
for p=1:size(P,2)
    GP = projectGroupsOntoPC(G, p);
    Sb = calcBetweenVariance(GP);
    Sw = calcWithinVariance(GP);
    varRatio(p) = Sb/Sw;
end
figure;
x = 1:size(P,2);
stem(x, varRatio, '-o', 'filled', 'MarkerFaceColor', 'y');
grid on; axis tight;
xlabel('PC');
ylabel('Ratio of variance S_b/S_w')
return;

%%%%%%
%
%
%%%%%%
function GP = projectGroupsOntoPC(G,ind)
for i=1:length(G)
    g = G{i};
    GP{i} = g(ind,:);
end
return;

%%%%%%
%
%
%%%%%%
function Sb = calcBetweenVariance(GP)
for i=1:length(GP)
    g = GP{i};
    mu(i) = mean(g,2);
end
Sb = var(mu);
return
%%%%%%
%
%
%%%%%%
function Sw = calcWithinVariance(GP)
Sw = 0;
for i=1:length(GP)
    g = GP{i};
    Sw(i) = var(g);
end
Sw = mean(Sw);
return


%%%%%
%
%
%%%%%
function [handles] = CalcProjVar(handles, details)
direc = [pwd, filesep, get(handles.AAM, 'modelDirec')];
[point_set_directory, image_set_name, ext, vers] = fileparts(direc);
[details, detailstr] = getDetails(direc);
templatename = get(handles.AAM, 'PointModelTemplate');
templatename = get(templatename, 'name');
[templatepath, templatename, ext, vers] = fileparts(templatename);
[groupindx,v] = listdlg('Name', 'Group Selection', 'PromptString','Select Base Group''s:', 'SelectionMode','multiple','ListString',{handles.groups.name}, 'ListSize', [300 160]);
if isempty(groupindx)
    return;
end
X = [];
for i=groupindx
    grp = handles.groups(i);
    X = [X, getGroupPointModels(grp, templatename)];
end
details = ProcrustesParamsDlg;
details = details(1:3);

[Xaligned, proc_param] = pmalign(X, details);
[Xm, P, b, pcaDat] = principle_component_analysis(Xaligned, 1);

M = length(pcaDat.b);
for i=1:M; str{i} = num2str(i); end;
[temp_Pcomponents,v] = listdlg('Name', 'PC Selection', 'PromptString','Select PC''s:', 'SelectionMode','multiple','ListString',str, 'ListSize', [300 160]);
if isempty(temp_Pcomponents)
    return;
end
A = pcaDat.P(:,temp_Pcomponents);
[groupindx,v] = listdlg('Name', 'Group Selection', 'PromptString','Select Group''s:', 'SelectionMode','multiple','ListString',{handles.groups.name}, 'ListSize', [300 160]);
fprintf('\n\n\n');

if size(Xaligned,2)==1
    C = 0;
else
    C = cov(Xaligned');
end
[V, L] = eig(C);
basetotalVar = trace(L);
fprintf('Model Data\n----------\n');
for a = 1:size(A,2)
    baseProj(:,a) = A(:,a)'*Xaligned;
    baseV(a) = var(baseProj(:,a));
    str = sprintf('\t PC: %d, var: %f, %%: %f\n', a, baseV(a), 100*(baseV(a)/basetotalVar));
    fprintf('%s', str);
end
fprintf('\nTotal var %f, total %%: %f\n\n\n', sum(baseV), 100*(sum(baseV)/basetotalVar));

fprintf('Group Data\n----------\n');
G = {};
for i=1:length(temp_Pcomponents)
   G{i+1, 1} = sprintf('PC %d', temp_Pcomponents(i)); 
end

indx = 2;
for i=groupindx
    grp = handles.groups(i);
    fprintf('Group: %s\n---------\n', grp.name);
    grpX = getGroupPointModels(grp, templatename);
    grpXaligned = pmalign(grpX, details);
    if size(grpXaligned,2)==1
        C = 0;
    else
        C = cov(grpXaligned');
    end
    [V, L] = eig(C);
    totalVar = trace(L);
    G{1 ,indx} = grp.name;
    for a = 1:size(A,2)
        grpProj = A(:,a)'*grpXaligned;

        grpV(a) = var(grpProj);
        P_t = ttest(baseProj(:,a), grpProj');
        str = sprintf('\t PC: %d, var: %f, %%: %f, P(t): %f\n', a, grpV(a), 100*(grpV(a)/totalVar), P_t);
        G{a+1 ,indx} = P_t;
        G{a+1 ,indx+1} = grpV(a);
        
        fprintf('%s', str);
    end
    projVar = sum(grpV);
    if totalVar == 0
        totalVar = 10000;
    end
    perc = 100*(projVar/totalVar);
    indx = indx + 2;
    fprintf('\tTotal var: %f, total %%: %f%%\n', projVar, perc);
end
fprintf('\n\n\n');

btnname=questdlg('Save to table?', 'Save', 'Yes','No','Yes');
if strcmp(lower(btnname), 'yes')
    saveT_TestTable(G);
end
return;
%%%%%%
%
%
%%%%%%
function saveT_TestTable(G)
[filename, pathname] = uiputfile('*.xls', 'Save as...');
if isequal(filename,0) || isequal(pathname,0)
    return;
end
[success,msg]=xlswrite([pathname, filesep, filename],G);
return
%%%%%%
%
%
%%%%%%
function [X] = getGroupPointModels(group, templatename)
N = length(group.elements);
X = [];
for i=1:N
    pmfilename = [pwd, filesep, 'PointModels', filesep,templatename, filesep, group.elements{i}];
    if exist(pmfilename)
        temp =load(pmfilename);
        if isfield(temp, 'pts')
            pts = temp.pts;
            pts = reshape(pts, 2, length(pts)/2);
            X(:,end+1) = pts(:);
        end
    end
end
return;
