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