Initial commit

1 files changed, 393 insertions, 0 deletions

diff --git a/noao/nproto/ace/skyfit.x b/noao/nproto/ace/skyfit.x
new file mode 100644
index 00000000..0b295e8e
--- /dev/null
+++ b/noao/nproto/ace/skyfit.x
@@ -0,0 +1,393 @@
+include	<imhdr.h>
+include	<math/curfit.h>
+include	<math/gsurfit.h>
+include	"skyfit.h"
+
+
+# SKY_FIT -- Fit sky surface.
+#
+# Compute a sky and/or sky sigma surface fit using a subset of the input
+# lines.  the input sky and sky sigma pointers are NULL.  The initial data
+# for the surface fit is measured at a subset of lines with any masked
+# pixels excluded.  Objects are removed by fitting a 1D curve to each line,
+# rejection points with large residuals and iterating until only sky is left.
+# The sky points are then accumulated for a 2D surface fit and the residuals
+# are added to a histogram.  The absolute deviations, scaled by 0.7979 to
+# convert to an gausian sigma, are accumulated for a sky sigma surface fit.
+# After all the sample lines are accumulated the surface fits are computed.
+# The histogram of residuals is then fit by a gaussian to estimate an
+# offset from the sky fit to the sky mode caused by unrejected object light.
+# The offset is applied to the sky surface.
+
+procedure sky_fit (par, dosky, dosig, im, bpm, expmap, skyname, signame,
+	skymap, sigmap, logfd)
+
+pointer	par			#U Sky parameters
+bool	dosky			#I Compute sky
+bool	dosig			#I Compute sigma
+pointer	im			#I Input image
+pointer	bpm			#I Input mask
+pointer	expmap			#I Exposure map
+char	skyname[ARB]		#I Sky map name
+char	signame[ARB]		#I Sigma map name
+pointer	skymap			#U Sky map
+pointer	sigmap			#U Sigma map
+int	logfd			#I Verbose?
+
+# Parameters
+real	step			# Line sample step
+int	lmin			# Minimum number of lines to fit
+int	func1d			# 1D fitting function
+int	func2d			# 2D fitting function
+int	xorder			# Sky fitting x order
+int	yorder			# Sky fitting y order
+int	xterms			# Sky fitting cross terms
+int	blk1d			# Block average
+real	hclip			# Sky fitting high sigma clip
+real	lclip			# Sky fitting low sigma clip
+int	niter			# Number of clipping iterations
+
+int	l1, l2
+int	i, j, c, l, n, nc, nl, nskyblk, ier
+real	res, sigma
+pointer	sp, x, y, z, r, a, x1, w1, w2, skydata, sigdata, expdata, w, ptr
+pointer	cvsky, cvsig, gssky, gssig
+
+pointer	imgl2r(), imgl2i(), map_opengs(), map_glr()
+bool	im_pmlne2()
+real	amedr()
+errchk	map_opengs, map_glr
+
+begin
+	if (!(dosky||dosig))
+	    return
+
+	# Set parameters.
+	if (par == NULL)
+	    call skf_pars ("open", "", par)
+	step = SKF_STEP(par)
+	lmin = SKF_LMIN(par)
+	xorder = SKF_XORDER(par)
+	yorder = SKF_YORDER(par)
+	xterms = SKF_XTERMS(par)
+	blk1d = SKF_BLK1D(par)
+	hclip = SKF_HCLIP(par)
+	lclip = SKF_LCLIP(par)
+	func1d = SKF_FUNC1D(par)
+	func2d = SKF_FUNC2D(par)
+	niter = SKF_NITER(par)
+
+	nc = IM_LEN(im,1)
+	nl = IM_LEN(im,2)
+	l1 = 1 + step / 2
+	l2 = nl - step / 2
+	step = real (l2-l1) / max (nint((l2-l1)/step),xorder+2,lmin)
+
+	if (logfd != NULL) {
+	    if (dosky && dosig)
+		call fprintf (logfd,
+		    "    Determine sky and sigma by surface fits:\n")
+	    else if (dosky)
+		call fprintf (logfd, "    Determine sky by surface fit:\n")
+	    else
+		call fprintf (logfd, "    Determine sigma by surface fit:\n")
+	    call fprintf (logfd,
+		"      start line = %d, end line = %d, step = %.1f\n")
+		call pargi (l1)
+		call pargi (l2)
+		call pargr (step)
+	    call fprintf (logfd,
+		"      xorder = %d, yorder = %d, xterms = %s\n")
+		call pargi (xorder)
+		call pargi (yorder)
+		switch (xterms) {
+		case GS_XNONE:
+		    call pargstr ("none")
+		case GS_XFULL:
+		    call pargstr ("full")
+		case GS_XHALF:
+		    call pargstr ("half")
+		}
+	    call fprintf (logfd, "      hclip = %g, lclip = %g\n")
+		call pargr (hclip)
+		call pargr (lclip)
+	}
+
+	# Allocate memory and initialize.
+	call smark (sp)
+	call salloc (x1, nc, TY_REAL)
+	call salloc (w1, nc, TY_REAL)
+	call salloc (w2, nc, TY_REAL)
+
+	nskyblk = nc / blk1d
+	call salloc (x, nskyblk, TY_REAL)
+	call salloc (y, nskyblk, TY_REAL)
+	call salloc (z, nskyblk, TY_REAL)
+	call salloc (r, nskyblk, TY_REAL)
+	call salloc (a, nskyblk, TY_REAL)
+	call salloc (skydata, nskyblk, TY_REAL)
+	call salloc (sigdata, nskyblk, TY_REAL)
+	if (expmap != NULL)
+	    call salloc (expdata, nskyblk, TY_REAL)
+
+	do c = 1, nc
+	    Memr[x1+c-1] = c
+	call amovkr (1., Memr[w1], nc)
+
+	# Initialize the 1D and 2D fitting pointers as needed.
+	if (dosky) {
+	    call cvinit (cvsky, func1d, xorder, Memr[x1],
+		Memr[x1+nc-1])
+	    call gsinit (gssky, func2d, xorder, yorder,
+		xterms, 1., real(nc), 1., real(nl))
+	}
+	if (dosig) {
+	    call cvinit (cvsig, CHEBYSHEV, 1, Memr[x1], Memr[x1+nc-1])
+	    call gsinit (gssig, GS_CHEBYSHEV, 1, 1, xterms,
+		1., real(nc), 1., real(nl))
+	}
+
+	# For each sample line find sky points by 1D fitting and sigma
+	# rejection and then accumulate 2D surface fitting points.
+	do j = 0, ARB {
+	    l = nint (l1 + j * step)
+	    if (l > l2)
+		break
+
+	    # Get input data and block average.
+	    if (bpm == NULL)
+		w = w1
+	    else if (!im_pmlne2 (bpm, l))
+		w = w1
+	    else {
+		w = imgl2i (bpm, l)
+		n = nc
+		do c = 0, nc-1 {
+		    if (Memi[w+c] != 0) {
+			Memr[w2+c] = 0
+			n = n - 1
+		    } else
+			Memr[w2+c] = Memr[w1+c]
+		}
+		w = w2
+		if (n < 10)
+		    next
+	    }
+
+	    # Block average.
+	    if (skymap != NULL) {
+		ptr = map_glr (skymap, l, READ_ONLY)
+		call blkavg1 (Memr[ptr], Memr[w], nc, Memr[skydata],
+		    nskyblk, blk1d)
+	    }
+	    if (expmap != NULL) {
+		ptr = map_glr (expmap, l, READ_ONLY)
+		call blkavg1 (Memr[ptr], Memr[w], nc, Memr[expdata],
+		    nskyblk, blk1d)
+	    }
+	    if (sigmap != NULL) {
+		ptr = map_glr (sigmap, l, READ_ONLY)
+		call blkavg1 (Memr[ptr], Memr[w], nc, Memr[sigdata],
+		    nskyblk, blk1d)
+		call adivkr (Memr[sigdata], sqrt(real(blk1d)), Memr[sigdata],
+		    nskyblk)
+		if (expmap != NULL)
+		    call expsigma (Memr[sigdata], Memr[expdata], nskyblk, 0)
+	    }
+	    call blkavg (Memr[x1], Memr[imgl2r(im,l)], Memr[w], nc,
+		Memr[x], Memr[z], Memr[w2], nskyblk, blk1d)
+	    w = w2
+
+	    # Iterate using line fitting.
+	    do i = 1, niter {
+
+		# Fit sky.
+		if (dosky) {
+		    call cvfit (cvsky, Memr[x], Memr[z], Memr[w], nskyblk,
+			WTS_USER, ier)
+		    if (ier == NO_DEG_FREEDOM)
+			call error (1, "Fitting error")
+		    call cvvector (cvsky, Memr[x], Memr[skydata], nskyblk)
+		}
+
+		# Compute residuals.
+		call asubr (Memr[z], Memr[skydata], Memr[r], nskyblk)
+
+		# Fit sky sigma.
+		if (dosig) {
+		    do c = 0, nskyblk-1
+			Memr[a+c] = abs(Memr[r+c]) / 0.7979
+		    if (expmap != NULL)
+			call expsigma (Memr[a], Memr[expdata], nskyblk, 1)
+		    if (i == 1)
+			call amovkr (amedr(Memr[a],nskyblk), Memr[sigdata],
+			    nskyblk)
+		    else {
+			call cvfit (cvsig, Memr[x], Memr[a], Memr[w], nskyblk,
+			    WTS_USER, ier)
+			if (ier == NO_DEG_FREEDOM)
+			    call error (1, "Fitting error")
+			call cvvector (cvsig, Memr[x], Memr[sigdata], nskyblk)
+		    }
+		    if (expmap != NULL)
+			call expsigma (Memr[sigdata], Memr[expdata], nskyblk, 0)
+		}
+
+		# Reject deviant points.
+		n = 0
+		do c = 0, nskyblk-1 {
+		    if (Memr[w+c] == 0.)
+			next
+		    res = Memr[r+c]
+		    sigma = Memr[sigdata+c]
+		    if (res > hclip * sigma || res < -lclip * sigma) {
+			Memr[w+c] = 0.
+			n = n + 1
+		    }
+		}
+		if (n == 0) {
+		    if (i == 1 && dosig) {
+			call cvfit (cvsig, Memr[x], Memr[a], Memr[w], nskyblk,
+			    WTS_USER, ier)
+			if (ier == NO_DEG_FREEDOM)
+			    call error (1, "Fitting error")
+		    }
+		    break
+		}
+	    }
+
+	    # Accumulate the sky data for the line.
+	    call amovkr (real(l), Memr[y], nskyblk)
+	    if (dosky && dosig) {
+		call amulkr (Memr[a], sqrt(real(blk1d)), Memr[a], nskyblk)
+		call gsacpts (gssky, Memr[x], Memr[y], Memr[z], Memr[w],
+		    nskyblk, WTS_USER)
+		call gsacpts (gssig, Memr[x], Memr[y], Memr[a], Memr[w],
+		    nskyblk, WTS_USER)
+	    } else if (dosky) {
+		call gsacpts (gssky, Memr[x], Memr[y], Memr[z], Memr[w],
+		    nskyblk, WTS_USER)
+	    } else {
+		call amulkr (Memr[a], sqrt(real(blk1d)), Memr[a], nskyblk)
+		call gsacpts (gssig, Memr[x], Memr[y], Memr[a],
+		    Memr[w], nskyblk, WTS_USER)
+	    }
+	}
+
+	# Compute the surface fits, store in header, and set output pointers.
+	if (dosky) {
+	    if (skymap != NULL)
+	        call map_close (skymap)
+	    call cvfree (cvsky)
+	    call gssolve (gssky, ier)
+	    if (ier == NO_DEG_FREEDOM)
+		call error (1, "Fitting error")
+	    if (skyname[1] != EOS)
+		call mgs_pgs (im, skyname, gssky)
+	    skydata = map_opengs (gssky, im); skymap = skydata
+	}
+	if (dosig) {
+	    if (sigmap != NULL)
+	        call map_close (sigmap)
+	    call cvfree (cvsig)
+	    call gssolve (gssig, ier)
+	    if (ier == NO_DEG_FREEDOM)
+		call error (1, "Fitting error")
+	    if (signame[1] != EOS)
+		call mgs_pgs (im, signame, gssig)
+	    sigdata = map_opengs (gssig, im); sigmap = sigdata
+	}
+
+	call sfree (sp)
+end
+
+
+procedure blkavg (xin, yin, win, nin, xout, yout, wout, nout, blksize)
+
+real	xin[nin]		#I Input values
+real	yin[nin]		#I Input values
+real	win[nin]		#I Input weights
+int	nin			#I Number of input values
+real	xout[nout]		#O Output values
+real	yout[nout]		#O Output values
+real	wout[nout]		#O Output weights
+int	nout			#O Number of output values
+int	blksize			#I Block size
+
+int	i, j, n, imax
+real	xavg, yavg, wsum, w
+
+begin
+	if (blksize == 1) {
+	    nout = nin
+	    call amovr (xin, xout, nout)
+	    call amovr (yin, yout, nout)
+	    call amovr (win, wout, nout)
+	    return
+	}
+
+	n = blksize
+	imax = nin - 2 * blksize + 1
+	nout = 0
+	for (i=1; i<=nin; ) {
+	    if (i > imax)
+		n = nin - i + 1
+	    xavg = 0.
+	    yavg = 0.
+	    wsum = 0.
+	    do j = 1, n {
+		w = win[i]
+		xavg = xavg + w * xin[i] 
+		yavg = yavg + w * yin[i] 
+		wsum = wsum + w
+		i = i + 1
+	    }
+	    if (wsum > 0.) {
+		nout = nout + 1
+		xout[nout] = xavg / wsum
+		yout[nout] = yavg / wsum
+		wout[nout] = wsum
+	    }
+	}
+end
+
+
+procedure blkavg1 (in, win, nin, out, nout, blksize)
+
+real	in[nin]			#I Input values
+real	win[nin]		#I Input weights
+int	nin			#I Number of input values
+real	out[nout]		#O Output values
+int	nout			#O Number of output values
+int	blksize			#I Block size
+
+int	i, j, n, imax
+real	avg, wsum, w
+
+begin
+	if (blksize == 1) {
+	    nout = nin
+	    call amovr (in, out, nout)
+	    return
+	}
+
+	n = blksize
+	imax = nin - 2 * blksize + 1
+	nout = 0
+	for (i=1; i<=nin; ) {
+	    if (i > imax)
+		n = nin - i + 1
+	    avg = 0.
+	    wsum = 0.
+	    do j = 1, n {
+		w = win[i]
+		avg = avg + w * in[i] 
+		wsum = wsum + w
+		i = i + 1
+	    }
+	    if (wsum > 0.) {
+		nout = nout + 1
+		out[nout] = avg / wsum
+	    }
+	}
+end