/*******************************************************************************
FIT -- Fit resolution templates to images with Poisson likelihood

A set of templates consisting of the PSF with different spatial scales,
two component composite templates of varying fraction of the PSF and
an extended component of the PSF at twice the scale, and three noise
templates are fit to determine the maximum Poisson likelihood.
The likelihood function is:

	L = sum (Qi - Ti - Qi * ln (Qi / Ti))

where the sum is over all pixels i in the image.  Qi is given by

	Qi = v0 + v1 * Di

where Di is the actual data value and v0 and v1 are parameters describing
the Poisson variance function.  Thus Qi is also the variance of the data
value.  Ti is the model template value given by

	Ti = I0 * ti + sky

where ti is the normalized and sky subtracted template, I0 is a scaling
determined for each object and sky is the sky of each object.

The fitting has the following steps:

	For each template determine I0 which maximizes L(ti)
	Find the template which has the maximum L

For efficiency not all templates are computed.  It is assumed that the
likelihood values vary smoothly for templates in which one parameter
varies, i.e. the scale parameter or the fraction parameter in the two
component templates and that the maximum can be found by a binary search.

There are special steps to allow DARK objects to be classified and to
deal with the situation where the Qi are zero or negative (i.e. sky <= 0.). 
*******************************************************************************/
#include	"focas1.h"

#define	NTMPS	104	/* Maximum number of templates */

static	int	DEBUG1 = 0;
static	int	DEBUG2 = 0;
static	int	DEBUG3 = 0;
static	int	DEBUG4 = 0;

static  struct	image    p;
static	int	param;
static	short   rct[4];
static	int	scale1, scale2, comp1, comp2, noise1, noise2, star, gal;
static	float	data[NYMAX][NXMAX], model[NTMPS][NYMAX][NXMAX];
static	float	params[NTMPS], prob[NTMPS], sky;

#define	ANINT(x)	((int) ((x < 0.) ? (x-0.5) : (x+0.5)))

/* FIT -- Fit resolution templates.  Return 1 if an error. */

fit (ob, area, type, center, v0, v1, tmpflag)

struct	objrec *ob;
struct	areas *area;
int	type;
int	center;
float	v0, v1;		/* Noise paramsters */
int	*tmpflag;

{
	int	i;
	float	p, fitmodel();

	/* Check object has been evaluated and is not saturated.
	/* Get the templates and return 1 if PSF is not found.
	/* Get the data and return 1 if error. */

	if (!(ob->eflgs & EVAL) || (ob->eflgs & PEAK))
	    return (0);
	if (getmodels (type, tmpflag) || getimage (ob, center))
	    return (1);

	/* Initialize */
	param = -1;
	for (i = 0; i < NTMPS; i++)
	    prob[i] = 1.;
	ob->eflgs &= ~CLSF;

	if (center > 0)
	    searchc (ob, area, center, star, v0, v1);
	else {
	    ob->fitxc = ANINT (ob->xc);
	    ob->fityc = ANINT (ob->yc);
	    getdata (ob, area, v0, v1);
	}

	/* Find the maximum likelihood over the scale and fraction templates.
	/* If the likelihood values are not well behaved for the binary
	/* search method evaluate all templates and set CLSF flag. */

	if (search1 (ob)) {
	    for (i = scale1; i <= scale2; i++)
	        p = fitmodel (ob, i);
	    ob->eflgs |= CLSF;
	}

	if (ob->eflgs & CLSF || search2 (ob)) {
	    for (i = comp1; i <= comp2; i++)
	        p = fitmodel (ob, i);
	    ob->eflgs |= CLSF;
	}

	/* Do noise templates */
	if (ob->eflgs & CLSF || params[param] < 0.8)
	    for (i = noise1; i <= noise2; i++)
	        p = fitmodel (ob, i);

	/* Set the catalog parameters for the maximum likelihood template */
	if ((param >= scale1) && (param <= scale2)) {
	    ob->scale = params[param];
	    ob->frac = 1.;
	} else if ((param >= comp1) && (param <= comp2)) {
	    ob->scale = params[gal];
	    ob->frac = params[param];
	} else if ((param >= noise1) && (param <= noise2)) {
	    ob->scale = params[param];
	    ob->frac = 1.;
	}

	return (0);
}	


