lib/piv_proc.c

   1 /* -*- Mode: C; indent-tabs-mode: nil; c-basic-offset: 4 c-style: "K&R" -*- */
   2
   3 /*
   4    libgpiv - library for Particle Image Velocimetry
   5
   6    Copyright (C) 2002, 2003, 2004, 2005, 2006, 2007, 2008
   7    Gerber van der Graaf
   8
   9    This file is part of libgpiv.
  10
  11    Libgpiv is free software; you can redistribute it and/or modify
  12    it under the terms of the GNU General Public License as published by
  13    the Free Software Foundation; either version 2, or (at your option)
  14    any later version.
  15
  16    This program is distributed in the hope that it will be useful,
  17    but WITHOUT ANY WARRANTY; without even the implied warranty of
  18    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  19    GNU General Public License for more details.
  20
  21    You should have received a copy of the GNU General Public License
  22    along with this program; if not, write to the Free Software Foundation,
  23    Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
  24
  25
  26
  27 ------------------------------------------------------------------------------
  28 FILENAME:                piv.c
  29 LIBRARY:                 libgpiv
  30 EXTERNAL FUNCTIONS:
  31                          gpiv_piv_count_pivdata_fromimage
  32                          gpiv_piv_select_int_point
  33                          gpiv_piv_interrogate_img
  34                          gpiv_piv_interrogate_ia
  35                          gpiv_piv_isizadapt
  36                          gpiv_piv_write_deformed_image
  37                          gpiv_piv_weight_kernel_1
  38                          gpiv_piv_weight_kernel_lin
  39                          gpiv_fread_fftw_wisdom
  40                          gpiv_fwrite_fftw_wisdom
  41
  42
  43 LAST MODIFICATION DATE:  $Id: piv.c,v 1.5 2008-09-25 13:19:53 gerber Exp $
  44 ---------------------------------------------------------------------------- */
  45
  46 #include <gpiv.h>
  47 #undef USE_MTRACE
  48 #ifdef USE_MTRACE
  49 #include <mcheck.h>
  50 #endif
  51 #include "my_utils.h"
  52
  53 #define FFTWISDOM "gpiv_fftwis" /* filename containg a string of wisdom for fft */
  54 #define FFTWISDOM_INV "gpiv_fftwis_inv" /* filename containg a string of wisdom for inverse fft */
  55
  56
  57
  58 /*-------------------------------------------------------------------------
  59  Levenbergh-Margardt parameter estimation constants */
  60 /* #define MA 4 */                      /* Number of parameters to be estimated by
  61                                    marquardt fit of cov peak */
  62 /* #define NPT 9 */                     /* Numbers of points to be fitted by mrqmin */
  63 /* #define SPREAD 0.001 */              /* multiplication factor for individual
  64                                    variances in measuremnts, i.e. covariance
  65                                    peak values */
  66 /* #undef NPT */
  67 /* #undef SPREAD */
  68
  69
  70
  71 #define NMIN_TO_INTERPOLATE 2
  72
  73 enum Position {
  74     LOWER,
  75     NEXT_LOWER,
  76     HIGHER,
  77     NEXT_HIGHER
  78 };
  79
  80 enum HorizontalPosition {
  81     WEST,
  82     WEST_WEST,
  83     EAST,
  84     EAST_EAST
  85 };
  86
  87 enum VerticalPosition {
  88     NORTH,
  89     NORTH_NORTH,
  90     SOUTH,
  91     SOUTH_SOUTH,
  92 };
  93
  94
  95 typedef struct __LoVarInterrogateImg LoVarII;
  96 /*
  97  * Local variables of gpiv_interrogate_img to distinguish from global or from
  98  * parameter list
  99  */
 100 struct __LoVarInterrogateImg
 101 {
 102     GpivImage *image;                /* local image that might be deformed */
 103     GpivPivPar *piv_par;             /* local piv parameters */
 104     GpivPivData *out_piv_data;       /* output piv data */
 105     GpivPivData *piv_data;           /* local piv data */
 106     gchar **err_msg_ar;              /* array of error messages for each thread */
 107 };
 108
 109 typedef struct __OpenMP Ompp;
 110 /*
 111  * OpenMP parallelization data
 112  */
 113 struct __OpenMP
 114 {
 115     int max_nr_thr;                  /* maximum number of threads possible */
 116     int nr_thr;                      /* current number of threads available to
 117                                         the program */
 118     int thr_id;                      /* thread identity */
 119 };
 120
 121 typedef struct __MessagePassingInterface Mpi;
 122 /*
 123  * Message Passing Interface data
 124  */
 125 struct __MessagePassingInterface
 126 {
 127     GpivPivData *piv_data;           /* scattered and gathered piv data to be
 128                                         used for parallel processing using MPI */
 129     /* gboolean scatter; */
 130     gint count;
 131     gint *counts;
 132     gint *displs;
 133     gint nprocs;
 134     gint rank;
 135 };
 136
 137 /*
 138  * Local functions prototypes
 139  */
 140 static Mpi *
 141 alloc_mpi(GpivPivData *piv_data);
 142
 143 static void
 144 free_mpi(Mpi *mpi);
 145
 146 static GpivPivData *
 147 alloc_pivdata_gridgen (const GpivImagePar *image_par,
 148                        const GpivPivPar *piv_par
 149                        );
 150
 151 static void
 152 report_progress (gint *progress_old,
 153                  gint index_y,
 154                  gint index_x,
 155                  GpivPivData *piv_data,
 156                  GpivPivPar *piv_par,
 157                  gint sweep,
 158                  gfloat cum_residu
 159                  );
 160
 161 static gboolean
 162 assign_img2intarr (gint ipoint_x,
 163                    gint ipoint_y,
 164                    guint16 **img_1,
 165                    guint16 **img_2,
 166                    gint int_size_f,
 167                    gint int_size_i,
 168                    gfloat **int_area_1,
 169                    gfloat **int_area_2,
 170                    gint pre_shift_row,
 171                    gint pre_shift_col,
 172                    gint nrows,
 173                    gint ncolumns,
 174                    gint int_size_0
 175                    );
 176
 177 static gboolean
 178 assign_img2intarr_central (gint ipoint_x,
 179                            gint ipoint_y,
 180                            guint16 **img_1,
 181                            guint16 **img_2,
 182                            gint int_size_f,
 183                            gint int_size_i,
 184                            gfloat **int_area_1,
 185                            gfloat **int_area_2,
 186                            gint pre_shift_row,
 187                            gint pre_shift_col,
 188                            gint nrows,
 189                            gint ncolumns,
 190                            gint int_size_0
 191                            );
 192
 193 static gboolean
 194 assign_img2intarr_forward (gint ipoint_x,
 195                            gint  ipoint_y,
 196                            guint16 **img_1,
 197                            guint16 **img_2,
 198                            gint int_size_f,
 199                            gint int_size_i,
 200                            gfloat **int_area_1,
 201                            gfloat **int_area_2,
 202                            gint pre_shift_row,
 203                            gint pre_shift_col,
 204                            gint nrows,
 205                            gint ncolumns,
 206                            gint int_size_0
 207                            );
 208
 209 static float
 210 int_mean_old (guint16 **img,
 211               int int_size,
 212               int int_size_i,
 213               int ipoint_x,
 214               int ipoint_y
 215               );
 216
 217 static gfloat
 218 int_mean (gfloat **int_area,
 219           gint int_size
 220           );
 221
 222 static gfloat
 223 int_range (gfloat **int_area,
 224            gint int_size
 225            );
 226
 227 static gboolean
 228 int_const (gfloat **int_area,
 229            const guint int_size
 230            );
 231
 232 static void
 233 cov_min_max (GpivCov *cov
 234              );
 235
 236 static void
 237 search_top (GpivCov *cov,
 238             gint peak_act,
 239             gint x_corr,
 240             gint sweep,
 241             gint i_skip_act,
 242             gint j_skip_act,
 243             float *z_max,
 244             long *i_max,
 245             long *j_max
 246             );
 247
 248 static char *
 249 cov_subtop (float **z,
 250             long *i_max,
 251             long *j_max,
 252             float *epsi_x,
 253             float *epsi_y,
 254             int ifit,
 255             int peak_act
 256             );
 257
 258 static int
 259 cov_top (GpivPivPar piv_par,
 260          GpivPivData * piv_data,
 261          int index_y,
 262          int index_x,
 263          GpivCov *cov,
 264          int x_corr,
 265          int ifit,
 266          int sweep,
 267          int last_sweep,
 268          int peak,
 269          int peak_act,
 270          int pre_shift_row_act,
 271          int pre_shift_col_act,
 272          int i_skip_act,
 273          int j_skip_act,
 274          gboolean *flag_subtop
 275          );
 276
 277 static
 278 void pack_cov (float **covariance,
 279                GpivCov *cov,
 280                int int_size_0
 281                );
 282
 283 static void
 284 piv_weight_kernel_lin (const guint int_size_0,
 285                        GpivCov *w_k
 286                        );
 287
 288 static void
 289 weight_cov (GpivCov *cov,
 290             GpivCov *w_k
 291             );
 292
 293 static gchar *
 294 filter_cov_spof (fftw_complex *a,
 295                  fftw_complex *b,
 296                  gint m,
 297                  gint n
 298                  );
 299
 300 static gchar *
 301 cova (const gboolean spof_filter,
 302       GpivFt *ft
 303       );
 304
 305 static gchar *
 306 ia_weight_gauss (gint int_size,
 307                  float **int_area
 308                  );
 309
 310 /*
 311  * Origined from piv_speed
 312  */
 313 static void
 314 nearest_point (gint i,
 315                gfloat x,
 316                gfloat point_x,
 317                gfloat *min,
 318                gint *index,
 319                gboolean *exist
 320                );
 321
 322 static gboolean
 323 nearest_index (enum Position pos,
 324                gint vlength,
 325                gfloat *src_point,
 326                gfloat dest_point,
 327                gint *index
 328                );
 329
 330 static gdouble
 331 bilinear_interpol (gdouble alpha_hor,
 332                    gdouble alpha_vert,
 333                    gdouble src_dx_nw,
 334                    gdouble src_dx_ne,
 335                    gdouble src_dx_sw,
 336                    gdouble src_dx_se
 337                    );
 338
 339 static void *
 340 intpol_facts (gfloat *src_point,
 341               gfloat *dest_point,
 342               gint src_vlength,
 343               gint dest_vlength,
 344               gint *index_l,
 345               gint *index_h,
 346               gint *index_ll,
 347               gint *index_hh,
 348               double *alpha
 349               );
 350
 351 static void
 352 dxdy_at_new_grid_block (const GpivPivData *piv_data_src,
 353                         GpivPivData *piv_data_dest,
 354                         gint expansion_factor,
 355                         gint smooth_window
 356                         );
 357
 358 static gchar *
 359 update_pivdata_imgdeform_zoff (const GpivImage *image,
 360                                GpivImage *lo_image,
 361                                const GpivPivPar *piv_par,
 362                                const GpivValidPar *valid_par,
 363                                GpivPivData *piv_data,
 364                                GpivPivData *lo_piv_data,
 365                                GpivFt *ft,
 366                                gfloat *cum_residu,
 367                                gboolean *cum_residu_reached,
 368                                gfloat *sum_dxdy,
 369                                gfloat *sum_dxdy_old,
 370                                gboolean isi_last,
 371                                gboolean grid_last,
 372                                gboolean sweep_last,
 373                                gboolean verbose
 374                                );
 375
 376 static LoVarII *
 377 piv_interrogate_img__init (const GpivImage *image,
 378                            const GpivPivPar *piv_par,
 379                            const GpivValidPar *valid_par,
 380                            Ompp *ompp,
 381                            Mpi *mpi,
 382                            gboolean *sweep_last,
 383                            const gboolean verbose
 384                            );
 385
 386 static gchar *
 387 piv_interrogate_img__finish (LoVarII *lo
 388                              );
 389
 390 static void
 391 piv_interrogate_img__err (LoVarII *lo,
 392                           GpivFt **ft,
 393                           const int nr_thr);
 394
 395 /*
 396  * Some MPI routines
 397  */
 398 #ifdef ENABLE_MPI
 399 static GpivPivData *
 400 piv_interrogate_img__scatterv_pivdata (GpivPivData *piv_data);
 401
 402 static GpivPivData *
 403 piv_interrogate_img__gatherv_pivdata (GpivPivData *lo_piv_data,
 404                                       GpivPivData *piv_data);
 405
 406 guint
 407 substr_noremain (guint n,
 408                  guint m);
 409
 410 #endif /* ENABLE_MPI */
 411
 412 /*
 413  * Public functions
 414  */
 415
 416 gchar *
 417 gpiv_piv_count_pivdata_fromimage (const GpivImagePar *image_par,
 418                                   const GpivPivPar *piv_par,
 419                                   guint *nx,
 420                                   guint *ny
 421                                   )
 422 /*-----------------------------------------------------------------------------
 423  *     Calculates the number of interrogation areas from the image sizes,
 424  *     pre-shift and area of interest
 425  *     NULL on success or error message on failure
 426  *---------------------------------------------------------------------------*/
 427 {
 428     gchar *err_msg = NULL;
 429     int row, column, row_1, column_1,
 430         pre_shift_row_max, pre_shift_col_max, count_x = 0, count_y = 0;
 431     int row_max, row_min, column_max, column_min;
 432
 433     int ncolumns = image_par->ncolumns;
 434     int nrows = image_par->nrows;
 435
 436     int int_geo = piv_par->int_geo;
 437     int row_start = piv_par->row_start;
 438     int row_end = piv_par->row_end;
 439     int col_start = piv_par->col_start;
 440     int col_end = piv_par->col_end;
 441     int int_line_col = piv_par->int_line_col;
 442     int int_line_col_start = piv_par->int_line_col_start;
 443     int int_line_col_end = piv_par->int_line_col_end;
 444     int int_line_row = piv_par->int_line_row;
 445     int int_line_row_start = piv_par->int_line_row_start;
 446     int int_line_row_end = piv_par->int_line_row_end;
 447     int int_point_col = piv_par->int_point_col;
 448     int int_point_row = piv_par->int_point_row;
 449     int int_size_f = piv_par->int_size_f;
 450     int int_size_i = piv_par->int_size_i;
 451     int int_shift = piv_par->int_shift;
 452     int pre_shift_row = piv_par->pre_shift_row;
 453     int pre_shift_col = piv_par->pre_shift_col;
 454
 455
 456     *nx = 0;
 457     *ny = 0;
 458
 459
 460     row_min = gpiv_min(-int_size_f / 2 + 1,
 461                        pre_shift_row - int_size_i / 2 + 1);
 462     column_min = gpiv_max(-int_size_f / 2 + 1,
 463                           pre_shift_col - int_size_i / 2 + 1);
 464     row_max = gpiv_max(int_size_f / 2, pre_shift_row + int_size_i / 2);
 465     column_max = gpiv_max(int_size_f / 2, pre_shift_col + int_size_i / 2);
 466
 467
 468     if (int_geo == GPIV_POINT) {
 469         *nx = 1;
 470         *ny = 1;
 471
 472
 473 /*
 474  * Counts number of Interrrogation Area for a single row
 475  */
 476     } else if (int_geo == GPIV_LINE_R) {
 477         if ((int_size_f - int_size_i) / 2 + pre_shift_col < 0) {
 478             column_1 = int_line_col_start -
 479                 ((int_size_f - int_size_i) / 2 +
 480                  pre_shift_col) + int_size_f / 2 - 1;
 481         } else {
 482                 column_1 = int_line_col_start + int_size_f / 2 - 1;
 483         }
 484
 485         for (column = column_1; column <= int_line_col_end - column_max;
 486              column += int_shift) {
 487             count_x++;
 488         }
 489
 490         *nx = count_x;
 491         *ny = 1;
 492
 493 /*
 494  * Counts number of Interrrogation Area for a single column
 495  */
 496     } else if (int_geo == GPIV_LINE_C) {
 497         if ((int_size_f - int_size_i) / 2 + pre_shift_row < 0) {
 498             row_1 = int_line_row_start -
 499                 ((int_size_f - int_size_i) / 2 +
 500                      pre_shift_row) + int_size_f / 2 - 1;
 501         } else {
 502             row_1 = int_line_row_start + int_size_f / 2 - 1;
 503         }
 504
 505         for (row = row_1; row <= int_line_row_end - row_max;
 506              row += int_shift) {
 507             count_y++;
 508         }
 509
 510         *ny = count_y;
 511         *nx = 1;
 512
 513
 514 /*
 515  * Counts number of Interrrogation Area for a Area Of Interest
 516  */
 517     } else if (int_geo == GPIV_AOI) {
 518         if ((int_size_f - int_size_i) / 2 + pre_shift_row < 0) {
 519             row_1 =
 520                 row_start - ((int_size_f - int_size_i) / 2 +
 521                              pre_shift_row) + int_size_f / 2 - 1;
 522         } else {
 523             row_1 = row_start + int_size_f / 2 - 1;
 524         }
 525         if ((int_size_f - int_size_i) / 2 + pre_shift_col < 0) {
 526             column_1 =
 527                 col_start - ((int_size_f - int_size_i) / 2 +
 528                              pre_shift_col) + int_size_f / 2 - 1;
 529         } else {
 530             column_1 = col_start + int_size_f / 2 - 1;
 531         }
 532
 533
 534         pre_shift_col_max = gpiv_max (pre_shift_col, 0);
 535         column_max =
 536             gpiv_max(int_size_f / 2, pre_shift_col + int_size_i / 2);
 537         pre_shift_row_max = gpiv_max (pre_shift_row, 0);
 538         row_max = gpiv_max (int_size_f / 2, pre_shift_row + int_size_i / 2);
 539
 540
 541         for (row = row_1; row + row_max <= row_end; row += int_shift) {
 542             for (column = column_1; column + column_max <= col_end;
 543                  column += int_shift) {
 544                 count_x++;
 545             }
 546             if (count_x > *nx)
 547                 *nx = count_x;
 548             count_x = 0;
 549             count_y++;
 550         }
 551         if (count_y > *ny)
 552           *ny = count_y;
 553     } else {
 554         err_msg = "gpiv_piv_count_pivdata_fromimage: should not arrive here";
 555         gpiv_warning ("%s", err_msg);
 556         return err_msg;
 557     }
 558
 559     if (*nx == 0 || *ny == 0) {
 560         err_msg = "gpiv_piv_count_pivdata_fromimage: line or AOI too small: nx=0 ny=0";
 561         gpiv_warning("gpiv_piv_count_pivdata_fromimage: line or AOI too small: nx = %d ny = %d\n",
 562                      *nx, *ny);
 563         return err_msg;
 564     }
 565
 566 #ifdef DEBUG
 567     g_message ("gpiv_piv_count_pivdata_fromimage:: 2 nx = %d, ny = %d", *nx, *ny);
 568 #endif
 569     return err_msg;
 570 }
 571
 572
 573
 574 GpivPivData *
 575 gpiv_piv_interrogate_img (const GpivImage *image,
 576                           const GpivPivPar *piv_par,
 577                           const GpivValidPar *valid_par,
 578                           const gboolean verbose
 579                           )
 580 /* ----------------------------------------------------------------------------
 581  * PIV interrogation of an image pair at an entire grid or single point
 582  *
 583  *     @param[in] image           image containing data and header info
 584  *     @param[in] piv_par         image evaluation parameters
 585  *     @param[in] valid_par       structure of PIV data validation parameters
 586  *     @param[out] verbose        prints progress of interrogation to stdout
 587  *     @return                    GpivPivData containing PIV estimators on succes
 588  *                                or NULL on failure
 589  */
 590 /*---------------------------------------------------------------------------*/
 591 {
 592     GpivPivData *piv_data = NULL;       /* piv data to be returned to caller */
 593     gboolean error = FALSE;             /* error flag */
 594     long int index_x = 0, index_y = 0;  /* array indices */
 595
 596     /*
 597      * Local variables with prefix lo_ to distinguish from global or from
 598      * parameter list
 599      */
 600     LoVarII *lo = NULL;                 /* Contains local variables */
 601     GpivFt **ft = NULL;                 /* arrays of stuctures to perform FFT */
 602     guint sweep = 1;                    /* itaration counter */
 603     gboolean grid_last = FALSE;         /* flag if final grid refinement has been
 604                                            reached */
 605     gboolean isi_last = FALSE;          /* flag if final interrogation area shift
 606                                            has been reached */
 607     gboolean cum_residu_reached = FALSE;/* flag if max. cumulative residu has
 608                                            been reached */
 609     gboolean sweep_last = FALSE;        /* perform the last iteration sweep */
 610     gboolean sweep_stop = FALSE;        /* stop the current iteration at the end */
 611     gfloat sum_dxdy = 0.0, sum_dxdy_old = 0.0;  /* */
 612     gfloat cum_residu = 914.6;          /* initial, large, arbitrary cumulative
 613                                            residu */
 614     guint progress_old = 0;             /* for monitoring calculation progress */
 615     gint i = 0;                         /* a counter */
 616
 617     /*
 618      * Variables for OMP and MPI
 619      */
 620     Ompp *ompp = g_new0(Ompp, 1);
 621     int max_nr_thr = -1;
 622     int nr_thr = -1;
 623     int thr_id = -1;
 624     Mpi *mpi = g_new0(Mpi,1);
 625
 626
 627     /*
 628      * Initializing all variables
 629      */
 630     lo = piv_interrogate_img__init(image, piv_par, valid_par, ompp, mpi,
 631                                    &sweep_last, verbose);
 632     max_nr_thr = ompp->max_nr_thr;
 633
 634     /*
 635      * Interrogates at single a point or at a grid, using advanced schemes
 636      */
 637     while (sweep <= GPIV_MAX_PIV_SWEEP && !sweep_stop) {
 638
 639         /*
 640          * Memory allocation of interrogation area's, covariance and FFT plans.
 641          */
 642         ft = g_new0 (GpivFt, ompp->max_nr_thr);
 643         for (i = 0; i < ompp->max_nr_thr; i++) {
 644             ft[i] = gpiv_alloc_ft (GPIV_ZEROPAD_FACT * lo->piv_par->int_size_i);
 645         }
 646
 647         /*
 648          * Interrogates a single interrogation area
 649          */
 650         if (lo->piv_par->int_geo == GPIV_POINT) {
 651
 652             if ((lo->err_msg_ar[0] =
 653                  gpiv_piv_interrogate_ia(0, 0, lo->image, lo->piv_par, sweep,
 654                                          sweep_last, lo->piv_data, ft[0]))
 655                 != NULL) {
 656                 piv_interrogate_img__err(lo, ft, 0);
 657                 return NULL;
 658             }
 659
 660         } else {
 661
 662             /*
 663              * Interrogates at a rectangular grid of points within the Area Of
 664              * Interest of the image
 665              */
 666
 667 #ifdef ENABLE_MPI
 668             /*
 669              * Scatter the PIV data over the rows to the different nodes.
 670              */
 671             lo->piv_data = piv_interrogate_img__scatterv_pivdata(lo->piv_data);
 672 #endif /* ENABLE_MPI */
 673
 674             /* grid_last, isi_last, cum_residu_reached, sum_dxdy, sum_dxdy_old */
 675 #pragma omp parallel default(none)                                             \
 676     private(thr_id, index_x)                                              \
 677     shared(nr_thr, max_nr_thr, index_y, error, sweep, sweep_stop,    \
 678            sweep_last, lo, cum_residu, progress_old, ft)
 679             {
 680 #ifdef ENABLE_OMP
 681                 nr_thr = omp_get_num_threads();
 682                 thr_id = omp_get_thread_num();
 683 #else /* ENABLE_OMP */
 684                 nr_thr = 1;
 685                 thr_id = 0;
 686 #endif /* ENABLE_OMP */
 687
 688 #pragma omp for schedule(dynamic, 1)
 689                 for (index_y = 0; index_y < lo->piv_data->ny; index_y++) {
 690                     for (index_x = 0; index_x < lo->piv_data->nx; index_x++) {
 691
 692 #ifdef DEBUG
 693                         if (rank == 0) {
 694                             g_message ("gpiv_piv_interrogate_img:: MPI: rank=%d",
 695                                        mpi->rank);
 696                         }
 697                         g_message("gpiv_piv_interrogate_img:: OMP: nr_thr = %d thr_id = %d i=%d j=%d",
 698                                   nr_thr, thr_id, index_y, index_x);
 699 #endif /* DEBUG */
 700
 701                         /*
 702                          * Interrogates a single interrogation area at the grid.
 703                          */
 704                         if ((lo->err_msg_ar[thr_id] =
 705                              gpiv_piv_interrogate_ia(index_y, index_x,
 706                                                      lo->image, lo->piv_par,
 707                                                      sweep, sweep_last,
 708                                                      lo->piv_data,
 709                                                      ft[thr_id]))
 710                             != NULL) {
 711                             error = TRUE;
 712 #ifdef ENABLE_MPI
 713                             MPI_Finalize();
 714 #endif
 715                         }
 716
 717                         /*
 718                          * Printing the progress of processing
 719                          */
 720                         if (verbose) {
 721 #pragma omp critical
 722                             report_progress(&progress_old, index_y, index_x,
 723                                             lo->piv_data, lo->piv_par, sweep,
 724                                             cum_residu);
 725                         }
 726                     }
 727                 }
 728             }
 729
 730 #ifdef ENABLE_MPI
 731             /*
 732              * Gather the scattered PIV data
 733              * and broadcasts the entire array to all nodes.
 734              */
 735             lo->piv_data =
 736                 piv_interrogate_img__gatherv_pivdata(lo->piv_data, lo->out_piv_data);
 737             gpiv_piv_mpi_bcast_pivdata (lo->piv_data);
 738 #endif
 739         }
 740
 741         /*
 742          * Error handling
 743          */
 744         if (error) {
 745             piv_interrogate_img__err(lo, ft, nr_thr);
 746             return NULL;
 747         }
 748
 749         if (sweep_last) {
 750             sweep_stop = TRUE;
 751         }
 752
 753         if (lo->piv_par->int_scheme == GPIV_IMG_DEFORM
 754             || lo->piv_par->int_scheme == GPIV_ZERO_OFF_FORWARD
 755             || lo->piv_par->int_scheme == GPIV_ZERO_OFF_CENTRAL) {
 756
 757             if((lo->err_msg_ar[thr_id] =
 758                 update_pivdata_imgdeform_zoff(image, lo->image, lo->piv_par,
 759                                               valid_par, lo->out_piv_data, lo->piv_data,
 760                                               ft[0], &cum_residu,
 761                                               &cum_residu_reached, &sum_dxdy,
 762                                               &sum_dxdy_old, isi_last,
 763                                               grid_last, sweep_last, verbose))
 764                != NULL) {
 765                 piv_interrogate_img__err(lo, ft, nr_thr);
 766                 return NULL;
 767             }
 768
 769         } else {
 770
 771             /*
 772              * Apply results to output piv_data
 773              */
 774
 775             gpiv_free_pivdata(lo->out_piv_data);
 776             lo->out_piv_data = gpiv_cp_pivdata(lo->piv_data);
 777             cum_residu_reached = TRUE;
 778         }
 779
 780         /*
 781          * De-allocating memory: other (smaller) sizes are eventually needed
 782          * for a next iteration sweep
 783          */
 784         for (i = 0; i < max_nr_thr; i++) {
 785             gpiv_free_ft(ft[i]);
 786         }
 787         g_free(ft);
 788
 789         /*
 790          * Adapt grid.
 791          * If final grid has been reached, grid_last will be set.
 792          */
 793         if (lo->piv_par->int_shift > piv_par->int_shift
 794             && !sweep_stop) {
 795             GpivPivData *pd = NULL;
 796
 797             pd = gpiv_piv_gridadapt(image->header, piv_par,  lo->piv_par,
 798                                     lo->out_piv_data, sweep, &grid_last);
 799             gpiv_free_pivdata(lo->out_piv_data);
 800             lo->out_piv_data = gpiv_cp_pivdata(pd);
 801             gpiv_free_pivdata(pd);
 802             gpiv_free_pivdata(lo->piv_data);
 803             lo->piv_data = gpiv_cp_pivdata(lo->out_piv_data);
 804
 805             if (lo->piv_par->int_scheme == GPIV_IMG_DEFORM) {
 806                 gpiv_0_pivdata(lo->piv_data);
 807             }
 808
 809         } else {
 810             grid_last = TRUE;
 811         }
 812
 813         /*
 814          *  Adapt interrogation area size.
 815          *  If final size has been reached, isi_last will be set.
 816          */
 817         gpiv_piv_isizadapt(piv_par, lo->piv_par, &isi_last);
 818
 819         /*
 820          * Test if all conditions have been reached
 821          */
 822         if (cum_residu_reached && isi_last && grid_last) {
 823             sweep_last = TRUE;
 824         }
 825
 826         sweep++;
 827     }
 828
 829     if (verbose) printf("\n");
 830
 831     /*
 832      * clean-up allocated memory, save existing fftw wisdom
 833      * and returns resulting PIV data to caller
 834      */
 835     piv_data = gpiv_cp_pivdata(lo->out_piv_data);
 836     piv_interrogate_img__finish(lo);
 837     return piv_data;
 838 }
 839
 840
 841
 842 gchar *
 843 gpiv_piv_interrogate_ia (const guint index_y,
 844                          const guint index_x,
 845                          const GpivImage *image,
 846                          const GpivPivPar *piv_par,
 847                          const guint sweep,
 848                          const guint last_sweep,
 849                          GpivPivData *piv_data,
 850                          GpivFt *ft
 851                          )
 852 /**----------------------------------------------------------------------------
 853  * Interrogates a single Interrogation Area
 854  */
 855 {
 856     gchar *err_msg = NULL;
 857
 858     guint ncolumns = image->header->ncolumns;
 859     guint nrows = image->header->nrows;
 860     gboolean x_corr = image->header->x_corr;
 861
 862     guint ifit = piv_par->ifit;
 863     guint int_size_f = piv_par->int_size_f;
 864     guint int_size_i = piv_par->int_size_i;
 865     gint peak = piv_par->peak;
 866     int pre_shift_row = piv_par->pre_shift_row;
 867     int pre_shift_col = piv_par->pre_shift_col;
 868     enum GpivIntScheme int_scheme = piv_par->int_scheme;
 869
 870     int return_val;
 871     int idum = gpiv_max((int_size_i - int_size_f) / 2, 0);
 872     int m = 0, n = 0;
 873     float intreg1_mean = 0.0, intreg2_mean = 0.0;
 874 /*     BUGFIX: gpiv_piv_interrogate_ia: disabled normalization I.A */
 875 #ifdef NORM_AI
 876     float intreg1_range = 0.0, intreg2_range = 0.0;
 877     float norm_fact = 0.0;
 878    guint img_top = (1 << image->header->depth) - 1;
 879 #endif
 880     int ipoint_x;
 881     int ipoint_y;
 882
 883     int pre_shift_row_act = 0, pre_shift_col_act = 0;
 884     int peak_act = 0, i_skip_act = 0, j_skip_act = 0;
 885
 886     gboolean flag_subtop = FALSE, flag_intar0 = FALSE, flag_accept = FALSE;
 887 /*
 888  * Interrogation area with zero padding
 889  */
 890     GpivCov *w_k = NULL;                /* covariance weighting kernel */
 891     int int_size_0 = GPIV_ZEROPAD_FACT * int_size_i;
 892     int j, k;
 893
 894
 895 /*
 896  * Checking for memory allocation of input variables
 897  */
 898     if ((err_msg = gpiv_check_alloc_ft(ft)) != NULL) {
 899         return (err_msg);
 900     }
 901
 902     if ((err_msg = gpiv_check_alloc_pivdata(piv_data)) != NULL) {
 903         return (err_msg);
 904     }
 905
 906     ipoint_x = (int) piv_data->point_x[index_y][index_x];
 907     ipoint_y = (int) piv_data->point_y[index_y][index_x];
 908
 909 /*
 910  * uses  piv values from previous estimation as pre-shifts and
 911  * searches closest Int. Area
 912  */
 913     if (int_scheme == GPIV_ZERO_OFF_FORWARD
 914         || int_scheme == GPIV_ZERO_OFF_CENTRAL) {
 915         pre_shift_col_act = piv_data->dx[index_y][index_x] + pre_shift_col;
 916         pre_shift_row_act = piv_data->dy[index_y][index_x] + pre_shift_row;
 917         piv_data->dx[index_y][index_x] = 0.0;
 918         piv_data->dy[index_y][index_x] = 0.0;
 919     } else {
 920         pre_shift_col_act = pre_shift_col;
 921         pre_shift_row_act = pre_shift_row;
 922     }
 923
 924     peak_act = peak;
 925
 926 /*
 927  * Assigning image data to the interrogation area's
 928  */
 929     if (int_scheme == GPIV_ZERO_OFF_CENTRAL ) {
 930         flag_accept =
 931             assign_img2intarr_central (ipoint_x, ipoint_y,
 932                                        image->frame1, image->frame2,
 933                                        int_size_f, int_size_i,
 934                                        ft->intreg1, ft->intreg2,
 935                                        pre_shift_row_act, pre_shift_col_act,
 936                                        nrows, ncolumns, int_size_0);
 937
 938     } else if (int_scheme == GPIV_ZERO_OFF_FORWARD) {
 939         flag_accept =
 940             assign_img2intarr_central (ipoint_x, ipoint_y,
 941                                        image->frame1, image->frame2,
 942                                        int_size_f, int_size_i,
 943                                        ft->intreg1, ft->intreg2,
 944                                        pre_shift_row_act, pre_shift_col_act,
 945                                        nrows, ncolumns, int_size_0);
 946     } else {
 947         flag_accept =
 948             assign_img2intarr (ipoint_x, ipoint_y,
 949                                image->frame1, image->frame2,
 950                                int_size_f, int_size_i,
 951                                ft->intreg1, ft->intreg2,
 952                                pre_shift_row_act, pre_shift_col_act,
 953                                nrows, ncolumns, int_size_0);
 954     }
 955
 956     if (flag_accept) {
 957
 958 /*
 959  * Weighting Interrogation Area with Gaussian function
 960  */
 961        if (piv_par->gauss_weight_ia) {
 962             if ((err_msg = ia_weight_gauss (int_size_f, ft->intreg1))
 963                 != NULL) {
 964                 return (err_msg);
 965             }
 966
 967             if ((err_msg = ia_weight_gauss (int_size_i, ft->intreg2))
 968                 != NULL) {
 969                 return (err_msg);
 970             }
 971         }
 972
 973 /*
 974  * The mean value of the image data within an interrogation area will be
 975  * subtracted from the data
 976  * BUXFIX: test on differences in estimator!
 977  */
 978         intreg1_mean = int_mean (ft->intreg1, int_size_f);
 979         intreg2_mean = int_mean (ft->intreg2, int_size_i);
 980 #ifdef NORM_AI
 981         intreg1_range = int_range (ft->intreg1, int_size_f);
 982         intreg2_range = int_range (ft->intreg2, int_size_i);
 983 #endif
 984
 985 /*
 986  * BUGFIX: this might be substituted by counting the number of particles within
 987  * Int. Area, as done in PTV
 988  */
 989         if (intreg1_mean == 0.0 || intreg2_mean == 0.0
 990             || int_const (ft->intreg1, int_size_f)
 991             || int_const (ft->intreg2, int_size_i)
 992             ) {
 993 /*             err_msg = "gpiv_piv_interrogate_ia: intreg1/2_mean = 0.0"; */
 994             flag_intar0 = TRUE;
 995 /*             return err_msg; */
 996
 997             }
 998
 999 #ifdef NORM_AI
1000         norm_fact = (gfloat) img_top / intreg1_range;
1001 #endif
1002         for (m = 0; m < int_size_f; m++) {
1003             for (n = 0; n < int_size_f; n++) {
1004                 ft->intreg1[m + idum ][n + idum ] -= intreg1_mean;
1005 #ifdef NORM_AI
1006                 ft->intreg1[m + idum ][n + idum ] *= norm_fact;
1007 #endif
1008             }
1009         }
1010
1011 #ifdef NORM_AI
1012         norm_fact = (gfloat) img_top / intreg2_range;
1013 #endif
1014         for (m = 0; m < int_size_i; m++) {
1015             for (n = 0; n < int_size_i; n++) {
1016                 ft->intreg2[m][n] -= intreg2_mean;
1017 #ifdef NORM_AI
1018                 ft->intreg2[m][n] *= norm_fact;
1019 #endif
1020             }
1021         }
1022
1023 /*
1024  * Calculate covariance function
1025  */
1026         if (!flag_intar0) {
1027             if ((err_msg = cova (piv_par->spof_filter, ft))
1028                 != NULL) {
1029                 gpiv_warning("%s", err_msg);
1030                 return err_msg;
1031             }
1032
1033 /*
1034  * Weighting covariance data with weight kernel
1035  */
1036             if (piv_par->int_scheme == GPIV_LK_WEIGHT) {
1037                 w_k = gpiv_alloc_cov (int_size_0, image->header->x_corr);
1038                 piv_weight_kernel_lin (int_size_0, w_k);
1039                 weight_cov (ft->cov, w_k);
1040                 gpiv_free_cov (w_k);
1041             }
1042
1043 /*
1044  * Searching maximum peak in covariance function
1045  */
1046             if ((return_val =
1047                  cov_top (*piv_par, piv_data, index_y, index_x, ft->cov, x_corr,
1048                           ifit, sweep, last_sweep, peak, peak_act,
1049                           pre_shift_row_act, pre_shift_col_act,
1050                           i_skip_act, j_skip_act, &flag_subtop)) != 0) {
1051                 err_msg = "gpiv_piv_interrogate_ia: Unable to call cov_top";
1052                 gpiv_warning("%s", err_msg);
1053                 return err_msg;
1054             }
1055
1056 /*
1057  * Writing values to piv_data
1058  */
1059             piv_data->dx[index_y][index_x] =
1060                 (double) pre_shift_col_act +
1061                 (double) ft->cov->top_x + ft->cov->subtop_x;
1062
1063             piv_data->dy[index_y][index_x] =
1064                 (double) pre_shift_row_act +
1065                 (double) ft->cov->top_y + ft->cov->subtop_y;
1066
1067 /*
1068  * Disabled as outliers are tested after each iteration
1069  */
1070 /*             if (last_sweep == 1) { */
1071                 piv_data->snr[index_y][index_x] = ft->cov->snr;
1072                 piv_data->peak_no[index_y][index_x] = peak_act;
1073 /*             } */
1074
1075         }
1076 /*
1077  * Check on validity of data
1078  */
1079         if (isnan(piv_data->dx[index_y][index_x]) != 0
1080             || isnan(piv_data->dy[index_y][index_x]) != 0
1081             || isnan(piv_data->snr[index_y][index_x]) != 0
1082             || flag_subtop
1083             || flag_intar0
1084             ) {
1085             piv_data->dx[index_y][index_x] = 0.0;
1086             piv_data->dy[index_y][index_x] = 0.0;
1087             piv_data->snr[index_y][index_x] = GPIV_SNR_NAN;
1088             piv_data->peak_no[index_y][index_x] = -1;
1089         }
1090
1091 /*
1092  * Uses old piv data and sets flag peak_no to -1 if:
1093  * for zero offsetting: 2nd Interrogation Area is out of image boundaries
1094  * for zero offsetting with central diff: one of the Interrogation Area's
1095  * is out of image  boundaries
1096  */
1097     } else {
1098         piv_data->dx[index_y][index_x] = piv_data->dx[index_y][index_x];
1099         piv_data->dy[index_y][index_x] = piv_data->dy[index_y][index_x];
1100         piv_data->snr[index_y][index_x] = piv_data->snr[index_y][index_x];
1101         piv_data->peak_no[index_y][index_x] = -1;
1102
1103     }
1104
1105
1106    return err_msg;
1107 }
1108
1109
1110
1111 void
1112 gpiv_piv_isizadapt (const GpivPivPar *piv_par_src,
1113                     GpivPivPar *piv_par_dest,
1114                     gboolean *isiz_last
1115                     )
1116 /*---------------------------------------------------------------------------*/
1117 /**
1118  *     Adjusts interrogation area sizes. For each interrogation sweep,
1119  *     (dest) int_size_i is halved, until it reaches (src)
1120  *     int_size_f. Then, isiz_last is set TRUE, which will avoid
1121  *     changing the interrogation sizes in next calls.
1122  *
1123  *     @param[in] piv_par_src        original parameters
1124  *     @param[out] piv_par_dest       actual parameters, to be modified during sweeps
1125  *     @param[out] isiz_last               flag for last interrogation sweep
1126  *     @return void
1127  */
1128 /*---------------------------------------------------------------------------*/
1129
1130 {
1131
1132 /*     if (piv_par_dest->int_size_i == piv_par_src->int_size_i) */
1133 /*         piv_par_dest->ad_int = 0; */
1134
1135 /*     if (piv_par_dest->ad_int == 1) { */
1136     if ((piv_par_dest->int_size_i) / 2 <= piv_par_src->int_size_f) {
1137         *isiz_last = TRUE;
1138         piv_par_dest->int_size_f =
1139             (piv_par_src->int_size_f - GPIV_DIFF_ISI);
1140         piv_par_dest->int_size_i =
1141             piv_par_src->int_size_f;
1142     } else {
1143         piv_par_dest->int_size_f = piv_par_dest->int_size_i / 2;
1144         piv_par_dest->int_size_i = piv_par_dest->int_size_i / 2;
1145     }
1146
1147 /*     } else if (piv_par_src->int_scheme == GPIV_ZERO_OFF_FORWARD */
1148 /*                || piv_par_src->int_scheme == GPIV_ZERO_OFF_CENTRAL */
1149 /*                || piv_par_src->int_scheme == GPIV_IMG_DEFORM */
1150 /*                ) { */
1151 /*         *isiz_last = TRUE; */
1152 /*         piv_par_dest->ifit = piv_par_src->ifit; */
1153 /*         piv_par_dest->int_size_f = (piv_par_src->int_size_f -  */
1154 /*                                          GPIV_DIFF_ISI); */
1155 /*         piv_par_dest->int_size_i = piv_par_src->int_size_f; */
1156 /*     } */
1157
1158 }
1159
1160
1161
1162 /* #define SAVE_TMP */
1163
1164 gchar *
1165 gpiv_piv_write_deformed_image (GpivImage *image
1166                                )
1167 /*-----------------------------------------------------------------------------
1168  * DESCRIPTION:
1169  *     Stores deformed image to file system with pre defined name to TMPDIR
1170  *     and prints message to stdout
1171  *
1172  * INPUTS:
1173  *     img1:                   first image of PIV image pair
1174  *     img2:                   second image of PIV image pair
1175  *     image_par:              image parameters to be stored in header
1176  *     verbose:                prints the storing to stdout
1177  *
1178  * OUTPUTS:
1179  *
1180  * RETURNS:
1181  *     char * to NULL on success or *err_msg on failure
1182  *---------------------------------------------------------------------------*/
1183 {
1184     gchar *err_msg = NULL;
1185     gchar *def_img;
1186     FILE *fp;
1187
1188
1189     def_img =  g_strdup_printf ("%s%s",  GPIV_DEFORMED_IMG_NAME,
1190                                 GPIV_EXT_PNG_IMAGE);
1191
1192 #ifdef SAVE_TMP
1193     if ((fp = my_utils_fopen_tmp (def_img, "wb")) == NULL) {
1194         err_msg = "gpiv_piv_write_deformed_image: Failure opening for output";
1195         return err_msg;
1196     }
1197
1198     g_message ("gpiv_piv_write_deformed_image: Writing deformed image to: %s",
1199                g_strdup_printf ("%s/%s/%s", tmp_dir, user_name, def_img));
1200 #else
1201     fp = fopen (def_img, "wb");
1202     g_message ("gpiv_piv_write_deformed_image: Writing deformed image to: %s",
1203                def_img);
1204 #endif
1205     if ((err_msg =
1206          gpiv_write_png_image (fp, image, FALSE)) != NULL) {
1207         fclose (fp);
1208         return err_msg;
1209     }
1210
1211     fclose(fp);
1212     g_free (def_img);
1213     return err_msg;
1214 }
1215
1216 #ifdef SAVE_TMP
1217 #undef SAVE_TMP
1218 #endif
1219
1220
1221
1222 static LoVarII *
1223 piv_interrogate_img__init (const GpivImage *image,
1224                            const GpivPivPar *piv_par,
1225                            const GpivValidPar *valid_par,
1226                            Ompp *ompp,
1227                            Mpi *mpi,
1228                            gboolean *sweep_last,
1229                            const gboolean verbose
1230                            )
1231 /* ----------------------------------------------------------------------------
1232  * Initializes variables for gpiv_interrogate_img
1233  */
1234 {
1235     gchar *err_msg = NULL;
1236     int i = 0;
1237
1238     LoVarII *lo = g_new0(LoVarII, 1);         /* struct containing local variables */
1239
1240
1241
1242 #ifdef ENABLE_OMP   /* allocate "thread array" for later use as xy[thr_id] */
1243     ompp->max_nr_thr = omp_get_max_threads();
1244 #else
1245     ompp->max_nr_thr = 1;
1246 #endif /* ENABLE_OMP */
1247
1248     /*
1249      * allocate err_msg_ar[] for each thread
1250      */
1251     lo->err_msg_ar = g_new0(gchar, ompp->max_nr_thr);
1252
1253     /*
1254      * Testing parameters on consistency and initializing derived
1255      * parameters/variables
1256      */
1257     if ((err_msg =
1258          gpiv_piv_testonly_parameters (image->header, piv_par)) != NULL) {
1259         gpiv_warning ("gpiv_piv_interrogate_img: %s", err_msg);
1260         return NULL;
1261     }
1262
1263     if ((err_msg =
1264          gpiv_valid_testonly_parameters (valid_par)) != NULL) {
1265         gpiv_warning ("gpiv_piv_interrogate_img: %s", err_msg);
1266         return NULL;
1267     }
1268
1269     /*
1270      * Local (actualized) parameters
1271      * Setting initial parameters and variables for adaptive grid and
1272      * Interrogation Area dimensions
1273      */
1274     lo->piv_par = gpiv_piv_cp_parameters (piv_par);
1275
1276     if (lo->piv_par->int_scheme == GPIV_ZERO_OFF_FORWARD
1277         || lo->piv_par->int_scheme == GPIV_ZERO_OFF_CENTRAL
1278         || lo->piv_par->int_scheme == GPIV_IMG_DEFORM
1279         || lo->piv_par->int_size_i > lo->piv_par->int_size_f) {
1280         lo->piv_par->int_size_f = lo->piv_par->int_size_i;
1281         *sweep_last = FALSE;
1282     } else {
1283         *sweep_last = TRUE;
1284     }
1285
1286     if (lo->piv_par->int_shift < lo->piv_par->int_size_i / GPIV_SHIFT_FACTOR) {
1287         lo->piv_par->int_shift = lo->piv_par->int_size_i / GPIV_SHIFT_FACTOR;
1288     }
1289
1290     /*
1291      * A copy of the image and PIV data are needed when image deformation is used.
1292      * To keep the algorithm simple (well, what's in the name :), copies are made
1293      * unconditionally.
1294      */
1295     lo->image = gpiv_cp_img(image);
1296     lo->out_piv_data = alloc_pivdata_gridgen(image->header, lo->piv_par);
1297
1298 #ifdef ENABLE_MPI
1299     MPI_Comm_size(MPI_COMM_WORLD, mpi->nprocs);
1300     MPI_Comm_rank(MPI_COMM_WORLD, mpi->rank);
1301 #endif /* ENABLE_MPI */
1302
1303     lo->piv_data = gpiv_cp_pivdata(lo->out_piv_data);
1304
1305 #ifdef DEBUG
1306     gpiv_write_pivdata(NULL, lo->piv_data, FALSE);
1307     fflush(stdout);
1308 #endif /* DEBUG */
1309
1310     gpiv_0_pivdata(lo->piv_data);
1311
1312     /*
1313      * Reads eventually existing fftw wisdom
1314      */
1315     gpiv_fread_fftw_wisdom(1);
1316     gpiv_fread_fftw_wisdom(-1);
1317
1318
1319     return lo;
1320 }
1321
1322
1323
1324 static gchar *
1325 piv_interrogate_img__finish (LoVarII *lo
1326                              )
1327 /* ----------------------------------------------------------------------------
1328  * Saves existing fftw wisdom
1329  * free allocated memory by gpiv_interrogate_img
1330  */
1331 {
1332     gchar *err_msg = NULL;
1333
1334
1335     fftw_cleanup();
1336     gpiv_fwrite_fftw_wisdom(1);
1337     gpiv_fwrite_fftw_wisdom(-1);
1338     gpiv_free_img(lo->image);
1339     gpiv_free_pivdata(lo->piv_data);
1340     gpiv_free_pivdata(lo->out_piv_data);
1341     g_free(*lo->err_msg_ar);
1342
1343     return err_msg;
1344 }
1345
1346
1347
1348 static void
1349 piv_interrogate_img__err (LoVarII *lo,
1350                           GpivFt **ft,
1351                           const int nr_thr)
1352 /**----------------------------------------------------------------------------
1353  * Prints error message (of specific thread in case of OpenMP)
1354  * Free up memory
1355  */
1356 {
1357     int i = 0;
1358
1359
1360     for (i = 0; i < nr_thr; i++) {
1361         gpiv_free_ft(ft[i]);
1362         if (lo->err_msg_ar[i] != NULL) {
1363             gpiv_warning("interrogate_img: thr_id = %d: %s",
1364                           i, lo->err_msg_ar[i]);
1365         }
1366     }
1367     g_free(ft);
1368     gpiv_free_pivdata(lo->out_piv_data);
1369     piv_interrogate_img__finish(lo);
1370
1371
1372     return;
1373 }
1374
1375
1376
1377 static void
1378 piv_weight_kernel_lin (const guint int_size_0,
1379                        GpivCov *w_k
1380                        )
1381 /*-----------------------------------------------------------------------------
1382  * DESCRIPTION:
1383  *      Calculate linear weight kernel values for covariance data
1384  *
1385  * INPUTS:
1386  *
1387  * OUTPUTS:
1388  *      w_k:            Structure containing weighting data
1389  *      int_size_0:     zero-padded interrogation size
1390  *
1391  * RETURNS:
1392  *
1393  *---------------------------------------------------------------------------*/
1394 {
1395     int i, j;
1396     int z_rnl = w_k->z_rnl, z_rnh = w_k->z_rnh, z_rpl = w_k->z_rpl,
1397         z_rph = w_k->z_rph;
1398     int z_cnl = w_k->z_cnl, z_cnh = w_k->z_cnh, z_cpl = w_k->z_cpl,
1399         z_cph = w_k->z_rph;
1400
1401
1402     g_return_if_fail (w_k != NULL);
1403
1404     for (i = z_rnl; i < z_rnh; i++) {
1405         for (j = z_cnl; j < z_cnh; j++) {
1406             w_k->z[i - int_size_0][j - int_size_0] =
1407                 (1 - abs(z_rnh - i) / (int_size_0 / 2.0))
1408                 * (1 - abs(z_cnh - j) / (int_size_0 / 2.0));
1409         }
1410
1411         for (j = z_cpl; j < z_cph; j++) {
1412             w_k->z[i - int_size_0][j] =
1413                 (1 - abs(z_rnh - i) / (int_size_0 / 2.0))
1414                 * (1 - abs(z_cpl - j) / (int_size_0 / 2.0));
1415         }
1416     }
1417
1418
1419     for (i = z_rpl; i < z_rph; i++) {
1420         for (j = z_cnl; j < z_cnh; j++) {
1421             w_k->z[i][j - int_size_0] =
1422                 (1 - abs(z_rpl - i) / (int_size_0 / 2.0))
1423                 * (1 - abs(z_cnh - j) / (int_size_0 / 2.0));
1424         }
1425
1426         for (j = z_cpl; j < z_cph; j++) {
1427             w_k->z[i][j] = (1 - abs(z_rpl - i) / (int_size_0 / 2.0))
1428                 * (1 - abs(z_cpl - j) / (int_size_0 / 2.0));
1429         }
1430     }
1431 }
1432
1433
1434
1435 void
1436 write_cov (int int_x,
1437            int int_y,
1438            int int_size,
1439            float **cov_area,
1440            int weight
1441            )
1442 /*-----------------------------------------------------------------------------
1443  * DESCRIPTION:
1444  *   Prints covariance data
1445  *
1446  * INPUTS:
1447  *   int_x:
1448  *   int_y
1449  *   cov_area:
1450  *   int_size,
1451  *
1452  * SOME MNENOSYNTACTICS OF LOCAL VARIABLES:
1453  *   cov_area:         name of array
1454  *   r:                row
1455  *   c:                column
1456  *   n:                negative displacements ; from 3/4 to 1 of int_size_0
1457  *   p:                positive displacements; from 0 to 1/4 of int_size_0
1458  *   l:                lowest index
1459  *   h:                highest index
1460  *---------------------------------------------------------------------------*/
1461 {
1462     int i, j;
1463     int cov_area_rnl, cov_area_rnh, cov_area_rpl, cov_area_rph,
1464         cov_area_cnl, cov_area_cnh, cov_area_cpl, cov_area_cph;
1465     float weight_kernel;
1466     int int_size_0 = GPIV_ZEROPAD_FACT * int_size;
1467
1468
1469     cov_area_rnl = 3.0 * (int_size_0) / 4 + 1;
1470     cov_area_rnh = int_size_0;
1471     cov_area_rpl = 0;
1472     cov_area_rph = int_size_0 / 4;
1473
1474     cov_area_cnl = 3.0 * (int_size_0) / 4 + 1;
1475     cov_area_cnh = int_size_0;
1476     cov_area_cpl = 0;
1477     cov_area_cph = int_size_0 / 4;
1478
1479
1480     for (i = cov_area_rnl; i < cov_area_rnh; i++) {
1481         for (j = cov_area_cnl; j < cov_area_cnh; j++) {
1482             if (weight == 1) {
1483                 weight_kernel =
1484                     (1 - abs(cov_area_rnh - i) / (int_size_0 / 2.0)) *
1485                   (1 - abs(cov_area_cnh - j) / (int_size_0 / 2.0));
1486             }
1487         }
1488         for (j = cov_area_cpl; j < cov_area_cph; j++) {
1489             if (weight == 1) {
1490                 weight_kernel =
1491                     (1 - abs(cov_area_rnh - i) / (int_size_0 / 2.0)) *
1492                   (1 - abs(cov_area_cpl - j) / (int_size_0 / 2.0));
1493             }
1494         }
1495     }
1496
1497
1498     for (i = cov_area_rpl; i < cov_area_rph; i++) {
1499         for (j = cov_area_cnl; j < cov_area_cnh; j++) {
1500             if (weight == 1) {
1501                 weight_kernel =
1502                     (1 - abs(cov_area_rpl - i) / (int_size_0 / 2.0)) *
1503                   (1 - abs(cov_area_cnh - j) / (int_size_0 / 2.0));
1504             }
1505         }
1506         for (j = cov_area_cpl; j < cov_area_cph; j++) {
1507             if (weight == 1) {
1508                 weight_kernel =
1509                     (1 - abs(cov_area_rpl - i) / (int_size_0 / 2.0)) *
1510                   (1 - abs(cov_area_cpl - j) / (int_size_0 / 2.0));
1511             }
1512         }
1513     }
1514 }
1515
1516
1517
1518 void
1519 gpiv_fread_fftw_wisdom (const gint dir
1520                         )
1521 /*-----------------------------------------------------------------------------
1522  * DESCRIPTION:
1523  *     Reads fftw wisdoms from file and stores into a string
1524  *
1525  * INPUTS:
1526  *     dir:            direction of fft; forward (+1) or inverse (-1)
1527  *
1528  * OUTPUTS:
1529  *
1530  * RETURNS:
1531  *     fftw_wisdom
1532  *---------------------------------------------------------------------------*/
1533 {
1534     gchar *fftw_filename;
1535     FILE *fp;
1536
1537
1538     g_return_if_fail (dir == 1 || dir == -1);
1539
1540 /*
1541  * Forward FFT or inverse FFT
1542  */
1543     if (dir == 1) {
1544         fftw_filename = g_strdup_printf ("%s", FFTWISDOM);
1545     } else {
1546         fftw_filename = g_strdup_printf ("%s", FFTWISDOM_INV);
1547     }
1548
1549     if ((fp = my_utils_fopen_tmp (fftw_filename, "r")) != NULL) {
1550         fftw_import_wisdom_from_file(fp);
1551         fclose(fp);
1552     }
1553
1554     g_free (fftw_filename);
1555 }
1556
1557
1558
1559 void
1560 gpiv_fwrite_fftw_wisdom (const gint dir
1561                          )
1562 /*-----------------------------------------------------------------------------
1563  * DESCRIPTION:
1564  *      Writes fftw wisdoms to a file
1565  *
1566  * INPUTS:
1567  *      dir:            direction of fft; forward (+1) or inverse (-1)
1568  *
1569  * OUTPUTS:
1570  *
1571  * RETURNS:
1572  *
1573  *---------------------------------------------------------------------------*/
1574 {
1575     gchar *fftw_filename;
1576     FILE *fp;
1577
1578     g_return_if_fail (dir == 1 || dir == -1);
1579
1580 /*
1581  * Forward FFT or inverse FFT
1582  */
1583     if (dir == 1) {
1584         fftw_filename = g_strdup_printf ("%s", FFTWISDOM);
1585     } else {
1586         fftw_filename = g_strdup_printf ("%s", FFTWISDOM_INV);
1587     }
1588
1589     if ((fp = my_utils_fopen_tmp (fftw_filename, "w")) != NULL) {
1590         fftw_export_wisdom_to_file(fp);
1591         fclose(fp);
1592     }
1593
1594     fftw_forget_wisdom();
1595     g_free (fftw_filename);
1596 }
1597
1598
1599 /*
1600  * Public functions, original from piv_speed
1601  */
1602
1603 gchar *
1604 gpiv_piv_dxdy_at_new_grid (const GpivPivData *piv_data_src,
1605                            GpivPivData *piv_data_dest
1606                            )
1607 /*---------------------------------------------------------------------------*/
1608 /**
1609  * calculates dx, dy of piv_data_dest from piv_data_src
1610  * by bi-linear interpolation of inner points with shifted knots
1611  * or extrapolation of outer lying points
1612  *
1613  *   dist_:      distance
1614  *   _n:         NORTH
1615  *   _s:         SOUTH
1616  *   _e:         EAST
1617  *   _w:         WEST
1618  *   _nn:        at NORTH of NORTH, etc.
1619  *
1620  *     @param[in] piv_data_src       input piv data
1621  *     @param[out] piv_data_dest      output piv data
1622  *     @return NULL on success or *err_msg on failure
1623  */
1624 /*---------------------------------------------------------------------------*/
1625 {
1626     char c_line[GPIV_MAX_LINES_C][GPIV_MAX_CHARS];
1627     char *err_msg = NULL;
1628
1629     gint *index_n, *index_s, *index_e, *index_w;
1630     gint *index_nn, *index_ss, *index_ee, *index_ww;
1631
1632     float *src_point_x = NULL, *dest_point_x = NULL;
1633     float *src_point_y = NULL, *dest_point_y = NULL;
1634     double *alpha_hor, *alpha_vert;
1635
1636     double epsi = 0.01;
1637     enum VerticalPosition vert_pos;
1638     enum HorizontalPosition hor_pos;
1639     gint i = 0, j = 0;
1640
1641     GpivPivData *gpd = NULL;
1642
1643
1644 /*
1645  * shift the knots of the grid for higher accuracies.
1646  * in order not to affect piv_data_src, a new PIV dataset will be copied
1647  */
1648 #ifdef DEBUG
1649     g_message ("gpiv_piv_dxdy_at_new_grid:: 0, src_nx = %d src_ny = %d dest_nx = %d dest_ny = %d",
1650                piv_data_src->nx, piv_data_src->ny,
1651                piv_data_dest->nx, piv_data_dest->ny);
1652 #endif
1653     gpd = gpiv_cp_pivdata (piv_data_src);
1654
1655
1656     if ((err_msg = gpiv_piv_shift_grid (gpd)) != NULL) {
1657         err_msg = "gpiv_piv_dxdy_at_new_grid: failing  gpiv_piv_shift_grid";
1658         g_warning ("%s", err_msg);
1659         return err_msg;
1660     }
1661
1662
1663     index_n = gpiv_ivector (piv_data_dest->ny);
1664     index_s = gpiv_ivector (piv_data_dest->ny);
1665     index_e = gpiv_ivector (piv_data_dest->nx);
1666     index_w = gpiv_ivector (piv_data_dest->nx);
1667     index_nn = gpiv_ivector (piv_data_dest->ny);
1668     index_ss = gpiv_ivector (piv_data_dest->ny);
1669     index_ee = gpiv_ivector (piv_data_dest->nx);
1670     index_ww = gpiv_ivector (piv_data_dest->nx);
1671
1672     alpha_vert = gpiv_dvector (piv_data_dest->ny);
1673     alpha_hor = gpiv_dvector (piv_data_dest->nx);
1674
1675     src_point_x = gpiv_vector (gpd->nx);
1676     src_point_y = gpiv_vector (gpd->ny);
1677     dest_point_x = gpiv_vector (piv_data_dest->nx);
1678     dest_point_y = gpiv_vector (piv_data_dest->ny);
1679
1680 /*
1681  * Calculate interpolation factors
1682  * in Horizontal direction
1683  */
1684 #ifdef DEBUG
1685     g_message ("gpiv_piv_dxdy_at_new_grid:: 1a, gpd_nx = %d gpd_ny = %d _ny",
1686                gpd->nx, gpd->ny);
1687 #endif
1688     if (gpd->nx >= NMIN_TO_INTERPOLATE) {
1689         for (i = 0, j = 0; j < gpd->nx; j++) {
1690             src_point_x[j] = gpd->point_x[i][j];
1691        }
1692
1693         for (i = 0, j = 0; j < piv_data_dest->nx; j++) {
1694             dest_point_x[j] = piv_data_dest->point_x[i][j];
1695         }
1696
1697         intpol_facts (src_point_x,
1698                       dest_point_x,
1699                       gpd->nx,
1700                       piv_data_dest->nx,
1701                       index_w,
1702                       index_e,
1703                       index_ww,
1704                       index_ee,
1705                       alpha_hor);
1706     } else {
1707         err_msg = "gpiv_piv_dxdy_at_new_grid: Not enough points in horizontal direction";
1708         return err_msg;
1709     }
1710
1711 /*
1712  * Calculate interpolation factors
1713  * in Vertical direction
1714  */
1715     if (gpd->ny >= NMIN_TO_INTERPOLATE) {
1716         for (i = 0, j = 0; i < gpd->ny; i++) {
1717             src_point_y[i] = gpd->point_y[i][j];
1718         }
1719
1720         for (i = 0, j = 0; i < piv_data_dest->ny; i++) {
1721             dest_point_y[i] = piv_data_dest->point_y[i][j];
1722         }
1723
1724         intpol_facts (src_point_y,
1725                       dest_point_y,
1726                       gpd->ny,
1727                       piv_data_dest->ny,
1728                       index_n,
1729                       index_s,
1730                       index_nn,
1731                       index_ss,
1732                       alpha_vert);
1733     } else {
1734         err_msg = "gpiv_piv_dxdy_at_new_grid: Not enough points in horizontal direction";
1735         return err_msg;
1736     }
1737
1738 /*
1739  * Calculate new displacements by bi-lineair interpolation
1740  */
1741     for (i = 0; i < piv_data_dest->ny; i++) {
1742         for (j = 0; j < piv_data_dest->nx; j++) {
1743             piv_data_dest->dx[i][j] = bilinear_interpol
1744                 (alpha_hor[j],
1745                  alpha_vert[i],
1746                  gpd->dx[index_n[i]][index_w[j]],
1747                  gpd->dx[index_n[i]][index_e[j]],
1748                  gpd->dx[index_s[i]][index_w[j]],
1749                  gpd->dx[index_s[i]][index_e[j]]);
1750
1751 #ifdef DEBUG2
1752 g_message ("piv_dxdy_at_new_grid:: alpha_hor[%d]=%f alpha_vert[%d]=%f dx_nw=%f dx_ne=%f dx_sw=%f dx_se=%f => dx=%f",
1753               j, alpha_hor[j], i, alpha_vert[i],
1754               gpd->dx[index_n[i]][index_w[j]],
1755               gpd->dx[index_n[i]][index_e[j]],
1756               gpd->dx[index_s[i]][index_w[j]],
1757               gpd->dx[index_s[i]][index_e[j]],
1758               piv_data_dest->dx[i][j]
1759               );
1760 #endif
1761
1762             piv_data_dest->dy[i][j] = bilinear_interpol
1763                 (alpha_hor[j],
1764                  alpha_vert[i],
1765                  gpd->dy[index_n[i]][index_w[j]],
1766                  gpd->dy[index_n[i]][index_e[j]],
1767                  gpd->dy[index_s[i]][index_w[j]],
1768                  gpd->dy[index_s[i]][index_e[j]]);
1769
1770 #ifdef DEBUG2
1771             g_message ("piv_dxdy_at_new_grid:: alpha_hor[%d]=%f alpha_vert[%d]=%f dy_nw=%f dy_ne=%f dy_sw=%f dy_se=%f => dy=%f",
1772               j, alpha_hor[j], i, alpha_vert[i],
1773               gpd->dy[index_n[i]][index_w[j]],
1774               gpd->dy[index_n[i]][index_e[j]],
1775               gpd->dy[index_s[i]][index_w[j]],
1776               gpd->dy[index_s[i]][index_e[j]],
1777               piv_data_dest->dy[i][j]
1778               );
1779 #endif
1780
1781         }
1782     }
1783
1784
1785     gpiv_free_ivector (index_n);
1786     gpiv_free_ivector (index_s);
1787     gpiv_free_ivector (index_e);
1788     gpiv_free_ivector (index_w);
1789     gpiv_free_ivector (index_nn);
1790     gpiv_free_ivector (index_ss);
1791     gpiv_free_ivector (index_ee);
1792     gpiv_free_ivector (index_ww);
1793
1794     gpiv_free_dvector (alpha_vert);
1795     gpiv_free_dvector (alpha_hor);
1796
1797     gpiv_free_vector (src_point_x);
1798     gpiv_free_vector (src_point_y);
1799     gpiv_free_vector (dest_point_x);
1800     gpiv_free_vector (dest_point_y);
1801
1802     gpiv_free_pivdata (gpd);
1803
1804     return err_msg;
1805 }
1806
1807
1808 gchar *
1809 gpiv_piv_shift_grid (GpivPivData *gpd_src
1810                      )
1811 /*---------------------------------------------------------------------------*/
1812 /**
1813  * shifts the knots of a 2-dimensional grid containing PIV data for improved
1814  * (bi-linear) interpolation
1815  *
1816  * See: T. Blu, P. Thevenaz, "Linear Interpolation Revitalized",
1817  * IEEE Trans. in Image Processing, vol13, no 5, May 2004
1818  *
1819  *     @param[in] piv_data_src   input piv data
1820  *     @return NULL on success or *err_msg on failure
1821  */
1822 /*---------------------------------------------------------------------------*/
1823 {
1824 #define SHIFT 0.2
1825
1826     char *err_msg = NULL;
1827     GpivPivData *h_gpd = NULL;
1828     GpivPivData *v_gpd = NULL;
1829     gfloat fact1 = -SHIFT / (1.0 - SHIFT);
1830     gfloat fact2 = 1.0 / (1 - SHIFT);
1831     gfloat **cx, **cy;
1832     gfloat delta = 0.0;
1833     gint i, j;
1834
1835
1836     delta = gpd_src->point_x[0][1] - gpd_src->point_x[0][0];
1837
1838 /*
1839  * Shift in horizontal (column-wise) direction
1840  */
1841     h_gpd = gpiv_alloc_pivdata (gpd_src->nx, gpd_src->ny);
1842
1843
1844     for (i = 0; i < gpd_src->ny; i++) {
1845         for (j = 0; j < gpd_src->nx; j++) {
1846 /*
1847  * Shift the knot (sample points)
1848  */
1849             h_gpd->point_x[i][j] = gpd_src->point_x[i][j] + SHIFT * delta;
1850             h_gpd->point_y[i][j] = gpd_src->point_y[i][j];
1851             if (j == 0) {
1852                 h_gpd->dx[i][j] = gpd_src->dx[i][j];
1853                 h_gpd->dy[i][j] = gpd_src->dy[i][j];
1854             } else {
1855 /*
1856  * Calculate value at shifted knot
1857  */
1858                 h_gpd->dx[i][j] = fact1 * h_gpd->dx[i][j-1] + fact2 *
1859                     gpd_src->dx[i][j];
1860                 h_gpd->dy[i][j] = fact1 * h_gpd->dy[i][j-1] + fact2 *
1861                     gpd_src->dy[i][j];
1862             }
1863         }
1864     }
1865
1866
1867 /*
1868  * Shift in vertical (row-wise) direction by using the horizontal shifted nodes
1869  */
1870     v_gpd = gpiv_alloc_pivdata (gpd_src->nx, gpd_src->ny);
1871
1872
1873     for (i = 0; i < gpd_src->ny; i++) {
1874         for (j = 0; j < gpd_src->nx; j++) {
1875             v_gpd->point_x[i][j] = h_gpd->point_x[i][j];
1876             v_gpd->point_y[i][j] = h_gpd->point_y[i][j] + SHIFT * delta;
1877             if (i == 0) {
1878                 v_gpd->dx[i][j] = h_gpd->dx[i][j];
1879                 v_gpd->dy[i][j] = h_gpd->dy[i][j];
1880             } else {
1881                 v_gpd->dx[i][j] = fact1 * v_gpd->dx[i-1][j] + fact2 *
1882                     h_gpd->dx[i][j];
1883                 v_gpd->dy[i][j] = fact1 * v_gpd->dy[i-1][j] + fact2 *
1884                     h_gpd->dy[i][j];
1885             }
1886         }
1887     }
1888
1889
1890 /*     gpiv_free_pivdata (gpd_src); */
1891 /*     gpd_src = gpiv_cp_pivdata (v_gpd); */
1892
1893     for (i = 0; i < gpd_src->ny; i++) {
1894         for (j = 0; j < gpd_src->nx; j++) {
1895             gpd_src->point_x[i][j] = v_gpd->point_x[i][j];
1896             gpd_src->point_y[i][j] = v_gpd->point_y[i][j];
1897             gpd_src->dx[i][j] = v_gpd->dx[i][j];
1898             gpd_src->dy[i][j] = v_gpd->dy[i][j];
1899             gpd_src->snr[i][j] = v_gpd->snr[i][j];
1900             gpd_src->peak_no[i][j] = v_gpd->peak_no[i][j];
1901         }
1902     }
1903
1904 /*     gpiv_write_pivdata (NULL, gpd_src, FALSE); */
1905
1906
1907     gpiv_free_pivdata (h_gpd);
1908     gpiv_free_pivdata (v_gpd);
1909     return err_msg;
1910 #undef SHIFT
1911 }
1912
1913
1914
1915 GpivPivData *
1916 gpiv_piv_gridgen (const guint nx,
1917                   const guint ny,
1918                   const GpivImagePar *image_par,
1919                   const GpivPivPar *piv_par
1920                   )
1921 /*-----------------------------------------------------------------------------
1922  * DESCRIPTION:
1923  *     Generates grid by Calculating the positions of interrogation areas
1924  *     Substitutes gpiv_piv_select_int_point
1925  *
1926  * INPUTS:
1927  *     nx               number of columns
1928  *     ny               number of rows
1929  *     @image_par:      structure of image parameters
1930  *     @piv_par:   structure of piv pivuation parameters
1931  *
1932  * OUTPUTS:
1933  *     @out_data:       output piv data from image analysis
1934  *
1935  * RETURNS:
1936  *     %char * NULL on success or *err_msg on failure
1937  *---------------------------------------------------------------------------*/
1938 {
1939     GpivPivData *piv_data = NULL;
1940     gchar *err_msg = NULL;
1941     int row, column, row_1, column_1, i, j;
1942     int row_max, row_min, column_max, column_min;
1943
1944     int ncolumns = image_par->ncolumns;
1945     int nrows = image_par->nrows;
1946
1947     int int_geo = piv_par->int_geo;
1948     int row_start = piv_par->row_start;
1949     int row_end = piv_par->row_end;
1950     int col_start = piv_par->col_start;
1951     int col_end= piv_par->col_end;
1952     int int_line_col = piv_par->int_line_col;
1953     int int_line_col_start = piv_par->int_line_col_start;
1954     int int_line_col_end = piv_par->int_line_col_end;
1955     int int_line_row = piv_par->int_line_row;
1956     int int_line_row_start = piv_par->int_line_row_start;
1957     int int_line_row_end = piv_par->int_line_row_end;
1958     int int_point_col = piv_par->int_point_col;
1959     int int_point_row = piv_par->int_point_row;
1960     int int_size_f = piv_par->int_size_f;
1961     int int_size_i = piv_par->int_size_i;
1962     int int_shift = piv_par->int_shift;
1963     int pre_shift_row = piv_par->pre_shift_row;
1964     int pre_shift_col = piv_par->pre_shift_col;
1965
1966
1967 /*     g_return_val_if_fail (piv_data->point_x != NULL, "piv_data->point_x != NULL"); */
1968 /*     g_return_val_if_fail (piv_data->point_y != NULL, "piv_data->point_y != NULL"); */
1969
1970     row_min = gpiv_min (-int_size_f / 2 + 1,
1971                         pre_shift_row - int_size_i / 2 + 1);
1972     column_min = gpiv_max (-int_size_f / 2 + 1,
1973                            pre_shift_col - int_size_i / 2 + 1);
1974     row_max = gpiv_max (int_size_f / 2, pre_shift_row + int_size_i / 2);
1975     column_max = gpiv_max (int_size_f / 2, pre_shift_col + int_size_i / 2);
1976
1977
1978 /*
1979  * Creates piv_data and centre points for one single interrogation area
1980  */
1981     piv_data = gpiv_alloc_pivdata (nx, ny);
1982
1983     if (int_geo == GPIV_POINT) {
1984         piv_data->point_y[0][0] = int_point_row;
1985         piv_data->point_x[0][0] = int_point_col;
1986
1987
1988 /*
1989  * Creates centre points for one single row
1990  */
1991     } else if (int_geo == GPIV_LINE_R) {
1992         if ((int_size_f - int_size_i) / 2 + pre_shift_col < 0) {
1993             column_1 = int_line_col_start -
1994                 ((int_size_f - int_size_i) / 2 +
1995                  pre_shift_col) + int_size_f / 2 - 1;
1996         } else {
1997             column_1 = int_line_col_start + int_size_f / 2 - 1;
1998         }
1999
2000         for (i = 0, j = 0, row = int_line_row, column = column_1;
2001              j < piv_data->nx; j++, column += int_shift) {
2002             piv_data->point_y[i][j] = row;
2003             piv_data->point_x[i][j] = column;
2004         }
2005
2006
2007 /*
2008  * Creates centre points for one single column
2009  */
2010     } else  if (int_geo == GPIV_LINE_C) {
2011         if ((int_size_f - int_size_i) / 2 + pre_shift_row < 0) {
2012             row_1 = int_line_row_start -
2013                 ((int_size_f - int_size_i) / 2 +
2014                  pre_shift_row) + int_size_f / 2 - 1;
2015         } else {
2016             row_1 = int_line_row_start + int_size_f / 2 - 1;
2017         }
2018
2019         for (i = 0, j = 0, column = int_line_col, row = row_1;
2020              i < piv_data->ny; i++, row += int_shift) {
2021             piv_data->point_y[i][j] = row;
2022             piv_data->point_x[i][j] = column;
2023         }
2024
2025
2026 /*
2027  * Creates an array of  centre points of the Interrrogation Area's:
2028  * int_ar_x and int_ar_y within the defined region of the image
2029  */
2030     } else if (int_geo == GPIV_AOI) {
2031         if ((int_size_f - int_size_i) / 2 + pre_shift_row < 0) {
2032             row_1 =
2033                 row_start - ((int_size_f - int_size_i) / 2 +
2034                              pre_shift_row) + int_size_f / 2 - 1;
2035         } else {
2036             row_1 = row_start + int_size_f / 2 - 1;
2037         }
2038
2039         if ((int_size_f - int_size_i) / 2 + pre_shift_col < 0) {
2040             column_1 =
2041                 col_start - ((int_size_f - int_size_i) / 2 +
2042                              pre_shift_col) + int_size_f / 2 - 1;
2043         } else {
2044             column_1 = col_start + int_size_f / 2 - 1;
2045         }
2046
2047         for (i = 0, row = row_1; i < ny; i++, row += int_shift) {
2048             for (j = 0, column = column_1; j < nx;
2049                  j++, column += int_shift) {
2050                 piv_data->point_y[i][j] = row;
2051                 piv_data->point_x[i][j] = column;
2052             }
2053         }
2054     } else {
2055         err_msg = "gpiv_piv_gridgen: should not arrive here";
2056         gpiv_warning ("%s", err_msg);
2057         return NULL;
2058     }
2059
2060
2061     return piv_data;
2062 }
2063
2064
2065
2066 GpivPivData *
2067 gpiv_piv_gridadapt (const GpivImagePar *image_par,
2068                     const GpivPivPar *piv_par_src,
2069                     GpivPivPar *piv_par_dest,
2070                     const GpivPivData *piv_data,
2071                     const guint sweep,
2072                     gboolean *grid_last
2073                     )
2074 /*-----------------------------------------------------------------------------
2075  * DESCRIPTION:
2076  *     Adjust grid nodes if zero_off or adaptive interrogation
2077  *     area has been used. This is performed by modifying int_shift equal
2078  *     to int_shift / GPIV_SHIFT_FACTOR , until it reaches (src)
2079  *     int_shift. Then, grid_last is set TRUE, which will avoid
2080  *     changing the interrogation shift in next calls and signal the
2081  *     (while loop in) the calling function.
2082  *
2083  * INPUTS:
2084  *     @image_par:                 image parameters
2085  *     @piv_par_src:               piv parameters
2086  *     @piv_data:                  input PIV data
2087  *     @sweep:                     interrogation sweep step
2088  *
2089  * OUTPUTS:
2090  *     @image_par:                image parameters
2091  *     @piv_par_dest:             modified piv parameters
2092  *     @grid_last:                flag if final grid refinement has been
2093  *                                reached
2094  *
2095  * RETURNS:
2096  *     piv_data:                 modified PIV data
2097  *---------------------------------------------------------------------------*/
2098 {
2099     GpivPivData *pd = NULL;
2100     gint local_int_shift, local_int_size_f, local_int_size_i;
2101     gint LOCAL_SHIFT_FACTOR = 2;
2102
2103     guint nx, ny;
2104
2105
2106     local_int_shift = piv_par_dest->int_shift / LOCAL_SHIFT_FACTOR;
2107     if (local_int_shift <= piv_par_src->int_shift) {
2108         *grid_last = TRUE;
2109     }
2110
2111     if (*grid_last == FALSE) {
2112 /*
2113  * - renew grid of PIV dataset
2114  * - calculate displacements at new grid points
2115  */
2116         piv_par_dest->int_shift = piv_par_dest->int_shift /
2117             LOCAL_SHIFT_FACTOR;
2118         gpiv_piv_count_pivdata_fromimage (image_par, piv_par_dest, &nx, &ny);
2119         pd = gpiv_piv_gridgen (nx, ny, image_par, piv_par_dest);
2120         gpiv_piv_dxdy_at_new_grid (piv_data, pd);
2121
2122     } else {
2123 /*
2124  * reset int_shift (and data positions) to the originally defined
2125  * parameter value.
2126  * For the last grid adaption, the number of interrogation area's may
2127  * not have been doubled perse, as int_size may be of
2128  * arbitrary quantity.
2129  */
2130
2131         piv_par_dest->int_shift =  piv_par_src->int_shift;
2132 /*
2133  * Set int_size_f and int_size_i of piv_par_dest temporarly to the
2134  * original settings, so that an identic grid will be constructued as
2135  * during the gpiv_gridgen call.
2136  */
2137         local_int_size_f =  piv_par_dest->int_size_f;
2138         local_int_size_i =  piv_par_dest->int_size_i;
2139         piv_par_dest->int_size_f = piv_par_src->int_size_f;
2140         piv_par_dest->int_size_i = piv_par_src->int_size_i;
2141         gpiv_piv_count_pivdata_fromimage (image_par, piv_par_dest, &nx, &ny);
2142         pd = gpiv_piv_gridgen (nx, ny, image_par, piv_par_dest);
2143         piv_par_dest->int_size_f = local_int_size_f;
2144         piv_par_dest->int_size_i = local_int_size_i;
2145
2146         gpiv_piv_dxdy_at_new_grid (piv_data, pd);
2147     }
2148
2149
2150     return pd;
2151 }
2152
2153
2154 /*
2155  * Local functions
2156  */
2157
2158 static GpivPivData *
2159 alloc_pivdata_gridgen (const GpivImagePar *image_par,
2160                        const GpivPivPar *piv_par
2161                        )
2162 /*-----------------------------------------------------------------------------
2163  * Determines the number of grid points, allocating memory for output
2164  * data and generates the grid
2165  */
2166 {
2167     GpivPivData *piv_data = NULL;
2168     gchar *err_msg = NULL;
2169     GpivPivPar *lo_piv_par = NULL;
2170     guint nx, ny;
2171
2172     if ((lo_piv_par = gpiv_piv_cp_parameters (piv_par)) == NULL) {
2173         gpiv_error ("LIBGPIV internal error: alloc_pivdata_gridgen: failing gpiv_piv_cp_parameters");
2174     }
2175
2176     if (piv_par->int_size_i > piv_par->int_size_f
2177         && piv_par->int_shift < piv_par->int_size_i / GPIV_SHIFT_FACTOR) {
2178         lo_piv_par->int_shift = lo_piv_par->int_size_i / GPIV_SHIFT_FACTOR;
2179     }
2180
2181     if ((err_msg =
2182          gpiv_piv_count_pivdata_fromimage (image_par, lo_piv_par, &nx, &ny))
2183         != NULL) {
2184         g_message ("LIBGPIV internal error: alloc_pivdata_gridgen: %s", err_msg);
2185         return NULL;
2186     }
2187
2188
2189     if ((piv_data =
2190          gpiv_piv_gridgen (nx, ny, image_par, lo_piv_par))
2191         == NULL) {
2192         g_message ("LIBGPIV internal error: alloc_pivdata_gridgen: failing gpiv_piv_gridgen");
2193         return NULL;
2194     }
2195
2196
2197     return piv_data;
2198 }
2199
2200
2201
2202 static void
2203 report_progress (gint *progress_old,
2204                  gint index_y,
2205                  gint index_x,
2206                  GpivPivData *piv_data,
2207                  GpivPivPar *piv_par,
2208                  gint sweep,
2209                  gfloat cum_residu
2210                  )
2211 /*-----------------------------------------------------------------------------
2212  * Printing the progress (between 0 and 100) of piv interrogation to stdout
2213  */
2214 {
2215     gint progress  = 100 * (index_y * piv_data->nx + index_x + 1) /
2216         (piv_data->nx * piv_data->ny);
2217
2218 #ifdef ENABLE_MPI
2219     gint rank, size;
2220 #endif
2221
2222
2223     if (progress != *progress_old) {
2224         *progress_old = progress;
2225
2226         if (index_y > 0 || index_x > 0)
2227             printf ("\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b");
2228
2229         if (piv_par->int_scheme == GPIV_ZERO_OFF_FORWARD
2230             || piv_par->int_scheme == GPIV_ZERO_OFF_CENTRAL
2231             || piv_par->int_scheme == GPIV_ZERO_OFF_CENTRAL
2232             || piv_par->int_scheme == GPIV_IMG_DEFORM
2233             || piv_par->int_size_i > piv_par->int_size_f) {
2234             printf
2235                 ("\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b"
2236                  "\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b"
2237                  "\b\b\b\b\b\b\b\b\b\b\b"
2238                  "\b\b\b\b\b\b\b\b\b\b\b\b\b\b"
2239                  );
2240
2241 #ifdef ENABLE_OMP
2242             printf ("\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b");
2243             printf ("nr_thr = %2d thr_id = %2d ",
2244                     omp_get_num_threads(), omp_get_thread_num());
2245 #endif
2246
2247 #ifdef ENABLE_MPI
2248             MPI_Comm_rank(MPI_COMM_WORLD, &rank);
2249             MPI_Comm_size(MPI_COMM_WORLD, &size);
2250             printf ("\b\b\b\b\b\b\b\b\b\b\b\b");
2251             printf ("rank %2d/%2d, ", rank, size);
2252 #endif
2253             printf ("sweep #%2d, int_size = %d int_shift = %d residu = %.3f: ",
2254                     sweep, piv_par->int_size_f, piv_par->int_shift,
2255                     cum_residu);
2256         }
2257
2258         printf ("%3d %%", progress);
2259         fflush (stdout);
2260     }
2261 }
2262
2263
2264
2265 static gboolean
2266 assign_img2intarr (gint ipoint_x,
2267                    gint ipoint_y,
2268                    guint16 **img_1,
2269                    guint16 **img_2,
2270                    gint int_size_f,
2271                    gint int_size_i,
2272                    gfloat **int_area_1,
2273                    gfloat **int_area_2,
2274                    gint pre_shift_row,
2275                    gint pre_shift_col,
2276                    gint nrows,
2277                    gint ncolumns,
2278                    gint int_size_0
2279                    )
2280 /*-----------------------------------------------------------------------------
2281  * Assigns image data to the interrogation area arrays in a straightforward way
2282  */
2283 {
2284     gint m, n;
2285     gint arg_int1_r = ipoint_y  - int_size_f / 2 + 1;
2286     gint arg_int1_c = ipoint_x  - int_size_f / 2 + 1;
2287     gint arg_int2_r = ipoint_y  - int_size_i / 2 + 1;
2288     gint arg_int2_c = ipoint_x  - int_size_i / 2 + 1;
2289
2290     gboolean flag_valid;
2291
2292
2293     assert (img_1[0] != NULL);
2294     assert (img_2[0] != NULL);
2295     assert (int_area_1[0] != NULL);
2296     assert (int_area_2[0] != NULL);
2297
2298 /*
2299  * Check if Interrogation Areas are within the image boundaries.
2300  * Principally arg_int1_r,c don't have to be tested as
2301  * int_size_i >= int_size_f, but has been kept to maintain generallity with the
2302  * other assign_img2intarr* functions
2303  */
2304     if ((arg_int1_r) >= 0
2305          && (arg_int1_r + int_size_f - 1) < nrows
2306          && (arg_int1_c) >= 0
2307          && (arg_int1_c + int_size_f - 1) < ncolumns
2308
2309          && (arg_int2_r) >= 0
2310          && (arg_int2_r + int_size_i - 1) < nrows
2311          && (arg_int2_c) >= 0
2312          && (arg_int2_c + int_size_i - 1) < ncolumns) {
2313         flag_valid = TRUE;
2314
2315     } else {
2316         flag_valid = FALSE;
2317     }
2318
2319
2320     if (flag_valid == TRUE) {
2321 /*
2322  * reset int_area_1, int_area_2 values
2323  */
2324         memset (int_area_1[0], 0.0, (sizeof(gfloat)) * int_size_0 * int_size_0);
2325         memset (int_area_2[0], 0.0, (sizeof(gfloat)) * int_size_0 * int_size_0);
2326
2327         for (m = 0; m < int_size_f; m++) {
2328             for (n = 0; n < int_size_f; n++) {
2329                 int_area_1[m][n] =
2330                     (float) img_1[m + arg_int1_r][n + arg_int1_c];
2331             }
2332         }
2333
2334         for (m = 0; m < int_size_i; m++) {
2335             for (n = 0; n < int_size_i; n++) {
2336                 int_area_2[m][n] =
2337                     (float) img_2[m + arg_int2_r][n + arg_int2_c];
2338             }
2339         }
2340     }
2341
2342
2343     return flag_valid;
2344 }
2345
2346
2347
2348 static gboolean
2349 assign_img2intarr_central (gint ipoint_x,
2350                            gint ipoint_y,
2351                            guint16 **img_1,
2352                            guint16 **img_2,
2353                            gint int_size_f,
2354                            gint int_size_i,
2355                            gfloat **int_area_1,
2356                            gfloat **int_area_2,
2357                            gint pre_shift_row,
2358                            gint pre_shift_col,
2359                            gint nrows,
2360                            gint ncolumns,
2361                            gint int_size_0
2362                            )
2363 /*-----------------------------------------------------------------------------
2364  * Assigns image data to the interrogation area arrays using the central
2365  * differential scheme
2366  */
2367 {
2368     gint m, n;
2369     gint idum = gpiv_max((int_size_i - int_size_f) / 2, 0);
2370     gint arg_int1_r = ipoint_y - int_size_f / 2 + 1 - pre_shift_row / 2 -
2371         pre_shift_row % 2;
2372     gint arg_int1_c = ipoint_x - int_size_f / 2 + 1 - pre_shift_col / 2 -
2373         pre_shift_col % 2;
2374     gint arg_int2_r = ipoint_y - int_size_i / 2 + 1 + pre_shift_row / 2;
2375     gint arg_int2_c = ipoint_x - int_size_i / 2 + 1 + pre_shift_col / 2;
2376     gboolean flag_valid;
2377
2378
2379     assert (img_1[0] != NULL);
2380     assert (img_2[0] != NULL);
2381     assert (int_area_1[0] != NULL);
2382     assert (int_area_2[0] != NULL);
2383 /*
2384  * Check if Interrogation Areas are within the image boundaries.
2385  */
2386      if ((arg_int1_r) >= 0
2387          && (arg_int1_r + int_size_f - 1) < nrows
2388          && (arg_int1_c) >= 0
2389          && (arg_int1_c + int_size_f - 1) < ncolumns
2390
2391          && (arg_int2_r) >= 0
2392          && (arg_int2_r + int_size_i - 1) < nrows
2393          && (arg_int2_c) >= 0
2394          && (arg_int2_c + int_size_i - 1) < ncolumns) {
2395          flag_valid = TRUE;
2396      } else {
2397          flag_valid = FALSE;
2398      }
2399
2400
2401     if (flag_valid == TRUE) {
2402 /*
2403  * reset int_area_1, int_area_2 values
2404  */
2405         memset(int_area_1[0], 0.0, (sizeof(gfloat)) * int_size_0 * int_size_0);
2406         memset(int_area_2[0], 0.0, (sizeof(gfloat)) * int_size_0 * int_size_0);
2407
2408         for (m = 0; m < int_size_f; m++) {
2409             for (n = 0; n < int_size_f; n++) {
2410                 int_area_1[m + idum][n + idum] =
2411                     (float) img_1[m + arg_int1_r][n + arg_int1_c];
2412             }
2413         }
2414
2415
2416         for (m = 0; m < int_size_i; m++) {
2417             for (n = 0; n < int_size_i; n++) {
2418                 int_area_2[m][n] =
2419                     (float) img_2[m + arg_int2_r][n + arg_int2_c];
2420             }
2421         }
2422
2423     }
2424
2425
2426     return flag_valid;
2427 }
2428
2429
2430
2431 static gboolean
2432 assign_img2intarr_forward (gint ipoint_x,
2433                            gint  ipoint_y,
2434                            guint16 **img_1,
2435                            guint16 **img_2,
2436                            gint int_size_f,
2437                            gint int_size_i,
2438                            gfloat **int_area_1,
2439                            gfloat **int_area_2,
2440                            gint pre_shift_row,
2441                            gint pre_shift_col,
2442                            gint nrows,
2443                            gint ncolumns,
2444                            gint int_size_0
2445                            )
2446 /*-----------------------------------------------------------------------------
2447  * Assigns image data to the interrogation area arrays for forward differential
2448  * scheme
2449  */
2450 {
2451     gint m, n;
2452     gint idum = gpiv_max((int_size_i - int_size_f) / 2, 0);
2453     gint arg_int1_r = ipoint_y - int_size_f / 2 + 1 + pre_shift_row + idum;
2454     gint arg_int1_c = ipoint_x - int_size_f / 2 + 1 + pre_shift_col + idum;
2455     gint arg_int2_r = ipoint_y - int_size_i / 2 + 1 + pre_shift_row;
2456     gint arg_int2_c = ipoint_x - int_size_i / 2 + 1 + pre_shift_col;
2457     gboolean flag_valid;
2458
2459
2460     assert (img_1[0] != NULL);
2461     assert (img_2[0] != NULL);
2462     assert (int_area_1[0] != NULL);
2463     assert (int_area_2[0] != NULL);
2464 /*
2465  * Check if Interrogation Areas are within the image boundaries.
2466  */
2467     if ((arg_int1_r) >= 0
2468         && (arg_int1_r + int_size_f - 1) < nrows
2469         && (arg_int1_c) >= 0
2470         && (arg_int1_c  + int_size_f - 1) < ncolumns
2471
2472         && (arg_int2_r) >= 0
2473         && (arg_int2_r + int_size_i - 1) < nrows
2474         && (arg_int2_c) >= 0
2475         && (arg_int2_c  + int_size_i - 1) < ncolumns) {
2476         flag_valid = TRUE;
2477     } else {
2478         flag_valid = FALSE;
2479     }
2480
2481
2482     if (flag_valid == TRUE) {
2483 /*
2484  * reset int_area_1, int_area_2 values
2485  */
2486         memset(int_area_1[0], 0.0,
2487                (sizeof(float)) * int_size_0 * int_size_0);
2488         memset(int_area_2[0], 0.0,
2489                (sizeof(float)) * int_size_0 * int_size_0);
2490
2491         for (m = 0; m < int_size_f; m++) {
2492             for (n = 0; n < int_size_f; n++) {
2493                 int_area_1[m + idum][n + idum] =
2494                     (float) img_1[m + arg_int1_r][n + arg_int1_c];
2495             }
2496         }
2497
2498
2499         for (m = 0; m < int_size_i; m++) {
2500             for (n = 0; n < int_size_i; n++) {
2501                 int_area_2[m][n] =
2502                     (float) img_2[m + arg_int2_r][n + arg_int2_c];
2503             }
2504         }
2505
2506     }
2507
2508
2509     return flag_valid;
2510 }
2511
2512
2513
2514 static float
2515 int_mean_old (guint16 **img,
2516               int int_size,
2517               int int_size_i,
2518               int ipoint_x,
2519               int ipoint_y
2520               )
2521 /* ----------------------------------------------------------------------------
2522  * calculates mean image value from which image data are taken
2523  */
2524 {
2525     int m = 0, n = 0, idum = gpiv_max((int_size_i - int_size) / 2, 0);
2526     int int_area_sum = 0;
2527     float mean;
2528
2529
2530     assert (img[0] != NULL);
2531
2532     for (m = 0; m < int_size; m++) {
2533         for (n = 0; n < int_size; n++) {
2534             int_area_sum +=
2535                 img[m + ipoint_y - int_size_i / 2 + 1 + idum]
2536                 [n + ipoint_x - int_size_i / 2 + 1 + idum];
2537         }
2538     }
2539
2540     mean = int_area_sum / (int_size * int_size);
2541
2542
2543     return mean;
2544 }
2545
2546
2547
2548 static gfloat
2549 int_mean (gfloat **int_area,
2550           gint int_size
2551           )
2552 /* ----------------------------------------------------------------------------
2553  * calculates mean value from interrogation area intensities
2554  */
2555 {
2556     int m = 0, n = 0;
2557     gfloat int_area_sum = 0;
2558     gfloat mean = 0.0;
2559
2560
2561     assert (int_area[0] != NULL);
2562
2563     for (m = 0; m < int_size; m++) {
2564         for (n = 0; n < int_size; n++) {
2565             int_area_sum += int_area[m][n];
2566         }
2567     }
2568
2569     mean = int_area_sum / (int_size * int_size);
2570
2571
2572     return mean;
2573 }
2574
2575
2576
2577 static gfloat
2578 int_range (gfloat **int_area,
2579           gint int_size
2580           )
2581 /* ----------------------------------------------------------------------------
2582  * calculates range of values from interrogation area intensities
2583  */
2584 {
2585     int m = 0, n = 0;
2586     gfloat int_area_sum = 0;
2587     gfloat min = 10.0e9, max = -10.0e9, range = 0.0;
2588
2589
2590     assert (int_area[0] != NULL);
2591
2592     for (m = 0; m < int_size; m++) {
2593         for (n = 0; n < int_size; n++) {
2594             if (int_area[m][n] > max) max = int_area[m][n];
2595             if (int_area[m][n] < min) min = int_area[m][n];
2596         }
2597     }
2598
2599     range = max - min;
2600
2601
2602     return range;
2603 }
2604
2605
2606
2607 static gboolean
2608 int_const (gfloat **int_area,
2609            const guint int_size
2610            )
2611 /* ----------------------------------------------------------------------------
2612  * Tests if all intesities values with an interrogation area are equal
2613  */
2614 {
2615     int m = 0, n = 0;
2616     gboolean flag = TRUE;
2617     gfloat val;
2618
2619
2620     assert (int_area[0] != NULL);
2621
2622     val = int_area[0][0];
2623     for (m = 1; m < int_size; m++) {
2624         for (n = 1; n < int_size; n++) {
2625             if (int_area[m][n] != val) flag = FALSE;
2626         }
2627     }
2628
2629
2630     return flag;
2631 }
2632
2633
2634
2635 static void
2636 cov_min_max (GpivCov *cov
2637              )
2638 /* ----------------------------------------------------------------------------
2639  * calculates minimum and maximum in cov
2640  */
2641 {
2642     gfloat max = -10.0e+9, min = 10.0e+9;
2643     gint z_rl = cov->z_rl, z_rh = cov->z_rh, z_cl = cov->z_cl,
2644         z_ch = cov->z_ch;
2645     gint i, j;
2646
2647
2648     for (i = z_rl + 1; i < z_rh - 1; i++) {
2649         for (j = z_cl + 1; j < z_ch - 1; j++) {
2650             if (cov->z[i][j] > max) max = cov->z[i][j];
2651             if (cov->z[i][j] < min) min = cov->z[i][j];
2652         }
2653     }
2654
2655     cov->min = min;
2656     cov->max = max;
2657 }
2658
2659
2660 static void
2661 search_top (GpivCov *cov,
2662             gint peak_act,
2663             gint x_corr,
2664             gint sweep,
2665             gint i_skip_act,
2666             gint j_skip_act,
2667             float *z_max,
2668             long *i_max,
2669             long *j_max
2670             )
2671 /* ----------------------------------------------------------------------------
2672  * searches top in cov. function
2673  */
2674 {
2675     gint h, i, j;
2676     gint z_rl = cov->z_rl, z_rh = cov->z_rh, z_cl = cov->z_cl,
2677         z_ch = cov->z_ch;
2678
2679
2680     for (h = 1; h <= peak_act + 1; h++) {
2681         z_max[h] = -1.0;
2682         i_max[h] = -2;
2683         j_max[h] = -3;
2684     }
2685
2686     for (h = 1; h <= peak_act; h++) {
2687         for (i = z_rl + 1; i < z_rh - 1; i++) {
2688             for (j = z_cl + 1; j < z_ch - 1; j++) {
2689
2690                 if (x_corr == 1 || (sweep == 0 || (i != i_skip_act ||
2691                                                    j != j_skip_act))) {
2692                     if (h == 1
2693                         || (h == 2
2694                             && (i != i_max[1] || j != j_max[1]))
2695                         || (h == 3
2696                             && (i != i_max[1] || j != j_max[1])
2697                             && (i != i_max[2] || j != j_max[2]))) {
2698                         if (cov->z[i][j] > z_max[h]) {
2699                             if ((cov->z[i][j] >= cov->z[i - 1][j]) &&
2700                                 (cov->z[i][j] >= cov->z[i + 1][j]) &&
2701                                 (cov->z[i][j] >= cov->z[i][j - 1]) &&
2702                                 (cov->z[i][j] >= cov->z[i][j + 1])) {
2703                                 z_max[h] = cov->z[i][j];
2704                                 i_max[h] = i;
2705                                 j_max[h] = j;
2706                             }
2707                         }
2708                     }
2709
2710                 }
2711             }
2712         }
2713     }
2714 }
2715
2716
2717
2718 static char *
2719 cov_subtop (float **z,
2720             long *i_max,
2721             long *j_max,
2722             float *epsi_x,
2723             float *epsi_y,
2724             int ifit,
2725             int peak_act
2726             )
2727 /*-----------------------------------------------------------------------------
2728  * Calculates particle displacements at sub-pixel level
2729  */
2730 {
2731     char *err_msg = NULL;
2732     double A_log, B_log, C_log;
2733     double min = 10e-3;
2734     gboolean flag = TRUE;
2735
2736
2737     if (ifit == GPIV_GAUSS) {
2738 /*
2739  * sub-pixel estimator by gaussian fit
2740  */
2741         if (z[i_max[peak_act]][j_max[peak_act]] > min) {
2742             C_log = log((double) z[i_max[peak_act]][j_max[peak_act]]);
2743         } else {
2744             flag = FALSE;
2745         }
2746
2747         if (flag && z[i_max[peak_act] - 1][j_max[peak_act]] > min) {
2748             A_log = log((double) z[i_max[peak_act] - 1][j_max[peak_act]]);
2749         } else {
2750             flag = FALSE;
2751         }
2752
2753         if (flag && z[i_max[peak_act] + 1][j_max[peak_act]] > min) {
2754             B_log = log((double) z[i_max[peak_act] + 1][j_max[peak_act]]);
2755         } else {
2756             flag = FALSE;
2757         }
2758
2759         if (flag && (2 * (A_log + B_log - 2 * C_log)) != 0.0) {
2760             *epsi_y = (A_log - B_log) / (2 * (A_log + B_log - 2 * C_log));
2761         } else {
2762             *epsi_y = 0.0;
2763             peak_act = -1;
2764             err_msg = "epsi_y = 0.0";
2765             flag = FALSE;
2766         }
2767
2768
2769         if (flag && z[i_max[peak_act]][j_max[peak_act] - 1] != 0.0) {
2770             A_log = log((double) z[i_max[peak_act]][j_max[peak_act] - 1]);
2771         } else {
2772             flag = FALSE;
2773         }
2774
2775         if (flag && z[i_max[peak_act]][j_max[peak_act] + 1] != 0.0) {
2776             B_log = log((double) z[i_max[peak_act]][j_max[peak_act] + 1]);
2777         } else {
2778             flag = FALSE;
2779         }
2780
2781         if (flag && (2 * (A_log + B_log - 2 * C_log)) != 0.0) {
2782             *epsi_x = (A_log - B_log) / (2 * (A_log + B_log - 2 * C_log));
2783         } else {
2784             *epsi_x = 0.0;
2785             peak_act = -1;
2786             err_msg = "epsi_x = 0.0";
2787             flag = FALSE;
2788         }
2789
2790
2791     } else if (ifit == GPIV_POWER) {
2792 /*
2793  * sub-pixel estimator by quadratic fit
2794  */
2795         *epsi_y = (z[i_max[peak_act] - 1][j_max[peak_act]] -
2796                    z[i_max[peak_act] + 1][j_max[peak_act]]) /
2797             (2 * (z[i_max[peak_act] - 1][j_max[peak_act]] +
2798                   z[i_max[peak_act] + 1][j_max[peak_act]] -
2799                   2 * z[i_max[peak_act]][j_max[peak_act]]));
2800
2801         *epsi_x = (z[i_max[peak_act]][j_max[peak_act] - 1] -
2802                    z[i_max[peak_act]][j_max[peak_act] + 1]) /
2803             (2 * (z[i_max[peak_act]][j_max[peak_act] - 1] +
2804                   z[i_max[peak_act]][j_max[peak_act] + 1] -
2805                   2 * z[i_max[peak_act]][j_max[peak_act]]));
2806
2807
2808     } else if (ifit == GPIV_GRAVITY) {
2809 /*
2810  * sub-pixel estimator by centre of gravity
2811  */
2812         *epsi_y = (z[i_max[peak_act] - 1][j_max[peak_act]] -
2813                    z[i_max[peak_act] + 1][j_max[peak_act]]) /
2814             (z[i_max[peak_act] - 1][j_max[peak_act]] +
2815              z[i_max[peak_act] + 1][j_max[peak_act]] +
2816              z[i_max[peak_act]][j_max[peak_act]]);
2817
2818         *epsi_x = (z[i_max[peak_act]][j_max[peak_act] - 1] -
2819                    z[i_max[peak_act]][j_max[peak_act] + 1]) /
2820             (z[i_max[peak_act]][j_max[peak_act] - 1] +
2821              z[i_max[peak_act]][j_max[peak_act] + 1] +
2822              z[i_max[peak_act]][j_max[peak_act]]);
2823
2824
2825     } else {
2826         err_msg = "LIBGPIV internal error: cov_subtop: invalid fit parameter";
2827         gpiv_warning("%s", err_msg);
2828         return err_msg;
2829     }
2830
2831
2832     return err_msg;
2833 }
2834
2835
2836
2837 static int
2838 cov_top (GpivPivPar piv_par,
2839          GpivPivData * piv_data,
2840          int index_y,
2841          int index_x,
2842          GpivCov *cov,
2843          int x_corr,
2844          int ifit,
2845          int sweep,
2846          int last_sweep,
2847          int peak,
2848          int peak_act,
2849          int pre_shift_row_act,
2850          int pre_shift_col_act,
2851          int i_skip_act,
2852          int j_skip_act,
2853          gboolean *flag_subtop
2854          )
2855 /* ----------------------------------------------------------------------------
2856  * detects location of peak and snr in correlation function
2857  */
2858 {
2859 #define INITIAL_MIN 9999.9
2860     char *err_msg = NULL;
2861     float z_min, *z_max, *z_max_next;
2862     int h, i, j, i_min, j_min;
2863     long *i_max, *j_max, *i_max_next, *j_max_next;
2864
2865     int z_rl = cov->z_rl, z_rh = cov->z_rh, z_cl = cov->z_cl, z_ch = cov->z_ch;
2866
2867     int ipoint_x = (int) piv_data->point_x[index_y][index_x];
2868     int ipoint_y = (int) piv_data->point_y[index_y][index_x];
2869 /*     float epsi_x = 0.0, epsi_y = 0.0; */
2870     gboolean flag_snr = TRUE;
2871     gint dim = peak_act;
2872
2873
2874     i_max = gpiv_nulvector_index(1, dim + 1);
2875     j_max = gpiv_nulvector_index(1, dim + 1);
2876     z_max = gpiv_vector_index(1, dim + 1);
2877     i_max_next = gpiv_nulvector_index(1, dim + 2);
2878    j_max_next = gpiv_nulvector_index(1, dim + 2);
2879    z_max_next = gpiv_vector_index(1, dim + 2);
2880
2881 /*
2882  * BUGFIX: CHECK!!
2883  * finding a local top within the interrogation region. In case of
2884  * autocorrelation, exclude the first max (normally at i = 0,j = 0 if no
2885  * pre-shifting has been used), by increasing peak_act with 1 during the first
2886  * iteration sweep, then call it skip_act
2887  */
2888
2889     if (sweep == 0 && x_corr == 0) {
2890         peak_act = peak + 1;
2891     } else {
2892         if (x_corr == 0)
2893             peak_act = peak;
2894     }
2895
2896     search_top (cov, peak_act, x_corr, sweep, i_skip_act, j_skip_act,
2897                 z_max, i_max, j_max);
2898
2899     for (h = 1; h <= peak_act + 1; h++) {
2900         if (z_max_next[h] == -1.0) {
2901             ifit = GPIV_NONE;
2902             flag_snr = FALSE;
2903         }
2904     }
2905
2906 /*
2907  * Define first max to be skipped if autocorr, eventually shift this
2908  * point with new pre-shifting values
2909  */
2910
2911
2912     if (x_corr == 0 && sweep == 0) {
2913         i_skip_act = i_max[1];
2914         j_skip_act = j_max[1];
2915     }
2916
2917 /* BUGFIX: don't calculate snr for the Challenge project */
2918 /*     flag_snr = FALSE;   */
2919
2920 /*
2921  * Search next higher peak for SNR calculation
2922  */
2923     if (flag_snr) {
2924         search_top (cov, peak_act + 1, x_corr, sweep, i_skip_act, j_skip_act,
2925                     z_max_next, i_max_next, j_max_next);
2926     }
2927
2928 /*
2929  * Check if second top has been found
2930  */
2931     for (h = 1; h <= peak_act + 1; h++) {
2932         if (z_max_next[h] == -1.0) {
2933             flag_snr = FALSE;
2934         }
2935     }
2936
2937
2938     if (flag_snr
2939         && cov->z[i_max_next[peak_act + 1]][j_max_next[peak_act + 1]] != 0.0) {
2940         cov->snr = (cov->z[i_max[peak_act]][j_max[peak_act]] - cov->min) /
2941         (cov->z[i_max_next[peak_act + 1]][j_max_next[peak_act + 1]] - cov->min);
2942     } else {
2943         cov->snr = 0.0;
2944         piv_data->snr[index_y][index_x] = cov->snr;
2945 /*         piv_data->peak_no[index_y][index_x] = -1; */
2946     }
2947
2948 /*
2949  * Searching of minimum around cov. peak_act and remove 'pedestal'
2950  */
2951     z_min = INITIAL_MIN;
2952     i_min = INITIAL_MIN;
2953     j_min = INITIAL_MIN;
2954     for (i = i_max[peak_act] - 1; i <= i_max[peak_act] + 1; i++) {
2955         for (j = j_max[peak_act] - 1; j <= j_max[peak_act] + 1; j++) {
2956             if ((i >= z_rl && i <= z_rh) && (j >= z_cl && j <= z_ch)) {
2957                 if (cov->z[i][j] < z_min) {
2958                     z_min = cov->z[i][j];
2959                     i_min = i;
2960                     j_min = j;
2961                 }
2962             }
2963         }
2964     }
2965
2966     if (z_min <= INITIAL_MIN) {
2967         for (i = i_max[peak_act] - 1; i <= i_max[peak_act] + 1; i++) {
2968             for (j = j_max[peak_act] - 1; j <= j_max[peak_act] + 1; j++) {
2969 /*           cov->z[i][j] = cov->z[i][j]-z_min; */
2970                 cov->z[i][j] = cov->z[i][j] - z_min + 0.1;
2971             }
2972         }
2973     } else {
2974         ifit = GPIV_NONE;
2975     }
2976
2977 /*
2978  * Calculate particle displacement at integer pixel numbers or at sub-pixel
2979  */
2980     if (ifit == GPIV_NONE) {
2981         cov->subtop_x = 0.0;
2982         cov->subtop_y = 0.0;
2983
2984     } else {
2985             if ((err_msg = cov_subtop (cov->z, i_max, j_max, &cov->subtop_x,
2986                                        &cov->subtop_y, ifit, peak_act))
2987                 != NULL) {
2988                 g_message("%s", err_msg);
2989                 *flag_subtop = TRUE;
2990             }
2991     }
2992
2993 /*
2994  * Writing maximuma  to cov
2995  */
2996     cov->top_y = i_max[peak_act];
2997     cov->top_x = j_max[peak_act];
2998
2999 /*
3000  * Free memeory
3001  */
3002     gpiv_free_nulvector_index(i_max, 1, dim + 1);
3003     gpiv_free_nulvector_index(j_max, 1, dim + 1);
3004     gpiv_free_vector_index(z_max, 1, dim + 1);
3005     gpiv_free_nulvector_index(i_max_next, 1, dim + 2);
3006     gpiv_free_nulvector_index(j_max_next, 1, dim + 2);
3007     gpiv_free_vector_index(z_max_next, 1, dim + 2);
3008
3009     return 0;
3010 }
3011
3012
3013
3014 static
3015 void pack_cov (float **covariance,
3016                GpivCov *cov,
3017                int int_size_0
3018                )
3019 /*-----------------------------------------------------------------------------
3020  * Packs the unordered covariance data in an ordered sequence when returning
3021  * from cova_nr
3022  */
3023 {
3024     int i, j;
3025     int z_rnl = cov->z_rnl, z_rnh = cov->z_rnh, z_rpl = cov->z_rpl,
3026         z_rph = cov->z_rph;
3027     int z_cnl = cov->z_cnl, z_cnh = cov->z_cnh, z_cpl = cov->z_cpl,
3028         z_cph = cov->z_rph;
3029
3030
3031     for (i = z_rnl; i < z_rnh; i++) {
3032         for (j = z_cnl; j < z_cnh; j++) {
3033             cov->z[i - int_size_0][j - int_size_0] = covariance[i][j];
3034         }
3035         for (j = z_cpl; j < z_cph; j++) {
3036             cov->z[i - int_size_0][j] = covariance[i][j];
3037         }
3038     }
3039
3040     for (i = z_rpl; i < z_rph; i++) {
3041         for (j = z_cnl; j < z_cnh; j++) {
3042             cov->z[i][j - int_size_0] = covariance[i][j];
3043         }
3044         for (j = z_cpl; j < z_cph; j++) {
3045             cov->z[i][j] = covariance[i][j];
3046         }
3047     }
3048 }
3049
3050
3051
3052 static void
3053 weight_cov (GpivCov *cov,
3054             GpivCov *w_k
3055             )
3056 /*-----------------------------------------------------------------------------
3057  * Corrects ordered covariance data with weighting kernel
3058  */
3059 {
3060     int i, j;
3061     int z_rl = w_k->z_rl, z_rh = w_k->z_rh;
3062     int z_cl = w_k->z_cl, z_ch = w_k->z_ch;
3063
3064
3065     if (w_k == NULL) {
3066         g_message("LIBGPIV internal error: weight_cov: w_k = NULL");
3067         return;
3068     }
3069
3070     for (i = z_rl; i < z_rh; i++) {
3071         for (j = z_cl; j < z_ch; j++) {
3072             cov->z[i][j] = cov->z[i][j] / w_k->z[i][j];
3073         }
3074     }
3075
3076 }
3077
3078
3079 static gchar *
3080 filter_cov_spof (fftw_complex *a,
3081                  fftw_complex *b,
3082                  gint m,
3083                  gint n
3084                  )
3085 /*-----------------------------------------------------------------------------
3086  * Applies Symmetric Phase Only filtering on the complex arrays a and b
3087  *
3088  * REFERENCES:
3089  *     M.P. Wernet, "Symmetric phase only filtering: a new paradigm for DPIV
3090  *     data processing", Meas. Sci. Technol, vol 16 (2005), pp 601 - 618
3091  */
3092 {
3093     gchar *err_msg = NULL;
3094     gfloat amplitude_a, amplitude_b;
3095     gint i, j, ij;
3096
3097
3098     /*         assert (a[0] != NULL); */
3099     /*         assert (b[0] != NULL); */
3100
3101     for (i = 0; i < m; i++) {
3102         for (j = 0; j < n / 2 + 1; j++) {
3103             ij = i * (n / 2 + 1) + j;
3104             amplitude_a = sqrt(a[ij][0] * a[ij][0] + a[ij][1] * a[ij][1]);
3105             amplitude_b = sqrt(b[ij][0] * b[ij][0] + b[ij][1] * b[ij][1]);
3106
3107             if (amplitude_a == 0.0 || amplitude_b == 0.0) {
3108                 a[ij][0] = 0.0;
3109                 a[ij][1] = 0.0;
3110                 b[ij][0] = 0.0;
3111                 b[ij][1] = 0.0;
3112             } else {
3113                 a[ij][0] /= sqrt(amplitude_a) * sqrt(amplitude_b);
3114                 a[ij][1] /= sqrt(amplitude_a) * sqrt(amplitude_b);
3115                 b[ij][0] /= sqrt(amplitude_a) * sqrt(amplitude_b);
3116                 b[ij][1] /= sqrt(amplitude_a) * sqrt(amplitude_b);
3117             }
3118         }
3119     }
3120
3121
3122     return err_msg;
3123 }
3124
3125
3126
3127 static gchar *
3128 cova (const gboolean spof_filter,
3129       GpivFt *ft
3130       )
3131 /*-----------------------------------------------------------------------------
3132  * Calculates the covariance function of intreg1 and intreg2 by
3133  * means of Fast Fourier Transforms.
3134  */
3135 {
3136     gchar *err_msg = NULL;
3137     int M = ft->size, N = ft->size;
3138     float covariance_max, covariance_min;
3139     gint i, j;
3140     gint ij = 0;
3141     float **covariance;
3142     gdouble scale = 1.0 / (M * N);
3143     double *l_re = &ft->re[0][0];
3144     fftw_complex *l_co = &ft->co[0][0];
3145     fftw_complex *l_A = &ft->A[0][0];
3146     fftw_complex *l_B = &ft->B[0][0];
3147
3148
3149 #ifdef DEBUG
3150     FILE *fp = my_utils_fopen_tmp (GPIV_LOGFILE, "w");
3151 #endif
3152 #ifdef USE_MTRACE
3153     mtrace();
3154 #endif
3155
3156     gpiv_check_alloc_ft (ft);
3157
3158     covariance = gpiv_matrix(M, N);
3159     memset(covariance[0], 0, (sizeof(float)) * M * N);
3160
3161     /*
3162      * FFT both interrogation areas
3163      */
3164 #ifdef DEBUG
3165     fprintf (fp, "New I.A:\n");
3166 #endif
3167     /* copying intreg1 to re[][] for first FFT */
3168     for (i = 0; i < M; ++i) {
3169         for (j = 0; j < N; ++j) {
3170             ij = i * N + j;
3171             l_re[ij] = (double) ft->intreg1[i][j];
3172 #ifdef DEBUG
3173             fprintf (fp, "cova:: intreg1[%d][%d] = %f  re[%d] = %f\n",
3174                      i, j, ft->intreg1[i][j],
3175                      i, j, l_re[ij]);
3176 #endif /* DEBUG */
3177         }
3178     }
3179
3180 #ifdef DEBUG
3181     fprintf (fp, "\n\n");
3182 #endif
3183
3184     fftw_execute(ft->plan_forwardFFT);
3185
3186     /*
3187      * save first FFT result in A[][]
3188      */
3189     for (i = 0; i < M; ++i) {
3190         for (j = 0; j < (N/2+1); ++j) {
3191             ij = i * (N/2+1) + j;
3192             l_A[ij][0] = l_co[ij][0];   /* real part */
3193             l_A[ij][1] = l_co[ij][1];   /* imaginary part */
3194         }
3195     }
3196
3197
3198     /*
3199      * copying intreg2 to re[][] for second FFT
3200      */
3201     for (i = 0; i < M; ++i) {
3202         for (j = 0; j < N; ++j) {
3203             ij = i * N + j;
3204             l_re[ij] = (double) ft->intreg2[i][j];
3205         }
3206     }
3207
3208     fftw_execute(ft->plan_forwardFFT);
3209
3210     /*
3211      * save second FFT result in B[][]
3212      */
3213     for (i = 0; i < M; ++i) {
3214         for (j = 0; j < (N/2+1); ++j) {
3215             ij = i * (N/2+1) + j;
3216             l_B[ij][0] = l_co[ij][0];   /* real part */
3217             l_B[ij][1] = l_co[ij][1];   /* imaginary part */
3218         }
3219     }
3220
3221
3222     if (spof_filter) {
3223         if ((err_msg = filter_cov_spof(l_A, l_B, M, N)) != NULL) {
3224             return (err_msg);
3225         }
3226     }
3227
3228     /*
3229      * B * conjugate(A) result in correct sign of displacements!
3230      */
3231     /* copying B x A* to co[][] */
3232     for (i = 0; i < M; ++i) {
3233         for (j = 0; j < N / 2 + 1; ++j) {
3234             ij = i * (N / 2 + 1) + j;
3235             l_co[ij][0] =
3236                 (l_B[ij][0] * l_A[ij][0] +
3237                  l_B[ij][1] * l_A[ij][1]) * scale;
3238             l_co[ij][1] =
3239                 (l_B[ij][1] * l_A[ij][0] -
3240                  l_B[ij][0] * l_A[ij][1]) * scale;
3241
3242 #ifdef DEBUG
3243             fprintf (fp, "cova:: A[%d]_re = %f A[%d]_im = %f  B[%d]_re = %f B[%d]_im = %f\n",
3244                      ij, l_A[ij][0], ij, l_A[ij][1],
3245                      ij, l_B[ij][0], ij, l_B[ij][1]
3246                      );
3247 #endif
3248         }
3249     }
3250 #ifdef DEBUG
3251     fprintf (fp, "\n\n");
3252 #endif
3253
3254     /*
3255      * inverse transform to get the covariance of intreg1 and intreg2;
3256      * executing reverse-FFT on co[][], result in re[][]
3257      */
3258     fftw_execute(ft->plan_reverseFFT);
3259
3260
3261     /*
3262      * Put the data back in a 2-dim array covariance[][]
3263      * copying re[][] to covariance[][]
3264      */
3265    for (i = 0; i < M; i++) {
3266         for (j = 0; j < N; j++) {
3267             ij = i * N + j;
3268             covariance[i][j] = (float) l_re[ij];
3269
3270         }
3271     }
3272 #ifdef DEBUG
3273     fprintf (fp, "\n\n");
3274 #endif
3275
3276     /*
3277      * normalisation => correlation function
3278      */
3279     /* using system limits from float.h here */
3280     covariance_max = FLT_MIN;
3281     covariance_min = FLT_MAX;
3282     for (i = 0; i < M; i++) {
3283         for (j = 0; j < N; j++) {
3284             if (covariance[i][j] > covariance_max)
3285                 covariance_max = covariance[i][j];
3286             if (covariance[i][j] < covariance_min)
3287                 covariance_min = covariance[i][j];
3288         }
3289     }
3290
3291     for (i = 0; i < M; i++) {
3292         for (j = 0; j < N; j++) {
3293             covariance[i][j] = covariance[i][j] / covariance_max;
3294         }
3295     }
3296
3297
3298     /*
3299      * Packing the unordered covariance data into the ordered array of
3300      * Covariance structure
3301      */
3302     pack_cov(covariance, ft->cov, M);
3303     /* BUGFIX: may be changed */
3304     cov_min_max(ft->cov);
3305     /*     cov->min = covariance_min; */
3306     /*     cov->max = covariance_max; */
3307
3308
3309     /*
3310      * free mems
3311      */
3312     gpiv_free_matrix (covariance);
3313
3314     /*
3315      * REMARK: fftw_cleanup really slows down!
3316      */
3317     /*     fftw_forget_wisdom(); */
3318     /*     fftw_cleanup(); */
3319
3320
3321 #ifdef DEBUG
3322     fclose(fp);
3323 #endif
3324 #ifdef USE_MTRACE
3325     muntrace();
3326 #endif
3327     return err_msg;
3328 }
3329
3330
3331
3332 static gchar *
3333 ia_weight_gauss (gint int_size,
3334                  float **int_area
3335                  )
3336 /**----------------------------------------------------------------------------
3337  * DESCRIPTION:
3338  *     Weights Interrogation Area's
3339  */
3340 {
3341     gchar *err_msg = NULL;
3342     gint i = 0, j = 0;
3343     gdouble weight;
3344
3345
3346     assert (int_area[0] != NULL);
3347
3348     for (i = 0; i < int_size; i++) {
3349         for (j = 0; j < int_size; j++) {
3350             weight = exp( -16.0 * (((double)i - (double)int_size / 2.0)
3351                                     * ((double)i - (double)int_size / 2.0)
3352                                     + ((double)j - (double)int_size / 2.0)
3353                                     * ((double)j - (double)int_size / 2.0))
3354                           / ((double)int_size * (double)int_size));
3355             int_area[i][j] *= weight;
3356         }
3357     }
3358
3359
3360     return err_msg;
3361 }
3362
3363
3364 /*
3365  * From piv_speed
3366  */
3367 static void
3368 nearest_point (gint i,
3369                gfloat x,
3370                gfloat point_x,
3371                gfloat *min,
3372                gint *index,
3373                gboolean *exist
3374                )
3375 /*-----------------------------------------------------------------------------
3376  * Test if current point_x is nearest from x
3377  */
3378 {
3379     gfloat dif = abs (x - point_x);
3380
3381     if (dif < *min) {
3382         *exist = TRUE;
3383         *min = dif;
3384         *index = i;
3385     }
3386
3387 }
3388
3389
3390 static gboolean
3391 nearest_index (enum Position pos,
3392                gint vlength,
3393                gfloat *src_point,
3394                gfloat dest_point,
3395                gint *index
3396                )
3397 /*-----------------------------------------------------------------------------
3398  * Search nearest index from piv_data belonging to point (x, y)
3399  * in horizontal direction
3400  * pos_x/y index left/right, top/bottom of point
3401  */
3402 {
3403     gint i, j;
3404     gfloat min = 10.0e4;
3405     gboolean exist = FALSE;
3406
3407     *index = 0;
3408     for (i = 0; i < vlength; i++) {
3409
3410         if (pos ==  LOWER
3411             && src_point[i] <= dest_point) {
3412             nearest_point (i, dest_point, src_point[i], &min,
3413                            index, &exist);
3414
3415         } else if (pos == NEXT_LOWER
3416                    &&  i > 0
3417                    && src_point[i - 1] < dest_point) {
3418
3419             nearest_point (i - 1, dest_point, src_point[i - 1], &min,
3420                            index, &exist);
3421
3422         } else if (pos == HIGHER
3423                    && src_point[i] > dest_point) {
3424             nearest_point (i, dest_point, src_point[i], &min, index, &exist);
3425
3426         } else if (pos ==  NEXT_HIGHER
3427                    && i <  vlength - 1
3428                    && src_point[i + 1] > dest_point) {
3429             nearest_point (i + 1, dest_point, src_point[i + 1],
3430                            &min,
3431                            index, &exist);
3432
3433         }
3434
3435     }
3436
3437     return exist;
3438 }
3439
3440
3441
3442 static gdouble
3443 bilinear_interpol (gdouble alpha_hor,
3444                    gdouble alpha_vert,
3445                    gdouble src_dx_nw,
3446                    gdouble src_dx_ne,
3447                    gdouble src_dx_sw,
3448                    gdouble src_dx_se
3449                    )
3450 /*-----------------------------------------------------------------------------
3451  * Bilinear interpolation of src_dx_*
3452  * _ne: NORTH_EAST
3453  * _nw: NORTH_WEST
3454  * _se: SOUTH_EAST
3455  * _SW: SOUTH_WEST
3456  */
3457 {
3458     gdouble dx, dx_n, dx_s;
3459
3460
3461     dx_n = (1.0 - alpha_hor) * src_dx_nw + alpha_hor * src_dx_ne;
3462     dx_s = (1.0 - alpha_hor) * src_dx_sw + alpha_hor * src_dx_se;
3463     dx =  (1.0 - alpha_vert) * dx_n + alpha_vert * dx_s;
3464
3465
3466     return dx;
3467 }
3468
3469
3470
3471 static void *
3472 intpol_facts (gfloat *src_point,
3473               gfloat *dest_point,
3474               gint src_vlength,
3475               gint dest_vlength,
3476               gint *index_l,
3477               gint *index_h,
3478               gint *index_ll,
3479               gint *index_hh,
3480               double *alpha
3481               )
3482 /*-----------------------------------------------------------------------------
3483  * calculates normalized interpolation factors for piv_data_src
3484  * Think of:
3485  *          _l (LOWER) is used for _w (WEST) or _n (NORTH),
3486  *          _h (HIGHER) is used for _e (EAST) or _s (SOUTH)
3487  *          _ll (NEXT_LOWER) is used for _ww (WEST_WEST) or _nn (NORTH_NORTH),
3488  *          _hh (NEXT_HIGHER) is used for _ee (EAST_EAST) or _ss (SOUTH_SOUTH)
3489  */
3490 {
3491     gboolean *exist_l, *exist_h, *exist_ll, *exist_hh;
3492     double *dist_l, *dist_h, *dist_ll, *dist_hh;
3493     enum Position pos;
3494     gint i;
3495
3496
3497     exist_l = gpiv_gbolvector (dest_vlength);
3498     exist_h = gpiv_gbolvector (dest_vlength);
3499     exist_ll = gpiv_gbolvector (dest_vlength);
3500     exist_hh = gpiv_gbolvector (dest_vlength);
3501
3502     dist_l = gpiv_dvector (dest_vlength);
3503     dist_h = gpiv_dvector (dest_vlength);
3504     dist_ll = gpiv_dvector (dest_vlength);
3505     dist_hh = gpiv_dvector (dest_vlength);
3506
3507 /*
3508  * Searching adjacent and next to adjacent points of dest_point in src_point
3509  * data array
3510  */
3511     for (i = 0; i < dest_vlength; i++) {
3512         pos = LOWER;
3513         exist_l[i] = FALSE;
3514         if (exist_l[i] =
3515             nearest_index (pos,
3516                            src_vlength,
3517                            src_point,
3518                            dest_point[i],
3519                            &index_l[i])) {
3520             dist_l[i] = dest_point[i] - src_point[index_l[i]];
3521         } else {
3522 /*
3523  * To be used for extrapolation in negative direction
3524  * by applying higher and next to higher points of src data
3525  */
3526             pos = NEXT_HIGHER;
3527             exist_hh[i] = FALSE;
3528             if (exist_hh[i] =
3529                 nearest_index (pos,
3530                                src_vlength,
3531                                src_point,
3532                                dest_point[i],
3533                                &index_hh[i])) {
3534                 dist_hh[i] = dest_point[i] - src_point[index_hh[i]];
3535             }
3536         }
3537
3538
3539
3540         pos = HIGHER;
3541         exist_h[i] = FALSE;
3542         if (exist_h[i] =
3543             nearest_index (pos,
3544                            src_vlength,
3545                            src_point,
3546                            dest_point[i],
3547                            &index_h[i])) {
3548             dist_h[i] = dest_point[i] - src_point[index_h[i]];
3549         } else {
3550
3551 /*
3552  * To be used for extrapolation in positive direction
3553  * by applying lower and next to lower points of src data
3554  */
3555             pos = NEXT_LOWER;
3556             exist_ll[i] = FALSE;
3557             index_ll[i] = 0;
3558             if (exist_ll[i] =
3559                 nearest_index (pos,
3560                                src_vlength,
3561                                src_point,
3562                                dest_point[i],
3563                                &index_ll[i])) {
3564                 dist_ll[i] = dest_point[i] - src_point[index_ll[i]];
3565             }
3566         }
3567
3568 /*
3569  * calculating of weight factors for inter- or extrapolation
3570  */
3571
3572         if (exist_l[i] && exist_h[i]) {
3573 /*
3574  * Inner point: bilinear interpolation
3575  */
3576             if (src_point[index_l[i]] != src_point[index_h[i]]) {
3577                 alpha[i] = dist_l[i] /
3578                     (src_point[index_h[i]] - src_point[index_l[i]]);
3579             } else {
3580                 alpha[i] = 1.0;
3581             }
3582
3583
3584         } else if (exist_l[i] && exist_ll[i] && !exist_h[i]) {
3585 /*
3586  * extrapolation from two lower values
3587  */
3588             if (src_point[index_ll[i]] != src_point[index_l[i]]) {
3589                 alpha[i] = dist_ll[i] /
3590                     (src_point[index_l[i]] - src_point[index_ll[i]]);
3591                 index_h[i] = index_l[i];
3592                 index_l[i] = index_ll[i];
3593             } else {
3594                 alpha[i] = 1.0;
3595             }
3596
3597
3598         } else if (!exist_l[i] && exist_h[i] && exist_hh[i]) {
3599 /*
3600  * extrapolation from two higher values
3601  */
3602             if (src_point[index_hh[i]] != src_point[index_h[i]]) {
3603                 alpha[i] = dist_h[i] /
3604                     (src_point[index_hh[i]] - src_point[index_h[i]]);
3605                 index_l[i] = index_h[i];
3606                 index_h[i] = index_hh[i];
3607             } else {
3608                 alpha[i] = 1.0;
3609             }
3610
3611
3612         } else {
3613             alpha[i] = 1.0;
3614         }
3615     }
3616
3617
3618     gpiv_free_gbolvector (exist_l);
3619     gpiv_free_gbolvector (exist_h);
3620     gpiv_free_gbolvector (exist_ll);
3621     gpiv_free_gbolvector (exist_hh);
3622
3623     gpiv_free_dvector (dist_l);
3624     gpiv_free_dvector (dist_h);
3625     gpiv_free_dvector (dist_ll);
3626     gpiv_free_dvector (dist_hh);
3627
3628 }
3629
3630
3631 static void
3632 dxdy_at_new_grid_block (const GpivPivData *piv_data_src,
3633                         GpivPivData *piv_data_dest,
3634                         gint expansion_factor,
3635                         gint smooth_window
3636                         )
3637 /*-----------------------------------------------------------------------------
3638  * Calculates displacements from old to new grid, that has been expanded by
3639  * factor 2 and avarages with smoothing window. Works only correct if all neighbours
3640  * at equal distances
3641  */
3642 {
3643     int i, j, k, l, ef = expansion_factor, sw = smooth_window;
3644     int count = 0;
3645     GpivPivData *pd = NULL;
3646
3647     pd = gpiv_alloc_pivdata (piv_data_dest->nx, piv_data_dest->ny);
3648
3649 /*
3650  * Copy blocks of 2x2 input data to pd
3651  */
3652     for (i = 0; i < piv_data_src->ny; i++) {
3653         for (j = 0; j < piv_data_src->nx; j++) {
3654             for (k = 0; k < 2; k++) {
3655                 if (ef * i+k < pd->ny) {
3656                     for (l = 0; l < 2; l++) {
3657                         if (ef * j+l < pd->nx) {
3658                             pd->dx[ef * i+k][ef * j+l] = piv_data_src->dx[i][j];
3659                             pd->dy[ef * i+k][ef * j+l] = piv_data_src->dy[i][j];
3660                         }
3661                     }
3662                 }
3663             }
3664         }
3665     }
3666
3667 /*
3668  * smooth the data
3669  */
3670     for (i = 0; i < piv_data_src->ny; i++) {
3671         for (j = 0; j < piv_data_src->nx; j++) {
3672             count = 0;
3673             for (k = -sw + 1; k < sw; k++) {
3674                 if (i + k > 0 && i + k < pd->ny) {
3675                     for (l = -sw + 1; l < sw; l++) {
3676                         if (j + l > 0 && j + l < pd->ny) {
3677                             count++;
3678                             piv_data_dest->dx[i][j] += pd->dx[i+k][j+l];
3679                             piv_data_dest->dy[i][j] += pd->dy[i+k][j+l];
3680                         }
3681                     }
3682                 }
3683             }
3684             piv_data_dest->dx[i][j] = piv_data_dest->dx[i][j] / (float)count;
3685             piv_data_dest->dy[i][j] = piv_data_dest->dx[i][j] / (float)count;
3686         }
3687     }
3688
3689     gpiv_free_pivdata (pd);
3690 }
3691
3692
3693
3694 static gchar *
3695 update_pivdata_imgdeform_zoff (const GpivImage *image,
3696                                GpivImage *lo_image,
3697                                const GpivPivPar *piv_par,
3698                                const GpivValidPar *valid_par,
3699                                GpivPivData *piv_data,
3700                                GpivPivData *lo_piv_data,
3701                                GpivFt *ft,
3702                                gfloat *cum_residu,
3703                                gboolean *cum_residu_reached,
3704                                gfloat *sum_dxdy,
3705                                gfloat *sum_dxdy_old,
3706                                gboolean isi_last,
3707                                gboolean grid_last,
3708                                gboolean sweep_last,
3709                                gboolean verbose
3710                                )
3711 /*-----------------------------------------------------------------------------
3712  * Validates and updates / renews pivdata and some other variables when image
3713  * deformation or zero-offset interrogation scheme is used
3714  */
3715 {
3716     gchar *err_msg = NULL;
3717
3718
3719     /*
3720      * Test on outliers
3721      */
3722      if ((err_msg =
3723          gpiv_valid_errvec (lo_piv_data, image, piv_par, valid_par, ft, TRUE))
3724         != NULL) {
3725         return err_msg;
3726     }
3727
3728
3729     if (piv_par->int_scheme == GPIV_IMG_DEFORM) {
3730
3731         /*
3732          * Update PIV estimators with those from the last interrogation
3733          * Resetting local PIV estimators for eventual next interrogation
3734          */
3735         if ((err_msg = gpiv_add_dxdy_pivdata (lo_piv_data, piv_data))
3736             != NULL) {
3737             return err_msg;
3738         }
3739         if ((err_msg = gpiv_0_pivdata (lo_piv_data))
3740             != NULL) {
3741             return err_msg;
3742         }
3743
3744         /*
3745          * Deform image with updated PIV estimators.
3746          * First, copy local from original image.
3747          * Deform with newly updated PIV estimators.
3748          * Eventually write deformed image.
3749          */
3750
3751         if ((err_msg = gpiv_cp_img_data (image, lo_image))
3752             != NULL) {
3753             return err_msg;
3754         }
3755
3756         if ((err_msg = gpiv_imgproc_deform (lo_image, piv_data))
3757             != NULL) {
3758             return err_msg;
3759         }
3760
3761 /* #define DEBUG */
3762 #ifdef DEBUG
3763         if (sweep_last && verbose) {
3764             printf ("\n");
3765             if ((err_msg =
3766                  my_utils_write_tmp_image (lo_image, GPIV_DEFORMED_IMG_NAME,
3767                                            "Writing deformed image to:"))
3768                 != NULL) {
3769                 return err_msg;
3770             }
3771         }
3772 #endif /* DEBUG */
3773 /* #undef DEBUG */
3774
3775     } else {
3776
3777         /*
3778          * Renew PIV estimators with those from the last interrogation
3779          */
3780         if ((err_msg = gpiv_0_pivdata (piv_data))
3781             != NULL) {
3782             return err_msg;
3783         }
3784         if ((err_msg = gpiv_add_dxdy_pivdata (lo_piv_data, piv_data))
3785             != NULL) {
3786             return err_msg;
3787         }
3788     }
3789
3790     /*
3791      * Checking the relative cumulative residu for convergence
3792      * if final residu has been reached, cum_residu_reached will be set.
3793      */
3794     if (isi_last && grid_last) {
3795         *sum_dxdy_old = *sum_dxdy;
3796         *sum_dxdy = 0.0;
3797         gpiv_sum_dxdy_pivdata (piv_data, sum_dxdy);
3798         *cum_residu = fabsf ((*sum_dxdy - *sum_dxdy_old) /
3799                              ((gfloat)piv_data->nx * (gfloat)piv_data->ny));
3800         if (*cum_residu < GPIV_CUM_RESIDU_MIN) {
3801             *cum_residu_reached = TRUE;
3802         }
3803     }
3804
3805
3806     return err_msg;
3807 }
3808
3809
3810
3811 #undef NMIN_TO_INTERPOLATE
3812
3813
3814 #ifdef ENABLE_MPI
3815 static GpivPivData *
3816 piv_interrogate_img__scatterv_pivdata(GpivPivData *piv_data)
3817 /*---------------------------------------------------------------------------------------
3818  * Scatter the piv_data equally over its rows.
3819  *
3820  * The number of rows in piv_data need not be a multiple of nprocs.
3821  * Therefore, the first (piv_data->ny)%nprocs  processes get
3822  * ceil(piv_data->ny/np/nprocs) rows, and the remaining processes get
3823  * floor(piv_data->ny)/nprocs) rows.
3824  */
3825 {
3826     GpivPivData *pd = NULL;
3827     Mpi *mpi = NULL;
3828     gint i;
3829
3830
3831     mpi = alloc_mpi(piv_data);
3832     MPI_Comm_size(MPI_COMM_WORLD, mpi->nprocs);
3833     MPI_Comm_rank(MPI_COMM_WORLD, mpi->rank);
3834
3835 #ifdef DEBUG
3836     if (rank ==0) g_message ("gpiv_piv_interrogate_img:: nx=%d ny=%d nprocs = %d",
3837                              piv_data->nx, piv_data->ny, mpi->nprocs);
3838 #endif
3839
3840     gpiv_free_pivdata (piv_data);
3841
3842     for (i = 0; i < mpi->nprocs; i++) {
3843         if (mpi->rank == i) pd = gpiv_alloc_pivdata(mpi->piv_data->nx, mpi->counts[i] / mpi->piv_data->nx);
3844     }
3845
3846     gpiv_piv_mpi_scatterv_pivdata (mpi->piv_data, pd, mpi->counts, mpi->displs);
3847
3848
3849     free_mpi(mpi);
3850     return pd;
3851 }
3852
3853
3854
3855 static GpivPivData *
3856 piv_interrogate_img__gatherv_pivdata(GpivPivData *lo_piv_data,
3857                                      GpivPivData *piv_data)
3858 /*---------------------------------------------------------------------------------------
3859  * Gathers the piv_data equally over its rows.
3860  * Counterpart of piv_interrogate_img__scatterv_pivdata
3861  *
3862  * The number of rows in piv_data need not be a multiple of nprocs.
3863  * Therefore, the first (piv_data->ny)%nprocs  processes get
3864  * ceil(piv_data->ny/np/nprocs) rows, and the remaining processes get
3865  * floor(piv_data->ny)/nprocs) rows.
3866  */
3867 {
3868     GpivPivData *pd = NULL;
3869     Mpi *mpi = NULL;
3870
3871
3872     mpi = alloc_mpi(piv_data);
3873     gpiv_piv_mpi_gatherv_pivdata (lo_piv_data, mpi->piv_data, mpi->counts,
3874                                   mpi->displs);
3875     pd = gpiv_cp_pivdata (mpi->piv_data);
3876
3877
3878     free_mpi(mpi);
3879     return pd;
3880 }
3881
3882
3883
3884
3885 static Mpi *
3886 alloc_mpi(GpivPivData *piv_data)
3887 {
3888     Mpi *mpi = g_new0(Mpi, 1);
3889
3890
3891     mpi->counts = gpiv_piv_mpi_compute_counts(piv_data->nx, piv_data->ny);
3892     mpi->displs = gpiv_piv_mpi_compute_displs(mpi->counts, piv_data->nx,
3893                                               piv_data->ny);
3894     mpi->piv_data = gpiv_cp_pivdata (piv_data);
3895
3896
3897     return mpi;
3898 }
3899
3900
3901
3902 static void
3903 free_mpi(Mpi *mpi)
3904 {
3905     gpiv_free_pivdata (mpi->piv_data);
3906     gpiv_free_ivector (mpi->counts);
3907     gpiv_free_ivector (mpi->displs);
3908 }
3909
3910
3911
3912 guint
3913 substr_noremain(guint n,
3914                 guint m)
3915 /*-------------------------------------------------------------------
3916  * Substracts 1 while remainder of n not equal to zero
3917  */
3918 {
3919     while (fmod((double) n, (double) m) != 0) {
3920         n--;
3921     }
3922
3923
3924     return (guint) n;
3925 }
3926
3927
3928
3929 #endif /* ENABLE_MPI */
3930 gchar *
3931 gpiv_piv_gnuplot                        (const gchar            *title,
3932                                         const gfloat            gnuplot_scale,
3933                                         const gchar             *GNUPLOT_DISPLAY_COLOR,
3934                                         const guint             GNUPLOT_DISPLAY_SIZE,
3935                                         const GpivImagePar      *image_par,
3936                                         const GpivPivPar        *piv_par,
3937                                         const GpivPivData       *piv_data
3938                                         )
3939 /*-----------------------------------------------------------------------------
3940  * DESCRIPTION:
3941  *      Plots piv data as vectors on screen with gnuplot
3942  *
3943  * INPUTS:
3944  *      fname:         file name containing plot
3945  *      title:         title of plot
3946  *      gnuplot_scale: vector scale
3947  *      GNUPLOT_DISPLAY_COLOR:  display color of window containing graph
3948  *      GNUPLOT_DISPLAY_SIZE:   display size of window containing graph
3949  *      image_par:      image parameters
3950  *      piv_par:   piv evaluation parameters
3951  *      piv_data:       piv data
3952  *      RCSID:          program name and version that interrogated the image
3953  *
3954  * OUTPUTS:
3955  *
3956  * RETURNS:
3957  *     NULL on success or *err_msg on failure
3958  *---------------------------------------------------------------------------*/
3959 {
3960     gchar *err_msg = NULL;
3961     FILE *fp_cmd;
3962     const gchar *tmp_dir = g_get_tmp_dir ();
3963     gchar *fname_loc = "gpiv_gnuplot.cmd";
3964     gchar command[GPIV_MAX_CHARS];
3965     gchar fname_cmd[GPIV_MAX_CHARS];
3966     gint i, j;
3967
3968
3969     snprintf (fname_cmd, GPIV_MAX_CHARS, "%s/%s", tmp_dir, fname_loc);
3970
3971     if ((fp_cmd = fopen (fname_cmd, "w")) == NULL)
3972         gpiv_error ("gpiv_piv_gnuplot: error: Failure opening %s for output",
3973                     fname_cmd);
3974
3975     fprintf (fp_cmd, "set xlabel \"x (pixels)\"");
3976     fprintf (fp_cmd, "\nset ylabel \"y (pixels)\"");
3977     fprintf (fp_cmd, "\nset title \"Piv of %s\" ", title);
3978
3979     for (i = 0; i < piv_data->ny; i++) {
3980         for (j = 0; j < piv_data->nx; j++) {
3981             fprintf (fp_cmd, "\nset arrow from %f, %f to %f, %f",
3982                      piv_data->point_x[i][j],
3983                      piv_data->point_y[i][j],
3984                      piv_data->point_x[i][j] + piv_data->dx[i][j]
3985                      * gnuplot_scale,
3986                      piv_data->point_y[i][j] + piv_data->dy[i][j]
3987                      * gnuplot_scale);
3988         }
3989     }
3990
3991     fprintf (fp_cmd, "\nset nokey");
3992     if (piv_par->int_geo == GPIV_AOI) {
3993         fprintf (fp_cmd, "\nplot [%d:%d] [%d:%d] %d",
3994                  piv_par->col_start, piv_par->col_end,
3995                  piv_par->row_start, piv_par->row_end,
3996                  piv_par->row_end);
3997     } else {
3998         fprintf (fp_cmd, "\nplot [%d:%d] [%d:%d] %d",
3999                  0, image_par->ncolumns,
4000                  0, image_par->nrows,
4001                  piv_par->row_end);
4002
4003     }
4004     fprintf (fp_cmd, "\npause -1 \"Hit return to exit\"");
4005     fprintf (fp_cmd, "\nquit");
4006     fclose (fp_cmd);
4007
4008     snprintf (command, GPIV_MAX_CHARS,
4009               "gnuplot -bg %s -geometry %dx%d %s",
4010               GNUPLOT_DISPLAY_COLOR, GNUPLOT_DISPLAY_SIZE, GNUPLOT_DISPLAY_SIZE,
4011               fname_cmd);
4012
4013     if (system (command) != 0) {
4014         err_msg = "gpiv_piv_gnuplot: could not exec shell command";
4015         gpiv_warning ("%s", err_msg);
4016         return err_msg;
4017     }
4018
4019
4020     return err_msg;
4021 }
4022
4023
4024
4025 gchar *
4026 gpiv_histo_gnuplot                      (const gchar            *fname_out,
4027                                         const gchar             *title,
4028                                         const gchar             *GNUPLOT_DISPLAY_COLOR,
4029                                         const gint              GNUPLOT_DISPLAY_SIZE
4030                                         )
4031 /*-----------------------------------------------------------------------------
4032  * DESCRIPTION:
4033  *      Plots histogram on screen with gnuplot
4034  *
4035  * INPUTS:
4036  *      fname_out:      output filename
4037  *      title:          plot title
4038  *      GNUPLOT_DISPLAY_COLOR:  display color of window containing graph
4039  *      GNUPLOT_DISPLAY_SIZE:   display size of window containing graph
4040  *
4041  * OUTPUTS:
4042  *
4043  * RETURNS:
4044  *     NULL on success or *err_msg on failure
4045  *---------------------------------------------------------------------------*/
4046 {
4047     gchar *err_msg = NULL;
4048     FILE *fp_cmd;
4049     const gchar *tmp_dir = g_get_tmp_dir ();
4050     gchar *fname_loc = "gpiv_gnuplot.cmd";
4051     gchar *function_name = "gpiv_histo_gnuplot";
4052     gchar fname_cmd[GPIV_MAX_CHARS];
4053     gchar command[GPIV_MAX_CHARS];
4054
4055
4056     snprintf (fname_cmd, GPIV_MAX_CHARS, "%s/%s", tmp_dir, fname_loc);
4057     if ((fp_cmd = fopen (fname_cmd,"w")) == NULL) {
4058         err_msg = "GPIV_HISTO_GNUPLOT: Failure opening for output";
4059         gpiv_warning ("%s", err_msg);
4060         return err_msg;
4061     }
4062
4063     fprintf (fp_cmd, "set xlabel \"residu (pixels)\"");
4064     fprintf (fp_cmd, "\nset ylabel \"frequency\"");
4065     fprintf (fp_cmd, "\nplot \"%s\" title \"%s\" with boxes",
4066              fname_out, title);
4067
4068     fprintf (fp_cmd, "\npause -1 \"Hit return to exit\"");
4069     fprintf (fp_cmd, "\nquit");
4070
4071     fclose (fp_cmd);
4072
4073
4074     snprintf (command, GPIV_MAX_CHARS, "gnuplot -bg %s -geometry %dx%d %s",
4075               GNUPLOT_DISPLAY_COLOR, GNUPLOT_DISPLAY_SIZE,
4076               GNUPLOT_DISPLAY_SIZE, fname_cmd);
4077
4078     if (system (command) != 0) {
4079         g_warning ("%s:%s could not exec shell command",
4080                    LIBNAME, function_name);
4081
4082         exit (1);
4083     }
4084
4085
4086     return err_msg;
4087 }
4088
4089
4090
4091 void
4092 gpiv_histo                              (const GpivPivData      *data,
4093                                         GpivBinData             *klass
4094                                         )
4095 /*-----------------------------------------------------------------------------
4096  * DESCRIPTION:
4097  *      Calculates histogram from GpivPivData (NOT from ScalarData!!)
4098  *
4099  * INPUTS:
4100  *      data:           Input data
4101  *
4102  * OUTPUTS:
4103  *      klass:          Output data
4104  *
4105  * RETURNS:
4106  *---------------------------------------------------------------------------*/
4107 {
4108     gint i, j, k;
4109     gint nx = data->nx, ny = data->ny, **peak_no = data->peak_no;
4110     float **snr = data->snr;
4111     float delta;
4112     float *bound = klass->bound, *centre = klass->centre;
4113     gint *count = klass->count, nbins = klass->nbins;
4114
4115
4116     g_return_if_fail (data->point_x != NULL); /* gpiv_error ("data->point_x not allocated"); */
4117     g_return_if_fail (data->point_y != NULL); /* gpiv_error ("data->point_y not allocated"); */
4118     g_return_if_fail (data->dx != NULL); /* gpiv_error ("data->dx not allocated"); */
4119     g_return_if_fail (data->dy != NULL); /* gpiv_error ("data->dy not allocated"); */
4120     g_return_if_fail (data->snr != NULL); /* gpiv_error ("data->snr not allocated"); */
4121     g_return_if_fail (data->peak_no != NULL); /* gpiv_error ("ata->peak_no not allocated"); */
4122
4123     g_return_if_fail (klass->count != NULL); /* gpiv_error ("klass->count not allocated"); */
4124     g_return_if_fail (klass->bound != NULL); /* gpiv_error ("klass->bound not allocated"); */
4125     g_return_if_fail (klass->centre != NULL); /* gpiv_error ("klass->centre not allocated"); */
4126
4127     klass->min = 10.0e+9, klass->max = -10.0e+9;
4128     /*
4129      * find min and max value
4130      */
4131     for (i = 0; i < ny; i++) {
4132         for (j = 0; j < nx; j++) {
4133             if (peak_no[i][j] != -1) {
4134
4135                 if (snr[i][j] < klass->min)
4136                     klass->min = snr[i][j];
4137                 if (snr[i][j] >= klass->max)
4138                     klass->max = snr[i][j];
4139
4140             }
4141         }
4142     }
4143
4144
4145
4146     /*
4147      * Calculating boundaries of bins
4148      */
4149     delta = (klass->max - klass->min) / nbins;
4150     for (i = 0; i < nbins; i++) {
4151         centre[i] = klass->min + delta / 2.0 + (float) i *delta;
4152         count[i] = 0;
4153     }
4154
4155     for (i = 0; i < nbins; i++) {
4156         bound[i] = klass->min + (float) i * delta;
4157     }
4158
4159
4160
4161     /*
4162      * Sorting of snr data in bins
4163      */
4164     for (i = 0; i < ny; i++) {
4165         for (j = 0; j < nx; j++) {
4166             if (peak_no[i][j] != -1) {
4167
4168                 for (k = 0; k < nbins; k++) {
4169                     if ((snr[i][j] > bound[k])
4170                         && (snr[i][j] <= bound[k] + delta)) {
4171                         count[k] = count[k] + 1;
4172                     }
4173
4174                 }
4175             }
4176         }
4177     }
4178 }
4179
4180
4181
4182 void
4183 gpiv_cumhisto                           (const GpivPivData      *data,
4184                                         GpivBinData             *klass
4185                                         )
4186 /*-----------------------------------------------------------------------------
4187  * DESCRIPTION:
4188  *      Calculates cumulative histogram from GpivPivData (NOT from ScalarData!!)
4189  *
4190  * INPUTS:
4191  *      data:           Input data
4192  *
4193  * OUTPUTS:
4194  *      klass:          Output data
4195  *
4196  * RETURNS:
4197  *---------------------------------------------------------------------------*/
4198 {
4199     gint i, j, k;
4200     gint nx = data->nx, ny = data->ny, **peak_no = data->peak_no;
4201     float **snr = data->snr;
4202     float delta;
4203     float *bound = klass->bound, *centre = klass->centre;
4204     gint *count = klass->count, nbins = klass->nbins;
4205
4206
4207     g_return_if_fail (data->point_x != NULL); /* gpiv_error ("data->point_x not allocated"); */
4208     g_return_if_fail (data->point_y != NULL); /* gpiv_error ("data->point_y not allocated"); */
4209     g_return_if_fail (data->dx != NULL); /* gpiv_error ("data->dx not allocated"); */
4210     g_return_if_fail (data->dy != NULL); /* gpiv_error ("data->dy not allocated"); */
4211     g_return_if_fail (data->snr != NULL); /* gpiv_error ("data->snr not allocated"); */
4212     g_return_if_fail (data->peak_no != NULL); /* gpiv_error ("ata->peak_no not allocated"); */
4213
4214     g_return_if_fail (klass->count != NULL); /* gpiv_error ("klass->count not allocated"); */
4215     g_return_if_fail (klass->bound != NULL); /* gpiv_error ("klass->bound not allocated"); */
4216     g_return_if_fail (klass->centre != NULL); /* gpiv_error ("klass->centre not allocated"); */
4217
4218     klass->min = 10e+9, klass->max = -10e+9;
4219     /*
4220      * find min and max value
4221      */
4222     for (i = 0; i < ny; i++) {
4223         for (j = 0; j < nx; j++) {
4224             if (peak_no[i][j] != -1) {
4225
4226                 if (snr[i][j] < klass->min)
4227                     klass->min = snr[i][j];
4228                 if (snr[i][j] >= klass->max)
4229                     klass->max = snr[i][j];
4230
4231             }
4232         }
4233     }
4234
4235
4236     /*
4237      * Calculating boundaries of bins
4238      */
4239     delta = (klass->max - klass->min) / nbins;
4240     for (i = 0; i < nbins; i++) {
4241         centre[i] = klass->min + delta / 2.0 + (float) i *delta;
4242         count[i] = 0;
4243     }
4244
4245     for (i = 0; i < nbins; i++) {
4246         bound[i] = klass->min + (float) i * delta;
4247     }
4248
4249
4250     /*
4251      * Sorting of snr data in bins
4252      */
4253     for (i = 0; i < ny; i++) {
4254         for (j = 0; j < nx; j++) {
4255             if (peak_no[i][j] != -1) {
4256
4257                 for (k = 0; k < nbins; k++) {
4258                     if (snr[i][j] <= bound[k] + delta) {
4259                         count[k] = count[k] + 1;
4260                     }
4261
4262                 }
4263             }
4264         }
4265     }
4266 }