Main Page   Data Structures   File List   Data Fields   Globals  

geomwarp.c

Go to the documentation of this file.
00001 #include        "sadie.h"
00002 
00003 /*-Copyright Information------------------------------------------------------*/
00004 /* Copyright (c) 1990 by the University of Arizona Digital Image Analysis Lab */
00005 /*----------------------------------------------------------------------------*/
00006 /*-General Information--------------------------------------------------------*/
00007 /*                                                                            */
00008 /*   This low-level function computes a cubic spline interpolation between    */
00009 /*   four datapoints.                                                         */
00010 /*                                                                            */
00011 /*----------------------------------------------------------------------------*/
00012 /*-Interface Information------------------------------------------------------*/
00013 static double
00014 CSPLINE (register double alpha, /*  I   Interpolation function parameter (alpha < 0). */
00015          double dx,             /*  I   Fractional location of interpolated datapoint.        */
00016          double g1,             /*  I   First datapoint.                                      */
00017          double g2,             /*  I   Second datapoint.                                     */
00018          double g3,             /*  I   Third datapoint.                                      */
00019          double g4              /*  I   Fourth datapoint.                                     */
00020 /*----------------------------------------------------------------------------*/
00021   )
00022 {
00023   register double dx1, dx2, dx3, dx4, value;
00024 
00025   if (TIMES)
00026     TIMING (T_CSPLINE);
00027 
00028   /* compute interpolation value */
00029   dx1 = 1.0 + dx;
00030   dx2 = dx;
00031   dx3 = 1.0 - dx;
00032   dx4 = 2.0 - dx;
00033   value =
00034     (-4.0 + 8.0 * dx1 - 5.0 * dx1 * dx1 + dx1 * dx1 * dx1) * alpha * g1 +
00035     (1.0 - (alpha + 3.0) * dx2 * dx2 + (alpha + 2.0) * dx2 * dx2 * dx2) * g2 +
00036     (1.0 - (alpha + 3.0) * dx3 * dx3 + (alpha + 2.0) * dx3 * dx3 * dx3) * g3 +
00037     (-4.0 + 8.0 * dx4 - 5.0 * dx4 * dx4 + dx4 * dx4 * dx4) * alpha * g4;
00038 
00039 the_end:
00040   if (TIMES)
00041     TIMING (T_EXIT);
00042   return (value);
00043 }
00044 
00045 /*-Copyright Information------------------------------------------------------*/
00046 /* Copyright (c) 1988 by the University of Arizona Digital Image Analysis Lab */
00047 /*----------------------------------------------------------------------------*/
00048 /*-General Information--------------------------------------------------------*/
00049 /*                                                                            */
00050 /*   This function geometrically warps a "transformation" image with          */
00051 /*   coordinate system (x',y') to a "reference" image with coordinate         */
00052 /*   system (x,y), using a power series polynomial distortion model.          */
00053 /*   The polynomial equations for the warp are:                               */
00054 /*                                                                            */
00055 /*      x  =    a1                      y  =    b1                            */
00056 /*            + a2 * x'                       + b2 * x'                       */
00057 /*            + a3 * y'                       + b3 * y'                       */
00058 /*            + a4 * x' * y'                  + b4 * x' * y'                  */
00059 /*            + a5 * x' * x'                  + b5 * x' * x'                  */
00060 /*            + a6 * y' * y'                  + b6 * y' * y'                  */
00061 /*                                                                            */
00062 /*----------------------------------------------------------------------------*/
00063 /*-Interface Information------------------------------------------------------*/
00064 void
00065 GEOMWARP (IMAGE * in,           /*  I   Pointer to the input image.                           */
00066           short nbnd,           /*  I   Number of bands in the output image.                  */
00067           short nlin,           /*  I   Number of lines in the output image.                  */
00068           short npix,           /*  I   Number of pixels/line in the output image.            */
00069           double *xc,           /*  I   Address of an array[nterms],                          */
00070           /*      containing the x-transformation coefficients.         */
00071           double *yc,           /*  I   Address of an array[nterms],                          */
00072           /*      containing the y-transformation coefficients.         */
00073           short nterms,         /*  I   Number of terms in the power series polynomial.       */
00074           double option,        /*  I   Switch, determines interpolation method:              */
00075           /*      < 0.0   -   Cubic interpolation (alpha = option).     */
00076           /*      = 0.0   -   Bilinear interpolation.                   */
00077           /*      > 0.0   -   Nearest neighbor interpolation.           */
00078           PIXEL glev,           /*  I   Graylevel to fill otherwise empty pixels.             */
00079           IMAGE ** out          /*  O   Address of a pointer to the output image.             */
00080 /*----------------------------------------------------------------------------*/
00081   )
00082 {
00083   register short i, j, k, ix, iy;
00084   char msg[SLEN];
00085   double x, y, gl1, gl2, gl3, gl4, pinc, psum;
00086 
00087   if (TIMES)
00088     TIMING (T_GEOMWARP);
00089   if (NAMES)
00090     {
00091       MESSAGE ('I', "");
00092       MESSAGE ('I', "GEOMWARP");
00093       MESSAGE ('I', "");
00094       sprintf (msg, " Input image:                  %s", in->text);
00095       MESSAGE ('I', msg);
00096       sprintf (msg, " Output image size: bands:     %d", nbnd);
00097       MESSAGE ('I', msg);
00098       sprintf (msg, "                    lines:     %d", nlin);
00099       MESSAGE ('I', msg);
00100       sprintf (msg, "                    pixels:    %d", npix);
00101       MESSAGE ('I', msg);
00102       sprintf (msg, " Fill graylevel:             %12.4e", glev);
00103       MESSAGE ('I', msg);
00104       sprintf (msg, " Number of polynomial terms:   %d", nterms);
00105       MESSAGE ('I', msg);
00106       sprintf (msg,
00107                " Transformation Coefficients:  x  =   %12.4e            y  =   %12.4e",
00108                xc[0], yc[0]);
00109       MESSAGE ('I', msg);
00110       sprintf (msg,
00111                "                                     +%12.4e * x'             +%12.4e * x'",
00112                xc[1], yc[1]);
00113       MESSAGE ('I', msg);
00114       sprintf (msg,
00115                "                                     +%12.4e * y'             +%12.4e * y'",
00116                xc[2], yc[2]);
00117       MESSAGE ('I', msg);
00118       if (nterms > 3)
00119         {
00120           sprintf (msg,
00121                    "                                     +%12.4e * x' * y'        +%12.4e * x' * y'",
00122                    xc[3], yc[3]);
00123           MESSAGE ('I', msg);
00124         }
00125       if (nterms > 4)
00126         {
00127           sprintf (msg,
00128                    "                                     +%12.4e * x' * x'        +%12.4e * x' * x'",
00129                    xc[4], yc[4]);
00130           MESSAGE ('I', msg);
00131         }
00132       if (nterms > 5)
00133         {
00134           sprintf (msg,
00135                    "                                     +%12.4e * y' * y'        +%12.4e * y' * y'",
00136                    xc[5], yc[5]);
00137           MESSAGE ('I', msg);
00138         }
00139       if (option < 0.0)
00140         {
00141           sprintf (msg, " Parametric cubic interpolation, alpha = %5.2f",
00142                    option);
00143           MESSAGE ('I', msg);
00144         }
00145       else if (option == 0.0)
00146         {
00147           MESSAGE ('I', " Bilinear interpolation");
00148         }
00149       else if (option > 0.0)
00150         {
00151           MESSAGE ('I', " Nearest-neighbor interpolation");
00152         }
00153       MESSAGE ('I', " ...............");
00154     }
00155 
00156   /* check input */
00157   if (!CHECKIMG (in))
00158     {
00159       MESSAGE ('E', " Can't identify image.");
00160       goto the_end;
00161     }
00162   else if (MINTERMS > nterms)
00163     {
00164       sprintf (msg, " Number of terms must be equal to or greater than %d.",
00165                MINTERMS);
00166       MESSAGE ('E', msg);
00167       goto the_end;
00168     }
00169   else if (nterms > MAXTERMS)
00170     {
00171       sprintf (msg, " Number of terms must be less than or equal to %d.",
00172                MAXTERMS);
00173       MESSAGE ('E', msg);
00174       goto the_end;
00175     }
00176 
00177   /* create image of appropriate size */
00178   if (!CHECKIMG (*out))
00179     GETMEM (nbnd, nlin, npix, out);
00180   if (!*out)
00181     goto the_end;
00182 
00183   /* initialize progress indicator */
00184   if (LINES && !PROGRESS (psum = 0.0))
00185     goto the_end;
00186   pinc = 1.0 / (double) nbnd / (double) nlin;
00187 
00188   /* cubic interpolation */
00189   if (option < 0.0)
00190     {
00191       for (i = 0; i < nbnd; i++)
00192         {
00193           for (j = 0; j < nlin; j++)
00194             {
00195               for (k = 0; k < npix; k++)
00196                 {
00197                   x = y = 0.0;
00198                   switch (nterms)
00199                     {
00200                     case 6:
00201                       x += xc[5] * (double) j *(double) j;
00202                       y += yc[5] * (double) j *(double) j;
00203                     case 5:
00204                       x += xc[4] * (double) k *(double) k;
00205                       y += yc[4] * (double) k *(double) k;
00206                     case 4:
00207                       x += xc[3] * (double) j *(double) k;
00208                       y += yc[3] * (double) j *(double) k;
00209                     case 3:
00210                       x += xc[0] + xc[1] * (double) k + xc[2] * (double) j;
00211                       y += yc[0] + yc[1] * (double) k + yc[2] * (double) j;
00212                     }
00213                   ix = (short) x;
00214                   iy = (short) y;
00215                   if (1 <= ix && ix <= in->npix - 3 && 1 <= iy
00216                       && iy <= in->nlin - 3)
00217                     {
00218                       gl1 =
00219                         CSPLINE (option, x - (double) ix,
00220                                  (double) in->data[i][iy - 1][ix - 1],
00221                                  (double) in->data[i][iy - 1][ix],
00222                                  (double) in->data[i][iy - 1][ix + 1],
00223                                  (double) in->data[i][iy - 1][ix + 2]);
00224                       gl2 =
00225                         CSPLINE (option, x - (double) ix,
00226                                  (double) in->data[i][iy][ix - 1],
00227                                  (double) in->data[i][iy][ix],
00228                                  (double) in->data[i][iy][ix + 1],
00229                                  (double) in->data[i][iy][ix + 2]);
00230                       gl3 =
00231                         CSPLINE (option, x - (double) ix,
00232                                  (double) in->data[i][iy + 1][ix - 1],
00233                                  (double) in->data[i][iy + 1][ix],
00234                                  (double) in->data[i][iy + 1][ix + 1],
00235                                  (double) in->data[i][iy + 1][ix + 2]);
00236                       gl4 =
00237                         CSPLINE (option, x - (double) ix,
00238                                  (double) in->data[i][iy + 2][ix - 1],
00239                                  (double) in->data[i][iy + 2][ix],
00240                                  (double) in->data[i][iy + 2][ix + 1],
00241                                  (double) in->data[i][iy + 2][ix + 2]);
00242                       (*out)->data[i][j][k] =
00243                         (PIXEL) CSPLINE (option, y - (double) iy, gl1, gl2,
00244                                          gl3, gl4);
00245                     }
00246                   else if (ix == x && 1 <= ix && ix <= in->npix - 2 && 1 <= iy
00247                            && iy <= in->nlin - 3)
00248                     {
00249                       (*out)->data[i][j][k] =
00250                         (PIXEL) CSPLINE (option, y - (double) iy,
00251                                          (double) in->data[i][iy - 1][ix],
00252                                          (double) in->data[i][iy][ix],
00253                                          (double) in->data[i][iy + 1][ix],
00254                                          (double) in->data[i][iy + 2][ix]);
00255                     }
00256                   else if (1 <= ix && ix <= in->npix - 3 && iy == y && 1 <= iy
00257                            && iy <= in->nlin - 2)
00258                     {
00259                       (*out)->data[i][j][k] =
00260                         (PIXEL) CSPLINE (option, x - (double) ix,
00261                                          (double) in->data[i][iy][ix - 1],
00262                                          (double) in->data[i][iy][ix],
00263                                          (double) in->data[i][iy][ix + 1],
00264                                          (double) in->data[i][iy][ix + 2]);
00265                     }
00266                   else if (ix == x && 1 <= ix && ix <= in->npix - 2 && iy == y
00267                            && 1 <= iy && iy <= in->nlin - 2)
00268                     {
00269                       (*out)->data[i][j][k] = in->data[i][iy][ix];
00270                     }
00271                   else
00272                     {
00273                       (*out)->data[i][j][k] = glev;
00274                     }
00275                 }
00276               if (LINES && !PROGRESS (psum += pinc))
00277                 goto the_end;
00278             }
00279         }
00280 
00281       /* bilinear interpolation */
00282     }
00283   else if (option == 0.0)
00284     {
00285       for (i = 0; i < nbnd; i++)
00286         {
00287           for (j = 0; j < nlin; j++)
00288             {
00289               for (k = 0; k < npix; k++)
00290                 {
00291                   x = y = 0.0;
00292                   switch (nterms)
00293                     {
00294                     case 6:
00295                       x += xc[5] * (double) j *(double) j;
00296                       y += yc[5] * (double) j *(double) j;
00297                     case 5:
00298                       x += xc[4] * (double) k *(double) k;
00299                       y += yc[4] * (double) k *(double) k;
00300                     case 4:
00301                       x += xc[3] * (double) j *(double) k;
00302                       y += yc[3] * (double) j *(double) k;
00303                     case 3:
00304                       x += xc[0] + xc[1] * (double) k + xc[2] * (double) j;
00305                       y += yc[0] + yc[1] * (double) k + yc[2] * (double) j;
00306                     }
00307                   ix = (short) x;
00308                   iy = (short) y;
00309                   if (0 <= ix && ix <= in->npix - 2 && 0 <= iy
00310                       && iy <= in->nlin - 2)
00311                     {
00312                       gl1 =
00313                         (double) in->data[i][iy][ix] + (x -
00314                                                         (double) ix) *
00315                         (double) (in->data[i][iy][ix + 1] -
00316                                   in->data[i][iy][ix]);
00317                       gl2 =
00318                         (double) in->data[i][iy + 1][ix] + (x -
00319                                                             (double) ix) *
00320                         (double) (in->data[i][iy + 1][ix + 1] -
00321                                   in->data[i][iy + 1][ix]);
00322                       (*out)->data[i][j][k] =
00323                         (PIXEL) (gl1 + (y - (double) iy) * (gl2 - gl1));
00324                     }
00325                   else if (ix == x && 0 <= ix && ix <= in->npix - 1 && 0 <= iy
00326                            && iy <= in->nlin - 2)
00327                     {
00328                       (*out)->data[i][j][k] =
00329                         (PIXEL) ((double) in->data[i][iy][ix] +
00330                                  (y -
00331                                   (double) iy) * (double) (in->data[i][iy +
00332                                                                        1][ix]
00333                                                            -
00334                                                            in->
00335                                                            data[i][iy][ix]));
00336                     }
00337                   else if (0 <= ix && ix <= in->npix - 2 && iy == y && 0 <= iy
00338                            && iy <= in->nlin - 1)
00339                     {
00340                       (*out)->data[i][j][k] =
00341                         (PIXEL) ((double) in->data[i][iy][ix] +
00342                                  (x -
00343                                   (double) ix) *
00344                                  (double) (in->data[i][iy][ix + 1] -
00345                                            in->data[i][iy][ix]));
00346                     }
00347                   else if (ix == x && 0 <= ix && ix <= in->npix - 1 && iy == y
00348                            && 0 <= iy && iy <= in->nlin - 1)
00349                     {
00350                       (*out)->data[i][j][k] = in->data[i][iy][ix];
00351                     }
00352                   else
00353                     {
00354                       (*out)->data[i][j][k] = glev;
00355                     }
00356                 }
00357               if (LINES && !PROGRESS (psum += pinc))
00358                 goto the_end;
00359             }
00360         }
00361 
00362       /* nearest neighbor interpolation */
00363     }
00364   else                          /* (option > 0.0) */
00365     {
00366       for (i = 0; i < nbnd; i++)
00367         {
00368           for (j = 0; j < nlin; j++)
00369             {
00370               for (k = 0; k < npix; k++)
00371                 {
00372                   x = y = 0.0;
00373                   switch (nterms)
00374                     {
00375                     case 6:
00376                       x += xc[5] * (double) j *(double) j;
00377                       y += yc[5] * (double) j *(double) j;
00378                     case 5:
00379                       x += xc[4] * (double) k *(double) k;
00380                       y += yc[4] * (double) k *(double) k;
00381                     case 4:
00382                       x += xc[3] * (double) j *(double) k;
00383                       y += yc[3] * (double) j *(double) k;
00384                     case 3:
00385                       x += xc[0] + xc[1] * (double) k + xc[2] * (double) j;
00386                       y += yc[0] + yc[1] * (double) k + yc[2] * (double) j;
00387                     }
00388                   ix = rnd (x);
00389                   iy = rnd (y);
00390                   if (0 <= ix && ix <= in->npix - 1 && 0 <= iy
00391                       && iy <= in->nlin - 1)
00392                     {
00393                       (*out)->data[i][j][k] = in->data[i][iy][ix];
00394                     }
00395                   else
00396                     {
00397                       (*out)->data[i][j][k] = glev;
00398                     }
00399                 }
00400               if (LINES && !PROGRESS (psum += pinc))
00401                 goto the_end;
00402             }
00403         }
00404     }
00405 
00406 the_end:
00407   if (TIMES)
00408     TIMING (T_EXIT);
00409 }
Generated on Sun May 18 15:36:10 2003 for tclSadie by doxygen1.3