%%%
%%% Lecture 8 (code mostly derived from lecture7.m)
%%%
clear;

%% General control variables
SHOW_EXAMPLE =0;


%% Image parameters
X_MIN = -5; X_MAX = 5; X_LEN = X_MAX-X_MIN;
Y_MIN = -5; Y_MAX = 5; Y_LEN = Y_MAX-Y_MIN;
IMG.SIZE = X_LEN*Y_LEN;
IMG.BOUNDARIES = [X_MIN X_MAX Y_MIN Y_MAX];


%%
%% GENERATE OBJECT MODEL AND OBJECT INSTANCE
%%
%% 
%% It is a Gaussian with initial parameters:
OBJ.mm = [-1 1 1 1 0 -1]';
N_PARAM = length(OBJ.mm);
OBJ.NF = N_PARAM/2; Y_INDEX=2*[1:OBJ.NF]; X_INDEX=Y_INDEX-1;
OBJ.SS = 0.1*eye(N_PARAM);
OBJ.LABELS = ['A'; 'B'; 'C'];
OBJ.P = [0.1 0.1 0.1]; %% Prob. of missing a feature
OBJ.PARAM = randnND(OBJ.mm,OBJ.SS,1);
OBJ.X = OBJ.PARAM(X_INDEX); OBJ.Y = OBJ.PARAM(Y_INDEX);

%% Show examples
if SHOW_EXAMPLE>0,
   N_EXAMPLES = 10;
   for ex=1:N_EXAMPLES,
      OBJ.PARAM = randnND(OBJ.mm,OBJ.SS,1);
      OBJ.X = OBJ.PARAM(X_INDEX); OBJ.Y = OBJ.PARAM(Y_INDEX);
      figure(1); clf; 
      patch(OBJ.X,OBJ.Y,'w'); hold on
      text(OBJ.X,OBJ.Y,OBJ.LABELS); hold off; axis(IMG.BOUNDARIES);
      title('Example of object'); pause;
   end;
end; %if SHOW_EXAMPLE


%%
%% GENERATE BACKGROUND MODEL AND FEATURES
%%
MAX_IMA_N = 20;
BGR.X = zeros(OBJ.NF,MAX_IMA_N); BGR.Y=BGR.X;
for ft=1:OBJ.NF, %% For each feature type generate
   				 %% randomly a model and then produce false alarms
	BGR.LAMBDA(ft) = 2+5*rand(1);
	BGR.N(ft) = round(BGR.LAMBDA(ft) + sqrt(BGR.LAMBDA(ft))*randn(1)); %% Approx. of Poisson
   BGR.N(ft) = min(max(BGR.N(ft),0),MAX_IMA_N-1);
   BGR.X(ft,1:BGR.N(ft)) = X_LEN*rand(1,BGR.N(ft))+X_MIN;
	BGR.Y(ft,1:BGR.N(ft)) = Y_LEN*rand(1,BGR.N(ft))+Y_MIN;
   BGR.LABEL(ft,:) = num2str(ft);
end;


%% Display both object and BACKGROUND FEATURES
figure(2); clf; hold on;
for ft=1:OBJ.NF, text(BGR.X(ft,1:BGR.N(ft)),BGR.Y(ft,1:BGR.N(ft)),BGR.LABEL(ft)); end; axis(IMG.BOUNDARIES);
hold on; text(OBJ.X,OBJ.Y,OBJ.LABELS,'FontSize',12); hold off;
pause


%%%
%%% CREATE IMAGE
%%%

IMA.X = BGR.X;
IMA.Y = BGR.Y;
for ft=1:OBJ.NF,
   IMA.N(ft) = BGR.N(ft)+1; %% How many features of each type have been detected
   IMA.X(ft,IMA.N(ft)) = OBJ.X(ft);
   IMA.Y(ft,IMA.N(ft)) = OBJ.Y(ft);
   %%IMA.LABEL(ft,:) = num2str(ft);
end;
IMA.LABEL = OBJ.LABELS;
%% Display image
figure(3); clf; hold on;
for ft=1:OBJ.NF,
   for n=1:IMA.N(ft),
      text(IMA.X(ft,n),IMA.Y(ft,n),['.' IMA.LABEL(ft) '_{' num2str(n) '}']);
   end;
