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