Date: Tue, 9 Nov 1999 15:45:00 -0600
Reply-To: "Blodgett, John G." <BlodgettJ@UMSYSTEM.EDU>
Sender: "SAS(r) Discussion" <SAS-L@LISTSERV.UGA.EDU>
From: "Blodgett, John G." <BlodgettJ@UMSYSTEM.EDU>
Subject: Re: find the max # of contiguous points in a plane
Content-Type: text/plain; charset="iso-8859-1"
-----Original Message-----
From: Fang Wan [mailto:fwan@UHMC.SUNYSB.EDU]
Sent: Monday, November 08, 1999 12:58 PM
To: SAS-L@LISTSERV.UGA.EDU
Subject: find the max # of contiguous points in a plane
I have a sas dataset with 66 variables xiyj (-5<= i <=5, 0<= j <=5)
which denote values at (i,j) in a plane, xiyj take values either 0 or 1,
my question is how to find the max number of points which have values 1
and are contiguous each other in the plane.
Note that two points are contiguous if they are adjacent by side or by
diagonal.
for example (1,1 ) and (1,2) are contiguous points and (1,2) and (2,3)
are contiguous points.
Any suggestions and ideas will be appreciated!
fang wan
suny at stony brook
========================================
I am attaching below a SAS program that I believe solves this problem. It
may need to be generalized a bit for different data but it was tested and
worked on the randomly-generated test data. If the code is too hard to
read here it can be viewed in its orginal format at
http://www.oseda.missouri.edu/cgi-bin/browse?/mscdc//sascode/contig_points.s
as .
title 'Contiguous Points in Space Problem';
*---this setup stored in d:\sas8\mysas\contig_points.sas --;
libname user 'd:\temp\contig_points_user';
*--create random array of 0/1 values using ranuni--;
data points;
do i=-5 to 5;
do j=0 to 5;
val=mod ( round(ranuni(0)*10) , 2); output;
end;
end;
keep val;
stop;
run;
*--Step to create the "first order neighbors" of each point in the array--*;
data nbrs1(keep=k k2) ref(keep=k nnbrs1 nbr1-nbr8 nbrstr);
set points;
length default=4 nbrstr $66; *--nbrstr is "neighbor string". The ith
character of nbrstr is either
blank or a hex digit, x, indicating that point i is an xth order
neighbor of current point k. For example,
a 2nd-order neighbor is a point that is not directly contiguous to pt
k but is directly contiguous to a
point that is directly contiguous (i.e. a 1st order neighbor). 3rd
order neighbors are direct
neighbors of 2nd order neighbors, etc.--;
k+1; *--key of the "current" point. There will be 66 values of k-;
array avals(66) aval1-aval66; retain aval1-aval66;
avals(k)=val;
if k < 66 then return;
*--all the point values are stored. Now create and output the 1st order
neighbors data-;
k=0;
do i=-5 to 5;
array nbrs(8) nbr1-nbr8; *--each point can have up to 8 immediate
neighbors---;
do j=0 to 5;
k+1;
val=avals(k);
nnbrs=0; nbrstr=' ';
if val then do; *--no need to check anything if the value at the
point is 0--;
substr(nbrstr,k,1)='1'; *--flag pt as contiguous to itself;
do i2=max(i-1,-5) to min(i,5); *--search row below, this row, row
above--;
k2=6*(i2+5); *--this will be the 1-66 index for the neighbor
point. We number from
the bottom row up. So row -5 contains pts 1-6,
row -4 contains pts
7-12, etc.---;
do j2=max(j-1,0) to min(j,5); *--search col left, this col, col
right--;
k2 + j2 +1; *--this becomes the pt number of the immediate
neighbor pt--;
if (i ne i2 or j ne j2) then if avals(k2) then do; *--if
directly contig pt has val of 1;
nnbrs+1; nbrs(nnbrs)=k2; *--count and store neighboring
pt number in array-;
substr(nbrstr,k2,1)='1'; *--set flag within 66-char string
to indicate this point is 1st order neighbor-;
output nbrs1;
end;
end; *--j2 loop--;
end; *--i2 loop--;
end; *--if val do group-;
do ix= nnbrs+1 to 8; nbrs(ix)=0; end; *--not really necessary but
looks better--;
nnbrs1=nnbrs; *--number of 1st order neighbors-;
output ref; *--initial version. Get updated in each cycle thru
&lvl loop in macro, below--;
end; *--j loop--;
end; *--i loop--;
stop;
run;
proc print data=nbrs1; *--have a look to make sure all looks reasonable--;
title 'nbrs1: 1st order neigbors';
run;
proc print data=ref; title 'Initial ref set'; run;
*--Idea now is to set up an iterative process where we look at ith-order
neighbors and create (i+1)th order ones.
We continue to do this until we detect no new neighbors at all during a
cycle, at which point we are done.
We get by with less than 10 cycles so we can use 1-character values to
indicate contiguity order. We really
should allow for that order number to go higher (by creating alpha codes
to represent values of 10+)--;
%macro dolvls;
%*----Just a big macro do loop to invoke the iterative
process--------------------------------*;
%do lvl=2 %to 9; %*--for some data we might need to check for more than 9
contiguity levels--;
%let plvl=%eval(&lvl -1); %*--"p" is for *p*revious level--;
%let lvlc=&lvl; %let plvlc=&plvl; %*--add code for hex digits if we
need over 9 levels--;
*---We have a data set with 66 observations, 1 per point in our 11 x 6
array. Each obs has the
index numbers of all points "contiguous" to that point. So if obs 24
has nnbrs=2, nbr1=23, nbr2=25
it means that point 24 is immediately contiguous with pts 23 and 25.
We can then go to obs 23 and
25 to get their immediate neighbors and do a "chain" to get all the
contig points.----*;
** options symbolgen;
data nbrs&lvl(keep=k k3 rename=(k3=k2))
newref(keep=k nnbrs1 nbr1-nbr8 nbrstr totnbrs lvl);
%*--the 2 sets created are:
1. the new neighbors set defining all &lvl-th level neighbors,
2. an updated ref set with info from previous cycle, which then gets
renamed
and reused in next cycle as input ref. Two new vars added that
were not in the
initial ref: totnbrs and lvl (total # of contiguous points and
the numeric contiguity level).
--*;
retain lvl "&lvl";
if _n_=1 then link loadref; *--load current ref set into temporary
arrays--;
array _nbrstr (66) $66 _temporary_;
array _nnbrs1(66) _temporary_;
array _nbrs1 (66,8) _temporary_;
array _n1 (8) nbr1-nbr8;
do until(last); *--use our own loop to read this set -- a tad bit more
efficient.--;
set nbrs&plvl end=last; *--gets variables containing the keys for the
current fixed point, k, and a
neighbor of that point at level &plvl, k2 --;
******* put k= k2= last=;
*--we have the 1st order neighbors, if any, of point k2 now in _nbr1
thru _nbr8--;
do _i_=1 to _nnbrs1(k2);
k3=_nbrs1(k2,_i_); *--retrieve _i_th first order neighbor of point
k2--;
*--check to see if this neighbor is already a lower order neighbor.
Add it at the
current level only if it has not already been assigned at a
lower level--;
if substr( _nbrstr(k),k3,1)=' ' then do;
output nbrs&lvl;
substr(_nbrstr(k),k3,1)="&lvlc"; *--insert code indicating this
is a current level nbr of
point k--;
_newnbrs+1; *--count new neighbors found during the entire
run--*;
end;
end;
end; *--loop to read nbrs&plvl;
*---Finished finding new nbrs but now we must rewrite the ref file
from the updated _nbrstr array---;
do _i_=1 to 66;
set ref;
nbrstr=_nbrstr(_i_);
totnbrs=length(compress(nbrstr)); *--total neighbor count based on
# non-blank chars in the neighbor string;
output newref;
end;
file log; put "***Processing for level &lvl neighbors completed. "
_newnbrs " new neighbors detected.";
call symput('newnbrs',left(put(_newnbrs,5.)));
stop;
loadref:
do _i_=1 to 66;
_obs=_i_;
set ref point=_obs;
_nbrstr(_i_)=nbrstr; _nnbrs1(_i_)=nnbrs1;
do _j_=1 to nnbrs1;
_nbrs1(_i_,_j_)=_n1(_j_);
end;
end;
file log;
put "****Ref data loaded for step &lvl ****";
return;
run;
proc print data=newref;
title2 "Listing of ref set after processing level &lvl";
run;
proc datasets library=user nolist; *--get ready for next cycle for renaming
newref to ref and ref to oldref;
age newref ref oldref ;
quit;
%if &newnbrs=0 %then %do;
%put *****NO more new neighbors found -- look for max value totnbrs on
the ref data set*****;
%goto aftrloop; %*--there should a %leave stmt for this--;
%end;
%end; %*-- the macro loop over lvl--;
%aftrloop:
%mend dolvls;
options mprint;
%dolvls
title2 'Look for the Answer In This Report';
proc means data=ref min max mean;
var totnbrs;
run;
John Blodgett
OSEDA (Office of Social & Economic Data Analysis)
626 Clark Hall -- University of Missouri
Columbia, MO 65211
PH: 573-884-2727 FX: 573-884-4635
e-mail: blodgettj@umsystem.edu
URL: http://www.oseda.missouri.edu/jgb/
| 