end; axis(IMG.BOUNDARIES);
hold off;
title('The image: input data to the algorithm');
pause




%%%
%%% RANK ALL THE N-plet HYPOTHESES
%%%
IMA.NZ = IMA.N + ones(1,OBJ.NF); %% Number of features including dummy ones

HYP.N = prod(IMA.NZ); %% Number of hypotheses, including the ones with zeros
HYP.i = [1:HYP.N]';

%%
%% Produce the indices for all the HYP.
%% These are stored in the matrix HYP.H where every row corresponds
%% to a hypothesis, and every column to a feature. So: the ith row
%% of HYP.H lists the F indices corresponding to the ith hypothesis.
%%
HYP.H = zeros(HYP.N,OBJ.NF); %% Allocate the space
hh = HYP.i-1; %% Auxhiliary vector for generating the cols of HYP.H
HYP.H(:,1) = floor(hh/prod(IMA.NZ(2:OBJ.NF)));
for ft=2:(OBJ.NF-1),
   HYP.H(:,ft) = mod( floor(hh/prod(IMA.NZ((ft+1):OBJ.NF))), IMA.NZ(ft));
end;
HYP.H(:,OBJ.NF) = mod(hh,IMA.NZ(OBJ.NF));

%%
%% Calculate score of each HYP.
%%
x_fg = zeros(OBJ.NF,1); y_fg = x_fg;
for h=HYP.i', %% Index through all hyp.
   hyp = HYP.H(h,:);    %%The vector of indices corresponding to the hypothesis
   d = (hyp ~= 0);	 %%d(i)=0 iff hyp(i)=0
   
   poisson_ratio = 1;
   for ft=1:OBJ.NF,
      if d(ft) ~=0,
      	x_fg(ft) = IMA.X(ft,hyp(ft)); y_fg(ft) = IMA.Y(ft,hyp(ft));
         pr(ft) = poisson(IMA.N(ft)-1, BGR.LAMBDA(ft))/poisson(IMA.N(ft), BGR.LAMBDA(ft));
      else %% if tau(ft)==0;
         x_fg(ft) = NaN; y_fg(ft) = NaN;
         pr(ft) = 1;
      end;
   end;
   poisson_ratio=prod(pr);
   
   param_fg = zeros(N_PARAM,1); param_fg(X_INDEX,1) = x_fg; param_fg(Y_INDEX,1) = y_fg;
   p_fg = gaussian(OBJ.mm,OBJ.SS,param_fg);
   p_bg = IMG.SIZE^(OBJ.NF-sum(d==0)); %% This is the uniform prob. dens. of the bg pts
   p_hyp = prod(OBJ.P.^(1-d))*prod((1-OBJ.P).^d)/prod(IMA.N.^d);	%% Prod. of N_i if d(i)=1 and 1 if d(i)=0
   HYP.R(h) = p_fg*p_bg*poisson_ratio*p_hyp;
end;

%% Explore and display the best hypotheses
[ps,hs] = sort(-HYP.R); ps = -ps;
figure(4); clf;
for hn=1:6,
   hyp = HYP.H(hs(hn),:); %% This is the hyp. being examined
   d = (hyp ~= 0);	 	%%d(i)=0 iff hyp(i)=0 - the non-zero-index parts of the hyp
   for ft=1:OBJ.NF,
      if d(ft) ~=0, x_fg(ft) = IMA.X(ft,hyp(ft)); y_fg(ft) = IMA.Y(ft,hyp(ft));
      else x_fg(ft) = NaN; y_fg(ft) = NaN;
      end;
   end;
   log_r = log10(HYP.R(hs(hn)));
   subplot(2,3,hn); 
   i_d = find(d);
   text(x_fg(i_d),y_fg(i_d),OBJ.LABELS(i_d));  axis(IMG.BOUNDARIES);
   title(['Hyp ' num2str(hs(hn)) '    log(p)=' num2str(log_r)]);
end;
pause;