include "wpdef.h"

.help median
.nf ----------------------------------------------------------------------------
              COMBINING IMAGES: MEDIAN ALGORITHM

The input images are combined by scaling and taking the median.  The exposure 
time of the output image is the scaled and weighted average of the input 
exposure times.  If some of the input images are real datatypes and the 
output image is short datatype there will be truncation.

PROCEDURES:

    MEDIAN	 -- Median of lines (no scaling).
    DQMEDIAN	 -- Scaled median of lines excluding bad pixels.
    SCMEDIAN	 -- Scaled median of lines.
    BIGSORT	 -- Sort by increasing value. Heapsort used for large arrays.
    SMALLSORT	 -- Sort by increasing value. Straight Insertion for small arrays.
.endhelp -----------------------------------------------------------------------

$for (silrdx)

#################################################################################
# MEDIAN --	Determine the median of image lines with no scaling.  This 	#
#		routine is based upon the `images.imcombine' package.  		#
#										#
#		Development version:	1/91	RAShaw				#

procedure median$t (data, median, nimages, npts)

# Calling arguments:
pointer		data[nimages]		# Input data line pointers
$if (datatype == sil)
real		median[npts]		# Output data line
$else
PIXEL		median[npts]		# Output data line
$endif
int		nimages			# Number of images to combine
int		npts			# Number of pixels per image line

# Local variables:
int		half			# Half the number of used images @pixel
int		i, j			# Loop indexs
$if (datatype == sil)
real		work[IMS_MAX]		# Work array
$else
PIXEL		work[IMS_MAX]		# Work array
$endif

begin

# Initialize working array. 
	do j = 1, npts {
	    do i = 1, nimages {
		work[i] = Mem$t[data[i]+j-1]
	    }

# Sort pixel values into increasing order.
	    switch (nimages) {
	    case 0, 1:
		return
	    case 2:
		;
	    case 3, 4, 5, 6, 7:
		call smallsort$t (work, nimages)
	    default:
		call bigsort$t (work, nimages)
	    }

# Select median value.  For an even number of elements compute the mean of the 
# two values @midpoint.  
	    half = nimages / 2
	    if (half*2 < nimages) 
		median[j] = work[half+1]
	    else
		median[j] = (work[half] + work[half+1]) / 2.
	}
end


#################################################################################
# DQMEDIAN --	Combine the images by scaling and taking the median, excluding 	#
#		bad pixels.  This routine is based upon the `images.imcombine' 	#
#		package.  							#
#										#
#		Development version:	1/91	RAShaw				#

procedure dqmedian$t (data, dqfdata, median, nimages, npts)

include "wpcom.h"

# Calling arguments:
pointer		data[nimages]		# Input data line pointers
int		dqfdata[nimages]	# Data Quality File pointers
$if (datatype == sil)
real		median[npts]		# Output data line
$else
PIXEL		median[npts]		# Output data line
$endif
int		nimages			# Number of images to combine
int		npts			# Number of pixels per image line

# Local variables:
int		bflag[IMS_MAX]		# 
int		half			# Half the number of used images @pixel
int		i, j			# Loop indexes
int		ncount			# Number of non-rejected images @pixel
$if (datatype == sil)
real		work[IMS_MAX]		# Scaled, non-flagged data
$else
PIXEL		work[IMS_MAX]		# Scaled, non-flagged data
$endif

begin
	do j = 1, npts {

# Select user-chosen Data Quality bits.
	    do i = 1, nimages
		bflag[i] = Memi[dqfdata[i]+j-1]
	    call aandki (bflag, BADBITS, bflag, npts)

# Initialize working array. 
	    ncount = 0
	    do i = 1, nimages {
		if (bflag[i] == 0) {
		    ncount       = ncount + 1
		    work[ncount] = Mem$t[data[i]+j-1] / SCALES[i] - ZEROS[i]
		} else
		    Mem$t[data[i]+j-1] = INDEF
	    }

# Sort pixel values into increasing order.
	    switch (ncount) {
	    case 0:
		median[j] = BLANK
		next
	    case 1, 2:
		;
	    case 3, 4, 5, 6, 7:
		call smallsort$t (work, ncount)
	    default:
		call bigsort$t (work, ncount)
	    }

# Select median value.  For an even number of elements compute the mean of the 
# two values @midpoint.  
	    half = ncount / 2
	    if (half*2 < ncount) 
		median[j] = work[half+1]
	    else
		median[j] = (work[half] + work[half+1]) / 2.
	}
end


#################################################################################
# SCMEDIAN --	Combine the images by scaling and taking the median, excluding 	#
#		bad pixels.  This routine is based upon the `images.imcombine' 	#
#		package.  							#
#										#
#		Development version:	1/91	RAShaw				#

