
function ris = training_ljmNL3 (classes, classfilemetadata, classdata, indx, N_GAUSS, comp, scale, tipo)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% ris = training_ljmNL3 (cl, indx, N_GAUSS, comp, scale, tipo)
%
% Date: 29 Dec 2003
% Description: training using averaged local jet features with third non
%   linearity. For each pixel a set of invariants are extracted and 
%   each of these invariants is splitted in positive negative and absolute
%   value (with respect to the background mean). The averaged values over
%   all image pixel are the feature vector. After dimension reduction
%   obtained thanks to FLD a MoG is estimated for all the data of each
%   class.
% USE:  cl classes to consider, ex. cl=1:12;
%       indx indexes of images to consider (the
%           same for all classes), i.e. it=1; indx{1}=0:499; etc.
%       N_GAUSS nr. of gaussians in MoG
%       comp nr. of dimensions in the final feature space
%       scale nr. of scale in wich invariants are computed; it is assumed
%           scale=3;
%       tipo define covariance matrix structure in MoG, ex. tipo = 'full'
%       ris variable storing the training models
% Author: Marc'Aurelio Ranzato
%
% modified doliver@mbari.org - Dec 05, 2004
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

fprintf(1,'TRAINING \n')

numclasses = length(classes);

n_tot_samp = 0;
for cc = 1:numclasses
    n_samples(cc) = length(indx{cc});
    n_tot_samp = n_tot_samp + n_samples(cc);
end
    
storet = zeros(n_tot_samp,(scale + 1)*9*3);

% *modified - total iterations
MAX_N_ITER = 1000;  % see gmmem

lc = numclasses;
index2 = 0;

for cc = 1:lc   % Load data   
    c=classdata{1,cc};  m=classfilemetadata{1,cc};
    storet(index2 + 1: index2 + n_samples(cc),:) = c(indx{cc}+1,:);
    resolution(cc,1) = mean(cell2mat(m(indx{cc}+1,2)));
    resolution(cc,2) = std(cell2mat(m(indx{cc}+1,2)));
    index2 = index2 + n_samples( cc );
end
                

%%% F L D %%%
fprintf(1,'Computing FLD \n')
mediatot = mean(storet);
storet = storet - repmat(mediatot,size(storet,1),1);
FLD = fisher_l_d(storet',comp,n_samples);
storet = (FLD'*storet')';   % on the 1st row the projection along FLD of 1st image features

% stima MoG
index2 = 0;
for jj = 1:lc
    
    st = storet(index2 + 1: index2 + n_samples(jj),:);
    
    fprintf(1,'computing MoG-%s \n', classes{jj});
    randn('state',0);
    mix = gmm(size(st,2), N_GAUSS{jj}, tipo{jj});
    option = zeros(1,14);
    option(1) = -1;
    option(3) = 1;
    option(5) = 1;
    option(14) = MAX_N_ITER;
    rand('state',0);
    randn('state',0);
    mix = gmminit(mix, st, option);
    mix = gmmem(mix, st, option);
    mog(jj).mix = mix;
    index2 = index2 + n_samples(jj);
    
end
    
%%ris.state =  s;
ris.mog = mog;
ris.FLD = FLD;
ris.mediatot = mediatot;
ris.resolution = resolution;
