%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Compute@DClassifiers.m
%
%  Copyright (c) 2012 Homayoun Bagherinia 
%
% This program is free software: you can redistribute it and%or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http:%%www.gnu.org%licenses/>.
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



% This function generates 18 filter parameters in 2-D space as follows:
% C: cell array of training data (for each image)
% output: output mat file containing 2-D classifiers parameters
% bgpath: background image path
% type: type of slope computations
%       1: eigenvector in 4-D
%       2: eigenvectors in 2-D
% intercept: 0: slope intersects the origin
%            1: slope calculated from centered data
% outliers: 0: do not remove saturated points 
%           >0: remove saturated points (with value > outliers) 
% a: offset e.g. [0 0 0];
% c: gain  e.g. [0.9489 1.1409 1.2712]
% g: GAMMA e.g. [0.5315 0.6322 0.5872]
% mr: margin ratio to expand the classifiers boundaries
% display: 0: no plots of classifiers shown
%          1: plot all classifiers


function Compute2DClassifiers(C, ParamStructure, output)

global SAT_VALUE;
global INTERCEPT;
global OUTLIER;
global MR;
global GAMMA;
global GAMMA_a;
global GAMMA_c;
global GAMMA_g;

bgpath = ParamStructure.BackgroundImagesFolder;
suffix = ParamStructure.BackgroundImageSuffix;
INTERCEPT = ParamStructure.Intercept;
OUTLIER = (ParamStructure.SaturatedPointRemoval > 0);
SAT_VALUE = ParamStructure.SaturatedPointRemoval;
MR = ParamStructure.MarginRatio;
display = ParamStructure.ScatterplotsDisplay;
GAMMA = ParamStructure.GammaCorrection;

if isfield(ParamStructure,'a')
    GAMMA_a = ParamStructure.a;
else
    GAMMA_a = [0 0 0];
end
if isfield(ParamStructure,'c')
    GAMMA_c = ParamStructure.c;
else
    GAMMA_c = [1 1 1];
end
if isfield(ParamStructure,'g')
    GAMMA_g = ParamStructure.g;
else
    GAMMA_g = [1 1 1];
end

close all;
Ct = C;
CC = [];
for k=1:size(C,1)
    Ct{k,1} = Ct{k,1} + randn(size(Ct{k,1}));
    CCt = Ct{k,1};
    CC = [CC CCt];
end;

if GAMMA
    CC(1:3:10,:) = invgamma(CC(1:3:10,:), GAMMA_a(1), GAMMA_c(1), GAMMA_g(1));
    CC(2:3:11,:) = invgamma(CC(2:3:11,:), GAMMA_a(2), GAMMA_c(2), GAMMA_g(2));
    CC(3:3:12,:) = invgamma(CC(3:3:12,:), GAMMA_a(3), GAMMA_c(3), GAMMA_g(3));
end

[m b1 b2 t_min t_max] = mintotaldistance(CC);

% find the false positive rate (FPR) of each filter
[fpr2Dfilters filterM]= filtersfprminpatch18(bgpath, suffix, m, b1, b2);

% plot all filters
if ( display )
    plotclassifiers(Ct, filterM);
end

% sort filters in increasing order of FPR
[B IX] = sort(fpr2Dfilters);

% save the filter parameters
save(output, 'm', 'b1', 'b2', 't_min', 't_max', 'IX', 'fpr2Dfilters');


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This function computes all 18 maximum distance D_ij in 2-D
% and a single maximum distance in 12-D
function [m b1 b2  t_min t_max] = mintotaldistance(C)

global SAT_VALUE;
global INTERCEPT;
global OUTLIER;
global MR;
global GAMMA;
global GAMMA_a;
global GAMMA_c;
global GAMMA_g;