procedure scmedian$t (data, median, nimages, npts)

include "wpcom.h"

# Calling arguments:
pointer		data[nimages]		# Input data line pointers
$if (datatype == sil)
real		median[npts]		# Output data line
$else
PIXEL		median[npts]		# Output data line
$endif
int		nimages			# Number of images to combine
int		npts			# Number of pixels per image line

# Local variables:
int		half			# Half the number of used images @pixel
int		i, j			# Loop indexes
$if (datatype == sil)
real		work[IMS_MAX]		# Scaled, non-flagged data
$else
PIXEL		work[IMS_MAX]		# Scaled, non-flagged data
$endif

begin

# Initialize working array.
	do j = 1, npts {
	    do i = 1, nimages {
		work[i] = Mem$t[data[i]+j-1] / SCALES[i] - ZEROS[i]
	    }

# Sort pixel values into increasing order.
	    switch (nimages) {
	    case 0, 1:
		return
	    case 2:
		;
	    case 3, 4, 5, 6, 7:
		call smallsort$t (work, nimages)
	    default:
		call bigsort$t (work, nimages)
	    }

# Select median value.  For an even number of elements compute the mean of the 
# two values @midpoint.  
	    half = nimages / 2
	    if (half*2 < nimages) 
		median[j] = work[half+1]
	    else
		median[j] = (work[half] + work[half+1]) / 2.
	}
end


#################################################################################
# BIGSORT --	Sort array "work" of length "n" into ascending numerical order 	#
#		using the Heapsort algorithm found in "Numerical Recipies".  	#
#		The array "work" is replaced on output by its sorted rearrange-	#
#		ment.  								#
#										#
#		Development version:	1/91	RAShaw				#

procedure bigsort$t (work, nc)

# Calling arguments:
$if (datatype == sil)
real		work[nc]	# Array of values
$else
PIXEL		work[nc]	# Array of values
$endif
int		nc		# Number of values to be sorted

# Local variables:
int		i, ir, j, m	# Dummy indexes
$if (datatype == sil)
real		temp		# Temporary value
$else
PIXEL		temp		# Temporary value
$endif

begin
	m  = nc / 2 + 1
	ir = nc

# The index "m" will be decremented from its initial value down to 1 during the 
# heap creation phase.  Once it reaches 1, the index "ir" will be decremented 
# from its initial value down to 1 during the heap selection phase. 
	repeat {
	    if (m > 1) {
		m    = m - 1
		temp = work[m]
	    } else {				# 
		temp     = work[ir]		# Clear a space @end of array & 
		work[ir] = work[1]		# retire top of heap into it.
		ir       = ir - 1		# 
		if (ir == 1) {			# Done with the last promotion?
		    work[1] = temp		# The lowest value
		    return
		}
	    }
	    i = m
	    j = m + m
	    while (j <= ir) {		# Sift "temp" down to its proper level
		if (j < ir) {
		    $if (datatype == x)
		    if (abs (work[j]) < abs (work[j+1]))
		    $else
		    if (work[j] < work[j+1])
		    $endif
			j = j + 1
		}
		$if (datatype == x)
		if (abs (temp) < abs (work[j])) {
		$else
		if (temp < work[j]) {		# Demote "temp"
		$endif
		    work[i] = work[j]
		    i = j
		    j = j + j
		} else				# Correct level for "temp" 
		    j = ir + 1			# Set "j" to terminate sift-down
	    }
	    work[i] = temp			# Put "temp" into its slot
	}
end


#################################################################################
# SMALLSORT --	Sort vector by increasing value.  This algorithm is based on 	#
#		the Straight Insertion routine ("PIKSRT") found in "Numerical 	#
#		Recipies", and is best for small vectors.  The array WORK is 	#
# 		replaced on output by its sorted rearrangement.  		#
#										#
#		Development version:	1/91	RAShaw				#

procedure smallsort$t (work, nc)

# Calling arguments:
$if (datatype == sil)
real		work[nc]		# Working array of data values
$else
PIXEL		work[nc]		# Working array of data values
$endif
int		nc			# Number of input images used per pixel.

# Local variables:
int		i, j			# Loop indexes
$if (datatype == sil)
real		temp			# Temporary value
$else
PIXEL		temp			# Temporary value
$endif

begin

# Pick out each element in turn.
	do j = 2, nc {
	    temp = work[j]

# Look for the place to insert it.
	    do i = j-1, 1, -1 {
		$if (datatype == x)
		if (abs(work[i]) <= abs(temp)) goto 10
		$else
		if (work[i] <= temp) goto 10
		$endif
		work[i+1] = work[i]
	    }
	    i = 0
10	    work[i+1] = temp
	}
end
$endfor