function [sep_measure,peak_order,transformation,decision_threshold] = compute_fisher(data,label,in_peak_order)

dimPC = 3;
number = size (data, 1);
dimgene = size(data, 2);
dimsel = dimgene;
dimPeak = length(in_peak_order);

%normgene = mean(abs(data), 1);

%[temp, normorder] = sort(normgene);
%idgene = normorder(dimgene - dimsel + 1 : dimgene);
%newdata = data(:, idgene);

newdata = data;
newdata = data;
id1 = (label == 1);
id2 = (label == 2);
data1 = newdata(id1, :);
data2 = newdata(id2, :);

Peakdata1 = data1(:, in_peak_order);
Peakdata2 = data2(:, in_peak_order);
Sw = (cov(Peakdata1) + cov(Peakdata2)) / number;
Sb = (mean(Peakdata1, 1) - mean(Peakdata2, 1))' * (mean(Peakdata1, 1) - mean(Peakdata2, 1));
[V, D] = eig(Sb, Sw);
[temp, PCorder] = sort(abs(diag(D)));
PCdata1 = Peakdata1 * V(:, PCorder(dimPeak - dimPC + 1 : dimPeak));
PCdata2 = Peakdata2 * V(:, PCorder(dimPeak - dimPC + 1 : dimPeak));

PCdata1=Peakdata1;
PCdata2=Peakdata2;

IOdisplay(number, PCdata1, PCdata2, label);

class1_mean = sum(PCdata1(:,3))/size(PCdata1,1);
class2_mean = sum(PCdata2(:,3))/size(PCdata2,1);
class1_std = std(PCdata1(:,3));
class2_std = std(PCdata2(:,3));
%decision_threshold.val = (class1_mean*class2_std+class2_mean*class1_std)/(class1_std+class2_std);
%decision_threshold.val = class1_mean*.5 +class2_mean*.5;
decision_threshold.val = class1_mean*(class1_std/(class1_std + class2_std)) +class2_mean*(class2_std/(class1_std + class2_std));
if(class1_mean > decision_threshold.val) decision_threshold.label =1;
else decision_threshold.label=2;
end

sep_measure = (mean(Peakdata1, 1) - mean(Peakdata2, 1)) * inv(Sw) * (mean(Peakdata1, 1) - mean(Peakdata2, 1))';
peak_order = in_peak_order;
transformation = V(:, PCorder(dimPeak - dimPC + 1 : dimPeak));