procedure borderstat (image,border)

file	image	=""	{prompt="image name"}
int	border	=1	{prompt="width of image border for statistics"}
int	npix	=0	{prompt="number of image pixels"}
real	mean	=0.	{prompt="average pixel value"}
real	midpt	=0.	{prompt="midpoint (median) pixel value"}
real	stddev	=0.	{prompt="standard deviation (sigma)"}
real	rmin	=0.	{prompt="minimum pixel value"}
real	rmax	=0.	{prompt="maximum pixel value"}
real	sum	=0.	{prompt="sum of all pixels"}
real	rmode	=0.	{prompt="mode of the pixel values"}
real	lower	=INDEF	{prompt="lower cutoff for pixel values"}
real	upper	=INDEF	{prompt="upper cutoff for pixel values"}
real	binwidth=0.1	{prompt="bin width of histogram in sigma"}
bool	domode	=no	{prompt="calculate mode of image pixels?"}
bool	verbose	=no	{prompt="print results to screen?"}

# @(#) task borderstat                Author: R.A. Jansen -- Feb 12 1996
# @(#)
# @(#) task borderstat runs task 'imgstat' and loads the statistical pa-
# @(#) rameters into different task parameters, that can be used outside
# @(#) the procedure (e.g. borderstat.mean). Note that the mode is actu-
# @(#) ally the mean of the modes of the four border sections.
# @(#)
# @(#) Apr 24 2001 -- Complete rewrite of this task: added support for
# @(#)                multi-extension / N-dim. images, got rid of ncol, 
# @(#)                nrow as required parameters, and uses output from
# @(#)                'imstat' directly, rather than calling 'imgstat',
# @(#)                so subsequent calls to these tasks won't conflict.

begin
	file	input
	string	sec1, sec2, sec3, sec4
	bool	domod
	int	nrow, ncol, x1, x2, y1, y2, rnpix, bd
	real	hi, lo, binw, rmean, rsum, rstddev, rrmode

	input = image
        bd = border
	domod = domode
        binw = binwidth
        hi = upper
        lo = lower

	imcopy (input, "_tmpbrdrstt.fits", verbose=no)

        imgets ("_tmpbrdrstt.fits", "naxis1") ; ncol = int(imgets.value)
        imgets ("_tmpbrdrstt.fits", "naxis1") ; ncol = int(imgets.value)
        imgets ("_tmpbrdrstt.fits", "naxis2") ; nrow = int(imgets.value)

        x1 = bd + 1
        x2 = ncol - bd
        y1 = bd + 1
        y2 = nrow - bd
        sec1 = "[1:"//str(bd)//",1:"//str(y2)//"]"
        sec2 = "[1:"//str(x2)//","//str(y2+1)//":"//str(nrow)//"]"
        sec3 = "["//str(x2+1)//":"//str(ncol)//","//str(y1)//":"//str(nrow)//"]"
        sec4 = "["//str(x1)//":"//str(ncol)//",1:"//str(bd)//"]"

        rnpix   = 0
        rmean   = 0.
        rsum    = 0.
        rstddev = 0.
        rmin    = 9999999.
        rmax    = -9999999.
        rrmode  = 0.

        if (domod) {
           imgstat ("_tmpbrdrstt.fits"//sec1, lower=lo, upper=hi, binwidth=binw, domode=yes)
           rnpix   = rnpix   + imgstat.npix
           rmean   = rmean   + imgstat.mean
           rsum    = rsum    + imgstat.sum
           rstddev = rstddev + imgstat.stddev
           rmin	   = min(rmin,imgstat.min)
           rmax    = max(rmax,imgstat.max)
           rrmode  = rrmode  + imgstat.rmode

           imgstat ("_tmpbrdrstt.fits"//sec2, lower=lo, upper=hi, binwidth=binw, domode=yes)
           rnpix   = rnpix   + imgstat.npix
           rmean   = rmean   + imgstat.mean
           rsum    = rsum    + imgstat.sum
           rstddev = rstddev + imgstat.stddev
           rmin	   = min(rmin,imgstat.min)
           rmax    = max(rmax,imgstat.max)
           rrmode  = rrmode  + imgstat.rmode

           imgstat ("_tmpbrdrstt.fits"//sec3, lower=lo, upper=hi, binwidth=binw, domode=yes)
           rnpix   = rnpix   + imgstat.npix
           rmean   = rmean   + imgstat.mean
           rsum    = rsum    + imgstat.sum
           rstddev = rstddev + imgstat.stddev
           rmin	   = min(rmin,imgstat.min)
           rmax    = max(rmax,imgstat.max)
           rrmode  = rrmode  + imgstat.rmode

           imgstat ("_tmpbrdrstt.fits"//sec4, lower=lo, upper=hi, binwidth=binw, domode=yes)
           rnpix   = rnpix + imgstat.npix
           rmean   = 0.25*(rmean + imgstat.mean)
           rsum    = rsum + imgstat.sum
           rstddev = 0.25*(rstddev + imgstat.stddev)
           rmin	   = min(rmin,imgstat.min)
           rmax    = max(rmax,imgstat.max)
           rrmode  = 0.25*(rrmode + imgstat.rmode)

	} else {
           imgstat ("_tmpbrdrstt.fits"//sec1)
           rnpix   = rnpix   + imgstat.npix
           rmean   = rmean   + imgstat.mean
           rsum    = rsum    + imgstat.sum
           rstddev = rstddev + imgstat.stddev
           rmin	   = min(rmin,imgstat.min)
           rmax    = max(rmax,imgstat.max)

           imgstat ("_tmpbrdrstt.fits"//sec2, domode=no)
           rnpix   = rnpix   + imgstat.npix
           rmean   = rmean   + imgstat.mean
           rsum    = rsum    + imgstat.sum
           rstddev = rstddev + imgstat.stddev
           rmin	   = min(rmin,imgstat.min)
           rmax    = max(rmax,imgstat.max)

           imgstat ("_tmpbrdrstt.fits"//sec3, domode=no)
           rnpix   = rnpix   + imgstat.npix
           rmean   = rmean   + imgstat.mean
           rsum    = rsum    + imgstat.sum
           rstddev = rstddev + imgstat.stddev
           rmin	   = min(rmin,imgstat.min)
           rmax    = max(rmax,imgstat.max)

           imgstat ("_tmpbrdrstt.fits"//sec4, domode=no)
           rnpix   = rnpix + imgstat.npix
           rmean   = 0.25*(rmean + imgstat.mean)
           rsum    = rsum + imgstat.sum
           rstddev = 0.25*(rstddev + imgstat.stddev)
           rmin	   = min(rmin,imgstat.min)
           rmax    = max(rmax,imgstat.max)
        }

        # Log results to the screen
        if (verbose) {
           print ("npix   mean   stddev     min     max    rmode    sum")
           print (rnpix, rmean, rstddev, rmin, rmax, rrmode, rsum, )
        }

	# Clean-up
	imdelete ("_tmpbrdrstt.fits", yes, verify=no, default_acti=yes)

	# Return results as task parameters
        npix   = rnpix
        mean   = rmean
        sum    = rsum
        stddev = rstddev
        rmode   = rrmode

end