function adapt_smooth,im,snr,rmax,snrim=snrim, expmap=expmap, thresh=thresh

	if n_params(0) eq 0 then begin
	  print,'smim = adapt_smooth(im,snr,rmax,snrm=snrim, expmap=expmap, thresh=thresh)'
	  print,'Adaptively smooth im until each pixel contains a minimum of
	  print,'snr^2 counts (resultant pixel will contain snr^2/n, n = no. of
	  print,'total pixels required in sum for pixel to have n counts'
	  print,'rmax = max distance to continue summing pixels'
	  print,'These parameters are analogous to the parameters for setplot rebin in XSPEC'
	  print,'snrim contains raw number of counts summed into each pixel (i.e., snr^2'
          print,'expmap contains optional exposure map to simultaneously smooth and divide'
          print,'  i.e., if units are secs, the output would be in cts/sec'
          print,'  only consider pixels where expmap > thresh (def. thresh=0)
	  retall
	endif

	sz = size(im)
	sm = fltarr(sz(1),sz(2))
	
	sn = snr^2
	snrim = im
;        emap = n_elements(expmap) gt 0  
        if n_elements(expmap) eq 0 then expmap = replicate(1l, sz(1), sz(2)) $
          else print,'Output will be divided by expmap'
        if n_elements(thresh) eq 0 then thresh = 0
	for j=0,sz(2)-1 do begin
         w = where(expmap(*,j) gt thresh, n)
         if n gt 0 then for i=min(w), max(w) do begin
;	 for i=0,sz(1)-1 do begin
	  c = long(im(i,j))
;	  n = 1.0
	  r = 1
          t = expmap(i,j)
          en = 0
	  while r lt rmax and c lt sn do begin
	    i1 = (i-r)>0
	    i2 = (i+r)<sz(1)-1
	    j1 = (j-r)>0
	    j2 = (j+r)<sz(2)-1
            temp = long(im(i1:i2, j1:j2))
; get no. of pixels in current sum region
;            if emap then begin
              tempe = expmap(i1:i2, j1:j2)
              we = where(tempe gt thresh, en)
              if en gt 0 then c = total(temp(we)) else c = 0l
;            endif else begin
;	      n = (i2-i1+1)*(j2-j1+1)
;	      c = total(temp)
;            endelse
;            n = float(n)
	    r = r+1
	  endwhile
;          t = 1.0
          if (en gt 0) then t = total(tempe(we))
; t = mean exposure of region
	  snrim(i,j) = c 	  
	  if t gt 0 then sm(i,j) = c/t
;	  print,'Pixel (',strn(i),',',strn(j),'), summing (',strn(i1),'-',$
;	      strn(i2),',',strn(j1),'-',strn(j2),'), n = ',strn(n),', c = ',$
;	      strn(c)
	 endfor
	 print,sz(2)-j
	endfor
	return,sm
	end