aboutsummaryrefslogtreecommitdiff
path: root/noao/digiphot/apphot/fitsky/apgauss.x
diff options
context:
space:
mode:
Diffstat (limited to 'noao/digiphot/apphot/fitsky/apgauss.x')
-rw-r--r--noao/digiphot/apphot/fitsky/apgauss.x296
1 files changed, 296 insertions, 0 deletions
diff --git a/noao/digiphot/apphot/fitsky/apgauss.x b/noao/digiphot/apphot/fitsky/apgauss.x
new file mode 100644
index 00000000..860301dc
--- /dev/null
+++ b/noao/digiphot/apphot/fitsky/apgauss.x
@@ -0,0 +1,296 @@
+include <mach.h>
+include <math/nlfit.h>
+include "../lib/fitsky.h"
+
+define TOL .001 # Fitting tolerance
+
+# AP_GAUSS -- Procedure to compute the peak, width and skew of the histogram
+# by fitting a skewed Gaussian function to the histogram.
+
+int procedure ap_gauss (skypix, coords, wgt, index, nskypix, snx, sny, maxfit,
+ k1, hwidth, binsize, smooth, losigma, hisigma, rgrow, maxiter,
+ sky_mode, sky_sigma, sky_skew, nsky, nsky_reject)
+
+real skypix[ARB] # array of unsorted sky pixels
+int coords[ARB] # array of coordinates for region growing
+real wgt[ARB] # weights for pixel rejection
+int index[ARB] # array of sort indices
+int nskypix # the number of sky pixels
+int snx, sny # the maximum dimensions of sky raster
+int maxfit # maximum number of iterations per fit
+real k1 # extent of the histogram in skysigma
+real hwidth # width of the histogram
+real binsize # the size of the histogram in sky sigma
+int smooth # smooth the histogram before fitting
+real losigma, hisigma # upper and lower sigma rejection criterion
+real rgrow # region growing radius in pixels
+int maxiter # maximum number of rejection cycles
+real sky_mode # computed sky value
+real sky_sigma # computed sigma of the sky pixels
+real sky_skew # skew of sky pixels
+int nsky # number of sky pixels used in fit
+int nsky_reject # number of sky pixels rejected
+
+double dsky, sumpx, sumsqpx, sumcbpx
+int i, j, nreject, nbins, nker, ier
+pointer sp, x, hgm, shgm, w
+real sky_mean, sky_msigma, dmin, dmax, hmin, hmax, dh, locut, hicut, cut
+real sky_zero
+int ap_grow_hist2(), aphigmr()
+real ap_asumr(), apmedr()
+
+begin
+ # Initialize.
+ nsky = nskypix
+ nsky_reject = 0
+ sky_mode = INDEFR
+ sky_sigma = INDEFR
+ sky_skew = INDEFR
+ if (nskypix <= 0)
+ return (AP_NOSKYAREA)
+
+ # Compute initial guess for sky statistics.
+ sky_zero = ap_asumr (skypix, index, nskypix) / nskypix
+ call ap_ialimr (skypix, index, nskypix, dmin, dmax)
+ call apfimoments (skypix, index, nskypix, sky_zero, sumpx, sumsqpx,
+ sumcbpx, sky_mean, sky_msigma, sky_skew)
+ sky_mean = apmedr (skypix, index, nskypix)
+ sky_mean = max (dmin, min (sky_mean, dmax))
+
+ # Compute the width and bin size of the sky histogram.
+ if (! IS_INDEFR(hwidth) && hwidth > 0.0) {
+ hmin = sky_mean - k1 * hwidth
+ hmax = sky_mean + k1 * hwidth
+ dh = binsize * hwidth
+ } else {
+ cut = min (sky_mean - dmin, dmax - sky_mean, k1 * sky_msigma)
+ hmin = sky_mean - cut
+ hmax = sky_mean + cut
+ dh = binsize * cut / k1
+ }
+
+ # Compute the number of histogram bins and width of smoothing kernel.
+ if (dh <= 0.0) {
+ nbins = 1
+ dh = 0.0
+ } else {
+ nbins = 2 * nint ((hmax - sky_mean) / dh) + 1
+ dh = (hmax - hmin) / (nbins - 1)
+ }
+
+ # Test for a valid histogram.
+ if (sky_msigma <= dh || dh <= 0.0 || k1 <= 0.0 || sky_msigma <= 0.0 ||
+ nbins < 4) {
+ sky_mode = sky_mean
+ sky_sigma = 0.0
+ sky_skew = 0.0
+ return (AP_NOHISTOGRAM)
+ }
+
+ # Allocate temporary working space.
+ call smark (sp)
+ call salloc (x, nbins, TY_REAL)
+ call salloc (hgm, nbins, TY_REAL)
+ call salloc (shgm, nbins, TY_REAL)
+ call salloc (w, nbins, TY_REAL)
+
+ # Compute the x array.
+ do i = 1, nbins
+ Memr[x+i-1] = i
+ call amapr (Memr[x], Memr[x], nbins, 1.0, real (nbins),
+ hmin + 0.5 * dh, hmax + 0.5 * dh)
+
+ # Accumulate the histogram.
+ call aclrr (Memr[hgm], nbins)
+ nsky_reject = nsky_reject + aphigmr (skypix, wgt, index, nskypix,
+ Memr[hgm], nbins, hmin, hmax)
+ nsky = nskypix - nsky_reject
+
+ # Do the initial rejection.
+ if (nsky_reject > 0) {
+ do i = 1, nskypix {
+ if (wgt[index[i]] <= 0.0) {
+ dsky = skypix[index[i]] - sky_zero
+ sumpx = sumpx - dsky
+ sumsqpx = sumsqpx - dsky ** 2
+ sumcbpx = sumcbpx - dsky ** 3
+ }
+ }
+ call apmoments (sumpx, sumsqpx, sumcbpx, nsky, sky_zero,
+ sky_mean, sky_msigma, sky_skew)
+ }
+
+ # Find the mode, sigma and skew of the histogram.
+ if (smooth == YES) {
+ nker = max (1, nint (sky_msigma / dh))
+ #call ap_lucy_smooth (Memr[hgm], Memr[shgm], nbins, nker, 2)
+ call ap_bsmooth (Memr[hgm], Memr[shgm], nbins, nker, 2)
+ call ap_hist_mode (Memr[x], Memr[shgm], Memr[w], nbins,
+ sky_mode, sky_sigma, sky_skew, maxfit, TOL, ier)
+ } else
+ call ap_hist_mode (Memr[x], Memr[hgm], Memr[w], nbins,
+ sky_mode, sky_sigma, sky_skew, maxfit, TOL, ier)
+ sky_mode = max (dmin, min (sky_mode, dmax))
+ if (ier != AP_OK) {
+ call sfree (sp)
+ return (ier)
+ }
+ if ((IS_INDEFR(losigma) && IS_INDEFR(hisigma)) || sky_sigma <= dh ||
+ maxiter < 1) {
+ call sfree (sp)
+ return (AP_OK)
+ }
+
+ # Fit histogram with pixel rejection and optional region growing.
+ do i = 1, maxiter {
+
+ # Compute the new rejection limits.
+ if (IS_INDEFR(losigma))
+ locut = -MAX_REAL
+ else
+ locut = sky_mode - losigma * sky_msigma
+ if (IS_INDEFR(hisigma))
+ hicut = MAX_REAL
+ else
+ hicut = sky_mode + hisigma * sky_msigma
+
+ # Detect and reject pixels.
+ nreject = 0
+ do j = 1, nskypix {
+ if (skypix[index[j]] >= locut && skypix[index[j]] <= hicut)
+ next
+ if (rgrow > 0.0)
+ nreject = nreject + ap_grow_hist2 (skypix, coords,
+ wgt, nskypix, sky_zero, index[j], snx, sny, Memr[hgm],
+ nbins, hmin, hmax, rgrow, sumpx, sumsqpx, sumcbpx)
+ else if (wgt[index[j]] > 0.0) {
+ call ap_hgmsub2 (Memr[hgm], nbins, hmin, hmax,
+ skypix[index[j]], sky_zero, sumpx, sumsqpx, sumcbpx)
+ wgt[index[j]] = 0.0
+ nreject = nreject + 1
+ }
+ }
+ if (nreject == 0)
+ break
+ nsky_reject = nsky_reject + nreject
+ nsky = nskypix - nsky_reject
+ if (nsky <= 0)
+ break
+ call apmoments (sumpx, sumsqpx, sumcbpx, nsky, sky_zero,
+ sky_mean, sky_msigma, sky_skew)
+
+ # Recompute mean, mode, sigma and skew.
+ if (smooth == YES) {
+ nker = max (1, nint (sky_msigma / dh))
+ #call ap_lucy_smooth (Memr[hgm], Memr[shgm], nbins, nker, 2)
+ call ap_bsmooth (Memr[hgm], Memr[shgm], nbins, nker, 2)
+ call ap_hist_mode (Memr[x], Memr[shgm], Memr[w], nbins,
+ sky_mode, sky_sigma, sky_skew, maxfit, TOL, ier)
+ } else
+ call ap_hist_mode (Memr[x], Memr[hgm], Memr[w], nbins,
+ sky_mode, sky_sigma, sky_skew, maxfit, TOL, ier)
+ sky_mode = max (dmin, min (sky_mode, dmax))
+ if (ier != AP_OK)
+ break
+ if (sky_sigma <= dh)
+ break
+ }
+
+ # Return appropriate error code.
+ call sfree (sp)
+ if (nsky == 0 || nsky_reject == nskypix) {
+ nsky = 0
+ nsky_reject = nskypix
+ sky_mode = INDEFR
+ sky_sigma = INDEFR
+ sky_skew = INDEFR
+ return (AP_NSKY_TOO_SMALL)
+ } else if (ier == AP_NSKY_TOO_SMALL) {
+ return (AP_NSKY_TOO_SMALL)
+ } else if (ier != AP_OK) {
+ return (ier)
+ } else
+ return (AP_OK)
+end
+
+
+# AP_HIST_MODE -- Procedure to fit the skewed Gaussian function to the histogram
+# of the sky pixels.
+
+procedure ap_hist_mode (x, hgm, w, nbins, sky_mode, sky_sigma, sky_skew,
+ maxiter, tol, ier)
+
+real x[ARB] # x array
+real hgm[ARB] # histogram
+real w[ARB] # weights
+int nbins # number of bins
+real sky_mode # sky value
+real sky_sigma # sky sigma
+real sky_skew # sky skew
+int maxiter # max number of iterations
+real tol # tolerances
+int ier # error code
+
+extern gausskew, dgausskew
+int i, imin, imax, np, fier
+pointer sp, list, fit, nl
+real p[4], dp[4], dummy1
+int locpr()
+
+begin
+ # Allocate working memory.
+ call smark (sp)
+ call salloc (list, 4, TY_INT)
+ call salloc (fit, nbins, TY_REAL)
+
+ # Initialize.
+ do i = 1, 4
+ Memi[list+i-1] = i
+
+ # Compute initial guesses for the parameters.
+ call ap_alimr (hgm, nbins, dummy1, p[1], imin, imax)
+ p[2] = x[imax]
+ #p[3] = max (sky_sigma ** 2, abs (x[2] - x[1]) ** 2)
+ p[3] = abs ((x[nbins] - x[1]) / 6.0) ** 2
+ p[4] = 0.0
+ np = 4
+
+ # Fit the histogram.
+ call nlinitr (nl, locpr (gausskew), locpr (dgausskew), p, dp, np,
+ Memi[list], 4, tol, maxiter)
+ call nlfitr (nl, x, hgm, w, nbins, 1, WTS_UNIFORM, fier)
+ call nlvectorr (nl, x, Memr[fit], nbins, 1)
+ call nlpgetr (nl, p, np)
+ call nlfreer (nl)
+
+ call sfree (sp)
+
+ # Return the appropriate error code.
+ ier = AP_OK
+ if (fier == NO_DEG_FREEDOM) {
+ sky_mode = INDEFR
+ sky_sigma = INDEFR
+ sky_skew = INDEFR
+ ier = AP_NSKY_TOO_SMALL
+ } else {
+ if (fier == SINGULAR)
+ ier = AP_SKY_SINGULAR
+ else if (fier == NOT_DONE)
+ ier = AP_SKY_NOCONVERGE
+ if (p[2] < x[1] || p[2] > x[nbins]) {
+ sky_mode = x[imax]
+ ier = AP_BADPARAMS
+ } else
+ sky_mode = p[2]
+ if (p[3] <= 0.0) {
+ sky_sigma = 0.0
+ sky_skew = 0.0
+ ier = AP_BADPARAMS
+ } else {
+ sky_sigma = sqrt (p[3])
+ sky_skew = 1.743875281 * abs (p[4]) ** (1.0 / 3.0) * sky_sigma
+ if (p[4] < 0.0)
+ sky_skew = - sky_skew
+ }
+ }
+end
generated by cgit (git 2.34.1) at 2025-09-20 07:18:26 -0400