 /*
  * Khoros: $Id: lvstats.c,v 1.3 1992/03/20 23:07:47 dkhoros Exp $
  */

#if !defined(lint) && !defined(SABER)
static char rcsid[] = "Khoros: $Id: lvstats.c,v 1.3 1992/03/20 23:07:47 dkhoros Exp $";
#endif

 /*
  * $Log: lvstats.c,v $
 * Revision 1.3  1992/03/20  23:07:47  dkhoros
 * VirtualPatch5
 *
  */

/*
 *----------------------------------------------------------------------
 *
 * Copyright 1991, University of New Mexico.  All rights reserved.
 * Permission to copy and modify this software and its documen-
 * tation only for internal use in your organization is hereby
 * granted, provided that this notice is retained thereon and
 * on all copies.  UNM makes no representations as to the sui-
 * tability and operability of this software for any purpose.
 * It is provided "as is" without express or implied warranty.
 * 
 * UNM DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
 * INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FIT-
 * NESS.  IN NO EVENT SHALL UNM BE LIABLE FOR ANY SPECIAL,
 * INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY OTHER DAMAGES WHAT-
 * SOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
 * IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS
 * ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR PER-
 * FORMANCE OF THIS SOFTWARE.
 * 
 * No other rights, including, for example, the right to redis-
 * tribute this software and its documentation or the right to
 * prepare derivative works, are granted unless specifically
 * provided in a separate license agreement.
 *---------------------------------------------------------------------
 */

#include "unmcopyright.h"        /* Copyright 1991 by UNM */

/*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>  <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 >>>>
 >>>>         File Name: lvstats.c
 >>>>
 >>>>      Program Name: vstats
 >>>>
 >>>> Date Last Updated: Fri Nov  8 15:26:50 1991 
 >>>>
 >>>>          Routines: lvstats - the library call for vstats
 >>>>
 >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>   <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<*/


#include "vinclude.h"


/* -library_includes */
#include "vipl/lvstats.h"
/* -library_includes_end */


/****************************************************************
*
* Routine Name: lvstats - library call for vstats
*
* Purpose:
*    
*    lvstats computes the mean, variance, standard deviation, RMS
*    level,  peaks  with  locations,  integral,  skew,  kurtosis,
*    entropy, and contrast for an input xviff image.
*    
*    

* Input:
*    
*    img1           (struct xvimage) - input xvimage structure
*    
*    img2           (struct xvimage) - gating or masking image
*    
*    mflag          (int) - a flag set (equal 1) if gating image  mask
*                   available
*    
*    process_dir    (int) - 0 = vector processing, 1 = band, and  2  =
*                   calculate one set of stats for entire file.
*    
*    

* Output:
*    
*    stats          structure loaded with the values computed from the
*                   image.
*    
*    

*
* Written By: Donna Koechner, John Rasure
*    
*    Audit Done  05-OCT-90  Donna Koechner
*    
*    06-OCT-90:    Fixed   fmin,fmax   initialization   bug    in
*    calculate_byte().  Donna Koechner
*    
*    Tue Mar 5 1991 - changed  method  of  calculating  variance,
*    stddev, skew and kurtosis.  Donna Koechner
*    
*    

****************************************************************/


/* -library_def */
int
lvstats(img1,img2,mflag,stats,process_dir)
struct xvimage  *img1,*img2;
int             mflag, process_dir;
struct lvstat   **stats;
/* -library_def_end */

/* -library_code */
{

  char *program = "lvstats";
  int calculate_byte();
  int calculate_short();
  int calculate_long();
  int calculate_float();
  int calculate_double();

  switch (img1->data_storage_type)
  {
    case VFF_TYP_BIT :
         (void)fprintf(stderr,"%s: Statistics unavailable for BIT images\n",
                       program);
          return(0);
          break;
    case VFF_TYP_1_BYTE : 
         if (! calculate_byte(img1, img2, stats, mflag, process_dir))
         {
           (void) fprintf(stderr,"%s: calculate_byte failed\n",program);
           return(0);
         }
         break;
    case VFF_TYP_2_BYTE :
         if (! calculate_short(img1, img2, stats, mflag, process_dir))
         {
           (void) fprintf(stderr,"%s: calculate_short failed\n",program);
           return(0);
         }
         break;
    case VFF_TYP_4_BYTE :
         if (! calculate_long(img1, img2, stats, mflag, process_dir))
         {
             (void) fprintf(stderr,"%s: calculate_long failed\n",program);
             return(0);
         }
         break;
    case VFF_TYP_FLOAT :
         if (! calculate_float(img1, img2, stats, mflag, process_dir))
         {
           (void) fprintf(stderr,"%s: calculate_float failed\n",program);
           return(0);
         }
         break;
    case VFF_TYP_DOUBLE :
         if (! calculate_double(img1, img2, stats, mflag, process_dir))
         {
           (void) fprintf(stderr,"%s: calculate_double failed\n",program);
           return(0);
         }
         break;
    case VFF_TYP_COMPLEX :
         (void) fprintf(stderr,
               "lvstats: Statistics unavailable on COMPLEX images\n");
         return(0);
         break;
    default:
         (void) fprintf(stderr,"lvstats: Unknown data storage type, %d\n",
                        img1->data_storage_type);
         return(0);
         break;
  }

  return(1); /* success */
}


/**************************************************************
*
* MODULE NAME: calculate_byte
*
*     PURPOSE: calculates the statistics for a BYTE image. The
*              stats structure stats will be filled out.
*
*       INPUT:  1. img1 -- first xvimage strcuture
*               2. img2 -- gating mask image
*               3. mflag -- a flag set (equal 1) if gating image mask
*                          available
*               4. process_dir -- 0 = vector processing, 1 = band,
*                                 and 2 = calculate one set of stats
*                                 for entire file.
*
*               The input images must be of the same data type
*               and size.
*
*      OUTPUT:  
*               The stat structure is loaded with the values computed
*               from the image.
*
* CALLED FROM: lvstats
*
* ROUTINES CALLED: dload_vector()
*
**************************************************************/

int calculate_byte(img1, img2, stats, mflag, process_dir)

int             process_dir, mflag;
struct xvimage  *img1, *img2;
struct lvstat   **stats;

{

  int           pmax,pmin,      /* positions of min and max values in vector */
                pos,neg,        /* number of pos and neg points in vector */
                dim1, dim2,     /* vector dimension of input/masking files */
                num_vects1,     /* number of vectors in input file */
                num_vects2,     /* number of vectors in masking file */
                tmp_band,
                tmp_col_size,
                i, j, n, k;

  unsigned int  hist[256];      /* holds histogram */

  unsigned char **c1,**c2;      /* image and masking vector arrays */

  double        skewness,       /* sum of cubed values in vector */
                kurtosis,       /* sum of values to the 4th pow in vector */
                variance,
                diff,
                stddev,
                sum,            /* sum of values in vector */
                rms,            /* sum of squared values in vector */
                pint,nint,      /* pos & neg integral sums */
                fmin,fmax,      /* min and max values in vector */
                norm_hist,
                entropy, contrast,
                dtmp,
                two = 2.0;

  struct lvstat *new_stats;

  char          **dload_vector(),
                *program = "calculate_byte";



  if ( (process_dir) > 2 || (process_dir < 0) )
  {
    (void)fprintf(stderr,"%s: unknown processing direction (%d)\n",
                  program, process_dir);
    return(0);
  }