/* SEARCHC - Search for maximum likelihood PSF center */

searchc (ob, area, center, k, v0, v1)
struct	objrec	*ob;
struct	areas   *area;
int	center;
int	k;
float	v0, v1;
{
	float	p, pmax;
	int	dx, dy, dxmax, dymax;
	float	fitmodel();

	pmax = -1e30;
	for (dy = -center; dy <= center; dy++) {
	    ob->fityc = ANINT (ob->yc + dy);
	    for (dx = -center; dx <= center; dx++) {
		if (dx == 0 && dy == 0)
		    continue;
	        ob->fitxc = ANINT (ob->xc + dx);
		getdata (ob, area, v0, v1);
		prob[k] = 1.;
		p = fitmodel (ob, k);
		if (p > pmax) {
		    dxmax = dx;
		    dymax = dy;
		    pmax = p;
		}
	    }
	}

        ob->fitxc = ANINT (ob->xc);
	ob->fityc = ANINT (ob->yc);
	getdata (ob, area, v0, v1);
	prob[k] = 1.;
	p = fitmodel (ob, k);

	if (p < pmax) {
	    ob->fitxc = ANINT (ob->xc + dxmax);
	    ob->fityc = ANINT (ob->yc + dymax);
	    prob[k] = pmax;
	    getdata (ob, area, v0, v1);
	}
}


/*SEARCH1 - Search for maximum likelihood over scale templates */
/* Return 1 if binary search fails. */

search1 (ob)
struct	objrec	*ob;
{
	int	k, kmin, kmax, ka, kb, kc;
	float	p, pmin, pmax, pa, pb, pc;
	float	fitmodel();

	/* Check if scale templates are defined */
	if (scale2 - scale1 + 1 < 1)
	    return (0);

	/* Do the end points and the PSF */
	kmin = scale1;
	kmax = scale2;
	ka = kmin + 1;
	kb = star;
	kc = kmax - 1;

	pmin = fitmodel (ob, kmin);
	pa = fitmodel (ob, ka);
	pb = fitmodel (ob, kb);
	pc = fitmodel (ob, kc);
	pmax = fitmodel (ob, kmax);

	/* Test for smoothness of likelihood function and return error */
	k = 0;
	if ((pmin - pa) * (pa - pb) < 0.)
	    k++;
	if ((pa - pb) * (pb - pc) < 0.)
	    k++;
	if ((pb - pc) * (pc - pmax) < 0.)
	    k++;
	if (k > 1)
	    return (1);

	if (pmin > pa)
	    return (0);
	if (pmax > pc)
	    return (0);

	if ((pa > pb) && (pb > pc)) {
	    pc = pb; kc = kb;
	    pb = pa; kb = ka;
	    pa = pmin; ka = kmin;
	}
	if ((pa < pb) && (pb < pc)) {
	    pa = pb; kb = ka;
	    pb = pc; kb = kc;
	    pc = pmax; kc = kmax;
	}

	/* Continue the binary search.  Return error if binary search fails. */
	for (;;) {
	    if (kb - ka > 1) {
		if (kb - ka == 2)
		    k = ka + 1;
		else if (pa < pb)
		    k = (ka + kb + kb) / 3;
		else
		    k = (ka + ka + kb) / 3;
		/* k = (ka + kb) / 2; */
	        p = fitmodel (ob, k);
	        if (p < pa)
		    return (1);
	        if (p > pb) {
		    pc = pb; kc = kb;
		    pb = p; kb = k;
	        } else {
		    pa = p; ka = k;
	        }
	    }
	    if (kc - kb > 1) {
		if (kc - kb == 2)
		    k = kb + 1;
		else if (pb < pc)
		    k = (kb + kc + kc) / 3;
		else
		    k = (kb + kb + kc) / 3;
	        /* k = (kb + kc) / 2; */
	        p = fitmodel (ob, k);
	        if (p < pc)
		    return (1);
	        if (p > pb) {
		    pa = pb; ka = kb;
		    pb = p; kb = k;
	        } else {
		    pc = p; kc = k;
	        }
	    }
	    if ((kb - ka < 2) && (kc - kb < 2))
	        break;
	}
	return (0);
}


