```Date: Tue, 9 Nov 1999 15:45:00 -0600 Reply-To: "Blodgett, John G." Sender: "SAS(r) Discussion" From: "Blodgett, John G." 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/ ```

Back to: Top of message | Previous page | Main SAS-L page