path: root/noao/digiphot/apphot/daofind/apfind.x

include <gset.h>
include <mach.h>
include <imhdr.h>
include "../lib/apphot.h"

# AP_FIND -- Detect images in the convolved image and then compute image
# characteristics using the original image.

int procedure ap_find (ap, im, cnv, out, id, ker2d, skip, nxk, nyk, skymode,
	threshold, relerr, emission, xsigsq, ysigsq, datamin, datamax,
	sharplo, sharphi, roundlo, roundhi, interactive, stid, mkdetections)

pointer ap			# the apphot descriptor
pointer	im			# pointer to the input image
pointer	cnv			# pointer to the output image
int	out			# the output file descriptor
pointer	id			# pointer to the display stream
real	ker2d[nxk,ARB]		# 2D Gaussian kernel
int	skip[nxk,ARB]		# 2D skip kernel
int	nxk, nyk		# dimensions of the kernel
real	skymode			# estimate of the sky
real	threshold		# threshold for image detection
real	relerr			# the relative error of the convolution kernel
int	emission		# emission features
real	xsigsq, ysigsq		# sigma of gaussian in x and y
real	datamin, datamax	# minimum and maximum good data values
real	sharplo, sharphi	# sharpness limits
real	roundlo,roundhi		# roundness parameter limits
int	interactive		# interactive mode
int	stid			# sequence number
int	mkdetections		# mark detections

int	inline, i, j, ncols, col1, col2, line1, line2, index, pos
int	xmiddle, ymiddle, nonzero, nobjs, nstars, ntotal
pointer	sp, bufptrs, imlbuf, cnvlbuf, imbuf, cnvbuf, cols
pointer	satur, sharp, round1, round2, x, y

int	ap_detect(), ap_test(), apstati()
pointer	imgs2r()
errchk	imgs2r()