/* SEARCH2 -- Search for maximum likelihood over composite templates */
/* Return 1 if binary search fails. */

search2 (ob)
struct	objrec	*ob;
{
	int	k, kmin, kmax, ka, kb, kc;
	float	p, pmin, pmax, pa, pb, pc;
	float	fitmodel();

	/* Check if fraction templates are defined */
	if (comp2 - comp1 + 1 < 1)
	    return (0);

	/* Test the endpoints and middle of range */
	kmin = star;
	kmax = gal;
	ka = comp1;
	kc = comp2;
	kb = (ka + kc) / 2;

	pmin = fitmodel (ob, kmin);
	pa = fitmodel (ob, ka);
	pb = fitmodel (ob, kb);
	pc = fitmodel (ob, kc);
	pmax = fitmodel (ob, kmax);

	/* Test for smoothness of likelihood function */
	k = 0;
	if ((pmin - pa) * (pa - pb) < 0.)
	    k++;
	if ((pa - pb) * (pb - pc) < 0.)
	    k++;
	if ((pb - pc) * (pc - pmax) < 0.)
	    k++;
	if (k > 1)
	    return (1);

	if (pmin > pa)
	    return (0);
	if (pmax > pc)
	    return (0);

	if ((pa > pb) && (pb > pc)) {
	    pc = pb; kc = kb;
	    pb = pa; kb = ka;
	    pa = pmin; ka = ka-1;
	}
	if ((pa < pb) && (pb < pc)) {
	    pa = pb; kb = ka;
	    pb = pc; kb = kc;
	    pc = pmax; kc = kc+1;
	}

	/* Continue the binary search.  Return error if binary search fails. */
	for (;;) {
	    if (kb - ka > 1) {
		if (kb - ka == 2)
		    k = ka + 1;
		else if (pa < pb)
		    k = (ka + kb + kb) / 3;
		else
		    k = (ka + ka + kb) / 3;
	        /* k = (ka + kb) / 2; */
	        p = fitmodel (ob, k);
	        if (p < pa)
		    return (1);
	        if (p > pb) {
		    pc = pb; kc = kb;
		    pb = p; kb = k;
	        } else {
		    pa = p; ka = k;
	        }
	    }
	    if (kc - kb > 1) {
		if (kc - kb == 2)
		    k = kb + 1;
		else if (pb < pc)
		    k = (kb + kc + kc) / 3;
		else
		    k = (kb + kb + kc) / 3;
	        /* k = (kb + kc) / 2; */
	        p = fitmodel (ob, k);
	        if (p < pc)
		    return (1);
	        if (p > pb) {
		    pa = pb; ka = kb;
		    pb = p; kb = k;
	        } else {
		    pc = p; kc = k;
	        }
	    }
	    if ((kb - ka < 2) && (kc - kb < 2))
	        break;
	}
	return (0);
}


/* FITMODEL -- Fit a model template to the data and return likelihood.
/* The intensity scale is determined iteratively to maximize the likelihood.
/* For efficiency the log is expanded for small values. */

