/* -*- Mode: C; indent-tabs-mode: nil; c-basic-offset: 4 c-style: "K&R" -*- */

/*---------------------------------------------------------------------------
   t-avg-img: calculates time-averaged values from a series if images

   Copyright (C) 2002, 2003, 2004 Gerber van der Graaf

   T-avg-img is free software; you can redistribute it and/or modify
   it under the terms of the GNU General Public License as published by
   the Free Software Foundation; either version 2, or (at your option)
   any later version.

   This program is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   GNU General Public License for more details.

   You should have received a copy of the GNU General Public License
   along with this program; if not, write to the Free Software Foundation,
   Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.  



TODO:

BUGS:

RCS version control
   $Log: t-avg-img.c,v $
   Revision 2.6  2006/01/31 14:18:03  gerber
   version 0.3.2

   Revision 2.5  2005/06/15 09:32:20  gerber
   check on memory leacks: change of internal structure; gpiv_alloc_pivdata and gpiv_free_pivdata in main optimizing by printing progress only if percentage has been changed

   Revision 2.4  2005/02/26 16:29:19  gerber
   parameter flags (parameter_logic) defined as gboolean, added scale_type to struct __GpivPostPar for inverse scaling

   Revision 2.3  2004/10/15 10:11:41  gerber
   GPIV_ and Gpiv prefix to defines and structure names of libgpiv

   Revision 2.2  2004/06/14 14:39:26  gerber
   addied image filtering and processing program

   Revision 2.1  2003/03/11 21:41:03  gerber
   GePeLled

   Revision 2.0  2002/12/03 12:44:40  gerber
   update version number

   Revision 1.1  2002/12/02 23:36:46  gerber
   init



------------------------------------------------------------------------*/

#include <stdlib.h>
#include <stdio.h>
#include <glib.h>
#include <gpiv.h>

#define PARFILE "t-avg-img.par"	/* Parameter file name */
#define USAGE "\
Usage: imgt-avg [-c int] [-dx] [-fb file] [-ff int] [fl int] [-h] [-p] [-r int] \n\
[-s] [-v] < stdin > stdout \n\
\n\
keys: \n\
\n\
-c  ncols:    number of columns  \n\
-fb file:     file base-name (without .number.r extension) instead of stdin \n\
              and stdout \n\
-ff N:        number of first image file \n\
-fl N:        number of last image file \n\
-fx:          prefix numbering to file base name \n\
-p:           print parameters to stdout \n\
-r nrows:     number of rows  \n\
-s:           subtract mean from input data  \n\
-v:           version number \n\
"

/*
 * For the moment: only from a single directory
 */
/* [-db dir] [-df int] [-dl int] \n\ */
/* db dir        directory base- name (without ending number) \n\ */
/* df N          number of first direcory \n\ */
/* dl N          number of first direcory \n\ */
/* dx             prefix numbering to directory base name \n\ */

#define HELP  "\
Calculates time-averaged intensities from a series of images for each pixel. \n\
Images should be numbered at the start or at the end of its name"

#define RCSID "$Id: t-avg-img.c,v 2.6 2006/01/31 14:18:03 gerber Exp $"

/*---------- Global variables ---------------------------------------------*/
gboolean fname_logic = FALSE;
gboolean print_par = FALSE;


/*
 * Function prototypes 
 */

void 
command_args(int argc, char *argv[], 
             char fname[GPIV_MAX_CHARS],
             GpivImagePar * image_par, 
             GpivImageProcPar * image_proc_par
             )