  if (process_dir == 2)  /* calculate one set of stats for entire file */
  {
    tmp_band = img1->num_data_bands;
    tmp_col_size = img1->col_size;
    img1->col_size *= tmp_band;
    img1->num_data_bands = 1;
    if ((c1 = (unsigned char **)dload_vector(img1, &num_vects1, &dim1, 
              DSP_BAND)) == NULL)
    {
      (void)fprintf(stderr,"%s: dload_vector failed\n",program);
      return(0);
    }
    free(img1->imagedata);

    if (mflag)
    {
      img2->col_size *= img2->num_data_bands;
      img2->num_data_bands = 1;
      if ((c2 = (unsigned char**)dload_vector(img2, &num_vects2, &dim2,
                DSP_BAND)) == NULL)
      {
        (void)fprintf(stderr,"%s: dload_vector failed\n",program);
        return(0);
      }
      if ( (num_vects1 != num_vects2) || (dim1 != dim2) )
      {
        (void)fprintf(stderr,"%s: incompatible mask size\n",program);
        return(0);
      }
      img2->num_data_bands = tmp_band;
      img2->col_size = tmp_col_size;
      free(img2->imagedata);
    }
    img1->num_data_bands = tmp_band;
    img1->col_size = tmp_col_size;
  }

  else  /* calculate stats on each band or vector separately */
  {
    if ((c1 = (unsigned char **)dload_vector(img1, &num_vects1, &dim1, 
              process_dir)) == NULL)
    {
      (void)fprintf(stderr,"%s: dload_vector failed\n",program);
      return(0);
    }
    free(img1->imagedata);
    if (mflag)
    {
      if ((c2 = (unsigned char**)dload_vector(img2, &num_vects2, &dim2,
                process_dir)) == NULL)
      {
        (void)fprintf(stderr,"%s: dload_vector failed\n",program);
        return(0);
      }
      if ( (num_vects1 != num_vects2) || (dim1 != dim2) )
      {
        (void)fprintf(stderr,"%s: incompatible mask size\n",program);
        return(0);
      }
      free(img2->imagedata);
    }
  }

  if (dim1 < 2) {
    (void)fprintf(stderr,"calculate_byte: will not calculate statistics ");
    (void)fprintf(stderr,"on less than 2 points\n");
    (void)fprintf(stderr,"Data has only %d points per vector.\n",dim1);
    return(0);
  }

  if (!mflag)
  {
    for (i=0; i<num_vects1; i++)
    {
      /* Zero out variables */
      /* pos = 0; */
      neg = 0;
      sum = 0.0;
      nint = 0.0;
      /* pint = 0.0; */
      rms = 0.0;
      variance = 0.0;
      entropy = 0.0;
      contrast = 0.0;
      skewness = 0.0;
      kurtosis = 0.0;
      fmin = (double) c1[i][0];
      fmax = (double) c1[i][0];
      pmin = 0;
      pmax = 0;
      bzero((char *) hist, 256 * sizeof(unsigned int));
      k = 0;
      /* calculate mean */
      for (j=0; j<dim1; j++) {
        sum += (double) c1[i][j];
        k++;
      }
      pint = sum;
      pos = k;
      new_stats = stats[i];
      new_stats->integ = (float)sum;
      sum /= (float)k;
      new_stats->mean = (float)sum;

      for (j=0; j<dim1; j++)
      {
        dtmp = (double) c1[i][j];
        hist[c1[i][j]]++;
        if (dtmp > fmax)
        {
          fmax = dtmp;
          pmax = j;
        }
        if (dtmp < fmin)
        {
          fmin = dtmp;
          pmin = j;
        }
        /* calculate variance, skew, and kurtosis */
        rms += dtmp * dtmp;
        diff = dtmp - sum;
        dtmp = diff * diff;
        variance += dtmp;
        dtmp *= diff;
        skewness += dtmp;
        dtmp *= diff;
        kurtosis += dtmp;
      }

      /* calculate entropy */
      for ( n = 0 ; n < 256; n++ )
      {
        if (hist[n] != 0) 
        {
          norm_hist = (double) hist[n]/k;
          entropy += norm_hist * (log10(norm_hist) / log10(two));
          contrast += (double) n * (double) n * norm_hist;
        }
      }
      variance /= (double)(k-1);
      stddev = sqrt(variance);
      if(variance != 0.0) {
        skewness /= ((double)k * pow(stddev,3.0));
        kurtosis = (kurtosis/((double)k * pow(variance,2.0))) - 3.0;
      }

      new_stats->rms = (float) sqrt(rms/(float)k);
      new_stats->pinteg = (float) pint;
      new_stats->ninteg = (float) nint;
      new_stats->pos = pos;
      new_stats->neg = neg;
      new_stats->var = (float) variance;
      new_stats->sdev = (float) stddev;
      new_stats->cpts = k;
      new_stats->fmax = (float) fmax;
      new_stats->maxy = (int)(pmax/img1->row_size);
      new_stats->maxx = pmax - (new_stats->maxy * img1->row_size);
      new_stats->maxp = pmax; 
      new_stats->fmin = (float) fmin;
      new_stats->miny = (int)(pmin/img1->row_size);
      new_stats->minx = pmin - (new_stats->miny * img1->row_size);
      new_stats->minp = pmin; 
      new_stats->entropy = (float) -entropy;
      new_stats->contrast = (float) contrast;
      new_stats->skewness = (float) skewness;
      new_stats->kurtosis = (float) kurtosis;
    }
  }

  else
  {
    for (i=0; i<num_vects1; i++)
    {
      /* Zero out variables */
      /* pos = 0; */
      neg = 0;
      sum = 0.0;
      nint = 0.0;
      /* pint = 0.0; */
      rms = 0.0;
      variance = 0.0;
      entropy = 0.0;
      contrast = 0.0;
      skewness = 0.0;
      kurtosis = 0.0;
      for(j=0; j<dim1; j++) {
        if(c2[i][j]) {
          fmin = (double) c1[i][j];
          fmax = (double) c1[i][j];
          pmin = j;
          pmax = j;
          j=dim1;
        }
      }
      bzero((char *) hist, 256 * sizeof(unsigned int));
      k = 0;
      /* calculate mean */
      for (j=0; j<dim1; j++) {
        if (c2[i][j]) {
          sum += (double) c1[i][j];
          k++;
        }
      }
      if (k < 2) {
        (void)fprintf(stderr,"calculate_byte: will not calculate statistics ");
        (void)fprintf(stderr,"on less than 2 points\n");
        (void)fprintf(stderr,"Data has only %d points per vector.\n",k);
        return(0);
      }
      pint = sum;
      pos = k;
      new_stats = stats[i];
      new_stats->integ = (float)sum;
      sum /= (float)k;
      new_stats->mean = (float)sum;

      for (j=0; j<dim1; j++)
      {
        if (c2[i][j])
        {
          dtmp = (double) c1[i][j];
          hist[c1[i][j]]++;
          if (dtmp > fmax)
          {
            fmax = dtmp;
            pmax = j;
          }
          if (dtmp < fmin)
          {
            fmin = dtmp;
            pmin = j;
          }
          /* calculate variance, skew, and kurtosis */
          rms += dtmp * dtmp;
          diff = dtmp - sum;
          dtmp = diff * diff;
          variance += dtmp;
          dtmp *= diff;
          skewness += dtmp;
          dtmp *= diff;
          kurtosis += dtmp;
        }
      }
      variance /= (double)(k-1);
      stddev = sqrt(variance);
      if(variance != 0.0) {
        skewness /= ((double)k * pow(stddev,3.0));
        kurtosis = (kurtosis/((double)k * pow(variance,2.0))) - 3.0;
      }

      /* calculate entropy, contrast */
      for ( n = 0 ; n < 256; n++ )
      {
        if (hist[n] != 0) 
        {
          norm_hist = (double) hist[n]/k;
          entropy += norm_hist * (log10(norm_hist) / log10(two));
          contrast += (double) n * (double) n * norm_hist;
        }
      }

      new_stats->rms = (float) sqrt(rms/(float)k);
      new_stats->pinteg = (float) pint;
      new_stats->ninteg = (float) nint;
      new_stats->pos = pos;
      new_stats->neg = neg;
      new_stats->var = (float) variance;
      new_stats->sdev = (float) stddev;
      new_stats->cpts = k;
      new_stats->fmax = (float) fmax;
      new_stats->maxy = (int)(pmax/img1->row_size);
      new_stats->maxx = pmax - (new_stats->maxy * img1->row_size);
      new_stats->maxp = pmax;   
      new_stats->fmin = (float) fmin;
      new_stats->miny = (int)(pmin/img1->row_size);
      new_stats->minx = pmin - (new_stats->miny * img1->row_size);
      new_stats->minp = pmin;   
      new_stats->entropy = (float) -entropy;
      new_stats->contrast = (float) contrast;
      new_stats->skewness = (float) skewness;
      new_stats->kurtosis = (float) kurtosis;
    }
  }