float	fitmodel (ob, k)
struct	objrec *ob;		/* Object */
int	k;			/* Template index */
{
	int	i, j, iter;
	float	Imodel, I0, max, sum1, sum2, sum3, sum4, ratio, ratio1;
	float	d, m, *dp, *mp;

	/* If already determined simply return the value. */
	if (prob[k] <= 0.)
	    return (prob[k]);

	/* Determine the initial intensity scale */
	sum1 = sum2 = 0.;
	for (i = rct[1]; i <= rct[3]; i++) {
	    dp = data[i] + rct[0];
	    mp = model[k][i] + rct[0];
	    for (j = rct[0]; j <= rct[2]; j++) {
		d = *dp;
		if (d > 0.) {
		    m = *mp;
		    sum1 += m * (d - sky);
		    sum2 += m * m * d;
		}
		dp++;
		mp++;
	    }
	}
	I0 = sky * sum1 / sum2;

	/* Iteratively adjust the intensity scale to maximize the likelihood. */
	max = -1e30;
	for (iter = 0; iter < 6; iter++) {
	    sum1 = sum2 = sum3 = sum4 = 0.;
	    for (i = rct[1]; i <= rct[3]; i++) {
	        dp = data[i] + rct[0];
	        mp = model[k][i] + rct[0];
	        for (j = rct[0]; j <= rct[2]; j++) {
		    d = *dp;
		    if (d > 0.) {
			m = *mp;
		        Imodel = I0 * m + sky;
			ratio = d / Imodel;
			if (ratio < 0.001)
			    ratio = 0.001;
			else if (ratio > 1000.)
			    ratio = 1000.;
			ratio1 = Imodel / d - 1.;
			if (ratio1 > -.02 && ratio1 < .02)
			    sum1 += -0.5 * d * ratio1 * ratio1;
			/*
			    sum1 += d * ratio1 * ratio1 * (-0.5 + ratio1 / 3.);
			*/
			else
			    sum1 += d * (1. - log (ratio)) - Imodel;
		        sum2 += m * (ratio - 1);
		        sum3 += ratio * m * m / Imodel;
		    }
		    dp++;
		    mp++;
	        }
	    }
	    if (sum1 > max)
	        max = sum1;
	    sum3 = sum2 / sum3;

	    if (DEBUG4)
	        fprintf (stderr, "%d: I0=%g dI0=%g prob=%g\n",
		    k, I0, sum3, sum1);

	    /* Convergence criteria */
	    if (sum3 / I0 < .001)
	        break;
	    I0 = I0 + sum3;
	}

	/* Update the maximum over all templates */
	prob[k] = max;
	if (param < 0 || max > prob[param])
	    param = k;
	    
	if (DEBUG3)
	    fprintf (stderr, "%d %g - %d %g\n", k, prob[k], param, prob[param]);
	return (max);
}


/* GETMODELS -- Get resolution templates */