begin
	# Set up useful line and column limits.
	ncols = IM_LEN(im,1) + nxk - 1
	col1 = 1 - nxk / 2
	col2 = IM_LEN(im,1) + nxk / 2
	line1 = 1 + nyk / 2
	line2 = IM_LEN(im,2) + nyk / 2
	xmiddle = 1 + nxk / 2
	ymiddle = 1 + nyk / 2

	# Compute find the number of defined elements in the kernel.
	nonzero = 0
	#skip[xmiddle,ymiddle] = NO
	do j = 1, nyk {
	    do i = 1, nxk {
		if (skip[i,j] == NO)
		    nonzero = nonzero + 1
	    }
	}
	skip[xmiddle,ymiddle] = YES
	nonzero = nonzero - 1

	# Set up a cylindrical buffers and some working space for
	# the detected images.
	call smark (sp)
	call salloc (bufptrs, nyk, TY_INT)
	call salloc (imbuf, nyk * ncols, TY_REAL)
	call salloc (cnvbuf, nyk * ncols, TY_REAL)
	call salloc (cols, ncols, TY_INT)
	call salloc (satur, ncols, TY_INT)
	call salloc (sharp, ncols, TY_REAL)
	call salloc (round1, ncols, TY_REAL)
	call salloc (round2, ncols, TY_REAL)
	call salloc (x, ncols, TY_REAL)
	call salloc (y, ncols, TY_REAL)

	# Read in the first nyk - 1 lines.
	pos = nyk
	do inline = 1 - nyk / 2, nyk / 2 {
	    imlbuf = imgs2r (im, col1, col2, inline, inline)
	    cnvlbuf = imgs2r (cnv, col1, col2, inline, inline)
	    if (emission == YES) {
	        call amovr (Memr[imlbuf], Memr[imbuf+(inline+ymiddle-2)*ncols],
	            ncols)
	        call amovr (Memr[cnvlbuf], Memr[cnvbuf+(inline+ymiddle-2)*
		    ncols], ncols)
	    } else {
	        call amulkr (Memr[imlbuf], -1.0, Memr[imbuf+(inline+ymiddle-2)*
		    ncols], ncols)
	        call amulkr (Memr[cnvlbuf], -1.0, Memr[cnvbuf+(inline+ymiddle-
		    2)* ncols], ncols)
	    }
	    Memi[bufptrs+pos-1] = pos - 1
	    pos = pos - 1
	}

	# Generate the starlist line by line.
	ntotal = 0
	pos = nyk
	do inline = line1, line2 {

	    # Setup the buffer pointer array.
	    do j = 2, nyk
		Memi[bufptrs+j-2] = Memi[bufptrs+j-1]
	    Memi[bufptrs+nyk-1] = pos
	    index = (pos - 1) * ncols

	    # Read in new image line.
	    imlbuf = imgs2r (im, col1, col2, inline, inline)
	    cnvlbuf = imgs2r (cnv, col1, col2, inline, inline)

	    # Copy new lines into cylindrical buffer.
	    if (emission == YES) {
	        call amovr (Memr[imlbuf], Memr[imbuf+index], ncols)
	        call amovr (Memr[cnvlbuf], Memr[cnvbuf+index], ncols)
	    } else {
	        call amulkr (Memr[imlbuf], -1.0, Memr[imbuf+index], ncols)
	        call amulkr (Memr[cnvlbuf], -1.0, Memr[cnvbuf+index], ncols)
	    }

	    # Detect stars in each image line. In order for a given pixel
	    # to be detected as an image the pixel must be above threshold
	    # and be greater than any other pixel within nsigma sigma.

	    # Increment the cylindrical buffer.
	    if (mod (pos, nyk) == 0)
		pos = 1
	    else
		pos = pos + 1

	    nobjs = ap_detect (Memr[cnvbuf], Memi[bufptrs], ncols, skip, nxk,
	        nyk, relerr * threshold, Memi[cols])
	    if (nobjs <= 0)
		next

	    # Compute the sharpness parameter.
	    call ap_sharp_round (Memr[imbuf], Memr[cnvbuf], Memi[bufptrs],
	        ncols, skip, nxk, nyk, Memi[cols], Memi[satur], Memr[round1],
		Memr[sharp], nobjs, nonzero, skymode, datamin, datamax) 

	    # Compute the roundness parameters.
	    call ap_xy_round (Memr[imbuf], Memi[bufptrs], ncols, ker2d, nxk,
	         nyk, Memi[cols], inline, Memr[round2], Memr[x],
		 Memr[y], nobjs, skymode, datamin, datamax, xsigsq, ysigsq)

	    # Test the image characeteristics of detected objects.
	    nstars = ap_test (Memi[cols], Memr[x], Memr[y], Memi[satur],
		Memr[round1], Memr[round2], Memr[sharp], nobjs, IM_LEN(im,1),
		IM_LEN(im,2), sharplo, sharphi, roundlo, roundhi)

	    # Mark the stars on the display.
	    if ((nstars > 0) && (interactive == YES) && (id != NULL) &&
	        (mkdetections == YES)) {
		call greactivate (id, 0)
		do j = 1, nstars {
		    #call ap_ltov (im, Memr[x+j-1], Memr[y+j-1], xc, yc, 1)
		    #call gmark (id, xc, yc, GM_PLUS, 1.0, 1.0)
		    call gmark (id, Memr[x+j-1], Memr[y+j-1], GM_PLUS, 1.0, 1.0)
		}
		call gdeactivate (id, 0)
	    }

	    switch (apstati (ap, WCSOUT)) {
	    case WCS_PHYSICAL:
		call ap_ltoo (ap, Memr[x], Memr[y], Memr[x], Memr[y], nstars)
	    case WCS_TV:
		call ap_ltov (im, Memr[x], Memr[y], Memr[x], Memr[y], nstars)
	    default:
		;
	    }

	    # Print results on the standard output.
	    if (interactive == YES)
	        call apstdout (Memr[cnvbuf], Memi[bufptrs], ncols, nyk,
	            Memi[cols], Memr[x], Memr[y], Memr[sharp], Memr[round1],
		    Memr[round2], nstars, ntotal, relerr * threshold)

	    # Save the results in the file.
	    call apdtfout (out, Memr[cnvbuf], Memi[bufptrs], ncols, nyk,
	        Memi[cols], Memr[x], Memr[y], Memr[sharp], Memr[round1],
		Memr[round2], nstars, ntotal, relerr * threshold, stid)


	    ntotal = ntotal + nstars

	}

	# Free space
	call sfree (sp)

	return (ntotal)