/* ----------------------------------------------------------------------------
 * Command line argument handling
 */
{
    char c;
    int argc_next;
    while (--argc > 0 && (*++argv)[0] == '-') {
	argc_next = 0;
/*
 * argc_next is set to 1 if the next cmd line argument has to be searched for; 
 * in case that the command line argument concerns more than one char or cmd 
 * line argument needs a parameter
 */
	while (argc_next == 0 && (c = *++argv[0])) {
	    switch (c) {
	    case 'v':
/*
 * Use Revision Control System (RCS) for version
 */
		printf("\n%s\n", RCSID);
		exit(0);
		break;
	    case 'h':
		printf("\n%s", RCSID);
		printf("\n%s", HELP);
		printf("\n%s", USAGE);
		exit(0);
		break;
/*
 * number of image columns
 */
	    case 'c':
		 image_par->ncolumns = atoi(*++argv);
		 image_par->ncolumns_logic = TRUE;
		 argc_next = 1;
		 --argc;
/*
 * direcory name and numbers
 */
	    case 'd':
	      if (strcmp("db", *argv) == 0) {
		strcpy(fname, *++argv);
		fname_logic = TRUE;
		argc_next = 1;
		--argc;
	      } else if (strcmp("df", *argv) == 0) {
		    image_proc_par->first_dir = atoi(*++argv);
		    image_proc_par->first_dir_logic = TRUE;
		    --argc;
		    argc_next = 1;
	      } else if (strcmp("dl", *argv) == 0) {
		    image_proc_par->last_dir = atoi(*++argv);
		    image_proc_par->last_dir_logic = TRUE;
		    --argc;
		    argc_next = 1;
	      } else if (strcmp("dx", *argv) == 0) {
		   image_proc_par->dir_prefix = 1;
	      }
		break;
/*
 * file name and numbers
 */
	    case 'f':
	      if (strcmp("fb", *argv) == 0) {
		strcpy(fname, *++argv);
		fname_logic = TRUE;
		argc_next = 1;
		--argc;
	      } else if (strcmp("ff", *argv) == 0) {
		    image_proc_par->first_file = atoi(*++argv);
		    image_proc_par->first_file_logic = TRUE;
		    --argc;
		    argc_next = 1;
	      } else if (strcmp("fl", *argv) == 0) {
		    image_proc_par->last_file = atoi(*++argv);
		    image_proc_par->last_file_logic = TRUE;
		    --argc;
		    argc_next = 1;
	      } else if (strcmp("fx", *argv) == 0) {
		   image_proc_par->file_prefix = 1;
	      }
		break;


	    case 'p':
		    print_par = TRUE;
		break;
/*
 * number of image rows
 */
	    case 'r':
		 image_par->nrows = atoi(*++argv);
		 image_par->nrows_logic = TRUE;
		 --argc;
		 argc_next = 1;
		 break;

/*
 * subtract mean from input data
 */
	    case 's':
		    image_proc_par->operator = GPIV_SUBTRACT;
		break;

	    default:
		gpiv_error("%s: unknown option: %s", RCSID, *argv);
		break;
	    }
	}
    }

    if (argc != 0) {
	gpiv_error("%s: unknown argument: %s", RCSID, *argv);
    }


}



void 
make_fname_out(char *fname_base, 
               char *fname_header, 
               char *fname_parameter,
               char *fname_out
               )
/*-----------------------------------------------------------------------------
 * generates filenames for output
 */
{

  char f_dum[GPIV_MAX_CHARS];
    if (fname_logic == FALSE) gpiv_error("%s: Filename has to be set", RCSID);

    gpiv_io_make_fname(fname_base, GPIV_EXT_HEADER, fname_header);
    if (print_par) printf("# Data parameter file: %s\n", fname_header);

    gpiv_io_make_fname(fname_base, GPIV_EXT_PAR, fname_parameter);
    if (print_par) printf("# Data parameter file: %s\n", fname_parameter);

    gpiv_io_make_fname(fname_base, GPIV_EXT_TA, f_dum);
    gpiv_io_make_fname(f_dum, GPIV_EXT_RAW_IMAGE, fname_out);
    if (print_par) printf("# Output file: %s\n", fname_out);

}



void 
make_fname_in(GpivImageProcPar image_proc_par, 
              char *fname_base, 
              int f_number, 
              char *fname_in
              )
/*-----------------------------------------------------------------------------
 * generates filenames for input
 */
{
  char f_dum[GPIV_MAX_CHARS];
    if (fname_logic == FALSE) gpiv_error("%s: Filename has to be set", RCSID);

    if (image_proc_par.file_prefix) {
        snprintf(f_dum, GPIV_MAX_CHARS, "%d%s", f_number, fname_base);
        gpiv_io_make_fname(f_dum, GPIV_EXT_RAW_IMAGE, fname_in);
    } else {
        snprintf(f_dum, GPIV_MAX_CHARS, "%s%d", fname_base, f_number);
        gpiv_io_make_fname(f_dum, GPIV_EXT_RAW_IMAGE, fname_in);
    }
      
    if (print_par) printf("# Input file: %s\n", fname_in);
     fprintf(stderr, "MAKE_FNAME_IN:: leaving\n");

}



void
gpiv_img_mean(GpivImagePar image_par, 
	      GpivImageProcPar image_proc_par, 
	      char *fname_base,
	      guint16 **img_mean
	      /*guint16 **img_sdev */
              )