getmodels (type, tmpflag)
int	type, *tmpflag;
{
	int	i, j, k, xc, yc, xi, yi;
	static int tmpflg = 0;
	float	x, y, scale, lum, v1, v2, v3, v4;
	double	a, b, c, d;

	/* Use tmpflag to mark if the templates have already been set */
	if (*tmpflag == 1)
	    return (tmpflg);

	/* Return error if the PSF is not known. */
	if ((sp.nx == 0) || (sp.ny == 0)) {
	    if (*tmpflag == 0)
		fprintf (stderr, "getmodels: Point spread function undefined.\n");
	    tmpflg = 1;
	    return (tmpflg);
	}

	tmpflg = 0;
	*tmpflag = 1;

	/* Define the templates */
	switch (type) {
	case 1:
	    scale1 = 0;
	    scale2 = 81;
	    comp1 = 82;
	    comp2 = 100;
	    noise1 = 101;
	    noise2 = 103;
	    star = 15;
	    gal = 65;
	    for (i = scale1; i < gal; i++)
		params[i] = 0.7 + (i - scale1) * 0.02;
	    for (i = gal; i < scale2; i++)
		params[i] = 2.0 + (i - gal) * 0.5;
	    params[scale2] = 100.;
	    for (i = comp1; i <= comp2; i++)
		params[i] = 0.05 + (i - comp1) * 0.05;
	    for (i = noise1; i <= noise2; i++)
		params[i] = 0.1 + (i - noise1) * 0.1;
	    break;
	case 2:
	    scale1 = 0;
	    scale2 = 81;
	    comp1 = 82;
	    comp2 = 81;
	    noise1 = 82;
	    noise2 = 84;
	    star = 15;
	    gal = 65;
	    for (i = scale1; i < gal; i++)
		params[i] = 0.7 + (i - scale1) * 0.02;
	    for (i = gal; i < scale2; i++)
		params[i] = 2.0 + (i - gal) * 0.5;
	    params[scale2] = 100.;
	    for (i = comp1; i <= comp2; i++)
		params[i] = 0.05 + (i - comp1) * 0.05;
	    for (i = noise1; i <= noise2; i++)
		params[i] = 0.1 + (i - noise1) * 0.1;
	    break;
	case 3:
	    scale1 = 0;
	    scale2 = 35;
	    comp1 = 36;
	    comp2 = 44;
	    noise1 = 45;
	    noise2 = 47;
	    star = 6;
	    gal = 26;
	    for (i = scale1; i < gal; i++)
		params[i] = 0.7 + (i - scale1) * 0.05;
	    for (i = gal; i < scale2; i++)
		params[i] = 2.0 + (i - gal) * 1.0;
	    params[scale2] = 100.;
	    for (i = comp1; i <= comp2; i++)
		params[i] = 0.1 + (i - comp1) * 0.1;
	    for (i = noise1; i <= noise2; i++)
		params[i] = 0.1 + (i - noise1) * 0.1;
	    break;
	case 4:
	    scale1 = 0;
	    scale2 = 35;
	    comp1 = 36;
	    comp2 = 35;
	    noise1 = 36;
	    noise2 = 38;
	    star = 6;
	    gal = 26;
	    for (i = scale1; i < gal; i++)
		params[i] = 0.7 + (i - scale1) * 0.05;
	    for (i = gal; i < scale2; i++)
		params[i] = 2.0 + (i - gal) * 1.0;
	    params[scale2] = 100.;
	    for (i = comp1; i <= comp2; i++)
		params[i] = 0.1 + (i - comp1) * 0.1;
	    for (i = noise1; i <= noise2; i++)
		params[i] = 0.1 + (i - noise1) * 0.1;
	    break;
	default:
	    fprintf (stderr, "getmodels: Unknown template type.\n");
	    tmpflg = 1;
	    return (tmpflg);
	}

	xc = sp.nx / 2;
	yc = sp.ny / 2;

	/* Make scale templates by fitting a gaussian and then evaluating it
	/* at different scales and linearly interpolating the residuals. */
	fitg (&a, &b, &c, &d);
	for (k = scale1; k <= scale2; k++) {
	    scale = params[k];
	    lum = 0.;
	    for (i = 0; i < sp.ny; i++) {
		for (j = 0; j < sp.nx; j++) {
		    x = (j - xc) / scale + xc;
		    y = (i - yc) / scale + yc;
		    if ((x<0.)||(x>=sp.nx-1.)||(y<0.)||(y>=sp.ny-1.)) {
			model[k][i][j] = 0.;
			continue;
		    }
		    xi = x;
		    yi = y;
		    model[k][i][j] = exp (a - (b * (x - xc) * (x - xc) +
			c * (y - yc) * (y - yc) + d * (x - xc) * (y - yc)));
		    v1 = template[yi][xi];
		    v2 = template[yi][xi + 1];
		    v3 = template[yi + 1][xi];
		    v4 = template[yi + 1][xi + 1];
		    model[k][i][j] += v1 * (1 - (x - xi)) * (1 - (y - yi)) +
		    v2 * (x - xi) * (1 - (y - yi)) +
		    v3 * (1 - (x - xi)) * (y - yi) +
		    v4 * (x - xi) * (y - yi);
		    lum += model[k][i][j];
		}
	    }
	    for (i = 0; i < sp.ny; i++)
		for (j = 0; j < sp.nx; j++)
		    model[k][i][j] /= lum;
	}

	/* Make composite templates */
	for (k = comp1; k <= comp2; k++) {
	    v1 = params[k];
	    v2 = 1.- v1;
	    for (i = 0; i < sp.ny; i++)
	        for (j = 0; j < sp.nx; j++)
		    model[k][i][j] = v1 * model[gal][i][j] +
			v2 * model[star][i][j];
	}

	/* Make noise templates */
	for (lum = 0., i = 0; i < sp.ny; i++)
	    for (j = 0; j < sp.nx; j++) {
		model[noise1][i][j] = 1.;
		lum += 1.;
	    }
	model[noise1][yc][xc] = 100 * lum;
	lum += 100.* lum;
	for (i = 0; i < sp.ny; i++)
	    for (j = 0; j < sp.nx; j++)
		model[noise1][i][j] /= lum;
	for (lum = 0., i = 0; i < sp.ny; i++)
	    for (j = 0; j < sp.nx; j++) {
		model[noise1+1][i][j] = 1.;
		lum += 1.;
	    }
	model[noise1+1][yc][xc] = 20 * lum;
	model[noise1+1][yc - 1][xc] = 20 * lum;
	model[noise1+1][yc][xc - 1] = 20 * lum;
	model[noise1+1][yc + 1][xc] = 20 * lum;
	model[noise1+1][yc][xc + 1] = 20 * lum;
	lum += 100.* lum;
	for (i = 0; i < sp.ny; i++)
	    for (j = 0; j < sp.nx; j++)
 		model[noise1+1][i][j] /= lum;
	for (lum = 0., i = 0; i < sp.ny; i++)
	    for (j = 0; j < sp.nx; j++) {
    		model[noise1+2][i][j] = 1.;
    		lum += 1.;
	    }
	model[noise1+2][yc][xc] = 11 * lum;
	model[noise1+2][yc - 1][xc] = 11 * lum;
	model[noise1+2][yc][xc - 1] = 11 * lum;
	model[noise1+2][yc + 1][xc] = 11 * lum;
	model[noise1+2][yc][xc + 1] = 11 * lum;
	model[noise1+2][yc - 1][xc - 1] = 11 * lum;
	model[noise1+2][yc - 1][xc + 1] = 11 * lum;
	model[noise1+2][yc + 1][xc - 1] = 11 * lum;
	model[noise1+2][yc + 1][xc + 1] = 11 * lum;
	lum += 9.* 11.* lum;
	for (i = 0; i < sp.ny; i++)
	    for (j = 0; j < sp.nx; j++)
		model[noise1+2][i][j] /= lum;

	return (tmpflg);
}


