#include	"focas1.h"

/* EVLSKY -- Main entry point for sky evaluation. */

evlsky (ob, area, pic, hstgrm, nhstgrm, step, start, skytype)
struct objrec  *ob;
struct areas   *area;
struct image   *pic;
int     hstgrm[];
int     nhstgrm;
PIXEL step, *start;
int     skytype;

{
    int     i, j, k, i1, i2, npix;
    PIXEL dmin, dmax, mean, mean1, mean2;

 /* Initialize */
    ob -> sbr = 0.;
    ob -> nsbr = 0;
    ob -> r1 = 0.;
    ob -> prob = 0.;
    ob -> ssbr = 0.;

    if (meanskyd (pic, &area -> d, &mean1) < 1)
	return (1);

    if (meanisphtd (pic, &area -> i, &mean2) < 1)
	return (1);;

    if (mean1 < mean2)
	mean = mean1;
    else
	mean = mean2;

    *start = mean - (nhstgrm / 2) * step;

 /* Determine the density histogram */

    if (skyhst (&area -> d, pic, hstgrm, nhstgrm, step, *start) < 1)
	return (1);

 /* Find the peak of the raw histogram */

    for (i = 0, k = 0; i < 512; i++)
	if (hstgrm[i] > k)
	    k = hstgrm[j = i];
    if (k == 0)
	return (1);
    ob -> r1 = film ((int) (*start + (j +.5) * step));

 /* Find the mean of the sky intensity for +- 2 sky sigma (density) * about the
    raw histogram peak */

    dmin = *start + (j +.5) * step - 2 * sp.skyhw;
    dmax = *start + (j +.5) * step + 2 * sp.skyhw;
    skystati (pic, &area -> d, dmin, dmax, ob);

 /* Find the peak of the histogram by convolution */

    ob -> prob = ob -> sbr;
    if (skytype == 1)
	evlsky1 (ob, hstgrm, nhstgrm, step, *start);

    return (0);
}

evlsky1 (ob, hstgrm, nhstgrm, step, start)
struct objrec  *ob;
int     hstgrm[];
int     nhstgrm;
PIXEL step, start;
{
    int     i, j;
    float   sig, dx, peak, convolution[512];
    static float    kernel[100];
    static int  nkernel = 100, kerflag = 1;

    if (kerflag) {
	sig = sp.skyhw / 2;
	for (i = 0; i < nkernel; i++) {
	    if (sig != 0)
		dx = i * step / sig;
	    else {
		printf("division by zero in evlsky.c\n");
		printf("Enter sky sigma for overide.\n");
		scanf ("%f",&sig);
		dx = i * step / sig;
	    }
	    if (dx > 3.) {
		nkernel = i;
		break;
	    }
	    kernel[i] = exp (-dx * dx / 2.);
	}
	kerflag = 0;
    }

    convolve (kernel, nkernel, hstgrm, convolution, nhstgrm);

    peak = 0.;
    for (i = nkernel; i < nhstgrm - nkernel; i++) {
	if (convolution[i] > peak) {
	    peak = convolution[i];
	    j = i;
	}
    }

    if ((j > 0) || (j < nhstgrm - 1))
	dx = 0.5 * (convolution[j + 1] - convolution[j - 1]) /
	    (2 * convolution[j] - convolution[j - 1] - convolution[j + 1]);

    ob -> sbr = film ((int) (start + (j + dx +.5) * step));
}

convolve (kernel, nkernel, datain, dataout, npts)
float   kernel[];
int     nkernel;
int     datain[];
float   dataout[];
int     npts;

{
    int     i, j;

    for (i = nkernel; i < npts - nkernel; i++) {
	dataout[i] = datain[i] * kernel[0];
	for (j = 1; j < nkernel; j++)
	    dataout[i] += (datain[i - j] + datain[i + j]) * kernel[j];
    }
}


meanskyd (pic, ar, mean)
struct image   *pic;
struct arblk   *ar;
PIXEL * mean;

