%
% Example on epage XXX
% 
% Problem on epage 517
% 
% Written by:
% -- 
% John L. Weatherwax                2008-02-20
% 
% email: wax@alum.mit.edu
% 
% Please send comments and especially bug reports to the
% above email address.
% 
%-----

clear all; 
close all; 
clc; 

addpath('../../Code/CSTool');

randn('seed',0); 
rand('seed',0); 

%
% First generate a large number of points according to a homogeneous Poisson process:
% 
% by making the variable "r" larger we can simulate more points ... 
% 

n_desired = 100; % <- the number of desired points
delta     = 0.1; % <- our inhibition width

% Set the lambda. 
lambda = 2;
%r = 5;
r = 10;
tol = 0;
i=1;
% Generate the exponential random variables.
while tol < pi*r^2
  x(i) = exprnd(1/lambda,1,1);
  tol = sum(x);
  i=i+1;
end
x(end)=[];
N = length(x);
% Get the coordinates for the angles.
th = 2*pi*rand(1,N);
R = zeros(1,N);
% Find the R_i.
for i = 1:N
  R(i) = sqrt(sum(x(1:i))/pi);
end
[Xc,Yc]=pol2cart(th,R);

fh=figure; plot( Xc, Yc, 'xr', 'markersize', 5 ); hold on; 
xlabel('X(i)'); ylabel('Y(i)'); 
title( 'the original homogeneous Poisson points' ); 

%
% Now "thin" these points [Xc,Yc] if they are "too" close
% 
inds      = randperm(N);              % <- permute [Xc,Yc] randomly 
Xc        = Xc(inds);   
Yc        = Yc(inds); 

dist      = pdist([Xc.',Yc.']);
DM        = squareform(dist); 
% these are the (ii,jj) pairs of the points (Xc(ii),Yc(ii)) and (Xc(jj),Yc(jj)) that are too close 
% ... one of these points must be deleted ... 
bad_inds  = find(DM(:)<=delta); 
[Ib,Jb]   = ind2sub(size(DM),bad_inds);
dum       = find( Ib==Jb ); 
Ib(dum)   = []; Jb(dum) = [];         % <- remove diagonal terms ... 
%
% this is a very computationally poor method for doing this ... 
% but a poor algorithm is better than none ... 
%
while( length(Ib)>0 )
  % remove one bad point: 
  Xc(Ib(1))=[]; Yc(Ib(1))=[]; 
  dist      = pdist([Xc.',Yc.']);
  DM        = squareform(dist); 
  % these are the (ii,jj) pairs of the points (Xc(ii),Yc(ii)) and (Xc(jj),Yc(jj)) that are too close 
  % ... one of these points must be deleted ... 
  bad_inds  = find(DM(:)<=delta); 
  [Ib,Jb]   = ind2sub(size(DM),bad_inds);
  dum       = find( Ib==Jb ); 
  Ib(dum)   = []; Jb(dum) = []; 
end

if( length(Xc) < n_desired )
  error('Not enough points left ... increase r and rerun'); 
end

% take a random selection of these remaining good points: 
r_perm = randperm(length(Xc));
r_inds = r_perm(1:n_desired);
Xc_thinned = Xc(r_inds); Yc_thinned = Yc(r_inds); 

figure(fh); plot( Xc_thinned, Yc_thinned, 'og', 'markersize', 8, 'markerfacecolor', 'g' ); 
title( 'original points with selected thinned points' ); 
saveas( gcf, '../../WriteUp/Graphics/Chapter12/prob_12_6_thinned', 'epsc' );

clear functions;