end


# AP_DETECT -- Detec 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 nsigma * sigma.

int procedure ap_detect (density, ptrs, ncols, skip, nxk, nyk, threshold, cols)

real	density[ncols, ARB]	# density array
int	ptrs[ARB]		# pointer array
int	ncols			# x dimesnsion of intensity buffer
int	skip[nxk,ARB]		# skip array
int	nxk, nyk		# size of convolution kernel
real	threshold		# density threshold
int	cols[ARB]		# column numbers of detected stars

int	i, j, k, kk, middle, nhalf, nobjs
define	nextpix_	11

begin
	middle = 1 + nyk / 2
	nhalf = nxk / 2

	# Loop over all the columns in an image line.
	nobjs = 0
	for (i = 1 + nhalf; i <= ncols - nhalf; ) {

	    # Test whether the density enhancement is above threshold.
	    if (density[i,ptrs[middle]] < threshold)
		goto nextpix_

	    # Test whether a given density enhancement is a local maximum.
	    do j = 1, nyk {
		kk = 1
	        do k = i - nhalf, i + nhalf {
		    if (skip[kk,j] == NO) {
		        if (density[i,ptrs[middle]] < density[k,ptrs[j]])
		           goto nextpix_
		    }
		    kk = kk + 1
	        }
	    }

	    # 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 + nhalf.
	    i = i + nhalf
nextpix_
	    # Work on the next pixel.
	    i = i + 1
	}

	return (nobjs)
end


# AP_SHARP_ROUND -- Compute an estimate of the roundness and sharpness of the
# detected objects. The roundness parameter is computed by comparing a measure
# of the bilateral symmetry with a measure of the four-fold symmetry. The
# sharpness parameter is defined as the ratio of the difference between the
# height of the central pixel and the mean of the surrounding pixels to the
# density enhancement of the central pixel.

procedure ap_sharp_round (data, density, ptrs, ncols, skip, nxk, nyk, cols,
        satur, round, sharps, nobjs, nonzero, skymode, datamin, datamax) 

real	data[ncols,ARB]		# image data
real	density[ncols,ARB]	# density enhancements
int	ptrs[ARB]		# buffer pointers
int	ncols			# length of data array
int	skip[nxk,ARB]		# 2D kernel
int	nxk, nyk		# size of convolution kernel
int	cols[ARB]		# array of columns
int	satur[ARB]		# array of saturated state parameters
real	round[ARB]		# array of roundness parameters
real	sharps[ARB]		# array of sharpness parameters
int	nobjs			# number of objects
int	nonzero			# number of nonzero kernel elements
real	skymode			# estimate of the sky mode
real	datamin, datamax	# minimum and maximum good data values

int	i, j, k, xmiddle, ymiddle, npixels, nhalf
real	pixval, midpix, temp, sharp, sum2, sum4

