function figure_2_pc_effects(modeldir)
%modeldir = 'C:\Documents and Settings\aih\My Documents\PostDoc\Shape_models\ModelsPrjs\Plants\Sandra\PRJ_Arabid_Antirrh\StatisticalModels\LE50Template_04_05_06\Set_1\ImageList_3';
if nargin<1
    modeldir = '';
end
if isempty(modeldir)
    modeldir = uigetdir(pwd, 'Pick your model directory (it must contain a mod_pmd.mat)');
end

pdmfilename = [modeldir, filesep, 'mod_pdm.mat'];
sfamfilename = [modeldir, filesep, 'mod_sfam.mat'];
aammodel = [modeldir, filesep, 'model.aam_dat'];
if ~exist(aammodel, 'file')
    fprintf('There was no model in this directory\n');
    return;
end
aammodel = load(aammodel, '-mat');
aammodel = aammodel.aam;

if ~exist(pdmfilename, 'file')
    fprintf('There was no model in this directory\n');
    return;
end
pdm = load(pdmfilename);
pdm = pdm.pdm;

b = pdm.pca.b;
str = cell(length(b),1);
for i=1:length(b); 
    str{i} = num2str(i); 
end;
[pcs,v] = listdlg('PromptString','Select PC''s to display:', 'SelectionMode','multiple', 'ListString',str);
if isempty(pcs)
    return;
end

figure; set(gcf, 'Color', [1 1 1]);
M = 3; N = 5;

template = get(aammodel, 'PointModelTemplate');
primaries = get(template, 'primaries');
edge = get(template, 'loops');
edge = edge{1};

 
usemodelmean=questdlg('Do you want to use the mean in the model?', 'Question', 'Yes','No','Yes');
if strcmpi(usemodelmean, 'yes')
    xm = pdm.Xm;
    xm = reshape(xm, 2, length(xm)/2);
else
    [filename, pathname] = uigetfile('*.mat', 'Pick an PM-file');
    if isequal(filename,0) || isequal(pathname,0)
        xm = pdm.Xm;
        xm = reshape(xm, 2, length(xm)/2);
        disp('User pressed cancel using mean shape from model')
    else
        disp(['User selected ', fullfile(pathname, filename)]);
        userxm = load([pathname, filesep, filename]);
        if ~isfield(userxm, 'pts');
            xm = pdm.Xm;
            xm = reshape(xm, 2, length(xm)/2);
            disp('Invalid pm file not pts field, using mean shape from model')
        else
            userxm = userxm.pts;
            if length(userxm)~=length(pdm.Xm)
                xm = pdm.Xm;
                xm = reshape(xm, 2, length(xm)/2);
                disp('Invalid pm file wrong number of pts, using mean shape from model')
            else
                xm = userxm;
                xm = reshape(xm, 2, length(xm)/2);
            end
        end
    end
end
offset = mean(xm,2);
xm(1,:)= xm(1,:) - offset(1);
xm(2,:)= xm(2,:) - offset(2);

xcol = [1 1 1];
ycol = [1 1 1];

prompt={'Enter vertical spacing: '};
numlines=1;
defaultanswer={'250'};
options.WindowStyle='normal';
options.Interpreter='tex';
vert=inputdlg(prompt,'',numlines,defaultanswer,options);
vert = str2double(vert{1});

%vert = input('Enter vertical spacing (i.e. 250): ');
if ~isnumeric(vert) || isempty(vert)
    vert = 250;
end
prompt={'Enter horizontal spacing: '};
numlines=1;
defaultanswer={'250'};
options.WindowStyle='normal';
options.Interpreter='tex';
horiz=inputdlg(prompt,'',numlines,defaultanswer,options);
horiz = str2double(horiz{1});
%horiz = input('Enter horizontal spacing (i.e. 420): ');
if ~isnumeric(horiz) || isempty(horiz)
    horiz = 420;
end
ah = gca;
set(ah, 'XColor', xcol, 'YColor', ycol);

fillmodel =questdlg('Do you want to flood fill your shape?', 'Question', 'Yes','No','Yes');
if strcmpi(fillmodel, 'yes')
    fillflag = 1;
else
    fillflag = 0;
end

convhullquest =questdlg('Do you want to use the convex hull for your shape?', 'Question', 'Yes','No','Yes');
if strcmpi(convhullquest, 'yes')
    convhullflag = 1;
