function Voxels = Lda(L, Y, Opts)
% Voxels = two most discriminant voxels
% L      = 1 x NumTrials condition labels vector
% Y      = NumVoxels x NumTrials data matrix

if nargin < 3, Opts = []; end

if ~isfield(Opts, 'Plot'),
    Opts.Plot = 1;
end;
if ~isfield(Opts, 'Legend'),
    Opts.Legend = [];
end;
if ~isfield(Opts, 'MarkerSize'),
    Opts.MarkerSize = 7;
end;
if ~isfield(Opts, 'FontSize'),
    Opts.FontSize = 8;
end;
if ~isfield(Opts, 'Voxels'),
    Opts.Voxels = [];
end;
if ~isfield(Opts, 'DistClasses'),
    Opts.DistClasses = [1 2];
end;

Z = []; 
NumClasses = size(unique(L),2);

%% Do it
if ~isempty(Opts.Voxels)
    Voxels = Opts.Voxels;
else
    cls = fitcdiscr(Y', L');
    % pick two most discriminant voxels btw classes A and B
    W = cls.Coeffs(Opts.DistClasses(1),Opts.DistClasses(2)).Linear;
    % sort voxels according to information content
    [d.W_sorted, W_order] = sort(W);
    
    % do LDA again with two of the most informative voxels only
    Voxels = [W_order(1) W_order(end)];
end

Yv = Y(Voxels,:);
Z = fitcdiscr(Yv', L');
% indeces for the Coeff matrix
[CoeffRows, CoeffCols] = find(tril(ones(NumClasses),-1));



%% plot
if Opts.Plot
    Axes = [min(Yv(1,:)) max(Yv(1,:)) min(Yv(2,:)) max(Yv(2,:))];
    Cols   = {'k', 'r', 'b', 'g'};
    Shapes = {'o', 'o', 'o', 'o'};
    SCols   = cell2mat(Cols(1:NumClasses));
    SShapes = cell2mat(Shapes(1:NumClasses));
    % border labels
    Borders = arrayfun(@(c,r) [num2str(r) '-' num2str(c)], CoeffRows, CoeffCols, 'UniformOutput', false);
    
    %figure(2)
    Sc = gscatter(Yv(1,:), Yv(2,:), L, SCols, SShapes, Opts.MarkerSize, 'off');
    hold on
    for c = 1:NumClasses
        Sc(c).MarkerFaceColor = Cols{c};
        
        K = Z.Coeffs(CoeffRows(c),CoeffCols(c)).Const; % First retrieve the coefficients for the linear
        B = Z.Coeffs(CoeffRows(c),CoeffCols(c)).Linear;% boundary between the 1st and 2nd classes
        
        % Vox_1ot the curve K + [x,y]*B  = 0.
        f = @(x1,x2) K + B(1)*x1 + B(2)*x2;
        h = ezplot(f, Axes);
        h.Color = Cols{c};
        h.LineWidth = 1;
    end
    
    if ~isempty(Opts.Legend)
        Legend = cellfun(@(l,b) l, Opts.Legend', Borders, 'UniformOutput', false);
        %legend(Legend, 'Location', 'northwest');
        legend(Legend, 'Location', 'best');
        %legend('boxoff');
    end
    
    set(gca, 'FontSize', Opts.FontSize);
    axis(Axes)
    title([]);
    xlabel([]);
    ylabel([]);
end
end

