function [inds,posmul]=ssp21sp(s,siz,fl) % [inds,posmul]=SSP21SP(s,siz,fl) % % Expresses the distance between a linearly indexed grid node and all the % other grid nodes as indices into the vector that contains the distances % between the first linear index of the grid and all the other nodes % % INPUT: % % s The linear index into a two-dimensional array % siz The size of the array (rows, columns) % fl 0 Runs a non-redundant full loop with some extra overhead [default] % 1 Runs a redundant full loop % % 2 Runs a short loop for all positions PAST the given one % % OUTPUT: % % inds The indices that map Dij(1,inds)=Dij(s,:) in the distance % matrix Dij that contains all pairwise distances between the % elements s and sp in the grid % posmul The relevant index positions and their multiplicity % % SEE ALSO: % % SSPDIST, XXPDIST % % Last modified by fjsimons-at-alum.mit.edu, 01/15/2014 defval('fl',1) % Keep time t=tic; % The size of the array M=siz(1); N=siz(2); MN=M*N; % A handful of special cases if s==1 inds=1:MN; posmul=[inds ; ones(size(inds))]; return end % Allocate index array inds=zeros(1,MN); %-[fl==3]*[s+1]); if s==MN inds=MN:-1:1; posmul=[inds ; ones(size(inds))]; return end % Only one remaining case after contraction of loops below is when max/min % don't commute which is when there is a NaN from 0/0, on the last grid row, % i.e. when a=-1 and d=0 it fails to fill out to the left since we picked % the max-min order to favor the right side... so prepare for that if rem(s-1,M)+1==M inds(s-M:-M:2)=[1:round([s-M-1]/M)]*M+1; end if s>1 && s0 && 1==3 disp(sprintf('%i out of %i are repeated hits, time taken %6.4fs',... redfax,MN,toc(t))) end end % Check that every element has been hit as it should diferm(sum(~~inds),MN) % If second output is requested if nargout==2 pos=unique(inds); posmul=hist(inds,pos); % And the reconstruction must be happening diferm(gamini(pos,posmul),sort(inds)) % Return both the position and its multiplicity posmul=[pos; posmul]; end