wolffd@0
|
1 function G = mk_2D_lattice_slow(nrows, ncols, wrap_around)
|
wolffd@0
|
2 % MK_2D_LATTICE Return adjacency matrix for 4-nearest neighbor connected 2D lattice
|
wolffd@0
|
3 % G = mk_2D_lattice(nrows, ncols, wrap_around)
|
wolffd@0
|
4 % G(k1, k2) = 1 iff k1=(i1,j1) is connected to k2=(i2,j2)
|
wolffd@0
|
5 %
|
wolffd@0
|
6 % If wrap_around = 1, we use toroidal boundary conditions (default = 0)
|
wolffd@0
|
7 %
|
wolffd@0
|
8 % Nodes are assumed numbered as in the following 3x3 lattice
|
wolffd@0
|
9 % 1 4 7
|
wolffd@0
|
10 % 2 5 8
|
wolffd@0
|
11 % 3 6 9
|
wolffd@0
|
12 %
|
wolffd@0
|
13 % e.g., G = mk_2D_lattice(3, 3, 0) returns
|
wolffd@0
|
14 % 0 1 0 1 0 0 0 0 0
|
wolffd@0
|
15 % 1 0 1 0 1 0 0 0 0
|
wolffd@0
|
16 % 0 1 0 0 0 1 0 0 0
|
wolffd@0
|
17 % 1 0 0 0 1 0 1 0 0
|
wolffd@0
|
18 % 0 1 0 1 0 1 0 1 0
|
wolffd@0
|
19 % 0 0 1 0 1 0 0 0 1
|
wolffd@0
|
20 % 0 0 0 1 0 0 0 1 0
|
wolffd@0
|
21 % 0 0 0 0 1 0 1 0 1
|
wolffd@0
|
22 % 0 0 0 0 0 1 0 1 0
|
wolffd@0
|
23 % so find(G(5,:)) = [2 4 6 8]
|
wolffd@0
|
24 % but find(G(1,:)) = [2 4]
|
wolffd@0
|
25 %
|
wolffd@0
|
26 % Using wrap around, G = mk_2D_lattice(3, 3, 1), we get
|
wolffd@0
|
27 % 0 1 1 1 0 0 1 0 0
|
wolffd@0
|
28 % 1 0 1 0 1 0 0 1 0
|
wolffd@0
|
29 % 1 1 0 0 0 1 0 0 1
|
wolffd@0
|
30 % 1 0 0 0 1 1 1 0 0
|
wolffd@0
|
31 % 0 1 0 1 0 1 0 1 0
|
wolffd@0
|
32 % 0 0 1 1 1 0 0 0 1
|
wolffd@0
|
33 % 1 0 0 1 0 0 0 1 1
|
wolffd@0
|
34 % 0 1 0 0 1 0 1 0 1
|
wolffd@0
|
35 % 0 0 1 0 0 1 1 1 0
|
wolffd@0
|
36 % so find(G(5,:)) = [2 4 6 8]
|
wolffd@0
|
37 % and find(G(1,:)) = [2 3 4 7]
|
wolffd@0
|
38
|
wolffd@0
|
39 if nargin < 3, wrap_around = 0; end
|
wolffd@0
|
40
|
wolffd@0
|
41 % M contains the number of each cell e.g.
|
wolffd@0
|
42 % 1 4 7
|
wolffd@0
|
43 % 2 5 8
|
wolffd@0
|
44 % 3 6 9
|
wolffd@0
|
45 % North neighbors (assuming wrap around) are
|
wolffd@0
|
46 % 3 6 9
|
wolffd@0
|
47 % 1 4 7
|
wolffd@0
|
48 % 2 5 8
|
wolffd@0
|
49 % Without wrap around, they are
|
wolffd@0
|
50 % 1 4 7
|
wolffd@0
|
51 % 1 4 7
|
wolffd@0
|
52 % 2 5 8
|
wolffd@0
|
53 % The first row is arbitrary, since pixels at the top have no north neighbor.
|
wolffd@0
|
54
|
wolffd@0
|
55 if nrows==1
|
wolffd@0
|
56 G = zeros(1, ncols);
|
wolffd@0
|
57 for i=1:ncols-1
|
wolffd@0
|
58 G(i,i+1) = 1;
|
wolffd@0
|
59 G(i+1,i) = 1;
|
wolffd@0
|
60 end
|
wolffd@0
|
61 if wrap_around
|
wolffd@0
|
62 G(1,ncols) = 1;
|
wolffd@0
|
63 G(ncols,1) = 1;
|
wolffd@0
|
64 end
|
wolffd@0
|
65 return;
|
wolffd@0
|
66 end
|
wolffd@0
|
67
|
wolffd@0
|
68
|
wolffd@0
|
69 npixels = nrows*ncols;
|
wolffd@0
|
70
|
wolffd@0
|
71 N = 1; E = 2; S = 3; W = 4;
|
wolffd@0
|
72 if wrap_around
|
wolffd@0
|
73 rows{N} = [nrows 1:nrows-1]; cols{N} = 1:ncols;
|
wolffd@0
|
74 rows{E} = 1:nrows; cols{E} = [2:ncols 1];
|
wolffd@0
|
75 rows{S} = [2:nrows 1]; cols{S} = 1:ncols;
|
wolffd@0
|
76 rows{W} = 1:nrows; cols{W} = [ncols 1:ncols-1];
|
wolffd@0
|
77 else
|
wolffd@0
|
78 rows{N} = [1 1:nrows-1]; cols{N} = 1:ncols;
|
wolffd@0
|
79 rows{E} = 1:nrows; cols{E} = [1 1:ncols-1];
|
wolffd@0
|
80 rows{S} = [2:nrows nrows]; cols{S} = 1:ncols;
|
wolffd@0
|
81 rows{W} = 1:nrows; cols{W} = [2:ncols ncols];
|
wolffd@0
|
82 end
|
wolffd@0
|
83
|
wolffd@0
|
84 M = reshape(1:npixels, [nrows ncols]);
|
wolffd@0
|
85 nbrs = cell(1, 4);
|
wolffd@0
|
86 for i=1:4
|
wolffd@0
|
87 nbrs{i} = M(rows{i}, cols{i});
|
wolffd@0
|
88 end
|
wolffd@0
|
89
|
wolffd@0
|
90
|
wolffd@0
|
91 G = zeros(npixels, npixels);
|
wolffd@0
|
92 if wrap_around
|
wolffd@0
|
93 for i=1:4
|
wolffd@0
|
94 if 0
|
wolffd@0
|
95 % naive
|
wolffd@0
|
96 for p=1:npixels
|
wolffd@0
|
97 G(p, nbrs{i}(p)) = 1;
|
wolffd@0
|
98 end
|
wolffd@0
|
99 else
|
wolffd@0
|
100 % vectorized
|
wolffd@0
|
101 ndx2 = sub2ind([npixels npixels], 1:npixels, nbrs{i}(:)');
|
wolffd@0
|
102 G(ndx2) = 1;
|
wolffd@0
|
103 end
|
wolffd@0
|
104 end
|
wolffd@0
|
105 else
|
wolffd@0
|
106 i = N;
|
wolffd@0
|
107 mask = ones(nrows, ncols);
|
wolffd@0
|
108 mask(1,:) = 0; % pixels in row 1 have no nbr to the north
|
wolffd@0
|
109 ndx = find(mask);
|
wolffd@0
|
110 ndx2 = sub2ind([npixels npixels], ndx, nbrs{i}(ndx));
|
wolffd@0
|
111 G(ndx2) = 1;
|
wolffd@0
|
112
|
wolffd@0
|
113 i = E;
|
wolffd@0
|
114 mask = ones(nrows, ncols);
|
wolffd@0
|
115 mask(:,ncols) = 0;
|
wolffd@0
|
116 ndx = find(mask);
|
wolffd@0
|
117 ndx2 = sub2ind([npixels npixels], ndx, nbrs{i}(ndx));
|
wolffd@0
|
118 G(ndx2) = 1;
|
wolffd@0
|
119
|
wolffd@0
|
120 i = S;
|
wolffd@0
|
121 mask = ones(nrows, ncols);
|
wolffd@0
|
122 mask(nrows,:)=0;
|
wolffd@0
|
123 ndx = find(mask);
|
wolffd@0
|
124 ndx2 = sub2ind([npixels npixels], ndx, nbrs{i}(ndx));
|
wolffd@0
|
125 G(ndx2) = 1;
|
wolffd@0
|
126
|
wolffd@0
|
127 i = W;
|
wolffd@0
|
128 mask = ones(nrows, ncols);
|
wolffd@0
|
129 mask(:,1)=0;
|
wolffd@0
|
130 ndx = find(mask);
|
wolffd@0
|
131 ndx2 = sub2ind([npixels npixels], ndx, nbrs{i}(ndx));
|
wolffd@0
|
132 G(ndx2) = 1;
|
wolffd@0
|
133 end
|
wolffd@0
|
134
|
wolffd@0
|
135 G = setdiag(G, 0);
|