m = [];
Cr = C(1,:);
Cr = [Cr; C(4,:)];
Cr = [Cr; C(7,:)];
Cr = [Cr; C(10,:)];
for i=1:3
    for j=i+1:4
       Ci = Cr(i,:); Cj = Cr(j,:);
       if OUTLIER
            maxidx = find(Ci<SAT_VALUE); 
            Ci = Ci(maxidx); Cj = Cj(maxidx);
            maxidx = find(Cj<SAT_VALUE); 
            Ci = Ci(maxidx); Cj = Cj(maxidx);
        end
       if INTERCEPT
           uu_ij = Cj / [Ci;ones(size(Ci))]; resid_ij = norm(uu_ij * [Ci;ones(size(Ci))] - Cj);
           uu_ji = Ci / [Cj;ones(size(Cj))]; resid_ji = norm(uu_ji * [Cj;ones(size(Cj))] - Ci);
       else
           uu_ij = Cj / Ci; resid_ij = norm(uu_ij * Ci - Cj);
           uu_ji = Ci / Cj; resid_ji = norm(uu_ji * Cj - Ci);
       end
       if resid_ij < resid_ji
           m = [m; uu_ij(1)];
       else
           m = [m; 1/uu_ji(1)];
       end
    end
end

Cg = C(2,:);
Cg = [Cg; C(5,:)];
Cg = [Cg; C(8,:)];
Cg = [Cg; C(11,:)];
for i=1:3
    for j=i+1:4
       Ci = Cg(i,:); Cj = Cg(j,:);
       if OUTLIER
            maxidx = find(Ci<SAT_VALUE); 
            Ci = Ci(maxidx); Cj = Cj(maxidx);
            maxidx = find(Cj<SAT_VALUE); 
            Ci = Ci(maxidx); Cj = Cj(maxidx);
       end
       if INTERCEPT
           uu_ij = Cj / [Ci;ones(size(Ci))]; resid_ij = norm(uu_ij * [Ci;ones(size(Ci))] - Cj);
           uu_ji = Ci / [Cj;ones(size(Cj))]; resid_ji = norm(uu_ji * [Cj;ones(size(Cj))] - Ci);
       else
           uu_ij = Cj / Ci; resid_ij = norm(uu_ij * Ci - Cj);
           uu_ji = Ci / Cj; resid_ji = norm(uu_ji * Cj - Ci);
       end
       if resid_ij < resid_ji
           m = [m; uu_ij(1)];
       else
           m = [m; 1/uu_ji(1)];
       end
    end
end

Cb = C(3,:);
Cb = [Cb; C(6,:)];
Cb = [Cb; C(9,:)];
Cb = [Cb; C(12,:)];
for i=1:3
    for j=i+1:4
       Ci = Cb(i,:); Cj = Cb(j,:);
       if OUTLIER
            maxidx = find(Ci<SAT_VALUE); 
            Ci = Ci(maxidx); Cj = Cj(maxidx);
            maxidx = find(Cj<SAT_VALUE); 
            Ci = Ci(maxidx); Cj = Cj(maxidx);
        end
       if INTERCEPT
           uu_ij = Cj / [Ci;ones(size(Ci))]; resid_ij = norm(uu_ij * [Ci;ones(size(Ci))] - Cj);
           uu_ji = Ci / [Cj;ones(size(Cj))]; resid_ji = norm(uu_ji * [Cj;ones(size(Cj))] - Ci);
       else
           uu_ij = Cj / Ci; resid_ij = norm(uu_ij * Ci - Cj);
           uu_ji = Ci / Cj; resid_ji = norm(uu_ji * Cj - Ci);
       end
       if resid_ij < resid_ji
           m = [m; uu_ij(1)];
       else
           m = [m; 1/uu_ji(1)];
       end
    end
end

m1 = [m(1) m(2) m(3) m(7) m(8) m(9) m(13) m(14) m(15) m(4) m(5) m(10) m(11) m(16) m(17) m(6) m(12) m(18)];
m = m1';
[b1,b2,t_min,t_max] = mandumaxdistance2D(C, m);
     
     

