function hullsort,xin,yin,xout=xout,yout=yout,xcen=xcen,ycen=ycen,\$ areap=areap,bycen=bycen,twait=twait,verbose=verbose, _extra=e ;+ ;function hullsort ; reorder a 2D list from outside in and return the indices of ; the reorderd list. At each stage, compute the convex hull ; and remove the "farthest" point from the hull centroid. ; ;syntax ; iout=hullsort(xin,yin,xout=xout,yout=yout,xcen=xcen,ycen=ycen,\$ ; areap=areap,/bycen,twait=twait,verbose=verbose) ; ;parameters ; xin [INPUT; required] 1D array of x-coordinate points ; yin [INPUT; required] 1D array of y-coordinate points ; * sizes must match ; * only useful if there are more than 4 points ; ;keywords ; xout [OUTPUT] reorderd XIN ; yout [OUTPUT] reorderd YIN ; xcen [OUTPUT] centroid x-coordinate computed at each step ; ycen [OUTPUT] centroid y-coordinate computed at each step ; areap [OUTPUT] area covered by the convex hull at each stage ; bycen [INPUT] if set, eliminates points from convex hull based ; on distance to center ; * default is to eliminate those farthest from all other ; points on the convex hull ; twait [INPUT] waiting time for replotting on screen [seconds] ; * only applies if VERBOSE>10 ; * default is 0.001 ; verbose [INPUT] controls chatter ; _extra [JUNK] here only to prevent crashing the program ; ;example ; .run hullsort ; ;subroutines ; QHULL (IDL built-in) ; KILROY ; AREAPOLY ; ;history ; vinay kashyap (MMXII.VII) ; added call to AREAPOLY() (VK; MMXII.X) ;- ; usage ok='ok' & np=n_params() & nx=n_elements(xin) & ny=n_elements(yin) if np lt 2 then ok='Insufficient parameters' else \$ if nx eq 0 then ok='XIN is undefined' else \$ if ny eq 0 then ok='YIN is undefined' else \$ if nx ne ny then ok='XIN and YIN are incompatible' else \$ if nx lt 4 then ok='too few points in the input' if ok ne 'ok' then begin print,'Usage: iout=hullsort(xin,yin,xout=xout,yout=yout,xcen=xcen,ycen=ycen,\$' print,' areap=areap,/bycen,twait=twait,verbose=verbose)' print,' reorder a 2D list from the outside in and return indices of reordered list' if np ne 0 then message,ok,/informational if nx gt 0 then return,lindgen(nx) else return,-1L endif ; keywords vv=0L & if keyword_set(verbose) then vv=long(verbose[0])>1L waitsec=0.001 & if keyword_set(twait) then waitsec=double(twait[0])>0 ; outputs iout=lonarr(nx) & jout=lindgen(nx) xx=xin[*] & yy=yin[*] & xout=xx & yout=yy xcen=fltarr(nx)+mean(xx,/double,/nan) & ycen=fltarr(ny)+mean(yy,/double,/nan) areap=fltarr(nx) for i=0L,nx-3L-1L do begin qhull,xx,yy,tr & jj=reform(tr[0,*]) & mm=n_elements(xx) ;avx=mean(xx[jj],/double,/nan) & avy=mean(yy[jj],/double,/nan) avx=mean(xx,/double,/nan) & avy=mean(yy,/double,/nan) dd=sqrt((xx[jj]-avx)^2+(yy[jj]-avy)^2) xcen[i:*]=avx & ycen[i:*]=avy areap[i:*]=areapoly(xx[jj],yy[jj]) if keyword_set(bycen) then begin jnk=max(dd,imx) kk=jj[imx] & iout[i]=jout[kk] & xout[i]=xx[kk] & yout[i]=yy[kk] endif else begin ndd=n_elements(dd) & dd2=dblarr(ndd,ndd) for j=0L,ndd-1L do dd2[j,*]=sqrt((xx[jj]-xx[jj[j]])^2+(yy[jj]-yy[jj[j]])^2) dd22=total(dd2,2)/double(ndd) ;dd2j=dblarr(ndd) & for j=0L,ndd-1L do begin & oj=where(dd2[j,*] gt 0) & dd2j[j]=min(dd2[j,oj]) & endfor ;dd2j=dblarr(ndd) & for j=0L,ndd-1L do begin & os=sort(dd2[j,*]) & dd2j[j]=mean(dd2[j,os[0:ndd/3]],/double) & endfor dd2j=dblarr(ndd) & for j=0L,ndd-1L do begin & os=sort(dd2[j,*]) & dd2j[j]=mean(dd2[j,os[0:1]],/double) & endfor ;for j=0L,ndd-1L do dd[j]=dd[j]+dd22[j]+dd2j[j] for j=0L,ndd-1L do dd[j]=dd2j[j] ;for j=0L,ndd-1L do dd[j]=dd[j]+mean(sqrt((xx[jj]-xx[jj[j]])^2+(yy[jj]-yy[jj[j]])^2),/double) ;for j=0L,ndd-1L do dd[j]=dd[j]+mean(sqrt((xx-xx[jj[j]])^2+(yy-yy[jj[j]])^2),/double) jnk=max(dd,imx) kk=jj[imx] & iout[i]=jout[kk] & xout[i]=xx[kk] & yout[i]=yy[kk] endelse if vv gt 10 then begin plot,xx,yy,psym=1,/xs,/ys,title=i,/nodata oplot,xx,yy,psym=1,col=100 if vv gt 50 then oplot,xx[jj],yy[jj],psym=4,col=200 if vv gt 100 then oplot,[xout[i]],[yout[i]],psym=4,col=150,symsize=2 wait,waitsec endif ;ii=lindgen(mm) & oo=where(ii ne jj[imx]) jout[kk]=-1 & o0=where(jout ge 0) & jout=jout[o0] xx=xx[o0] & yy=yy[o0] if vv gt 0 and vv lt 100 then \$ if i*vv eq 1000*long((i*vv)/1000) then kilroy endfor xout[i:*]=xx & yout[i:*]=yy & iout[i:*]=jout[o0] if vv gt 1000 then stop,xout[i],yout[i] return,iout end ; example run if not keyword_set(npt) then npt=500L if not keyword_set(verbose) then verbose=101 if not keyword_set(twait) then twait=0.01 ;xin=randomn(seed,npt)*3+2 & yin=randomn(seed,npt)*6+4 xin=randomn(seed,npt/2)*3+2 & yin=randomn(seed,npt/2)*6+2 xin=[xin,randomn(seed,npt-npt/2)*1+4] yin=[yin,randomn(seed,npt-npt/2)*2+4] iout=hullsort(xin,yin,xout=xout,yout=yout,xcen=xcen,ycen=ycen,verbose=verbose,twait=twait) icol=byte((lindgen(npt)/float(npt))*192.)+63 if !d.name eq 'X' then window,0 & plot,xin,yin,xtitle='X',ytitle='Y',/xs,/ys,psym=3,title='all points and centroid evolution' for i=1L,npt-1L do oplot,[xcen[i-1L],xcen[i]],[ycen[i-1L],ycen[i]],col=icol[i] if !d.name eq 'X' then window,2 & plot,xin,yin,/nodata,xtitle='X',ytitle='Y',/xs,/ys,psym=3,title='points color coded by EE' for i=0L,npt-1L do plots,xout[i],yout[i],col=icol[i],psym=1 ; calling sequence print,'' jnk=hullsort() print,'' end