# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.

include	"../icombine.h"

$for (sird)
# IC_MEDIAN -- Median of lines

procedure ic_median$t (d, n, npts, doblank, median)

pointer	d[ARB]			# Input data line pointers
int	n[npts]			# Number of good pixels
int	npts			# Number of output points per line
int	doblank			# Set blank values?
$if (datatype == sil)
real	median[npts]		# Median
$else
PIXEL	median[npts]		# Median
$endif

int	i, j, k, j1, j2, n1, lo, up, lo1, up1
bool	even
$if (datatype == silx)
real	val1, val2, val3
$else
PIXEL	val1, val2, val3
$endif
PIXEL	temp, wtemp
$if (datatype == x)
real	abs_temp
$endif

include	"../icombine.com"

begin
	# If no data return after possibly setting blank values.
	if (dflag == D_NONE) {
	    if (doblank == YES) {
		do i = 1, npts
		    median[i]= blank
	    }
	    return
	}

	# If the data were previously sorted then directly compute the median.
	if (mclip) {
	    if (dflag == D_ALL) {
		n1 = n[1]
		even = (mod (n1, 2) == 0)
		j1 = n1 / 2 + 1
		j2 = n1 / 2
		do i = 1, npts {
		    k = i - 1
		    if (even) {
			val1 = Mem$t[d[j1]+k]
			val2 = Mem$t[d[j2]+k]
			median[i] = (val1 + val2) / 2.
		    } else
			median[i] = Mem$t[d[j1]+k]
		}
	    } else {
		do i = 1, npts {
		    k = i - 1
		    n1 = n[i]
		    if (n1 > 0) {
			j1 = n1 / 2 + 1
			if (mod (n1, 2) == 0) {
			    j2 = n1 / 2
			    val1 = Mem$t[d[j1]+k]
			    val2 = Mem$t[d[j2]+k]
			    median[i] = (val1 + val2) / 2.
			} else
			    median[i] = Mem$t[d[j1]+k]
		    } else if (doblank == YES)
			median[i] = blank
		}
	    }
	    return
	}

	# Compute the median.
	do i = 1, npts {
	    k = i - 1
	    n1 = n[i]

	    # If there are more than 3 points use Wirth algorithm.  This
	    # is the same as vops$amed.gx except for an even number of
	    # points it selects the middle two and averages.
	    if (n1 > 3) {
		lo = 1
		up = n1
		j  = max (lo, min (up, (up+1)/2))

		while (lo < up) {
		    if (! (lo < up))
			break

		    temp = Mem$t[d[j]+k];  lo1 = lo;  up1 = up
		    $if (datatype == x)
		    abs_temp = abs (temp)
		    $endif

		    repeat {
			$if (datatype == x)
			while (abs (Mem$t[d[lo1]+k]) < abs_temp)
			$else
			while (Mem$t[d[lo1]+k] < temp)
			$endif
			    lo1 = lo1 + 1
			$if (datatype == x)
			while (abs_temp < abs (Mem$t[d[up1]+k]))
			$else
			while (temp < Mem$t[d[up1]+k])
			$endif
			    up1 = up1 - 1
			if (lo1 <= up1) {
			    wtemp = Mem$t[d[lo1]+k]
			    Mem$t[d[lo1]+k] = Mem$t[d[up1]+k]
			    Mem$t[d[up1]+k] = wtemp
			    lo1 = lo1 + 1;  up1 = up1 - 1
			}
		    } until (lo1 > up1)

		    if (up1 < j)
			lo = lo1
		    if (j < lo1)
			up = up1
		}

		median[i] = Mem$t[d[j]+k]

		if (mod (n1,2) == 0) {
		    lo = 1
		    up = n1
		    j  = max (lo, min (up, (up+1)/2)+1)

		    while (lo < up) {
			if (! (lo < up))
			    break

			temp = Mem$t[d[j]+k];  lo1 = lo;  up1 = up
			$if (datatype == x)
			abs_temp = abs (temp)
			$endif

			repeat {
			    $if (datatype == x)
			    while (abs (Mem$t[d[lo1]+k]) < abs_temp)
			    $else
			    while (Mem$t[d[lo1]+k] < temp)
			    $endif
				lo1 = lo1 + 1
			    $if (datatype == x)
			    while (abs_temp < abs (Mem$t[d[up1]+k]))
			    $else
			    while (temp < Mem$t[d[up1]+k])
			    $endif
				up1 = up1 - 1
			    if (lo1 <= up1) {
				wtemp = Mem$t[d[lo1]+k]
				Mem$t[d[lo1]+k] = Mem$t[d[up1]+k]
				Mem$t[d[up1]+k] = wtemp
				lo1 = lo1 + 1;  up1 = up1 - 1
			    }
			} until (lo1 > up1)

			if (up1 < j)
			    lo = lo1
			if (j < lo1)
			    up = up1
		    }
		    median[i] = (median[i] + Mem$t[d[j]+k]) / 2
		}

	    # If 3 points find the median directly.
	    } else if (n1 == 3) {
		$if (datatype == x)
	        val1 = abs (Mem$t[d[1]+k])
	        val2 = abs (Mem$t[d[2]+k])
	        val3 = abs (Mem$t[d[3]+k])
	        if (val1 < val2) {
		    if (val2 < val3)		# abc
		        median[i] = Mem$t[d[2]+k]
		    else if (val1 < val3)	# acb
		        median[i] = Mem$t[d[3]+k]
		    else			# cab
		        median[i] = Mem$t[d[1]+k]
	        } else {
		    if (val2 > val3)		# cba
		        median[i] = Mem$t[d[2]+k]
		    else if (val1 > val3)	# bca
		        median[i] = Mem$t[d[3]+k]
		    else			# bac
		        median[i] = Mem$t[d[1]+k]
	        }
		$else
	        val1 = Mem$t[d[1]+k]
	        val2 = Mem$t[d[2]+k]
	        val3 = Mem$t[d[3]+k]
	        if (val1 < val2) {
		    if (val2 < val3)		# abc
		        median[i] = val2
		    else if (val1 < val3)	# acb
		        median[i] = val3
		    else			# cab
		        median[i] = val1
	        } else {
		    if (val2 > val3)		# cba
		        median[i] = val2
		    else if (val1 > val3)	# bca
		        median[i] = val3
		    else			# bac
		        median[i] = val1
	        }
		$endif

	    # If 2 points average.
	    } else if (n1 == 2) {
		val1 = Mem$t[d[1]+k]
		val2 = Mem$t[d[2]+k]
		median[i] = (val1 + val2) / 2

	    # If 1 point return the value.
	    } else if (n1 == 1)
		median[i] = Mem$t[d[1]+k]

	    # If no points return with a possibly blank value.
	    else if (doblank == YES)
		median[i] = blank
	}
end
$endfor