aboutsummaryrefslogtreecommitdiff
path: root/noao/digiphot/apphot/fitsky/apfitsky.x
diff options
context:
space:
mode:
Diffstat (limited to 'noao/digiphot/apphot/fitsky/apfitsky.x')
-rw-r--r--noao/digiphot/apphot/fitsky/apfitsky.x452
1 files changed, 452 insertions, 0 deletions
diff --git a/noao/digiphot/apphot/fitsky/apfitsky.x b/noao/digiphot/apphot/fitsky/apfitsky.x
new file mode 100644
index 00000000..1c105b79
--- /dev/null
+++ b/noao/digiphot/apphot/fitsky/apfitsky.x
@@ -0,0 +1,452 @@
+include "../lib/apphotdef.h"
+include "../lib/apphot.h"
+include "../lib/noisedef.h"
+include "../lib/fitskydef.h"
+include "../lib/fitsky.h"
+
+# APFITSKY -- Procedure to the compute the sky value in an annular region
+# around a given position in the IRAF image.
+
+int procedure apfitsky (ap, im, wx, wy, sd, gd)
+
+pointer ap # pointer to the apphot structure
+pointer im # pointer to the IRAF image
+real wx # object x coordinate
+real wy # object y coordinate
+int sd # pointer to input text file containing sky values
+pointer gd # pointer to graphics stream
+
+int ier, nclip, nsky, ilo, ihi
+pointer sky, nse, gt
+real x, y
+int apskybuf(), ap_mode(), ap_centroid(), ap_histplot(), ap_readsky()
+int ap_median(), ap_radplot(), ap_gauss(), ap_lgsky(), ap_crosscor()
+int ap_mean(), ap_clip()
+pointer ap_gtinit()
+
+begin
+ # Initialize.
+ sky = AP_PSKY(ap)
+ nse = AP_NOISE(ap)
+ AP_SXCUR(sky) = wx
+ AP_SYCUR(sky) = wy
+ if (IS_INDEFR(wx) || IS_INDEFR(wy)) {
+ AP_OSXCUR(sky) = INDEFR
+ AP_OSYCUR(sky) = INDEFR
+ } else {
+ switch (AP_WCSOUT(ap)) {
+ case WCS_WORLD, WCS_PHYSICAL:
+ call ap_ltoo (ap, wx, wy, AP_OSXCUR(sky), AP_OSYCUR(sky), 1)
+ case WCS_TV:
+ call ap_ltov (im, wx, wy, AP_OSXCUR(sky), AP_OSYCUR(sky), 1)
+ default:
+ AP_OSXCUR(sky) = wx
+ AP_OSYCUR(sky) = wy
+ }
+ }
+
+ AP_SKY_MODE(sky) = INDEFR
+ AP_SKY_SIG(sky) = INDEFR
+ AP_SKY_SKEW(sky) = INDEFR
+ AP_NSKY(sky) = 0
+ AP_NSKY_REJECT(sky) = 0
+ if (IS_INDEFR(wx) || IS_INDEFR(wy))
+ return (AP_NOSKYAREA)
+
+ switch (AP_SKYFUNCTION(sky)) {
+
+ case AP_MEAN:
+
+ # Fetch the sky pixels.
+ ier = apskybuf (ap, im, wx, wy)
+ if (ier != AP_OK)
+ return (ier)
+
+ # Initialze the weights.
+ call amovkr (1.0, Memr[AP_SWGT(sky)], AP_NSKYPIX(sky))
+
+ # Clip the data.
+ if (AP_SLOCLIP(sky) > 0.0 || AP_SHICLIP(sky) > 0.0) {
+ nclip = ap_clip (Memr[AP_SKYPIX(sky)], Memi[AP_INDEX(sky)],
+ AP_NSKYPIX(sky), AP_SLOCLIP(sky), AP_SHICLIP(sky), ilo,
+ ihi)
+ if (nclip >= AP_NSKYPIX(sky))
+ return (AP_NSKY_TOO_SMALL)
+ nsky = AP_NSKYPIX(sky) - nclip
+ } else {
+ nclip = 0
+ call ap_index (Memi[AP_INDEX(sky)], AP_NSKYPIX(sky))
+ ilo = 1
+ nsky = AP_NSKYPIX(sky)
+ }
+
+ # Compute the mean of the sky pixel distribution with pixel
+ # rejection and region growing.
+ ier = ap_mean (Memr[AP_SKYPIX(sky)], Memi[AP_COORDS(sky)],
+ Memr[AP_SWGT(sky)], Memi[AP_INDEX(sky)+ilo-1],
+ nsky, AP_SNX(sky), AP_SNY(sky), AP_SLOREJECT(sky),
+ AP_SHIREJECT(sky), AP_RGROW(sky) * AP_SCALE(ap),
+ AP_SNREJECT(sky), AP_SKY_MODE(sky), AP_SKY_SIG(sky),
+ AP_SKY_SKEW(sky), AP_NSKY(sky), AP_NSKY_REJECT(sky))
+ AP_NSKY_REJECT(sky) = AP_NBADSKYPIX(sky) + nclip +
+ AP_NSKY_REJECT(sky)
+
+ return (ier)
+
+ case AP_MEDIAN:
+
+ # Fetch the sky pixels.
+ ier = apskybuf (ap, im, wx, wy)
+ if (ier != AP_OK)
+ return (ier)
+
+ # Initialze the weights.
+ call amovkr (1.0, Memr[AP_SWGT(sky)], AP_NSKYPIX(sky))
+
+ # Clip the data.
+ if (AP_SLOCLIP(sky) > 0.0 || AP_SHICLIP(sky) > 0.0) {
+ nclip = ap_clip (Memr[AP_SKYPIX(sky)], Memi[AP_INDEX(sky)],
+ AP_NSKYPIX(sky), AP_SLOCLIP(sky), AP_SHICLIP(sky), ilo,
+ ihi)
+ if (nclip >= AP_NSKYPIX(sky))
+ return (AP_NSKY_TOO_SMALL)
+ nsky = AP_NSKYPIX(sky) - nclip
+ } else {
+ nclip = 0
+ call apqsort (Memr[AP_SKYPIX(sky)], Memi[AP_INDEX(sky)],
+ Memi[AP_INDEX(sky)], AP_NSKYPIX(sky))
+ ilo = 1
+ nsky = AP_NSKYPIX(sky)
+ }
+
+ # Compute the median of the sky pixel distribution with pixel
+ # rejection and region growing.
+ ier = ap_median (Memr[AP_SKYPIX(sky)], Memi[AP_COORDS(sky)],
+ Memr[AP_SWGT(sky)], Memi[AP_INDEX(sky)+ilo-1], nsky,
+ AP_SNX(sky), AP_SNY(sky), AP_SLOREJECT(sky),
+ AP_SHIREJECT(sky), AP_RGROW(sky) * AP_SCALE(ap),
+ AP_SNREJECT(sky), AP_SKY_MODE(sky), AP_SKY_SIG(sky),
+ AP_SKY_SKEW(sky), AP_NSKY(sky), AP_NSKY_REJECT(sky))
+ AP_NSKY_REJECT(sky) = AP_NBADSKYPIX(sky) + nclip +
+ AP_NSKY_REJECT(sky)
+
+ return (ier)
+
+ case AP_MODE:
+
+ # Fetch the sky pixels.
+ ier = apskybuf (ap, im, wx, wy)
+ if (ier != AP_OK)
+ return (ier)
+
+ # Initialze the weights.
+ call amovkr (1.0, Memr[AP_SWGT(sky)], AP_NSKYPIX(sky))
+
+ # Clip the data.
+ if (AP_SLOCLIP(sky) > 0.0 || AP_SHICLIP(sky) > 0.0) {
+ nclip = ap_clip (Memr[AP_SKYPIX(sky)], Memi[AP_INDEX(sky)],
+ AP_NSKYPIX(sky), AP_SLOCLIP(sky), AP_SHICLIP(sky), ilo,
+ ihi)
+ if (nclip >= AP_NSKYPIX(sky))
+ return (AP_NSKY_TOO_SMALL)
+ nsky = AP_NSKYPIX(sky) - nclip
+ } else {
+ nclip = 0
+ call apqsort (Memr[AP_SKYPIX(sky)], Memi[AP_INDEX(sky)],
+ Memi[AP_INDEX(sky)], AP_NSKYPIX(sky))
+ ilo = 1
+ nsky = AP_NSKYPIX(sky)
+ }
+
+ # Compute the median of the sky pixel distribution with pixel
+ # rejection and region growing.
+ ier = ap_mode (Memr[AP_SKYPIX(sky)], Memi[AP_COORDS(sky)],
+ Memr[AP_SWGT(sky)], Memi[AP_INDEX(sky)+ilo-1], nsky,
+ AP_SNX(sky), AP_SNY(sky), AP_SLOREJECT(sky),
+ AP_SHIREJECT(sky), AP_RGROW(sky) * AP_SCALE(ap),
+ AP_SNREJECT(sky), AP_SKY_MODE(sky), AP_SKY_SIG(sky),
+ AP_SKY_SKEW(sky), AP_NSKY(sky), AP_NSKY_REJECT(sky))
+ AP_NSKY_REJECT(sky) = AP_NBADSKYPIX(sky) + nclip +
+ AP_NSKY_REJECT(sky)
+
+ return (ier)
+
+ case AP_CENTROID:
+
+ # Fetch the sky pixels.
+ ier = apskybuf (ap, im, wx, wy)
+ if (ier != AP_OK)
+ return (ier)
+
+ # Initialze the weights.
+ call amovkr (1.0, Memr[AP_SWGT(sky)], AP_NSKYPIX(sky))
+
+ # Clip the data.
+ if (AP_SLOCLIP(sky) > 0.0 || AP_SHICLIP(sky) > 0.0) {
+ nclip = ap_clip (Memr[AP_SKYPIX(sky)], Memi[AP_INDEX(sky)],
+ AP_NSKYPIX(sky), AP_SLOCLIP(sky), AP_SHICLIP(sky), ilo,
+ ihi)
+ if (nclip >= AP_NSKYPIX(sky))
+ return (AP_NSKY_TOO_SMALL)
+ nsky = AP_NSKYPIX(sky) - nclip
+ } else {
+ nclip = 0
+ call ap_index (Memi[AP_INDEX(sky)], AP_NSKYPIX(sky))
+ ilo = 1
+ nsky = AP_NSKYPIX(sky)
+ }
+
+ # Compute the sky value by performing a moment analysis of the
+ # sky pixel histogram.
+ ier = ap_centroid (Memr[AP_SKYPIX(sky)], Memi[AP_COORDS(sky)],
+ Memr[AP_SWGT(sky)], Memi[AP_INDEX(sky)+ilo-1], nsky,
+ AP_SNX(sky), AP_SNY(sky), AP_K1(sky), INDEFR,
+ AP_BINSIZE(sky), AP_SMOOTH(sky), AP_SLOREJECT(sky),
+ AP_SHIREJECT(sky), AP_RGROW(sky) * AP_SCALE(ap),
+ AP_SMAXITER(sky), AP_SKY_MODE(sky), AP_SKY_SIG(sky),
+ AP_SKY_SKEW(sky), AP_NSKY(sky), AP_NSKY_REJECT(sky))
+ AP_NSKY_REJECT(sky) = AP_NBADSKYPIX(sky) + nclip +
+ AP_NSKY_REJECT(sky)
+
+ return (ier)
+
+ case AP_CONSTANT:
+
+ # Set the sky value to a constant.
+ AP_SKY_MODE(sky) = AP_SKYBACKGROUND(sky)
+ AP_SKY_SIG(sky) = AP_SKYSIGMA(nse)
+ AP_SKY_SKEW(sky) = INDEFR
+ AP_NSKY(sky) = 0
+ AP_NSKY_REJECT(sky) = 0
+ return (AP_OK)
+
+ case AP_SKYFILE:
+
+ # Read the sky values from a file.
+ if (sd == NULL)
+ return (AP_NOSKYFILE)
+ ier = ap_readsky (sd, x, y, AP_SKY_MODE(sky), AP_SKY_SIG(sky),
+ AP_SKY_SKEW(sky), AP_NSKY(sky), AP_NSKY_REJECT(sky))
+ if (ier == EOF)
+ return (AP_EOFSKYFILE)
+ else if (ier != 7)
+ return (AP_BADSKYSCAN)
+ else if (AP_NSKY(sky) <= 0)
+ return (AP_NOSKYAREA)
+ else
+ return (AP_OK)
+
+ case AP_RADPLOT:
+
+ # Check the status of the graphics stream.
+ if (gd == NULL)
+ return (AP_NOGRAPHICS)
+
+ # Fetch the sky pixels.
+ ier = apskybuf (ap, im, wx, wy)
+ if (ier != AP_OK)
+ return (ier)
+
+ # Clip the data.
+ if (AP_SLOCLIP(sky) > 0.0 || AP_SHICLIP(sky) > 0.0) {
+ nclip = ap_clip (Memr[AP_SKYPIX(sky)], Memi[AP_INDEX(sky)],
+ AP_NSKYPIX(sky), AP_SLOCLIP(sky), AP_SHICLIP(sky), ilo,
+ ihi)
+ if (nclip >= AP_NSKYPIX(sky))
+ return (AP_NSKY_TOO_SMALL)
+ nsky = AP_NSKYPIX(sky) - nclip
+ } else {
+ nclip = 0
+ call ap_index (Memi[AP_INDEX(sky)], AP_NSKYPIX(sky))
+ ilo = 1
+ nsky = AP_NSKYPIX(sky)
+ }
+
+ # Mark the sky level on the radial profile plot.
+ call gactivate (gd, 0)
+ gt = ap_gtinit (AP_IMROOT(ap), wx, wy)
+ ier = ap_radplot (gd, gt, Memr[AP_SKYPIX(sky)],
+ Memi[AP_COORDS(sky)], Memi[AP_INDEX(sky)+ilo-1], nsky,
+ AP_SXC(sky), AP_SYC(sky), AP_SNX(sky), AP_SNY(sky),
+ AP_SCALE(ap), AP_SKY_MODE(sky), AP_SKY_SKEW(sky),
+ AP_SKY_SIG(sky), AP_NSKY(sky), AP_NSKY_REJECT(sky))
+ AP_NSKY_REJECT(sky) = AP_NBADSKYPIX(sky) + nclip +
+ AP_NSKY_REJECT(sky)
+ call ap_gtfree (gt)
+ call gdeactivate (gd, 0)
+
+ return (ier)
+
+ case AP_HISTPLOT:
+
+ # Check the status of the graphics stream.
+ if (gd == NULL)
+ return (AP_NOGRAPHICS)
+
+ # Fetch the sky pixels.
+ ier = apskybuf (ap, im, wx, wy)
+ if (ier != AP_OK)
+ return (ier)
+
+ # Initialze the weights.
+ call amovkr (1.0, Memr[AP_SWGT(sky)], AP_NSKYPIX(sky))
+
+ # Clip the data.
+ if (AP_SLOCLIP(sky) > 0.0 || AP_SHICLIP(sky) > 0.0) {
+ nclip = ap_clip (Memr[AP_SKYPIX(sky)], Memi[AP_INDEX(sky)],
+ AP_NSKYPIX(sky), AP_SLOCLIP(sky), AP_SHICLIP(sky),
+ ilo, ihi)
+ if (nclip >= AP_NSKYPIX(sky))
+ return (AP_NSKY_TOO_SMALL)
+ nsky = AP_NSKYPIX(sky) - nclip
+ } else {
+ nclip = 0
+ call ap_index (Memi[AP_INDEX(sky)], AP_NSKYPIX(sky))
+ ilo = 1
+ nsky = AP_NSKYPIX(sky)
+ }
+
+ # Mark the peak of the histogram on the histogram plot.
+ #call gactivate (gd, 0)
+ gt = ap_gtinit (AP_IMROOT(ap), wx, wy)
+ ier = ap_histplot (gd, gt, Memr[AP_SKYPIX(sky)],
+ Memr[AP_SWGT(sky)], Memi[AP_INDEX(sky)+ilo-1], nsky,
+ AP_K1(sky), INDEFR, AP_BINSIZE(sky), AP_SMOOTH(sky),
+ AP_SKY_MODE(sky), AP_SKY_SIG(sky), AP_SKY_SKEW(sky),
+ AP_NSKY(sky), AP_NSKY_REJECT(sky))
+ AP_NSKY_REJECT(sky) = AP_NBADSKYPIX(sky) + nclip +
+ AP_NSKY_REJECT(sky)
+ call ap_gtfree (gt)
+ #call gdeactivate (gd, 0)
+
+ return (ier)
+
+ case AP_OFILT:
+
+ # Fetch the sky pixels.
+ ier = apskybuf (ap, im, wx, wy)
+ if (ier != AP_OK)
+ return (ier)
+
+ # Initialze the weights.
+ call amovkr (1.0, Memr[AP_SWGT(sky)], AP_NSKYPIX(sky))
+
+ # Clip the data.
+ if (AP_SLOCLIP(sky) > 0.0 || AP_SHICLIP(sky) > 0.0) {
+ nclip = ap_clip (Memr[AP_SKYPIX(sky)], Memi[AP_INDEX(sky)],
+ AP_NSKYPIX(sky), AP_SLOCLIP(sky), AP_SHICLIP(sky),
+ ilo, ihi)
+ if (nclip >= AP_NSKYPIX(sky))
+ return (AP_NSKY_TOO_SMALL)
+ nsky = AP_NSKYPIX(sky) - nclip
+ } else {
+ nclip = 0
+ call ap_index (Memi[AP_INDEX(sky)], AP_NSKYPIX(sky))
+ ilo = 1
+ nsky = AP_NSKYPIX(sky)
+ }
+
+ # Compute the sky value using the histogram of the sky pixels
+ # and a variation of the optimal filtering technique.
+ ier = ap_lgsky (Memr[AP_SKYPIX(sky)], Memi[AP_COORDS(sky)],
+ Memr[AP_SWGT(sky)], Memi[AP_INDEX(sky)+ilo-1], nsky,
+ AP_SNX(sky), AP_SNY(sky), AP_K1(sky), INDEFR,
+ AP_BINSIZE(sky), AP_SMOOTH(sky), AP_SLOREJECT(sky),
+ AP_SHIREJECT(sky), AP_RGROW(sky) * AP_SCALE(ap),
+ AP_SMAXITER(sky), AP_SKY_MODE(sky), AP_SKY_SIG(sky),
+ AP_SKY_SKEW(sky), AP_NSKY(sky), AP_NSKY_REJECT(sky))
+ AP_NSKY_REJECT(sky) = AP_NBADSKYPIX(sky) + nclip +
+ AP_NSKY_REJECT(sky)
+
+ return (ier)
+
+ case AP_GAUSS:
+
+ # Fetch the sky pixels.
+ ier = apskybuf (ap, im, wx, wy)
+ if (ier != AP_OK)
+ return (ier)
+
+ # Initialze the weights.
+ call amovkr (1.0, Memr[AP_SWGT(sky)], AP_NSKYPIX(sky))
+
+ # Clip the data.
+ if (AP_SLOCLIP(sky) > 0.0 || AP_SHICLIP(sky) > 0.0) {
+ nclip = ap_clip (Memr[AP_SKYPIX(sky)], Memi[AP_INDEX(sky)],
+ AP_NSKYPIX(sky), AP_SLOCLIP(sky), AP_SHICLIP(sky),
+ ilo, ihi)
+ if (nclip >= AP_NSKYPIX(sky))
+ return (AP_NSKY_TOO_SMALL)
+ nsky = AP_NSKYPIX(sky) - nclip
+ } else {
+ nclip = 0
+ call ap_index (Memi[AP_INDEX(sky)], AP_NSKYPIX(sky))
+ ilo = 1
+ nsky = AP_NSKYPIX(sky)
+ }
+
+ # Compute the sky value by a fitting a skewed Gaussian function
+ # to the sky pixel histogram.
+ ier = ap_gauss (Memr[AP_SKYPIX(sky)], Memi[AP_COORDS(sky)],
+ Memr[AP_SWGT(sky)], Memi[AP_INDEX(sky)+ilo-1], nsky,
+ AP_SNX(sky), AP_SNY(sky), AP_SMAXITER(sky), AP_K1(sky),
+ INDEFR, AP_BINSIZE(sky), AP_SMOOTH(sky),
+ AP_SLOREJECT(sky), AP_SHIREJECT(sky), AP_RGROW(sky) *
+ AP_SCALE(ap), AP_SNREJECT(sky), AP_SKY_MODE(sky),
+ AP_SKY_SIG(sky), AP_SKY_SKEW(sky), AP_NSKY(sky),
+ AP_NSKY_REJECT(sky))
+ AP_NSKY_REJECT(sky) = AP_NBADSKYPIX(sky) + nclip +
+ AP_NSKY_REJECT(sky)
+
+ return (ier)
+
+ case AP_CROSSCOR:
+
+ # Fetch the sky pixels.
+ ier = apskybuf (ap, im, wx, wy)
+ if (ier != AP_OK)
+ return (ier)
+
+ # Initialze the weights.
+ call amovkr (1.0, Memr[AP_SWGT(sky)], AP_NSKYPIX(sky))
+
+ # Clip the data.
+ if (AP_SLOCLIP(sky) > 0.0 || AP_SHICLIP(sky) > 0.0) {
+ nclip = ap_clip (Memr[AP_SKYPIX(sky)], Memi[AP_INDEX(sky)],
+ AP_NSKYPIX(sky), AP_SLOCLIP(sky), AP_SHICLIP(sky),
+ ilo, ihi)
+ if (nclip >= AP_NSKYPIX(sky))
+ return (AP_NSKY_TOO_SMALL)
+ nsky = AP_NSKYPIX(sky) - nclip
+ } else {
+ nclip = 0
+ call ap_index (Memi[AP_INDEX(sky)], AP_NSKYPIX(sky))
+ ilo = 1
+ nsky = AP_NSKYPIX(sky)
+ }
+
+ # Fit the sky value by computing the cross-correlation
+ # function of the histogram and an estimate of the noise
+ # distribution.
+ ier = ap_crosscor (Memr[AP_SKYPIX(sky)], Memi[AP_COORDS(sky)],
+ Memr[AP_SWGT(sky)], Memi[AP_INDEX(sky)+ilo-1], nsky,
+ AP_SNX(sky), AP_SNY(sky), AP_K1(sky), INDEFR,
+ AP_BINSIZE(sky), AP_SMOOTH(sky), AP_SLOREJECT(sky),
+ AP_SHIREJECT(sky), AP_RGROW(sky) * AP_SCALE(ap),
+ AP_SMAXITER(sky), AP_SKY_MODE(sky), AP_SKY_SIG(sky),
+ AP_SKY_SKEW(sky), AP_NSKY(sky), AP_NSKY_REJECT(sky))
+ AP_NSKY_REJECT(sky) = AP_NBADSKYPIX(sky) + nclip +
+ AP_NSKY_REJECT(sky)
+
+ return (ier)
+
+ default:
+
+ AP_SKY_MODE(sky) = INDEFR
+ AP_SKY_SIG(sky) = INDEFR
+ AP_SKY_SKEW(sky) = INDEFR
+ AP_NSKY(sky) = AP_NSKYPIX(sky)
+ AP_NSKY_REJECT(sky) = 0
+ return (AP_OK)
+ }
+end
generated by cgit (git 2.34.1) at 2025-09-20 09:24:50 -0400