;
;+
; NAME:
;
;       REJECT_1D
;
; PURPOSE:
;
;		Clean a 1-D array from the points out of a sigma threshold from
;		their neigborhood.
;		The outliers are deleted from the resulting array.
;
; CATEGORY:
;
;       Array processing.
;
; CALLING SEQUENCE:
;
;		REJECT_1D, x,y, x_ext, y_ext, box, threshold, iter, rejected=rejected, median=median
;
; INPUTS:
;
;		y: 1-D array to clean
;		x: 1-D array of the x coordinates of y.
;		   The x array must be monotonically increasing
;
;		box: width of the neigborhood to which compare each point.
;
;		threshold: sigma threshold for the clipping
;
;		iter: number of iterations to repeat the sigma-clipping
;
; KEYWORD PARAMETERS:
;
;		MEDIAN: if set, the median of the neigborhood is used instead of the mean
;				to compute the standard deviation.
;
;		REJECTED: output vector with the indexes of the rejected points
;
; OUTPUTS:
;
;		y_ext: 1-D array cleaned
;		x_ext: 1-D array of the x coordinates of y_ext
;
;
; CALLS:
;
;		The ROBUST_SIGMA of the ASTROLIB is called.
;
; MODIFICATION HISTORY:
;
;       Feb 2004 - 	Gianluca Li Causi, INAF - Rome Astronomical Observatory
;					licausi@mporzio.astro.it
;					http://www.mporzio.astro.it/~licausi/
;-

PRO reject_1d, x,y, x_ext, y_ext, box, soglia, iter, rejected=rejected, median=median
box2 = fix(box/2.)
x_ext = reform(x)
y_ext = reform(y)

np = n_elements(x_ext)
rejected = bytarr(np) + 1.
index0 = indgen(np)

FOR it = 0, iter-1 DO BEGIN			;ITERAZIONI

	np = n_elements(x_ext)
	reject = bytarr(np)


	y_ext0 = y_ext
	y_ext = [reverse(y_ext[0:box2]), y_ext, reverse(y_ext[np-box2: np-1])]
	FOR j = 0L, np-1 DO BEGIN

		je = j + box2
		loc_piece = [y_ext[je-box2:je-1],y_ext[je+1:je+box2]]
		IF keyword_set(median) THEN $
			loc_median = Med_Neigborhood(loc_piece)	$		;mediana locale senza il punto centrale
		ELSE $

			loc_median = Med_Neigborhood(loc_piece, /MEAN)	;media locale senza il punto centrale

		sigma = robust_sigma(loc_piece)
		dev = ABS((y_ext0[j] - loc_median[n_elements(loc_median)/2])/sigma)

		IF dev GE soglia THEN reject[j] = 1

	ENDFOR

	index = where(reject EQ 0)
	x_ext = x_ext[index]
	y_ext = y_ext0[index]
	index0 = index0[index]

ENDFOR

rejected[index0] = 0


END