function [points,edges] = lattice(X,Y,connect)%Function [points,edges]=lattice(X,Y,connect) generates a 2D, %"4-connected" Cartesian lattice with dimensions X and Y. %%Inputs:    X/Y - The dimensions of the lattice (i.e., the number of %               nodes will equal X*Y)%           connect - Optional iInteger indicating the connectivity of %               the lattice%               0 - 4-connected (Default)%               1 - 8-connected%               K>1 - Radially connected with radius = K%%Outputs:   points - Nx2 lattice nodes%           edges - Mx2 lattice edges%%%6/5/03 - Leo Grady% Copyright (C) 2002, 2003 Leo Grady <lgrady@cns.bu.edu>%   Computer Vision and Computational Neuroscience Lab%   Department of Cognitive and Neural Systems%   Boston University%   Boston, MA  02215%% This program is free software; you can redistribute it and/or% modify it under the terms of the GNU General Public License% as published by the Free Software Foundation; either version 2% of the License, or (at your option) any later version.%% This program is distributed in the hope that it will be useful,% but WITHOUT ANY WARRANTY; without even the implied warranty of% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the% GNU General Public License for more details.%% You should have received a copy of the GNU General Public License% along with this program; if not, write to the Free Software% Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.%% Date - $Id: lattice.m,v 1.2 2003/08/21 17:29:29 lgrady Exp $%========================================================================%%Check for single pointif X*Y==1    points=[1;1];    edges=[];    returnend%Read connect argumentif nargin < 3    connect=0;endif connect < 2    %Generate points    rangex=0:(X-1);    rangey=0:(Y-1);    [x y]=meshgrid(rangey,rangex);    points=[x(:),y(:)];    N=X*Y;    %Connect points    edges=[[1:N]',[(1:N)+1]'];    edges=[[edges(:,1);[1:N]'],[edges(:,2);[1:N]'+X]];    if connect == 1 %8-connected        border=1:N;        border1=find(mod(border,X)-1)';        border2=find(mod(border,X))';        edges=[edges;[border1,border1+X-1;border2,border2+X+1]];    end        excluded=find((edges(:,1)>N)|(edges(:,1)<1)|(edges(:,2)>N)| ...        (edges(:,2)<1));    edges([excluded;[X:X:((Y-1)*X)]'],:)=[]; else %Radial connectivity    %Generate points    X=X+2*connect;    Y=Y+2*connect;    rangex=0:(X-1);    rangey=0:(Y-1);    [x y]=meshgrid(rangey,rangex);    points=[x(:),y(:)];    N=X*Y;    %Generate adjacency matrix    distances=(points(:,1)-floor(Y/2)).^2 + (points(:,2)-floor(X/2)).^2;    inside=find(distances<=(connect^2));    adjVec=sparse(1,inside,1,1,N);    midPoint=find(points(:,1)==floor(Y/2) & points(:,2)==floor(X/2));    %Adjust vector to be centered around diagonal    adjVec=adjVec([midPoint:N,1:(midPoint-1)]);    adjVec(1)=0; %Remove self-connection    W=circulant(adjVec);    %Remove phantom points    index=1:N;    tmp=1:X:N;    tmp=tmp'*ones(1,connect)+ones(length(tmp),1)*(0:(connect-1));    tmp2=tmp+X-connect;    index([tmp(:)',tmp2(:)',1:(X*connect),N+[-(X*connect-1):0]])=[];    edges=adjtoedges(W(index,index));    points=points(index,:);end