begin
	# Loop over the detected objects.
	nhalf = min (nxk / 2, nyk / 2)
	xmiddle = 1 + nxk / 2 
	ymiddle = 1 + nyk / 2
	do i = 1, nobjs {

	    # Compute the first estimate of roundness.
	    sum2 = 0.0
	    sum4 = 0.0
	    do k = 0, nhalf {
		do j = 1, nhalf {
		    sum2 = sum2 +
		        density[cols[i]-k,ptrs[ymiddle-j]] +
		        density[cols[i]+k,ptrs[ymiddle+j]] -
		        density[cols[i]-j,ptrs[ymiddle+k]] -
		        density[cols[i]+j,ptrs[ymiddle-k]]
		    sum4 = sum4 +
		        abs (density[cols[i]-k,ptrs[ymiddle-j]]) +
		        abs (density[cols[i]+k,ptrs[ymiddle+j]]) +
		        abs (density[cols[i]-j,ptrs[ymiddle+k]]) +
		        abs (density[cols[i]+j,ptrs[ymiddle-k]])
		}
	    }
	    if (sum2 == 0.0)
		round[i] = 0.0
	    else if (sum4 <= 0.0)
		round[i] = INDEFR
	    else
	        round[i] = 2.0 * sum2 / sum4

	    satur[i] = NO

	    # Eliminate the sharpness test if the central pixel is bad.
	    midpix = data[cols[i],ptrs[ymiddle]]
	    if (midpix > datamax) {
		satur[i] = YES
		sharps[i] = INDEFR
		next
	    }
	    if (midpix < datamin) {
		sharps[i] = INDEFR
		next
	    }

	    # Accumulate the sharpness statistic.
	    sharp = 0.0
	    npixels = nonzero
	    do j = 1, nyk {
		temp = 0.0
	        do k = 1, nxk {
		    if (skip[k,j] == YES)
			next
		    pixval = data[cols[i]-xmiddle+k,ptrs[j]]
		    if (pixval > datamax) {
			satur[i] = YES
			npixels = npixels - 1
		    } else if (pixval < datamin) {
			npixels = npixels - 1
		    } else {
		        temp = temp + (pixval - skymode)
		    }
		}
		sharp = sharp + temp
	    }

	    # Compute the sharpness statistic.
	    if (density[cols[i],ptrs[ymiddle]] <= 0.0 || npixels <= 0)
		sharps[i] = INDEFR
	    else
	        sharps[i] = (midpix - skymode - sharp / real (npixels)) /
		    density[cols[i],ptrs[ymiddle]]

	}
end


# AP_XY_ROUND -- Estimate the x-y centers and the roundness of the detected
# objects. The height of the equivalent Gaussian function in x and y is fit by
# least squares to the marginal distribution of the image data. If either
# of these of these heights is negative set the roundess characteristic to
# -MAX_REAL, otherwise compute a roundness characteristic. At the same
# time setup the necessary sums for computing the first order corection
# to the centroid of the gaussian profile.

procedure ap_xy_round (data, ptrs, ncols, ker2d, nxk, nyk, cols, inline,
	rounds, x, y, nobjs, skymode, datamin, datamax, xsigsq, ysigsq)

real	data[ncols,ARB]		# density enhancements
int	ptrs[ARB]		# buffer pointers
int	ncols			# number of columns in cylindrical buffer
real	ker2d[nxk,ARB]		# the gaussian convolution kernel
int	nxk			# size of kernel in x
int	nyk			# size of kernel in y
int	cols[ARB]		# the input positions	
int	inline			# the input image line
real	rounds[ARB]		# array of sharpness parameters
real	x[ARB]			# output x coords
real	y[ARB]			# output y coords
int	nobjs			# number of objects
real	skymode			# estimate of the sky mode
real	datamin, datamax	# minium and maximum data values
real	xsigsq, ysigsq		# x-y gaussian sigma squared

int	i, j, k, xmiddle, ymiddle, n
real	sumgd, sumgsq, sumg, sumd, sumdx, dgdx, sdgdx, sdgdxsq, sddgdx, sgdgdx
real	pixval, p, sg, sd, wt, hx, hy, dx, dy, skylvl, xhalf, yhalf