  for (i=0; i<num_vects1; i++)
  {
    free(c1[i]);
    if(mflag) free(c2[i]);
  }
  free(c1);
  if(mflag) free(c2);

  return(1);
}

/**************************************************************
*
* MODULE NAME: calculate_short
*
*     PURPOSE: calculates the statistics for a SHORT image. The
*              stats structure stats will be filled out.
*
*       INPUT:  1. img1 -- first xvimage strcuture
*               2. img2 -- gating mask image
*               3. mflag -- a flag set (equal 1) if gating image mask
*                          available
*               4. process_dir -- 0 = vector processing, 1 = band,
*                                 and 2 = calculate one set of stats 
*                                 for entire file.
*
*               The input images must be of the same data type
*               and size.
*
*      OUTPUT:  
*               The stat structure is loaded with the values computed
*               from the image.
*
* CALLED FROM: lvstats
*
* ROUTINES CALLED: dload_vector()
**************************************************************/

int calculate_short(img1, img2, stats, mflag, process_dir)

struct xvimage  *img1,*img2;
int             mflag, process_dir;
struct lvstat   **stats;

{
  int           pmax,pmin,      /* positions of min and max values in vector */
                pos,neg,        /* number of pos and neg points in vector */
                dim1, dim2,     /* vector dimension of input/masking files */
                num_vects1,     /* number of vectors in input file */
                num_vects2,     /* number of vectors in masking file */
                tmp_band,
                tmp_col_size,
                i, j, k;

  short         **s1, **s2;     /* image and masking vector arrays */

  double        skewness,       /* sum of cubed values in vector */
                kurtosis,       /* sum of values to the 4th pow in vector */
                variance,
                diff,
                stddev,
                sum,            /* sum of values in vector */
                rms,            /* sum of squared values in vector */
                pint,nint,      /* pos & neg integral sums */
                fmin,fmax,      /* min and max values in vector */
                entropy, contrast,
                dtmp;

  struct lvstat *new_stats;

  char          **dload_vector(),
                *program = "calculate_short";


  if ( (process_dir) > 2 || (process_dir < 0) )
  {
    (void)fprintf(stderr,"%s: unknown processing direction (%d)\n",
                  program, process_dir);
    return(0);
  }

  if (process_dir == 2)  /* calculate one set of stats for entire file */
  {
    tmp_band = img1->num_data_bands;
    tmp_col_size = img1->col_size;
    img1->col_size *= tmp_band;
    img1->num_data_bands = 1;
    if ((s1 = (short **)dload_vector(img1, &num_vects1, &dim1, DSP_BAND))
        == NULL)
    {
      (void)fprintf(stderr,"%s: dload_vector failed \n",program);
      return(0);
    }
    free(img1->imagedata);

    if (mflag)
    {
      img2->col_size *= img2->num_data_bands;
      img2->num_data_bands = 1;
      if ((s2 = (short **)dload_vector(img2, &num_vects2, &dim2, DSP_BAND))
          == NULL)
      {
        (void)fprintf(stderr,"%s: dload_vector failed\n",program);
        return(0);
      }
      if ( (num_vects1 != num_vects2) || (dim1 != dim2) )
      {
        (void)fprintf(stderr,"%s: incompatible mask size\n",program);
        return(0);
      }
      img2->num_data_bands = tmp_band;
      img2->col_size = tmp_col_size;
      free(img2->imagedata);
    }
    img1->num_data_bands = tmp_band;
    img1->col_size = tmp_col_size;
  }

  else  /* calculate stats on each band or vector separately */
  {
    if ((s1 = (short **)dload_vector(img1, &num_vects1, &dim1, process_dir))
        == NULL)
    {
      (void)fprintf(stderr,"%s: dload_vector failed\n",program);
      return(0);
    }
    free(img1->imagedata);
    if (mflag)
    {
      if ((s2 = (short **)dload_vector(img2, &num_vects2, &dim2, process_dir))
          == NULL)
      {
        (void)fprintf(stderr,"%s: dload_vector failed\n",program);
        return(0);
      }
      if ( (num_vects1 != num_vects2) || (dim1 != dim2) )
      {
        (void)fprintf(stderr,"%s: incompatible mask size\n",program);
        return(0);
      }
      free(img2->imagedata);
    }
  }

  if (dim1 < 2) {
    (void)fprintf(stderr,"calculate_short: will not calculate statistics ");
    (void)fprintf(stderr,"on less than 2 points\n");
    (void)fprintf(stderr,"Data has only %d points per vector.\n",dim1);
    return(0);
  }

