wolffd@0: function G = mk_2D_lattice_slow(nrows, ncols, wrap_around) 
wolffd@0: % MK_2D_LATTICE Return adjacency matrix for 4-nearest neighbor connected 2D lattice
wolffd@0: % G = mk_2D_lattice(nrows, ncols, wrap_around)
wolffd@0: % G(k1, k2) = 1 iff k1=(i1,j1) is connected to k2=(i2,j2)
wolffd@0: %
wolffd@0: % If wrap_around = 1, we use toroidal boundary conditions (default = 0)
wolffd@0: %
wolffd@0: % Nodes are assumed numbered as in the following 3x3 lattice
wolffd@0: %   1 4 7
wolffd@0: %   2 5 8 
wolffd@0: %   3 6 9
wolffd@0: %
wolffd@0: % e.g., G = mk_2D_lattice(3, 3, 0) returns
wolffd@0: %   0 1 0 1 0 0 0 0 0 
wolffd@0: %   1 0 1 0 1 0 0 0 0 
wolffd@0: %   0 1 0 0 0 1 0 0 0 
wolffd@0: %   1 0 0 0 1 0 1 0 0 
wolffd@0: %   0 1 0 1 0 1 0 1 0 
wolffd@0: %   0 0 1 0 1 0 0 0 1 
wolffd@0: %   0 0 0 1 0 0 0 1 0 
wolffd@0: %   0 0 0 0 1 0 1 0 1 
wolffd@0: %   0 0 0 0 0 1 0 1 0 
wolffd@0: % so find(G(5,:)) = [2 4 6 8] 
wolffd@0: % but find(G(1,:)) = [2 4]
wolffd@0: %
wolffd@0: % Using wrap around, G = mk_2D_lattice(3, 3, 1), we get
wolffd@0: %   0 1 1 1 0 0 1 0 0 
wolffd@0: %   1 0 1 0 1 0 0 1 0 
wolffd@0: %   1 1 0 0 0 1 0 0 1 
wolffd@0: %   1 0 0 0 1 1 1 0 0 
wolffd@0: %   0 1 0 1 0 1 0 1 0 
wolffd@0: %   0 0 1 1 1 0 0 0 1 
wolffd@0: %   1 0 0 1 0 0 0 1 1 
wolffd@0: %   0 1 0 0 1 0 1 0 1 
wolffd@0: %   0 0 1 0 0 1 1 1 0 
wolffd@0: % so find(G(5,:)) = [2 4 6 8] 
wolffd@0: % and find(G(1,:)) = [2 3 4 7]
wolffd@0: 
wolffd@0: if nargin < 3, wrap_around = 0; end
wolffd@0: 
wolffd@0: % M contains the number of each cell e.g.
wolffd@0: %   1 4 7
wolffd@0: %   2 5 8 
wolffd@0: %   3 6 9
wolffd@0: % North neighbors (assuming wrap around) are
wolffd@0: %   3 6 9
wolffd@0: %   1 4 7
wolffd@0: %   2 5 8
wolffd@0: % Without wrap around, they are
wolffd@0: %   1 4 7
wolffd@0: %   1 4 7
wolffd@0: %   2 5 8
wolffd@0: % The first row is arbitrary, since pixels at the top have no north neighbor.
wolffd@0: 
wolffd@0: if nrows==1
wolffd@0:   G = zeros(1, ncols);
wolffd@0:   for i=1:ncols-1
wolffd@0:     G(i,i+1) = 1;
wolffd@0:     G(i+1,i) = 1;
wolffd@0:   end
wolffd@0:   if wrap_around
wolffd@0:     G(1,ncols) = 1;
wolffd@0:     G(ncols,1) = 1;
wolffd@0:   end
wolffd@0:   return;
wolffd@0: end
wolffd@0: 
wolffd@0:   
wolffd@0: npixels = nrows*ncols;
wolffd@0: 
wolffd@0: N = 1; E = 2; S = 3; W = 4;
wolffd@0: if wrap_around
wolffd@0:   rows{N} = [nrows 1:nrows-1]; cols{N} = 1:ncols;
wolffd@0:   rows{E} = 1:nrows; cols{E} = [2:ncols 1];
wolffd@0:   rows{S} = [2:nrows 1]; cols{S} = 1:ncols;
wolffd@0:   rows{W} = 1:nrows; cols{W} = [ncols 1:ncols-1];
wolffd@0: else
wolffd@0:   rows{N} = [1 1:nrows-1]; cols{N} = 1:ncols;
wolffd@0:   rows{E} = 1:nrows; cols{E} = [1 1:ncols-1];
wolffd@0:   rows{S} = [2:nrows nrows]; cols{S} = 1:ncols;
wolffd@0:   rows{W} = 1:nrows; cols{W} = [2:ncols ncols];
wolffd@0: end
wolffd@0: 
wolffd@0: M = reshape(1:npixels, [nrows ncols]);
wolffd@0: nbrs = cell(1, 4);
wolffd@0: for i=1:4
wolffd@0:   nbrs{i} = M(rows{i}, cols{i});
wolffd@0: end
wolffd@0: 
wolffd@0: 
wolffd@0: G = zeros(npixels, npixels);
wolffd@0: if wrap_around
wolffd@0:   for i=1:4
wolffd@0:     if 0
wolffd@0:       % naive
wolffd@0:       for p=1:npixels
wolffd@0: 	G(p, nbrs{i}(p)) = 1;
wolffd@0:       end
wolffd@0:     else
wolffd@0:       % vectorized
wolffd@0:       ndx2 = sub2ind([npixels npixels], 1:npixels, nbrs{i}(:)');
wolffd@0:       G(ndx2) = 1;
wolffd@0:     end
wolffd@0:   end
wolffd@0: else
wolffd@0:   i = N;
wolffd@0:   mask = ones(nrows, ncols);
wolffd@0:   mask(1,:) = 0; % pixels in row 1 have no nbr to the north
wolffd@0:   ndx = find(mask);
wolffd@0:   ndx2 = sub2ind([npixels npixels], ndx, nbrs{i}(ndx));
wolffd@0:   G(ndx2) = 1;
wolffd@0: 
wolffd@0:   i = E;
wolffd@0:   mask = ones(nrows, ncols);
wolffd@0:   mask(:,ncols) = 0;
wolffd@0:   ndx = find(mask);
wolffd@0:   ndx2 = sub2ind([npixels npixels], ndx, nbrs{i}(ndx));
wolffd@0:   G(ndx2) = 1;
wolffd@0: 
wolffd@0:   i = S;
wolffd@0:   mask = ones(nrows, ncols);
wolffd@0:   mask(nrows,:)=0;
wolffd@0:   ndx = find(mask);
wolffd@0:   ndx2 = sub2ind([npixels npixels], ndx, nbrs{i}(ndx));
wolffd@0:   G(ndx2) = 1;
wolffd@0:   
wolffd@0:   i = W;
wolffd@0:   mask = ones(nrows, ncols);
wolffd@0:   mask(:,1)=0;
wolffd@0:   ndx = find(mask);
wolffd@0:   ndx2 = sub2ind([npixels npixels], ndx, nbrs{i}(ndx));
wolffd@0:   G(ndx2) = 1;
wolffd@0: end
wolffd@0: 
wolffd@0: G = setdiag(G, 0);