diff options
Diffstat (limited to 'noao/digiphot/apphot/fitsky/apcentroid.x')
| -rw-r--r-- | noao/digiphot/apphot/fitsky/apcentroid.x | 244 |
1 files changed, 244 insertions, 0 deletions
diff --git a/noao/digiphot/apphot/fitsky/apcentroid.x b/noao/digiphot/apphot/fitsky/apcentroid.x new file mode 100644 index 00000000..310db85d --- /dev/null +++ b/noao/digiphot/apphot/fitsky/apcentroid.x @@ -0,0 +1,244 @@ +include <mach.h> +include "../lib/fitsky.h" + +# AP_CENTROID -- Procedure to find the mode, width and skew of the sky +# distribution by computing the moments of the histogram. + +int procedure ap_centroid (skypix, coords, wgt, index, nskypix, snx, sny, k1, + hwidth, binsize, smooth, losigma, hisigma, rgrow, maxiter, sky_mode, + sky_sigma, sky_skew, nsky, nsky_reject) + +real skypix[ARB] # array of sky pixels +int coords[ARB] # array of coordinates for regions growing +real wgt[ARB] # array of weights for rejection +int index[ARB] # array of sorted indices +int nskypix # the number of sky pixels +int snx, sny # the maximum dimensions of sky raster +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 k-sigma rejection limits +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 nreject, nbins, nker, ier, i, j +pointer sp, hgm, shgm +real dmin, dmax, sky_mean, hmin, hmax, dh, locut, hicut, cut +real sky_zero +int ap_grow_hist2(), aphigmr() +real ap_asumr(), apmedr() + +begin + # Intialize. + nsky = nskypix + nsky_reject = 0 + sky_mode = INDEFR + sky_sigma = INDEFR + sky_skew = INDEFR + if (nskypix <= 0) + return (AP_NOSKYAREA) + + # Compute the histogram width and binsize. + 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_sigma, sky_skew) + sky_mean = apmedr (skypix, index, nskypix) + sky_mean = max (dmin, min (sky_mean, dmax)) + 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_sigma) + hmin = sky_mean - cut + hmax = sky_mean + cut + dh = binsize * cut / k1 + } + + # Compute the number of histogram bins and the histogram resolution. + if (dh <= 0.0) { + nbins = 1 + dh = 0.0 + } else { + nbins = 2 * nint ((hmax - sky_mean) / dh) + 1 + dh = (hmax - hmin) / (nbins - 1) + } + + # Check for a valid histogram. + if (nbins < 2 || k1 <= 0.0 || sky_sigma < dh || dh <= 0.0 || + sky_sigma <= 0.0) { + sky_mode = sky_mean + sky_sigma = 0.0 + sky_skew = 0.0 + return (AP_NOHISTOGRAM) + } + + # Allocate temporary space. + call smark (sp) + call salloc (hgm, nbins, TY_REAL) + call salloc (shgm, nbins, TY_REAL) + + # 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_sigma, sky_skew) + } + + # Fit the mode, sigma an skew of the histogram. + if (smooth == YES) { + nker = max (1, nint (sky_sigma / dh)) + #call ap_lucy_smooth (Memr[hgm], Memr[shgm], nbins, nker, 2) + call ap_bsmooth (Memr[hgm], Memr[shgm], nbins, nker, 2) + call ap_imode (Memr[shgm], nbins, hmin, hmax, YES, sky_mode, ier) + } else + call ap_imode (Memr[hgm], nbins, hmin, hmax, NO, sky_mode, ier) + sky_mode = max (dmin, min (sky_mode, dmax)) + 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 new histogram limits. + if (IS_INDEFR(losigma)) + locut = -MAX_REAL + else + locut = sky_mode - losigma * sky_sigma + if (IS_INDEFR(hisigma)) + hicut = MAX_REAL + else + hicut = sky_mode + hisigma * sky_sigma + + # 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 + + # Recompute the moments. + call apmoments (sumpx, sumsqpx, sumcbpx, nsky, sky_zero, sky_mean, + sky_sigma, sky_skew) + + # Recompute the histogram peak. + if (smooth == YES) { + nker = max (1, nint (sky_sigma / dh)) + #call ap_lucy_smooth (Memr[hgm], Memr[shgm], nbins, nker, 2) + call ap_bsmooth (Memr[hgm], Memr[shgm], nbins, nker, 2) + call ap_imode (Memr[shgm], nbins, hmin, hmax, YES, sky_mode, + ier) + } else + call ap_imode (Memr[hgm], nbins, hmin, hmax, NO, sky_mode, ier) + sky_mode = max (dmin, min (sky_mode, dmax)) + if (sky_sigma <= dh || ier != AP_OK) + break + } + + # Return the error codes. + call sfree (sp) + if (nsky == 0 || nsky_reject == nskypix || ier == AP_NOHISTOGRAM) { + nsky = 0 + nsky_reject = nskypix + sky_mode = INDEFR + sky_sigma = INDEFR + sky_skew = INDEFR + return (AP_NSKY_TOO_SMALL) + } else + return (AP_OK) +end + + +# AP_IMODE -- Procedure to compute the 1st, 2nd and third moments of the +# histogram of sky pixels. + +procedure ap_imode (hgm, nbins, z1, z2, smooth, sky_mode, ier) + +real hgm[ARB] # histogram +int nbins # number of bins +real z1, z2 # min and max of the histogram +int smooth # is the histogram smoothed +real sky_mode # mode of histogram +int ier # error code + +int i, noccup +double sumi, sumix, x, dh +real hmean, dz +real asumr() + +begin + # Initialize the sums. + sumi = 0.0 + sumix = 0.0 + + # Compute a continuum level. + if (smooth == NO) + hmean = asumr (hgm, nbins) / nbins + else { + call alimr (hgm, nbins, dz, hmean) + hmean = 2.0 * hmean / 3.0 + } + + # Accumulate the sums. + noccup = 1 + dz = (z2 - z1) / (nbins - 1) + x = z1 + 0.5 * dz + do i = 1, nbins { + dh = hgm[i] - hmean + if (dh > 0.0d0) { + sumi = sumi + dh + sumix = sumix + dh * x + noccup = noccup + 1 + } + x = x + dz + } + + # Compute the sky mode, sigma and skew. + if (sumi > 0.0) { + sky_mode = sumix / sumi + ier = AP_OK + } else { + sky_mode = INDEFR + ier = AP_NOHISTOGRAM + } +end |