  if (!mflag)
  {
    for (i = 0; i < num_vects1; i++)
    {
      /* Zero out variables */
      pos = 0;
      neg = 0;
      sum = 0.0;
      nint = 0.0;
      pint = 0.0;
      rms = 0.0;
      variance = 0.0;
      entropy = 0.0;
      contrast = 0.0;
      skewness = 0.0;
      kurtosis = 0.0;
      fmin = (double) s1[i][0];
      fmax = (double) s1[i][0];
      pmin = 0;
      pmax = 0;
      k = 0;
      /* calculate mean */
      for (j=0; j<dim1; j++) {
        sum += (double) s1[i][j];
        k++;
      }
      new_stats = stats[i];
      new_stats->integ = (float)sum;
      sum /= (float)k;
      new_stats->mean = (float)sum;

      for (j=0; j<dim1; j++)
      {
        dtmp = (double) s1[i][j];
        /* calculate pos & neg intregrals & num points, & min/max */
        if (dtmp >= 0) { pint += dtmp; pos++; }
        if (dtmp <  0) { nint += dtmp; neg++; }
        if (dtmp > fmax)
        {
          fmax = dtmp;
          pmax = j;
        }
        if (dtmp < fmin)
        {
          fmin = dtmp;
          pmin = j;
        }
        /* calculate variance, skew, and kurtosis */
        rms += dtmp * dtmp;
        diff = dtmp - sum;
        dtmp = diff * diff;
        variance += dtmp;
        dtmp *= diff;
        skewness += dtmp;
        dtmp *= diff;
        kurtosis += dtmp;
      }
      variance /= (double)(k-1);
      stddev = sqrt(variance);
      if(variance != 0.0) {
        skewness /= ((double)k * pow(stddev,3.0));
        kurtosis = (kurtosis/((double)k * pow(variance,2.0))) - 3.0;
      }

      new_stats->rms = (float) sqrt(rms/(float)k);
      new_stats->pinteg = (float)pint;
      new_stats->ninteg = (float)nint;
      new_stats->pos = pos;
      new_stats->neg = neg;
      new_stats->var = (float) variance;
      new_stats->sdev = (float) stddev;
      new_stats->cpts = k;
      new_stats->fmax = (float) fmax;
      new_stats->maxy = (int)(pmax/img1->row_size);
      new_stats->maxx = pmax - (new_stats->maxy * img1->row_size);
      new_stats->maxp = pmax;   
      new_stats->fmin = (float) fmin;
      new_stats->miny = (int)(pmin/img1->row_size);
      new_stats->minx = pmin - (new_stats->miny * img1->row_size);
      new_stats->minp = pmin;   
      new_stats->skewness = (float) skewness;
      new_stats->kurtosis = (float) kurtosis;
      new_stats->entropy = entropy;
      new_stats->contrast = contrast;
    }
  }

  else  /* use masking image */
  {
    for (i = 0; i < num_vects1; i++)
    {
      /* Zero out variables */
      pos = 0;
      neg = 0;
      sum = 0.0;
      nint = 0.0;
      pint = 0.0;
      rms = 0.0;
      variance = 0.0;
      entropy = 0.0;
      contrast = 0.0;
      skewness = 0.0;
      kurtosis = 0.0;
      for(j=0; j<dim1; j++) {
        if(s2[i][j]) {
          fmin = (double) s1[i][j];
          fmax = (double) s1[i][j];
          pmin = j;
          pmax = j;
          j=dim1;
        }
      }
      k = 0;
      /* calculate mean */
      for (j=0; j<dim1; j++) {
        if (s2[i][j]) {
          sum += (double) s1[i][j];
          k++;
        }
      }
      if (k < 2) {
        (void)fprintf(stderr,"calculate_short: will not calculate statistics ");
        (void)fprintf(stderr,"on less than 2 points\n");
        (void)fprintf(stderr,"Data has only %d points per vector.\n",k);
        return(0);
      }
      new_stats = stats[i];
      new_stats->integ = (float)sum;
      sum /= (float)k;
      new_stats->mean = (float)sum;

      for (j=0; j<dim1; j++)
      {
        if (s2[i][j])
        {
          dtmp = (double) s1[i][j];
          /* calculate pos, neg integrals, num points, & min/max */
          if (dtmp >= 0) { pint += dtmp; pos++; }
          if (dtmp <  0) { nint += dtmp; neg++; }
          if (dtmp > fmax)
          {
            fmax = dtmp;
            pmax = j;
          }
          if (dtmp < fmin)
          {
            fmin = dtmp;
            pmin = j;
          }
          /* calculate variance, skew, and kurtosis */
          rms += dtmp * dtmp;
          diff = dtmp - sum;
          dtmp = diff * diff;
          variance += dtmp;
          dtmp *= diff;
          skewness += dtmp;
          dtmp *= diff;
          kurtosis += dtmp;
        }
      }
      variance /= (double)(k-1);
      stddev = sqrt(variance);
      if(variance != 0.0) {
        skewness /= ((double)k * pow(stddev,3.0));
        kurtosis = (kurtosis/((double)k * pow(variance,2.0))) - 3.0;
      }

      new_stats->rms = (float) sqrt(rms/(float)k);
      new_stats->pinteg = (float) pint;
      new_stats->ninteg = (float) nint;
      new_stats->pos = pos;
      new_stats->neg = neg;
      new_stats->var = (float) variance;
      new_stats->sdev = (float) stddev;
      new_stats->cpts = k;
      new_stats->fmax = (float) fmax;
      new_stats->maxy = (int)(pmax/img1->row_size);
      new_stats->maxx = pmax - (new_stats->maxy * img1->row_size);
      new_stats->maxp = pmax;   
      new_stats->fmin = (float) fmin;
      new_stats->miny = (int)(pmin/img1->row_size);
      new_stats->minx = pmin - (new_stats->miny * img1->row_size);
      new_stats->minp = pmin;   
      new_stats->skewness = (float) skewness;
      new_stats->kurtosis = (float) kurtosis;
      new_stats->entropy = entropy;
      new_stats->contrast = contrast;
    }
  }

  for (i=0; i<num_vects1; i++)
  {
    free(s1[i]);
    if(mflag) free(s2[i]);
  }
  free(s1);
  if(mflag) free(s2);

  return(1);
}

/**************************************************************
*
* MODULE NAME: calculate_long
*
*     PURPOSE: calculates the statistics for a LONG image. The
*              stats structure stats will be filled out.
*
*       INPUT:  1. img1 -- first xvimage strcuture
*               2. img2 -- gating mask image
*               3. mflag -- a flag set (equal 1) if gating image mask
*                          available
*               4. process_dir -- 0 = vector processing, 1 = band,
*                                 and 2 = calculate one set of stats
*                                 for entire file.
*
*               The input images must be of the same data type
*               and size.
*
*      OUTPUT:  
*               The stat structure is loaded with the values computed
*               from the image.
*
* CALLED FROM: lvstats
*
* ROUTINES CALLED: dload_vector()
**************************************************************/

int calculate_long(img1, img2, stats, mflag, process_dir)

struct xvimage  *img1,*img2;
int             mflag, process_dir;
struct lvstat   **stats;

{
  int           pmax,pmin,      /* positions of min and max values in vector */
                pos,neg,        /* number of pos and neg points in vector */
                dim1, dim2,     /* vector dimension of input/masking files */
                num_vects1,     /* number of vectors in input file */
                num_vects2,     /* number of vectors in masking file */
                tmp_band,
                tmp_col_size,
                i, j, k;

  int           **i1, **i2;     /* image and masking vector arrays */

  double        skewness,       /* sum of cubed values in vector */
                kurtosis,       /* sum of values to the 4th pow in vector */
                variance,
                diff,
                stddev,
                sum,            /* sum of values in vector */
                rms,            /* sum of squared values in vector */
                pint,nint,      /* pos & neg integral sums */
                fmin,fmax,      /* min and max values in vector */
                entropy, contrast,
                dtmp;

  struct lvstat *new_stats;

  char          **dload_vector(),
                *program = "calculate_long";


  if ( (process_dir) > 2 || (process_dir < 0) )
  {
    (void)fprintf(stderr,"%s: unknown processing direction (%d)\n",
                  program, process_dir);
    return(0);
  }

