/* @(#)cstat.c	17.1.1.1 (ES0-DMD) 01/25/02 17:50:05 */
/*===========================================================================
  Copyright (C) 1995 European Southern Observatory (ESO)
 
  This program is free software; you can redistribute it and/or 
  modify it under the terms of the GNU General Public License as 
  published by the Free Software Foundation; either version 2 of 
  the License, or (at your option) any later version.
 
  This program is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU General Public License for more details.
 
  You should have received a copy of the GNU General Public 
  License along with this program; if not, write to the Free 
  Software Foundation, Inc., 675 Massachusetss Ave, Cambridge, 
  MA 02139, USA.
 
  Corresponding concerning ESO-MIDAS should be addressed as follows:
	Internet e-mail: midas@eso.org
	Postal address: European Southern Observatory
			Data Management Division 
			Karl-Schwarzschild-Strasse 2
			D 85748 Garching bei Muenchen 
			GERMANY
===========================================================================*/

#include <math.h>
#define    MM              400       /* Max. no. of sample areas        */
#define    MPS             256       /* Max. no. of sample per area     */
#define    MNL              40       /* Max. no. of line read           */

int cstat(image,npix,sparam)
float *image;
float sparam[8];
int npix[2];
{
int ns,nps,nxd,nyd,nxf,nyf,nls,nx,ny,nsx,nsy;
int i,j, ix, iy, n, ioff, imm, mpw, npw, mft ;
float tlim,qlim, *pf;
double  fv, sv, tv, qv, vv, dv, fnps, fac;
double  fm[MM], sm[MM], tm[MM], qm[MM], sd;
ns = 9;
nps = 1600;
tlim = 0.25;
qlim = 0.5;
nx = npix[0]; ny = npix[1];
nxd = nyd = 100;
/*nxd = nx/ns; nyd = ny/ns;*/
nsx = nx/nxd;
nsy = ny/nxd;
nxf = nxd/2; nyf = nyd/2;
nls = nps / MNL; fnps = (double) (nls * MNL);
for (i=0; i<MM; i++) fm[i] = sm[i] = tm[i] = qm[i] = 0.0;
for (i=0; i<8; i++) sparam[i] = 0.0;

ioff =  (nyf-1)*nx;
imm = 0; mpw = 0;
for (iy=0; iy<nsy; iy++) {
     for (ix=nxf; ix<nx; ix += nxd) {
        fv = sv = tv = qv = 0.0;
        for (n=0; n<MNL; n++) {
           pf = &image[ioff + ix+n*nx];
           for (i=0; i<nls; i++) {
               vv = *pf++;
               fv +=vv;
               }
           }
        fv /=fnps; fm[imm] = fv;
        for (n=0; n<MNL; n++) {
           pf = &image[ioff+ix+n*nx];
           for (i=0; i<nls; i++) {
              vv = (*pf++ -fv); dv = vv*vv;
              sv += dv;
              tv += vv * dv;
              qv += dv * dv;
            }
         }
        fm[imm] = fv;
        sv /= fnps; sd = sv;
        sd = (sd<0.0) ? 0.0 : sqrt(sd); sm[imm] = sd;
        tv /= fnps;
        qv /=fnps;
        tv = tv/(sd*sd*sd); tm[imm] = tv;
        qv = qv/(sd*sd*sd*sd); qm[imm] = qv;
        imm++;
      }
     ioff += nx*nyd;
   }

/* compute mean values and save them in descriptor  */

fv = sv = tv = qv = 0.0; 
for (i=0; i<imm; i++) {
   fv += fm[i]; sv += sm[i]; tv += tm[i]; qv += qm[i];
}
fac = 1.0/(double) imm;
sparam[0] = fac*fv;  sparam[1] = fac*sv;
sparam[2] = fac*tv;  sparam[3] = fac*qv;
fv = sv = 0.0; tv = qv = fm[0]; n = 0;
for (i=0; i<imm; i++) {
     if (fabs(tm[i])<tlim && fabs(qm[i]-3.0)<qlim) {
        fv += fm[i]; sv += sm[i]; n++;
        if (qv<fm[i]) qv = fm[i];
        if (fm[i]<tv) tv = fm[i];
      }
}
fac = 1.0/(double) n;
sparam[4] = fac*fv;  sparam[5] = fac*sv;
sparam[6] = tv;  sparam[7] = qv;
/*printf("%d, %9.3f, %9.3f, %9.3f, %9.3f\n",n, sparam[0],sparam[1],sparam[2],sparam[3]); 
printf("%9.3f, %9.3f, %9.3f, %9.3f\n",sparam[4],sparam[5],sparam[6],sparam[7]); 
*/
}