%%%%%
function [b1 b2 t_min t_max] = mandumaxdistance2D(C, m)

global SAT_VALUE;
global INTERCEPT;
global OUTLIER;
global MR;
global GAMMA;
global GAMMA_a;
global GAMMA_c;
global GAMMA_g;

b1 = -inf*ones(size(m)); b2 = -b1;

siz = size(C);
dc = 0;

t_min=[]; t_max = [];

for j=1:12
    for k=j+3:3:12
        dc = dc+1;
        for i=1:siz(2)
            c = [C(j,i); C(k,i)]; 
            if (~OUTLIER) || all(c < SAT_VALUE)
                aux = c(2) - m(dc) * c(1);
                if aux > b1(dc)
                    b1(dc) = aux;
                end
                 if aux < b2(dc)
                    b2(dc) = aux;
                 end
            end
        end
        
        bb = (b1(dc) + b2(dc)) / 2;
        b1(dc) = bb + MR * (b1(dc) - bb);
        b2(dc) = bb - MR * (bb - b2(dc));
        
        aux_min(1) = max([0,min(C(j,:)) / MR]);
        aux_min(2) = max([0, min(C(k,:)) / MR]);
        t_min = [t_min; round(aux_min)];
        
        aux_max(1) = min([255, max(C(j,:)) * MR]);
        aux_max(2) = min([255, max(C(k,:)) * MR]);        
        t_max = [t_max; round(aux_max)];
       
    end
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% finds the false positive rate of all 18 filters
function [fprarray filterM]= filtersfprminpatch18(folder, suffix,m, b1, b2)

global GAMMA;
global GAMMA_a;
global GAMMA_c;
global GAMMA_g;

%imfile = dir(fullfile(folder, '*png'));
imfile = dir(fullfile(folder, strcat('*', suffix)));

totalfilter = length(m);

filterM = Filters18(m, b1, b2);

fprarray = zeros(1, totalfilter);

inminsize = size(imread(strcat(folder, imfile(1).name)));

probwidth = 12;
probwith2 = probwidth/2;
% the probes are taken from every 3 pixels for faster computation
allimprobes = round((inminsize(1)-probwidth-1)/3)*round((inminsize(2)-probwidth-1)/3)*size(imfile,1);

for l=1:size(imfile,1)
    imin = double(imread(strcat(folder, imfile(l).name)));             % load current image
    if GAMMA
        imin(:,:,1) = invgamma(imin(:,:,1), GAMMA_a(1), GAMMA_c(1), GAMMA_g(1));
        imin(:,:,2) = invgamma(imin(:,:,2), GAMMA_a(2), GAMMA_c(2), GAMMA_g(2));
        imin(:,:,3) = invgamma(imin(:,:,3), GAMMA_a(3), GAMMA_c(3), GAMMA_g(3));
    end
    
    for i=1:18
        pcounter = detecttarget_2D_1(imin, filterM(i,:), probwith2);
        fprarray(1,i) = fprarray(1,i) + pcounter;
    end;
end;

fprarray = fprarray ./ allimprobes;

%%%%%%%
function filterM = Filters18(m, b1, b2)

