/*******************************************************************************
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 (Ti/Vi * (Di - Ti - Di * ln (Di / Ti)))

where the sum is over all pixels i in the image.
The Di are the observed data values and the Ti are the model template
values given by

	Ti = I0 * ti + fitsky

where ti is the normalized and sky subtracted template, I0 is a scaling
determined for each object and fitsky is the sky of each object.  The
Vi are the expected variance given by

	Vi = v0 + Ti

where v0 is the readout noise variance in data units.

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. fitsky <= 0.). 
*******************************************************************************/
#include	"focas1.h"

/* The following is to avoid name conflicts in the IRAF version */
#define	DEBUG1		NDEBUG1
#define	DEBUG2		NDEBUG2
#define	DEBUG3		NDEBUG3
#define	DEBUG4		NDEBUG4
#define	fit		newfit
#define	searchc		nsearchc
#define	search1		nsearch1
#define	search2		nsearch2
#define	fitmodel	nfitmodel
#define	getmodels	ngetmodels
#define	getimage	ngetimage
#define	getdata		ngetdata
#define	fitg		nfitg

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

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

int	tmpflag;

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

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

fit (ob, area, type, center, clsf, minrad, maxrad, rn, skysub)

struct	objrec *ob;
struct	areas *area;
int	type;
int	center;
int	clsf;			/* Check all templates? */
float	minrad, maxrad;		/* Minimum and maximum fitting radius */
float	rn, skysub;		/* Noise paramsters */

{
	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) || getimage (ob, center))
	    return (1);

	/* Initialize */
	v0 = rn * rn;
	sky0 = skysub;
	minr = minrad;
	maxr = maxrad;
	param = -1;
	for (i = 0; i < NTMPS; i++)
	    prob[i] = 1.;
	if (clsf)
	    ob->eflgs &= CLSF;
	else
	    ob->eflgs &= ~CLSF;

	if (center > 0)
	    searchc (ob, area, center, star);
	else {
	    ob->fitxc = (int) (ob->xc + 0.5);
	    ob->fityc = (int) (ob->yc + 0.5);
	    getdata (ob, area);
	}

	/* 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 (ob->eflgs & CLSF || 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)
struct	objrec	*ob;
struct	areas   *area;
int	center;
int	k;
{
	float	p, pmax;
	int	dx, dy, dxmax, dymax;
	float	fitmodel();

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

        ob->fitxc = (int) (ob->xc + 0.5);
	ob->fityc = (int) (ob->yc + 0.5);
	getdata (ob, area);
	prob[k] = 1.;
	p = fitmodel (ob, k);

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


/*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	I0, max;
	float	*dp, *mp;
	float	d, m, t, v, dt, mv, ln, f0, f1, f2, w0, w1, w2, p0, p1, p2;
	float	sum1, sum2, sum3, sum4, sum5, sum6, sum7, sum8, sum9;

	/* 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 - fitsky);
		    sum2 += m * m * d;
		}
		dp++;
		mp++;
	    }
	}
	I0 = fitsky * sum1 / sum2;

	/* Iteratively adjust the intensity scale to maximize the likelihood. */
	max = -1e30;
	for (iter = 0; iter < 20; iter++) {
	    if (v0 == 0.) {
		sum1 = sum2 = sum3 = 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;
			    t = I0 * m + fitsky;
			    dt = d / t;

			    ln = 1 / dt  - 1;
			    ln *= ln;
			    if (ln < .0004)
				f0 = -0.5 * d * ln;
			    else
				f0 = d - t - d * log (dt);
			    f1 = m * (dt - 1);
			    f2 = m * m * dt / t; 

			    sum1 += f0;
			    sum2 += f1;
			    sum3 += f2;
			}
			dp++;
			mp++;
		    }
		}
		p0 = sum1;
		p1 = sum2;
		p2 = sum3;
	    } else {
		sum1=sum2=sum3=sum4=sum5=sum6=sum7=sum8=sum9=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;
			    t = I0 * m + fitsky;
			    dt = d / t;
			    v = v0 + t;

			    ln = 1 / dt - 1;
			    ln *= ln;
			    if (ln < .0004)
				f0 = -0.5 * d * ln;
			    else
				f0 = d - t - d * log (dt);
			    f1 = m * (dt - 1);
			    f2 = m * m * dt / t;

			    mv = m / v;
			    w0 = t / v;
			    w1 = -mv * (w0 - 1);
			    w2 = -2 * mv * w1;

			    sum1 += w0 * f0;
			    sum2 += w0 * f1;
			    sum3 += w0 * f2;
			    sum4 += w1 * f0;
			    sum5 += w1 * f1;
			    sum6 += w2 * f0;
			    sum7 += w0;
			    sum8 += w1;
			    sum9 += w2;
			}
			dp++;
			mp++;
		    }
		}
		p0 = sum1 / sum7;
		p1 = (sum2 + sum4 - p0 * sum8) / sum7;
		p2 = (sum3 + 2*sum5 + sum6 - 2 * p1 * sum8 - p0 * sum9) / sum7;
	    }
	    if (p0 > max)
	        max = p0;
	    p1 = p1 / p2;

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

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

	/* 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)
int	type;
{
	int	i, j, k, xc, yc, xi, yi;
	static 	int tmpflg = 1;
	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)
	    return (tmpflg);
	tmpflag = 1;

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

	/* 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);
		    if (model[k][i][j] < 0.)
			model[k][i][j] = 0.;
		    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];
		    if (model[k][i][j] < 0.)
			model[k][i][j] = 0.;
		}
	}

	/* 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;

	tmpflg = 0;
	return (tmpflg);
}


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

getimage (ob, center)

struct	objrec  *ob;
int	center;

{
	rct[0] = (int) (lpcoords (ob->xc,1,2) + 0.5) - (sp.nx + 2 * center) / 2;
	rct[1] = (int) (lpcoords (ob->yc,2,2) + 0.5) - (sp.ny + 2 * center) / 2;
	rct[2] = (int) (lpcoords (ob->xc,1,2) + 0.5) + (sp.nx + 2 * center) / 2;
	rct[3] = (int) (lpcoords (ob->yc,2,2) + 0.5) + (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.  Allow for DARK objects.
/* Set the data array to zero outside the fitting region. */

