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

/*------------------------------------------------------------------------------
 
  rr - rr calculates the mean particle displacement for (Digital) Particle
       Image Velocimetry (DPIV) by means of FFT. 

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

   This program 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.  

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

#include <stdlib.h>
#include <fftw3.h>
#include <gpiv.h>
#undef USE_MTRACE
#ifdef USE_MTRACE
#include <mcheck.h>
#endif

#define PARFILE "gpivrc"	/* Parameter file name */

/* static float sqrarg; */
/* #define SQR(a) ((sqrarg=(a)) == 0.0 ? 0.0 : sqrarg*sqrarg) */
/* #define SWAP(a,b) tempr=(a);(a)=(b);(b)=tempr */

#define USAGE "\
Usage: gpiv_rr [--c_dif] [--no_cdif] [--cf int] [--cl int] [--cp int] \n\
          [-g] [--gauss][-h | --help] [--ia_size_i int] [--ia_size_f int] \n\
          [--ia_shift int] [--ischeme int] [--ifit 0/1/2/3] [--linec int int int] \n\
          [--liner int int int] [-p | --print] [--peak 1/2/3] [--p_piv] \n\
          [--point int int] [--rf int] [--rl int] [--rp int] [--spof] \n\
          [-s float] [-v | --version] [--val_r int] [--val_s] [--val_t float] \n\
          [filename] < stdin > stdout \n\
\n\
keys: \n\
--cf COL:              first column COL of area to interrogate \n\
--cl COL:              last column COL of area to interrogate \n\
--cp COLUMNS:          pre-shift in x-direction (COLUMNS) \n\
-g:                    graphical visualization with gnuplot (needs -f) \n\
--gauss:               Gauss weighting of interrogation area \n\
-h | --help:           this on-line help \n\
--ia_size_i SIZE:      size of initial interrogation area (>= isi1) \n\
--ia_size_f SIZE:      size of final interrogation area \n\
--ia_shift SHIFT:      shift of adjacent interrogation area \n\
--ischeme 0/1/2/3/4:   interrogation scheme: no correction (0), linear \n\
                       weighting (1), zero offset (2), zero offset with central \n\
                       differential (3), image deformation (4) \n\
--ifit 0/1/2/3:        sub-pixel interpolation type: none (0), Gauss (1), \n\
                       Parabolic (2) or Centre of Gravity (3). \n\
--linec C RF RL:       selects a vertical line at column R to interrogate \n\
                       from row RF to row RL \n\
--liner R CF CL:       selects an horizontal line at row R to interrogate \n\
                       from column CF to column CL \n\
-p | --print:          print parameters and other info to stdout \n\
--peak #:              find maximum of #-th covariance peak \n\
--p_piv:               prints piv results to stdout, even if -f has been used \n\
--point X Y:           select one single point in the image to interrogate \n\
--rf ROW:              first ROW of area to interrogate \n\
--rl ROW:              last ROW of area to interrogate \n\
--rp ROWS:             pre-shift in y-direction (ROWS) \n\
--spof:                symmetric phase only filtering \n\
-s S:                  scale factor for graphic output with gnuplot \n\
-v | --version:        version number \n\
--val_r N:             validation; residu type calculated from: snr (0), median \n\
                       (1) or normalized median (2) \n\
--val_s N:             validation; substitution of erroneous vector by: nothing \n\
                       (0), local mean from the surroundings (1), \n\
                       the median of the surroundings (2), next highest \n\
                       correlation peak (3) (needs -f) \n\
--val_t F:             validation; threshhold of residus to be accepted \n\
filename:              full image filename. Overrides stdin and stdout.  \n\
                       If stdin and stdout are used, input is expected \n\
                       to be a .png formatted image \n\
"

#ifdef DEBUG
#define USAGE_DEBUG "\
Developers version also contains: \n\
rr:	  [-p_main int] \n\
\n\
keys: \n\
-p_'function' int:     prints data to be generated in the function; the \n\
                       higher the number, the more detailed the output.  \n\
                       For int = 10, rr will exit at the end of the function \n\