filterM = zeros(length(m),5);
filterM(:,1) = m;
filterM(:,2) = b1;
filterM(:,3) = b2;
filterM(1,4) = 1;   filterM(1,5) = 4;       % Green(r), White(r)
filterM(2,4) = 1;   filterM(2,5) = 7;       % Green(r), Black(r)
filterM(3,4) = 1;   filterM(3,5) = 10;      % Green(r), Red(r)
filterM(4,4) = 2;   filterM(4,5) = 5;       % Green(g), White(g)
filterM(5,4) = 2;   filterM(5,5) = 8;       % Green(g), Black(g)
filterM(6,4) = 2;   filterM(6,5) = 11;      % Green(g), Red(g)
filterM(7,4) = 3;  filterM(7,5) = 6;        % Green(b), White(b)
filterM(8,4) = 3;  filterM(8,5) = 9;        % Green(b), Black(b)
filterM(9,4) = 3;  filterM(9,5) = 12;       % Green(b), Red(b)
filterM(10,4) = 4;   filterM(10,5) = 7;     % White(r), Black(r)
filterM(11,4) = 4;   filterM(11,5) = 10;    % White(r), Red(r)
filterM(12,4) = 5;  filterM(12,5) = 8;      % White(g), Black(g)
filterM(13,4) = 5;  filterM(13,5) = 11;     % White(g), Red(g)
filterM(14,4) = 6;  filterM(14,5) = 9;      % White(b), Black(b)
filterM(15,4) = 6;  filterM(15,5) = 12;     % White(b), Red(b)
filterM(16,4) = 7;   filterM(16,5) = 10;    % Black(r), Red(r)
filterM(17,4) = 8;  filterM(17,5) = 11;     % Black(g), Red(g)
filterM(18,4) = 9;  filterM(18,5) = 12;     % Black(b), Red(b)

%%%%%%%%
function fpcounter = detecttarget_2D_1(imin, filter, width)
imsize = size(imin);
imR = imin(:,:,1);
imG = imin(:,:,2);
imB = imin(:,:,3);

fpcounter = 0;
C = zeros(12,1);
for i=width+1:3:imsize(1)-width-1
    for j=width+1:3:imsize(2)-width-1
        
        % the current probe
        C(1) = imR(i-width,j-width);    
        C(2) = imG(i-width,j-width);    
        C(3) = imB(i-width,j-width);    
        C(4) = imR(i-width,j+width);   
        C(5) = imG(i-width,j+width);   
        C(6) = imB(i-width,j+width);   
        C(7) = imR(i+width,j+width);    
        C(8) = imG(i+width,j+width);
        C(9) = imB(i+width,j+width);
        C(10) = imR(i+width,j-width);  
        C(11) = imG(i+width,j-width);   
        C(12) = imB(i+width,j-width);
        
        fpcounter = filter_2D(C, filter, fpcounter);
    end;
end;
   
%%%%%
function fpcounter = filter_2D(C, filter, fpcounter)
    diff = C(filter(1,5)) - filter(1,1)*C(filter(1,4));
    if ( diff > filter(1,2) || diff < -filter(1,3) )
        return;
    end;
    fpcounter = fpcounter + 1;

%%%%%
function I = invgamma(I, GAMMA_a, GAMMA_c, GAMMA_g)
GAMMA_a = GAMMA_a/255;
I = double(I)./255;
I = round((((I-GAMMA_a)./GAMMA_c).^(1/GAMMA_g))*255);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function plotclassifiers(C, filter)

global GAMMA;
global GAMMA_a;
global GAMMA_c;
global GAMMA_g;


