wolffd@0: function G = mk_2D_lattice(nrows, ncols, con)
wolffd@0: % MK_2D_LATTICE Return adjacency matrix for nearest neighbor connected 2D lattice
wolffd@0: % G = mk_2D_lattice(nrows, ncols, con)
wolffd@0: % G(k1, k2) = 1 iff k1=(i1,j1) is a neighbor of k2=(i2,j2) 
wolffd@0: % (Two pixels are neighbors if their Euclidean distance is less than r.)
wolffd@0: % Default connectivity = 4.
wolffd@0: %
wolffd@0: % WE ASSUME NO WRAP AROUND. 
wolffd@0: %
wolffd@0: % This is the neighborhood as a function of con:
wolffd@0: %
wolffd@0: %  con=4,r=1  con=8,r=sqrt(2)   con=12,r=2   con=24,r=sqrt(8)
wolffd@0: %  nn            2nd order        4th order
wolffd@0: %                                  x         x x x x x
wolffd@0: %    x          x x x            x x x       x x x x x
wolffd@0: %  x o x        x o x          x x o x x     x x o x x
wolffd@0: %    x          x x x            x x x       x x x x x
wolffd@0: %                                  x         x x x x x
wolffd@0: %
wolffd@0: % Examples:
wolffd@0: % Consider a 3x4 grid
wolffd@0: %  1  4  7  10
wolffd@0: %  2  5  8  11
wolffd@0: %  3  6  9  12
wolffd@0: %  
wolffd@0: % 4-connected:
wolffd@0: % G=mk_2D_lattice(3,4,4);
wolffd@0: % find(G(1,:)) = [2 4]
wolffd@0: % find(G(5,:)) = [2 4 6 8]
wolffd@0: %
wolffd@0: % 8-connected:
wolffd@0: % G=mk_2D_lattice(3,4,8);
wolffd@0: % find(G(1,:)) = [2 4 5]
wolffd@0: % find(G(5,:)) = [1 2 3 4 6 7 8 9]
wolffd@0: 
wolffd@0: % meshgrid trick due to Temu Gautama (temu@neuro.kuleuven.ac.be)
wolffd@0: 
wolffd@0: if nargin < 3, con = 4; end
wolffd@0: 
wolffd@0: switch con,
wolffd@0:  case 4, r = 1;
wolffd@0:  case 8, r = sqrt(2);
wolffd@0:  case 12, r = 2;
wolffd@0:  case 24, r = sqrt(8);
wolffd@0:  otherwise, error(['unrecognized connectivity ' num2str(con)])
wolffd@0: end
wolffd@0: 
wolffd@0: 
wolffd@0: npixels = nrows*ncols;
wolffd@0: 
wolffd@0: [x y]=meshgrid(1:ncols, 1:nrows);
wolffd@0: M = [x(:) y(:)];
wolffd@0: M1 = repmat(reshape(M',[1 2 npixels]),[npixels 1 1]);
wolffd@0: M2 = repmat(M,[1 1 npixels]);
wolffd@0: %D = squeeze(sum(abs(M1-M2),2)); % Manhattan distance
wolffd@0: M3 = M1-M2;
wolffd@0: D = sqrt(squeeze(M3(:,1,:)) .^2 + squeeze(M3(:,2,:)) .^2); % Euclidean distance
wolffd@0: G = reshape(D <= r,npixels,npixels);
wolffd@0: G = setdiag(G, 0);