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