aboutsummaryrefslogtreecommitdiff
path: root/pkg/images/imcoords/src/sffind.x
diff options
context:
space:
mode:
Diffstat (limited to 'pkg/images/imcoords/src/sffind.x')
-rw-r--r--pkg/images/imcoords/src/sffind.x739
1 files changed, 739 insertions, 0 deletions
diff --git a/pkg/images/imcoords/src/sffind.x b/pkg/images/imcoords/src/sffind.x
new file mode 100644
index 00000000..367893e5
--- /dev/null
+++ b/pkg/images/imcoords/src/sffind.x
@@ -0,0 +1,739 @@
+include <error.h>
+include <mach.h>
+include <imhdr.h>
+include <imset.h>
+include <fset.h>
+include <math.h>
+include "starfind.h"
+
+
+# SF_FIND -- Find stars in an image using a pattern matching technique and
+# a circularly symmetric Gaussian pattern.
+
+procedure sf_find (im, out, sf, nxblock, nyblock, wcs, wxformat, wyformat,
+ boundary, constant, verbose)
+
+pointer im #I pointer to the input image
+int out #I the output file descriptor
+pointer sf #I pointer to the apphot structure
+int nxblock #I the x dimension blocking factor
+int nyblock #I the y dimension blocking factor
+char wcs[ARB] #I the world coordinate system
+char wxformat[ARB] #I the x axis world coordinate format
+char wyformat[ARB] #I the y axis world coordinate format
+int boundary #I type of boundary extension
+real constant #I constant for constant boundary extension
+int verbose #I verbose switch
+
+int i, j, fwidth, swidth, norm
+int l1, l2, c1, c2, ncols, nlines, nxb, nyb, nstars, stid
+pointer sp, gker2d, ngker2d, skip, fmtstr, twxformat, twyformat
+pointer imbuf, denbuf, str, mw, ct
+real sigma, nsigma, a, b, c, f, gsums[LEN_GAUSS], relerr, dmin, dmax
+real maglo, maghi
+
+bool streq()
+int sf_stfind()
+pointer mw_openim(), mw_sctran()
+real sf_egkernel()
+errchk mw_openim(), mw_sctran(), mw_gattrs()
+
+begin
+ # Allocate working space.
+ call smark (sp)
+ call salloc (twxformat, SZ_FNAME, TY_CHAR)
+ call salloc (twyformat, SZ_FNAME, TY_CHAR)
+ call salloc (fmtstr, SZ_LINE, TY_CHAR)
+ call salloc (str, SZ_LINE, TY_CHAR)
+
+ # Compute the parameters of the Gaussian kernel.
+ sigma = HWHM_TO_SIGMA * SF_HWHMPSF(sf)
+ nsigma = SF_FRADIUS(sf) / HWHM_TO_SIGMA
+ call sf_egparams (sigma, 1.0, 0.0, nsigma, a, b, c, f, fwidth, fwidth)
+
+ # Compute the separation parameter
+ swidth = max (2, int (SF_SEPMIN(sf) * SF_HWHMPSF(sf) + 0.5))
+
+ # Compute the minimum and maximum pixel values.
+ if (IS_INDEFR(SF_DATAMIN(sf)) && IS_INDEFR(SF_DATAMAX(sf))) {
+ norm = YES
+ dmin = -MAX_REAL
+ dmax = MAX_REAL
+ } else {
+ norm = NO
+ if (IS_INDEFR(SF_DATAMIN(sf)))
+ dmin = -MAX_REAL
+ else
+ dmin = SF_DATAMIN(sf)
+ if (IS_INDEFR(SF_DATAMAX(sf)))
+ dmax = MAX_REAL
+ else
+ dmax = SF_DATAMAX(sf)
+ }
+
+ # Compute the magnitude limits
+ if (IS_INDEFR(SF_MAGLO(sf)))
+ maglo = -MAX_REAL
+ else
+ maglo = SF_MAGLO(sf)
+ if (IS_INDEFR(SF_MAGHI(sf)))
+ maghi = MAX_REAL
+ else
+ maghi = SF_MAGHI(sf)
+
+ # Open the image WCS.
+ if (wcs[1] == EOS) {
+ mw = NULL
+ ct = NULL
+ } else {
+ iferr {
+ mw = mw_openim (im)
+ } then {
+ call erract (EA_WARN)
+ mw = NULL
+ ct = NULL
+ } else {
+ iferr {
+ ct = mw_sctran (mw, "logical", wcs, 03B)
+ } then {
+ call erract (EA_WARN)
+ ct = NULL
+ call mw_close (mw)
+ mw = NULL
+ }
+ }
+ }
+
+ # Set the WCS formats.
+ if (ct == NULL)
+ call strcpy (wxformat, Memc[twxformat], SZ_FNAME)
+ else if (wxformat[1] == EOS) {
+ if (mw != NULL) {
+ iferr (call mw_gwattrs (mw, 1, "format", Memc[twxformat],
+ SZ_FNAME)) {
+ if (streq (wcs, "world"))
+ call strcpy ("%11.8g", Memc[twxformat], SZ_FNAME)
+ else
+ call strcpy ("%9.3f", Memc[twxformat], SZ_FNAME)
+ }
+ } else
+ call strcpy ("%9.3f", Memc[twxformat], SZ_FNAME)
+ } else
+ call strcpy (wxformat, Memc[twxformat], SZ_FNAME)
+ if (ct == NULL)
+ call strcpy (wyformat, Memc[twyformat], SZ_FNAME)
+ else if (wyformat[1] == EOS) {
+ if (mw != NULL) {
+ iferr (call mw_gwattrs (mw, 2, "format", Memc[twyformat],
+ SZ_FNAME)) {
+ if (streq (wcs, "world"))
+ call strcpy ("%11.8g", Memc[twyformat], SZ_FNAME)
+ else
+ call strcpy ("%9.3f", Memc[twyformat], SZ_FNAME)
+ }
+ } else
+ call strcpy ("%9.3f", Memc[twyformat], SZ_FNAME)
+ } else
+ call strcpy (wyformat, Memc[twyformat], SZ_FNAME)
+
+ # Create the output format string.
+ call sprintf (Memc[fmtstr],
+ SZ_LINE, " %s %s %s %s %s %s %s %s %s %s %s\n")
+ call pargstr ("%9.3f")
+ call pargstr ("%9.3f")
+ call pargstr (Memc[twxformat])
+ call pargstr (Memc[twyformat])
+ call pargstr ("%7.2f")
+ call pargstr ("%6d")
+ call pargstr ("%6.2f")
+ call pargstr ("%6.3f")
+ call pargstr ("%6.1f")
+ call pargstr ("%7.3f")
+ call pargstr ("%6d")
+
+ # Set up the image boundary extension characteristics.
+ call imseti (im, IM_TYBNDRY, boundary)
+ call imseti (im, IM_NBNDRYPIX, 1 + fwidth / 2 + swidth)
+ if (boundary == BT_CONSTANT)
+ call imsetr (im, IM_BNDRYPIXVAL, constant)
+
+ # Set up the blocking factor.
+ # Compute the magnitude limits
+ if (IS_INDEFI(nxblock))
+ nxb = IM_LEN(im,1)
+ else
+ nxb = nxblock
+ if (IS_INDEFI(nyblock))
+ nyb = IM_LEN(im,2)
+ else
+ nyb = nyblock
+
+ # Print the detection criteria on the standard output.
+ if (verbose == YES) {
+ call fstats (out, F_FILENAME, Memc[str], SZ_LINE)
+ call printf ("\nImage: %s Output: %s\n")
+ call pargstr (IM_HDRFILE(im))
+ call pargstr (Memc[str])
+ call printf ("Detection Parameters\n")
+ call printf (
+ " Hwhmpsf: %0.3f (pixels) Threshold: %g (ADU) Npixmin: %d\n")
+ call pargr (SF_HWHMPSF(sf))
+ call pargr (SF_THRESHOLD(sf))
+ call pargi (SF_NPIXMIN(sf))
+ call printf (" Datamin: %g (ADU) Datamax: %g (ADU)\n")
+ call pargr (SF_DATAMIN(sf))
+ call pargr (SF_DATAMAX(sf))
+ call printf (" Fradius: %0.3f (HWHM) Sepmin: %0.3f (HWHM)\n\n")
+ call pargr (SF_FRADIUS(sf))
+ call pargr (SF_SEPMIN(sf))
+ }
+
+ if (out != NULL) {
+ call fstats (out, F_FILENAME, Memc[str], SZ_LINE)
+ call fprintf (out, "\n# Image: %s Output: %s\n")
+ call pargstr (IM_HDRFILE(im))
+ call pargstr (Memc[str])
+ call fprintf (out, "# Detection Parameters\n")
+ call fprintf (out,
+ "# Hwhmpsf: %0.3f (pixels) Threshold: %g (ADU) Npixmin: %d\n")
+ call pargr (SF_HWHMPSF(sf))
+ call pargr (SF_THRESHOLD(sf))
+ call pargi (SF_NPIXMIN(sf))
+ call fprintf (out, "# Datamin: %g (ADU) Datamax: %g (ADU)\n")
+ call pargr (SF_DATAMIN(sf))
+ call pargr (SF_DATAMAX(sf))
+ call fprintf (out, "# Fradius: %g (HWHM) Sepmin: %g (HWHM)\n")
+ call pargr (SF_FRADIUS(sf))
+ call pargr (SF_SEPMIN(sf))
+ call fprintf (out, "# Selection Parameters\n")
+ call pargi (SF_NPIXMIN(sf))
+ call fprintf (out, "# Maglo: %0.3f Maghi: %0.3f\n")
+ call pargr (SF_MAGLO(sf))
+ call pargr (SF_MAGHI(sf))
+ call fprintf (out, "# Roundlo: %0.3f Roundhi: %0.3f\n")
+ call pargr (SF_ROUNDLO(sf))
+ call pargr (SF_ROUNDHI(sf))
+ call fprintf (out, "# Sharplo: %0.3f Sharphi: %0.3f\n")
+ call pargr (SF_SHARPLO(sf))
+ call pargr (SF_SHARPHI(sf))
+ call fprintf (out, "# Columns\n")
+ call fprintf (out, "# 1: X 2: Y \n")
+ if (ct == NULL) {
+ call fprintf (out, "# 3: Mag 4: Area\n")
+ call fprintf (out, "# 5: Hwhm 6: Roundness\n")
+ call fprintf (out, "# 7: Pa 8: Sharpness\n\n")
+ } else {
+ call fprintf (out, "# 3: Wx 4: Wy \n")
+ call fprintf (out, "# 5: Mag 6: Area\n")
+ call fprintf (out, "# 7: Hwhm 8: Roundness\n")
+ call fprintf (out, "# 9: Pa 10: Sharpness\n\n")
+ }
+ }
+
+ # Process the image block by block.
+ stid = 1
+ nstars = 0
+ do j = 1, IM_LEN(im,2), nyb {
+
+ l1 = j
+ l2 = min (IM_LEN(im,2), j + nyb - 1)
+ nlines = l2 - l1 + 1 + 2 * (fwidth / 2 + swidth)
+
+ do i = 1, IM_LEN(im,1), nxb {
+
+ # Allocate space for the convolution kernel.
+ call malloc (gker2d, fwidth * fwidth, TY_REAL)
+ call malloc (ngker2d, fwidth * fwidth, TY_REAL)
+ call malloc (skip, fwidth * fwidth, TY_INT)
+
+ # Allocate space for the data and the convolution.
+ c1 = i
+ c2 = min (IM_LEN(im,1), i + nxb - 1)
+ ncols = c2 - c1 + 1 + 2 * (fwidth / 2 + swidth)
+ call malloc (imbuf, ncols * nlines, TY_REAL)
+ call malloc (denbuf, ncols * nlines, TY_REAL)
+
+ # Compute the convolution kernels.
+ relerr = sf_egkernel (Memr[gker2d], Memr[ngker2d], Memi[skip],
+ fwidth, fwidth, gsums, a, b, c, f)
+
+ # Do the convolution.
+ if (norm == YES)
+ call sf_fconvolve (im, c1, c2, l1, l2, swidth, Memr[imbuf],
+ Memr[denbuf], ncols, nlines, Memr[ngker2d], Memi[skip],
+ fwidth, fwidth)
+ else
+ call sf_gconvolve (im, c1, c2, l1, l2, swidth, Memr[imbuf],
+ Memr[denbuf], ncols, nlines, Memr[gker2d], Memi[skip],
+ fwidth, fwidth, gsums, dmin, dmax)
+
+ # Find the stars.
+ nstars = sf_stfind (out, Memr[imbuf], Memr[denbuf], ncols,
+ nlines, c1, c2, l1, l2, swidth, Memi[skip], fwidth,
+ fwidth, SF_HWHMPSF(sf), SF_THRESHOLD(sf), dmin, dmax,
+ ct, SF_NPIXMIN(sf), maglo, maghi, SF_ROUNDLO(sf),
+ SF_ROUNDHI(sf), SF_SHARPLO(sf), SF_SHARPHI(sf),
+ Memc[fmtstr], stid, verbose)
+
+ # Increment the sequence number.
+ stid = stid + nstars
+
+ # Free the memory.
+ call mfree (imbuf, TY_REAL)
+ call mfree (denbuf, TY_REAL)
+ call mfree (gker2d, TY_REAL)
+ call mfree (ngker2d, TY_REAL)
+ call mfree (skip, TY_INT)
+ }
+ }
+
+ # Print out the selection parameters.
+ if (verbose == YES) {
+ call printf ("\nSelection Parameters\n")
+ call printf ( " Maglo: %0.3f Maghi: %0.3f\n")
+ call pargr (SF_MAGLO(sf))
+ call pargr (SF_MAGHI(sf))
+ call printf ( " Roundlo: %0.3f Roundhi: %0.3f\n")
+ call pargr (SF_ROUNDLO(sf))
+ call pargr (SF_ROUNDHI(sf))
+ call printf ( " Sharplo: %0.3f Sharphi: %0.3f\n")
+ call pargr (SF_SHARPLO(sf))
+ call pargr (SF_SHARPHI(sf))
+ }
+
+ if (mw != NULL) {
+ call mw_ctfree (ct)
+ call mw_close (mw)
+ }
+ call sfree (sp)
+end
+
+
+# SF_STFIND -- Detect images in the convolved image and then compute image
+# characteristics using the original image.
+
+int procedure sf_stfind (out, imbuf, denbuf, ncols, nlines, c1, c2, l1, l2,
+ sepmin, skip, nxk, nyk, hwhmpsf, threshold, datamin, datamax,
+ ct, nmin, maglo, maghi, roundlo, roundhi, sharplo, sharphi,
+ fmtstr, stid, verbose)
+
+int out #I the output file descriptor
+real imbuf[ncols,nlines] #I the input data buffer
+real denbuf[ncols,nlines] #I the input density enhancements buffer
+int ncols, nlines #I the dimensions of the input buffers
+int c1, c2 #I the image columns limits
+int l1, l2 #I the image lines limits
+int sepmin #I the minimum object separation
+int skip[nxk,ARB] #I the pixel fitting array
+int nxk, nyk #I the dimensions of the fitting array
+real hwhmpsf #I the HWHM of the PSF in pixels
+real threshold #I the threshold for object detection
+real datamin, datamax #I the minimum and maximum good data values
+pointer ct #I the coordinate transformation pointer
+int nmin #I the minimum number of good object pixels
+real maglo,maghi #I the magnitude estimate limits
+real roundlo,roundhi #I the ellipticity estimate limits
+real sharplo, sharphi #I the sharpness estimate limits
+char fmtstr[ARB] #I the format string
+int stid #U the object sequence number
+int verbose #I verbose mode
+
+int line1, line2, inline, xmiddle, ymiddle, ntotal, nobjs, nstars
+pointer sp, cols, sharp, x, y, ellip, theta, npix, mag, size
+int sf_detect(), sf_test()
+
+begin
+ # Set up useful line and column limits.
+ line1 = 1 + sepmin + nyk / 2
+ line2 = nlines - sepmin - nyk / 2
+ xmiddle = 1 + nxk / 2
+ ymiddle = 1 + nyk / 2
+
+ # Set up a cylindrical buffers and some working space for
+ # the detected images.
+ call smark (sp)
+ call salloc (cols, ncols, TY_INT)
+ call salloc (x, ncols, TY_REAL)
+ call salloc (y, ncols, TY_REAL)
+ call salloc (mag, ncols, TY_REAL)
+ call salloc (npix, ncols, TY_INT)
+ call salloc (size, ncols, TY_REAL)
+ call salloc (ellip, ncols, TY_REAL)
+ call salloc (theta, ncols, TY_REAL)
+ call salloc (sharp, ncols, TY_REAL)
+
+ # Generate the starlist line by line.
+ ntotal = 0
+ do inline = line1, line2 {
+
+ # Detect local maximum in the density enhancement buffer.
+ nobjs = sf_detect (denbuf[1,inline-nyk/2-sepmin], ncols, sepmin,
+ nxk, nyk, threshold, Memi[cols])
+ if (nobjs <= 0)
+ next
+
+ # Do not skip the middle pixel in the moments computation.
+ call sf_moments (imbuf[1,inline-nyk/2], denbuf[1,inline-nyk/2],
+ ncols, skip, nxk, nyk, Memi[cols], Memr[x], Memr[y],
+ Memi[npix], Memr[mag], Memr[size], Memr[ellip], Memr[theta],
+ Memr[sharp], nobjs, datamin, datamax, threshold, hwhmpsf,
+ real (-sepmin - nxk / 2 + c1 - 1), real (inline - sepmin -
+ nyk + l1 - 1))
+
+ # Test the image characeteristics of detected objects.
+ nstars = sf_test (Memi[cols], Memr[x], Memr[y], Memi[npix],
+ Memr[mag], Memr[size], Memr[ellip], Memr[theta], Memr[sharp],
+ nobjs, real (c1 - 0.5), real (c2 + 0.5), real (l1 - 0.5),
+ real (l2 + 0.5), nmin, maglo, maghi, roundlo, roundhi,
+ sharplo, sharphi)
+
+ # Print results on the standard output.
+ if (verbose == YES)
+ call sf_write (STDOUT, Memi[cols], Memr[x], Memr[y],
+ Memr[mag], Memi[npix], Memr[size], Memr[ellip],
+ Memr[theta], Memr[sharp], nstars, ct, fmtstr,
+ ntotal + stid)
+
+ # Save the results in the file.
+ call sf_write (out, Memi[cols], Memr[x], Memr[y], Memr[mag],
+ Memi[npix], Memr[size], Memr[ellip], Memr[theta],
+ Memr[sharp], nstars, ct, fmtstr, ntotal + stid)
+
+ ntotal = ntotal + nstars
+
+ }
+
+ # Free space
+ call sfree (sp)
+
+ return (ntotal)
+end
+
+
+# SF_DETECT -- Detect stellar objects in an image line. In order to be
+# detected as a star the candidate object must be above threshold and have
+# a maximum pixel value greater than any pixels within sepmin pixels.
+
+int procedure sf_detect (density, ncols, sepmin, nxk, nyk, threshold, cols)
+
+real density[ncols, ARB] #I the input density enhancements array
+int ncols #I the x dimension of the input array
+int sepmin #I the minimum separation in pixels
+int nxk, nyk #I size of the fitting area
+real threshold #I density threshold
+int cols[ARB] #O column numbers of detected stars
+
+int i, j, k, ymiddle, nxhalf, nyhalf, ny, b2, nobjs, rj2, r2
+define nextpix_ 11
+
+begin
+ ymiddle = 1 + nyk / 2 + sepmin
+ nxhalf = nxk / 2
+ nyhalf = nyk / 2
+ ny = 2 * sepmin + 1
+ b2 = sepmin ** 2
+
+ # Loop over all the columns in an image line.
+ nobjs = 0
+ for (i = 1 + nxhalf + sepmin; i <= ncols - nxhalf - sepmin; ) {
+
+ # Test whether the density enhancement is above threshold.
+ if (density[i,ymiddle] < threshold)
+ goto nextpix_
+
+ # Test whether a given density enhancement satisfies the
+ # separation criterion.
+ do j = 1, ny {
+ rj2 = (j - sepmin - 1) ** 2
+ do k = i - sepmin, i + sepmin {
+ r2 = (i - k) ** 2 + rj2
+ if (r2 <= b2) {
+ if (density[i,ymiddle] < density[k,j+nyhalf])
+ goto nextpix_
+ }
+ }
+ }
+
+ # Add the detected object to the list.
+ nobjs = nobjs + 1
+ cols[nobjs] = i
+
+ # If a local maximum is detected there can be no need to
+ # check pixels in this row between i and i + sepmin.
+ i = i + sepmin
+nextpix_
+ # Work on the next pixel.
+ i = i + 1
+ }
+
+ return (nobjs)
+end
+
+
+# SF_MOMENTS -- Perform a moments analysis on the dectected objects.
+
+procedure sf_moments (data, den, ncols, skip, nxk, nyk, cols, x, y,
+ npix, mag, size, ellip, theta, sharp, nobjs, datamin, datamax,
+ threshold, hwhmpsf, xoff, yoff)
+
+real data[ncols,ARB] #I the input data array
+real den[ncols,ARB] #I the input density enhancements array
+int ncols #I the x dimension of the input buffer
+int skip[nxk,ARB] #I the input fitting array
+int nxk, nyk #I the dimensions of the fitting array
+int cols[ARB] #I the input initial positions
+real x[ARB] #O the output x coordinates
+real y[ARB] #O the output y coordinates
+int npix[ARB] #O the output area in number of pixels
+real mag[ARB] #O the output magnitude estimates
+real size[ARB] #O the output size estimates
+real ellip[ARB] #O the output ellipticity estimates
+real theta[ARB] #O the output position angle estimates
+real sharp[ARB] #O the output sharpness estimates
+int nobjs #I the number of objects
+real datamin, datamax #I the minium and maximum good data values
+real threshold #I threshold for moments computation
+real hwhmpsf #I the HWHM of the PSF
+real xoff, yoff #I the x and y coordinate offsets
+
+int i, j, k, xmiddle, ymiddle, sumn
+double pixval, sumix, sumiy, sumi, sumixx, sumixy, sumiyy, r2, dx, dy, diff
+double mean
+
+begin
+ # Initialize
+ xmiddle = 1 + nxk / 2
+ ymiddle = 1 + nyk / 2
+
+ # Compute the pixel sum, number of pixels, and the x and y centers.
+ do i = 1, nobjs {
+
+ # Estimate the background using the input data and the
+ # best fitting Gaussian amplitude
+ sumn = 0
+ sumi = 0.0
+ do j = 1, nyk {
+ do k = 1, nxk {
+ if (skip[k,j] == NO)
+ next
+ pixval = data[cols[i]-xmiddle+k,j]
+ if (pixval < datamin || pixval > datamax)
+ next
+ sumi = sumi + pixval
+ sumn = sumn + 1
+ }
+ }
+ if (sumn <= 0)
+ mean = data[cols[i],ymiddle] - den[cols[i],ymiddle]
+ else
+ mean = sumi / sumn
+
+ # Compute the first order moments.
+ sumi = 0.0
+ sumn = 0
+ sumix = 0.0d0
+ sumiy = 0.0d0
+ do j = 1, nyk {
+ do k = 1, nxk {
+ if (skip[k,j] == YES)
+ next
+ pixval = data[cols[i]-xmiddle+k,j]
+ if (pixval < datamin || pixval > datamax)
+ next
+ pixval = pixval - mean
+ if (pixval <= 0.0)
+ next
+ sumi = sumi + pixval
+ sumix = sumix + (cols[i] - xmiddle + k) * pixval
+ sumiy = sumiy + j * pixval
+ sumn = sumn + 1
+ }
+
+ }
+
+ # Use the first order moments to estimate the positions
+ # magnitude, area, and amplitude of the object.
+ if (sumi <= 0.0) {
+ x[i] = cols[i]
+ y[i] = (1.0 + nyk) / 2.0
+ mag[i] = INDEFR
+ npix[i] = 0
+ } else {
+ x[i] = sumix / sumi
+ y[i] = sumiy / sumi
+ mag[i] = -2.5 * log10 (sumi)
+ npix[i] = sumn
+ }
+
+ # Compute the second order central moments using the results of
+ # the first order moment analysis.
+ sumixx = 0.0d0
+ sumiyy = 0.0d0
+ sumixy = 0.0d0
+ do j = 1, nyk {
+ dy = j - y[i]
+ do k = 1, nxk {
+ if (skip[k,j] == YES)
+ next
+ pixval = data[cols[i]-xmiddle+k,j]
+ if (pixval < datamin || pixval > datamax)
+ next
+ pixval = pixval - mean
+ if (pixval <= 0.0)
+ next
+ dx = cols[i] - xmiddle + k - x[i]
+ sumixx = sumixx + pixval * dx ** 2
+ sumixy = sumixy + pixval * dx * dy
+ sumiyy = sumiyy + pixval * dy ** 2
+ }
+ }
+
+ # Use the second order central moments to estimate the size,
+ # ellipticity, position angle, and sharpness of the objects.
+ if (sumi <= 0.0) {
+ size[i] = 0.0
+ ellip[i] = 0.0
+ theta[i] = 0.0
+ sharp[i] = INDEFR
+ } else {
+ sumixx = sumixx / sumi
+ sumixy = sumixy / sumi
+ sumiyy = sumiyy / sumi
+ r2 = sumixx + sumiyy
+ if (r2 <= 0.0) {
+ size[i] = 0.0
+ ellip[i] = 0.0
+ theta[i] = 0.0
+ sharp[i] = INDEFR
+ } else {
+ size[i] = sqrt (LN_2 * r2)
+ sharp[i] = size[i] / hwhmpsf
+ diff = sumixx - sumiyy
+ ellip[i] = sqrt (diff ** 2 + 4.0d0 * sumixy ** 2) / r2
+ if (diff == 0.0d0 && sumixy == 0.0d0)
+ theta[i] = 0.0
+ else
+ theta[i] = RADTODEG (0.5d0 * atan2 (2.0d0 * sumixy,
+ diff))
+ if (theta[i] < 0.0)
+ theta[i] = theta[i] + 180.0
+ }
+ }
+
+ # Convert the computed coordinates to the image system.
+ x[i] = x[i] + xoff
+ y[i] = y[i] + yoff
+ }
+end
+
+
+# SF_TEST -- Check that the detected objects are in the image, contain
+# enough pixels above background to be measurable objects, and are within
+# the specified magnitude, roundness and sharpness range.
+
+int procedure sf_test (cols, x, y, npix, mag, size, ellip, theta, sharps,
+ nobjs, c1, c2, l1, l2, nmin, maglo, maghi, roundlo, roundhi,
+ sharplo, sharphi)
+
+int cols[ARB] #U the column ids of detected object
+real x[ARB] #U the x position estimates
+real y[ARB] #U the y positions estimates
+int npix[ARB] #U the area estimates
+real mag[ARB] #U the magnitude estimates
+real size[ARB] #U the size estimates
+real ellip[ARB] #U the ellipticity estimates
+real theta[ARB] #U the position angle estimates
+real sharps[ARB] #U sharpness estimates
+int nobjs #I the number of detected objects
+real c1, c2 #I the image column limits
+real l1, l2 #I the image line limits
+int nmin #I the minimum area
+real maglo, maghi #I the magnitude limits
+real roundlo, roundhi #I the roundness limits
+real sharplo, sharphi #I the sharpness limits
+
+int i, nstars
+
+begin
+ # Loop over the detected objects.
+ nstars = 0
+ do i = 1, nobjs {
+
+ if (x[i] < c1 || x[i] > c2)
+ next
+ if (y[i] < l1 || y[i] > l2)
+ next
+ if (npix[i] < nmin)
+ next
+ if (mag[i] < maglo || mag[i] > maghi)
+ next
+ if (ellip[i] < roundlo || ellip[i] > roundhi)
+ next
+ if (! IS_INDEFR(sharps[i]) && (sharps[i] < sharplo ||
+ sharps[i] > sharphi))
+ next
+
+ # Add object to the list.
+ nstars = nstars + 1
+ cols[nstars] = cols[i]
+ x[nstars] = x[i]
+ y[nstars] = y[i]
+ mag[nstars] = mag[i]
+ npix[nstars] = npix[i]
+ size[nstars] = size[i]
+ ellip[nstars] = ellip[i]
+ theta[nstars] = theta[i]
+ sharps[nstars] = sharps[i]
+ }
+
+ return (nstars)
+end
+
+
+# SF_WRITE -- Write the results to the output file.
+
+procedure sf_write (fd, cols, x, y, mag, npix, size, ellip, theta, sharp,
+ nstars, ct, fmtstr, stid)
+
+int fd #I the output file descriptor
+int cols[ARB] #I column numbers
+real x[ARB] #I xcoords
+real y[ARB] #I y coords
+real mag[ARB] #I magnitudes
+int npix[ARB] #I number of pixels
+real size[ARB] #I object sizes
+real ellip[ARB] #I ellipticities
+real theta[ARB] #I position angles
+real sharp[ARB] #I sharpnesses
+int nstars #I number of detected stars in the line
+pointer ct #I coordinate transformation
+char fmtstr[ARB] #I the output format string
+int stid #I output file sequence number
+
+double lx, ly, wx, wy
+int i
+
+begin
+ if (fd == NULL)
+ return
+
+ do i = 1, nstars {
+ call fprintf (fd, fmtstr)
+ call pargr (x[i])
+ call pargr (y[i])
+ if (ct != NULL) {
+ lx = x[i]
+ ly = y[i]
+ call mw_c2trand (ct, lx, ly, wx, wy)
+ call pargd (wx)
+ call pargd (wy)
+ }
+ call pargr (mag[i])
+ call pargi (npix[i])
+ call pargr (size[i])
+ call pargr (ellip[i])
+ call pargr (theta[i])
+ call pargr (sharp[i])
+ call pargi (stid + i - 1)
+ }
+end
generated by cgit (git 2.34.1) at 2025-09-20 05:25:43 -0400