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