begin
	xhalf = real (nxk / 2) + 0.5
	yhalf = real (nyk / 2) + 0.5
	xmiddle = 1 + nxk / 2
	ymiddle = 1 + nyk / 2

	# Loop over the detected objects.
	do i = 1, nobjs {

	    # Initialize the x fit.
	    sumgd = 0.0
	    sumgsq = 0.0
	    sumg = 0.0
	    sumd = 0.0
	    sumdx = 0.0
	    sdgdx = 0.0
	    sdgdxsq = 0.0
	    sddgdx = 0.0
	    sgdgdx = 0.0
	    p = 0.0
	    n = 0

	    # Compute the sums required for the x fit.
	    do k = 1, nxk {

		sg = 0.0
		sd = 0.0
		do j = 1, nyk {
		    wt = real (ymiddle - abs (j - ymiddle))
		    pixval = data[cols[i]-xmiddle+k,ptrs[j]]
		    if (pixval < datamin || pixval > datamax)
			next
		    sd = sd + (pixval - skymode) * wt
		    sg = sg + ker2d[k,j] * wt
		}

	        if (sg <= 0.0)
		    next
	        wt = real (xmiddle - abs (k - xmiddle))
	        sumgd = sumgd + wt * sg * sd
	        sumgsq = sumgsq + wt * sg ** 2
	        sumg = sumg + wt * sg
	        sumd = sumd + wt * sd
		sumdx = sumdx + wt * sd * (xmiddle - k)
	        p = p + wt
	        n = n + 1
	        dgdx = sg * (xmiddle - k)
	        sdgdxsq = sdgdxsq + wt * dgdx ** 2
	        sdgdx = sdgdx + wt * dgdx
	        sddgdx = sddgdx + wt * sd * dgdx
	        sgdgdx = sgdgdx + wt * sg * dgdx
	    }

	    # Need at least three points to estimate the x height, position
	    # and local sky brightness of the star.

	    if (n <= 2 || p <= 0.0) {
		x[i] = INDEFR
		y[i] = INDEFR
		rounds[i] = INDEFR
		next
	    }

	    # Solve for the height of the best-fitting gaussian to the
	    # xmarginal. Reject the star if the height is non-positive.

	    hx = sumgsq - (sumg ** 2) / p
	    if (hx <= 0.0) {
		x[i] = INDEFR
		y[i] = INDEFR
		rounds[i] = INDEFR
		next
	    }
	    hx = (sumgd - sumg * sumd / p) / hx
	    if (hx <= 0.0) {
		x[i] = INDEFR
		y[i] = INDEFR
		rounds[i] = INDEFR
		next
	    }

	    # Solve for the new x centroid.
	    skylvl = (sumd - hx * sumg) / p
	    dx = (sgdgdx - (sddgdx - sdgdx * (hx * sumg + skylvl * p))) /
		(hx * sdgdxsq / xsigsq)
	    if (abs (dx) > xhalf) {
		if (sumd == 0.0)
		    dx = 0.0
		else
		    dx = sumdx / sumd 
		if (abs (dx) > xhalf)
		    dx = 0.0
	    }
	    x[i] = (cols[i] - xmiddle + 1) + dx 

	    # Initialize y fit.
	    sumgd = 0.0
	    sumgsq = 0.0
	    sumg = 0.0
	    sumd = 0.0
	    sumdx = 0.0
	    sdgdx = 0.0
	    sdgdxsq = 0.0
	    sddgdx = 0.0
	    sgdgdx = 0.0
	    p = 0.0
	    n = 0

	    do j = 1, nyk {
		sg = 0.0
		sd = 0.0
		do k = 1, nxk {
		    wt = real (xmiddle - abs (k - xmiddle))
		    pixval = data[cols[i]-xmiddle+k,ptrs[j]]
		    if (pixval < datamin || pixval > datamax)
			next
		    sd = sd + (pixval - skymode) * wt
		    sg = sg + ker2d[k,j] * wt
		}
	        if (sg <= 0.0)
		    next
	        wt = real (ymiddle - abs (j - ymiddle))
	        sumgd = sumgd + wt * sg * sd
	        sumgsq = sumgsq + wt * sg ** 2
	        sumg = sumg + wt * sg
	        sumd = sumd + wt * sd
		sumdx = sumdx + wt * sd * (j - ymiddle)
	        p = p + wt
	        n = n + 1
	        dgdx = sg * (ymiddle - j)
	        sdgdx = sdgdx + wt * dgdx
	        sdgdxsq = sdgdxsq + wt * dgdx ** 2
	        sddgdx = sddgdx + wt * sd * dgdx
	        sgdgdx = sgdgdx + wt * sg * dgdx
	    }

	    # Need at least three points to estimate the y height, position
	    # and local sky brightness of the star.

	    if (n <= 2 || p <= 0.0) {
		x[i] = INDEFR
		y[i] = INDEFR
		rounds[i] = INDEFR
		next
	    }

	    # Solve for the height of the best-fitting gaussian to the
	    # y marginal. Reject the star if the height is non-positive.

	    hy = sumgsq - (sumg ** 2) / p 
	    if (hy <= 0.0) {
		x[i] = INDEFR
		y[i] = INDEFR
		rounds[i] = INDEFR
		next
	    }
	    hy = (sumgd - sumg * sumd / p) / (sumgsq - (sumg ** 2) / p)
	    if (hy <= 0.0) {
		x[i] = INDEFR
		y[i] = INDEFR
		rounds[i] = INDEFR
		next
	    }

	    # Solve for the new x centroid.
	    skylvl = (sumd - hy * sumg) / p
	    dy = (sgdgdx - (sddgdx - sdgdx * (hy * sumg + skylvl * p))) /
		(hy * sdgdxsq / ysigsq)
	    if (abs (dy) > yhalf) {
		if (sumd == 0.0)
		    dy = 0.0
		else
		    dy = sumdx / sumd 
		if (abs (dy) > yhalf)
		    dy = 0.0
	    }
	    y[i] = (inline - ymiddle + 1) + dy

	    # Compute the roundness.
	    rounds[i] = 2.0 * (hx - hy) / (hx + hy)
	}
