OK, not much interest in this... I came up with a solution that works for me. However, I'm not sure if this is the right way of doing that. And it might need some tweaking...
x and y are coordinates, footprint is the desired spatial bin size.
;generate distance table
n = n_elements(x)
d=sqrt( (rebin(transpose(x),n,n,/SAMPLE)-rebin(x,n,n,/SAMPLE))^2 + $
(rebin(transpose(y),n,n,/SAMPLE)-rebin(y,n,n,/SAMPLE))^2 )
iclose=INTARR(n,n)
igroup=LONARR(n)-1
g=0L
hh=HISTOGRAM(d, REVERSE_INDICES=ri, BINSIZE=round(footprint), min=0, NBINS=2)
for i=0, N_ELEMENTS(hh)-1 do begin
if i eq 0 then k=1 else k=0 ; 1 is for points within footprint
ii=ri[ri[i]:ri[i+1]-1]
iclose[ii]=k
endfor
repeat begin
; find remaining points
hh=HISTOGRAM(igroup, REVERSE_INDICES=ri, BINSIZE=1, min=-1, NBINS=1)
remaining=ri[ri[0]:ri[0+1]-1] ;points that are not grouped yet
;choose new startpoint
nclose=total(iclose,1,/INTEGER)
tt=min(nclose[remaining],damin) ;next start point is the one with least close neighbors
;now check to which points this one is close
hh=HISTOGRAM(iclose[remaining[damin],*], REVERSE_INDICES=ri, BINSIZE=1, min=0, NBINS=2)
ii=remaining[ri[ri[1]:ri[1+1]-1]] ;index to close neighbors
;and check which of these neighbors has max number of close neighbors and small distane to them
tt=total(iclose[ii, *], 2) * 1d/(1+sqrt(VARIANCE(x[iclose[ii]])+VARIANCE(y[iclose[ii]])))
tt=max(tt,damax)
;ii[damax] is new center point
;store indices for this bin
hh=HISTOGRAM(iclose[ii[damax],*], REVERSE_INDICES=ri, BINSIZE=1, min=0, NBINS=2)
ii=ri[ri[1]:ri[1+1]-1]
igroup[ii]=g & g++
;remove indices from the iclose (symmetric) array
iclose[ii,*]=0 & iclose[*,ii]=0
endrep until total(iclose) eq 0 ;stop when no close points are left
;do cleanup
;...
;... ?
;make bins
hh=HISTOGRAM(igroup, REVERSE_INDICES=ri)
nr=N_ELEMENTS(hh)
xr=FLTARR(nr) & yr=FLTARR(nr) & zr=FLTARR(nr)
for i=0, N_ELEMENTS(hh)-1 do begin
ii=ri[ri[i]:ri[i+1]-1]
xr[i]=median(x[ii])
yr[i]=median(y[ii])
zr[i]=mean(z[ii])
endfor
|