else
    convhullflag = 0;
end
prompt={'Specify the line width:'};
name='Line Width';
numlines=1;
defaultanswer={'1'};

answer=inputdlg(prompt,name,numlines,defaultanswer);
if isempty(answer)
    linewidth = 1;
else
    linewidth = str2double(answer{1});
end



for j=1:length(pcs)
    %plotshape(xm(:), primaries, edge, [0*horiz (j-1)*vert], ah, [],
    %fillflag);PAUL
    plotshape(xm(:), primaries, edge, [0*horiz (j-1)*vert], ah, [], fillflag,[0 0 0],linewidth);
end

PCNUM = repmat(pcs, 2, 1);
PCNUM = PCNUM(:);

OFFSET = [];
for i=1:length(pcs)
    OFFSET = [OFFSET; horiz (i-1)*vert; -horiz (i-1)*vert];
end
for i=[2 3]
    for j=1:length(pcs)
        %plotshape(xm(:), primaries, edge, [i*horiz (j-1)*vert], ah,[],
        %0);paul
        plotshape(xm(:), primaries, edge, [0*horiz (j-1)*vert], ah, [], fillflag,[0 0 0],linewidth);
    end
end
prompt={'Enter SD Range: '};
numlines=1;
defaultanswer={'2'};
options.WindowStyle='normal';
options.Interpreter='tex';
SDLIMIT=inputdlg(prompt,'',numlines,defaultanswer,options);
SDLIMIT = str2double(SDLIMIT{1});
%SDLIMIT = input('Enter SD range (i.e. 2): ');
if ~isnumeric(SDLIMIT) || isempty(SDLIMIT)
    SDLIMIT = 2;
end
SD = ones(1,length(pcs));
SD = [SD; -1*SD];
SD = SD(:);
SD = SDLIMIT*SD;
bcol = [0 0 1];
rcol = [1 0 0];
cols = [];
for i=1:2:size(OFFSET,1)
    cols = [cols; rcol; bcol];
end
for i=1:size(OFFSET,1)
    b = SD(i)*sqrt(pdm.pca.b(PCNUM(i)));
    x = xm(:) + pdm.pca.P(:,PCNUM(i))*b;
    x = reshape(x, 2, length(x)/2);
    x = x - mean(x,2)*ones(1, size(x,2));
    plotshape(x(:), primaries, edge, OFFSET(i,:), ah, [], fillflag, cols(i,:));
    plotshape(x(:), primaries, edge, [abs(OFFSET(i,1)), OFFSET(i,2)].*[2 1], ah,[],0, cols(i,:));
    if convhullflag
        p = x;
        K = convhull(p(1,:), p(2,:));
        p = p(:,K);
        a = polyarea(p(1,:), p(2,:));
        ma = polyarea(xm(1,K), xm(2,K));
    else
        a = polyarea(x(1,:), x(2,:));
        ma = polyarea(xm(1,:), xm(2,:));
        
    end
    x = x./sqrt(a);
    x = x*sqrt(ma);
    x = normalizeShapeSize(xm, x); 
    plotshape(x(:), primaries, edge, [abs(OFFSET(i,1)), OFFSET(i,2)].*[3 1], ah,[],0, cols(i,:));
end
xcoef = -2;

%% calculate the area of effect for SD
%number of pc's is
numPCs = size(pdm.pca.b,1);
%for each pc loop and calculate the areas for the + and - SD values
for i = 1:numPCs
    b = SDLIMIT*sqrt(pdm.pca.b(i));
    x = xm(:) + pdm.pca.P(:,i)*b;
    x = reshape(x, 2, length(x)/2);
    x = x - mean(x,2)*ones(1, size(x,2));
    %possSDareas is array of areas of sd
    possSDareas(i) = polyarea(x(1,:), x(2,:));
    
    b = -SDLIMIT*sqrt(pdm.pca.b(i));
    x = xm(:) + pdm.pca.P(:,i)*b;
    x = reshape(x, 2, length(x)/2);
    x = x - mean(x,2)*ones(1, size(x,2));
    negSDareas(i) = polyarea(x(1,:), x(2,:));    