  if (process_dir == 2)  /* calculate one set of stats for entire file */
  {
    tmp_band = img1->num_data_bands;
    tmp_col_size = img1->col_size;
    img1->col_size *= tmp_band;
    img1->num_data_bands = 1;
    if ((i1 = (int **)dload_vector(img1, &num_vects1, &dim1, DSP_BAND))
        == NULL)
    {
      (void)fprintf(stderr,"%s: dload_vector failed \n",program);
      return(0);
    }
    free(img1->imagedata);

    if (mflag)
    {
      img2->col_size *= img2->num_data_bands;
      img2->num_data_bands = 1;
      if ((i2 = (int **)dload_vector(img2, &num_vects2, &dim2, DSP_BAND))
          == NULL)
      {
        (void)fprintf(stderr,"%s: dload_vector failed\n",program);
        return(0);
      }
      if ( (num_vects1 != num_vects2) || (dim1 != dim2) )
      {
        (void)fprintf(stderr,"%s: incompatible mask size\n",program);
        return(0);
      }
      img2->num_data_bands = tmp_band;
      img2->col_size = tmp_col_size;
      free(img2->imagedata);
    }
    img1->num_data_bands = tmp_band;
    img1->col_size = tmp_col_size;
  }

  else  /* calculate stats on each band or vector separately */
  {
    if ((i1 = (int **)dload_vector(img1, &num_vects1, &dim1, process_dir))
        == NULL)
    {
      (void)fprintf(stderr,"%s: dload_vector failed\n",program);
      return(0);
    }
    free(img1->imagedata);
    if (mflag)
    {
      if ((i2 = (int **)dload_vector(img2, &num_vects2, &dim2, process_dir))
          == NULL)
      {
        (void)fprintf(stderr,"%s: dload_vector failed\n",program);
        return(0);
      }
      if ( (num_vects1 != num_vects2) || (dim1 != dim2) )
      {
        (void)fprintf(stderr,"%s: incompatible mask size\n",program);
        return(0);
      }
      free(img2->imagedata);
    }
  }

  if (dim1 < 2) {
    (void)fprintf(stderr,"calculate_long: will not calculate statistics ");
    (void)fprintf(stderr,"on less than 2 points\n");
    (void)fprintf(stderr,"Data has only %d points per vector.\n",dim1);
    return(0);
  }

  if (!mflag)
  {
    for (i = 0; i < num_vects1; i++)
    {
      /* Zero out variables */
      pos = 0;
      neg = 0;
      sum = 0.0;
      nint = 0.0;
      pint = 0.0;
      rms = 0.0;
      variance = 0.0;
      entropy = 0.0;
      contrast = 0.0;
      skewness = 0.0;
      kurtosis = 0.0;
      fmin = (double) i1[i][0];
      fmax = (double) i1[i][0];
      pmin = 0;
      pmax = 0;
      k = 0;
      /* calculate mean */
      for (j=0; j<dim1; j++) {
        sum += (double) i1[i][j];
        k++;
      }
      new_stats = stats[i];
      new_stats->integ = (float)sum;
      sum /= (float)k;
      new_stats->mean = (float)sum;

      for (j=0; j<dim1; j++)
      {
        dtmp = (double) i1[i][j];
        /* calculate pos & neg intregrals & num points, & min/max */
        if (dtmp >= 0) { pint += dtmp; pos++; }
        if (dtmp <  0) { nint += dtmp; neg++; }
        if (dtmp > fmax)
        {
          fmax = dtmp;
          pmax = j;
        }
        if (dtmp < fmin)
        {
          fmin = dtmp;
          pmin = j;
        }
        /* calculate variance, skew, and kurtosis */
        rms += dtmp * dtmp;
        diff = dtmp - sum;
        dtmp = diff * diff;
        variance += dtmp;
        dtmp *= diff;
        skewness += dtmp;
        dtmp *= diff;
        kurtosis += dtmp;
      }
      variance /= (double)(k-1);
      stddev = sqrt(variance);
      if(variance != 0.0) {
        skewness /= ((double)k * pow(stddev,3.0));
        kurtosis = (kurtosis/((double)k * pow(variance,2.0))) - 3.0;
      }

      new_stats->rms = (float) sqrt(rms/(float)k);
      new_stats->pinteg = (float)pint;
      new_stats->ninteg = (float)nint;
      new_stats->pos = pos;
      new_stats->neg = neg;
      new_stats->var = (float) variance;
      new_stats->sdev = (float) stddev;
      new_stats->cpts = k;
      new_stats->fmax = (float) fmax;
      new_stats->maxy = (int)(pmax/img1->row_size);
      new_stats->maxx = pmax - (new_stats->maxy * img1->row_size);
      new_stats->maxp = pmax;   
      new_stats->fmin = (float) fmin;
      new_stats->miny = (int)(pmin/img1->row_size);
      new_stats->minx = pmin - (new_stats->miny * img1->row_size);
      new_stats->minp = pmin;   
      new_stats->skewness = (float) skewness;
      new_stats->kurtosis = (float) kurtosis;
      new_stats->entropy = entropy;
      new_stats->contrast = contrast;
    }
  }

  else  /* use masking image */
  {
    for (i = 0; i < num_vects1; i++)
    {
      /* Zero out variables */
      pos = 0;
      neg = 0;
      sum = 0.0;
      nint = 0.0;
      pint = 0.0;
      rms = 0.0;
      variance = 0.0;
      entropy = 0.0;
      contrast = 0.0;
      skewness = 0.0;
      kurtosis = 0.0;
      for(j=0; j<dim1; j++) {
        if(i2[i][j]) {
          fmin = (double) i1[i][j];
          fmax = (double) i1[i][j];
          pmin = j;
          pmax = j;
          j=dim1;
        }
      }
      k = 0;
      /* calculate mean */
      for (j=0; j<dim1; j++) {
        if (i2[i][j]) {
          sum += (double) i1[i][j];
          k++;
        }
      }
      if (k < 2) {
        (void)fprintf(stderr,"calculate_long: will not calculate statistics ");
        (void)fprintf(stderr,"on less than 2 points\n");
        (void)fprintf(stderr,"Data has only %d points per vector.\n",k);
        return(0);
      }
      new_stats = stats[i];
      new_stats->integ = (float)sum;
      sum /= (float)k;
      new_stats->mean = (float)sum;

      for (j=0; j<dim1; j++)
      {
        if (i2[i][j])
        {
          dtmp = (double) i1[i][j];
          /* calculate pos, neg integrals, num points, & min/max */
          if (dtmp >= 0) { pint += dtmp; pos++; }
          if (dtmp <  0) { nint += dtmp; neg++; }
          if (dtmp > fmax)
          {
            fmax = dtmp;
            pmax = j;
          }
          if (dtmp < fmin)
          {
            fmin = dtmp;
            pmin = j;
          }
          /* calculate variance, skew, and kurtosis */
          rms += dtmp * dtmp;
          diff = dtmp - sum;
          dtmp = diff * diff;
          variance += dtmp;
          dtmp *= diff;
          skewness += dtmp;
          dtmp *= diff;
          kurtosis += dtmp;
        }
      }
      variance /= (double)(k-1);
      stddev = sqrt(variance);
      if(variance != 0.0) {
        skewness /= ((double)k * pow(stddev,3.0));
        kurtosis = (kurtosis/((double)k * pow(variance,2.0))) - 3.0;
      }

      new_stats->rms = (float) sqrt(rms/(float)k);
      new_stats->pinteg = (float) pint;
      new_stats->ninteg = (float) nint;
      new_stats->pos = pos;
      new_stats->neg = neg;
      new_stats->var = (float) variance;
      new_stats->sdev = (float) stddev;
      new_stats->cpts = k;
      new_stats->fmax = (float) fmax;
      new_stats->maxy = (int)(pmax/img1->row_size);
      new_stats->maxx = pmax - (new_stats->maxy * img1->row_size);
      new_stats->maxp = pmax;   
      new_stats->fmin = (float) fmin;
      new_stats->miny = (int)(pmin/img1->row_size);
      new_stats->minx = pmin - (new_stats->miny * img1->row_size);
      new_stats->minp = pmin;   
      new_stats->skewness = (float) skewness;
      new_stats->kurtosis = (float) kurtosis;
      new_stats->entropy = entropy;
      new_stats->contrast = contrast;
    }
  }

