/*----------------------------------------------------------------------PROGRAM: contrastmap.cDATE:    6/28/94AUTHOR:  Baback Moghaddam, baback@media.mit.edu------------------------------------------------------------------------   This routine reads in an ASCII file of the format   filename     .   .   .   and contrast normalizes each image according to the following   1. compute the (mean,std.dev.) of (possibly non-zero) pixels   2. map the image into zero-mean unit-sigma    3. output image is this normalized image multiplied by a range      factor and offset by an pedestal factor   4. if specified, the output will be written as uchar   It outputs the transformed files in an output directory   with the same names   ---------------------------------------------------------------------- */#include <stdio.h>#include <stdlib.h>#include <string.h>#include <math.h>#include <float.h>#include "util.h"#include "io.h"int contrast(float **image_in, float **image_out, 	     int nrow, int ncol,	     float pedestal, float scale,	     float *mean, float *sigma);/* ----------- Command-Line Parsing Stuff ------- */extern int optind;extern char *optarg;char *progname;          /* used to store the name of the program  */char comline[256];       /* used to store the entire command line  */#define OPTIONS "i:l:o:s:p:fz"char *usage = "\t-i indir -l list -o outdir \n\\t\t\t-s scale -p pedestal [-f] [-z]\n";char *help="\Contrast Normalization\n\n\-i indir    \t input directory\n\-l list     \t ASCII list of filenames (one per line)\n\-o outdir   \t output directory\n\-s scale    \t scale factor (multiplies the normalized image)\n\-p pedestal \t pedestal (offset for scaled normalized image)\n\-f          \t write float output image (default is same as input)\n\-z          \t computes statistics over non-zero portion of input\n";/* ----- global command-line defaults ----- */int out_float_flag = -1;  /* -1 = unspecified */int nonzero_flag = 0;#define MAX_CHARS  256/*----------------------------------------------------------------------*/main(int argc, char **argv){   register int i,j,k,l,ii,jj;  int f, c, nframe, nfeatures, sets, bytes_pixel;  int nrow, ncol, M, N;  char command[MAX_CHARS],indir[MAX_CHARS],listfile[MAX_CHARS], \    infile[MAX_CHARS], modelfile[MAX_CHARS],outdir[MAX_CHARS], \      filename[MAX_CHARS], line[MAX_CHARS];  float **image_in, **image_out;  unsigned char **uchar_image;  float fval1, fval2, fval;  FILE *fp, *fp2;         /* for output values dump */    float scale;  float pedestal;  float mean;  float sigma;  int npixels;  /* required input flags */    int errflag   = 0;  int inflag    = 0;  int listflag  = 0;  int outflag   = 0;  int s_flag    = 0;  int p_flag    = 0;  progname = argv[0];  for (i=0; i<argc; i++)    strcat(comline, argv[i]),strcat(comline, " ");      /* ----------------------  Command Line Parse ------------------------ */    while ((c = getopt(argc, argv, OPTIONS)) != EOF)    switch (c) {          case 'i':      strcpy(indir, optarg);      inflag = 1;      break;    case 'l':      strcpy(listfile, optarg);      listflag = 1;      break;          case 'o':      strcpy(outdir, optarg);      outflag = 1;      break;    case 's':      scale = atof(optarg);      s_flag = 1;      break;    case 'p':      pedestal = atof(optarg);      p_flag = 1;      break;    case 'f':      out_float_flag = 1;      break;    case 'z':      nonzero_flag = 1;      break;    case '?':      errflag = 1;      break;    }    /* command line error check */    if (errflag || !inflag || !listflag || !outflag || !s_flag || !p_flag) {    fprintf(stderr,"\nUSAGE: %s %s\n%s\n", progname, usage, help);    exit(1);  }  /* ----  read indir descriptor file -------- */    read_descriptor(indir, &nframe, &sets, &bytes_pixel, &ncol, &nrow);  if (sets>1)     myerror("Input files must be single-set DAT files!");  if (out_float_flag == -1)    out_float_flag = (bytes_pixel==4 ? 1:0);  /* ---- set up matrices  ---- */  image_in   = matrix(1, nrow, 1, ncol);  image_out  = matrix(1, nrow, 1, ncol);  if (out_float_flag == 0)    uchar_image = cmatrix(1, nrow, 1, ncol);  /* ---- loop over input list file and warp ------- */  if ((fp = fopen(listfile, "r")) == NULL) {    fprintf(stderr,"ERROR: Could not open input file %s \n\n", listfile);    exit(1);  }  nframe = 0;  while (strlen(fgets(line, MAX_CHARS, fp))) {    if (strncmp(line, "#", 1) != 0 && strlen(line)>1) {      nframe++;          /* --- read filename --- */            sscanf(line, "%s", infile);      /* --- read file --- */            sprintf(filename,"%s/%s", indir, infile);      if (bytes_pixel == 4)	read_RAW_float(filename, image_in, nrow, ncol);      else	read_RAW(filename, image_in, nrow, ncol);      fprintf(stdout, "Read %d-by-%d %s file %s\n",	      nrow, ncol, (bytes_pixel==4 ? "float":"uchar"), infile);                  /* --- normalize the contrast --- */            npixels = contrast(image_in, image_out, nrow, ncol, 	       pedestal, scale, &mean, &sigma);      fprintf(stdout,"Npixels = %d\nMean    = %f\nSigma   = %f\n",	      npixels, mean, sigma);                  /* --- write output image --- */            if (out_float_flag == 0) {	for (i=1; i<=nrow; i++) {	  for (j=1; j<=ncol; j++) {	    if (image_out[i][j]>255)	      uchar_image[i][j] = 255;	    else if (image_out[i][j]<0)	      uchar_image[i][j] = 0;	    else 	      uchar_image[i][j] = (int) (image_out[i][j] + 0.5);	  }	}      }      sprintf(filename,"%s/%s", outdir, infile);      if (out_float_flag == 1)	write_RAW_float(filename, image_out, nrow, ncol);      else	write_RAW(filename, uchar_image, nrow, ncol);      fprintf(stdout,"Wrote %d-by-%d %s output file %s\n\n", 	      nrow, ncol, (out_float_flag==1 ? "float":"uchar"),	      filename);                }  }  /* --- write output dir descriptor --- */    write_descriptor(outdir, nframe, ncol, nrow,		   ((bytes_pixel==1 || out_float_flag==1) ? 1:4),		   comline);  free_matrix(image_in, 1, nrow, 1, ncol);  free_matrix(image_out, 1, nrow, 1, ncol);  if (out_float_flag == 0)    free_cmatrix(uchar_image, 1, nrow, 1, ncol);  fclose(fp);  return 0;}/* ====================================================================== */int contrast(float **image_in, float **image_out, 	      int nrow, int ncol,	      float pedestal, float scale,	      float *mean, float *sigma){  register int i,j;  float sum;  float my_mean, my_sigma;  int count;  /* --- compute statistics ---- */  sum = count = 0;  for (i=1; i<=nrow; i++)    for (j=1; j<=ncol; j++) {      if (nonzero_flag == 1) {	if (image_in[i][j] != 0) {	  sum += image_in[i][j];	  count ++;	}      }      else {	sum += image_in[i][j];	count ++;      }    }  my_mean = sum/count;  sum = count = 0;  for (i=1; i<=nrow; i++)    for (j=1; j<=ncol; j++) {      if (nonzero_flag == 1) {	if (image_in[i][j] != 0) {	  sum += SQR(image_in[i][j] - my_mean);	  count ++;	}      }      else {	sum += SQR(image_in[i][j] - my_mean);	count ++;      }    }  my_sigma = sqrt(sum/(count-1));  /* --- remap input image ---- */  for (i=1; i<=nrow; i++)    for (j=1; j<=ncol; j++)       if (nonzero_flag == 1) {	if (image_in[i][j] != 0)	  image_out[i][j] = pedestal + 	    scale * (image_in[i][j]-my_mean)/my_sigma;	else	  image_out[i][j] = 0;      }  *mean = my_mean;  *sigma = my_sigma;  return count;}