"
#endif

#define HELP "\
rr interrogates an image (pair) in order to obtain displacements of \n\
particles for (Digital) Particle Image Velocimetry (PIV). \n\
"
#define RCSID "$Id: rr.c,v 3.26 2008-09-25 13:18:48 gerber Exp $"

/*
 * Global parameters and variables
 */
gboolean use_stdin_stdout = FALSE;
gboolean verbose = FALSE;


#define GNUPLOT_DISPLAY_SIZE 500
#define GNUPLOT_DISPLAY_COLOR "LightBlue"

int gnuplot = 0;
float gnuplot_scale = 1.0;
gboolean gnuplot__set = FALSE, gnuplot_scale__set = FALSE;



static void 
command_args (int argc, 
              char *argv[], 
              char fname[GPIV_MAX_CHARS],
              GpivPivPar *piv_par,
              GpivValidPar *valid_par
              )
/*-----------------------------------------------------------------------------
 * Command line argument handling
 */
{
    char c;
    int argc_next;


    /*
     * Processing of command line arguments
     */
    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) {

                /*
                 * graphic output with gnuplot
                 */
            case 'g':	       
                gnuplot = 1;
                gnuplot__set = TRUE;
                /* 		 argc_next = 1; */
                break;
                
            case 'h':
                printf ("%s\n", argv[0]);
                printf ("%s\n", HELP);
                printf ("%s", USAGE);
#ifdef DEBUG
                printf ("\n%s", USAGE_DEBUG);
#endif
                exit (0);

                /*
                 * print paramaters to stdout
                 */
            case 'p':
                verbose = TRUE;
                break;
                
                /*
                 * scaling for graphic output with gnuplot
                 */
            case 's':	       
                gnuplot_scale = atof (*++argv);
                gnuplot_scale__set = TRUE;
                --argc;
                argc_next = 1;
                break;
                                
                /*
                 * Print version and exits
                 * Use Revision Control System (RCS) for version
                 */
            case 'v':
                printf ("%s\n", RCSID);
                exit (0);
                
                
                /*
                 * long option keys
                 */
            case '-':
                if (strcmp ("-help", *argv) == 0) {
                    printf ("\n%s", argv[0]);
                    printf ("\n%s", HELP);
                    printf ("\n%s", USAGE);
                    exit (0);
                } else if (strcmp ("-print", *argv) == 0) {
                    verbose = TRUE;
                } else if (strcmp ("-version", *argv) == 0) {
                    printf ("%s\n", RCSID);
                    exit (0);

                    /*
                     * Rr.Col_start
                     */
                } else if (strcmp ("-cf", *argv) == 0) {
                    piv_par->col_start = atoi (*++argv);
                    piv_par->col_start__set = TRUE;
                    piv_par->int_geo = GPIV_AOI;
                    piv_par->int_geo__set = TRUE;
                    --argc;
                    argc_next = 1;

                    /*
                     * overrides Rr.Col_end
                     */
                } else if (strcmp ("-cl", *argv) == 0) {
                    piv_par->col_end = atoi (*++argv);
                    piv_par->col_end__set = TRUE;
                    piv_par->int_geo = GPIV_AOI;
                    piv_par->int_geo__set = TRUE;
                    --argc;
                    argc_next = 1;

                    /*
                     * pre-shift c (col)
                     */
                } else if (strcmp ("-cp", *argv) == 0) {
                    piv_par->pre_shift_col = atoi (*++argv);
                    piv_par->pre_shift_col__set = TRUE;
                    --argc;
                    argc_next = 1;

                    /*
                     * interrogation area weighting with gauss function
                     */
                } else if (strcmp ("-gauss", *argv) == 0) {
                    piv_par->gauss_weight_ia = TRUE;
                    piv_par->gauss_weight_ia__set = TRUE;

                    /*
                     * overrides Rr.int_size_f
                     */
                } else if (strcmp ("-ia_size_f", *argv) == 0) {
                    piv_par->int_size_f = atoi (*++argv);
                    piv_par->int_size_f__set = TRUE;
                    --argc;
                    argc_next = 1;
                                        
                    /*
                     * overrides Rr.int_size_i
                     */
                } else if (strcmp ("-ia_size_i", *argv) == 0) {
                    piv_par->int_size_i = atoi (*++argv);
                    piv_par->int_size_i__set = TRUE;
                    --argc;
                    argc_next = 1;

                    /*
                     * overrides Rr.int_shift
                     */
                } else if (strcmp ("-ia_shift", *argv) == 0) {
                    piv_par->int_shift = atoi (*++argv);
                    piv_par->int_shift__set = TRUE;
                    --argc;
                    argc_next = 1;
                    
                    /*
                     * overrides Rr.Int_scheme
                     */
                } else if (strcmp ("-ischeme", *argv) == 0) {
                    piv_par->int_scheme = atoi (*++argv);
                    if (piv_par->int_scheme != 0 
                        && piv_par->int_scheme != 1 
                        && piv_par->int_scheme != 2
                        && piv_par->int_scheme != 3 
                        && piv_par->int_scheme != 4) {
                        gpiv_error ("%s: invalid value of Int_scheme (1, 2, 3 or 4)", 
                                    argv[0]);
                    }
                    piv_par->int_scheme__set = TRUE;
                    --argc;
                    argc_next = 1;

                    /*
                     * overrides Rr.Ifit
                     */
                } else if (strcmp ("-ifit", *argv) == 0) {
                    piv_par->ifit = atoi (*++argv);
                    if (piv_par->ifit != -1 && piv_par->ifit != 0 && 
                        piv_par->ifit != 1 && piv_par->ifit != 2
                        && piv_par->ifit != 3) {
                        gpiv_error ("%s: invalid value of Ifit (1, 2 or 3)", 
                                    argv[0]);
                    }
                    piv_par->ifit__set = TRUE;
                    --argc;
                    argc_next = 1;

                    /*
                     * define interrogate at a line
                     */
                } else if ((strcmp ("-linec", *argv) == 0)) {
                    piv_par->int_line_col = atoi (*++argv);
                    piv_par->int_line_row_start = atoi (*++argv);
                    piv_par->int_line_row_end =  atoi (*++argv);
                    piv_par->int_geo = GPIV_LINE_C;
                    piv_par->int_line_col__set = TRUE;
                    piv_par->int_line_row_start__set = TRUE;
                    piv_par->int_line_row_end__set = TRUE;
                    if (piv_par->int_geo__set == TRUE) {
                        gpiv_error ("%s: only a horizontal line"
                                    " may be defined\n"
                                    " not a vertical, a point or an area of"
                                    " interest as well\n",
                                    argv[0]);
                    } else {
                        piv_par->int_geo__set = TRUE;
                    }
                    --argc;
                    --argc;
                    --argc;
                    argc_next = 1;
                } else if ((strcmp ("-liner", *argv) == 0)) {
                    piv_par->int_line_row = atoi (*++argv);
                    piv_par->int_line_col_start = atoi (*++argv);
                    piv_par->int_line_col_end =  atoi (*++argv);
                    piv_par->int_geo = GPIV_LINE_R;
                    piv_par->int_line_row__set = TRUE;
                    piv_par->int_line_col_start__set = TRUE;
                    piv_par->int_line_col_end__set = TRUE;
                    if (piv_par->int_geo__set == TRUE) {
                        gpiv_error ("%s: a vertical horizontal line"
                                    " may be defined\n"
                                    " not an horizontal, a point or an area of"
                                    " interest as well\n",
                                    argv[0]);
                    } else {
                        piv_par->int_geo__set = TRUE;
                    }
                    --argc;
                    --argc;
                    --argc;
                    argc_next = 1;
                    
                    /*
                     * cov peak # to fit
                     */
                } else if (strcmp ("-peak", *argv) == 0) {
                    piv_par->peak = atoi (*++argv);
                    --argc;
                    piv_par->peak__set = TRUE;
                    argc_next = 1;
                    
                    /*
                     * # piv values to stdout
                     */
                } else if (strcmp ("-p_piv", *argv) == 0) {
                    piv_par->print_piv = 1;
                    piv_par->print_piv__set = TRUE;
                    argc_next = 1;
                    
                    /*
                     * interrogation at point(row, col)
                     */
                } else if (strcmp ("-point", *argv) == 0) {
                    piv_par->int_point_col = atoi (*++argv);
                    piv_par->int_point_row = atoi (*++argv);
                    piv_par->int_point_col__set = TRUE;
                    piv_par->int_point_row__set = TRUE;
                    piv_par->int_geo = GPIV_POINT;
                    if (piv_par->int_geo__set == TRUE) {
                        gpiv_error ("%s: only a point may be defined\n"
                                    " not an horizontal, a vertical line"
                                    " or an area of interest as well\n",
                                    argv[0]);
                    } else {
                        piv_par->int_geo__set = TRUE;
                    }
                    --argc;
                    --argc;
                    argc_next = 1;
                    
                    /*
                     * overrides Rr.Row_start
                     */
                } else if (strcmp ("-rf", *argv) == 0) {
                    piv_par->row_start = atoi (*++argv);
                    piv_par->row_start__set = TRUE;
                    piv_par->int_geo = GPIV_AOI;
                    piv_par->int_geo__set = TRUE;
                    --argc;
                    argc_next = 1;

                    /*
                     * overrides Rr.Row_end
                     */
		} else if (strcmp ("-rl", *argv) == 0) {
		    piv_par->row_end = atoi (*++argv);
		    piv_par->row_end__set = TRUE;
                    piv_par->int_geo = GPIV_AOI;
                    piv_par->int_geo__set = TRUE;
		    --argc;
		    argc_next = 1;
                    
                    /*
                     * pre-shift r (row)
                     */
		} else if (strcmp ("-rp", *argv) == 0) {
		    piv_par->pre_shift_row = atoi (*++argv);
		    piv_par->pre_shift_row__set = TRUE;
		    --argc;
		    argc_next = 1;

                    /*
                     * Symmetric phase only filtering
                     */
		} else if (strcmp ("-spof", *argv) == 0) {
                    piv_par->spof_filter = TRUE;
                    piv_par->spof_filter__set = TRUE;

                    /*
                     * validation parameter: residu type
                     */
                } else if (strcmp ("-val_r", *argv) == 0) {
                    valid_par->residu_type = atoi (*++argv);
                    valid_par->residu_type__set = TRUE;
                    --argc;
                    argc_next = 1;
                    if (valid_par->residu_type != 0 
                        && valid_par->residu_type != 1 
                        && valid_par->residu_type != 2) {
                        gpiv_error ("%s: invalid value of VALID.Residu_type (0, 1, or 2)", 
                                    argv[0]);
                    }
                    
                    /*
                     * validation parameter: substitution type
                     */
                } else if (strcmp ("-val_s", *argv) == 0) {
                    valid_par->subst_type = atoi (*++argv);
                    valid_par->subst_type__set = TRUE;
                    --argc;
                    argc_next = 1;
                    if (valid_par->subst_type != 0 
                        && valid_par->subst_type != 1 
                        && valid_par->subst_type != 2
                        && valid_par->subst_type != 3) {
                        gpiv_error ("%s: invalid value of VALID.Subst_type (0, 1, 2 or 3)", 
                                    argv[0]);
                    }

                    /*
                     * validation parameter: threshold value
                     */
                } else if (strcmp ("-val_t", *argv) == 0) {
                    valid_par->residu_max = atof (*++argv);
                    valid_par->residu_max__set = TRUE;
                    --argc;
                    argc_next = 1;

                } else {
		    gpiv_error ("%s: unknown option: %s", argv[0], *argv);
		}
		argc_next = 1;
		break;

	    default:
		gpiv_error ("%s: unknown option: %s", argv[0], *argv);
		break;

	    }
    }


    /*
     * Check if filename or stdin /stdout is used
     */
    if (argc == 1) {
        use_stdin_stdout = FALSE;
        strcpy(fname, argv[argc - 1]);
    } else if (argc == 0) {
        use_stdin_stdout = TRUE;
        verbose = FALSE;
        if (gnuplot == 1) {
            gpiv_error ("%s: an image filename for input has to be used "
                        "in combination with 'gnuplot'", 
                        argv[0]);
        }
    } else {
        gpiv_error("\n%s: unknown argument: %s: %s", argv[0], *argv, USAGE);
    }



}