color = ['+r'; '+g'; '+b'; '+c'; '+m'; '+y'; '+k'; ....
         'or'; 'og'; 'ob'; 'oc'; 'om'; 'oy'; 'ok'; ...
         '*r'; '*g'; '*b'; '*c'; '*m'; '*y'; '*k'; ...
         'xr'; 'xg'; 'xb'; 'xc'; 'xm'; 'xy'; 'xk'; ...
         '.r'; '.g'; '.b'; '.c'; '.m'; '.y'; '.k'; ...
         '^r'; '^g'; '^b'; '^c'; '^m'; '^y'; '^k'; ...
         'vr'; 'vg'; 'vb'; 'vc'; 'vm'; 'vy'; 'vk'; ...
         '>r'; '>g'; '>b'; '>c'; '>m'; '>y'; '>k'; ...
         '<r'; '<g'; '<b'; '<c'; '<m'; '<y'; '<k'; ...
         'sr'; 'sg'; 'sb'; 'sc'; 'sm'; 'sy'; 'sk'; ...
         'dr'; 'dg'; 'db'; 'dc'; 'dm'; 'dy'; 'dk'; ...
         'pr'; 'pg'; 'pb'; 'pc'; 'pm'; 'py'; 'pk'; ...
         'hr'; 'hg'; 'hb'; 'hc'; 'hm'; 'hy'; 'hk'; ...
         '+r'; '+g'; '+b'; '+c'; '+m'; '+y'; '+k'; ....
         'or'; 'og'; 'ob'; 'oc'; 'om'; 'oy'; 'ok'; ...
         '*r'; '*g'; '*b'; '*c'; '*m'; '*y'; '*k'; ...
         'xr'; 'xg'; 'xb'; 'xc'; 'xm'; 'xy'; 'xk'; ...
         '.r'; '.g'; '.b'; '.c'; '.m'; '.y'; '.k'; ...
         '^r'; '^g'; '^b'; '^c'; '^m'; '^y'; '^k'; ...
         'vr'; 'vg'; 'vb'; 'vc'; 'vm'; 'vy'; 'vk'; ...
         '>r'; '>g'; '>b'; '>c'; '>m'; '>y'; '>k'; ...
         '<r'; '<g'; '<b'; '<c'; '<m'; '<y'; '<k'; ...
         'sr'; 'sg'; 'sb'; 'sc'; 'sm'; 'sy'; 'sk'; ...
         'dr'; 'dg'; 'db'; 'dc'; 'dm'; 'dy'; 'dk'; ...
         'pr'; 'pg'; 'pb'; 'pc'; 'pm'; 'py'; 'pk'; ...
         'hr'; 'hg'; 'hb'; 'hc'; 'hm'; 'hy'; 'hk'];

for i=1:size(filter,1)
    x1 = [0:255];
    maxC1X = 0;
    maxC1Y = 0;
    figure;
    hold on;
    for k=1:size(C,1)
        C1 = C{k,1};
        if GAMMA
            C1(1:3:10,:) = invgamma(C1(1:3:10,:), GAMMA_a(1), GAMMA_c(1), GAMMA_g(1));
            C1(2:3:11,:) = invgamma(C1(2:3:11,:), GAMMA_a(2), GAMMA_c(2), GAMMA_g(2));
            C1(3:3:12,:) = invgamma(C1(3:3:12,:), GAMMA_a(3), GAMMA_c(3), GAMMA_g(3));
        end
        plot(C1(filter(i,4),:),C1(filter(i,5),:),color(k,:),'MarkerSize', 3);
        maxC1X = max(maxC1X, max(C1(filter(i,4))));
        maxC1Y = max(maxC1Y, max(C1(filter(i,5))));
    end
    if ( GAMMA_a ~= 0 )
        x1 = [0:maxC1X];
    end
    
    plot(x1, filter(i,1)*x1+filter(i,2), 'b', x1, filter(i,1)*x1+filter(i,3), 'b');
    hold off
    axis([0 255 0 255])
    if ( GAMMA_a ~= 0 )
        axis([0 maxC1X 0 maxC1Y])
    end
    axis equal;
    xlabel(chname(filter(i,4)));
    ylabel(chname(filter(i,5)));
    grid on;
%    hgsave(strcat(chname(filter(i,4)),'_',chname(filter(i,5)),'.fig'));
end;

%%%%%%
function name = chname(id)
switch id
    case 1 
        name='G_r'; return;
    case 2 
        name='G_g'; return;
    case 3 
        name='G_b'; return;
    case 4 
        name='W_r'; return;
    case 5 
        name='W_g'; return;
    case 6 
        name='W_b'; return;
    case 7 
        name='O_r'; return;
    case 8 
        name='O_g'; return;
    case 9 
        name='O_b'; return;
    case 10 
        name='B_r'; return;
    case 11 
        name='B_g'; return;
    case 12 
        name='B_b'; return;
    otherwise
        name=num2str(id); return;
end;