/* GETIMAGE -- Get the image raster. */

getimage (ob, center)

struct	objrec  *ob;
int	center;

{
	rct[0] = ANINT (lpcoords (ob->xc,1,2)) - (sp.nx+2*center) / 2;
	rct[1] = ANINT (lpcoords (ob->yc,2,2)) - (sp.ny+2*center) / 2;
	rct[2] = ANINT (lpcoords (ob->xc,1,2)) + (sp.nx+2*center) / 2;
	rct[3] = ANINT (lpcoords (ob->yc,2,2)) + (sp.ny+2*center) / 2;
	if (rct[0] < 0)
	    rct[0] = 0;
	if (rct[1] < 0)
	    rct[1] = 0;
	if (rct[2] > pthdr.naxis1 - 1)
	    rct[2] = pthdr.naxis1 - 1;
	if (rct[3] > pthdr.naxis2 - 1)
	    rct[3] = pthdr.naxis2 - 1;
	if (rdimg (pfd, rct, &p) < 0)
	    return (1);
	return (0);
}


/* GETDATA -- Get the data to be fit.  Convert to "noise equivalent quanta"
/* and allow DARK objects.  Set the data array to zero outside the fitting
/* region. */

getdata (ob, a, v0, v1)

struct	objrec  *ob;
struct	areas   *a;
float	v0, v1;