static gchar *
make_fname (FILE * fp_par_out, 
            char *fname_in, 
            char *fname_out,
            GpivPivPar *piv_par
            )
/*-----------------------------------------------------------------------------
 * generates filenames
 */
{
    gchar *err_msg = NULL;
    gchar *fname_base = NULL;

    if (fname_in == NULL ) {
        err_msg = "make_fname: \"fname_in == NULL\"";
        return (err_msg);
    }
 
    /*
     * Stripping filename
     */
    fname_base = g_strdup(fname_in);
    strtok(fname_base, ".");

    /*
     * filenames for output PIV data
     */
    if (fname_out != NULL ) {
        gpiv_io_make_fname (fname_base, GPIV_EXT_PIV, fname_out);
        if (piv_par->print_piv == 0) {
           if (verbose)
                printf ("\n# outputfile is: %s\n", fname_out);
            fprintf (fp_par_out, "\n# outputfile is: %s\n", fname_out);
        }
    }

    g_free (fname_base);
    return (err_msg);
}


#ifdef ENABLE_MPI

static GpivImage *
alloc_image_struct(void)
/*-----------------------------------------------------------------------------
 * Allocates GpivImage structure
 */
{
    GpivImage *image = g_new0 (GpivImage, 1);
    GpivImagePar *image_par = g_new0 (GpivImagePar, 1);;


    image->header = image_par;


    return image;
}