  for (i=0; i<num_vects1; i++)
  {
    free(i1[i]);
    if(mflag) free(i2[i]);
  }
  free(i1);
  if(mflag) free(i2);

  return(1);
}

/**************************************************************
*
* MODULE NAME: calculate_float
*
*     PURPOSE: calculates the statistics for a FLOAT image. The
*              stats structure stats will be filled out.
*
*       INPUT:  1. img1 -- first xvimage strcuture
*               2. img2 -- gating mask image
*               3. mflag -- a flag set (equal 1) if gating image mask
*                          available
*               4. process_dir -- 0 = vector processing, 1 = band,
*                                 and 2 = calculate one set of stats
*                                 for entire file.
*
*               The input images must be of the same data type
*               and size.
*
*      OUTPUT:  
*               The stat structure is loaded with the values computed
*               from the image.
*
* CALLED FROM: lvstats
*
* ROUTINES CALLED: dload_vector() 
**************************************************************/

int calculate_float(img1, img2, stats, mflag, process_dir)

struct xvimage  *img1,*img2;
int             mflag, process_dir;
struct lvstat   **stats;

{
  int           pmax,pmin,      /* positions of min and max values in vector */
                pos,neg,        /* number of pos and neg points in vector */
                dim1, dim2,     /* vector dimension of input/masking files */
                num_vects1,     /* number of vectors in input file */ 
                num_vects2,     /* number of vectors in masking file */
                tmp_band,
                tmp_col_size,
                i, j, k;

  float         **f1, **f2;     /* image and masking vector arrays */

  double        skewness,       /* sum of cubed values in vector */
                kurtosis,       /* sum of values to the 4th pow in vector */
                variance,
                diff,
                stddev,
                sum,            /* sum of values in vector */
                rms,            /* sum of squared values in vector */
                pint,nint,      /* pos & neg integral sums */
                fmin,fmax,      /* min and max values in vector */
                entropy, contrast, 
                dtmp;

  struct lvstat *new_stats;

  char          **dload_vector(),
                *program = "calculate_float";


  if ( (process_dir) > 2 || (process_dir < 0) ) {
    (void)fprintf(stderr,"%s: unknown processing direction (%d)\n",
                  program, process_dir);
    return(0);
  }

  if (process_dir == 2)  /* calculate one set of stats for entire file */
  {
    tmp_band = img1->num_data_bands;
    tmp_col_size = img1->col_size;
    img1->col_size *= tmp_band;
    img1->num_data_bands = 1;
    if ((f1 = (float **)dload_vector(img1, &num_vects1, &dim1, DSP_BAND))
        == NULL)
    {
      (void)fprintf(stderr,"%s: dload_vector failed \n",program);
      return(0);
    }
    free(img1->imagedata);

    if (mflag)
    {
      img2->col_size *= img2->num_data_bands;
      img2->num_data_bands = 1;
      if ((f2 = (float **)dload_vector(img2, &num_vects2, &dim2, DSP_BAND))
          == NULL)
      {
        (void)fprintf(stderr,"%s: dload_vector failed\n",program);
        return(0);
      }
      if ( (num_vects1 != num_vects2) || (dim1 != dim2) )
      {
        (void)fprintf(stderr,"%s: incompatible mask size\n",program);
        return(0);
      }
      img2->num_data_bands = tmp_band;
      img2->col_size = tmp_col_size;
      free(img2->imagedata);
    }
    img1->num_data_bands = tmp_band;
    img1->col_size = tmp_col_size;
  }

  else  /* calculate stats on each band or vector separately */
  {
    if ((f1 = (float **)dload_vector(img1, &num_vects1, &dim1, process_dir))
        == NULL)
    {
      (void)fprintf(stderr,"%s: dload_vector failed\n",program);
      return(0);
    }
    free(img1->imagedata);
    if (mflag)
    {
      if ((f2 = (float **)dload_vector(img2, &num_vects2, &dim2, process_dir))
          == NULL)
      {
        (void)fprintf(stderr,"%s: dload_vector failed\n",program);
        return(0);
      }
      if ( (num_vects1 != num_vects2) || (dim1 != dim2) )
      {
        (void)fprintf(stderr,"%s: incompatible mask size\n",program);
        return(0);
      }
      free(img2->imagedata);
    }
  }

  if (dim1 < 2) {
    (void)fprintf(stderr,"calculate_float: will not calculate statistics ");
    (void)fprintf(stderr,"on less than 2 points\n");
    (void)fprintf(stderr,"Data has only %d points per vector.\n",dim1);
    return(0);
  }

  if (!mflag)
  {
    for (i = 0; i < num_vects1; i++)
    {
      /* Zero out variables */
      pos = 0;
      neg = 0;
      sum = 0.0;
      nint = 0.0;
      pint = 0.0;
      rms = 0.0;
      variance = 0.0;
      entropy = 0.0;
      contrast = 0.0;
      skewness = 0.0;
      kurtosis = 0.0;
      fmin = (double) f1[i][0];
      fmax = (double) f1[i][0];
      pmin = 0;
      pmax = 0;
      k = 0;
      /* calculate mean */
      for (j=0; j<dim1; j++) {
        sum += (double) f1[i][j];
        k++;
      }
      new_stats = stats[i];
      new_stats->integ = (float)sum;
      sum /= (float)k;
      new_stats->mean = (float)sum;

      for (j=0; j<dim1; j++)
      {
        dtmp = (double) f1[i][j];
        /* calculate pos & neg intregrals & num points, & min/max */
        if (dtmp >= 0) { pint += dtmp; pos++; }
        if (dtmp <  0) { nint += dtmp; neg++; }
        if (dtmp > fmax)
        {
          fmax = dtmp;
          pmax = j;
        }
        if (dtmp < fmin)
        {
          fmin = dtmp;
          pmin = j;
        }
        /* calculate variance, skew, and kurtosis */
        rms += dtmp * dtmp;
        diff = dtmp - sum;
        dtmp = diff * diff;
        variance += dtmp;
        dtmp *= diff;
        skewness += dtmp;
        dtmp *= diff;
        kurtosis += dtmp;
      }
      variance /= (double)(k-1);
      stddev = sqrt(variance);
      if(variance != 0.0) {
        skewness /= ((double)k * pow(stddev,3.0));
        kurtosis = (kurtosis/((double)k * pow(variance,2.0))) - 3.0;
      }

      new_stats->rms = (float) sqrt(rms/(double)k);
      new_stats->pinteg = (float)pint;
      new_stats->ninteg = (float)nint;
      new_stats->pos = pos;
      new_stats->neg = neg;
      new_stats->var = (float) variance;
      new_stats->sdev = (float) stddev;
      new_stats->cpts = k;
      new_stats->fmax = (float) fmax;
      new_stats->maxy = (int)(pmax/img1->row_size);
      new_stats->maxx = pmax - (new_stats->maxy * img1->row_size);
      new_stats->maxp = pmax;   
      new_stats->fmin = (float) fmin;
      new_stats->miny = (int)(pmin/img1->row_size);
      new_stats->minx = pmin - (new_stats->miny * img1->row_size);
      new_stats->minp = pmin;   
      new_stats->skewness = (float) skewness;
      new_stats->kurtosis = (float) kurtosis;
      new_stats->entropy = entropy;
      new_stats->contrast = contrast;
    }
  }


