-
Notifications
You must be signed in to change notification settings - Fork 3
/
xxpdist.m
170 lines (160 loc) · 5.67 KB
/
xxpdist.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
function dxxp=xxpdist(MDX,NDXP,meth)
% dxxp=XXPDIST(MDX,NDXP,meth)
% dxxp=XXPDIST(MDX,[],meth)
%
% Calculates all pairwise distances between the entries of a
% M-dimensional vector of D-dimensional positions and a different
% N-dimensional vector of D-dimensional positions
%
% INPUT:
%
% MDX An MxD dimensional array with M vectors in D dimensions,
% for example, the M positions [x(:) y(:)] in D=2
% OR: an Mx1 dimensional vector [x(:)] of x-coordinates
% defining ONE dimension of a grid when also NDXP is thus defined
% NDXP An NxD dimensional array with N vectors in D dimensions,
% for example, the N different positions [xp(:) yp(:)] in D=2
% ---> If this is empty, the graph is compared to itself
% OR: an Nx1 dimensional vector [y(:)] of y-coordinates
% meth 1 using something I pieced together [default unless 3 possible]
% 2 my original way of doing things in D=2
% 3 Matlab's built-in way of doing things (needs STATS toolbox)
% 4 as found on http://tinyurl.com/xxpdist
%
% OUTPUT:
%
% dxxp An MxN array with the distances between every pair of points
%
% EXAMPLE:
%
% xxpdist('demo1') % Produces no errors if all goes well
%
%% Compute one gridded set of distances from the first corner point
% x=[0:11]*30; y=[0:9]*40;
% [MDX,NDXP]=meshgrid(x(:),y(:));
% dxxp1=xxpdist([MDX(:) NDXP(:)],[MDX(1) NDXP(1)]);
% imagesc(reshape(dxxp1,length(y),length(x))); axis xy
% dxxp7=xxpdist([MDX(:) NDXP(:)],[MDX(7) NDXP(7)]);
%
% SEE ALSO:
%
% KERNELC2D, LOCALIZATION2D, SSPDIST, SSP21SP
%
% Last modified by fjsimons-at-alum.mit.edu, 01/15/2014
% Not tested for complex input, but the hope is that it works
% For non-integers there might be some rounding error
% Initialize
defval('MDX',[])
defval('NDXP',[])
if ~isstr(MDX)
if prod(size(MDX))==length(MDX) && prod(size(NDXP))==length(NDXP)
% This means that both together define one D=2 grid of coordinates
% Obviously could easily build in D=3 case with ZZ but then again
[MDX,NDXP]=meshgrid(MDX(:),NDXP(:));
% Unwrap and pass through again
MDX=[MDX(:) NDXP(:)];
NDXP=[];
end
% If the STATS toolbox is available, use that one
if help('stats'); defval('meth',3); end
% Prepare to do it
defval('meth',1)
if ~isempty(NDXP)
% The numbers of dimensions, D, must be the same of course
if ~all(size(MDX,2)==size(NDXP,2))
error('The number of dimensions must be the same in both arrays')
end
end
t=tic;
% Do it
switch meth
case 1
% We expand the distance as the square root of the dot product
% Matlab's DOT function isn't cutting it for us
% Matlab's BSXFUN function lets us expand the singleton dimensions
if ~isempty(NDXP)
dxxp=sqrt(bsxfun(@plus,sum([MDX.*conj(MDX)],2),sum([NDXP.*conj(NDXP)],2)')...
-2*(MDX*conj(NDXP')));
else
dxxp=-2*(MDX*conj(MDX'));
dxxp=sqrt(dxxp-bsxfun(@plus,diag(dxxp),diag(dxxp)')/2);
end
case 2
if isempty(NDXP)
% Can't make this any shorter if it's a grid onto itself, so duplicate
NDXP=MDX;
end
if size(MDX,2)~=2 || size(NDXP,2)~=2
error(sprintf('Method %s only works in 2 dimensions',meth))
end
% Obviously could easily build in D=3 case with ZZ but then again
[XX,XXP]=meshgrid(NDXP(:,1),MDX(:,1));
[YY,YYP]=meshgrid(NDXP(:,2),MDX(:,2));
% Could do this with BSXFUN(@minus as in PDIST
dxxp=sqrt((XX-XXP).^2+(YY-YYP).^2);
case 3
if ~isempty(NDXP)
dxxp=pdist2(MDX,NDXP,'euclidean');
else
% The next line takes only about half the time since it does half the work
dxxp=pdist(MDX,'euclidean');
% Compare this with TRILOS and TRILOSI... this doubles the time again
dxxp=squareform(dxxp);
% The indexing is, apparently D(I,J)=DD((I-1)*(M-I/2)+J-I)
end
case 4
% Is the DOT function in Matlab somehow more efficient than case 1??
if ~isempty(NDXP)
dxxp=sqrt(bsxfun(@plus,dot(MDX,MDX,2),dot(NDXP,NDXP,2)')-2*(MDX*conj(NDXP')));
else
dxxp=sqrt(bsxfun(@plus,dot(MDX,MDX,2),dot(MDX,MDX,2)')-2*(MDX*conj(MDX')));
end
otherwise
error('Specify a valid method: 1, 2, 3 or 4')
end
% Report speed
disp(sprintf('%s took %f seconds (method %i)',...
upper(mfilename),toc(t),meth))
% Provide output
varns={dxxp};
varargout=varns(1:nargout);
elseif strcmp(MDX,'demo1')
x=[0:63]*30; y=[0:127]*40;
% One set of 2-dimensional coordinates yet to be made
dxxp11=xxpdist(x,y,1);
dxxp21=xxpdist(x,y,2);
dxxp31=xxpdist(x,y,3);
dxxp41=xxpdist(x,y,4);
difer(dxxp11-dxxp21,[],[],NaN)
difer(dxxp11-dxxp31,[],[],NaN)
difer(dxxp11-dxxp41,[],[],NaN)
difer(dxxp21-dxxp31,[],[],NaN)
difer(dxxp31-dxxp41,[],[],NaN)
[X,Y]=meshgrid(x+1,y+2);
xp=[0:31]*10; yp=[0:101]*20; [XP,YP]=meshgrid(xp,yp);
% Two sets of 2-dimensional coordinates fully defined
dxxp1=xxpdist([X(:) Y(:)],[XP(:) YP(:)],1);
dxxp2=xxpdist([X(:) Y(:)],[XP(:) YP(:)],2);
dxxp3=xxpdist([X(:) Y(:)],[XP(:) YP(:)],3);
dxxp4=xxpdist([X(:) Y(:)],[XP(:) YP(:)],4);
difer(dxxp1-dxxp2,[],[],NaN)
difer(dxxp1-dxxp3,[],[],NaN)
difer(dxxp1-dxxp4,[],[],NaN)
difer(dxxp2-dxxp3,[],[],NaN)
difer(dxxp3-dxxp4,[],[],NaN)
% One set of 2-dimensional coordinates already defined
dxxp1=xxpdist([X(:) Y(:)],[],1);
dxxp2=xxpdist([X(:) Y(:)],[],2);
dxxp3=xxpdist([X(:) Y(:)],[],3);
dxxp4=xxpdist([X(:) Y(:)],[],4);
difer(dxxp1-dxxp2,[],[],NaN)
difer(dxxp1-dxxp3,[],[],NaN)
difer(dxxp1-dxxp4,[],[],NaN)
difer(dxxp2-dxxp3,[],[],NaN)
difer(dxxp3-dxxp4,[],[],NaN)
% Check that the first set and the last set are identical
difer(dxxp1-dxxp11,[],[],NaN)
difer(dxxp2-dxxp21,[],[],NaN)
difer(dxxp3-dxxp31,[],[],NaN)
difer(dxxp4-dxxp41,[],[],NaN)
end