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