Repatch (from linux) of OSX IRAF

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