  else  /* use masking image */
  {
    for (i = 0; i < num_vects1; i++)
    {
      /* Zero out variables */
      pos = 0;
      neg = 0;
      sum = 0.0;
      nint = 0.0;
      pint = 0.0;
      rms = 0.0;
      variance = 0.0;
      entropy = 0.0;
      contrast = 0.0;
      skewness = 0.0;
      kurtosis = 0.0;
      for(j=0; j<dim1; j++) {
        if(f2[i][j]) {
          fmin = (double) f1[i][j];
          fmax = (double) f1[i][j];
          pmin = j;
          pmax = j;
          j=dim1;
        }
      }
      k = 0;
      /* calculate mean */
      for (j=0; j<dim1; j++) {
        if (f2[i][j]) {
          sum += (double) f1[i][j];
          k++;
        }
      }
      if (k < 2) {
        (void)fprintf(stderr,"calculate_float: will not calculate statistics ");
        (void)fprintf(stderr,"on less than 2 points\n");
        (void)fprintf(stderr,"Data has only %d points per vector.\n",k);
        return(0);
      }
      new_stats = stats[i];
      new_stats->integ = (float)sum;
      sum /= (float)k;
      new_stats->mean = (float)sum;

      for (j=0; j<dim1; j++)
      {
        if (f2[i][j])
        {
          dtmp = (double) f1[i][j];
          /* calculate pos, neg integrals, num points, & min/max */
          if (dtmp >= 0) { pint += dtmp; pos++; }
          if (dtmp <  0) { nint += dtmp; neg++; }
          if (dtmp > fmax)
          {
            fmax = dtmp;
            pmax = j;
          }
          if (dtmp < fmin)
          {
            fmin = dtmp;
            pmin = j;
          }
          /* calculate variance, skew, and kurtosis */
          rms += dtmp * dtmp;
          diff = dtmp - sum;
          dtmp = diff * diff;
          variance += dtmp;
          dtmp *= diff;
          skewness += dtmp;
          dtmp *= diff;
          kurtosis += dtmp;
        }
      }
      variance /= (double)(k-1);
      stddev = sqrt(variance);
      if(variance != 0.0) {
        skewness /= ((double)k * pow(stddev,3.0));
        kurtosis = (kurtosis/((double)k * pow(variance,2.0))) - 3.0;
      }

      new_stats->rms = (float) sqrt(rms/(double)k);
      new_stats->pinteg = (float) pint;
      new_stats->ninteg = (float) nint;
      new_stats->pos = pos;
      new_stats->neg = neg;
      new_stats->var = (float) variance;
      new_stats->sdev = (float) stddev;
      new_stats->cpts = k;
      new_stats->fmax = (float) fmax;
      new_stats->maxy = (int)(pmax/img1->row_size);
      new_stats->maxx = pmax - (new_stats->maxy * img1->row_size);
      new_stats->maxp = pmax;   
      new_stats->fmin = (float) fmin;
      new_stats->miny = (int)(pmin/img1->row_size);
      new_stats->minx = pmin - (new_stats->miny * img1->row_size);
      new_stats->minp = pmin;   
      new_stats->skewness = (float) skewness;
      new_stats->kurtosis = (float) kurtosis;
      new_stats->entropy = entropy;
      new_stats->contrast = contrast;
    }
  }

  for (i=0; i<num_vects1; i++)
  {
    free(f1[i]);
    if(mflag) free(f2[i]);
  }
  free(f1);
  if(mflag) free(f2);

  return(1);
}

/**************************************************************
*
* MODULE NAME: calculate_double
*
*     PURPOSE: calculates the statistics for a DOUBLE image. The
*              stats structure stats will be filled out.
*
*       INPUT:  1. img1 -- first xvimage strcuture
*               2. img2 -- gating mask image
*               3. mflag -- a flag set (equal 1) if gating image mask
*                          available
*               4. process_dir -- 0 = vector processing, 1 = band,
*                                 and 2 = calculate one set of stats
*                                 for entire file.
*
*               The input images must be of the same data type
*               and size.
*
*      OUTPUT:  
*               The stat structure is loaded with the values computed
*               from the image.
*
* CALLED FROM: lvstats
*
* ROUTINES CALLED: dload_vector() 
**************************************************************/

int calculate_double(img1, img2, stats, mflag, process_dir)

struct xvimage  *img1,*img2;
int             mflag, process_dir;
struct lvstat   **stats;

{
  int           pmax,pmin,      /* positions of min and max values in vector */
                pos,neg,        /* number of pos and neg points in vector */
                dim1, dim2,     /* vector dimension of input/masking files */
                num_vects1,     /* number of vectors in input file */ 
                num_vects2,     /* number of vectors in masking file */
                tmp_band,
                tmp_col_size,
                i, j, k;

  double        **d1, **d2;     /* image and masking vector arrays */

  double        skewness,       /* sum of cubed values in vector */
                kurtosis,       /* sum of values to the 4th pow in vector */
                variance,
                diff,
                stddev,
                sum,            /* sum of values in vector */
                rms,            /* sum of squared values in vector */
                pint,nint,      /* pos & neg integral sums */
                fmin,fmax,      /* min and max values in vector */
                entropy, contrast, 
                dtmp;

  struct lvstat *new_stats;

  char          **dload_vector(),
                *program = "calculate_float";


  if ( (process_dir) > 2 || (process_dir < 0) ) {
    (void)fprintf(stderr,"%s: unknown processing direction (%d)\n",
                  program, process_dir);
    return(0);
  }

  if (process_dir == 2)  /* calculate one set of stats for entire file */
  {
    tmp_band = img1->num_data_bands;
    tmp_col_size = img1->col_size;
    img1->col_size *= tmp_band;
    img1->num_data_bands = 1;
    if ((d1 = (double **)dload_vector(img1, &num_vects1, &dim1, DSP_BAND))
        == NULL)
    {
      (void)fprintf(stderr,"%s: dload_vector failed \n",program);
      return(0);
    }
    free(img1->imagedata);

    if (mflag)
    {
      img2->col_size *= img2->num_data_bands;
      img2->num_data_bands = 1;
      if ((d2 = (double **)dload_vector(img2, &num_vects2, &dim2, DSP_BAND))
          == NULL)
      {
        (void)fprintf(stderr,"%s: dload_vector failed\n",program);
        return(0);
      }
      if ( (num_vects1 != num_vects2) || (dim1 != dim2) )
      {
        (void)fprintf(stderr,"%s: incompatible mask size\n",program);
        return(0);
      }
      img2->num_data_bands = tmp_band;
      img2->col_size = tmp_col_size;
      free(img2->imagedata);
    }
    img1->num_data_bands = tmp_band;
    img1->col_size = tmp_col_size;
  }

