include <imhdr.h>

define	SZ_KEYW		11	# size of buffer for header keyword

# fsquares -- compute sum of squares
# This task computes the sum of squares and normalizes by the square of
# the sum of the values as a measure of the sharpness of an image.
# The output from fextract is the input to this routine.
#
# Phil Hodge,  7-Feb-1990  Task created.

procedure fsquares()

pointer input			# scr for name of input image
pointer output			# scr for name of output ascii file
int	box1			# separation betw different stars in a slice
int	box2			# separation betw slices at diff focus positions
#--
pointer sp
pointer rownum			# scratch for list of row numbers
pointer im			# pointer to imhdr struct
pointer v			# pointer to input data
char	keyword[SZ_KEYW]	# keyword for getting refoc position
real	refoc			# refoc mirror position
real	sharp			# sharpness:  sum of squares / sum of values
int	fd			# fd for output file
int	numimages		# number of focus positions
int	naxis1			# width of input image
int	y1, y2			# limits of image section on second axis
int	k			# loop index for focus position
pointer immap(), imgs2r()
int	open(), imgeti(), imaccf()
int	findrow()
real	imgetr()

begin
	call smark (sp)
	call salloc (input, SZ_FNAME, TY_CHAR)
	call salloc (output, SZ_FNAME, TY_CHAR)
	call clgstr ("input", Memc[input], SZ_FNAME)
	call clgstr ("output", Memc[output], SZ_FNAME)

	im = immap (Memc[input], READ_ONLY, NULL)
	fd = open (Memc[output], NEW_FILE, TEXT_FILE)

	# Get the list of numbers of rows that have been filled in.
	call salloc (rownum, SZ_LINE, TY_CHAR)
	call imgstr (im, "rownum", Memc[rownum], SZ_LINE)

	# fextract puts different stars along the horizontal direction
	# and different focus positions (images input to fextract) in
	# the vertical direction.  The slice we are going to take from
	# the input image is of size naxis1 x box2.
	naxis1 = IM_LEN(im,1)
	box1 = imgeti (im, "box1")
	box2 = imgeti (im, "box2")
	numimages = IM_LEN(im,2) / box2		# should divide exactly

	do k = 1, numimages {		# do for each focus position

	    # Has this row been filled in?  If not, skip it.
	    if (findrow (Memc[rownum], k) == NO)
		next

	    # Get the refoc mirror position from the input header.
	    call sprintf (keyword, SZ_KEYW, "refoc%d")
		call pargi (k)
	    if (imaccf (im, keyword) == YES)
		refoc = imgetr (im, keyword)
	    else
		refoc = 0.

	    # Get the data (a horizontal slice of size naxis1 x box2)
	    # for all stars for a given focus position.
	    y1 = (k - 1) * box2 + 1
	    y2 = k * box2
	    v = imgs2r (im, 1, naxis1, y1, y2)

	    # Compute the sharpness of this image section.
	    call fs_sumsq (Memr[v], naxis1, box2, sharp)

	    # Write the refoc position and sharpness measure to the file.
	    call fprintf (fd, "%15.7g %15.7g\n")
		call pargr (refoc)
		call pargr (sharp)
	    call flush (fd)
	}

	call close (fd)
	call imunmap (im)
	call sfree (sp)
end


# fs_sumsq -- compute sum of squares
# This routine computes a measure of the sharpness of an image:
# the sum of squares of the values divided by the square of the sum of values.

procedure fs_sumsq (v, naxis1, box2, sharp)

real	v[naxis1,box2]		# i: input section
int	naxis1, box2		# i: size of input
real	sharp			# o: normalized sum of squares
#--
real	sumsq			# for accumulating sum of squares
real	sum			# for accumulating sum of values
int	j, k			# loop indexes

begin
	sum = 0.
	sumsq = 0.
	do k = 1, box2 {
	    do j = 1, naxis1 {
		sum = sum + v[j,k]
		sumsq = sumsq + v[j,k] * v[j,k]
	    }
	}

	if (sum == 0.)
	    sharp = 0.
	else
	    sharp = sumsq / (sum * sum)
end