-
Notifications
You must be signed in to change notification settings - Fork 20
/
dists.m
192 lines (183 loc) · 6.68 KB
/
dists.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
function D = dists(X1,X2,p,e)
%DISTS Calculates distances between vectors of points.
% D = dists(X1,X2,p,e)
% X1 = n x d matrix of n d-dimensional points
% X2 = m x d matrix of m d-dimensional points
% D = n x m matrix of distances
% = (n-1) x 1 vector of distances between X1 points, if X2 = []
% p = 2, Euclidean (default): D(i,j) = sqrt(sum((X1(i,:) - X2(j,:))^2))
% = 1, rectilinear: D(i,j) = sum(abs(X1(i,:) - X2(j,:))
% = Inf, Chebychev dist: D(i,j) = max(abs(X1(i,:) - X2(j,:))
% = (1 Inf), lp norm: D(i,j) = sum(abs(X1(i,:) - X2(j,:))^p)^(1/p)
% = 'rad', great circle distance in radians of a sphere
% (where X1 and X2 are decimal degree latitudes and longitudes)
% = 'mi' or 'sm', great circle distance in statute miles on the earth
% = 'km', great circle distance in kilometers on the earth
% e = epsilon for hyperboloid approximation gradient estimation
% = 0 (default); no error checking if any non-empty 'e' input
% ~= 0 => general lp used for rect., Cheb., and p outside [1,2]
%
% Great circle distances are calculated using the Haversine Formula (from R.W.
% Sinnott, "Virtues of the Haversine", Sky and Telescope, vol. 68, no. 2, 1984
% p. 159, as reported in "http://www.census.gov/cgi-bin/geo/gisfaq?Q5.1")
% Copyright (c) 1998 by Michael G. Kay
% Matlog Version 1.0 Apr-3-98
% Input Error Checking *******************************************************
if nargin == 4 & ~isempty(e) % No error checking is 'e' input
[n,d] = size(X1);
m = size(X2,1);
else % Do error checking
error(nargchk(1,4,nargin));
e = 0; % 'e' default
[n,d] = size(X1);
if n == 0 | ~isreal(X1)
error('X1 must be non-empty real matrix');
end
if nargin < 2 | isempty(X2) % Calc intra-seq dist b/w X1 points
m = 0; % X2 default
if n < 2
error('X1 must have more than 1 point');
end
else % Calc dist b/w X1 and X2 points
[m,dX2] = size(X2);
if m == 0 | ~isreal(X2)
error('X2 must be non-empty real matrix');
end
if d ~= dX2
error('Rows of X1 and X2 must have same dimensions');
end
end
if nargin < 3 | isempty(p)
p = 2; % 'p' default
elseif ~ischar(p) % lp distances
if length(p(:)) ~= 1 | ~isreal(p)
error('''p'' must be a real scalar number');
end
elseif ischar(p) % Great circle distances
p = lower(p);
if d ~= 2
error('Points must be 2-dimensional for great-circle distances');
end
if ~any(strcmp(p,{'rad','mi','sm','km'}))
error('''p'' must be either ''rad,'' ''mi,'' ''sm,'' or ''km''');
end
else
error('''p'' not valid value');
end
end
% End (Input Error Checking) ***********************************************
% Interchange if X2 is the only 1 point
intrchg = 0;
if n > 1 & m == 1
tmp = X2; X2 = X1; X1 = tmp;
m = n;n = 1;
intrchg = 1;
end
% 1-dimensional points
if d == 1
if e == 0
if m ~= 0
D = abs(X1(:,ones(1,m)) - X2(:,ones(1,n))');
else
D = abs(X1(1:n-1) - X1(2:n))'; % X1 intra-seq. dist.
end
else
if m ~= 0
D = sqrt((X1(:,ones(1,m)) - X2(:,ones(1,n))').^2 + e);
else
D = sqrt((X1(1:n-1) - X1(2:n)).^2 + e)';
end
end
% X1 only 1 point or intra-seq dist
elseif n == 1 | m == 0
if n == 1 % Expand X1 to match X2
X1 = X1(ones(1,m),:);
n = m;
else % X1 intra-seq. dist.
X2 = X1(2:n,:); % X2 = ending points
n = n - 1;
X1 = X1(1:n,:); % X1 = beginning points
end
if p == 2 % Euclidean distance
D = sqrt(sum(((X1 - X2).^2 + e)'));
elseif ischar(p) % Great-circle distance
X1 = pi*X1/180;X2 = pi*X2/180;
D = 2*asin(min(1,sqrt(sin((X1(:,1) - X2(:,1))/2).^2 + ...
cos(X1(:,1)).*cos(X2(:,1)).* ...
sin((X1(:,2) - X2(:,2))/2).^2)))';
elseif p == 1 & e == 0 % Rectilinear distance
D = sum(abs(X1 - X2)');
elseif (p >= 1 & p <= 2) | (e ~= 0 & p > 0) % General lp distance
D = sum((((X1 - X2).^2 + e).^(p/2))').^(1/p);
elseif p == Inf & e == 0 % Chebychev distance
D = max(abs(X1 - X2)');
else % Otherwise
D = zeros(1,n);
for j = 1:n
D(j) = norm(X1(j,:) - X2(j,:),p);
end
end
% X1 and X2 are 2-dimensional points
elseif d == 2
if p == 2 % Euclidean distance
D = sqrt((X1(:, ones(1,m)) - X2(:, ones(1,n))').^2 + e + ...
(X1(:,2*ones(1,m)) - X2(:,2*ones(1,n))').^2 + e);
elseif ischar(p) % Great-circle distance
X1 = pi*X1/180;X2 = pi*X2/180;
cosX1lat = cos(X1(:,1));cosX2lat = cos(X2(:,1));
D = 2*asin(min(1,sqrt(...
sin((X1(:,ones(1,m)) - X2(:,ones(1,n))')/2).^2 + ...
cosX1lat(:,ones(1,m),1).*cosX2lat(:,ones(1,n))'.* ...
sin((X1(:,2*ones(1,m)) - X2(:,2*ones(1,n))')/2).^2)));
elseif p == 1 & e == 0 % Rectilinear distance
D = abs(X1(:,ones(1,m)) - X2(:,ones(1,n))') + ...
abs(X1(:,2*ones(1,m)) - X2(:,2*ones(1,n))');
elseif (p >= 1 & p <= 2) | (e ~= 0 & p > 0) % General lp distance
D = (((X1(:, ones(1,m)) - X2(:, ones(1,n))').^2 + e).^(p/2) + ...
((X1(:,2*ones(1,m)) - X2(:,2*ones(1,n))').^2 + e).^(p/2)).^(1/p);
elseif p == Inf & e == 0 % Chebychev distance
D = max(abs(X1(:,ones(1,m)) - X2(:,ones(1,n))'),...
abs(X1(:,2*ones(1,m)) - X2(:,2*ones(1,n))'));
else % Otherwise
D = zeros(n,m);
for i = 1:n
for j = 1:m
D(i,j) = norm(X1(i,:) - X2(j,:),p);
end
end
end
% X1 and X2 are 3 or more dim. point
else
if p == 2 % Euclidean distance
D = sqrt(sum((repmat(reshape(X1,[n 1 k]),1,m) - ...
repmat(reshape(X2,[1 m k]),n,1)).^2 + e,3));
elseif p == 1 & e == 0 % Rectilinear distance
D = sum(abs(repmat(reshape(X1,[n 1 k]),1,m) - ...
repmat(reshape(X2,[1 m k]),n,1)),3);
elseif (p >= 1 & p <= 2) | (e ~= 0 & p > 0) % General lp distance
D = sum(((repmat(reshape(X1,[n 1 k]),1,m) - ...
repmat(reshape(X2,[1 m k]),n,1)).^2 + e).^(p/2),3).^(1/p);
elseif p == Inf & e == 0 % Chebychev distance
D = max(abs(repmat(reshape(X1,[n 1 k]),1,m) - ...
repmat(reshape(X2,[1 m k]),n,1)),[],3);
else % Otherwise
D = zeros(n,m);
for i = 1:n
for j = 1:m
D(i,j) = norm(X1(i,:) - X2(j,:),p);
end
end
end
end
% Transpose D if X2 was interchanged or X1 intra-seq distances
if intrchg == 1 | m == 0
D = D.';
end
% Convert 'rad' to 'km' or 'mi' (or 'sm')
if ischar(p) & ~strcmp(p,'rad')
if strcmp(p,'km')
D = (6378.388 - 21.476*abs(sin(mean([X1(:,1);X2(:,1)]))))*D;
else % 'mi' or 'sm'
D = (3963.34 - 13.35*abs(sin(mean([X1(:,1);X2(:,1)]))))*D;
end
end