  else  /* calculate stats on each band or vector separately */
  {
    if ((d1 = (double **)dload_vector(img1, &num_vects1, &dim1, process_dir))
        == NULL)
    {
      (void)fprintf(stderr,"%s: dload_vector failed\n",program);
      return(0);
    }
    free(img1->imagedata);
    if (mflag)
    {
      if ((d2 = (double **)dload_vector(img2, &num_vects2, &dim2, process_dir))
          == NULL)
      {
        (void)fprintf(stderr,"%s: dload_vector failed\n",program);
        return(0);
      }
      if ( (num_vects1 != num_vects2) || (dim1 != dim2) )
      {
        (void)fprintf(stderr,"%s: incompatible mask size\n",program);
        return(0);
      }
      free(img2->imagedata);
    }
  }

  if (dim1 < 2) {
    (void)fprintf(stderr,"calculate_double: will not calculate statistics ");
    (void)fprintf(stderr,"on less than 2 points\n");
    (void)fprintf(stderr,"Data has only %d points per vector.\n",dim1);
    return(0);
  }

  if (!mflag)
  {
    for (i = 0; i < num_vects1; i++)
    {
      /* Zero out variables */
      pos = 0;
      neg = 0;
      sum = 0.0;
      nint = 0.0;
      pint = 0.0;
      rms = 0.0;
      variance = 0.0;
      entropy = 0.0;
      contrast = 0.0;
      skewness = 0.0;
      kurtosis = 0.0;
      fmin = d1[i][0];
      fmax = d1[i][0];
      pmin = 0;
      pmax = 0;
      k = 0;
      /* calculate mean */
      for (j=0; j<dim1; j++) {
        sum += d1[i][j];
        k++;
      }
      new_stats = stats[i];
      new_stats->integ = (float)sum;
      sum /= (float)k;
      new_stats->mean = (float)sum;

      for (j=0; j<dim1; j++)
      {
        dtmp = d1[i][j];
        /* calculate pos & neg intregrals & num points, & min/max */
        if (dtmp >= 0) { pint += dtmp; pos++; }
        if (dtmp <  0) { nint += dtmp; neg++; }
        if (dtmp > fmax)
        {
          fmax = dtmp;
          pmax = j;
        }
        if (dtmp < fmin)
        {
          fmin = dtmp;
          pmin = j;
        }
        /* calculate variance, skew, and kurtosis */
        rms += dtmp * dtmp;
        diff = dtmp - sum;
        dtmp = diff * diff;
        variance += dtmp;
        dtmp *= diff;
        skewness += dtmp;
        dtmp *= diff;
        kurtosis += dtmp;
      }
      variance /= (double)(k-1);
      stddev = sqrt(variance);
      if(variance != 0.0) {
        skewness /= ((double)k * pow(stddev,3.0));
        kurtosis = (kurtosis/((double)k * pow(variance,2.0))) - 3.0;
      }

      new_stats->rms = (float) sqrt(rms/(double)k);
      new_stats->pinteg = (float)pint;
      new_stats->ninteg = (float)nint;
      new_stats->pos = pos;
      new_stats->neg = neg;
      new_stats->var = (float) variance;
      new_stats->sdev = (float) stddev;
      new_stats->cpts = k;
      new_stats->fmax = (float) fmax;
      new_stats->maxy = (int)(pmax/img1->row_size);
      new_stats->maxx = pmax - (new_stats->maxy * img1->row_size);
      new_stats->maxp = pmax;   
      new_stats->fmin = (float) fmin;
      new_stats->miny = (int)(pmin/img1->row_size);
      new_stats->minx = pmin - (new_stats->miny * img1->row_size);
      new_stats->minp = pmin;   
      new_stats->skewness = (float) skewness;
      new_stats->kurtosis = (float) kurtosis;
      new_stats->entropy = entropy;
      new_stats->contrast = contrast;
    }
  }


  else  /* use masking image */
  {
    for (i = 0; i < num_vects1; i++)
    {
      /* Zero out variables */
      pos = 0;
      neg = 0;
      sum = 0.0;
      nint = 0.0;
      pint = 0.0;
      rms = 0.0;
      variance = 0.0;
      entropy = 0.0;
      contrast = 0.0;
      skewness = 0.0;
      kurtosis = 0.0;
      for(j=0; j<dim1; j++) {
        if(d2[i][j]) {
          fmin = d1[i][j];
          fmax = d1[i][j];
          pmin = j;
          pmax = j;
          j=dim1;
        }
      }
      k = 0;
      /* calculate mean */
      for (j=0; j<dim1; j++) {
        if (d2[i][j]) {
          sum += d1[i][j];
          k++;
        }
      }
      if (k < 2) {
        (void)fprintf(stderr,"calculate_double: will not calculate statistics");
        (void)fprintf(stderr," on less than 2 points\n");
        (void)fprintf(stderr,"Data has only %d points per vector.\n",k);
        return(0);
      }
      new_stats = stats[i];
      new_stats->integ = (float)sum;
      sum /= (float)k;
      new_stats->mean = (float)sum;

      for (j=0; j<dim1; j++)
      {
        if (d2[i][j])
        {
          dtmp = d1[i][j];
          /* calculate pos, neg integrals, num points, & min/max */
          if (dtmp >= 0) { pint += dtmp; pos++; }
          if (dtmp <  0) { nint += dtmp; neg++; }
          if (dtmp > fmax)
          {
            fmax = dtmp;
            pmax = j;
          }
          if (dtmp < fmin)
          {
            fmin = dtmp;
            pmin = j;
          }
          /* calculate variance, skew, and kurtosis */
          rms += dtmp * dtmp;
          diff = dtmp - sum;
          dtmp = diff * diff;
          variance += dtmp;
          dtmp *= diff;
          skewness += dtmp;
          dtmp *= diff;
          kurtosis += dtmp;
        }
      }
      variance /= (double)(k-1);
      stddev = sqrt(variance);
      if(variance != 0.0) {
        skewness /= ((double)k * pow(stddev,3.0));
        kurtosis = (kurtosis/((double)k * pow(variance,2.0))) - 3.0;
      }

      new_stats->rms = (float) sqrt(rms/(double)k);
      new_stats->pinteg = (float) pint;
      new_stats->ninteg = (float) nint;
      new_stats->pos = pos;
      new_stats->neg = neg;
      new_stats->var = (float) variance;
      new_stats->sdev = (float) stddev;
      new_stats->cpts = k;
      new_stats->fmax = (float) fmax;
      new_stats->maxy = (int)(pmax/img1->row_size);
      new_stats->maxx = pmax - (new_stats->maxy * img1->row_size);
      new_stats->maxp = pmax;   
      new_stats->fmin = (float) fmin;
      new_stats->miny = (int)(pmin/img1->row_size);
      new_stats->minx = pmin - (new_stats->miny * img1->row_size);
      new_stats->minp = pmin;   
      new_stats->skewness = (float) skewness;
      new_stats->kurtosis = (float) kurtosis;
      new_stats->entropy = entropy;
      new_stats->contrast = contrast;
    }
  }

  for (i=0; i<num_vects1; i++)
  {
    free(d1[i]);
    if(mflag) free(d2[i]);
  }
  free(d1);
  if(mflag) free(d2);

  return(1);
}

/* -library_code_end */
