lib/imgproc_deform.c

   1 /* -*- Mode: C; indent-tabs-mode: nil; c-basic-offset: 4 c-style: "K&R" -*- */
   2
   3 /*
   4    libgpiv - library for Particle Image Velocimetry
   5
   6    Copyright (C) 2002, 2003, 2004, 2005, 2006, 2007, 2008
   7    Gerber van der Graaf
   8
   9    This file is part of libgpiv.
  10    Libgpiv 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.
  14
  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.
  19
  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.
  23
  24
  25
  26 -------------------------------------------------------------------------------
  27 FILENAME:                imgproc_deform.c
  28 LIBRARY:                 libgpiv:
  29 EXTERNAL FUNCTIONS:      gpiv_imgproc_deform
  30
  31 LOCAL FUNCTIONS:         compare_float
  32
  33 LAST MODIFICATION DATE:  $Id: imgproc_deform.c,v 1.7 2008-09-25 13:19:53 gerber Exp $
  34
  35 -------------------------------------------------------------------------------
  36
  37
  38   This software is provided by P. Thévenaz and based on the following paper:
  39
  40   P. Thévenaz, T. Blu, M. Unser, "Interpolation Revisited," IEEE
  41   Transactions on Medical Imaging, vol. 19, no. 7, pp. 739-758, July 2000.
  42
  43   EPFL/STI/IOA/BIG
  44   Philippe Thevenaz
  45   Bldg. BM-Ecublens 4.137
  46   CH-1015 Lausanne
  47   phone (CET):  +41(21)693.51.61
  48   fax:                  +41(21)693.37.01
  49   RFC-822:              philippe.thevenaz@epfl.ch
  50   X-400:                /C=ch/A=400net/P=switch/O=epfl/S=thevenaz/G=philippe/
  51   URL:                  http://bigwww.epfl.ch/
  52 ------------------------------------------------------------------------------*/
  53
  54 #include        <float.h>
  55 /* #include     <math.h> */
  56 /* #include     <stddef.h> */
  57 /* #include     <stdio.h> */
  58 /* #include     <stdlib.h> */
  59
  60 #include <stddef.h>
  61 #include <gpiv.h>
  62
  63 /*
  64  * Prototyping of static functions
  65  */
  66 static void
  67 ConvertToInterpolationCoefficients
  68                                 (double c[],            /* input samples --> output coefficients */
  69                                 long    DataLength,     /* number of samples or coefficients */
  70                                 double  z[],            /* poles */
  71                                 long    NbPoles,        /* number of poles */
  72                                 double  Tolerance       /* admissible relative error */
  73                                 );
  74
  75 static void
  76 GetColumn                       (float  *Image,         /* input image array */
  77                                 long    Width,          /* width of the image */
  78                                 long    x,              /* x coordinate of the selected line */
  79                                 double  Line[],         /* output linear array */
  80                                 long    Height          /* length of the line and height of the image */
  81                                 );
  82
  83 static void
  84 GetRow                          (float  *Image,         /* input image array */
  85                                 long    y,              /* y coordinate of the selected line */
  86                                 double  Line[],         /* output linear array */
  87                                 long    Width           /* length of the line and width of the image */
  88                                 );
  89
  90 static double
  91 InitialCausalCoefficient        (double c[],            /* coefficients */
  92                                 long    DataLength,     /* number of coefficients */
  93                                 double  z,              /* actual pole */
  94                                 double  Tolerance       /* admissible relative error */
  95                                 );
  96
  97 static double
  98 InitialAntiCausalCoefficient    (double c[],            /* coefficients */
  99                                 long    DataLength,     /* number of samples or coefficients */
 100                                 double  z               /* actual pole */
 101                                 );
 102
 103 static void
 104 PutColumn                       (float  *Image,         /* output image array */
 105                                 long    Width,          /* width of the image */
 106                                 long    x,              /* x coordinate of the selected line */
 107                                 double  Line[],         /* input linear array */
 108                                 long    Height          /* length of the line and height of the image */
 109                                 );
 110
 111 static void
 112 PutRow                          (float  *Image,         /* output image array */
 113                                 long    y,              /* y coordinate of the selected line */
 114                                 double  Line[],         /* input linear array */
 115                                 long    Width           /* length of the line and width of the image */
 116                                 );
 117
 118
 119
 120 static int
 121 SamplesToCoefficients           (float  *Image,         /* in-place processing */
 122                                 long    Width,          /* width of the image */
 123                                 long    Height,         /* height of the image */
 124                                 long    SplineDegree    /* degree of the spline model */
 125                                 );
 126
 127 static gchar *
 128 imgproc_deform_frame_from_pivdata(guint16 **img,
 129                                   const GpivImagePar *image_par,
 130                                   const GpivPivData *piv_data,
 131                                   const gdouble magnitude
 132                                   );
 133
 134 static int
 135 imgproc_deform_frame            (guint16                **img,
 136                                  const GpivImagePar     *image_par,
 137                                  const GpivPivData      *piv_data
 138                                 );
 139
 140 static double
 141 InterpolatedValue               (float  *Bcoeff,        /* input B-spline array of coefficients */
 142                                 const long      Width,          /* width of the image */
 143                                 const long      Height,         /* height of the image */
 144                                 const double    x,              /* x coordinate where to interpolate */
 145                                 const double    y,              /* y coordinate where to interpolate */
 146                                 const long      SplineDegree    /* degree of the spline model */
 147                                 );
 148
 149 /*
 150  * Public functions
 151  */
 152 gchar *
 153 gpiv_imgproc_deform (GpivImage *image,
 154                      const GpivPivData *piv_data
 155                      )
 156 /*-----------------------------------------------------------------------------
 157  * DESCRIPTION:
 158  *     Image shifting and deformation routine for a single exposed, double
 159  *     frame PIV image pair with magnitude of PIV estimations at each pixel
 160  *     Deforms first frame halfway forward and second frame halfway backward.
 161  *
 162  *---------------------------------------------------------------------------*/
 163 {
 164     gchar *err_msg = NULL;
 165
 166
 167     if (image->header->x_corr__set == FALSE
 168         || image->header->x_corr == FALSE) {
 169         err_msg = "gpiv_imgproc_deform: only to be used for cross correlation";
 170         return (err_msg);
 171     }
 172
 173     if (image->frame1[0] == NULL) {
 174         err_msg = "gpiv_imgproc_deform: image->frame1 == NULL";
 175         return (err_msg);
 176     }
 177
 178     if (image->frame2[0] == NULL) {
 179         err_msg = "gpiv_imgproc_deform: image->frame2 == NULL";
 180         return (err_msg);
 181     }
 182
 183     /*
 184      * Deform first image frame with half shift forwards
 185      * and second frame with half shift backwards.
 186      */
 187     /* OK: results approved */
 188 #pragma omp parallel sections
 189     {
 190 #pragma omp section
 191         {
 192             if ((err_msg =
 193                  imgproc_deform_frame_from_pivdata(image->frame1, image->header,
 194                                                    piv_data, 0.5))
 195                 != NULL) {
 196 #ifdef ENABLE_OMP
 197                 exit;
 198 #else
 199                 return (err_msg);
 200 #endif
 201             }
 202         }
 203
 204 #pragma omp section
 205         {
 206             if ((err_msg =
 207                  imgproc_deform_frame_from_pivdata(image->frame2, image->header,
 208                                                    piv_data, -0.5))
 209                 != NULL) {
 210 #ifdef ENABLE_OMP
 211                 exit;
 212 #else
 213                 return (err_msg);
 214 #endif
 215             }
 216         }
 217     }
 218
 219
 220     return (err_msg);
 221 }
 222
 223
 224
 225 /*
 226  * Definition of static functions
 227  */
 228
 229 static void
 230 ConvertToInterpolationCoefficients (
 231                                     double      c[],            /* input samples --> output coefficients */
 232                                     long        DataLength,     /* number of samples or coefficients */
 233                                     double      z[],            /* poles */
 234                                     long        NbPoles,        /* number of poles */
 235                                     double      Tolerance       /* admissible relative error */
 236                                     )
 237 /*--------------------------------------------------------------------------*/
 238 {
 239         double  Lambda = 1.0;
 240         long    n, k;
 241         /* special case required by mirror boundaries */
 242         if (DataLength == 1L) {
 243                 return;
 244         }
 245         /* compute the overall gain */
 246         for (k = 0L; k < NbPoles; k++) {
 247                 Lambda = Lambda * (1.0 - z[k]) * (1.0 - 1.0 / z[k]);
 248         }
 249         /* apply the gain */
 250         for (n = 0L; n < DataLength; n++) {
 251                 c[n] *= Lambda;
 252         }
 253         /* loop over all poles */
 254         for (k = 0L; k < NbPoles; k++) {
 255                 /* causal initialization */
 256                 c[0] = InitialCausalCoefficient(c, DataLength, z[k], Tolerance);
 257                 /* causal recursion */
 258                 for (n = 1L; n < DataLength; n++) {
 259                         c[n] += z[k] * c[n - 1L];
 260                 }
 261                 /* anticausal initialization */
 262                 c[DataLength - 1L] = InitialAntiCausalCoefficient(c, DataLength, z[k]);
 263                 /* anticausal recursion */
 264                 for (n = DataLength - 2L; 0 <= n; n--) {
 265                         c[n] = z[k] * (c[n + 1L] - c[n]);
 266                 }
 267         }
 268 }
 269
 270
 271
 272 static double
 273 InitialCausalCoefficient (
 274                           double        c[],    /* coefficients */
 275                           long          DataLength,     /* number of coefficients */
 276                           double        z,      /* actual pole */
 277                           double        Tolerance       /* admissible relative error */
 278                           )
 279 /*--------------------------------------------------------------------------*/
 280 {
 281         double  Sum, zn, z2n, iz;
 282         long    n, Horizon;
 283         /* this initialization corresponds to mirror boundaries */
 284         Horizon = DataLength;
 285         if (Tolerance > 0.0) {
 286                 Horizon = (long)ceil(log(Tolerance) / log(fabs(z)));
 287         }
 288         if (Horizon < DataLength) {
 289                 /* accelerated loop */
 290                 zn = z;
 291                 Sum = c[0];
 292                 for (n = 1L; n < Horizon; n++) {
 293                         Sum += zn * c[n];
 294                         zn *= z;
 295                 }
 296                 return(Sum);
 297         }
 298         else {
 299                 /* full loop */
 300                 zn = z;
 301                 iz = 1.0 / z;
 302                 z2n = pow(z, (double)(DataLength - 1L));
 303                 Sum = c[0] + z2n * c[DataLength - 1L];
 304                 z2n *= z2n * iz;
 305                 for (n = 1L; n <= DataLength - 2L; n++) {
 306                         Sum += (zn + z2n) * c[n];
 307                         zn *= z;
 308                         z2n *= iz;
 309                 }
 310                 return(Sum / (1.0 - zn * zn));
 311         }
 312 }
 313
 314
 315
 316 static void
 317 GetColumn (
 318            float        *Image,         /* input image array */
 319            long         Width,          /* width of the image */
 320            long         x,              /* x coordinate of the selected line */
 321            double       Line[],         /* output linear array */
 322            long         Height          /* length of the line */
 323            )
 324 /*--------------------------------------------------------------------------*/
 325 { /* begin GetColumn */
 326         long    y;
 327         Image = Image + (ptrdiff_t)x;
 328         for (y = 0L; y < Height; y++) {
 329                 Line[y] = (double)*Image;
 330                 Image += (ptrdiff_t)Width;
 331         }
 332 } /* end GetColumn */
 333
 334
 335
 336 static void
 337 GetRow (
 338         float   *Image,         /* input image array */
 339         long    y,              /* y coordinate of the selected line */
 340         double  Line[],         /* output linear array */
 341         long    Width           /* length of the line */
 342         )
 343 /*--------------------------------------------------------------------------*/
 344 {
 345         long    x;
 346         Image = Image + (ptrdiff_t)(y * Width);
 347         for (x = 0L; x < Width; x++) {
 348                 Line[x] = (double)*Image++;
 349         }
 350 }
 351
 352
 353
 354 static double
 355 InitialAntiCausalCoefficient (
 356                               double    c[],            /* coefficients */
 357                               long      DataLength,     /* number of samples or coefficients */
 358                               double    z               /* actual pole */
 359                               )
 360 /*--------------------------------------------------------------------------*/
 361 {
 362         /* this initialization corresponds to mirror boundaries */
 363         return((z / (z * z - 1.0)) * (z * c[DataLength - 2L] + c[DataLength - 1L]));
 364 }
 365
 366
 367
 368 static void
 369 PutColumn (
 370            float        *Image,         /* output image array */
 371            long Width,                  /* width of the image */
 372            long x,                      /* x coordinate of the selected line */
 373            double       Line[],         /* input linear array */
 374            long Height                  /* length of the line and height of the image */
 375            )
 376 /*--------------------------------------------------------------------------*/
 377 {
 378         long    y;
 379         Image = Image + (ptrdiff_t)x;
 380         for (y = 0L; y < Height; y++) {
 381                 *Image = (float)Line[y];
 382                 Image += (ptrdiff_t)Width;
 383         }
 384 }
 385
 386
 387
 388 static void
 389 PutRow (
 390         float   *Image,         /* output image array */
 391         long    y,              /* y coordinate of the selected line */
 392         double  Line[],         /* input linear array */
 393         long    Width           /* length of the line and width of the image */
 394         )
 395 /*--------------------------------------------------------------------------*/
 396 {
 397         long    x;
 398         Image = Image + (ptrdiff_t)(y * Width);
 399         for (x = 0L; x < Width; x++) {
 400                 *Image++ = (float)Line[x];
 401         }
 402 }
 403
 404
 405
 406 static int
 407 SamplesToCoefficients (
 408                        float    *Image,         /* in-place processing */
 409                        long     Width,          /* width of the image */
 410                        long     Height,         /* height of the image */
 411                        long     SplineDegree    /* degree of the spline model */
 412                        )
 413 /*--------------------------------------------------------------------------*/
 414 {
 415         double  *Line;
 416         double  Pole[2];
 417         long    NbPoles;
 418         long    x, y;
 419         /* recover the poles from a lookup table */
 420
 421         switch (SplineDegree) {
 422                 case 2L:
 423                         NbPoles = 1L;
 424                         Pole[0] = sqrt(8.0) - 3.0;
 425                         break;
 426                 case 3L:
 427                         NbPoles = 1L;
 428                         Pole[0] = sqrt(3.0) - 2.0;
 429                         break;
 430                 case 4L:
 431                         NbPoles = 2L;
 432                         Pole[0] = sqrt(664.0 - sqrt(438976.0)) + sqrt(304.0) - 19.0;
 433                         Pole[1] = sqrt(664.0 + sqrt(438976.0)) - sqrt(304.0) - 19.0;
 434                         break;
 435                 case 5L:
 436                         NbPoles = 2L;
 437                         Pole[0] = sqrt(135.0 / 2.0 - sqrt(17745.0 / 4.0)) + sqrt(105.0 / 4.0)
 438                                 - 13.0 / 2.0;
 439                         Pole[1] = sqrt(135.0 / 2.0 + sqrt(17745.0 / 4.0)) - sqrt(105.0 / 4.0)
 440                                 - 13.0 / 2.0;
 441                         break;
 442                 default:
 443                         printf("Invalid spline degree\n");
 444                         return(1);
 445         }
 446
 447 /*
 448  * convert the image samples into interpolation coefficients
 449  * in-place separable process, along x
 450  */
 451         Line = (double *)malloc((size_t)(Width * (long)sizeof(double)));
 452         if (Line == (double *)NULL) {
 453                 printf("Row allocation failed\n");
 454                 return(1);
 455         }
 456
 457
 458 /* BUGFIX SamplesToCoefficients: OMP causes error */
 459 /* KEEP OUT! #pragma omp parallel for */
 460         for (y = 0L; y < Height; y++) {
 461                 GetRow(Image, y, Line, Width);
 462                 ConvertToInterpolationCoefficients(Line, Width, Pole, NbPoles,
 463                                                    DBL_EPSILON);
 464                 PutRow(Image, y, Line, Width);
 465         }
 466
 467
 468         free(Line);
 469         /* in-place separable process, along y */
 470         Line = (double *)malloc((size_t)(Height * (long)sizeof(double)));
 471         if (Line == (double *)NULL) {
 472                 printf("Column allocation failed\n");
 473                 return(1);
 474         }
 475
 476
 477 /* BUGFIX SamplesToCoefficients: OMP causes error */
 478 /* KEEP OUT! #pragma omp parallel for */
 479         for (x = 0L; x < Width; x++) {
 480                 GetColumn(Image, Width, x, Line, Height);
 481                 ConvertToInterpolationCoefficients(Line, Height, Pole, NbPoles, DBL_EPSILON);
 482                 PutColumn(Image, Width, x, Line, Height);
 483         }
 484
 485
 486         free(Line);
 487         return(0);
 488 }
 489
 490
 491
 492 static double
 493 InterpolatedValue(
 494                   float *Bcoeff,        /* input B-spline array of coefficients */
 495                   const long   Width,           /* width of the image */
 496                   const long   Height,  /* height of the image */
 497                   const double x,               /* x coordinate where to interpolate */
 498                   const double y,               /* y coordinate where to interpolate */
 499                   const long SplineDegree     /* degree of the spline model */
 500                   )
 501 /*--------------------------------------------------------------------------*/
 502 {
 503     float       *p;
 504     double      xWeight[6], yWeight[6];
 505     double      interpolated;
 506     double      w, w2, w4, t, t0, t1;
 507     long        xIndex[6], yIndex[6];
 508     long        Width2 = 2L * Width - 2L, Height2 = 2L * Height - 2L;
 509     long        i, j, k;
 510
 511 /*
 512  * compute the interpolation indexes
 513  */
 514     if (SplineDegree & 1L) {
 515         i = (long)floor(x) - SplineDegree / 2L;
 516         j = (long)floor(y) - SplineDegree / 2L;
 517         for (k = 0L; k <= SplineDegree; k++) {
 518             xIndex[k] = i++;
 519             yIndex[k] = j++;
 520         }
 521     } else {
 522         i = (long)floor(x + 0.5) - SplineDegree / 2L;
 523         j = (long)floor(y + 0.5) - SplineDegree / 2L;
 524         for (k = 0L; k <= SplineDegree; k++) {
 525             xIndex[k] = i++;
 526             yIndex[k] = j++;
 527         }
 528     }
 529
 530 /*
 531  * compute the interpolation weights
 532  */
 533     switch (SplineDegree) {
 534     case 2L:
 535 /*
 536  * x
 537  */
 538         w = x - (double)xIndex[1];
 539         xWeight[1] = 3.0 / 4.0 - w * w;
 540         xWeight[2] = (1.0 / 2.0) * (w - xWeight[1] + 1.0);
 541         xWeight[0] = 1.0 - xWeight[1] - xWeight[2];
 542 /*
 543  * y
 544  */
 545         w = y - (double)yIndex[1];
 546         yWeight[1] = 3.0 / 4.0 - w * w;
 547         yWeight[2] = (1.0 / 2.0) * (w - yWeight[1] + 1.0);
 548         yWeight[0] = 1.0 - yWeight[1] - yWeight[2];
 549         break;
 550
 551     case 3L:
 552 /*
 553  * x
 554  */
 555         w = x - (double)xIndex[1];
 556         xWeight[3] = (1.0 / 6.0) * w * w * w;
 557         xWeight[0] = (1.0 / 6.0) + (1.0 / 2.0) * w * (w - 1.0) - xWeight[3];
 558         xWeight[2] = w + xWeight[0] - 2.0 * xWeight[3];
 559         xWeight[1] = 1.0 - xWeight[0] - xWeight[2]
 560             - xWeight[3];
 561 /*
 562  * y
 563  */
 564         w = y - (double)yIndex[1];
 565         yWeight[3] = (1.0 / 6.0) * w * w * w;
 566         yWeight[0] = (1.0 / 6.0) + (1.0 / 2.0) * w * (w - 1.0) - yWeight[3];
 567         yWeight[2] = w + yWeight[0] - 2.0 * yWeight[3];
 568         yWeight[1] = 1.0 - yWeight[0] - yWeight[2]
 569             - yWeight[3];
 570         break;
 571
 572     case 4L:
 573 /*
 574  * x
 575  */
 576         w = x - (double)xIndex[2];
 577         w2 = w * w;
 578         t = (1.0 / 6.0) * w2;
 579         xWeight[0] = 1.0 / 2.0 - w;
 580         xWeight[0] *= xWeight[0];
 581         xWeight[0] *= (1.0 / 24.0) * xWeight[0];
 582         t0 = w * (t - 11.0 / 24.0);
 583         t1 = 19.0 / 96.0 + w2 * (1.0 / 4.0 - t);
 584         xWeight[1] = t1 + t0;
 585         xWeight[3] = t1 - t0;
 586         xWeight[4] = xWeight[0] + t0 + (1.0 / 2.0) * w;
 587         xWeight[2] = 1.0 - xWeight[0] - xWeight[1] - xWeight[3] - xWeight[4];
 588 /*
 589  * y
 590  */
 591         w = y - (double)yIndex[2];
 592         w2 = w * w;
 593         t = (1.0 / 6.0) * w2;
 594         yWeight[0] = 1.0 / 2.0 - w;
 595         yWeight[0] *= yWeight[0];
 596         yWeight[0] *= (1.0 / 24.0) * yWeight[0];
 597         t0 = w * (t - 11.0 / 24.0);
 598         t1 = 19.0 / 96.0 + w2 * (1.0 / 4.0 - t);
 599         yWeight[1] = t1 + t0;
 600         yWeight[3] = t1 - t0;
 601         yWeight[4] = yWeight[0] + t0 + (1.0 / 2.0) * w;
 602         yWeight[2] = 1.0 - yWeight[0] - yWeight[1] - yWeight[3] - yWeight[4];
 603         break;
 604
 605     case 5L:
 606 /*
 607  * x
 608  */
 609         w = x - (double)xIndex[2];
 610         w2 = w * w;
 611         xWeight[5] = (1.0 / 120.0) * w * w2 * w2;
 612         w2 -= w;
 613         w4 = w2 * w2;
 614         w -= 1.0 / 2.0;
 615         t = w2 * (w2 - 3.0);
 616         xWeight[0] = (1.0 / 24.0) * (1.0 / 5.0 + w2 + w4) - xWeight[5];
 617         t0 = (1.0 / 24.0) * (w2 * (w2 - 5.0) + 46.0 / 5.0);
 618         t1 = (-1.0 / 12.0) * w * (t + 4.0);
 619         xWeight[2] = t0 + t1;
 620         xWeight[3] = t0 - t1;
 621         t0 = (1.0 / 16.0) * (9.0 / 5.0 - t);
 622         t1 = (1.0 / 24.0) * w * (w4 - w2 - 5.0);
 623         xWeight[1] = t0 + t1;
 624         xWeight[4] = t0 - t1;
 625 /*
 626  * y
 627  */
 628         w = y - (double)yIndex[2];
 629         w2 = w * w;
 630         yWeight[5] = (1.0 / 120.0) * w * w2 * w2;
 631         w2 -= w;
 632         w4 = w2 * w2;
 633         w -= 1.0 / 2.0;
 634         t = w2 * (w2 - 3.0);
 635         yWeight[0] = (1.0 / 24.0) * (1.0 / 5.0 + w2 + w4) - yWeight[5];
 636         t0 = (1.0 / 24.0) * (w2 * (w2 - 5.0) + 46.0 / 5.0);
 637         t1 = (-1.0 / 12.0) * w * (t + 4.0);
 638         yWeight[2] = t0 + t1;
 639         yWeight[3] = t0 - t1;
 640         t0 = (1.0 / 16.0) * (9.0 / 5.0 - t);
 641         t1 = (1.0 / 24.0) * w * (w4 - w2 - 5.0);
 642         yWeight[1] = t0 + t1;
 643         yWeight[4] = t0 - t1;
 644         break;
 645
 646     default:
 647         printf("Invalid spline degree\n");
 648         return(0.0);
 649     }
 650
 651 /*
 652  * apply the mirror boundary conditions
 653  */
 654     for (k = 0L; k <= SplineDegree; k++) {
 655         xIndex[k] = (Width == 1L) ? (0L) : ((xIndex[k] < 0L) ?
 656                                             (-xIndex[k] - Width2
 657                                              * ((-xIndex[k]) / Width2))
 658                                             : (xIndex[k] - Width2
 659                                                * (xIndex[k] / Width2)));
 660         if (Width <= xIndex[k]) {
 661             xIndex[k] = Width2 - xIndex[k];
 662         }
 663         yIndex[k] = (Height == 1L) ? (0L) : ((yIndex[k] < 0L) ?
 664                                              (-yIndex[k] - Height2
 665                                               * ((-yIndex[k]) / Height2))
 666                                              : (yIndex[k] - Height2
 667                                                 * (yIndex[k] / Height2)));
 668         if (Height <= yIndex[k]) {
 669             yIndex[k] = Height2 - yIndex[k];
 670         }
 671     }
 672
 673 /*
 674  * perform interpolation
 675  */
 676     interpolated = 0.0;
 677     for (j = 0L; j <= SplineDegree; j++) {
 678         p = Bcoeff + (ptrdiff_t)(yIndex[j] * Width);
 679         w = 0.0;
 680         for (i = 0L; i <= SplineDegree; i++) {
 681             w += xWeight[i] * p[xIndex[i]];
 682         }
 683         interpolated += yWeight[j] * w;
 684     }
 685     return(interpolated);
 686 }
 687
 688
 689
 690 static gchar *
 691 imgproc_deform_frame_from_pivdata(guint16 **img,
 692                                   const GpivImagePar *image_par,
 693                                   const GpivPivData *piv_data,
 694                                   const gdouble magnitude)
 695 /*-----------------------------------------------------------------------------
 696  */
 697 {
 698     gchar *err_msg = NULL;
 699     gint i, j;
 700     GpivPivData *pd = NULL;
 701     gint k, l, count = 0, window = 4;
 702     gfloat dy_sum, dx_sum;
 703
 704
 705     if (img[0] == NULL) {
 706         err_msg = "imgproc_deform_frame_from_pivdata: image->frame1 == NULL";
 707         return (err_msg);
 708     }
 709
 710     /*
 711      * Interpolate piv data at each pixel point of the images.
 712      */
 713
 714     pd = gpiv_alloc_pivdata (image_par->ncolumns,
 715                              image_par->nrows);
 716 /* OK: results approved */
 717 #pragma omp parallel for shared(pd) private(i, j)
 718     for (i = 0; i < pd->ny; i++) {
 719         for (j = 0; j < pd->nx; j++) {
 720             pd->point_x[i][j] = (float) j;
 721             pd->point_y[i][j] = (float) i;
 722         }
 723     }
 724     gpiv_piv_dxdy_at_new_grid (piv_data, pd);
 725
 726 #undef SMOOTH
 727 #ifdef SMOOTH
 728     /*
 729      * Smoothing
 730      */
 731     for (i = 0; i < pd->ny; i++) {
 732         for (j = 0; j < pd->nx; j++) {
 733             for (k = -window / 2; k < window / 2; k++) {
 734                 if (i + k > 0 && i + k < pd->ny) {
 735                     for (l = -window / 2; l < window / 2; l++) {
 736                         if (j + l > 0 && j + l < pd->nx) {
 737                             dy_sum += pd->dy[i+k][j+l];
 738                             dx_sum += pd->dx[i+k][j+l];
 739                             count++;
 740                         }
 741                     }
 742                 }
 743             }
 744             count = 0;
 745         }
 746     }
 747 #endif /* SMOOTH */
 748
 749 /* OK: results approved */
 750 #pragma omp parallel for shared(pd) private(i, j)
 751     for (i = 0; i < pd->ny; i++) {
 752         for (j = 0; j < pd->nx; j++) {
 753             pd->dx[i][j] = magnitude * pd->dx[i][j];
 754             pd->dy[i][j] = magnitude * pd->dy[i][j];
 755         }
 756     }
 757
 758     imgproc_deform_frame (img, image_par, pd);
 759
 760
 761     gpiv_free_pivdata (pd);
 762 #ifdef SMOOTH
 763 #undefine SMOOTH
 764 #endif
 765     return (err_msg);
 766 }
 767
 768
 769 static int
 770 imgproc_deform_frame (guint16 **img,
 771                       const GpivImagePar *image_par,
 772                       const GpivPivData *piv_data
 773                       )
 774 /*-----------------------------------------------------------------------------
 775  * DESCRIPTION:
 776  *      Image deformation/interpolation routine for a single image frame
 777  *      Origins from demo.c
 778  *
 779  * INPUTS:
 780  *      piv_data:       PIV data structure which defines shift/deformation at
 781  *                      each image pixel
 782  *      img:            input image frame
 783  *
 784  * OUTPUTS:
 785  *      img:            deformed image frame
 786  *
 787  * RETURNS:
 788  *
 789  *---------------------------------------------------------------------------*/
 790 {
 791     float       *ImageRasterArray = NULL, *OutputImage = NULL;
 792     float       *p;
 793     double      x1, y1;
 794     long        Width = image_par->ncolumns, Height = image_par->nrows;
 795     long        SplineDegree = 5;
 796     long        x, y;
 797     int         Masking = 1;
 798     int         Error;
 799
 800     guint img_top = (1 << image_par->depth) - 1;
 801     unsigned char *Line;
 802     gint i, j;
 803
 804 /*     double   a11, a12, a21, a22; */
 805 /*     double   x0, y0, x1, y1; */
 806 /*     double   xOrigin = 0.0, yOrigin = 0.0; */
 807 /*     double   Angle = 1.0, xShift = 0.0, yShift = 0.0; */
 808
 809     long        rounded;
 810
 811
 812     ImageRasterArray = gpiv_vector (image_par->nrows * image_par->ncolumns);
 813     p = ImageRasterArray;
 814
 815 /* OK: results approved */
 816 #pragma omp parallel for shared(image_par, img, p) private(i, j)
 817     for (i = 0; i < image_par->nrows; i++) {
 818         for (j = 0; j < image_par->ncolumns; j++) {
 819             p[i*image_par->ncolumns + j] = (float) img[i][j];
 820 /*             *p++ = (float) img[i][j]; */
 821         }
 822     }
 823
 824
 825 /*
 826  * convert between a representation based on image samples
 827  * and a representation based on image B-spline coefficients
 828  */
 829     if (Error = SamplesToCoefficients (ImageRasterArray, Width, Height,
 830                                        SplineDegree)) {
 831         gpiv_free_vector (ImageRasterArray);
 832         g_warning("LIBGPIV internal error: imgproc_deform_frame: Change of basis failed");
 833         return(1);
 834     }
 835
 836 /*
 837  * visit all pixels of the output image and assign their value
 838  */
 839     OutputImage = gpiv_vector (image_par->nrows * image_par->ncolumns);
 840     p = OutputImage;
 841
 842 /* BUGFIX imgproc_deform_frame: OMP causes differences in results */
 843 /* BUGFIX imgproc_deform_frame: When disabled OMP in imgproc_deform */
 844 /* KEEP OUT! #pragma omp parallel for private(x1, y1) */
 845     for (i = 0; i < image_par->nrows; i++) {
 846         for (j = 0; j < image_par->ncolumns; j++) {
 847             x1 = (double) j - (double) piv_data->dx[i][j];
 848             y1 = (double) i - (double) piv_data->dy[i][j];
 849             if (Masking) {
 850                 if ((x1 <= -0.5)
 851                     || (((double)Width - 0.5) <= x1)
 852                     || (y1 <= -0.5)
 853                     || (((double)Height - 0.5) <= y1)) {
 854 /*                     *p++ = 0.0F; */
 855                     p[i*image_par->ncolumns + j] = 0.0F;
 856                 } else {
 857 /*                     *p++ = */
 858                     p[i*image_par->ncolumns + j] =
 859                         (float)InterpolatedValue (ImageRasterArray,
 860                                                   Width, Height,
 861                                                   x1, y1, SplineDegree);
 862                 }
 863             } else {
 864 /*                 *p++ = */
 865                 p[i*image_par->ncolumns + j] =
 866                     (float)InterpolatedValue (ImageRasterArray,
 867                                               Width, Height,
 868                                               x1, y1, SplineDegree);
 869             }
 870
 871         }
 872     }
 873
 874 /*
 875  * Store back into img
 876  */
 877
 878     p = OutputImage;
 879 /* OK: results approved */
 880 #pragma omp parallel for private(i, j, rounded)
 881     for (i = 0; i < image_par->nrows; i++) {
 882         for (j = 0; j < image_par->ncolumns; j++) {
 883             rounded = (long)floor ((double) p[i*image_par->ncolumns + j] + 0.5);
 884             if (rounded < 0.0) rounded = 0.0;
 885             if (rounded > (gfloat) img_top /* 255.0 */) rounded = (gfloat) img_top /* 255.0 */;
 886             img[i][j] =  (guint16) rounded;
 887         }
 888     }
 889
 890
 891
 892     gpiv_free_vector (ImageRasterArray);
 893     gpiv_free_vector (OutputImage);
 894     return(0);
 895 }
 896
 897
 898