Initial commit

1 files changed, 64 insertions, 0 deletions

diff --git a/noao/digiphot/apphot/fitsky/aprgrow.x b/noao/digiphot/apphot/fitsky/aprgrow.x
new file mode 100644
index 00000000..902b91d6
--- /dev/null
+++ b/noao/digiphot/apphot/fitsky/aprgrow.x
@@ -0,0 +1,64 @@
+# AP_GROW_REGIONS -- Perform region growing around a rejected pixel.
+
+int procedure ap_grow_regions (skypix, coords, wgt, nskypix, sky_zero,
+	index, snx, sny, rgrow, sumpx, sumsqpx, sumcbpx)
+
+real	skypix[ARB]		# sky pixels
+int	coords[ARB]		# sky cordinates
+real	wgt[ARB]		# weights
+int	nskypix			# total number of sky pixels
+real	sky_zero		# sky zero point for moment analysis
+int	index			# index of pixel to be rejected
+int	snx, sny		# size of sky subraster
+real	rgrow			# region growing radius
+double	sumpx, sumsqpx, sumcbpx	# sum and sum of squares of sky pixels
+
+double	dsky
+int	j, k, ixc, iyc, xmin, xmax, ymin, ymax, cstart, c, nreject
+real	rgrow2, r2, d
+
+begin
+	# Find the center of the region to be rejected.
+	ixc = mod (coords[index], snx)
+	if (ixc == 0)
+	    ixc = snx
+	iyc = (coords[index] - ixc) / snx + 1
+
+	# Define the region to be searched.
+	rgrow2 = rgrow ** 2
+	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
+
+	# Reject the pixels.
+	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) {
+		    dsky = skypix[cstart] - sky_zero
+		    sumpx = sumpx - dsky
+		    sumsqpx = sumsqpx - dsky ** 2
+		    sumcbpx = sumcbpx - dsky ** 3
+		    nreject = nreject + 1
+		    wgt[cstart] = 0.0
+		}
+	    }
+	}
+
+	return (nreject)
+end