{
	int	nx, ny, dx, dy, i, j, i1, j1, i2, j2;
	int	xcp, ycp, ip, jp, i1p, j1p, i2p, j2p;
	int     k, kp, xl, xr;
	float	s;
	PIXEL	*cp;

	/* Define the data array */
	nx = sp.nx;
	xcp = ANINT (lpcoords (ob->fitxc, 1, 2)) - p.x0;
	dx = nx / 2 - xcp;
	j1p = 0 - dx;
	j2p = nx - 1 - dx;
	j1p = (j1p < 0) ? 0 : j1p;
	j2p = (j2p >= p.dx) ? p.dx-1 : j2p;
	j1 = j1p + dx;
	j2 = j2p + dx;

	ny = sp.ny;
	ycp = ANINT (lpcoords (ob->fityc, 2, 2)) - p.y0;
	dy = ny / 2 - ycp;
	i1p = 0 - dy;
	i2p = ny - 1 - dy;
	i1p = (i1p < 0) ? 0 : i1p;
	i2p = (i2p >= p.dy) ? p.dy-1 : i2p;
	i1 = i1p + dy;
	i2 = i2p + dy;

	/* Set the sky level in NEQ.  Set a minimum sky at 10 * sky sigma. */
	/* For DARK objects flip the data in intensity. */

	sky = v0 + v1 * ob->sbr;
	if (sky < v1 * 10. * sp.skyhw)
	    sky = v1 * 10. * sp.skyhw;
	if (ob->eflgs & DARK)
	    s = -v1;
	else
	    s = v1;

	/* Set data array from image raster.  Allow for outside image raster */ 
	/* image by filling with zero.  Set points outside the isophote to */
	/* zero except if within 4 pixels of the center. */
	/* Convert to NEQ. */

	for (i = ny - 1; i > i2; i--)
	    for (j = 0; j < nx; j++)
		data[i][j] = 0.;
	for (ip = i2p; ip >= i1p; i--, ip--) {
	    for (j = nx - 1; j > j2; j--)
		data[i][j] = 0.;
	    jp = j2p;
	    cp = p.pic + ip * p.dx;
	    kp = a->d.ychst - p.y0 - ip;
	    if ((kp >= 0) && (kp < a->d.nlines)) {
		k = kp == 0 ? 0 : a->d.nx[kp-1];
	        for (; k < a->d.nx[kp]; k++) {
		    xl = a->d.xla[k] - p.x0;
		    xr = a->d.xra[k] - p.x0;
		    for (; jp > xr && jp >= j1p; j--, jp--)
			if ((abs (ip - ycp) < 4) && (abs (jp - xcp) < 4))
			    data[i][j] = s * (film (cp[jp]) - ob->sbr) + sky;
			else
			    data[i][j] = 0;
		    for (; jp >= xl && jp >= j1p; j--, jp--)
			data[i][j] = s * (film (cp[jp]) - ob->sbr) + sky;
		}
	    }
	    for (; jp >= j1p; j--, jp--)
		if ((abs (ip - ycp) < 4) && (abs (jp - xcp) < 4))
		    data[i][j] = s * (film (cp[jp]) - ob->sbr) + sky;
		else
		    data[i][j] = 0.;
	    for (; j >=0; j--)
		data[i][j] = 0.;
	}
	for (; i >= 0; i--)
	    for (j = 0; j < nx; j++)
		data[i][j] = 0.;

	/* Find the limits of the bounding box */

	rct[0] = nx; rct[1] = ny; rct[2] = 0; rct[3] = 0;
	for (i = 0; i < ny; i++)
	    for (j = 0; j < nx; j++)
	        if (data[i][j] > 0.) {
	            if (rct[0] > j) rct[0] = j;
	            if (rct[1] > i) rct[1] = i;
	            if (rct[3] < i) rct[3] = i;
		    break;
		}
	for (i = rct[1]; i <= rct[3]; i++)
	    for (j = nx-1; j >= 0; j--)
		if (data[i][j] > 0.) {
	    	    if (rct[2] < j) rct[2] = j;
		    break;
		}

	if (DEBUG2)
	    for (i = 0; i < ny; i++) {
		for (j = 0; j < nx; j++)
		    fprintf (stderr, "%4.0f ", data[i][j]);
		fprintf (stderr, "\n");
	    }
}

/* FITG -- Fit a gaussian to the PSF.  This is only for the purposes of
/* interpolation.  The PSF is assumed to be centered and the center position
/* is not determined.  The goal is to fit and remove the steep part of the
/* PSF so that we can use a combination of functional evaluation and linear
/* interpolation. */

