function [P, b] = max_var_ratio(modelDirec, P, b, indices, pmdirectory)
% function [P, b] = max_var_ratio(modelDirec, P, b)
%
% Dr. A. I. Hanna (2006)
X = [];
groups_file = [modelDirec, filesep, 'groups.mat'];
if exist(groups_file)
    groups = load(groups_file);
    groups = groups.groups;
    str = {groups.name};
    [s,v] = listdlg('PromptString','Select a group to project:', 'SelectionMode','single','ListString',str)
    if v == 0
        return;
    end
    elements = groups(s(1)).elements;
    indx = 1;
    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 = load(name);
            pts = pts.pts;
            pts = reshape(pts, 2, length(pts)/2);
            pts = pts(:, indices);
            X(:,indx) = pts(:)';
            indx = indx + 1;

        end

    end
end

indx = 1;
for e =1:length(filenames)
    name = filenames{e};
    name = name(1:length(name)-4);
    %pts = load(['PointModels', filesep, templatename, filesep,  name, '_aligned']);
    pts = load(fullfile(modelDirec,[name,'_aligned']));
    pts = pts.pts;
    pts = reshape(pts, 2, length(pts)/2);
    pts = pts(:, indices);
    modelX(:, indx) = pts(:);
    indx = indx+1;
end



for a = 1:size(P,2)
    lv(a,:) =  P(:,a)'*X;
   lv(a, :) = lv(a,:) - mean(lv(a,:));
   modelv(a,:) = P(:,a)'*modelX;
   modelv(a,:) = modelv(a,:) - mean(modelv(a,:));
   r(a) = var(modelv(a,:))/ var(lv(a,:));

   %all(a) =  cov(P(:,a)'*[L,  R]);
   %r(a) = all(a)/(lv(a) + rv(a));
end
[sr, ind] = sort(r, 'descend');
%P = P(:, ind);
%b = b(ind);
N = 3;
%M = ceil(size(P,2)/N);
M = 4
figure(1); clf;
indx = 1;
for PC_to_draw = ind
    leo = P(:,PC_to_draw)'*X;
    rem = P(:,PC_to_draw)'*modelX;
    leo = leo - mean(leo);
    rem = rem - mean(rem);
    
    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, '-');
    plot(rem, zeros(length(rem),1), 'ro', 'MarkerFaceColor', 'b', 'MarkerEdgeColor', 'b');
    plot(gxr, gyr, 'r-');
    title(sprintf('PC: %d, ratio: %f\n', PC_to_draw, sr(PC_to_draw)));
    grid on;
    %set(gca, 'Visible', 'off');
    axis tight;
end
return;