/*-----------------------------------------------------------------------------
*/
{
     guint16 **img_in1 = NULL, **img_in2 = NULL;
     float **img_sum = NULL;
/*        guint16 **img_sumsqr = NULL;; */
     int /* i, */ j, k, l;
/*      int first_dir = image_proc_par.first_dir; */
/*      int last_dir = image_proc_par.last_dir; */
     int first_file = image_proc_par.first_file;
     int last_file = image_proc_par.last_file;

     int nrows = image_par.nrows;
     int ncolumns = image_par.ncolumns;
     char fname_in[GPIV_MAX_CHARS];

     if (print_par) printf("# Calculating mean from input data\n");
     img_sum = gpiv_matrix(image_par.nrows, image_par.ncolumns);
     memset(img_sum[0], 0, (sizeof(float)) * nrows * ncolumns);
     img_in1 = gpiv_alloc_img(image_par);
     img_in2 = gpiv_alloc_img(image_par);
     
/*
 * For the moment: only from single directory and only mean value
 */
/*      img_sumsqr = gpiv_alloc_img(image_par); */
/*     for (i = first_dir; i <= last_dir; i++) { */

     for (j = first_file; j <= last_file; j++) {
	 memset(img_in1[0], 0, (sizeof(guint16)) * nrows * ncolumns);
	 memset(img_in2[0], 0, (sizeof(guint16)) * nrows * ncolumns);
         if (fname_logic == TRUE) {
             make_fname_in(image_proc_par, fname_base, j, fname_in);
             gpiv_fread_image(fname_in, img_in1, img_in2, image_par);
         } else {
             gpiv_read_image(img_in1, img_in2, image_par);
         }
         
	 for (k = 0; k < nrows; k++) {
             for (l = 0; l < ncolumns; l++) {
                 img_sum[k][l] += 
                     ((float)img_in1[k][l] + (float)img_in1[k][l])/2.0;
             }
	 }
     }
    
     for (k = 0; k < nrows; k++) {
	 for (l = 0; l < ncolumns; l++) {
             float mean = 0.0;
             mean = img_sum[k][l] / (float) (last_file - first_file + 1);
             img_mean[k][l] = (guint16) mean;
	 }
     }
     
     gpiv_free_img(img_in1, image_par);
     gpiv_free_img(img_in2, image_par);
     gpiv_free_matrix(img_sum);
     
}




void
gpiv_img_subtract(GpivImagePar image_par, 
                  GpivImageProcPar image_proc_par, 
                  char *fname_base,
                  guint16 **img_subtr
                  )
/*-----------------------------------------------------------------------------
*/
{
    char f_dum1[GPIV_MAX_CHARS] = "";
    char f_dum2[GPIV_MAX_CHARS] = "";
    char fname_out[GPIV_MAX_CHARS] = "";

    int j, k, l;
     guint16 **img_in1 = NULL, **img_in2 = NULL;
     int first_file = image_proc_par.first_file;
     int last_file = image_proc_par.last_file;

     int nrows = image_par.nrows;
     int ncolumns = image_par.ncolumns;
     char fname_in[GPIV_MAX_CHARS];

     if (print_par) printf("# Subtracting mean from input data\n");
     img_in1 = gpiv_alloc_img(image_par);
     img_in2 = gpiv_alloc_img(image_par);
     for (j = first_file; j <= last_file; j++) {
	 memset(img_in1[0], 0, (sizeof(guint16)) * nrows * ncolumns);
	 memset(img_in2[0], 0, (sizeof(guint16)) * nrows * ncolumns);
         if (fname_logic == TRUE) {
             make_fname_in(image_proc_par, fname_base, j, fname_in);
             gpiv_fread_image(fname_in, img_in1, img_in2, image_par);
         } else {
             gpiv_read_image(img_in1, img_in2, image_par);
         }
         
	 for (k = 0; k < nrows; k++) {
	      for (l = 0; l < ncolumns; l++) {
		   img_in1[k][l] -= img_subtr[k][l];
		   img_in2[k][l] -= img_subtr[k][l];
	      }
	 }

         if (fname_logic == TRUE) {
             snprintf(f_dum1, GPIV_MAX_CHARS, "%s%d", fname_base, j);
             gpiv_io_make_fname(f_dum1, GPIV_EXT_TA, f_dum2);
             gpiv_io_make_fname(f_dum2, GPIV_EXT_RAW_IMAGE, fname_out);
             gpiv_fwrite_image(fname_out, img_in1, img_in2, image_par);
             if (print_par) printf("# Output file: %s\n", fname_out);
             snprintf(f_dum1, GPIV_MAX_CHARS, " ");;
             snprintf(f_dum2, GPIV_MAX_CHARS, " ");;
             snprintf(fname_out, GPIV_MAX_CHARS, " ");;
         } else {
             gpiv_write_image(img_in1, img_in2, image_par);
         }

    }
    
    gpiv_free_img(img_in1, image_par);
    gpiv_free_img(img_in2, image_par);

}



int 
main(int argc, 
     char *argv[]
     )
