; $Id: image_background.pro, v 1.1 Apr 2000 e.d. $
;
;+
; NAME:
;	IMAGE_BACKGROUND
;
; PURPOSE:
;	Compute the intensity distribution of the background in a given image.
;	Local measurements of the background are obtained after partitioning the
;	input array; the sampled distribution is then magnified onto the same
;	grid of the original image by interpolation.
;	This method is especially useful when the image background is not
;	spatially uniform.
;
; CATEGORY:
;	Signal processing.
;
; CALLING SEQUENCE:
;	Result = IMAGE_BACKGROUND(Image, Step)
;
; INPUTS:
;	Image:	2D array
;
;	Step:	Size of sub-regions to measure local background
;
; KEYWORD PARAMETERS:
;	SKY_MEDBOX:	The array of local measurements is smoothed by median
;		filtering. Use the keyword SKY_MEDBOX to fix a box size for median
;		smoothing. In general, the larger the Step size the smaller should
;		be the box for median filtering. The default is SKY_MEDBOX = 3,
;		which represents also the minimum box size.
;
;	CUBIC:	Set this keyword to magnify the sampled background with the
;		cubic convolution interpolation method implemented in the library
;		routine CONGRID. The default is to use bilinear interpolation as;
;		implemented in the IDL intrinsic function REBIN.
;
; OUTPUTS:
;	Result:	2D array, having the same size as the input Image.
;
; RESTRICTIONS:
;	1) Bad pixels are assumed to have been corrected.
;	2) The background may be over-estimated in the presence of strong
;	sources, even though the present method seems to be less biased
;	than median smoothing.
;
; PROCEDURE:
;	The input Image is partitioned into sub-frames defined by a rectangular
;	grid. The mesh points coincide with the centers of the corresponding
;	regions. For each sub-image a local background estimate is obtained by
;	means of the IDL library routine SKY, written by W.Landsman and adapted
;	from the DAOPHOT routine of the same name. If SKY doesn't converge, the
;	local background is estimated as the median value of the sub-image.
;	Each local background estimate is associated to the mesh point
;	corresponding to the sub-image under examination. The sampled version of
;	the background is smoothed by median filtering and magnified onto the
;	same grid of the input Image by means of interpolation.
;
; MODIFICATION HISTORY:
; 	Written by:	Emiliano Diolaiti, August 1999.
;	Updates:
;	1) Corrected edge effect due to 0-padding
;	   (Emiliano Diolaiti, October 1999).
;	2) Fixed bug on small arrays (Emiliano Diolaiti, April 2000).
;	3) Introduced error handler to handle run-time error generated by call
;	   to SKY (Emiliano Diolaiti, April 2000).
;-

FUNCTION image_background, image, sampling_step, $
						   SKY_MEDBOX = med_box, CUBIC = cubic

	on_error, 2
	siz = size52(image, /DIM)  &  sx = siz[0]  &  sy = siz[1]
	step = round(sampling_step)
	if  n_elements(med_box) ne 0  then  box = med_box  else  box = 3
	box = box > 3
	; Pad the image with 0s
	sx = (sx / step + 1) * step  &  sy = (sy / step + 1) * step
	ima = extend_array(image, sx, sy, OFFSET = off)
	lo = off  &  up = lo + siz - 1
	; Define the sampling grid
	nx = sx / step  &  ny = sy / step
	x = sampling_grid(nx, step, lx0, ux0)  &  ux0 = ux0 < sx
	y = sampling_grid(ny, step, ly0, uy0)  &  uy0 = uy0 < sy
	lx = round(x - step / 2.) > lo[0]  &  ux = round(x + step / 2.) < up[0]
	ly = round(y - step / 2.) > lo[1]  &  uy = round(y + step / 2.) < up[1]
	; Background sampling
	b = fltarr(nx, ny)
	for  i = 0L, ny - 1  do  for  j = 0L, nx - 1  do begin
	   catch, error
	   if  error ne 0  then begin
	      local_back = median(ima[lx[j]:ux[j],ly[i]:uy[i]], /EVEN)
	      catch, /CANCEL
	   endif else begin
	      sky, /SILENT, ima[lx[j]:ux[j],ly[i]:uy[i]], local_back, stdev
	      if  stdev lt 0  then $
	         local_back = median(ima[lx[j]:ux[j],ly[i]:uy[i]], /EVEN)
	   endelse
	   b[j,i] = local_back
	endfor
	; Sampled background smoothing
	b = median_filter(b, box)
	; Background interpolation
	siz = size52(b, /DIM) * step
	if  n_elements(b) gt 1  then begin
	   if  keyword_set(cubic)  then $
	      b = congrid(b, siz[0], siz[1], CUBIC = -0.5)  else $
	      b = rebin(b, siz[0], siz[1])
	endif else  b = make_array(siz, siz, VALUE = b)
	b = b[lo[0]:up[0],lo[1]:up[1]]
	return, b
end