Merge ../gpivtools-omp
[gpivtools.git] / src / post / aint.c
blob8c86978c71be1c42770d68c83211b87c9800c33b
1 /* -*- Mode: C; indent-tabs-mode: nil; c-basic-offset: 4 c-style: "K&R" -*- */
3 /*---------------------------------------------------------------------------
5 aint - Calculates mean image intensity at each interrogation area
7 Copyright (C) 2006, 2007, 2008
8 Gerber van der Graaf <gerber_graaf@users.sourceforge.net
10 This program is free software; you can redistribute it and/or modify
11 it under the terms of the GNU General Public License as published by
12 the Free Software Foundation; either version 2, or (at your option)
13 any later version.
15 This program is distributed in the hope that it will be useful,
16 but WITHOUT ANY WARRANTY; without even the implied warranty of
17 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18 GNU General Public License for more details.
20 You should have received a copy of the GNU General Public License
21 along with this program; if not, write to the Free Software Foundation,
22 Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
24 ------------------------------------------------------------------------*/
26 #include <stdio.h>
27 #include <stdlib.h>
28 #include <glib.h>
29 #include <gpiv.h>
30 #include "config.h"
31 #include "git-rev.h"
33 #define PARFILE "gpivrc" /* Parameter file name */
34 #define USAGE "\
35 Usage: gpiv_aint [ -h | --help] [-p | --print] [-s | --second]\n\
36 [-V | --verbose] [-v | --version] filename \n\
37 \n\
38 keys: \n\
39 -h | --help: this on-line help \n\
40 -p | --print: print parameters to stdout \n\
41 -s | --second: obtain averages from second image frame \n\
42 -V | --verbose: prints data of interrogation area's and resulting mean \n\
43 -v | --version: version number \n\
44 filename: image file name \n\
47 #define HELP "\
48 gpiv_aint calculates mean image intensity at each interrogation area"
51 * Extension of filename means: "Interrogation Area Intensity"
53 #define AINT_EXT ".ain"
56 * Global variables
58 gboolean second_image = FALSE;
59 gboolean print_par = FALSE, verbose = FALSE;
61 #ifdef DEBUG
63 * Variables for development version
65 int print_main=0;
67 #endif
70 static gboolean
71 assign_img2intarr (gint ipoint_x,
72 gint ipoint_y,
73 guint16 **img_1,
74 guint16 **img_2,
75 gint int_size_f,
76 gint int_size_i,
77 gfloat **int_area_1,
78 gfloat **int_area_2,
79 gint pre_shift_row,
80 gint pre_shift_col,
81 gint nrows,
82 gint ncolumns,
83 gint int_size_0
86 static gboolean
87 assign_img2intarr_central (gint ipoint_x,
88 gint ipoint_y,
89 guint16 **img_1,
90 guint16 **img_2,
91 gint int_size_f,
92 gint int_size_i,
93 gfloat **int_area_1,
94 gfloat **int_area_2,
95 gint pre_shift_row,
96 gint pre_shift_col,
97 gint nrows,
98 gint ncolumns,
99 gint int_size_0
102 static gfloat
103 int_mean (gfloat **int_area,
104 gint int_size
107 static void
108 command_args (int argc,
109 char *argv[],
110 char fname[GPIV_MAX_CHARS]
113 static char *
114 make_fname (char *fname,
115 char *fname_par,
116 char *fname_piv,
117 char *fname_out
120 static GpivScalarData *
121 post_ia_intensity (GpivImage *image,
122 GpivPivData *piv_data,
123 GpivPivPar *piv_par
126 static GpivScalarData *
127 post_ia_intensity (GpivImage *image,
128 GpivPivData *piv_data,
129 GpivPivPar *piv_par
131 /* ----------------------------------------------------------------------------
132 * Calculates mean image intensity for each I.A.
135 GpivScalarData *sc_data = NULL;
137 gint i = 0, j = 0;
138 gfloat **int_area_1;
139 gfloat **int_area_2;
140 GpivImage *loc_image = NULL;
142 gboolean flag_accept = FALSE;
143 gint pre_shift_col_act = 0, pre_shift_row_act = 0;
146 assert (image->frame1[0] != NULL);
147 assert (image->frame2[0] != NULL);
149 assert (piv_data->point_x != NULL);
150 assert (piv_data->point_y != NULL);
151 assert (piv_data->dx != NULL);
152 assert (piv_data->dy != NULL);
153 assert (piv_data->snr != NULL);
154 assert (piv_data->peak_no != NULL);
156 int_area_1 = gpiv_matrix (piv_par->int_size_f, piv_par->int_size_f);
157 int_area_2 = gpiv_matrix (piv_par->int_size_f, piv_par->int_size_f);
160 sc_data = gpiv_alloc_scdata (piv_data->nx, piv_data->ny);
163 * Deform entire image
165 if (piv_par->int_scheme == GPIV_IMG_DEFORM) {
166 loc_image = gpiv_cp_img (image);
167 gpiv_imgproc_deform (loc_image, piv_data);
168 } else {
169 loc_image->frame1 = image->frame1;
170 loc_image->frame2 = image->frame2;
174 * For each I.A. obtain image data
175 * Source: gpiv_piv_interr_reg
177 if (piv_par->int_scheme == GPIV_ZERO_OFF_CENTRAL
178 || piv_par->int_scheme == GPIV_ZERO_OFF_FORWARD
181 #ifdef DEBUG
182 g_message("post_ia_intensity:: assign_img2intarr_central");
183 #endif /* DEBUG */
185 flag_accept = FALSE;
186 pre_shift_col_act = piv_data->dx[i][j] + piv_par->pre_shift_col;
187 pre_shift_row_act = piv_data->dy[i][j] + piv_par->pre_shift_row;
188 for (i = 0; i < piv_data->ny; i++) {
189 for (j = 0; j < piv_data->nx; j++) {
190 flag_accept =
191 assign_img2intarr_central ((int) piv_data->point_x[i][j],
192 (int) piv_data->point_y[i][j],
193 loc_image->frame1,
194 loc_image->frame2,
195 piv_par->int_size_f,
196 piv_par->int_size_f,
197 int_area_1,
198 int_area_2,
199 pre_shift_row_act,
200 pre_shift_col_act,
201 image->header->nrows,
202 image->header->ncolumns,
203 piv_par->int_size_f);
205 * Calculate average value and apply to out_data
207 sc_data->point_x[i][j] = piv_data->point_x[i][j];
208 sc_data->point_y[i][j] = piv_data->point_y[i][j];
209 if (flag_accept) {
210 sc_data->flag[i][j] = 1;
211 if (second_image) {
212 sc_data->scalar[i][j] = int_mean(int_area_2,
213 piv_par->int_size_f);
214 } else {
215 sc_data->scalar[i][j] = int_mean(int_area_1,
216 piv_par->int_size_f);
217 if (verbose) printf("x=%d y=%d mean = %f \n\n",
218 (int)sc_data->point_x[i][j],
219 (int)sc_data->point_x[i][j],
220 sc_data->scalar[i][j]);
222 } else {
223 sc_data->flag[i][j] = 0;
224 sc_data->scalar[i][j] = 0.0;
228 } else {
230 #ifdef DEBUG
231 g_message("post_ia_intensity:: assign_img2intarr");
232 #endif /* DEBUG */
233 flag_accept = FALSE;
234 pre_shift_col_act = piv_par->pre_shift_col;
235 pre_shift_row_act = piv_par->pre_shift_row;
236 for (i = 0; i < piv_data->ny; i++) {
237 for (j = 0; j < piv_data->nx; j++) {
238 flag_accept =
239 assign_img2intarr((int) piv_data->point_x[i][j],
240 (int) piv_data->point_y[i][j],
241 loc_image->frame1,
242 loc_image->frame2,
243 piv_par->int_size_f,
244 piv_par->int_size_f,
245 int_area_1,
246 int_area_2,
247 pre_shift_row_act,
248 pre_shift_col_act,
249 image->header->nrows,
250 image->header->ncolumns,
251 piv_par->int_size_f);
253 * Calculate average value and apply to out_data
255 sc_data->point_x[i][j] = piv_data->point_x[i][j];
256 sc_data->point_y[i][j] = piv_data->point_y[i][j];
257 if (flag_accept) {
258 sc_data->flag[i][j] = 1;
259 if (second_image) {
260 sc_data->scalar[i][j] = int_mean(int_area_2,
261 piv_par->int_size_f);
262 } else {
263 sc_data->scalar[i][j] = int_mean(int_area_1,
264 piv_par->int_size_f);
265 if (verbose) printf("x=%d y=%d mean = %f \n\n",
266 (int)sc_data->point_x[i][j],
267 (int)sc_data->point_x[i][j],
268 sc_data->scalar[i][j]);
270 } else {
271 sc_data->flag[i][j] = 0;
272 sc_data->scalar[i][j] = 0.0;
279 * Free allocated memory
281 gpiv_free_matrix(int_area_1);
282 gpiv_free_matrix(int_area_2);
283 if (piv_par->int_scheme == GPIV_IMG_DEFORM) {
284 gpiv_free_img (loc_image);
287 return sc_data;
291 static gboolean
292 assign_img2intarr(gint ipoint_x,
293 gint ipoint_y,
294 guint16 **img_1,
295 guint16 **img_2,
296 gint int_size_f,
297 gint int_size_i,
298 gfloat **int_area_1,
299 gfloat **int_area_2,
300 gint pre_shift_row,
301 gint pre_shift_col,
302 gint nrows,
303 gint ncolumns,
304 gint int_size_0
306 /*-----------------------------------------------------------------------------
307 * Assigns image data to the interrogation area arrays straightforward
310 gint m, n;
311 gint arg_int1_r = ipoint_y - int_size_f / 2 + 1;
312 gint arg_int1_c = ipoint_x - int_size_f / 2 + 1;
313 gint arg_int2_r = ipoint_y - int_size_i / 2 + 1;
314 gint arg_int2_c = ipoint_x - int_size_i / 2 + 1;
316 gboolean flag_valid;
318 #ifdef DEBUG
319 g_message("assign_img2intarr:: 0 arg_int1_r = %d arg_int1_c = %d",
320 arg_int1_r, arg_int1_c);
321 g_message("assign_img2intarr:: 0 arg_int2_r = %d arg_int2_c = %d",
322 arg_int2_r, arg_int2_c);
323 #endif /* DEBUG */
325 * Check if Interrogation Areas are within the image boundaries.
326 * Principally arg_int1_r,c don't have to be tested as
327 * int_size_i >= int_size_f, but has been kept to maintain generallity with the
328 * other assign_img2intarr* functions
330 if ((arg_int1_r) >= 0
331 && (arg_int1_r + int_size_f - 1) < nrows
332 && (arg_int1_c) >= 0
333 && (arg_int1_c + int_size_f - 1) < ncolumns
335 && (arg_int2_r) >= 0
336 && (arg_int2_r + int_size_i - 1) < nrows
337 && (arg_int2_c) >= 0
338 && (arg_int2_c + int_size_i - 1) < ncolumns) {
339 flag_valid = TRUE;
340 } else {
341 flag_valid = FALSE;
345 #ifdef DEBUG
346 g_message("assign_img2intarr:: 1 flag_valid=%d", flag_valid);
347 #endif /* DEBUG */
348 if (flag_valid == TRUE) {
350 * reset int_area_1, int_area_2 values
352 memset(int_area_1[0], 0.0, (sizeof(gfloat)) * int_size_0 * int_size_0);
353 memset(int_area_2[0], 0.0, (sizeof(gfloat)) * int_size_0 * int_size_0);
355 if (verbose) printf("assign_img2intarr: x = %d y = %d",
356 ipoint_x, ipoint_y);
357 for (m = 0; m < int_size_f; m++) {
358 if (verbose) printf("\n");
359 for (n = 0; n < int_size_f; n++) {
360 int_area_1[m][n] =
361 (float) img_1[m + arg_int1_r][n + arg_int1_c];
362 if (verbose) printf("%d ", (int) int_area_1[m][n]);
365 if (verbose) printf("\n\n");
367 for (m = 0; m < int_size_i; m++) {
368 for (n = 0; n < int_size_i; n++) {
369 int_area_2[m][n] =
370 (float) img_2[m + arg_int2_r][n + arg_int2_c];
375 #ifdef DEBUG
376 g_message("assign_img2intarr:: END");
377 #endif /* DEBUG */
378 return flag_valid;
383 static gboolean
384 assign_img2intarr_central(gint ipoint_x,
385 gint ipoint_y,
386 guint16 **img_1,
387 guint16 **img_2,
388 gint int_size_f,
389 gint int_size_i,
390 gfloat **int_area_1,
391 gfloat **int_area_2,
392 gint pre_shift_row,
393 gint pre_shift_col,
394 gint nrows,
395 gint ncolumns,
396 gint int_size_0
398 /*-----------------------------------------------------------------------------
399 * Assigns image data to the interrogation area arrays for central differential
400 * scheme
403 gint m, n;
404 gint idum = gpiv_max((int_size_i - int_size_f) / 2, 0);
405 gint arg_int1_r = ipoint_y - int_size_f / 2 + 1 - pre_shift_row / 2 -
406 pre_shift_row % 2;
407 gint arg_int1_c = ipoint_x - int_size_f / 2 + 1 - pre_shift_col / 2 -
408 pre_shift_col % 2;
409 gint arg_int2_r = ipoint_y - int_size_i / 2 + 1 + pre_shift_row / 2;
410 gint arg_int2_c = ipoint_x - int_size_i / 2 + 1 + pre_shift_col / 2;
411 gboolean flag_valid;
413 #ifdef DEBUG
414 g_message("assign_img2intarr_central:: 0 arg_int1_r = %d arg_int1_c = %d",
415 arg_int1_r, arg_int1_c);
416 g_message("assign_img2intarr_central:: 0 arg_int2_r = %d arg_int2_c = %d",
417 arg_int2_r, arg_int2_c);
418 #endif /* DEBUG */
420 * Check if Interrogation Areas are within the image boundaries.
422 if ((arg_int1_r) >= 0
423 && (arg_int1_r + int_size_f - 1) < nrows
424 && (arg_int1_c) >= 0
425 && (arg_int1_c + int_size_f - 1) < ncolumns
427 && (arg_int2_r) >= 0
428 && (arg_int2_r + int_size_i - 1) < nrows
429 && (arg_int2_c) >= 0
430 && (arg_int2_c + int_size_i - 1) < ncolumns) {
431 flag_valid = TRUE;
432 } else {
433 flag_valid = FALSE;
437 #ifdef DEBUG
438 g_message("assign_img2intarr_central:: 1 flag_valid=%d", flag_valid);
439 #endif /* DEBUG */
440 if (flag_valid == TRUE) {
442 * reset int_area_1, int_area_2 values
444 memset(int_area_1[0], 0.0, (sizeof(gfloat)) * int_size_0 * int_size_0);
445 memset(int_area_2[0], 0.0, (sizeof(gfloat)) * int_size_0 * int_size_0);
447 if (verbose) printf("assign_img2intarr_central: x = %d y = %d",
448 ipoint_x, ipoint_y);
449 for (m = 0; m < int_size_f; m++) {
450 if (verbose) printf("\n");
451 for (n = 0; n < int_size_f; n++) {
452 int_area_1[m + idum][n + idum] =
453 (float) img_1[m + arg_int1_r][n + arg_int1_c];
454 if (verbose) printf(" %d ", (int) int_area_1[m][n]);
457 if (verbose) printf("\n\n");
459 for (m = 0; m < int_size_i; m++) {
460 for (n = 0; n < int_size_i; n++) {
461 int_area_2[m][n] =
462 (float) img_2[m + arg_int2_r][n + arg_int2_c];
468 #ifdef DEBUG
469 g_message("assign_img2intarr_central:: END");
470 #endif /* DEBUG */
471 return flag_valid;
476 static gfloat
477 int_mean(gfloat **int_area,
478 gint int_size
480 /* ----------------------------------------------------------------------------
481 * calculates mean value from interrogation area
484 int m = 0, n = 0;
485 gfloat int_area_sum = 0;
486 gfloat mean = 0.0;
488 for (m = 0; m < int_size; m++) {
489 for (n = 0; n < int_size; n++) {
490 int_area_sum += int_area[m][n];
493 mean = int_area_sum / (int_size * int_size);
495 return mean;
500 static void
501 command_args (int argc,
502 char *argv[],
503 char fname[GPIV_MAX_CHARS]
505 /* ----------------------------------------------------------------------------
506 * Command line argument handling
509 char c;
510 int argc_next;
513 while (--argc > 0 && (*++argv)[0] == '-') {
514 argc_next=0;
516 * argc_next is set to 1 if the next cmd line argument has to be searched for;
517 * in case that the command line argument concerns more than one char or cmd
518 * line argument needs a parameter
520 while (argc_next == 0 && (c = *++argv[0])) {
521 switch (c) {
523 case 'v':
525 * Use Git revision control system
527 #ifdef GIT_HASH
528 printf ("git hash: %s\n", GIT_REV);
529 #else
530 printf ("version: %s\n", GPIVTOOLS_VERSION);
531 #endif
532 exit (0);
533 break;
535 case 'h':
536 printf ("%s\n", argv[0]);
537 printf ("%s\n", HELP);
538 printf ("%s\n", USAGE);
539 exit (0);
540 break;
542 case 'p':
543 print_par = TRUE;
544 break;
546 case 's':
547 second_image = TRUE;
548 break;
550 case 'V':
551 verbose = TRUE;
552 break;
555 * long option keys
557 case '-':
558 if (strcmp ("-help", *argv) == 0) {
559 printf ("\n%s", argv[0]);
560 printf ("\n%s", HELP);
561 printf ("\n%s", USAGE);
562 exit (0);
563 } else if (strcmp ("-print", *argv) == 0) {
564 print_par = TRUE;
565 } else if (strcmp ("-second", *argv) == 0) {
566 second_image = TRUE;
567 } else if (strcmp ("-verbose", *argv) == 0) {
568 verbose = TRUE;
569 } else if (strcmp ("-version", *argv) == 0) {
570 #ifdef GIT_HASH
571 printf ("git hash: %s\n", GIT_REV);
572 #else
573 printf ("version: %s\n", GPIVTOOLS_VERSION);
574 #endif
575 exit (0);
576 } else {
577 gpiv_error ("%s: unknown option: %s", argv[0], *argv);
579 argc_next = 1;
580 break;
582 default:
583 fprintf (stderr, USAGE);
584 exit (1);
585 break;
591 if (argc == 1) {
592 strcpy (fname, *argv);
593 } else {
594 gpiv_error ("%s: %s", argv[0], USAGE);
601 static gchar *
602 make_fname (char *fname,
603 char *fname_par,
604 char *fname_piv,
605 char *fname_out
607 /* ---------------------------------------------------------------------------
608 * function to generate filenames
611 gchar *err_msg = NULL;
612 gchar *fname_base = NULL;
614 if (fname == NULL ) {
615 err_msg = "make_fname: \"fname == NULL\"";
616 return (err_msg);
618 fname_base = g_strdup (fname);
619 strtok (fname_base, ".");
621 gpiv_io_make_fname (fname_base, GPIV_EXT_PAR, fname_par);
622 if (print_par) printf("# Parameter file: %s\n",fname_par);
624 gpiv_io_make_fname (fname_base, GPIV_EXT_PIV, fname_piv);
625 if (print_par) printf("# PIV data file: %s\n",fname_piv);
627 gpiv_io_make_fname (fname_base, AINT_EXT, fname_out);
628 if (print_par) printf("# Output file: %s\n",fname_out);
630 g_free (fname_base);
631 return (err_msg);
637 main (int argc,
638 char *argv[]
640 /* ----------------------------------------------------------------------------
641 * Start of the main program
644 gchar * err_msg = NULL;
645 FILE *fp, *fp_par;
646 gchar fname[GPIV_MAX_CHARS], fname_out[GPIV_MAX_CHARS],
647 fname_par[GPIV_MAX_CHARS], fname_piv[GPIV_MAX_CHARS];
649 GpivPivPar *piv_par = g_new0 (GpivPivPar, 1);
650 GpivPostPar *post_par = g_new0 (GpivPostPar, 1);
652 GpivImage *image = NULL;
653 GpivPivData *piv_in_data = NULL;
654 GpivScalarData *sc_out_data = NULL;
658 gpiv_piv_parameters_set (piv_par, FALSE);
659 gpiv_post_parameters_set (post_par, FALSE);
662 command_args (argc, argv, fname);
666 * Generating proper filenames
668 if ((err_msg =
669 make_fname (fname, fname_par, fname_piv, fname_out))
670 != NULL) {
671 gpiv_error ("%s: Failure calling make_fname", argv[0]);
676 * Prints command line parameters to par-file
678 if ((fp_par = fopen (fname_par, "a")) == NULL) {
679 gpiv_error ("\n%s: failure opening %s for input",
680 argv[0], fname_par);
682 #ifdef GIT_HASH
683 fprintf (fp_par, "\n\n# Software: %s\n# git hash: %s\n# Command line options:\n",
684 argv[0], GIT_REV);
685 #else
686 fprintf (fp_par, "\n\n# Software: %s\n# version: %s\n# Command line options:\n",
687 argv[0], GPIVTOOLS_VERSION);
688 #endif
691 * Reading parametes from PARFILE and (system) resources (and writing to data
692 * par-file)
695 gpiv_scan_parameter (GPIV_PIVPAR_KEY, PARFILE, piv_par, print_par);
696 if ((err_msg =
697 gpiv_scan_resourcefiles (GPIV_PIVPAR_KEY, piv_par, print_par))
698 != NULL) gpiv_error ("%s: %s", argv[0], err_msg);
699 gpiv_piv_print_parameters (fp_par, piv_par);
701 gpiv_scan_parameter (GPIV_POSTPAR_KEY, PARFILE, post_par, print_par);
702 if ((err_msg =
703 gpiv_scan_resourcefiles (GPIV_POSTPAR_KEY, post_par, print_par))
704 != NULL) gpiv_error ("%s: %s", argv[0], err_msg);
705 gpiv_post_print_parameters (fp_par, post_par);
708 if ((err_msg =
709 gpiv_post_check_parameters_read (post_par, NULL))
710 != NULL) gpiv_error ("%s: %s", argv[0], err_msg);
713 * Read input image
715 if ((image = gpiv_fread_image (fname)) == NULL) {
716 gpiv_error ("%s: gpiv_fread_image_ni", argv[0]);
721 * Read PIV data for input
723 if ((fp = fopen (fname_piv, "rb")) == NULL) {
724 guint nx, ny;
726 fprintf (fp_par, "# WARNING: no PIV data found; particle displacements set to 0");
727 if (print_par) gpiv_warning ("Particle displacements set to 0");
729 if ((err_msg =
730 gpiv_piv_count_pivdata_fromimage (image->header, piv_par, &nx, &ny))
731 != NULL) gpiv_error ("%s: %s", argv[0], err_msg);
733 if ((piv_in_data =
734 gpiv_piv_gridgen (nx, ny, image->header, piv_par))
735 == NULL) gpiv_error ("%: %s", argv[0], err_msg);
737 if ((err_msg =
738 gpiv_0_pivdata (piv_in_data))
739 != NULL) gpiv_error (": %s", argv[0], err_msg);
741 } else {
742 if ((piv_in_data = gpiv_read_pivdata (fp))
743 == NULL) gpiv_error ("%s: %s", argv[0], err_msg);
744 fclose(fp);
747 fclose(fp_par);
751 * Here the function calls of the post-processing; intensity of
752 * Interrogation Area's
754 if (print_par == TRUE) printf ("\n");
756 if ((sc_out_data =
757 post_ia_intensity (image, piv_in_data, piv_par))
758 == NULL) gpiv_error ("%s: Failure %s", argv[0], err_msg);
761 * And writing to output
763 if ((fp = fopen (fname_out, "wb")) == NULL) {
764 gpiv_error ("%s: failing opening %s", argv[0], fname_out);
766 if ((err_msg =
767 gpiv_write_scdata (fp, sc_out_data, TRUE))
768 != NULL) gpiv_error ("%s: %s", argv[0], err_msg);
769 fclose (fp);
772 * Freeing allocated memory
774 gpiv_free_pivdata (piv_in_data);
775 gpiv_free_img (image);
776 if (print_par == TRUE) printf ("\n");
777 exit (0);