/* ----------------------------------------------------------------------------
 * main routine to calculates time-averaged image intensity
 */
{
    FILE *fp_par_dat;
    char fname_base[GPIV_MAX_CHARS], 
        fname_header[GPIV_MAX_CHARS],
        fname_out[GPIV_MAX_CHARS],
	fname_parameter[GPIV_MAX_CHARS];
    guint16 **img_mean = NULL;
    GpivImageProcPar image_proc_par;
    GpivImagePar image_par;


/*
 * Image header info initialization
 */
    image_par.ncolumns_logic = FALSE;
    image_par.nrows_logic = FALSE;
    image_par.depth_logic = FALSE;
    image_par.x_corr_logic = FALSE;

/*
 * Image processing parameter initialization
 */
    image_proc_par.first_dir = 0;
    image_proc_par.last_dir = 0;
    image_proc_par.dir_prefix = 0;    
    image_proc_par.first_file = 0;
    image_proc_par.last_file = 0;
    image_proc_par.file_prefix = 0;
    image_proc_par.operator = 0;

    image_proc_par.first_dir_logic = FALSE;
    image_proc_par.last_dir_logic = FALSE;
    image_proc_par.dir_prefix_logic = FALSE;     
    image_proc_par.first_file_logic = FALSE;
    image_proc_par.last_file_logic = FALSE;
    image_proc_par.file_prefix_logic = FALSE;
    image_proc_par.operator_logic = FALSE;

    command_args(argc, argv, fname_base, &image_par, &image_proc_par);
    if (print_par) {
        printf("# %s\n# Command line options:\n", RCSID);
        gpiv_img_print_header(image_par);
        gpiv_imgproc_print_parameters(image_proc_par);
    }


    if (fname_logic == TRUE) {
/* 
 * Generating proper filenames 
 */
	make_fname_out(fname_base, fname_header, fname_parameter, fname_out);
	if (print_par) 
            printf("\n# Parameters written to: %s", fname_parameter);


/* 
 * Prints command line parameters to par-file 
 */
	if ((fp_par_dat = fopen(fname_parameter, "a")) == NULL) {
	    gpiv_error("%s: failure opening %s for input",
                       RCSID, fname_parameter);
	}
	fprintf(fp_par_dat, "\n\n# %s\n# Command line options:\n", RCSID);
	gpiv_imgproc_fprint_parameters(fp_par_dat, image_proc_par);


/*
 * Reads image header data from file.h and/or resource files 
 * if not overridden by the commandline options
 */
        gpiv_scan_parameter( "", fname_header, &image_par, print_par);
        gpiv_scan_resourcefiles(GPIV_IMAGE_PAR_KEY, &image_par, print_par);
        gpiv_img_fprint_header(fp_par_dat, image_par);


/*
 *  Reading parametes from PARFILE (and writing to data par-file)
 */
        gpiv_scan_parameter(GPIV_IMAGEPROC_PAR_KEY, PARFILE, &image_proc_par, 
			    print_par);
        gpiv_scan_resourcefiles(GPIV_IMAGEPROC_PAR_KEY, &image_proc_par, 
                                print_par);
        gpiv_imgproc_fprint_parameters(fp_par_dat, image_proc_par);
	fclose(fp_par_dat);


    } else {
/*
 * use stdin, stdout
 */
        gpiv_scan_resourcefiles("", &image_par, print_par);
        gpiv_scan_parameter(GPIV_IMAGEPROC_PAR_KEY, PARFILE, &image_proc_par, 
			    print_par);
        gpiv_scan_resourcefiles(GPIV_IMAGEPROC_PAR_KEY, &image_proc_par, 
                                print_par);
    }


/*
 * Check parameters on correct values and adjust belonging  variables
 */



/* 
 * Now the parameters are known, data  memory can be allocated
 */
    img_mean = gpiv_alloc_img(image_par);
 /*    imgsdev = gpiv_alloc_img(image_par); */


/*
 * Here the function calls of the post-processing; flipping x and y data etc.
 * Inclusive reading of the PIV data its own functions.
 */

    gpiv_img_mean(image_par, image_proc_par, fname_base, img_mean
                  /* , img_sdev */
                  );

    if (image_proc_par.operator == GPIV_SUBTRACT)
        gpiv_img_subtract(image_par, image_proc_par, fname_base, img_mean);


/* 
 * And writing to output 
 */
    if (fname_logic == TRUE) {
	 gpiv_fwrite_image(fname_out, img_mean, img_mean, image_par);
    } else {
	 gpiv_write_image(img_mean, img_mean, image_par);
    }


/*
 * Freeing allocated memory of matrices
 */
    gpiv_free_img(img_mean, image_par);
/*     gpiv_free_img(img_sdev, image_par); */

    return 0;
}