{
    float   sum;
    int     nsum;
    int     i, j, k, l, xl, xr, xi, yi;
    PIXEL * cp;

    sum = 0.;
    nsum = 0;
    for (i = pic -> dy - 1, l = 0; i >= 0; i--) {
	yi = ar -> ychst - pic -> y0 - i;
	cp = pic -> pic + (i + 1) * pic -> dx - 1;
	if ((yi >= 0) && (yi < ar -> nlines)) {
	    for (xi = pic -> dx - 1; l < ar -> nx[yi]; l++) {
		xl = ar -> xla[l] - pic -> x0;
		xr = ar -> xra[l] - pic -> x0;
		for (; xi > xr; xi--) {
		    if (xi >= pic -> dx - sp.skwdth) {
			sum += *cp--;
			nsum++;
		    }
		    else
			cp--;
		}
		cp -= xi - xl + 1;
		xi = xl - 1;
	    }
	    for (; xi >= 0; xi--) {
		if (xi < sp.skwdth) {
		    sum += *cp--;
		    nsum++;
		}
		else
		    cp--;
	    }
	}
	else
	    if ((i < sp.skwdth) || (i >= pic -> dy - sp.skwdth)) {
		for (xi = pic -> dx - 1; xi >= 0; xi--) {
		    sum += *cp--;
		    nsum++;
		}
	    }
	    else {
		for (xi = pic -> dx - 1; xi >= 0; xi--) {
		    if ((xi < sp.skwdth) || (xi >= pic -> dx - sp.skwdth)) {
			sum += *cp--;
			nsum++;
		    }
		    else
			cp--;
		}
	    }
    }

    if (nsum > 0)
	*mean = sum / nsum;

    return (nsum);
}


/*Isophotal average*/

meanisphtd (pic, ar, mean)
struct image   *pic;
struct arblk   *ar;
PIXEL * mean;

{
    float   sum;
    int     nsum;
    int     i, j, k, l, xl, xr, xi, yi;
    PIXEL * cp;

    sum = 0.;
    nsum = 0;
    for (i = pic -> dy - 1, l = 0; i >= 0; i--) {
	yi = ar -> ychst - pic -> y0 - i;
	if ((yi >= 0) && (yi < ar -> nlines)) {
	    for (; l < ar -> nx[yi]; l++) {
		xl = ar -> xla[l] - pic -> x0;
		xr = ar -> xra[l] - pic -> x0;
		cp = pic -> pic + i * pic -> dx + xl;
		for (j = xl; j <= xr; j++) {
		    sum += *cp++;
		    nsum++;
		}
	    }
	}
    }

    if (nsum > 0)
	*mean = sum / nsum;

    return (nsum);
}


skystati (pic, ar, dmin, dmax, ob)
struct image   *pic;
struct arblk   *ar;
PIXEL dmin, dmax;
struct objrec  *ob;

{
    float   intensity, sumi, sumi2;
    int     nsum;
    int     i, j, k, l, xl, xr, xi, yi;
    PIXEL density, *cp;

    sumi = sumi2 = 0.;
    nsum = 0;
    for (i = pic -> dy - 1, l = 0; i >= 0; i--) {
	yi = ar -> ychst - pic -> y0 - i;
	cp = pic -> pic + (i + 1) * pic -> dx - 1;
	if ((yi >= 0) && (yi < ar -> nlines)) {
	    for (xi = pic -> dx - 1; l < ar -> nx[yi]; l++) {
		xl = ar -> xla[l] - pic -> x0;
		xr = ar -> xra[l] - pic -> x0;
		for (; xi > xr; xi--) {
		    density = *cp--;
		    if (xi >= pic -> dx - sp.skwdth) {
			if ((density >= dmin) && (density <= dmax)) {
			    intensity = film (density);
			    sumi += intensity;
			    sumi2 += intensity * intensity;
			    nsum++;
			}
		    }
		}
		cp -= xi - xl + 1;
		xi = xl - 1;
	    }
	    for (; xi >= 0; xi--) {
		density = *cp--;
		if ((xi < sp.skwdth) && (density >= dmin) && (density <= dmax)) {
		    intensity = film (density);
		    sumi += intensity;
		    sumi2 += intensity * intensity;
		    nsum++;
		}
	    }
	}
	else
	    if ((i < sp.skwdth) || (i >= pic -> dy - sp.skwdth)) {
		for (xi = pic -> dx - 1; xi >= 0; xi--) {
		    density = *cp--;
		    if ((density >= dmin) && (density <= dmax)) {
			intensity = film (density);
			sumi += intensity;
			sumi2 += intensity * intensity;
			nsum++;
		    }
		}
	    }
	    else {
		for (xi = pic -> dx - 1; xi >= 0; xi--) {
		    density = *cp--;
		    if ((xi < sp.skwdth) || (xi >= pic -> dx - sp.skwdth)) {
			if ((density >= dmin) && (density <= dmax)) {
			    intensity = film (density);
			    sumi += intensity;
			    sumi2 += intensity * intensity;
			    nsum++;
			}
		    }
		}
	    }
    }

    if (nsum > 0) {
	ob -> nsbr = nsum;
	ob -> sbr = sumi / nsum;
	ob -> ssbr = sumi2 * nsum - sumi * sumi;
	if (ob -> ssbr > 0.) {
	    ob -> ssbr = sumi2 * nsum - sumi * sumi;
	    ob -> ssbr = ob -> ssbr > 0. ? ob -> ssbr : 0.;
	    ob -> ssbr = sqrt (ob -> ssbr) / nsum;
	} else
	    ob -> ssbr = 0.;
    }
}