getdata (ob, a)

struct	objrec  *ob;
struct	areas   *a;

{
	int	dx, dy, x0, y0, i, j, i1, j1, i2, j2;
	int	xcp, ycp, ip, jp, i1p, j1p, i2p, j2p;
	int     k, kp, xl, xr;
	float	s, r, fastr();
	PIXEL	*cp;

	/* Define the data array */
	dx = sp.nx;
	dy = sp.ny;
	x0 = (int) (lpcoords (ob->fitxc, 1, 2) + 0.5) - dx / 2;
	y0 = (int) (lpcoords (ob->fityc, 2, 2) + 0.5) - dy / 2;
	i1 = p.y0 - y0;
	j1 = p.x0 - x0;
	i2 = p.y0 + p.dy - 1 - y0;
	j2 = p.x0 + p.dx - 1 - x0;

	xcp = (int) (lpcoords (ob->fitxc, 1, 2) + 0.5) - p.x0;
	ycp = (int) (lpcoords (ob->fityc, 2, 2) + 0.5) - p.y0;
	i1p = y0 - p.y0;
	j1p = x0 - p.x0;
	i2p = y0 + dy - 1 - p.y0;
	j2p = x0 + dx - 1 - p.x0;

	/* Set a minimum sky at 10 * sky sigma. */

	fitsky = sky0 + ob->sbr;
	if (fitsky < 10 * sp.skyhw)
	    fitsky = 10 * sp.skyhw;

	/* 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 minr pixels of the center. */

	for (i = dy - 1; i > i2; i--)
	    for (j = 0; j < dx; j++)
		data[i][j] = 0.;
	if (ob->eflgs & DARK) {
	    s = fitsky + ob->sbr;
	    for (ip = i2p; i >= i1 && ip >= i1p; i--, ip--) {
		for (jp = j2p, j = dx - 1; j > j2; j--, jp--)
		    data[i][j] = 0.;
		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 >= 0; j--, jp--) {
			    r = fastr (jp-xcp, ip-ycp);
			    if (r <= minr)
				data[i][j] = s - film (cp[jp]);
			    else
				data[i][j] = 0;
			}
			for (; jp >= xl && jp >= j1p && j >= j1; j--, jp--) {
			    r = fastr (jp-xcp, ip-ycp);
			    if (r <= maxr)
				data[i][j] = s - film (cp[jp]);
			    else
				data[i][j] = 0;
			}
		    }
		}
		for (; jp >= j1p && j >= j1; j--, jp--) {
		    r = fastr (jp-xcp, ip-ycp);
		    if (r <= minr)
			data[i][j] = s - film (cp[jp]);
		    else
			data[i][j] = 0.;
		}
		for (; j >=0; j--)
		    data[i][j] = 0.;
	    }
	} else {
	    s = fitsky - ob->sbr;
	    for (ip = i2p; i >= i1 && ip >= i1p; i--, ip--) {
		for (jp = j2p, j = dx - 1; j > j2; j--, jp--)
		    data[i][j] = 0.;
		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 >= 0; j--, jp--) {
			    r = fastr (jp-xcp, ip-ycp);
			    if (r <= minr)
				data[i][j] = s + film (cp[jp]);
			    else
				data[i][j] = 0;
			}
			for (; jp >= xl && jp >= j1p && j >= j1; j--, jp--) {
			    r = fastr (jp-xcp, ip-ycp);
			    if (r <= maxr)
				data[i][j] = s + film (cp[jp]);
			    else
				data[i][j] = 0;
			}
		    }
		}
		for (; jp >= j1p && j >= j1; j--, jp--) {
		    r = fastr (jp-xcp, ip-ycp);
		    if (r <= minr)
			data[i][j] = s + film (cp[jp]);
		    else
			data[i][j] = 0.;
		}
		for (; j >=0; j--)
		    data[i][j] = 0.;
	    }
	}
	for (; i >= 0; i--)
	    for (j = 0; j < dx; j++)
		data[i][j] = 0.;

	/* Find the limits of the bounding box */

	rct[0] = dx; rct[1] = dy; rct[2] = 0; rct[3] = 0;
	for (i = 0; i < dy; i++)
	    for (j = 0; j < dx; 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 = dx-1; j >= 0; j--)
		if (data[i][j] > 0.) {
	    	    if (rct[2] < j) rct[2] = j;
		    break;
		}

	if (DEBUG2)
	    for (i = 0; i < dy; i++) {
		for (j = 0; j < dx; 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. */

/* 6/7/93: Negative template values are set to zero. */

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++) {
		if (template[i][j] < 0.)
		    template[i][j] = 0.;
		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;
	    }
	}
}