Initial commit

1 files changed, 184 insertions, 0 deletions

diff --git a/noao/digiphot/apphot/fitsky/apmedian.x b/noao/digiphot/apphot/fitsky/apmedian.x
new file mode 100644
index 00000000..5998492c
--- /dev/null
+++ b/noao/digiphot/apphot/fitsky/apmedian.x
@@ -0,0 +1,184 @@
+include <mach.h>
+include "../lib/fitsky.h"
+
+define	MEDCUT	0.025
+
+# AP_MEDIAN -- Procedure to calculate the median of the sky pixel array.
+
+int procedure ap_median (skypix, coords, wgt, index, nskypix, snx, sny, losigma,
+	hisigma, rgrow, maxiter, sky_med, sky_sigma, sky_skew, nsky,
+	nsky_reject)
+
+real	skypix[ARB]		# unsorted array of sky pixels
+int	coords[ARB]		# coordinate array for regions growing
+real	wgt[ARB]		# array of weights for rejections
+int	index[ARB]		# sky pixel indices in sort order
+int	nskypix			# total number of sky pixels
+int	snx, sny		# dimensions of the sky subraster
+real	losigma, hisigma	# upper and lower sigma for rejection
+real	rgrow			# radius of region growing
+int	maxiter			# maximum number of cycles of rejection
+real	sky_med			# computed sky value
+real	sky_sigma		# the computed sky sigma
+real	sky_skew		# skewness of sky distribution
+int	nsky			# the number of sky pixels used
+int	nsky_reject		# the number of sky pixels rejected
+
+double  dsky, sumpx, sumsqpx, sumcbpx
+int	i, j, ilo, ihi, il, ih, med, medcut
+real	sky_zero, sky_mean, locut, hicut, dmin, dmax
+int	ap_grow_regions(), apimed()
+real	apsmed(), apwsmed(), ap_asumr()
+
+begin
+	# Intialize.
+	nsky = nskypix
+	nsky_reject = 0
+	sky_med = INDEFR
+	sky_sigma = INDEFR
+	sky_skew = INDEFR
+	if (nskypix <= 0)
+	    return (AP_NOSKYAREA)
+
+	# Sort the sky pixels and compute the median, mean, sigma and skew.
+	# MEDCUT tries to correct for quantization effects
+	call ap_ialimr (skypix, index, nskypix, dmin, dmax)
+	sky_zero = ap_asumr (skypix, index, nskypix) / nskypix
+	medcut = nint (MEDCUT * real (nskypix))
+	sky_med = apsmed (skypix, index, nskypix, medcut)
+	sky_med = max (dmin, min (sky_med, dmax))
+	call apfimoments (skypix, index, nskypix, sky_zero, sumpx, sumsqpx,
+	    sumcbpx, sky_mean, sky_sigma, sky_skew)
+	sky_mean = max (dmin, min (sky_mean, dmax))
+	if (maxiter < 1 || (IS_INDEFR(losigma) && IS_INDEFR(hisigma)) ||
+	    sky_sigma <= 0.0)
+	    return (AP_OK)
+	#call printf (
+	    #"mean=%g med=%g sigma=%g nsky=%d nrej=%d dmin=%g dmax=%g\n")
+	    #call pargr (sky_mean)
+	    #call pargr (sky_med)
+	    #call pargr (sky_sigma)
+	    #call pargi (nsky)
+	    #call pargi (nsky_reject)
+	    #call pargr (dmin)
+	    #call pargr (dmax)
+
+	# Reject points outside losigma * sky_sigma and hisigma * sky_sigma
+	# of the median.
+	ilo = 1
+	ihi = nskypix
+	do i = 1, maxiter {
+
+	    # Compute the new rejection limits.
+	    if (IS_INDEFR(losigma))
+		locut = max (-MAX_REAL, dmin)
+	    else if (i == 1)
+		locut = sky_med - min (sky_med - dmin, dmax - sky_med,
+		    losigma * sky_sigma)
+	    else
+	        locut = sky_med - losigma * sky_sigma
+	    if (IS_INDEFR(hisigma))
+		hicut = min (MAX_REAL, dmax)
+	    else if (i == 1)
+		hicut = sky_med + min (dmax - sky_med, sky_med - dmin,
+		    hisigma * sky_sigma)
+	    else
+	        hicut = sky_med + hisigma * sky_sigma
+	    #call printf ("    locut=%g hicut=%g\n")
+		#call pargr (locut)
+		#call pargr (hicut)
+
+	    # Detect pixels to be rejected.
+	    for (il = ilo; il <= nskypix; il = il + 1) {
+		if (skypix[index[il]] >= locut)
+		    break
+	    }
+	    for (ih = ihi; ih >= 1; ih = ih - 1) {
+		if (skypix[index[ih]] <= hicut)
+		    break
+	    }
+	    if (il == ilo && ih == ihi)
+		break
+
+	    # Reject pixels with optional region growing.
+	    if (rgrow > 0.0) {
+
+		# Reject low side pixels with region growing.
+	        do j = ilo, il - 1 
+		    nsky_reject = nsky_reject + ap_grow_regions (skypix, coords,
+		        wgt, nskypix, sky_zero, index[j], snx, sny, rgrow,
+			sumpx, sumsqpx, sumcbpx)
+
+		# Reject high side pixels with region growing.
+		do j = ih + 1, ihi
+		    nsky_reject = nsky_reject + ap_grow_regions (skypix, coords,
+		        wgt, nskypix, sky_zero, index[j], snx, sny, rgrow,
+			sumpx, sumsqpx, sumcbpx)
+
+		# Recompute the median.
+		# and the median pixel
+		nsky = nskypix - nsky_reject
+		med = apimed (wgt, index, il, ih, (nsky + 1) / 2)
+
+
+	    } else {
+
+		# Recompute the number of sky pixels, the number of rejected
+		# pixels and median pixel
+		nsky_reject = nsky_reject + (il - ilo) + (ihi - ih)
+		nsky = nskypix - nsky_reject
+		med = (ih + il) / 2
+
+		# Reject pixels on the low side.
+		do j = ilo, il - 1 {
+		    dsky = skypix[index[j]] - sky_zero
+		    sumpx = sumpx - dsky
+		    sumsqpx = sumsqpx - dsky ** 2
+		    sumcbpx = sumcbpx - dsky ** 3
+		}
+
+		# Reject pixels on the high side.
+		do j = ih + 1, ihi {
+		    dsky = skypix[index[j]] - sky_zero
+		    sumpx = sumpx - dsky
+		    sumsqpx = sumsqpx - dsky ** 2
+		    sumcbpx = sumcbpx - dsky ** 3
+		}
+	    }
+	    if (nsky <= 0)
+		break
+
+	    # Recompute the median, sigma and skew.
+	    medcut = nint (MEDCUT * real (nsky))
+	    sky_med = apwsmed (skypix, index, wgt, nskypix, med, medcut)
+	    sky_med = max (dmin, min (sky_med, dmax))
+	    call apmoments (sumpx, sumsqpx, sumcbpx, nsky, sky_zero,
+	        sky_mean, sky_sigma, sky_skew)
+	    sky_mean = max (dmin, min (sky_mean, dmax))
+	    #call printf (
+	    #"    mean=%g med=%g sigma=%g nsky=%d nrej=%d dmin=%g dmax=%g\n")
+	        #call pargr (sky_mean)
+	        #call pargr (sky_med)
+	        #call pargr (sky_sigma)
+	        #call pargi (nsky)
+	        #call pargi (nsky_reject)
+	        #call pargr (dmin)
+	        #call pargr (dmax)
+
+	    if (sky_sigma <= 0.0)
+		break
+	    ilo = il
+	    ihi = ih
+	}
+
+	# Return an appropriate error code.
+	if (nsky == 0 || nsky_reject == nskypix) {
+	    nsky = 0
+	    nsky_reject = nskypix
+	    sky_med = INDEFR
+	    sky_sigma = INDEFR
+	    sky_skew = INDEFR
+	    return (AP_NSKY_TOO_SMALL)
+	} else
+	    return (AP_OK)
+end