%%%
%%% Lecture 7 - Spring 2000
%%% (code mostly derived from lecture8-99.m)
%%% (c) California Institute of Technology
%%% Pietro Perona, Apr. 10 1999 and Apr. 19 2000
%%%

clear;

%%
%% General control variables -- Play around with them to test
%% the software and the algorithm.
%%
SHOW_EXAMPLE =1;
SHAPE_VARIABILITY = 0.3; %% In display units
MAX_IMA_N = 20; %% Maximum number of false alarms per feature type
MIN_LAMBDA = 2; %% Minimum Poisson lambda for each feature type
MAX_LAMBDA = 4; %% Maximum Poisson lambda for each feature type


%%
%% Image size parameters; used by code that
%% generates the test images.
%%
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.LABELS = ['A'; 'B'; 'C'; 'D'; 'E']; %% Labels of the parts - used in displaying data
OBJ.P = [0.9 0.9 0.9 0.9 0.9]; %% Prob. of missing a feature
OBJ.mm = [-1 1 1 1 1 -1 -1 -1 0 0]';  %% Mean shape
N_PARAM = length(OBJ.mm);   %% Number of shape parameters
OBJ.NF = N_PARAM/2;         %% Number of parts (each has 2 coordinates)
Y_INDEX=2*[1:OBJ.NF]; X_INDEX=Y_INDEX-1; %% Indices for accessing parameters
OBJ.SS = SHAPE_VARIABILITY^2*eye(N_PARAM);  %% Variance matrix of shape

%%
%% Show examples in order to demonstrate the variabiliy
%% of the object class
%%
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(1);
   end;
end; %if SHOW_EXAMPLE


%% Generate an instance of the object - for the purpose
%% of testing the detection code.
OBJ.PARAM = randnND(OBJ.mm,OBJ.SS,1);
OBJ.X = OBJ.PARAM(X_INDEX); OBJ.Y = OBJ.PARAM(Y_INDEX);

%%
%% GENERATE BACKGROUND MODEL AND FEATURES
%%
BGR.X = zeros(OBJ.NF,MAX_IMA_N); BGR.Y=BGR.X;
%% For each feature type generate
%% randomly a model and then produce false alarms
for ft=1:OBJ.NF,
        %% Generate the Poisson Lambda at random within the bounds
	BGR.LAMBDA(ft) = MIN_LAMBDA+(MAX_LAMBDA-MIN_LAMBDA)*rand(1);
	%% Given lambda generage n. of false alarms with approx. of Poisson
	BGR.N(ft) = round(BGR.LAMBDA(ft) + sqrt(BGR.LAMBDA(ft))*randn(1));
	%% Clip n. of false alarms to avoid less than zero or more than array size
	BGR.N(ft) = min(max(BGR.N(ft),0),MAX_IMA_N-1);
	%% Generate the false alarms
	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,:) = OBJ.LABELS(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);
text(OBJ.X,OBJ.Y,OBJ.LABELS,'FontSize',24); hold off;
pause(1);


%%%
%%% CREATE IMAGE
%%% This code sets up the data structure describing the data that
%%% were collected from the 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(1);




%%%
%%% 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. - This code is somewhat
%% difficult to understand, but easy in concept.
%% 
%% The hyp 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.i is the list number of each hyp - not terribly important
%%
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; %% Initialize to 1 in case no part was detected i.e. d===0
   for ft=1:OBJ.NF,
     if d(ft) ~=0, %% If the feature corresponding to the part ft was detected
	%% Pick out the coordinates of the foreground feature of type ft
	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));
	pr(ft) = IMA.N(ft) / BGR.LAMBDA(ft);
      else
	%% If the feature was not detected then the coordinates are undefined
	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);
   %% This is the uniform prob. dens. of the bg pts in the current hyp
   %% divided by the uniform prob. density of the bg pts in the null hyp
   p_bg = IMG.SIZE^sum(d); %%% Old clunky-looking code: (OBJ.NF-sum(d==0)); 
   %% Prod. of N_i if d(i)=1 and 1 if d(i)=0 i.e. how many such hyp. are possible
   p_hyp = 1/prod(IMA.N.^d);
   p_d = prod(OBJ.P.^d)*prod((1-OBJ.P).^(1-d));
   HYP.R(h) = p_fg * p_bg * poisson_ratio * p_hyp * p_d;
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;