/*
 *  Template program for computing features 
 *
 *  Feature program fpstat computes pixel mean and standard deviation
 *   This version accepts all main pixel types
 *
 */
#include "VisXV4.h"    /* VisionX structure include file       */
#include "Vutil.h"     /* VisionX utility header files         */


VXparam_t par[] =      /* command line structure               */
{
{    "if=",    0, " input file,  vptemp: evaluate fpstat "},
{    "of=",    0, " output file "},
{     0,       0,    0}
};
#define  INIMAGE   par[0].val
#define  OUTCSV    par[1].val

int fpstat ( float *, VisXimage_t * );

int main(int argc, char** argv)
{
Vfstruct (im);
float result [3];
FILE *outfile;

VXparse(&argc, &argv, par);  /* parse the command line         */

  outfile = stdout;
  if (OUTCSV) { outfile = fopen (OUTCSV,"w"); }
  fprintf (outfile, "Count,Mean,Std.Dev\n");

  while ( Vfread( &im, INIMAGE) ) {
    switch ( im.type ) {
         case VX_PBYTE:
         case VX_PSHORT:
         case VX_PINT:
         case VX_PFLOAT:
         case VX_PDOUBLE:
	                  break;
   	 default: fprintf (stderr, "error: image type not supported\n");
	 exit (1);

    }
    fpstat (result, &im );
    fprintf (outfile, "%5.3f,%5.3f,%5.3f\n",result[0],result[1],result[2]);
  }

  if ( OUTCSV) { fclose (outfile); }
  exit(0);
}

/* feature subprogram  */

/* macro to accumulate stats */
#define CSTAT(type)                                 \
    for ( y= im->ylo; y <= im->yhi; y++) {           \
       for ( x= im->xlo; x <= im->xhi; x++) {        \
          sum += im->type[y][x];                     \
          sum2 += im->type[y][x] * im->type[y][x];   \
	  count++;                                   \
       }                                             \
    }

int fpstat ( float * result, VisXimage_t *im ) {
int x, y;
float sum, sum2, mean, sdev;
int count;
    count = 0;
    sum = 0.0;
    sum2 = 0.0;
    switch ( im->type ) {
        case VX_PBYTE:   CSTAT(u) break;
        case VX_PSHORT:  CSTAT(s) break;
        case VX_PINT:    CSTAT(i) break;
        case VX_PFLOAT:  CSTAT(f) break;
        case VX_PDOUBLE: CSTAT(d) break;
    }
    
    mean = sum/count;
    result[0] = count;  
    result[1] = mean;  
    result[2] = sqrt(sum2/count - mean * mean);  
    return (1);
}