end


# AP_TEST -- Test the characteristic of the detected images for roundness
# and sharpness.

int procedure ap_test (cols, x, y, satur, round1, round2, sharps, nobjs,
	ncols, nlines, sharplo, sharphi, roundlo, roundhi)

int	cols[ARB]			# col IDS of detected images
real	x[ARB]				# x positions
real	y[ARB]				# y positions
int	satur[ARB]			# saturation condition
real	round1[ARB]			# first roundness parameters
real	round2[ARB]			# second roundness parameters
real	sharps[ARB]			# sharpness parameters
int	nobjs				# number of objects
int	ncols, nlines			# size of the input image
real	sharplo, sharphi		# sharpness parameters
real	roundlo, roundhi		# roundness parameters

int	i, nstars

begin
	# Loop over the detected objects.
	nstars = 0
	do i = 1, nobjs {

	    # Compute the sharpness statistic
	    if (! IS_INDEFR(sharps[i]) && (sharps[i] < sharplo ||
	        sharps[i] > sharphi))
		next
	    if (IS_INDEFR(round1[i]) || round1[i] < roundlo ||
	        round1[i] > roundhi)
		next
	    if (satur[i] == NO) {
	        if (IS_INDEFR(round2[i]) || round2[i] < roundlo ||
	            round2[i] > roundhi)
		next
	    }
	    if (IS_INDEFR(x[i]) || x[i] < 0.5 || x[i] > (ncols+0.5))
		next
	    if (IS_INDEFR(y[i]) || y[i] < 0.5 || y[i] > (nlines+0.5))
		next

	    # Add object to the list.
	    nstars = nstars + 1
	    cols[nstars] = cols[i]
	    x[nstars] = x[i]
	    y[nstars] = y[i]
	    sharps[nstars] = sharps[i]
	    round1[nstars] = round1[i]
	    round2[nstars] = round2[i]
	}

	return (nstars)
end