static GpivPivPar *
alloc_pivpar_struct(void)
/*-----------------------------------------------------------------------------
 * Allocates GpivPivPar structure
 */
{
    GpivPivPar *piv_par = g_new0 (GpivPivPar, 1);
    return piv_par;
}


static GpivValidPar *
alloc_validpar_struct(void)
/*-----------------------------------------------------------------------------
 * Allocates GpivValidPar structure
 */
{
    GpivValidPar *valid_par = g_new0 (GpivValidPar, 1);
    return valid_par;
}

#endif /* ENABLE_MPI */


int 
main (int argc, 
      char *argv[]
      )
/*-----------------------------------------------------------------------------
 * Main routine of rr to calculate estimators of particle image displacements
 * from PIV images 
 */
{
    char *err_msg = NULL;
    FILE *fp_out = NULL, *fp_par_out = NULL;
    char fname_in[GPIV_MAX_CHARS],
        fname_out[GPIV_MAX_CHARS],
        fname_par[GPIV_MAX_CHARS];

    GpivImage *image = NULL;
    GpivPivData *piv_data = NULL;
    GpivPivPar *piv_par = g_new0 (GpivPivPar, 1);
    GpivValidPar *valid_par = g_new0 (GpivValidPar, 1);
    
#ifdef _OPENMP
    int np, nprocs, num_threads;
#endif /* _OPENMP */

#ifdef ENABLE_MPI
    /*     gint argc_mpi; */
    /*     gchar *argv_mpi[10]; */
    int nprocs, rank = 0, namelen, root = 0;
    char processor_name[MPI_MAX_PROCESSOR_NAME];
#endif /* ENABLE_MPI */
        
#ifdef ENABLE_MPI
    /* 
     * Initializing MPI
     */
    MPI_Init(&argc, &argv);
    MPI_Comm_rank(MPI_COMM_WORLD, &rank);
    MPI_Comm_size(MPI_COMM_WORLD, &nprocs);
    MPI_Get_processor_name(processor_name, &namelen);

/*     gpiv_warning("main:: 0 rank=%d nprocs=%d on %s", rank, nprocs, processor_name); */
    /*     argc_mpi = argc - 4; */
    /*     for (i = 0; i<argc_mpi; i++) { */
    /*         gpiv_warning("main:: rank=%d argv[%d]=%s", rank, i, argv[i]); */
    /*         argv_mpi[i] = argv[i]; */
    /*     } */
#endif /* ENABLE_MPI */

    /* 
     * Initializing parameters
     */
    gpiv_piv_parameters_set (piv_par, FALSE);
    gpiv_valid_parameters_set (valid_par, FALSE);

#ifdef _OPENMP
    /* 
     * Initializing OMP
     */
    np = omp_get_num_procs();
    omp_set_num_threads(np);
    num_threads = omp_get_max_threads();
#endif /* _OPENMP */

#ifdef ENABLE_MPI
    if (rank == 0) {
#endif /* ENABLE_MPI */

        command_args (argc, 
                      argv, 
                      fname_in, 
                      piv_par,
                      valid_par);

        if (verbose) {
            printf ("\n# %s\n# Command line options:\n", RCSID);
            gpiv_piv_print_parameters (stdout, piv_par);
            gpiv_valid_print_parameters (stdout, valid_par);
        }

        /*
         * Reads piv and valid parameters from PARFILE and/or resource files 
         * if not overridden by the commandline options.
         * Check if all parameters have been read. Else, set as libgpiv defaults
         */
        gpiv_scan_parameter (GPIV_PIVPAR_KEY, PARFILE, piv_par, verbose);
        gpiv_scan_resourcefiles (GPIV_PIVPAR_KEY, piv_par, verbose);
        if ((err_msg = 
             gpiv_piv_check_parameters_read (piv_par, NULL))
            != NULL) gpiv_warning ("%s: %s", argv[0], err_msg);

        gpiv_scan_parameter (GPIV_VALIDPAR_KEY, PARFILE, valid_par, verbose);
        gpiv_scan_resourcefiles (GPIV_VALIDPAR_KEY, valid_par, verbose);
        if ((err_msg = 
             gpiv_valid_check_parameters_read (valid_par, NULL))
            != NULL) gpiv_warning ("%s: %s", argv[0], err_msg);

       /*
         * Creates file names if not stdin / stdout are used
         * and save parameters
         */
        if (use_stdin_stdout == FALSE) {
            gchar *fname_base = g_strdup (fname_in);
            strtok(fname_base, ".");
            gpiv_io_make_fname (fname_base, GPIV_EXT_PAR, fname_par);
            g_free (fname_base);

            if ((fp_par_out = fopen (fname_par, "w")) == NULL) {
                gpiv_error ("%s error: failure opening %s for input",
                            argv[0], fname_par);
            }
            if ((err_msg = 
                 make_fname (fp_par_out, fname_in, fname_out, piv_par))
                != NULL) {
                gpiv_error ("%s: %s\n", argv[0], err_msg);
            }
       

            /*
             * Write parameters to parameterfile
             */
            gpiv_piv_print_parameters (fp_par_out, piv_par);
            gpiv_valid_print_parameters (fp_par_out, valid_par);
            fclose (fp_par_out);
        }


        /*
         * reads the image
         */

#ifdef USE_MTRACE
        /*     fprintf(stderr, "Calling interrogate with mtrace()\n"); */
        mtrace();
#endif
        if (use_stdin_stdout) {
#ifdef DEBUG
            gpiv_warning ("gpiv_rr:: using stdin");
#endif
            if ((image = gpiv_read_image (stdin)) == NULL) {
                gpiv_error ("%s: gpiv_fread_image_ni", argv[0]);
            }
        } else {
#ifdef DEBUG
            gpiv_warning ("rr:: using fp, fname = %s", fname_in);
#endif
            if ((image = gpiv_fread_image (fname_in)) == NULL) {
                gpiv_error ("%s: gpiv_fread_image", argv[0]);
            }
        }

#ifdef ENABLE_MPI
    } else {
        /*
         * allocates the memory of structures for threads of which rank != 0
         */
        image = alloc_image_struct();
        piv_par = alloc_pivpar_struct();
        valid_par = alloc_validpar_struct();

    }


    if (MPI_Bcast(&verbose, 1, MPI_INT, root, MPI_COMM_WORLD) 
        != MPI_SUCCESS) {
        gpiv_error ("gpiv_rr: An error ocurred when calling MPI_Bcast");
    }
    gpiv_img_mpi_bcast_image (image, TRUE);
    if ((err_msg = gpiv_piv_mpi_bcast_pivpar (piv_par)) != NULL) {
        gpiv_error (err_msg);
    }
    gpiv_valid_mpi_bcast_validpar (valid_par);

#endif /* ENABLE_MPI */
#ifdef DEBUG
    gpiv_warning ("gpiv_rr:: calling gpiv_piv_interrogate_img");
#endif


/*
 * Interrogating image
 */
    if ((piv_data = gpiv_piv_interrogate_img (image, piv_par, valid_par, verbose))
        == NULL) {
        gpiv_error ("%s: gpiv_piv_interrogate_img failed", argv[0]);
    }


#ifdef ENABLE_MPI
    if (rank == 0) {
#endif
        /*
         * Displays PIV data
         */
        if (gnuplot) gpiv_piv_gnuplot (fname_in, gnuplot_scale, 
                                       GNUPLOT_DISPLAY_COLOR, GNUPLOT_DISPLAY_SIZE,
                                       image->header, piv_par, piv_data);

        /*
         * Add comments to PIV data and saving
         */
        piv_data->comment = g_strdup_printf ("# Software: %s\n", RCSID);
        piv_data->comment = gpiv_add_datetime_to_comment (piv_data->comment);
        piv_data->comment = g_strconcat (piv_data->comment, 
                                         "# Data type: Particle Image Velocities\n", 
                                         NULL);


        if ((piv_par->print_piv == TRUE && use_stdin_stdout == FALSE)
            || use_stdin_stdout == TRUE) {
            fp_out = stdout;
            if ((err_msg = 
                 gpiv_write_pivdata (fp_out, piv_data, FALSE))
                != NULL) gpiv_error ("%s: %s", argv[0], err_msg);

        } else {
            if ((fp_out = fopen (fname_out, "wb")) == NULL) {
                gpiv_error ("%s: failing fopen %s", argv[0], fname_out);
            }
            if ((err_msg = 
                 gpiv_write_pivdata (fp_out, piv_data, FALSE))
                != NULL) gpiv_error ("%s: %s", argv[0], err_msg);
            fclose (fp_out);
        
        }
#ifdef ENABLE_MPI
    }
#endif

    /*
     * Free memory of images and data
     */
    gpiv_free_pivdata (piv_data);
    if (image->frame1 != NULL && image->frame2 != NULL) {
        gpiv_free_img (image);
    }


#ifdef USE_MTRACE
    muntrace();
#endif

/* #endif */ /* __DISABLE */
#ifdef ENABLE_MPI
#ifdef DEBUG
        gpiv_warning("main:: rank=%d: Calling MPI_Finalize()", rank);
#endif /* DEBUG */
        MPI_Finalize();
#endif /* ENABLE+MPI */
        exit (0);
}