src/evaluate/rr.c

   1 /* -*- Mode: C; indent-tabs-mode: nil; c-basic-offset: 4 c-style: "K&R" -*- */
   2
   3 /*------------------------------------------------------------------------------
   4
   5   rr - rr calculates the mean particle displacement for (Digital) Particle
   6        Image Velocimetry (DPIV) by means of FFT.
   7
   8    Copyright (C) 2002, 2003, 2004, 2005, 2006, 2007, 2008
   9    Gerber van der Graaf <gerber_graaf@users.sourceforge.net
  10
  11    This program is free software; you can redistribute it and/or modify
  12    it under the terms of the GNU General Public License as published by
  13    the Free Software Foundation; either version 2, or (at your option)
  14    any later version.
  15
  16    This program is distributed in the hope that it will be useful,
  17    but WITHOUT ANY WARRANTY; without even the implied warranty of
  18    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  19    GNU General Public License for more details.
  20
  21    You should have received a copy of the GNU General Public License
  22    along with this program; if not, write to the Free Software Foundation,
  23    Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
  24
  25
  26  --------------------------------------------------------------------------- */
  27
  28 #include <stdlib.h>
  29 #include <fftw3.h>
  30 #include <gpiv.h>
  31 #undef USE_MTRACE
  32 #ifdef USE_MTRACE
  33 #include <mcheck.h>
  34 #endif
  35
  36 #define PARFILE "gpivrc"        /* Parameter file name */
  37
  38 /* static float sqrarg; */
  39 /* #define SQR(a) ((sqrarg=(a)) == 0.0 ? 0.0 : sqrarg*sqrarg) */
  40 /* #define SWAP(a,b) tempr=(a);(a)=(b);(b)=tempr */
  41
  42 #define USAGE "\
  43 Usage: gpiv_rr [--c_dif] [--no_cdif] [--cf int] [--cl int] [--cp int] \n\
  44           [-g] [--gauss][-h | --help] [--ia_size_i int] [--ia_size_f int] \n\
  45           [--ia_shift int] [--ischeme int] [--ifit 0/1/2/3] [--linec int int int] \n\
  46           [--liner int int int] [-p | --print] [--peak 1/2/3] [--p_piv] \n\
  47           [--point int int] [--rf int] [--rl int] [--rp int] [--spof] \n\
  48           [-s float] [-v | --version] [--val_r int] [--val_s] [--val_t float] \n\
  49           [filename] < stdin > stdout \n\
  50 \n\
  51 keys: \n\
  52 --cf COL:              first column COL of area to interrogate \n\
  53 --cl COL:              last column COL of area to interrogate \n\
  54 --cp COLUMNS:          pre-shift in x-direction (COLUMNS) \n\
  55 -g:                    graphical visualization with gnuplot (needs -f) \n\
  56 --gauss:               Gauss weighting of interrogation area \n\
  57 -h | --help:           this on-line help \n\
  58 --ia_size_i SIZE:      size of initial interrogation area (>= isi1) \n\
  59 --ia_size_f SIZE:      size of final interrogation area \n\
  60 --ia_shift SHIFT:      shift of adjacent interrogation area \n\
  61 --ischeme 0/1/2/3/4:   interrogation scheme: no correction (0), linear \n\
  62                        weighting (1), zero offset (2), zero offset with central \n\
  63                        differential (3), image deformation (4) \n\
  64 --ifit 0/1/2/3:        sub-pixel interpolation type: none (0), Gauss (1), \n\
  65                        Parabolic (2) or Centre of Gravity (3). \n\
  66 --linec C RF RL:       selects a vertical line at column R to interrogate \n\
  67                        from row RF to row RL \n\
  68 --liner R CF CL:       selects an horizontal line at row R to interrogate \n\
  69                        from column CF to column CL \n\
  70 -p | --print:          print parameters and other info to stdout \n\
  71 --peak #:              find maximum of #-th covariance peak \n\
  72 --p_piv:               prints piv results to stdout, even if -f has been used \n\
  73 --point X Y:           select one single point in the image to interrogate \n\
  74 --rf ROW:              first ROW of area to interrogate \n\
  75 --rl ROW:              last ROW of area to interrogate \n\
  76 --rp ROWS:             pre-shift in y-direction (ROWS) \n\
  77 --spof:                symmetric phase only filtering \n\
  78 -s S:                  scale factor for graphic output with gnuplot \n\
  79 -v | --version:        version number \n\
  80 --val_r N:             validation; residu type calculated from: snr (0), median \n\
  81                        (1) or normalized median (2) \n\
  82 --val_s N:             validation; substitution of erroneous vector by: nothing \n\
  83                        (0), local mean from the surroundings (1), \n\
  84                        the median of the surroundings (2), next highest \n\
  85                        correlation peak (3) (needs -f) \n\
  86 --val_t F:             validation; threshhold of residus to be accepted \n\
  87 filename:              full image filename. Overrides stdin and stdout.  \n\
  88                        If stdin and stdout are used, input is expected \n\
  89                        to be a .png formatted image \n\
  90 "
  91
  92 #ifdef DEBUG
  93 #define USAGE_DEBUG "\
  94 Developers version also contains: \n\
  95 rr:       [-p_main int] \n\
  96 \n\
  97 keys: \n\
  98 -p_'function' int:     prints data to be generated in the function; the \n\
  99                        higher the number, the more detailed the output.  \n\
 100                        For int = 10, rr will exit at the end of the function \n\
 101 "
 102 #endif
 103
 104 #define HELP "\
 105 rr interrogates an image (pair) in order to obtain displacements of \n\
 106 particles for (Digital) Particle Image Velocimetry (PIV). \n\
 107 "
 108 #define RCSID "$Id: rr.c,v 3.26 2008-09-25 13:18:48 gerber Exp $"
 109
 110 /*
 111  * Global parameters and variables
 112  */
 113 gboolean use_stdin_stdout = FALSE;
 114 gboolean verbose = FALSE;
 115
 116
 117 #define GNUPLOT_DISPLAY_SIZE 500
 118 #define GNUPLOT_DISPLAY_COLOR "LightBlue"
 119
 120 int gnuplot = 0;
 121 float gnuplot_scale = 1.0;
 122 gboolean gnuplot__set = FALSE, gnuplot_scale__set = FALSE;
 123
 124
 125
 126 static void
 127 command_args (int argc,
 128               char *argv[],
 129               char fname[GPIV_MAX_CHARS],
 130               GpivPivPar *piv_par,
 131               GpivValidPar *valid_par
 132               )
 133 /*-----------------------------------------------------------------------------
 134  * Command line argument handling
 135  */
 136 {
 137     char c;
 138     int argc_next;
 139
 140
 141     /*
 142      * Processing of command line arguments
 143      */
 144     while (--argc > 0 && (*++argv)[0] == '-') {
 145         argc_next = 0;
 146         /* argc_next is set to 1 if the next cmd
 147          * line argument has to be searched for;
 148          * in case that the command line argument
 149          * concerns more than one char or cmd
 150          * line argument needs a parameter
 151          */
 152         while (argc_next == 0 && (c = *++argv[0]))
 153             switch (c) {
 154
 155                 /*
 156                  * graphic output with gnuplot
 157                  */
 158             case 'g':
 159                 gnuplot = 1;
 160                 gnuplot__set = TRUE;
 161                 /*               argc_next = 1; */
 162                 break;
 163
 164             case 'h':
 165                 printf ("%s\n", argv[0]);
 166                 printf ("%s\n", HELP);
 167                 printf ("%s", USAGE);
 168 #ifdef DEBUG
 169                 printf ("\n%s", USAGE_DEBUG);
 170 #endif
 171                 exit (0);
 172
 173                 /*
 174                  * print paramaters to stdout
 175                  */
 176             case 'p':
 177                 verbose = TRUE;
 178                 break;
 179
 180                 /*
 181                  * scaling for graphic output with gnuplot
 182                  */
 183             case 's':
 184                 gnuplot_scale = atof (*++argv);
 185                 gnuplot_scale__set = TRUE;
 186                 --argc;
 187                 argc_next = 1;
 188                 break;
 189
 190                 /*
 191                  * Print version and exits
 192                  * Use Revision Control System (RCS) for version
 193                  */
 194             case 'v':
 195                 printf ("%s\n", RCSID);
 196                 exit (0);
 197
 198
 199                 /*
 200                  * long option keys
 201                  */
 202             case '-':
 203                 if (strcmp ("-help", *argv) == 0) {
 204                     printf ("\n%s", argv[0]);
 205                     printf ("\n%s", HELP);
 206                     printf ("\n%s", USAGE);
 207                     exit (0);
 208                 } else if (strcmp ("-print", *argv) == 0) {
 209                     verbose = TRUE;
 210                 } else if (strcmp ("-version", *argv) == 0) {
 211                     printf ("%s\n", RCSID);
 212                     exit (0);
 213
 214                     /*
 215                      * Rr.Col_start
 216                      */
 217                 } else if (strcmp ("-cf", *argv) == 0) {
 218                     piv_par->col_start = atoi (*++argv);
 219                     piv_par->col_start__set = TRUE;
 220                     piv_par->int_geo = GPIV_AOI;
 221                     piv_par->int_geo__set = TRUE;
 222                     --argc;
 223                     argc_next = 1;
 224
 225                     /*
 226                      * overrides Rr.Col_end
 227                      */
 228                 } else if (strcmp ("-cl", *argv) == 0) {
 229                     piv_par->col_end = atoi (*++argv);
 230                     piv_par->col_end__set = TRUE;
 231                     piv_par->int_geo = GPIV_AOI;
 232                     piv_par->int_geo__set = TRUE;
 233                     --argc;
 234                     argc_next = 1;
 235
 236                     /*
 237                      * pre-shift c (col)
 238                      */
 239                 } else if (strcmp ("-cp", *argv) == 0) {
 240                     piv_par->pre_shift_col = atoi (*++argv);
 241                     piv_par->pre_shift_col__set = TRUE;
 242                     --argc;
 243                     argc_next = 1;
 244
 245                     /*
 246                      * interrogation area weighting with gauss function
 247                      */
 248                 } else if (strcmp ("-gauss", *argv) == 0) {
 249                     piv_par->gauss_weight_ia = TRUE;
 250                     piv_par->gauss_weight_ia__set = TRUE;
 251
 252                     /*
 253                      * overrides Rr.int_size_f
 254                      */
 255                 } else if (strcmp ("-ia_size_f", *argv) == 0) {
 256                     piv_par->int_size_f = atoi (*++argv);
 257                     piv_par->int_size_f__set = TRUE;
 258                     --argc;
 259                     argc_next = 1;
 260
 261                     /*
 262                      * overrides Rr.int_size_i
 263                      */
 264                 } else if (strcmp ("-ia_size_i", *argv) == 0) {
 265                     piv_par->int_size_i = atoi (*++argv);
 266                     piv_par->int_size_i__set = TRUE;
 267                     --argc;
 268                     argc_next = 1;
 269
 270                     /*
 271                      * overrides Rr.int_shift
 272                      */
 273                 } else if (strcmp ("-ia_shift", *argv) == 0) {
 274                     piv_par->int_shift = atoi (*++argv);
 275                     piv_par->int_shift__set = TRUE;
 276                     --argc;
 277                     argc_next = 1;
 278
 279                     /*
 280                      * overrides Rr.Int_scheme
 281                      */
 282                 } else if (strcmp ("-ischeme", *argv) == 0) {
 283                     piv_par->int_scheme = atoi (*++argv);
 284                     if (piv_par->int_scheme != 0
 285                         && piv_par->int_scheme != 1
 286                         && piv_par->int_scheme != 2
 287                         && piv_par->int_scheme != 3
 288                         && piv_par->int_scheme != 4) {
 289                         gpiv_error ("%s: invalid value of Int_scheme (1, 2, 3 or 4)",
 290                                     argv[0]);
 291                     }
 292                     piv_par->int_scheme__set = TRUE;
 293                     --argc;
 294                     argc_next = 1;
 295
 296                     /*
 297                      * overrides Rr.Ifit
 298                      */
 299                 } else if (strcmp ("-ifit", *argv) == 0) {
 300                     piv_par->ifit = atoi (*++argv);
 301                     if (piv_par->ifit != -1 && piv_par->ifit != 0 &&
 302                         piv_par->ifit != 1 && piv_par->ifit != 2
 303                         && piv_par->ifit != 3) {
 304                         gpiv_error ("%s: invalid value of Ifit (1, 2 or 3)",
 305                                     argv[0]);
 306                     }
 307                     piv_par->ifit__set = TRUE;
 308                     --argc;
 309                     argc_next = 1;
 310
 311                     /*
 312                      * define interrogate at a line
 313                      */
 314                 } else if ((strcmp ("-linec", *argv) == 0)) {
 315                     piv_par->int_line_col = atoi (*++argv);
 316                     piv_par->int_line_row_start = atoi (*++argv);
 317                     piv_par->int_line_row_end =  atoi (*++argv);
 318                     piv_par->int_geo = GPIV_LINE_C;
 319                     piv_par->int_line_col__set = TRUE;
 320                     piv_par->int_line_row_start__set = TRUE;
 321                     piv_par->int_line_row_end__set = TRUE;
 322                     if (piv_par->int_geo__set == TRUE) {
 323                         gpiv_error ("%s: only a horizontal line"
 324                                     " may be defined\n"
 325                                     " not a vertical, a point or an area of"
 326                                     " interest as well\n",
 327                                     argv[0]);
 328                     } else {
 329                         piv_par->int_geo__set = TRUE;
 330                     }
 331                     --argc;
 332                     --argc;
 333                     --argc;
 334                     argc_next = 1;
 335                 } else if ((strcmp ("-liner", *argv) == 0)) {
 336                     piv_par->int_line_row = atoi (*++argv);
 337                     piv_par->int_line_col_start = atoi (*++argv);
 338                     piv_par->int_line_col_end =  atoi (*++argv);
 339                     piv_par->int_geo = GPIV_LINE_R;
 340                     piv_par->int_line_row__set = TRUE;
 341                     piv_par->int_line_col_start__set = TRUE;
 342                     piv_par->int_line_col_end__set = TRUE;
 343                     if (piv_par->int_geo__set == TRUE) {
 344                         gpiv_error ("%s: a vertical horizontal line"
 345                                     " may be defined\n"
 346                                     " not an horizontal, a point or an area of"
 347                                     " interest as well\n",
 348                                     argv[0]);
 349                     } else {
 350                         piv_par->int_geo__set = TRUE;
 351                     }
 352                     --argc;
 353                     --argc;
 354                     --argc;
 355                     argc_next = 1;
 356
 357                     /*
 358                      * cov peak # to fit
 359                      */
 360                 } else if (strcmp ("-peak", *argv) == 0) {
 361                     piv_par->peak = atoi (*++argv);
 362                     --argc;
 363                     piv_par->peak__set = TRUE;
 364                     argc_next = 1;
 365
 366                     /*
 367                      * # piv values to stdout
 368                      */
 369                 } else if (strcmp ("-p_piv", *argv) == 0) {
 370                     piv_par->print_piv = 1;
 371                     piv_par->print_piv__set = TRUE;
 372                     argc_next = 1;
 373
 374                     /*
 375                      * interrogation at point(row, col)
 376                      */
 377                 } else if (strcmp ("-point", *argv) == 0) {
 378                     piv_par->int_point_col = atoi (*++argv);
 379                     piv_par->int_point_row = atoi (*++argv);
 380                     piv_par->int_point_col__set = TRUE;
 381                     piv_par->int_point_row__set = TRUE;
 382                     piv_par->int_geo = GPIV_POINT;
 383                     if (piv_par->int_geo__set == TRUE) {
 384                         gpiv_error ("%s: only a point may be defined\n"
 385                                     " not an horizontal, a vertical line"
 386                                     " or an area of interest as well\n",
 387                                     argv[0]);
 388                     } else {
 389                         piv_par->int_geo__set = TRUE;
 390                     }
 391                     --argc;
 392                     --argc;
 393                     argc_next = 1;
 394
 395                     /*
 396                      * overrides Rr.Row_start
 397                      */
 398                 } else if (strcmp ("-rf", *argv) == 0) {
 399                     piv_par->row_start = atoi (*++argv);
 400                     piv_par->row_start__set = TRUE;
 401                     piv_par->int_geo = GPIV_AOI;
 402                     piv_par->int_geo__set = TRUE;
 403                     --argc;
 404                     argc_next = 1;
 405
 406                     /*
 407                      * overrides Rr.Row_end
 408                      */
 409                 } else if (strcmp ("-rl", *argv) == 0) {
 410                     piv_par->row_end = atoi (*++argv);
 411                     piv_par->row_end__set = TRUE;
 412                     piv_par->int_geo = GPIV_AOI;
 413                     piv_par->int_geo__set = TRUE;
 414                     --argc;
 415                     argc_next = 1;
 416
 417                     /*
 418                      * pre-shift r (row)
 419                      */
 420                 } else if (strcmp ("-rp", *argv) == 0) {
 421                     piv_par->pre_shift_row = atoi (*++argv);
 422                     piv_par->pre_shift_row__set = TRUE;
 423                     --argc;
 424                     argc_next = 1;
 425
 426                     /*
 427                      * Symmetric phase only filtering
 428                      */
 429                 } else if (strcmp ("-spof", *argv) == 0) {
 430                     piv_par->spof_filter = TRUE;
 431                     piv_par->spof_filter__set = TRUE;
 432
 433                     /*
 434                      * validation parameter: residu type
 435                      */
 436                 } else if (strcmp ("-val_r", *argv) == 0) {
 437                     valid_par->residu_type = atoi (*++argv);
 438                     valid_par->residu_type__set = TRUE;
 439                     --argc;
 440                     argc_next = 1;
 441                     if (valid_par->residu_type != 0
 442                         && valid_par->residu_type != 1
 443                         && valid_par->residu_type != 2) {
 444                         gpiv_error ("%s: invalid value of VALID.Residu_type (0, 1, or 2)",
 445                                     argv[0]);
 446                     }
 447
 448                     /*
 449                      * validation parameter: substitution type
 450                      */
 451                 } else if (strcmp ("-val_s", *argv) == 0) {
 452                     valid_par->subst_type = atoi (*++argv);
 453                     valid_par->subst_type__set = TRUE;
 454                     --argc;
 455                     argc_next = 1;
 456                     if (valid_par->subst_type != 0
 457                         && valid_par->subst_type != 1
 458                         && valid_par->subst_type != 2
 459                         && valid_par->subst_type != 3) {
 460                         gpiv_error ("%s: invalid value of VALID.Subst_type (0, 1, 2 or 3)",
 461                                     argv[0]);
 462                     }
 463
 464                     /*
 465                      * validation parameter: threshold value
 466                      */
 467                 } else if (strcmp ("-val_t", *argv) == 0) {
 468                     valid_par->residu_max = atof (*++argv);
 469                     valid_par->residu_max__set = TRUE;
 470                     --argc;
 471                     argc_next = 1;
 472
 473                 } else {
 474                     gpiv_error ("%s: unknown option: %s", argv[0], *argv);
 475                 }
 476                 argc_next = 1;
 477                 break;
 478
 479             default:
 480                 gpiv_error ("%s: unknown option: %s", argv[0], *argv);
 481                 break;
 482
 483             }
 484     }
 485
 486
 487     /*
 488      * Check if filename or stdin /stdout is used
 489      */
 490     if (argc == 1) {
 491         use_stdin_stdout = FALSE;
 492         strcpy(fname, argv[argc - 1]);
 493     } else if (argc == 0) {
 494         use_stdin_stdout = TRUE;
 495         verbose = FALSE;
 496         if (gnuplot == 1) {
 497             gpiv_error ("%s: an image filename for input has to be used "
 498                         "in combination with 'gnuplot'",
 499                         argv[0]);
 500         }
 501     } else {
 502         gpiv_error("\n%s: unknown argument: %s: %s", argv[0], *argv, USAGE);
 503     }
 504
 505
 506
 507 }
 508
 509
 510
 511 static gchar *
 512 make_fname (FILE * fp_par_out,
 513             char *fname_in,
 514             char *fname_out,
 515             GpivPivPar *piv_par
 516             )
 517 /*-----------------------------------------------------------------------------
 518  * generates filenames
 519  */
 520 {
 521     gchar *err_msg = NULL;
 522     gchar *fname_base = NULL;
 523
 524     if (fname_in == NULL ) {
 525         err_msg = "make_fname: \"fname_in == NULL\"";
 526         return (err_msg);
 527     }
 528
 529     /*
 530      * Stripping filename
 531      */
 532     fname_base = g_strdup(fname_in);
 533     strtok(fname_base, ".");
 534
 535     /*
 536      * filenames for output PIV data
 537      */
 538     if (fname_out != NULL ) {
 539         gpiv_io_make_fname (fname_base, GPIV_EXT_PIV, fname_out);
 540         if (piv_par->print_piv == 0) {
 541            if (verbose)
 542                 printf ("\n# outputfile is: %s\n", fname_out);
 543             fprintf (fp_par_out, "\n# outputfile is: %s\n", fname_out);
 544         }
 545     }
 546
 547     g_free (fname_base);
 548     return (err_msg);
 549 }
 550
 551
 552 #ifdef ENABLE_MPI
 553
 554 static GpivImage *
 555 alloc_image_struct(void)
 556 /*-----------------------------------------------------------------------------
 557  * Allocates GpivImage structure
 558  */
 559 {
 560     GpivImage *image = g_new0 (GpivImage, 1);
 561     GpivImagePar *image_par = g_new0 (GpivImagePar, 1);;
 562
 563
 564     image->header = image_par;
 565
 566
 567     return image;
 568 }
 569
 570
 571 static GpivPivPar *
 572 alloc_pivpar_struct(void)
 573 /*-----------------------------------------------------------------------------
 574  * Allocates GpivPivPar structure
 575  */
 576 {
 577     GpivPivPar *piv_par = g_new0 (GpivPivPar, 1);
 578     return piv_par;
 579 }
 580
 581
 582 static GpivValidPar *
 583 alloc_validpar_struct(void)
 584 /*-----------------------------------------------------------------------------
 585  * Allocates GpivValidPar structure
 586  */
 587 {
 588     GpivValidPar *valid_par = g_new0 (GpivValidPar, 1);
 589     return valid_par;
 590 }
 591
 592 #endif /* ENABLE_MPI */
 593
 594
 595 int
 596 main (int argc,
 597       char *argv[]
 598       )
 599 /*-----------------------------------------------------------------------------
 600  * Main routine of rr to calculate estimators of particle image displacements
 601  * from PIV images
 602  */
 603 {
 604     char *err_msg = NULL;
 605     FILE *fp_out = NULL, *fp_par_out = NULL;
 606     char fname_in[GPIV_MAX_CHARS],
 607         fname_out[GPIV_MAX_CHARS],
 608         fname_par[GPIV_MAX_CHARS];
 609
 610     GpivImage *image = NULL;
 611     GpivPivData *piv_data = NULL;
 612     GpivPivPar *piv_par = g_new0 (GpivPivPar, 1);
 613     GpivValidPar *valid_par = g_new0 (GpivValidPar, 1);
 614
 615 #ifdef _OPENMP
 616     int np, nprocs, num_threads;
 617 #endif /* _OPENMP */
 618
 619 #ifdef ENABLE_MPI
 620     /*     gint argc_mpi; */
 621     /*     gchar *argv_mpi[10]; */
 622     int nprocs, rank = 0, namelen, root = 0;
 623     char processor_name[MPI_MAX_PROCESSOR_NAME];
 624 #endif /* ENABLE_MPI */
 625
 626 #ifdef ENABLE_MPI
 627     /*
 628      * Initializing MPI
 629      */
 630     MPI_Init(&argc, &argv);
 631     MPI_Comm_rank(MPI_COMM_WORLD, &rank);
 632     MPI_Comm_size(MPI_COMM_WORLD, &nprocs);
 633     MPI_Get_processor_name(processor_name, &namelen);
 634
 635 /*     gpiv_warning("main:: 0 rank=%d nprocs=%d on %s", rank, nprocs, processor_name); */
 636     /*     argc_mpi = argc - 4; */
 637     /*     for (i = 0; i<argc_mpi; i++) { */
 638     /*         gpiv_warning("main:: rank=%d argv[%d]=%s", rank, i, argv[i]); */
 639     /*         argv_mpi[i] = argv[i]; */
 640     /*     } */
 641 #endif /* ENABLE_MPI */
 642
 643     /*
 644      * Initializing parameters
 645      */
 646     gpiv_piv_parameters_set (piv_par, FALSE);
 647     gpiv_valid_parameters_set (valid_par, FALSE);
 648
 649 #ifdef _OPENMP
 650     /*
 651      * Initializing OMP
 652      */
 653     np = omp_get_num_procs();
 654     omp_set_num_threads(np);
 655     num_threads = omp_get_max_threads();
 656 #endif /* _OPENMP */
 657
 658 #ifdef ENABLE_MPI
 659     if (rank == 0) {
 660 #endif /* ENABLE_MPI */
 661
 662         command_args (argc,
 663                       argv,
 664                       fname_in,
 665                       piv_par,
 666                       valid_par);
 667
 668         if (verbose) {
 669             printf ("\n# %s\n# Command line options:\n", RCSID);
 670             gpiv_piv_print_parameters (stdout, piv_par);
 671             gpiv_valid_print_parameters (stdout, valid_par);
 672         }
 673
 674         /*
 675          * Reads piv and valid parameters from PARFILE and/or resource files
 676          * if not overridden by the commandline options.
 677          * Check if all parameters have been read. Else, set as libgpiv defaults
 678          */
 679         gpiv_scan_parameter (GPIV_PIVPAR_KEY, PARFILE, piv_par, verbose);
 680         gpiv_scan_resourcefiles (GPIV_PIVPAR_KEY, piv_par, verbose);
 681         if ((err_msg =
 682              gpiv_piv_check_parameters_read (piv_par, NULL))
 683             != NULL) gpiv_warning ("%s: %s", argv[0], err_msg);
 684
 685         gpiv_scan_parameter (GPIV_VALIDPAR_KEY, PARFILE, valid_par, verbose);
 686         gpiv_scan_resourcefiles (GPIV_VALIDPAR_KEY, valid_par, verbose);
 687         if ((err_msg =
 688              gpiv_valid_check_parameters_read (valid_par, NULL))
 689             != NULL) gpiv_warning ("%s: %s", argv[0], err_msg);
 690
 691        /*
 692          * Creates file names if not stdin / stdout are used
 693          * and save parameters
 694          */
 695         if (use_stdin_stdout == FALSE) {
 696             gchar *fname_base = g_strdup (fname_in);
 697             strtok(fname_base, ".");
 698             gpiv_io_make_fname (fname_base, GPIV_EXT_PAR, fname_par);
 699             g_free (fname_base);
 700
 701             if ((fp_par_out = fopen (fname_par, "w")) == NULL) {
 702                 gpiv_error ("%s error: failure opening %s for input",
 703                             argv[0], fname_par);
 704             }
 705             if ((err_msg =
 706                  make_fname (fp_par_out, fname_in, fname_out, piv_par))
 707                 != NULL) {
 708                 gpiv_error ("%s: %s\n", argv[0], err_msg);
 709             }
 710
 711
 712             /*
 713              * Write parameters to parameterfile
 714              */
 715             gpiv_piv_print_parameters (fp_par_out, piv_par);
 716             gpiv_valid_print_parameters (fp_par_out, valid_par);
 717             fclose (fp_par_out);
 718         }
 719
 720
 721         /*
 722          * reads the image
 723          */
 724
 725 #ifdef USE_MTRACE
 726         /*     fprintf(stderr, "Calling interrogate with mtrace()\n"); */
 727         mtrace();
 728 #endif
 729         if (use_stdin_stdout) {
 730 #ifdef DEBUG
 731             gpiv_warning ("gpiv_rr:: using stdin");
 732 #endif
 733             if ((image = gpiv_read_image (stdin)) == NULL) {
 734                 gpiv_error ("%s: gpiv_fread_image_ni", argv[0]);
 735             }
 736         } else {
 737 #ifdef DEBUG
 738             gpiv_warning ("rr:: using fp, fname = %s", fname_in);
 739 #endif
 740             if ((image = gpiv_fread_image (fname_in)) == NULL) {
 741                 gpiv_error ("%s: gpiv_fread_image", argv[0]);
 742             }
 743         }
 744
 745 #ifdef ENABLE_MPI
 746     } else {
 747         /*
 748          * allocates the memory of structures for threads of which rank != 0
 749          */
 750         image = alloc_image_struct();
 751         piv_par = alloc_pivpar_struct();
 752         valid_par = alloc_validpar_struct();
 753
 754     }
 755
 756
 757     if (MPI_Bcast(&verbose, 1, MPI_INT, root, MPI_COMM_WORLD)
 758         != MPI_SUCCESS) {
 759         gpiv_error ("gpiv_rr: An error ocurred when calling MPI_Bcast");
 760     }
 761     gpiv_img_mpi_bcast_image (image, TRUE);
 762     if ((err_msg = gpiv_piv_mpi_bcast_pivpar (piv_par)) != NULL) {
 763         gpiv_error (err_msg);
 764     }
 765     gpiv_valid_mpi_bcast_validpar (valid_par);
 766
 767 #endif /* ENABLE_MPI */
 768 #ifdef DEBUG
 769     gpiv_warning ("gpiv_rr:: calling gpiv_piv_interrogate_img");
 770 #endif
 771
 772
 773 /*
 774  * Interrogating image
 775  */
 776     if ((piv_data = gpiv_piv_interrogate_img (image, piv_par, valid_par, verbose))
 777         == NULL) {
 778         gpiv_error ("%s: gpiv_piv_interrogate_img failed", argv[0]);
 779     }
 780
 781
 782 #ifdef ENABLE_MPI
 783     if (rank == 0) {
 784 #endif
 785         /*
 786          * Displays PIV data
 787          */
 788         if (gnuplot) gpiv_piv_gnuplot (fname_in, gnuplot_scale,
 789                                        GNUPLOT_DISPLAY_COLOR, GNUPLOT_DISPLAY_SIZE,
 790                                        image->header, piv_par, piv_data);
 791
 792         /*
 793          * Add comments to PIV data and saving
 794          */
 795         piv_data->comment = g_strdup_printf ("# Software: %s\n", RCSID);
 796         piv_data->comment = gpiv_add_datetime_to_comment (piv_data->comment);
 797         piv_data->comment = g_strconcat (piv_data->comment,
 798                                          "# Data type: Particle Image Velocities\n",
 799                                          NULL);
 800
 801
 802         if ((piv_par->print_piv == TRUE && use_stdin_stdout == FALSE)
 803             || use_stdin_stdout == TRUE) {
 804             fp_out = stdout;
 805             if ((err_msg =
 806                  gpiv_write_pivdata (fp_out, piv_data, FALSE))
 807                 != NULL) gpiv_error ("%s: %s", argv[0], err_msg);
 808
 809         } else {
 810             if ((fp_out = fopen (fname_out, "wb")) == NULL) {
 811                 gpiv_error ("%s: failing fopen %s", argv[0], fname_out);
 812             }
 813             if ((err_msg =
 814                  gpiv_write_pivdata (fp_out, piv_data, FALSE))
 815                 != NULL) gpiv_error ("%s: %s", argv[0], err_msg);
 816             fclose (fp_out);
 817
 818         }
 819 #ifdef ENABLE_MPI
 820     }
 821 #endif
 822
 823     /*
 824      * Free memory of images and data
 825      */
 826     gpiv_free_pivdata (piv_data);
 827     if (image->frame1 != NULL && image->frame2 != NULL) {
 828         gpiv_free_img (image);
 829     }
 830
 831
 832 #ifdef USE_MTRACE
 833     muntrace();
 834 #endif
 835
 836 /* #endif */ /* __DISABLE */
 837 #ifdef ENABLE_MPI
 838 #ifdef DEBUG
 839         gpiv_warning("main:: rank=%d: Calling MPI_Finalize()", rank);
 840 #endif /* DEBUG */
 841         MPI_Finalize();
 842 #endif /* ENABLE+MPI */
 843         exit (0);
 844 }