include	"../gcombine.h"

$for (sird)
# G_MEDIAN -- Median of lines
#
# CYZhang 6 May, 1994 

procedure g_median$t (data, id, nimages, n, npts, median)

# Calling arguments
pointer	data[ARB]		# Input data line pointers
pointer	id[ARB]			# IDs pointers
int	nimages			# Number of images
int	n[npts]			# Number of retained pixels
int	npts			# Number of output points per line
$if (datatype == sil)
real	median[npts]		# Median
$else
PIXEL	median[npts]		# Median
$endif

# Local variables
pointer	sp, work
int	i, j, n1

include	"../gcombine.com"

begin
	
	call smark (sp)
	call salloc (work, nimages, TY_PIXEL)
	call gc_sort$t (data, Mem$t[work], n, npts)
	do i = 1, npts {
	    n1 = n[i]
	    if (n1 > 0) {
		do j = 1, n1
	    	    Mem$t[work+j-1] = Mem$t[data[j]+i-1]
		call g_med$t (Mem$t[work], n1, median[i])
	    } else 
		median[i] = BLANK
	}
	call sfree (sp)
end

# G_MED -- Median of an array
#
# CYZhang April 19, 1994

procedure g_med$t (a, n, med)

PIXEL	a[n]
int	n, n2
$if (datatype == sil)
real	med, low, high
$else
PIXEL	med, low, high
$endif

include	"../gcombine.com"

begin
	if (n == 0) {
	    med = BLANK
	    return
	}
	else if (n == 1) {
	    med = a[1]
	    return
	} else if (n == 2) {
	    low = a[1]
	    high = a[2]
	    med = (low + high) / 2.
	    return
	} 

	# Median
	if (n >= 3) {
	    n2 = 1 + n / 2
	    if (mod (n, 2) == 0) {
	        low = a[n2-1]
	        high = a[n2]
	        med = (low + high) / 2.
	    } else
	        med = a[n2]
	}
end
$endfor