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