src/piveval_interrogate.c

   1 /* -*- Mode: C; indent-tabs-mode: nil; c-basic-offset: 4 c-style: "K&R" -*- */
   2
   3 /*----------------------------------------------------------------------
   4
   5   gpiv - Graphic program for Particle Image Velocimetry, based on gtk/gnome
   6           libraries.
   7
   8    Copyright (C) 2006, 2007, 2008
   9    Gerber van der Graaf
  10
  11    This file is part of gpiv.
  12
  13    Gpiv is free software; you can redistribute it and/or modify
  14    it under the terms of the GNU General Public License as published by
  15    the Free Software Foundation; either version 2, or (at your option)
  16    any later version.
  17
  18    This program is distributed in the hope that it will be useful,
  19    but WITHOUT ANY WARRANTY; without even the implied warranty of
  20    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  21    GNU General Public License for more details.
  22
  23    You should have received a copy of the GNU General Public License
  24    along with this program; if not, write to the Free Software Foundation,
  25    Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
  26
  27 ----------------------------------------------------------------------*/
  28
  29 /* $Log: piveval_interrogate.c,v $
  30 /* Revision 1.4  2008-10-09 15:19:03  gerber
  31 /* Release 0.6.0 using parallelized code
  32 /*
  33 /* Revision 1.3  2008-10-09 14:43:37  gerber
  34 /* paralellized with OMP and MPI
  35 /*
  36 /* Revision 1.2  2008-09-25 13:32:22  gerber
  37 /* Adapted for use on cluster (using MPI/OMP) parallelised gpiv_rr from gpivtools)
  38 /*
  39 /* Revision 1.1  2008-09-16 10:17:36  gerber
  40 /* added piveval_interrogate routines
  41 /*
  42  */
  43
  44 #include "gpiv_gui.h"
  45 /* #include "utils.h" */
  46 #include "piveval.h"
  47 #include "piveval_interrogate.h"
  48
  49 #undef USE_MTRACE
  50 #ifdef USE_MTRACE
  51 #include <mcheck.h>
  52 #endif
  53
  54
  55 static void
  56 report_progress (const guint index_y,
  57                  const guint index_x,
  58                  gdouble *progress_value_prev,
  59                  const GpivPivData *piv_data,
  60                  const GpivPivPar *piv_par,
  61                  const guint sweep,
  62                  const gfloat cum_residu,
  63                  const GpivConsole *gpiv,
  64                  const Display *disp
  65                  );
  66
  67 static GpivPivData *
  68 alloc_pivdata_gridgen (const GpivImagePar *image_par,
  69                        const GpivPivPar *piv_par
  70                        );
  71
  72 static gchar *
  73 update_pivdata_imgdeform_zoff (const GpivImage *image,
  74                                GpivImage *lo_image,
  75                                const GpivPivPar *piv_par,
  76                                const GpivValidPar *valid_par,
  77                                GpivPivData *piv_data,
  78                                GpivPivData *lo_piv_data,
  79                                gfloat *cum_residu,
  80                                gboolean *cum_residu_reached,
  81                                gfloat *sum_dxdy,
  82                                gfloat *sum_dxdy_old,
  83                                gboolean isi_last,
  84                                gboolean grid_last,
  85                                gboolean sweep_last,
  86                                gboolean verbose
  87                                );
  88
  89 #ifdef ENABLE_MPI
  90 static GpivPivData *
  91 exec_piv_mpi (GpivImage *image,
  92               GpivPivPar *piv_par,
  93               GpivValidPar *valid_par,
  94               GpivConsole *gpiv
  95               );
  96 #endif /* ENABLE_MPI */
  97
  98 /*
  99  * Program-wide public piv interrogation functions
 100  */
 101
 102
 103 void
 104 exec_piv (GpivConsole *gpiv
 105          )
 106 /*-----------------------------------------------------------------------------
 107  * Performs the execution of image interrogation for PIV
 108  */
 109 {
 110     char *err_msg = NULL;
 111     char message[2 * GPIV_MAX_CHARS];
 112     guint nx = 0, ny = 0;
 113
 114
 115     if (display_act == NULL
 116         || display_act->img->exist_img == FALSE
 117         || cancel_process) {
 118         err_msg =
 119             _("At first, open an image. \n"
 120               "Than we'll further see what will happen.");
 121         g_warning (err_msg);
 122         warning_gpiv (err_msg);
 123         return;
 124     }
 125
 126 /*
 127  * Free memory of pivdata and clean the display from its vectors
 128  */
 129     if (display_act->pida->exist_piv
 130 /*             && (m_select != SINGLE_AREA_MS */
 131 /*                 || m_select != DRAG_AREA_MS) */
 132         ) {
 133         destroy_all_vectors (display_act->pida);
 134         gpiv_free_pivdata (display_act->pida->piv_data);
 135         display_act->pida->exist_piv = FALSE;
 136         display_act->pida->averaged_piv = FALSE;
 137         if (display_act->pida->scaled_piv) {
 138             gpiv_free_pivdata (display_act->pida->piv_data_scaled);
 139             display_act->pida->scaled_piv = FALSE;
 140         }
 141     }
 142
 143 /*
 144  * Free eventually existing memory of vor_data and cleanup the display
 145  * as they do not belong to the piv data anymore
 146  */
 147     free_post_bufmems (display_act);
 148
 149     exec_process = TRUE;
 150
 151 /*
 152  * Set mouse selection to None for using correct AOI
 153  */
 154 /*         if (m_select != NO_MS) { */
 155 /*             gtk_toggle_button_set_state (GTK_TOGGLE_BUTTON */
 156 /*                                          (gpiv->piveval->radiobutton_mouse_1),  */
 157 /*                                          TRUE); */
 158 /*         } */
 159
 160 /*
 161  * Setting interrogation scheme to GPIV_ZERO_OFF_CENTRAL if image deformation
 162  * is impossible with too small (initial or final) grid.
 163  */
 164     if (gl_piv_par->int_scheme == GPIV_IMG_DEFORM) {
 165         gpiv_piv_count_pivdata_fromimage (display_act->img->image->header,
 166                                           gl_piv_par, &nx, &ny);
 167         if (nx  < 2 || ny < 2) {
 168             g_snprintf (message, 2 * GPIV_MAX_CHARS,
 169                         _("Image deformation is impossibe with grid of nx = %d ny =%d.\n\
 170 Setting Interrogation scheme to Central difference.\n                   \
 171 This will be reset automatically."),
 172                         nx, ny);
 173             warning_gpiv ("%s", message);
 174             gtk_toggle_button_set_state (GTK_TOGGLE_BUTTON
 175                                          (gpiv->piveval->radiobutton_centraldiff),
 176                                          TRUE);
 177             int_scheme_autochanged = TRUE;
 178         }
 179     }
 180
 181 /*
 182  * checking parameters and interrogating image by local function
 183  */
 184     display_act->pida->piv_par = gpiv_piv_cp_parameters (gl_piv_par);
 185     if ((err_msg =
 186          gpiv_piv_testadjust_parameters (display_act->img->image->header,
 187                                          display_act->pida->piv_par))
 188         != NULL) {
 189         warning_gpiv ("%s", err_msg);
 190         return;
 191     }
 192
 193     display_act->pida->valid_par = gpiv_valid_cp_parameters (gl_valid_par);
 194     gpiv_valid_print_parameters (stdout, display_act->pida->valid_par);
 195     if ((err_msg =
 196          gpiv_valid_testadjust_parameters (display_act->pida->valid_par))
 197         != NULL) {
 198         warning_gpiv ("%s", err_msg);
 199         return;
 200     }
 201     gpiv_valid_print_parameters (stdout, display_act->pida->valid_par);
 202
 203 #ifdef USE_LIBGPIV_INTERR
 204 #undef USE_LIBGPIV_INTERR
 205 #endif
 206
 207     if ((display_act->pida->piv_data =
 208 #ifdef ENABLE_MPI
 209 /* Calling extern tool gpiv_rr from gpivtools with MPI enabled */
 210          exec_piv_mpi(display_act->img->image,
 211                       display_act->pida->piv_par,
 212                       display_act->pida->valid_par,
 213                       gpiv)
 214
 215 #elif USE_LIBGPIV_INTERR
 216 /* Libgpiv's function for image interrogation */
 217          gpiv_piv_interrogate_img (display_act->img->image,
 218                                    display_act->pida->piv_par,
 219                                    display_act->pida->valid_par,
 220                                    TRUE)
 221 #else
 222 /* Gpiv's function for image interrogation */
 223          interrogate_img (display_act->img->image,
 224                           display_act->pida->piv_par,
 225                           display_act->pida->valid_par,
 226                           gpiv)
 227 #endif
 228          ) == NULL) {
 229         warning_gpiv ("exec_piv: failing interrogate_img");
 230         return;
 231     }
 232
 233 #ifdef USE_LIBGPIV_INTERR
 234 #undef USE_LIBGPIV_INTERR
 235 #endif
 236
 237     display_act->pida->exist_piv = TRUE;
 238
 239 /*
 240  *  Adding some comment to piv_data
 241  */
 242     display_act->pida->piv_data->comment =
 243         g_strdup_printf ("# Software: %s %s\n", PACKAGE, VERSION);
 244     display_act->pida->piv_data->comment =
 245         gpiv_add_datetime_to_comment (display_act->pida->piv_data->comment);
 246     display_act->pida->piv_data->comment =
 247         g_strconcat (display_act->pida->piv_data->comment,
 248                      "# Data type: Particle Image Velocities\n", NULL);
 249
 250 /*
 251  * Drawing and displaying PIV vectors
 252  */
 253     if (gl_piv_par->int_geo == GPIV_POINT
 254         || m_select == SINGLE_AREA_MS
 255         || m_select == SINGLE_POINT_MS
 256         || m_select == DRAG_MS) {
 257     } else {
 258         if (display_act->display_piv) {
 259             create_all_vectors (display_act->pida);
 260         }
 261         display_act->display_piv = TRUE;
 262     }
 263
 264 /*
 265  * Some settings for displaying features
 266  * Update vectors to correct for colors/gray-scale
 267  */
 268     display_act->pida->scaled_piv = FALSE;
 269     display_act->pida->saved_piv = FALSE;
 270     display_act->pida->averaged_piv = FALSE;
 271     display_act->pida->exist_cov = TRUE;
 272     update_all_vectors (display_act->pida);
 273     exec_process = FALSE;
 274
 275 /*
 276  * Resetting interrogation scheme
 277  */
 278     if (int_scheme_autochanged) {
 279         int_scheme_autochanged = FALSE;
 280         gtk_toggle_button_set_state (GTK_TOGGLE_BUTTON
 281                                      (gpiv->piveval->radiobutton_imgdeform),
 282                                      TRUE);
 283     }
 284
 285     gtk_progress_set_value (GTK_PROGRESS (gnome_appbar_get_progress
 286                                           (GNOME_APPBAR (gpiv->appbar))),
 287                             0.0);
 288
 289 }
 290
 291
 292 GpivPivData *
 293 interrogate_img (const GpivImage *image,
 294                  const GpivPivPar *piv_par,
 295                  const GpivValidPar *valid_par,
 296                  GpivConsole *gpiv
 297                  )
 298 /* ----------------------------------------------------------------------------
 299  * PIV interrogation of an image pair at an entire grid or single point
 300  * Similar to gpiv_piv_interr_img, but adapted for the GUI
 301  */
 302 {
 303     GpivPivData *piv_data = NULL;       /* piv data to be returned */
 304     gchar *err_msg = NULL;              /* error message */
 305     guint index_x = 0, index_y = 0;     /* array indices */
 306
 307     /*
 308      * Local variables with prefix lo_ to distinguish from global or from
 309      * parameter list
 310      */
 311     GpivImage *lo_image = NULL;         /* local image that might be deformed */
 312     GpivPivData *lo_piv_data = NULL;    /* local piv data */
 313     GpivPivPar *lo_piv_par = NULL;      /* local piv parameters */
 314
 315     gfloat **intreg1 = display_act->pida->intreg1; /* first interrogation area */
 316     gfloat **intreg2 = display_act->pida->intreg2; /* second interrogation area */
 317     guint int_size_0;                   /* zero-padded interrogation area size */
 318
 319     GpivCov *cov = NULL;                /* covariance */
 320     guint sweep = 1;                    /* itaration counter */
 321     gboolean grid_last = FALSE;         /* flag if final grid refinement has been reached */
 322     gboolean isi_last = FALSE;          /* flag if final interrogation area shift has been reached */
 323     gboolean cum_residu_reached = FALSE;/* flag if max. cumulative residu has been reached */
 324     gboolean sweep_last = FALSE;        /* perform the last iteration sweep */
 325     gboolean sweep_stop = FALSE;        /* stop the current iteration at the end */
 326     gfloat sum_dxdy = 0.0, sum_dxdy_old = 0.0;  /* */
 327     gfloat cum_residu = 914.6;          /* initial, large, arbitrary cumulative residu */
 328     gdouble progress = 0.0;             /* for monitoring calculation progress */
 329
 330
 331     /*
 332      * Testing parameters on consistency and initializing derived
 333      * parameters/variables
 334      */
 335     if ((err_msg =
 336          gpiv_piv_testonly_parameters (image->header, piv_par))
 337         != NULL) {
 338         warning_gpiv ("%s", err_msg);
 339         return NULL;
 340     }
 341
 342     if ((err_msg =
 343          gpiv_valid_testonly_parameters (valid_par))
 344         != NULL) {
 345         warning_gpiv ("%s", err_msg);
 346         return NULL;
 347     }
 348
 349     /*
 350      * Local (actualized) parameters
 351      * Setting initial parameters and variables for adaptive grid and
 352      * Interrogation Area dimensions
 353      */
 354     lo_piv_par = gpiv_piv_cp_parameters (piv_par);
 355
 356     if (lo_piv_par->int_scheme == GPIV_ZERO_OFF_FORWARD
 357         || lo_piv_par->int_scheme == GPIV_ZERO_OFF_CENTRAL
 358         || lo_piv_par->int_scheme == GPIV_IMG_DEFORM
 359         || lo_piv_par->int_size_i > lo_piv_par->int_size_f) {
 360         lo_piv_par->int_size_f = lo_piv_par->int_size_i;
 361         sweep_last = FALSE;
 362     } else {
 363         sweep_last = TRUE;
 364     }
 365
 366     if (lo_piv_par->int_shift < lo_piv_par->int_size_i / GPIV_SHIFT_FACTOR) {
 367         lo_piv_par->int_shift = lo_piv_par->int_size_i / GPIV_SHIFT_FACTOR;
 368     }
 369
 370     /*
 371      * A copy of the image and PIV data are needed when image deformation is used.
 372      * To keep the algorithm simple, copies are made unconditionally.
 373      */
 374     lo_image = gpiv_cp_img (image);
 375     piv_data = alloc_pivdata_gridgen (image->header, lo_piv_par);
 376     lo_piv_data = gpiv_cp_pivdata (piv_data);
 377     gpiv_0_pivdata (lo_piv_data);
 378
 379     /*
 380      * Reads eventually existing fftw wisdom
 381      */
 382     gpiv_fread_fftw_wisdom(1);
 383     gpiv_fread_fftw_wisdom(-1);
 384
 385 #ifdef USE_MTRACE
 386     mtrace();
 387 #endif
 388
 389     while (sweep <= GPIV_MAX_PIV_SWEEP
 390            && !sweep_stop
 391            && !cancel_process) {
 392
 393         /*
 394          * Memory allocation of interrogation area's and covariance.
 395          * These memory chunks are allocated here to optimize calculation
 396          * speed and for eventually monitoring their contents.
 397          */
 398         int_size_0 = GPIV_ZEROPAD_FACT * lo_piv_par->int_size_i;
 399         intreg1 = gpiv_matrix (int_size_0, int_size_0);
 400         intreg2 = gpiv_matrix (int_size_0, int_size_0);
 401         cov = gpiv_alloc_cov (int_size_0, image->header->x_corr);
 402         display_act->pida->cov = cov;
 403
 404         /*
 405          * Interrogates a single interrogation area
 406          */
 407         if (m_select != SINGLE_AREA_MS
 408             && m_select != DRAG_MS) {
 409             destroy_all_vectors (display_act->pida);
 410             gnome_canvas_update_now (GNOME_CANVAS (display_act->canvas));
 411         }
 412
 413         if (lo_piv_par->int_geo == GPIV_POINT
 414             || m_select == SINGLE_AREA_MS
 415             || m_select == SINGLE_POINT_MS
 416             || m_select == DRAG_MS) {
 417
 418             if ((err_msg =
 419                  gpiv_piv_interrogate_ia (m_select_index_y,
 420                                           m_select_index_x,
 421                                           lo_image,
 422                                           lo_piv_par,
 423                                           sweep,
 424                                           sweep_last,
 425                                           intreg1,
 426                                           intreg2,
 427                                           cov,
 428                                           lo_piv_data
 429                                           ))
 430                 != NULL) {
 431                 gpiv_free_img (lo_image);
 432                 gpiv_free_pivdata (lo_piv_data);
 433                 gpiv_free_pivdata (piv_data);
 434                 gpiv_free_matrix (intreg1);
 435                 gpiv_free_matrix (intreg2);
 436                 gpiv_free_cov (cov);
 437                 error_gpiv ("interrogate_img: %s", err_msg);
 438             }
 439
 440             /*
 441              * display piv values, draw interrogation areas and
 442              * covariance function
 443              */
 444             if (gpiv_var->piv_disproc == TRUE) {
 445                 display_piv_vector (m_select_index_y,
 446                                     m_select_index_x,
 447                                     piv_data,
 448                                     gpiv->piveval);
 449                 display_img_intreg1 (intreg1,
 450                                      lo_piv_par->int_size_i,
 451                                      gpiv->piveval);
 452                 display_img_intreg2 (intreg2,
 453                                      lo_piv_par->int_size_i,
 454                                      gpiv->piveval);
 455                 display_img_cov (cov,
 456                                  lo_piv_par->int_size_i,
 457                                  gpiv->piveval);
 458             }
 459
 460
 461         } else {
 462
 463             /*
 464              * Interrogates at a rectangular grid of points within the Area Of
 465              * Interest of the image
 466              */
 467             for (index_y = 0; index_y < lo_piv_data->ny; index_y++) {
 468                 for (index_x = 0; index_x < lo_piv_data->nx; index_x++) {
 469                     if (cancel_process) break;
 470
 471                     /*
 472                      * Interrogates a single interrogation area.
 473                      */
 474                     if ((err_msg =
 475                          gpiv_piv_interrogate_ia (index_y,
 476                                                   index_x,
 477                                                   lo_image,
 478                                                   lo_piv_par,
 479                                                   sweep,
 480                                                   sweep_last,
 481                                                   intreg1,
 482                                                   intreg2,
 483                                                   cov,
 484                                                   lo_piv_data
 485                                                   ))
 486                         != NULL) {
 487                         gpiv_free_img (lo_image);
 488                         gpiv_free_pivdata (lo_piv_data);
 489                         gpiv_free_pivdata (piv_data);
 490                         gpiv_free_matrix (intreg1);
 491                         gpiv_free_matrix (intreg2);
 492                         gpiv_free_cov (cov);
 493                         error_gpiv ("interrogate_img: %s", err_msg);
 494                     }
 495
 496 #ifdef DEBUG
 497                     gpiv_warning("interrogate_img:: back from gpiv_piv_interr_ia: sweep=%d x[%d][%d]=%f y[%d][%d]=%f dx[%d][%d]=%f dy[%d][%d]=%f snr=%f p_no=%d",
 498                                  sweep,
 499                                  index_y, index_x, lo_piv_data->point_x[index_y][index_x],
 500                                  index_y, index_x, lo_piv_data->point_y[index_y][index_x],
 501                                  index_y, index_x, lo_piv_data->dx[index_y][index_x],
 502                                  index_y, index_x, lo_piv_data->dy[index_y][index_x],
 503                                  lo_piv_data->snr[index_y][index_x],
 504                                  lo_piv_data->peak_no[index_y][index_x]
 505                                  );
 506 #endif
 507                     /*
 508                      * Printing the progress of processing
 509                      */
 510                     report_progress (index_y,
 511                                      index_x,
 512                                      &progress,
 513                                      piv_data,
 514                                      lo_piv_par,
 515                                      sweep,
 516                                      cum_residu,
 517                                      gpiv,
 518                                      display_act
 519                                      );
 520                     /*
 521                      * Draw interrogation areas, covariance function,
 522                      * display piv vector to monitor the process.
 523                      * Includes report_rogress to include
 524                      * as an argument in gpiv_piv_interr_img
 525                      */
 526                     if (gpiv_var->piv_disproc == TRUE) {
 527                         display_piv_vector (index_y,
 528                                             index_x,
 529                                             piv_data,
 530                                             gpiv->piveval);
 531                         display_img_intreg1 (intreg1,
 532                                              lo_piv_par->int_size_i,
 533                                              gpiv->piveval);
 534                         display_img_intreg2 (intreg2,
 535                                              lo_piv_par->int_size_i,
 536                                              gpiv->piveval);
 537                         display_img_cov (cov,
 538                                          lo_piv_par->int_size_i,
 539                                          gpiv->piveval);
 540                     }
 541
 542
 543                 }
 544             }
 545
 546         }
 547
 548         /*
 549          * De-allocating memory: other (smaller) sizes are eventually needed
 550          * for a next iteration sweep
 551          */
 552         gpiv_free_matrix (intreg1);
 553         gpiv_free_matrix (intreg2);
 554         gpiv_free_cov (cov);
 555
 556         if (sweep_last) {
 557             sweep_stop = TRUE;
 558         }
 559
 560         if ((lo_piv_par->int_scheme == GPIV_IMG_DEFORM
 561              || lo_piv_par->int_scheme == GPIV_ZERO_OFF_FORWARD
 562              || lo_piv_par->int_scheme == GPIV_ZERO_OFF_CENTRAL)
 563             /* BUGFIX: crashes with single point */
 564             && (lo_piv_par->int_geo != GPIV_POINT
 565                 && m_select != SINGLE_AREA_MS
 566                 && m_select != SINGLE_POINT_MS
 567                 && m_select != DRAG_MS)
 568             ) {
 569
 570             if ((err_msg =
 571                  update_pivdata_imgdeform_zoff (image, lo_image,
 572                                                 lo_piv_par, valid_par,
 573                                                 piv_data, lo_piv_data,
 574                                                 &cum_residu,
 575                                                 &cum_residu_reached,
 576                                                 &sum_dxdy, &sum_dxdy_old,
 577                                                 isi_last, grid_last,
 578                                                 sweep_last, gpiv_par->verbose))
 579                 != NULL) {
 580                 g_warning ("GPIV_PIV_INTERR_IMG: %s", err_msg);
 581                 gpiv_free_img (lo_image);
 582                 gpiv_free_pivdata (lo_piv_data);
 583                 gpiv_free_pivdata (piv_data);
 584                 return NULL;
 585             }
 586
 587         } else {
 588
 589             /*
 590              * Apply results to output piv_data
 591              */
 592             gpiv_free_pivdata (piv_data);
 593             piv_data = gpiv_cp_pivdata (lo_piv_data);
 594             cum_residu_reached = TRUE;
 595         }
 596
 597         /*
 598          * Adapt grid.
 599          * If final grid has been reached, grid_last will be set.
 600          */
 601         if (lo_piv_par->int_shift > piv_par->int_shift
 602             && !sweep_stop) {
 603             GpivPivData *pd = NULL;
 604
 605             pd = gpiv_piv_gridadapt (image->header,
 606                                      piv_par,
 607                                      lo_piv_par,
 608                                      piv_data,
 609                                      sweep,
 610                                      &grid_last);
 611             gpiv_free_pivdata (piv_data);
 612             piv_data = gpiv_cp_pivdata (pd);
 613             gpiv_free_pivdata (pd);
 614
 615             gpiv_free_pivdata (lo_piv_data);
 616             lo_piv_data = gpiv_cp_pivdata (piv_data);
 617
 618             if (lo_piv_par->int_scheme == GPIV_IMG_DEFORM) {
 619                 gpiv_0_pivdata (lo_piv_data);
 620             }
 621
 622         } else {
 623             grid_last = TRUE;
 624         }
 625
 626         /*
 627          *  Adapt interrogation area size.
 628          *  If final size has been reached, isi_last will be set.
 629          */
 630         gpiv_piv_isizadapt (piv_par,
 631                             lo_piv_par,
 632                             &isi_last);
 633
 634         if (cum_residu_reached && isi_last && grid_last) {
 635             sweep_last = TRUE;
 636 /*             if (lo_piv_par->int_scheme == GPIV_ZERO_OFF_FORWARD */
 637 /*                 || lo_piv_par->int_scheme == GPIV_ZERO_OFF_CENTRAL) { */
 638 /*                 lo_piv_par->ifit = piv_par->ifit; */
 639 /*             } */
 640         }
 641         sweep++;
 642     }
 643
 644     /*
 645      * Writes existing fftw wisdom
 646      */
 647     gpiv_fwrite_fftw_wisdom (1);
 648     gpiv_fwrite_fftw_wisdom (-1);
 649     fftw_forget_wisdom();
 650     fftw_cleanup ();
 651
 652 #ifdef USE_MTRACE
 653     muntrace();
 654 #endif
 655     gpiv_free_img (lo_image);
 656     gpiv_free_pivdata (lo_piv_data);
 657
 658
 659     return piv_data;
 660 }
 661
 662
 663
 664 /*
 665  * Private piv interrogation functions
 666  */
 667 static GpivPivData *
 668 alloc_pivdata_gridgen (const GpivImagePar *image_par,
 669                        const GpivPivPar *piv_par
 670                        )
 671 /*-----------------------------------------------------------------------------
 672  * Determining the number of grid points, allocating memory for output data
 673  * and calculate locations of Interrogation Area's. If ad_int is enabled, a
 674  * course grid is started with and adapted for subsequent sweeps.
 675  ----------------------------------------------------------------------------*/
 676 {
 677     GpivPivData *piv_data = NULL;
 678     gchar *err_msg = NULL;
 679     GpivPivPar *lo_piv_par = NULL;
 680     guint nx, ny;
 681
 682
 683     if ((lo_piv_par = gpiv_piv_cp_parameters (piv_par)) == NULL) {
 684         gpiv_error ("alloc_pivdata_gridgen: failing gpiv_piv_cp_parameters");
 685     }
 686
 687     if (piv_par->int_size_i > piv_par->int_size_f
 688         && piv_par->int_shift < piv_par->int_size_i / GPIV_SHIFT_FACTOR) {
 689         lo_piv_par->int_shift = lo_piv_par->int_size_i / GPIV_SHIFT_FACTOR;
 690     }
 691
 692     if ((err_msg =
 693          gpiv_piv_count_pivdata_fromimage (image_par, lo_piv_par, &nx, &ny))
 694         != NULL) {
 695         warning_gpiv ("%s", err_msg);
 696         return NULL;
 697     }
 698
 699
 700     if ((piv_data = gpiv_piv_gridgen (nx, ny, image_par, lo_piv_par))
 701         == NULL) error_gpiv ("%s: %s", RCSID, err_msg);
 702
 703
 704     return piv_data;
 705 }
 706
 707
 708
 709 static void
 710 report_progress (const guint index_y,
 711                  const guint index_x,
 712                  gdouble *progress_value_prev,
 713                  const GpivPivData *piv_data,
 714                  const GpivPivPar *piv_par,
 715                  const guint sweep,
 716                  const gfloat cum_residu,
 717                  const GpivConsole *gpiv,
 718                  const Display *disp
 719                  )
 720 /*-----------------------------------------------------------------------------
 721  * Calculates progress of interrogation processing and other variables
 722  * and prints to message bar of the console if progress value has been changed
 723  */
 724 {
 725     gchar progress_string[GPIV_MAX_CHARS];
 726     gdouble progress_value = 100 * (index_y * piv_data->nx + index_x +1) /
 727         (piv_data->nx * piv_data->ny);
 728
 729
 730     if (progress_value != *progress_value_prev) {
 731         *progress_value_prev = progress_value;
 732         if (piv_par->int_scheme == GPIV_ZERO_OFF_FORWARD
 733             || piv_par->int_scheme == GPIV_ZERO_OFF_CENTRAL
 734             || piv_par->int_scheme == GPIV_IMG_DEFORM) {
 735             if (piv_par->int_size_i > piv_par->int_size_f /* substitutes ad_int == TRUE */
 736                 ) {
 737
 738                 g_snprintf (progress_string, GPIV_MAX_CHARS,
 739                             "Interrogating image #%d:"
 740                             " sweep #%d"
 741                             " size=%d"
 742                             " shift=%d"
 743                             " residu=%.3f"
 744                             ,
 745                             disp->id,
 746                             sweep,
 747                             piv_par->int_size_i,
 748                             piv_par->int_shift,
 749                             cum_residu
 750                             );
 751
 752             } else {
 753
 754                 g_snprintf (progress_string, GPIV_MAX_CHARS,
 755                             "Interrogating image #%d:"
 756                             " sweep #%d"
 757                             " shift=%d"
 758                             " residu=%.3f"
 759                             ,
 760                             disp->id,
 761                             sweep,
 762                             piv_par->int_shift,
 763                             cum_residu
 764                             );
 765
 766             }
 767
 768
 769         } else {
 770
 771             g_snprintf (progress_string, GPIV_MAX_CHARS,
 772                         "Interrogating image #%d: "
 773                         " shift=%d "
 774                         ,
 775                         disp->id,
 776                         piv_par/* _dest */ ->int_shift
 777                         );
 778         }
 779
 780         while (g_main_context_iteration (NULL, FALSE));
 781         gtk_progress_set_value (GTK_PROGRESS (gnome_appbar_get_progress
 782                                               (GNOME_APPBAR (gpiv->appbar))),
 783                                 progress_value);
 784
 785         gnome_appbar_set_status (GNOME_APPBAR(gpiv->appbar),
 786                                  progress_string);
 787     }
 788
 789 }
 790
 791
 792 static gchar *
 793 update_pivdata_imgdeform_zoff (const GpivImage *image,
 794                                GpivImage *lo_image,
 795                                const GpivPivPar *piv_par,
 796                                const GpivValidPar *valid_par,
 797                                GpivPivData *piv_data,
 798                                GpivPivData *lo_piv_data,
 799                                gfloat *cum_residu,
 800                                gboolean *cum_residu_reached,
 801                                gfloat *sum_dxdy,
 802                                gfloat *sum_dxdy_old,
 803                                gboolean isi_last,
 804                                gboolean grid_last,
 805                                gboolean sweep_last,
 806                                gboolean verbose
 807                                )
 808 /*-----------------------------------------------------------------------------
 809  * Validates and updates / renews pivdata and some other variables when image
 810  * deformation or zero-offset interrogation scheme is used
 811  */
 812 {
 813
 814     gchar *err_msg = NULL;
 815
 816
 817     /*
 818      * Test on outliers
 819      */
 820     if ((err_msg =
 821          gpiv_valid_errvec (lo_piv_data,
 822                             image,
 823                             piv_par,
 824                             valid_par,
 825                             TRUE))
 826         != NULL) {
 827         return err_msg;
 828     }
 829
 830
 831     if (piv_par->int_scheme == GPIV_IMG_DEFORM) {
 832
 833         /*
 834          * Update PIV estimators with those from the last interrogation
 835          * Resetting local PIV estimators for eventual next interrogation
 836          */
 837         if ((err_msg = gpiv_add_dxdy_pivdata (lo_piv_data, piv_data))
 838             != NULL) {
 839             return err_msg;
 840         }
 841
 842         if ((err_msg = gpiv_0_pivdata (lo_piv_data))
 843             != NULL) {
 844             return err_msg;
 845         }
 846
 847         /*
 848          * Deform image with updated PIV estimators.
 849          * First, copy local from original image.
 850          * Deform with newly updated PIV estimators.
 851          * Eventually write deformed image.
 852          */
 853
 854         if ((err_msg = gpiv_cp_img_data (image, lo_image))
 855             != NULL) {
 856             return err_msg;
 857         }
 858
 859         if ((err_msg = gpiv_imgproc_deform (lo_image, piv_data))
 860             != NULL) {
 861             return err_msg;
 862         }
 863 /* #define DEBUG */
 864 #ifdef DEBUG
 865         if (sweep_last && verbose) {
 866             printf ("\n");
 867             if ((err_msg =
 868                  gpiv_piv_write_deformed_image (lo_image))
 869                 != NULL) {
 870                 return err_msg;
 871             }
 872         }
 873 #endif /* DEBUG */
 874 /* #undef DEBUG */
 875
 876     } else {
 877
 878         /*
 879          * Renew PIV estimators with those from the last interrogation
 880          */
 881         if ((err_msg = gpiv_0_pivdata (piv_data))
 882             != NULL) {
 883             return err_msg;
 884         }
 885         if ((err_msg = gpiv_add_dxdy_pivdata (lo_piv_data, piv_data))
 886             != NULL) {
 887             return err_msg;
 888         }
 889     }
 890
 891     /*
 892      * Checking the relative cumulative residu for convergence
 893      * if final residu has been reached, cum_residu_reached will be set.
 894      */
 895     if (isi_last && grid_last) {
 896         *sum_dxdy_old = *sum_dxdy;
 897         *sum_dxdy = 0.0;
 898         gpiv_sum_dxdy_pivdata (piv_data, sum_dxdy);
 899         *cum_residu = fabsf ((*sum_dxdy - *sum_dxdy_old) /
 900                             ((gfloat)piv_data->nx * (gfloat)piv_data->ny));
 901         if (*cum_residu < GPIV_CUM_RESIDU_MIN) {
 902             *cum_residu_reached = TRUE;
 903         }
 904     }
 905
 906
 907     return err_msg;
 908 }
 909
 910
 911 #ifdef ENABLE_MPI
 912 static GpivPivData *
 913 exec_piv_mpi (GpivImage *image,
 914               GpivPivPar *piv_par,
 915               GpivValidPar *valid_par,
 916               GpivConsole *gpiv
 917               )
 918 /*-----------------------------------------------------------------------------
 919  * Interrogates an image on a distributed memory (Beowulf) cluster
 920  * using gpiv_rr from the gpivtools package
 921  */
 922 {
 923     GpivPivData *pd = NULL;
 924     FILE *fp = NULL;
 925     gchar *command = NULL;
 926     const gchar      *tmp_dir = g_get_tmp_dir ();
 927     const gchar      *user_name = g_get_user_name ();
 928     const gchar      *basename = "gpiv_mpi";
 929
 930     gchar progress_string[GPIV_MAX_CHARS];
 931
 932
 933 /*
 934  * prints to message console bar
 935  */
 936     g_snprintf (progress_string, GPIV_MAX_CHARS,
 937                 "Interrogating image #%d using %s %d gpiv_rr",
 938                 display_act->id, MPIRUN_CMD, gpiv_par->mpi_nodes);
 939     gnome_appbar_set_status (GNOME_APPBAR(gpiv->appbar),
 940                              progress_string);
 941 /*     gtk_progress_set_value (GTK_PULSE? (gnome_appbar_get_progress */
 942 /*                                           (GNOME_APPBAR (gpiv->appbar))),  */
 943 /*                             progress_value); */
 944
 945     /*
 946      * Save image at tmp/
 947      */
 948     gchar            *image_name =
 949         g_strdup_printf ("%s/%s/%s.png" ,tmp_dir,
 950                          user_name, basename);
 951     gchar            *pivdata_name =
 952         g_strdup_printf ("%s/%s/%s.piv" ,tmp_dir,
 953                          user_name, basename);
 954     gchar            *par_name =
 955         g_strdup_printf ("%s/%s/%s.par" ,tmp_dir,
 956                          user_name, basename);
 957
 958     if (gpiv_par->verbose) g_message ("exec_piv_mpi: image name = %s, piv name = %s",
 959                image_name, pivdata_name);
 960
 961     if ((fp = fopen (image_name, "w")) == NULL) {
 962         warning_gpiv ("exec_piv_mpi: failing fopen %s", image_name);
 963     }
 964     gpiv_write_png_image (fp, image, FALSE);
 965     fclose (fp);
 966
 967     /*
 968      * Run parallellised gpiv_rr
 969      */
 970     command =
 971         g_strdup_printf ("%s %d gpiv_rr -p \
 972 --cf %d --cl %d --cp %d \
 973 --rf %d --rl %d --rp %d \
 974 --ia_size_i %d --ia_size_f %d \
 975 --ia_shift %d --ischeme %d \
 976 --ifit %d --peak %d \
 977 --val_r %d --val_s %d --val_t %f",
 978                          MPIRUN_CMD, gpiv_par->mpi_nodes,
 979                          piv_par->col_start, piv_par->col_end, piv_par->pre_shift_col,
 980                          piv_par->row_start, piv_par->row_end, piv_par->pre_shift_row,
 981                          piv_par->int_size_i, piv_par->int_size_f,
 982                          piv_par->int_shift, piv_par->int_scheme,
 983                          piv_par->ifit, piv_par->peak,
 984                          valid_par->residu_type, valid_par->subst_type, valid_par->residu_max);
 985
 986     if (piv_par->spof_filter) {
 987         command = g_strconcat (command, " --spof", NULL);
 988     }
 989     if (gpiv_par->verbose) {
 990         command = g_strconcat (command, " -p", NULL);
 991     }
 992     command = g_strconcat (command, " ", image_name, NULL);
 993     if (gpiv_par->verbose) g_message ("exec_piv_mpi: command = %s", command);
 994
 995     if (system (command) != 0) {
 996         error_gpiv ("exec_piv_mpi: could not exec shell command:\n %s", command);
 997     }
 998     g_free (command);
 999
1000     /*
1001      * Obtain PIV-data from tmp/
1002      */
1003     pd = gpiv_fread_pivdata(pivdata_name);
1004
1005     /*
1006      * Cleanup image and data from tmp/
1007      */
1008     command = g_strdup_printf ("rm %s %s %s", image_name, pivdata_name, par_name);
1009     if (system (command) != 0) {
1010         error_gpiv ("exec_piv_mpi: could not exec shell command:\n %s", command);
1011     }
1012
1013     g_free (command);
1014     g_free (image_name);
1015     g_free (pivdata_name);
1016     return pd;
1017 }
1018 #endif /* ENABLE_MPI */