include <imhdr.h>
include "gsky.h"

define	MAX_BINS	40000	#maximum number of bins in the histogram


# gsky_sky  --  Computes the sky for a set of images.
#
#
# 12-Jun-1996   Adapted from xcrrej to account for flagged pixels 
#               on input images (I.Busko)


procedure gsky_sky (ipin, ipmask, nfiles, nmasks, dim_x, dim_y, 
                    lower, upper, dqpat, skyval)

# input:
pointer	ipin[ARB]	# input images' descriptors
pointer	ipmask[ARB]	# input masks' descriptors
int	nfiles		# number of images
int	nmasks		# number of masks
int	dim_x, dim_y	# dimensions of the input image
real	lower, upper, 	# lower and upper limits of usable data
short	dqpat		# data quality bits

# output:
real	skyval[ARB]	# the sky level

# local:
pointer	histgrm		# pointer to the histogram
int	nbins		# number of bins
real	hwidth		# histogram resolution
real	hmin		# minimum histogram value
real	hmax		# maximum histogram value
char	text[SZ_LINE]
int	range
short	szero
pointer	sp, v, buf, buf2, buf_msk
int	i, k, npix
int     imgnlr()
short   imgnls(), ands()
#===============================================================================
begin
	szero = short (0)

	# first try to use upper and lower as boundary of the histogram, and
	# use the bin size of 1 (DN)
	range = nint(upper) - nint(lower)
	if (range > 1 && range < MAX_BINS) {
	    hmin = real (nint(lower))
	    hmax = real (nint(upper))
	    hwidth = 1.
	    nbins = range + 1

	# if the range is too large
	} else {
	    call sprintf (text, SZ_LINE, 
		    "Data range too large, must be less than %d")
	 	call pargi (MAX_BINS)
	    call error (1, text)
	}

	# set up the histogram and clean-pixel arrays
	call smark (sp)
	call salloc (v, IM_MAXDIM, TY_LONG)
	call malloc (histgrm, nbins, TY_INT)
	if (nmasks > 0)
	    call malloc (buf2, dim_x, TY_REAL)

	do k = 1, nfiles {

	    # reset the line-counter
            call amovkl (long(1), Meml[v], IM_MAXDIM)

	    # reset the histogram
	    do i = 1, nbins
		Memi[histgrm+i-1] = 0

	    # construct the histogram line by line
            while (imgnlr (ipin[k], buf, Meml[v]) != EOF) {

	        # read mask
	        if (nmasks > 0) {
	            i = imgnls (ipmask[k], buf_msk, Meml[v])

	            # discard flagged pixels
	            npix = 0
	            do i = 0, dim_x - 1 {
	                if (ands (dqpat, Mems[buf_msk+i]) == szero) {
	                    Memr[buf2+npix] = Memr[buf+i]
	                    npix = npix + 1
	                }
	            }

	            # build histogram from cleaned-up data
                    call ahgmr (Memr[buf2],npix,Memi[histgrm],nbins,hmin,hmax)

	        } else 

	            # build histogram from raw data
                    call ahgmr (Memr[buf],dim_x,Memi[histgrm],nbins,hmin,hmax)
	    }
	
	    # calculate the mode from the histogram
	    call gsky_mode (Memi[histgrm],nbins,hwidth,hmin,hmax, skyval[k])
	}

	if (nmasks > 0)
	    call mfree (buf2, TY_REAL)
	call mfree (histgrm, TY_INT)
	call sfree (sp)
end