Function Cluster,x,y,nx,ny,NREGIONS=nregions,THRES=threshold,PLOT=iplot ; ; Cluster takes two 1D integer arrays with X and Y coodinates of ; pixels somehow selected in a 2D (Nx,Ny) array and identifies ; clusters of pixels. X and Y should be Y-sorted, that is X should run faster ; while Y should increase monotonously. If that is not the case, it could be ; achieved with the following commends: ; i=sort(Y*2L^15+X) & Y=Y(i) & X=X(i) ; A cluster is defined so that every pixel in a cluster is adjacent to at ; list one other pixel from the same cluster (their X and Y coordinates do ; not differ by more then 1) and all the pixels that are adjacent to any ; pixel in a cluster belong ; to the same cluster. ; The function returns a 1D integer array of the same size as X and Y ; that contains cluster number for each pixel. Non-cluster members ; are marked with 0. ; Optional parameters: ; nregions - (output) contain the number of identified clusters also ; counting non-cluster pixels (if any) as a separate cluster ; thres - (input, default: 1) consider clusters with "threshold" ; or smaller number of members to be non-clustered (index eq 0) ; plot - (input) if set color pixels as clusters are identified ; (a bit slower), otherwise print % of completion ; ; History: 17-July-2000 N.Piskunov wrote the version optimized for ; clusters oriented preferentially along rows or ; columns. ; 21-July-2000 N.Piskunov modified to handle arbitrary oriented ; clusters in the optimal way. It is slower than the ; original version by about 10% for spectral orders ; which are nearly horizontal/vertical. ; if(n_params() lt 4) then begin print,'Syntax: ind=Cluster(x,y,nx,ny[,NREGIONS=n[,THRES=t[,/PLOT]]]' print,'where: x and y are 1D arrays with column and row numbers of pixels' print,' to be analized for clustering,' print,' nx and ny are the dimensions of the original image,' print,' n is the number of detected clusters+1' print,' t is the threshold for the size of a cluster' print,' /PLOT displays clusters as they are detected' print,'Cluster returns an array of the same size as x and y. The value' print,'of index associates pixel with the cluster number. Pixels not' print,'associated with clusters are marked with index 0.' return,0L endif if(n_elements(x) ne n_elements(y)) then begin print,'Cluster: X and Y arrays should have the same size' return,0L endif if(max(x) ge nx) then begin print,'Cluster: X cannot exceed NX-1' return,0L endif if(max(y) ge ny) then begin print,'Cluster: Y cannot exceed NY-1' return,0L endif if(keyword_set(iplot)) then plot,x,y,xs=3,ys=3,psym=3 n=n_elements(x) ; Number of pixels in objects index=intarr(n) ; Created index array ; The initial X and Y are Y-sorted, which ; means that X runs faster while Y is ; increasing monotonously. ix_sort=sort( x*2L^15+y) ; Pointers to X-sorted arrays iy_sort=sort((y*2L^15+x)(ix_sort)) ; Pointers to Y-sorted array from X-sorted ; arrays. This means that: ; x(ix_sort(iy_sort)) eq x and ; y(ix_sort(iy_sort)) eq y ; The double index cunstructed above allows us to find the nearest pixels ; directly without using "where" which is costly in large arrays. max_branch=1 ; Number of found branches branches=[0L,0L] ; First/last pixel of each branch iper=0 for i=0L,n-1L do begin ; Loop through pixels if(not keyword_set(iplot) and i*10L/n gt iper) then begin iper=i*10L/n ; % complete print,iper*10,'% done' endif j=i ; Find pixels nearest to i in the same column j1=iy_sort(i)-1L j2=iy_sort(i)+1L if(j1 ge 0L) then j=[ix_sort(j1),j] if(j2 lt n-1L) then j=[j,ix_sort(j2)] if(i gt 0L) then begin j=[i-1L,j] ; Find pixels nearest to i in the previous column j1=iy_sort(i-1L)-1L j2=iy_sort(i-1L)+1L if(j1 ge 0L) then j=[ix_sort(j1),j] if(j2 lt n-1L) then j=[j,ix_sort(j2)] endif if(i lt n-1L) then begin j=[j,i+1L] ; Find pixels nearest to i in the next column j1=iy_sort(i+1L)-1L j2=iy_sort(i+1L)+1L if(j1 ge 0L) then j=[ix_sort(j1),j] if(j2 lt n-1L) then j=[j,ix_sort(j2)] endif j=j(sort(j)) ; Find immediate neighours. Where searches in 9 pixels or less. j=j(where(x(j) ge x(i)-1 and x(j) le x(i)+1 and \$ y(j) ge y(i)-1 and y(j) le y(i)+1, nj)) nmax=max(j) if(nj gt 1) then begin ; Check if this pixel has neighbours jj=where(index(j) gt 0, njj) ; Check for existing branches if(njj eq 0) then begin ; None found (only unmarked pixels) free_branch=where(branches(0,*) eq -1, nfree) if(nfree gt 0) then begin free_branch=free_branch(0) branches(*,free_branch)=minmax(j) index(j)=free_branch endif else begin branches=[[branches],[minmax(j)]] ; First/last pixel of the current branch index(j)=max_branch ; Start a new branch max_branch=max_branch+1 ; Increment branch number endelse endif else begin ; This is part of the existing branch curr_index=index(j(jj)) if(min(curr_index) eq max(curr_index)) then begin ind=curr_index(0) index(j)=ind branches(0,ind)=j(0) branches(1,ind) endif else begin curr_branch=min(curr_index) curr_indices=curr_index(sort(curr_index)) curr_indices=curr_indices(uniq(curr_indices)) ; Paint all relevant branches with the same index for jjj=1,n_elements(curr_indices)-1 do begin ind=curr_indices(jjj) ; Next non-zero index i1=branches(0,ind) ; First pixel of the ind branch i2=branches(1,ind) ; Last pixel of the ind branch ii=where(index(i1:i2) eq ind)+i1 ; Find all pixels in this branch index(ii)=curr_branch ; and merge branches if needed branches(0,curr_branch)=i1branches(1,curr_branch) endfor ; Remove the branches that were merged to a lower index (curr_branch) for jjj=1,n_elements(curr_indices)-1 do begin ind=curr_indices(jjj) branches(*,ind)=[-1,-1] endfor ; Index unmarked pixels and verify they do not change the pointers index(j)=curr_branch ; Mark immediate neighbours branches(0,curr_branch)=j(0) branches(1,curr_branch) endelse endelse if(keyword_set(iplot)) then oplot,x(j),y(j),psym=3,col=(index(i)+1) mod 128 endif endfor if(not keyword_set(iplot)) then print,100,'% done' if(keyword_set(iplot)) then plot,x,y,xs=3,ys=3,psym=3,/nodata ind=0 if(not keyword_set(threshold)) then threshold=1L for i=0,max_branch-1 do begin i1=branches(0,i) ; First pixel of the ind branch i2=branches(1,i) ; Last pixel of the ind branch if(i1 ge 0) then begin j=where(index(i1:i2) eq i, nj)+i1 if(nj gt threshold) then begin ind=ind+1 index(j)=ind if(keyword_set(iplot)) then oplot,x(j),y(j),psym=3,col=(ind mod 128) endif else if(nj gt 0) then begin index(j)=0 endif endif endfor nregions=ind+1 return,index end