end
%calculate the area abs differences between positive and negative SD's for
%each PC
areaDiff = abs(possSDareas - negSDareas);
%%
for i=1:length(pcs)
    text(xcoef*horiz, (i-1)*vert, ['PC', num2str(pcs(i))], 'HorizontalAlignment', 'center');
end

text(-1.3*horiz, -0.5*vert, sprintf('-%f SD', SDLIMIT), 'HorizontalAlignment', 'center');
text(0, -0.5*vert, 'Mean', 'HorizontalAlignment', 'center');
text(horiz, -0.5*vert, sprintf('%f SD', SDLIMIT), 'HorizontalAlignment', 'center');
text(2*horiz, -0.6*vert, sprintf('Non\nNormalised'), 'HorizontalAlignment', 'center');
text(3*horiz, -0.6*vert, sprintf('Area\nNormalised'), 'HorizontalAlignment', 'center');
text(4*horiz, -0.5*vert, sprintf('Variance'), 'HorizontalAlignment', 'center');
text(5*horiz, -0.6*vert, sprintf('Area\nEffect'), 'HorizontalAlignment', 'center');

for i=1:length(pcs)
    text(4*horiz, (i-1)*vert, sprintf('%2.2f %%\n', 100*(pdm.pca.b(pcs(i))/sum(pdm.pca.b))), 'HorizontalAlignment', 'center');
end

%display area of effect
for i=1:length(pcs)
    text(5*horiz, (i-1)*vert, sprintf('%2.2f %%\n', 100*(areaDiff(i)/sum(areaDiff))), 'HorizontalAlignment', 'center');
end
% text(4*horiz, vert, sprintf('%2.2f %%\n', 100*(pdm.pca.b(1)/sum(pdm.pca.b))), 'HorizontalAlignment', 'center');
% text(4*horiz, 0, sprintf('%2.2f %%\n', 100*(pdm.pca.b(2)/sum(pdm.pca.b))), 'HorizontalAlignment', 'center');
% text(4*horiz, -vert, sprintf('%2.2f %%\n', 100*(pdm.pca.b(3)/sum(pdm.pca.b))), 'HorizontalAlignment', 'center');

axis image tight ij;

%axis([-1000 2500 -400 500]);

%axis fiddles


%%%%%
%
%
%%%%%

function plotshape(pts, pri, edge, offset, axish, fillcol, fillflag, edgecol, linewidth)
[s, o] = isAugmentedProject;
pts = reshape(pts, 2, length(pts)/2);
[pts] = splitPoints(pts', s, o);
pts = pts';
pts = pts(:);
facecol = [1 1 1];
priedgecol = [0 0 0];
prifacecol = [0 0 0];

if nargin<8
    edgecol = [0 0 0];
end
if nargin<9
    linewidth = 1;
end
pts = reshape(pts, 2, length(pts)/2);
pts(1,:) = pts(1,:) + offset(1);
pts(2,:) = pts(2,:) + offset(2);
hold(axish, 'on');
if fillflag
    fillcol = [0.7529    0.7529    0.7529];
    fill(pts(1,:),pts(2,:),fillcol)
end 
for i=1:length(edge)
    xp = [pts(1,edge(i,1)) pts(1,edge(i,2))];
    yp = [pts(2,edge(i,1)) pts(2,edge(i,2))];
    plot(axish, xp, yp, '-', 'Color', edgecol, 'LineWidth', linewidth);
end

function new_pts = normalizeShapeSize(xm, pts)
[s, o] = isAugmentedProject;
xm_P = partitionPts(xm', s);
P = partitionPts(pts', s);
new_pts = [];
for i=1:length(xm_P)
    new_pts = cat(1, new_pts, normalize_size(P{i},  shape_size(xm_P{i})));
end
new_pts = new_pts';
return;

function pts = normalize_size(pts, new_s)
mu = mean(pts);
pts = pts - ones(size(pts,1),1)*mu;
s = shape_size(pts);
pts = pts.*(new_s/s);
pts = pts + ones(size(pts,1),1)*mu;
return

function s = shape_size(pts)
pts = pts - ones(size(pts,1),1)*mean(pts);
s = sqrt(sum(sum(pts.^2, 2)));

function P = partitionPts(pts, segpoint)
segpoint = [0; segpoint; size(pts,1)];
P = {};
for i=2:length(segpoint)
    P{i-1} = pts(segpoint(i-1)+1:segpoint(i),:);
end