sigrange.pro

basic median filter simulation

	FUNCTION SIGRANGE,ARRAY,FRACTION=FRACTION,MISSING=MISSING,RANGE=RANGE
;+
; NAME: 
;	SIGRANGE()
; PURPOSE: 
;	Selects the most significant data range in an image.
; EXPLANATION: 
;	Selects out the most significant range in the data to be used in 
;	displaying images.  The histogram of ARRAY is used to select the most
;	significant range.      Useful for scaling an image display.
; CALLING SEQUENCE: 
;	OUTPUT = SIGRANGE( ARRAY )
; INPUTS: 
;	ARRAY	 = Array to take most significant range of.
; OPTIONAL INPUTS: 
;	None.
; OUTPUTS: 
;	The function returns an array where values above and below the
;	selected range are set equal to the maximum and minimum of the
;	range respectively.
; OPTIONAL INPUT KEYWORDS: 
;	FRACTION = Fraction of data to consider most significant.
;		   Defaults to 0.99
;	MISSING	 = Value used to flag missing points.  Data points with this
;		   value are not considered or changed.
; OPTIONAL OUTPUT KEYWORD
;	RANGE    = 2 element vector, giving the range (minimum and maxmimum) 
;		used
;
; NOTES:
;       If the image array contains more than 10,000 points then SIGRANGE() 
;       uses random indexing of a subset of the points to determine the range
;       (for speed).    Thus identical calls to SIGRANGE() might not yield
;       identical results (although they should be very close).     
; RESTRICTIONS: 
;	ARRAY must have more than two points.  Fraction must be greater than 0 
;	and less than 1.
;
;	SIGRANGE was originally part of the SERTS image display package.   
;	Other routines from this package are available at 
;
;	http://sohowww.nascom.nasa.gov/solarsoft/gen/idl/image/
;
;	Note that this version of SIGRANGE does not include the non-standard 
;	system variables used in the SERTS package.
; REVISION HISTORY: 
;	Version 1, William Thompson, GSFC, 12 May 1993.
;		Incorporated into CDS library.
;	Version 2, William Thompson, GSFC, 25 May 1993.
;		Changed call to HISTOGRAM to be compatible with OpenVMS/ALPHA
;       Version 3, CDP, RAL, Add RANGE keyword.  16-Apr-96
;	Version 4, William Thompson, GSFC, 17 April 1996
;		Corrected some problems when range is too high.
;	Version 5, 13-Jan-1998, William Thompson, GSFC
;		Use random numbers to improve statistics when only using a
;		fraction of the array.
;	Version 6, 06-Mar-1998, William Thompson, GSFC
;		Change default to 0.99
;-
;
	IF N_ELEMENTS(FRACTION) NE 1 THEN FRACTION = 0.99
	IF N_ELEMENTS(ARRAY) LE 2 THEN BEGIN
	    MESSAGE, /CONTINUE, 'Not enough points to form histogram'
	    RETURN, ARRAY
	END ELSE IF (FRACTION LE 0) OR (FRACTION GE 1) THEN BEGIN
	    MESSAGE, /CONTINUE, 'Fraction must be GT 0 and LT 1'
	    RETURN, ARRAY
	ENDIF
;
;  To speed up the process, work on a reduced version of ARRAY.
;
	IF N_ELEMENTS(ARRAY) LT 10000 THEN ATEMP0 = ARRAY ELSE BEGIN
	    NN = 1000 > (N_ELEMENTS(ARRAY) / 25) < 100000
	    ATEMP0 = ARRAY[N_ELEMENTS(ARRAY)*RANDOMU(SEED,NN)]
	ENDELSE
;
;  Get the total range of the data, excluding any missing points.
;
        IF N_ELEMENTS(MISSING) EQ 1 THEN BEGIN
                W = WHERE(ATEMP0 NE MISSING, COUNT)
                IF COUNT GT 0 THEN ATEMP0 = ATEMP0(W)
	ENDIF 
	N_TOTAL = N_ELEMENTS(ATEMP0)
	AMAX = 1.*MAX(ATEMP0)
	AMIN = 1.*MIN(ATEMP0)
	IF AMIN EQ AMAX THEN GOTO, EXIT_POINT
;
;  Set up some initial parameters for the reiteration.
;
	ATEMP = ATEMP0
	DELTA = 0
;
;  Form the histogram, and calculate an array expressing the fraction of points
;  that fall within or below the given bin.
;
FIND_RANGE:
	LAST_DELTA = DELTA
	X = AMIN  +  FINDGEN(1001) * (AMAX - AMIN) / 1000.
	H = HISTOGRAM(LONG((ATEMP-AMIN)*1000./(AMAX - AMIN)))
	FOR I = 1,N_ELEMENTS(H)-1 DO H[I] = H[I] + H[I-1]
	H = H / FLOAT(N_TOTAL)
;
;  Estimate the endpoints corresponding to the specified range, and calculate
;  the values at these endpoints.  Limit the array to be within these values.
;
	IMIN = (MIN( WHERE( H GT ((1. - FRACTION) / 2.) )) - 1) > 0
	IMAX =  MIN( WHERE( H GT ((1. + FRACTION) / 2.) ))
	IF IMAX LT 0 THEN IMAX = 1000
	AMIN = X[IMIN]
	AMAX = X[IMAX]
;
;  If the calculated range is zero, then use 2% of the full range of the data.
;
	IF AMAX EQ AMIN THEN BEGIN
		BMAX = MAX(ATEMP0, MIN=BMIN)
		AMAX = MAX(ATEMP0(WHERE(ATEMP0 LE (AMAX + 0.01*(BMAX-BMIN)))))
		AMIN = MIN(ATEMP0(WHERE(ATEMP0 GE (AMIN - 0.01*(BMAX-BMIN)))))
	ENDIF
;
;  If the range calculated has changed by more than 5% from the last iteration,
;  the reiterate.
;
	ATEMP = AMIN > ATEMP0 < AMAX
	DELTA = AMAX - AMIN
	RATIO = (DELTA - LAST_DELTA) / (DELTA + LAST_DELTA)
	IF ABS(RATIO) GT 0.05 THEN GOTO, FIND_RANGE
;
;  If a missing pixel flag value was passed, then reset those points to the
;  flag value.  Return the adjusted array.
;
EXIT_POINT:
	ATEMP = AMIN > ARRAY < AMAX
	IF N_ELEMENTS(MISSING) EQ 1 THEN BEGIN
		WW = WHERE(ARRAY EQ MISSING,N_MISSING)
		IF N_MISSING GT 0 THEN ATEMP[WW] = MISSING
	ENDIF
        RANGE = [AMIN,AMAX]
	RETURN, ATEMP
	END