Initial commit

1 files changed, 136 insertions, 0 deletions

diff --git a/noao/digiphot/apphot/fitsky/apgrowhist.x b/noao/digiphot/apphot/fitsky/apgrowhist.x
new file mode 100644
index 00000000..6f01a789
--- /dev/null
+++ b/noao/digiphot/apphot/fitsky/apgrowhist.x
@@ -0,0 +1,136 @@
+# AP_GROW_HIST -- Procedure to reject pixels with region growing.
+
+int procedure ap_grow_hist (skypix, coords, wgt, nskypix, index, snx, sny, hgm,
+    nbins, z1, z2, rgrow)
+
+real	skypix[ARB]		# sky pixels
+int	coords[ARB]		# sky pixel coordinates
+real	wgt[ARB]		# array of weights
+int	nskypix			# number of sky pixels
+int	index			# index of pixel to be rejected
+int	snx, sny		# size of sky subraster
+real	hgm[ARB]		# histogram
+int	nbins			# number of bins
+real	z1, z2			# value of first and last histogram bin
+real	rgrow			# region growing radius
+
+int	j, k, ixc, iyc, ymin, ymax, xmin, xmax, nreject, cstart, c, bin
+real	dh, r2, rgrow2, d
+
+begin
+	# Find the x and y coordinates of the pixel to be rejected.
+	ixc = mod (coords[index], snx)
+	if (ixc == 0)
+	    ixc = snx
+	iyc = (coords[index] - ixc) / snx + 1
+
+	# Find the coordinate space to be used for regions growing.
+	ymin = max (1, int (iyc - rgrow))
+	ymax = min (sny, int (iyc + rgrow))
+	xmin = max (1, int (ixc - rgrow))
+	xmax = min (snx, int (ixc + rgrow))
+	if (ymin <= iyc)
+	    cstart = min (nskypix, max (1, index - int (rgrow) + snx *
+	        (ymin - iyc)))
+	else
+	    cstart = index
+
+	# Perform the region growing.
+	dh = real (nbins - 1) / (z2 - z1)
+	rgrow2 = rgrow ** 2
+	nreject = 0
+	do j = ymin, ymax {
+	    d = rgrow2 - (j - iyc) ** 2
+	    if (d <= 0.0)
+		d = 0.0
+	    else
+		d = sqrt (d)
+	    do k = max (xmin, int (ixc - d)), min (xmax, int (ixc + d)) {
+		c = k + (j - 1) * snx
+		while (coords[cstart] < c && cstart < nskypix)
+		    cstart = cstart + 1
+		r2 = (k - ixc) ** 2 + (j - iyc) ** 2
+		if (r2 <= rgrow2 && c == coords[cstart] && wgt[cstart] > 0.0) {
+		    nreject = nreject + 1
+		    wgt[cstart] = 0.0
+		    if (skypix[cstart] >= z1 && skypix[cstart] <= z2) {
+		        bin = int ((skypix[cstart] - z1) * dh) + 1
+		        hgm[bin] = hgm[bin] - 1.0
+		    }
+		}
+	    }
+	}
+
+	return (nreject)
+end
+
+
+# AP_GROW_HIST2 -- Procedure to reject pixels with region growing.
+
+int procedure ap_grow_hist2 (skypix, coords, wgt, nskypix, sky_zero, index,
+    snx, sny, hgm, nbins, z1, z2, rgrow, sumpx, sumsqpx, sumcbpx)
+
+real	skypix[ARB]		# sky pixels
+int	coords[ARB]		# coordinates
+real	wgt[ARB]		# array of weights
+int	nskypix			# number of sky pixels
+real	sky_zero		# the sky zero point for moment analysis
+int	index			# index of pixel to be rejected
+int	snx, sny		# size of sky subraster
+real	hgm[ARB]		# histogram
+int	nbins			# number of bins
+real	z1, z2			# value of first and last histogram bin
+real	rgrow			# region growing radius
+double	sumpx			# sum of sky values
+double	sumsqpx			# sum of sky values squared
+double	sumcbpx			# sum of sky values cubed
+
+double	dsky
+int	j, k, ixc, iyc, ymin, ymax, xmin, xmax, nreject, cstart, c, bin
+real	dh, r2, rgrow2, d
+
+begin
+	# Find the coordinates of the region growing center.
+	ixc = mod (coords[index], snx)
+	if (ixc == 0)
+	    ixc = snx
+	iyc = (coords[index] - ixc) / snx + 1
+
+	ymin = max (1, int (iyc - rgrow))
+	ymax = min (sny, int (iyc + rgrow))
+	xmin = max (1, int (ixc - rgrow))
+	xmax = min (snx, int (ixc + rgrow))
+	dh = real (nbins - 1) / (z2 - z1)
+	if (ymin <= iyc)
+	    cstart = min (nskypix, max (1, index - int (rgrow) + snx *
+	        (ymin - iyc)))
+	else
+	    cstart = index
+
+	# Perform the region growing.
+	nreject = 0
+	rgrow2 = rgrow ** 2
+	do j = ymin, ymax {
+	    d = sqrt (rgrow2 - (j - iyc) ** 2)
+	    do k = max (xmin, int (ixc - d)), min (xmax, int (ixc + d)) {
+		c = k + (j - 1) * snx
+		while (coords[cstart] < c && cstart < nskypix)
+		    cstart = cstart + 1
+		r2 = (k - ixc) ** 2 + (j - iyc) ** 2
+		if (r2 <= rgrow2 && c == coords[cstart] && wgt[cstart] > 0.0) {
+		    nreject = nreject + 1
+		    wgt[cstart] = 0.0
+		    dsky = skypix[cstart] - sky_zero
+		    sumpx = sumpx - dsky
+		    sumsqpx = sumsqpx - dsky ** 2
+		    sumcbpx = sumcbpx - dsky ** 3
+		    if (skypix[cstart] >= z1 && skypix[cstart] <= z2) {
+		        bin = int ((skypix[cstart] - z1) * dh) + 1
+		        hgm[bin] = hgm[bin] - 1.0
+		    }
+		}
+	    }
+	}
+
+	return (nreject)
+end