%% Geo-502 Lab 
% Matlab's distance related functions
% Phaedon Kyriakidis

%% Matlab functions for distance computation
% pdist  -> pairwise distances between 1 set of points
% pdist2 -> pairwise distances between 2 sets of points (R2010a)
% squareform -> transform pdist output, of length (N*N-N)/2, to square matrix

%% Create some scattered point locations (random coordinates)
rng(1234) % 1234 is the seed
x1 = rand(10,1)*100; y1 = rand(10,1)*100;
x2 = rand(20,1)*100; y2 = rand(20,1)*100;

%% Display point locations
plot(x1,y1,'*')
hold on
plot(x2,y2,'ro')
hold off

%% Add text with point id's
id1 = num2str((1:10)'); % convert numbers to character strings
id2 = num2str((1:20)');
h=text(x1,y1,id1);
set(h,'color','b'); 
h=text(x2,y2,id2);
set(h,'color','r');

%% Create grid node coordinates
hold on
xG = linspace(5,95,10); % put 10 points equally spaced from 5 to 95
% same as xG = 5:10:95
yG = linspace(5,95,10);
[X,Y] = ndgrid(xG,yG); 
%[X,Y] = medhgrid(xG,yG);
size(X)
xG = X(:); yG = Y(:); % make them long column vectors
plot(xG,yG,'k.');
axis([0 100 0 100]); % display limits [xmin xmax ymin ymax]
idG = num2str((1:100)'); % text with id's
text(xG,yG,idG);
hold off

%% Matrix of pairwise distances between points of a single set
D = squareform(pdist([x1 y1]));
size(D) % dimensions of array D
size(pdist([x1 y1])) % pdist outputs upper triangular part in a vector form
(10*10-10)/2
% Display all entries of matrix D as an image
imagesc(D); colorbar
% Symmetric matrix with zeros along the diagonal:
[D(1,2) D(2,1) D(4,5) D(5,4)]
[D(1,1) D(2,2) D(10,10)]
%%
% Row (or column) #5 of matrix |D|, for example, contains the distances
% between event location with coordinates (x(5),y(5)) and all other 
% event locations

%% Pairwise distance between 2 sets of points
D12 = pdist2([x1 y1],[x2 y2]);
size(D12)
imagesc(D12); colorbar
%axis equal
%axis tight

%% DEFINING A CENTRAL POINT (ONE OF EXISTING EVENTS)
% The objective here is to define a "central" existing location in the 
% point pattern, or that event (one from the N in the point pattern) for
% which the total (Euclidean) distance from all other events is minimum.

%%
% We are looking for that event location, or for that row of matrix |D|,
% whose element sum (total event-to-all-other-events distance) is minimum.
% The next step would be therefore to compute the row-sum of matrix |D|:
sumDrow = sum(D,2); % 2 indicates the 2nd dimension (columns) of a matrix
size(sumDrow) % size of output
hist(sumDrow,30); % histogram of total event-to-all-other-event distances
xlabel('distance'); ylabel('frequency');
title('Histogram of total distances')
%%
% Last, we are looking for the index of the smallest element in array 
% |sumDrow|, which corresponds to the minimum event-to-all-other-events
% distance:
[minTotDist,indexMin] = min(sumDrow)
%%
% where |minTotDist| is the smallest element (distance) of array |sumDrow|,
% and |indexMin| is the index of that element. In other words, |indexMin|
% is the integer we are looking for:
x1Centr = x1(indexMin); y1Centr = y1(indexMin);
plot(x1,y1,'b*',x1Centr,y1Centr,'ro');

%% ARBITRARY CENTRAL LOCATION (NOT ONE OF EXISTING EVENTS)
% First, compute the distances beteween all candidate points (those on the
% 10x10 grid) and the scattered events, and then take the minimum

%% Distances between grid nodes and scattered points
DG1 = pdist2([xG yG],[x1 y1]);
size(DG1)
DG2 = pdist2([xG yG],[x2 y2]);
size(DG2)

%% Central point of same RED set searching on a 10x10 grid
sumDG2 = sum(DG2,2); % row sum of DG2 
[dMin2,idxMin2] = min(sumDG2);
plot(xG,yG,'k.');
axis([0 100 0 100]); % display limits [xmin xmax ymin ymax]
hold on
plot(x2,y2,'ro')
plot(xG(idxMin2),yG(idxMin2),'g*')
hold off

%% Search for arbitrary central point (not necessarily one of existing points)
xGzoom = linspace(45,65,10); % put 10 points equally spaced from 5 to 95
yGzoom = linspace(45,65,10);
[Xzoom,Yzoom] = ndgrid(xGzoom,yGzoom); 
xGzoom = Xzoom(:); yGzoom = Yzoom(:); % make them long column vectors
plot(xGzoom,yGzoom,'k.');
axis([0 100 0 100]); % display limits [xmin xmax ymin ymax]
idG = num2str((1:100)'); % text with id's
text(xGzoom,yGzoom,idG);