aboutsummaryrefslogtreecommitdiff
path: root/noao/imred/dtoi/spotlist.x
diff options
context:
space:
mode:
Diffstat (limited to 'noao/imred/dtoi/spotlist.x')
-rw-r--r--noao/imred/dtoi/spotlist.x395
1 files changed, 395 insertions, 0 deletions
diff --git a/noao/imred/dtoi/spotlist.x b/noao/imred/dtoi/spotlist.x
new file mode 100644
index 00000000..8f4ffeee
--- /dev/null
+++ b/noao/imred/dtoi/spotlist.x
@@ -0,0 +1,395 @@
+include <error.h>
+include <imhdr.h>
+include <mach.h>
+include <fset.h>
+
+# T_SPOTLIST -- calculates the densities and standard deviations of calibration
+# spots read in as IRAF images. Entries for density, sdev, spot number
+# and number of pixels used in the average are made in the output database.
+# These values are entered for all spots and the fog level. The
+# fog level is calculated but not subtracted in this procedure.
+
+procedure t_spotlist ()
+
+pointer db, sp, den, ngpix, sdev, dloge, option, npix
+int spot_fd, fog_fd, fngpix, nspots, maxad, fnpix
+real scale, sigma, fsdev, fog
+double maxden
+
+pointer ddb_map()
+int imtopenp(), imtlen(), clgeti(), strncmp()
+real clgetr()
+
+begin
+ call smark (sp)
+ call salloc (dloge, SZ_FNAME, TY_CHAR)
+ call salloc (option, SZ_FNAME, TY_CHAR)
+
+ # Get parameters and open file name templates
+ spot_fd = imtopenp ("spots")
+ fog_fd = imtopenp ("fogs")
+ call clgstr ("database", Memc[dloge], SZ_FNAME)
+ scale = clgetr ("scale")
+ maxad = clgeti ("maxad")
+ sigma = clgetr ("sigma")
+ call clgstr ("option", Memc[option], SZ_FNAME)
+
+ maxden = maxad * double (scale)
+
+ # Allocate space for density, standard deviation and ngpix arrays
+ nspots = imtlen (spot_fd)
+ call salloc ( den, nspots, TY_REAL)
+ call salloc (sdev, nspots, TY_REAL)
+ call salloc (ngpix, nspots, TY_INT)
+ call salloc ( npix, nspots, TY_INT)
+
+ # Calculate densities depending on algorithm option. The
+ # number of saturated pixels per spot is also calculated now.
+
+ if (strncmp (Memc[option], "median", 3) == 0)
+ call hd_median (spot_fd, Memr[den], Memr[sdev], Memi[ngpix],
+ nspots, scale, Memi[npix])
+ else
+ call hd_mean (spot_fd, Memr[den], Memr[sdev], Memi[ngpix],
+ nspots, scale, sigma, Memi[npix])
+
+ # Calculate fog level and count saturated pixels
+ call hd_fogcalc (fog_fd, fog, fsdev, fngpix, scale, sigma,
+ Memc[option], fnpix)
+
+ # Now print results to stdout
+ call hd_printit (Memr[den], Memr[sdev], Memi[npix], Memi[ngpix],
+ fog, fsdev, fnpix, fngpix, nspots)
+
+ # Open output database file and write spot information
+ db = ddb_map (Memc[dloge], APPEND)
+ call hd_wspotdb (db, Memr[den], Memr[sdev], Memi[ngpix], nspots)
+
+ # Write fog information to database as single record
+ call ddb_prec (db, "fog")
+ call ddb_putr (db, "density", fog)
+ call ddb_putr (db, "sdev", fsdev)
+ call ddb_puti (db, "ngpix", fngpix)
+
+ # Scale info gets written to database also (very precisely!)
+ call ddb_prec (db, "common")
+ call ddb_putr (db, "scale", scale)
+ call ddb_putd (db, "maxden", maxden)
+ call ddb_pstr (db, "option", Memc[option])
+
+ call ddb_unmap (db)
+
+ call imtclose (spot_fd)
+ call imtclose (fog_fd)
+ call sfree (sp)
+end
+
+
+# HD_MEAN -- Calculate mean density of calibration spots.
+
+procedure hd_mean (spot_fd, den, sdev, ngpix, nspots, scale, sigma, npix)
+
+int spot_fd # File descriptor for list of spots
+real den[ARB] # Mean density values - filled on return
+real sdev[ARB] # Standard deviation array - filled on return
+int ngpix[ARB] # Number of unrejected pixels - filled on return
+int nspots # Number of spots in list
+real scale # Scale for voltage to density conversion
+real sigma # Rejection criteria set by user
+int npix[ARB] # Number of pixels per spot
+
+pointer im, spot, sp, pix
+int i, junk, ncols, nlines
+pointer immap(), imgs2r()
+int imtgetim(), hd_aravr()
+errchk imgs2r, amulkr, hd_aravr
+
+begin
+ call smark (sp)
+ call salloc (spot, SZ_FNAME, TY_CHAR)
+
+ # Loop over spot rasters. Calculate density and standard deviation.
+ for (i = 1; i <= nspots; i = i + 1) {
+ junk = imtgetim (spot_fd, Memc[spot], SZ_FNAME)
+ iferr (im = immap (Memc[spot], READ_ONLY, 0)) {
+ call erract (EA_WARN)
+ next
+ }
+
+ ncols = IM_LEN(im,1)
+ nlines = IM_LEN(im,2)
+ npix[i] = ncols * nlines
+
+ # For all pixels in the image, scale the a/d value to density and
+ # calculate the mean value, using a mean rejection algorithm.
+
+ pix = imgs2r (im, 1, ncols, 1, nlines)
+ call amulkr (Memr[pix], scale, Memr[pix], npix[i])
+ ngpix[i] = hd_aravr (Memr[pix], npix[i], den[i], sdev[i], sigma)
+ call imunmap (im)
+ }
+
+ call sfree (sp)
+end
+
+
+# HD_MEDIAN -- Calculate median density of calibration spots.
+
+procedure hd_median (spot_fd, den, sdev, ngpix, nspots, scale, npix)
+
+int spot_fd # File descriptor for list of spots
+real den[ARB] # Mean density values - filled on return
+real sdev[ARB] # Standard deviation of input spots
+int ngpix[ARB] # Number of pixels not rejected
+int nspots # Number of spots in list
+real scale # Scale for voltage to density conversion
+int npix[ARB] # Number of pixels per spot
+
+pointer im, spot, sp, pix
+int i, junk, ncols, nlines
+real junk_mean
+pointer immap(), imgs2r()
+int imtgetim()
+real amedr()
+errchk imgs2r, amulkr, amedr
+
+begin
+ call smark (sp)
+ call salloc (spot, SZ_FNAME, TY_CHAR)
+
+ # Loop over spot rasters. Calculate density and standard deviation.
+ for (i = 1; i <= nspots; i = i + 1) {
+ junk = imtgetim (spot_fd, Memc[spot], SZ_FNAME)
+ iferr (im = immap (Memc[spot], READ_ONLY, 0)) {
+ call erract (EA_WARN)
+ next
+ }
+
+ ncols = IM_LEN(im,1)
+ nlines = IM_LEN(im,2)
+ npix[i] = ncols * nlines
+
+ # For all pixels in the image, scale the a/d value to density and
+ # calculate the median value. For the user's information, the
+ # sigma is also calculated in a call to aavgr.
+
+ pix = imgs2r (im, 1, ncols, 1, nlines)
+ call amulkr (Memr[pix], scale, Memr[pix], npix[i])
+ den[i] = amedr (Memr[pix], npix[i])
+ ngpix[i] = npix[i]
+ call aavgr (Memr[pix], npix[i], junk_mean, sdev[i])
+ call imunmap (im)
+ }
+
+ call sfree (sp)
+end
+
+
+# HD_FOGCALC -- Calculate fog density.
+
+procedure hd_fogcalc (fog_fd, fog, fsdev, fngpix, scale, sigma, option, nfpix)
+
+int fog_fd # File descriptor for list of fog spots
+real fog # Mean fog value - returned
+real fsdev # Standard deviation - returned
+int fngpix # Number of pixels used - returned
+real scale # Voltage to density scaling factor
+real sigma # Rejection criteria - set by user
+char option[ARB] # User's choice of mean/median algorithm
+int nfpix # Total number of fog pixels
+
+pointer pix, im, sp, fogfile, ptr
+int nfog, maxfog, i, junk, ncols, nlines, npix, total_pix, op
+real junk_mean
+pointer immap(), imgs2r()
+int imtlen(), imtgetim(), hd_aravr(), strncmp()
+real amedr()
+errchk calloc, imgs2r, aaddr, amulkr, hd_aravr
+
+begin
+ call smark (sp)
+ call salloc (fogfile, SZ_FNAME, TY_CHAR)
+
+ pix = NULL
+ total_pix = 0
+ op = 0
+ nfog = imtlen (fog_fd)
+ maxfog = nfog
+ do i = 1, maxfog {
+ junk = imtgetim (fog_fd, Memc[fogfile], SZ_FNAME)
+ iferr (im = immap (Memc[fogfile], READ_ONLY, 0)) {
+ call erract (EA_WARN)
+ nfog = nfog - 1
+ next
+ }
+
+ ncols = IM_LEN(im,1)
+ nlines = IM_LEN(im,2)
+ npix = ncols * nlines
+ total_pix = total_pix + npix
+
+ if (pix == NULL)
+ # Initialize space for accumulating pixels
+ call calloc (pix, npix, TY_REAL)
+ else
+ # Increase space for accumulating pixels
+ call realloc (pix, total_pix, TY_REAL)
+
+ # Build up vector of fog pixels
+ ptr = imgs2r (im, 1, ncols, 1, nlines)
+ call amovr (Memr[ptr], Memr[pix+op], npix)
+ op = op + npix
+
+ call imunmap (im)
+ }
+
+ # Scale values to density and calculate fog and std deviation
+
+ if (nfog > 0) {
+# 7 Sept 1989, S. Rooke: in Suzanne's absence, made following bugfix after
+# bug reported by Steve Majewski that fog values are off by 1/n images where
+# multiple fog images are used in a single run. The total_pix already contains
+# the sum of all pixel values, so the fog pixel values should not be divided
+# by nfog. This should be verified by Suzanne on her return, and these comments
+# removed.
+# call amulkr (Memr[pix], scale / real (nfog), Memr[pix], total_pix)
+ call amulkr (Memr[pix], scale, Memr[pix], total_pix)
+ if (strncmp (option, "median", 3) == 0) {
+ fog = amedr (Memr[pix], total_pix)
+ fngpix = total_pix
+ call aavgr (Memr[pix], total_pix, junk_mean, fsdev)
+ } else
+ fngpix = hd_aravr (Memr[pix], total_pix, fog, fsdev, sigma)
+ } else {
+ fog = 0.0
+ fsdev = 0.0
+ fngpix = 0
+ }
+ nfpix = total_pix
+
+ call mfree (pix, TY_REAL)
+ call sfree (sp)
+end
+
+
+# HD_WSPOTDB -- Write spot information to database file. Values are first
+# sorted in order of increasing density.
+
+procedure hd_wspotdb (db, density, sdev, ngpix, nspots)
+
+pointer db # Pointer to database
+real density[ARB] # Array of densities
+real sdev [ARB] # Array of standard deviations
+int ngpix [ARB] # Array of npix used in calculations
+int nspots # Number of density spots
+
+begin
+ if (density[1] > density[nspots]) {
+ # Need to reorder arrays
+ call hd_reorderr (density, nspots)
+ call hd_reorderr ( sdev, nspots)
+ call hd_reorderi ( ngpix, nspots)
+ }
+
+ call ddb_ptime (db)
+
+ # Density record
+ call ddb_prec (db, "density")
+ call ddb_par (db, "den_val", density, nspots)
+
+ # Standard deviation of density is written to a record
+ call ddb_prec (db, "standard deviation")
+ call ddb_par (db, "sdev_val", sdev, nspots)
+
+ # Record for npix_used
+ call ddb_prec (db, "ngpix")
+ call ddb_pai (db, "npix_val", ngpix, nspots)
+end
+
+
+# HD_REORDERR - Flip order of real array in place.
+
+procedure hd_reorderr (array, nvals)
+
+real array[ARB] # Real array to be reordered
+int nvals # Number of elements in array
+
+pointer sp, tmp
+int i
+
+begin
+ call smark (sp)
+ call salloc (tmp, nvals, TY_REAL)
+
+ call amovr (array, Memr[tmp], nvals)
+ do i = 1, nvals
+ array[i] = Memr[tmp+nvals-i]
+
+ call sfree (sp)
+end
+
+
+# HD_REORDERI -- Flip order of integer array in place.
+
+procedure hd_reorderi (array, nvals)
+
+int array[ARB] # Integer array to be ordered
+int nvals # Number of elements in array
+
+pointer sp, tmp
+int i
+
+begin
+ call smark (sp)
+ call salloc (tmp, nvals, TY_INT)
+
+ call amovi (array, Memi[tmp], nvals)
+ do i = 1, nvals
+ array[i] = Memi[tmp+nvals-i]
+
+ call sfree (sp)
+end
+
+
+# HD_PRINTIT -- Neatly printout all the accumulated information.
+
+procedure hd_printit (den, sdev, npix, ngpix, fog, fsdev, fnpix, fngpix, nspots)
+
+real den[ARB] # density array
+real sdev[ARB] # std deviation array
+int npix[ARB] # npix array
+int ngpix[ARB] # ngoodpix array
+real fog # fog value
+real fsdev # std deviation of fog
+int fnpix # npix in fog
+int fngpix # ngoodpix in fog
+int nspots # number of spots
+int i
+
+begin
+
+ call fseti (STDOUT, F_FLUSHNL, YES)
+
+ call printf ("\n # Number of P")
+ call printf ("ixels\n")
+ call printf ("# Spot Number Density Std Deviation Total Used Rej")
+ call printf ("ected \n")
+
+ do i = 1, nspots {
+ call printf (" %d %17t%.4f %27t%.4f %43t%d %5d %6d \n")
+ call pargi (i)
+ call pargr (den[i])
+ call pargr (sdev[i])
+ call pargi (npix[i])
+ call pargi (ngpix[i])
+ call pargi (npix[i] - ngpix[i])
+ }
+
+ call printf (" FOG %17t%.4f %27t%.4f %43t%d %5d %6d \n")
+ call pargr (fog)
+ call pargr (fsdev)
+ call pargi (fnpix)
+ call pargi (fngpix)
+ call pargi (fnpix - fngpix)
+
+end
generated by cgit (git 2.34.1) at 2025-09-20 09:04:11 -0400