function [ii jj] = sparse_adj_matrix(sz, r, p) % % Construct sparse adjacency matrix (provides ii and jj indices into the % matrix) % % Usage: % [ii jj] = sparse_adj_matrix(sz, r, p) % % inputs: % sz - grid size (determine the number of variables n=prod(sz), and the % geometry) % r - the radius around each point for which edges are formed % p - in what p-norm to measure the r-ball, can be 1,2 or 'inf' % % outputs % ii, jj - linear indices into adjacency matrix (for each pair (m,n) % there is also the pair (n,m)) % % How to construct the adjacency matrix? % >> A = sparse(ii, jj, ones(1,numel(ii)), prod(sz), prod(sz)); % % % Example: % >> [ii jj] = sparse_adj_matrix([10 20], 1, inf); % construct indices for 200x200 adjacency matrix for 8-connect graph over a % grid of 10x20 nodes. % To visualize the graph: % >> [r c]=ndgrid(1:10,1:20); % >> A = sparse(ii, jj, ones(1,numel(ii)), 200, 200);; % >> gplot(A, [r(:) c(:)]); % % % % Copyright (c) Bagon Shai % Department of Computer Science and Applied Mathmatics % Wiezmann Institute of Science % http://www.wisdom.weizmann.ac.il/ % % Permission is hereby granted, free of charge, to any person obtaining a copy % of this software and associated documentation files (the "Software"), to deal % in the Software without restriction, subject to the following conditions: % % The above copyright notice and this permission notice shall be included in % all copies or substantial portions of the Software. % % The Software is provided "as is", without warranty of any kind. % % Sep. 2010 % % number of variables n = prod(sz); % number of dimensions ndim = numel(sz); tovec = @(x) x(:); N=cell(ndim,1); I=cell(ndim,1); % construct the neighborhood fr=floor(r); for di=1:ndim N{di}=-fr:fr; I{di}=1:sz(di); end [N{1:ndim}]=ndgrid(N{:}); [I{1:ndim}]=ndgrid(I{:}); N = cellfun(tovec, N, 'UniformOutput',false); N=[N{:}]; I = cellfun(tovec, I, 'UniformOutput',false); I=[I{:}]; % compute N radius according to p switch lower(p) case {'1','l1',1} R = sum(abs(N),2); case {'2','l2',2} R = sum(N.*N,2); r=r*r; case {'inf',inf} R = max(abs(N),[],2); otherwise error('sparse_adj_matrix:norm_type','Unknown norm p (should be either 1,2 or inf'); end N = N(R<=r+eps,:); % "to"-index (not linear indices) ti = bsxfun(@plus, permute(I,[1 3 2]), permute(N, [3 1 2])); sel = all(ti >= 1, 3) & all( bsxfun(@le, ti, permute(sz, [1 3 2])), 3); csz = cumprod([1 sz(1:(ndim-1))]); jj = sum( bsxfun(@times, ti-1, permute(csz, [1 3 2])), 3)+1; % convert to linear indices ii = repmat( (1:n)', [1 size(jj,2)]); jj = jj(sel(:)); ii = ii(sel(:));