/*      FILE: /mxtools/src/mximstat/img_stat.c
 *   PURPOSE: 
 *    AUTHOR: Kenneth J. Mighell (mighell@noao.edu)
 *  LANGUAGE: ANSI C
 *      DATE: 2001NOV04
 * COPYRIGHT: (C) 2001 Assoc. of Universities for Research in Astronomy Inc.
 */

#include <mx.h>
#include "inc.h"

int inc_ImageS_Stat_i2(
  inc_ParS *ParS,
  inc_ImageS *ImageS
) {

  char mxfunc[] = "inc_ImageS_Stat_i2";
  int mxstatus = 0;
  int status = 1;

  int x;
  int nx;
  int y;
  int ny;
  int npx;
  double pixel;
  double max;
  int maxx;
  int maxy;
  double min;
  int minx;
  int miny;
  double sum;
  double sumsq;
  double delta;
  double adev;
  double mean;
  double sigma;
  double var;
  double rmean;
  int rmean_npx;
  double rmedian;
  double radev;
  double *tv;
  double rsigma;

  nx = ImageS->nxi;
  ny = ImageS->nyi;
  npx = ImageS->npxi;

  max = ImageS->matrixd[0][0];
  maxx = 0;
  maxy = 0;
  min = max;
  minx = maxx;
  miny = maxy;
  sum = 0;
  sumsq = 0;

  for (y=0; y<ny; ++y) { 
  for (x=0; x<nx; ++x) {

    pixel = ImageS->matrixd[y][x];

    sum += pixel;
    sumsq += pixel*pixel;

    if (pixel>max) {
      max = pixel;
      maxx = x;
      maxy = y;
    }

    if (pixel<min) {
      min = pixel;
      minx = x;
      miny = y;
    }

  }
  }
  mean = sum/npx;

  adev = 0;
  var = 0;
  for (y=0; y<ny; ++y) { 
  for (x=0; x<nx; ++x) {
    pixel = ImageS->matrixd[y][x];
    delta = pixel - mean;
    adev += fabs(delta);
    var += delta*delta;
  }
  }
  adev /= npx;
  var /= (npx-1);
  sigma = sqrt(var);
  
  tv = (double *)calloc( npx, sizeof(double));
  VecdRobustMdn_v6(
    ImageS->vectord,
    npx,
    tv,
    &rmean_npx,
    &rmean,
    &radev
  );
  rsigma = radev/0.78;

  ParS->omean = mean;
  ParS->osigma = sigma;
  ParS->oadev = adev;
  ParS->ormean = rmean;
  ParS->orsigma = rsigma;
  ParS->oradev = radev;
  ParS->ornpix = rmean_npx;
  ParS->omax = max;
  ParS->omaxx = maxx+1; /* C -> IRAF notation */
  ParS->omaxy = maxy+1; /* C -> IRAF notation */
  ParS->omin = min;
  ParS->ominx = minx+1; /* C -> IRAF notation */
  ParS->ominy = miny+1; /* C -> IRAF notation */
  ParS->onpix = npx;
  ParS->osum = sum;
  ParS->osumsq = sumsq;

ok:
  status = 0;
  goto bye;

mx_error1:
  mxp_tmpmsg_init_f0();
mx_error2:
  mxp_errmsg_append_f3 (mxfunc, mxstatus, MX.tmpmsg);
  goto bye;

bye:
  return(status);

}
/* end-of-file */