/* * * * * * * * * * * * * * * * * * *  * * * * * * * * * * * * * * * * * *
 *                                                                        *
 * rs_8_16:    Calculate statistics on BIL datasets.                      *
 *                                                                        *
 *         -   Input parameters (at command line)are:                     *
 *                 -        Dataset filename                              *
 *                 -        Rows                                          *
 *                 -        Columns                                       *
 *                 -        Bands         and                             *
 *                 -        Header size                                   *
 *         -   Output parameters (on output) for each band are:           *
 *                 -        Minimum                                       *
 *                 -        Maximum                                       *
 *                 -        Average                                       *
 *                 -        Median        and                             *
 *                 -        Std. Deviation                                *
 *         -   Works on 8 or 16 bit unsigned integer datasets             *
 *                                                                        *
 * * * * * * * * * * * * * * * * * * * *  * * * * * * * * * * * * * * * * *
 *                                                                        *
 * rs_8_16  by  Christos Iossifides (c)  January 2003 NaTUReS Lab.        *
 *                                                                        *
 * * * * * * * * * * * * * * * * * * *  * * * * * * * * * * * * * * * * * */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

int main(int argc, char *argv[]){
  FILE *fp;
  int row, col, ban, hea, half, sum, d, pixelsize;
  int *mini, *maxi, *histo, *med;
  unsigned char *p;
  unsigned short int *pp;
  double *s, *ss;
  register int r, b, c;

  if(argc!=6){
    puts("Usage: >rs_8_16 filename rows columns bands header.");
    exit(1);
    }
  if((fp=fopen(argv[1],"rb"))==NULL){
    printf("rs_8_16 reports: Bad filename <%s>.\n",argv[1]);
    exit(1);
    }
  row=atoi(argv[2]);
  col=atoi(argv[3]);
  ban=atoi(argv[4]);
  hea=atoi(argv[5]);
  if (row<1 ||col<1|| ban<1|| hea<0|| hea>=row*col*ban*sizeof(short)){
    printf("rs_8_16 reports: Bad data <%d %d %d %d>.\n",row,col,ban,hea);
    exit(1);
    }
  if (fseek(fp,0,SEEK_END)){
    puts("rs_8_16 reports: Bad seek.");
    exit(1);
    }
  pixelsize=(ftell(fp)-hea)/row/col/ban;
  if (fseek(fp,0,SEEK_SET)){
    puts("rs_8_16 reports: Bad seek.");
    exit(1);
    }
  if (pixelsize<1||pixelsize>2){
    printf("rs_8_16 reports: Bad pixel size <%d %d %d %d -> %d>.\n",row,col,ban,hea,pixelsize);
    exit(1);
    }
  if(hea&&fseek(fp,hea,SEEK_SET)){
    puts("rs_8_16 reports: Bad seek.");
    exit(1);
    }
  if((pp=(unsigned short int *)p=(unsigned char *)malloc(col*ban*pixelsize))==NULL){
    puts("rs_8_16 reports: Bad memory allocation.");
    exit(1);
    }
  if((mini=malloc(ban*sizeof(int)))==NULL){
    puts("rs_8_16 reports: Bad memory allocation.");
    exit(1);
    }
  if((maxi=malloc(ban*sizeof(int)))==NULL){
    puts("rs_8_16 reports: Bad memory allocation.");
    exit(1);
    }
  if((med=malloc(ban*sizeof(int)))==NULL){
    puts("rs_8_16 reports: Bad memory allocation.");
    exit(1);
    }
  if((s=malloc(ban*sizeof(double)))==NULL){
    puts("rs_8_16 reports: Bad memory allocation.");
    exit(1);
    }
  if((ss=malloc(ban*sizeof(double)))==NULL){
    puts("rs_8_16 reports: Bad memory allocation.");
    exit(1);
    }
  if((histo=calloc(ban*65535,sizeof(int)))==NULL){
    puts("rs_8_16 reports: Bad memory allocation.");
    exit(1);
    }
  for (b=0;b<ban;b++){
    mini[b]=65536;
    maxi[b]=0;
    ss[b]=s[b]=0.0;
    }
  for (r=0;r<row;r++){
    if(fread(p,pixelsize,col*ban,fp)!=col*ban){
      puts("rs_8_16 reports: EOF encountered.");
      exit(1);
      }
    for (b=0;b<ban;b++)
      for (c=0;c<col;c++){
        s[b]+=d=(int)(pixelsize-1)?pp[b*col+c]:p[b*col+c];
        ss[b]+=d*d;
        histo[b*65535+d]++;
        if(mini[b]>(int)d)
          mini[b]=(int)d;
        if(maxi[b]<(int)d)
          maxi[b]=(int)d;
        }
    }

  fclose(fp);  
  half=row*col/2;
  for (b=0;b<ban;b++){
    sum=0;
    for (c=mini[b];sum<half;r=c++)
      sum+=histo[b*65535+c];
    med[b]=r;
    }
  printf("\n Band | Minimum | Maximum |  Average |  Median | Std. Deviation");
  printf("\n------+---------+---------+----------+---------+----------------");
  for (b=0;b<ban;b++)
    printf("\n %4d | %7d | %7d | %8.3lf | %7d |  %12.8lf",
               b+1, mini[b], maxi[b], s[b]/row/col, med[b], sqrt(ss[b]/row/col-s[b]*s[b]/row/row/col/col));

  printf("\n\nRS_8_16 by Ch Iossif (c) Jan 2003 NaTUReS Lab.\n");

  return 0;
}