fitg (A, B, C, D)
double *A, *B, *C, *D;
{
	int	i, j, xc, yc;
	double	x2, y2, xy, w, z, I0, model;
	double	a[4][5], b[3][4], c[2][3], d[1][2];

	for (i=0; i<4; i++)
	    for (j=0; j<5; j++)
		a[i][j] = 0.;

	xc = sp.nx / 2;
	yc = sp.ny / 2;
	I0 = template[yc][xc];
	for (i=0; i<sp.ny; i++) {
	    for (j=0; j<sp.nx; j++) {
		z = template[i][j];
		if (z > 0.1 * I0) {
		    w = z * z;
		    z = log (z);
		    x2 = (j - xc) * (j - xc);
		    y2 = (i - yc) * (i - yc);
		    xy = (j - xc) * (i - yc);

		    a[0][0] += w * z;
		    a[0][1] -= w;
		    a[0][2] += w * x2;
		    a[0][3] += w * y2;
		    a[0][4] += w * xy;

		    a[1][0] += w * z * x2;
		    a[1][1] -= w * x2;
		    a[1][2] += w * x2 * x2;
		    a[1][3] += w * y2 * x2;
		    a[1][4] += w * xy * x2;

		    a[2][0] += w * z * y2;
		    a[2][1] -= w * y2;
		    a[2][2] += w * x2 * y2;
		    a[2][3] += w * y2 * y2;
		    a[2][4] += w * xy * y2;

		    a[3][0] += w * z * xy;
		    a[3][1] -= w * xy;
		    a[3][2] += w * x2 * xy;
		    a[3][3] += w * y2 * xy;
		    a[3][4] += w * xy * xy;
		}
	    }
	}

	b[0][0] = a[0][0] - a[3][0] * a[0][4] / a[3][4];
	b[0][1] = a[0][1] - a[3][1] * a[0][4] / a[3][4];
	b[0][2] = a[0][2] - a[3][2] * a[0][4] / a[3][4];
	b[0][3] = a[0][3] - a[3][3] * a[0][4] / a[3][4];
	b[1][0] = a[1][0] - a[3][0] * a[1][4] / a[3][4];
	b[1][1] = a[1][1] - a[3][1] * a[1][4] / a[3][4];
	b[1][2] = a[1][2] - a[3][2] * a[1][4] / a[3][4];
	b[1][3] = a[1][3] - a[3][3] * a[1][4] / a[3][4];
	b[2][0] = a[2][0] - a[3][0] * a[2][4] / a[3][4];
	b[2][1] = a[2][1] - a[3][1] * a[2][4] / a[3][4];
	b[2][2] = a[2][2] - a[3][2] * a[2][4] / a[3][4];
	b[2][3] = a[2][3] - a[3][3] * a[2][4] / a[3][4];

	c[0][0] = b[0][0] - b[2][0] * b[0][3] / b[2][3];
	c[0][1] = b[0][1] - b[2][1] * b[0][3] / b[2][3];
	c[0][2] = b[0][2] - b[2][2] * b[0][3] / b[2][3];
	c[1][0] = b[1][0] - b[2][0] * b[1][3] / b[2][3];
	c[1][1] = b[1][1] - b[2][1] * b[1][3] / b[2][3];
	c[1][2] = b[1][2] - b[2][2] * b[1][3] / b[2][3];

	d[0][0] = c[0][0] - c[1][0] * c[0][2] / c[1][2];
	d[0][1] = c[0][1] - c[1][1] * c[0][2] / c[1][2];

	*A = -d[0][0] / d[0][1];
	*B = -(c[1][0] + *A * c[1][1]) / c[1][2];
	*C = -(b[2][0] + *A * b[2][1] + *B * b[2][2])/ b[2][3];
	*D = -(a[3][0] + *A * a[3][1] + *B * a[3][2] + *C * a[3][3])/ a[3][4];

	if (DEBUG1) {
	    for (i = 0; i < sp.ny; i++) {
		for (j = 0; j < sp.nx; j++)
		    fprintf (stderr, "%5.0f ", 1000. * template[i][j]/ I0);
		fprintf (stderr, "\n");
	    }
	    fprintf (stderr, "\n");
	    for (i = 0; i < sp.ny; i++) {
		for (j = 0; j < sp.nx; j++) {
		    x2 = (j - xc) * (j - xc);
		    y2 = (i - yc) * (i - yc);
		    xy = (j - xc) * (i - yc);
		    model = 1000. / I0 *
			exp (*A - (*B * x2 + *C * y2 + *D * xy));
		    fprintf (stderr, "%5.0f ", model);
		}
		fprintf (stderr, "\n");
	    }
	    fprintf (stderr, "\n");
	    for (i = 0; i < sp.ny; i++) {
		for (j = 0; j < sp.nx; j++) {
		    z = 1000. * template[i][j] / I0;
		    x2 = (j - xc) * (j - xc);
		    y2 = (i - yc) * (i - yc);
		    xy = (j - xc) * (i - yc);
		    model = 1000. / I0 *
			exp (*A - (*B * x2 + *C * y2 + *D * xy));
		    fprintf (stderr, "%5.0f ", z - model);
		}
		fprintf (stderr, "\n");
	    }
	}
	for (i = 0; i < sp.ny; i++) {
	    for (j = 0; j < sp.nx; j++) {
		x2 = (j - xc) * (j - xc);
		y2 = (i - yc) * (i - yc);
		xy = (j - xc) * (i - yc);
		model = exp (*A - (*B * x2 + *C * y2 + *D * xy));
		template[i][j] -= model;
	    }
	}
}