geomwarp.c Source 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 CSPLINE (
00014 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 ) { register double dx1, dx2, dx3, dx4, value;
00022 
00023     if (TIMES) TIMING(T_CSPLINE);
00024 
00025     /* compute interpolation value */
00026     dx1 = 1.0 + dx;
00027     dx2 =       dx;
00028     dx3 = 1.0 - dx;
00029     dx4 = 2.0 - dx;
00030     value = (-4.0 + 8.0*dx1 - 5.0*dx1*dx1 + dx1*dx1*dx1) * alpha  * g1  +
00031             (1.0 - (alpha+3.0)*dx2*dx2 + (alpha+2.0)*dx2*dx2*dx2) * g2  +
00032             (1.0 - (alpha+3.0)*dx3*dx3 + (alpha+2.0)*dx3*dx3*dx3) * g3  +
00033             (-4.0 + 8.0*dx4 - 5.0*dx4*dx4 + dx4*dx4*dx4) * alpha  * g4;
00034 
00035     the_end:
00036     if (TIMES) TIMING(T_EXIT);
00037     return(value);
00038 }
00039 
00040 /*-Copyright Information------------------------------------------------------*/
00041 /* Copyright (c) 1988 by the University of Arizona Digital Image Analysis Lab */
00042 /*----------------------------------------------------------------------------*/
00043 /*-General Information--------------------------------------------------------*/
00044 /*                                                                            */
00045 /*   This function geometrically warps a "transformation" image with          */
00046 /*   coordinate system (x',y') to a "reference" image with coordinate         */
00047 /*   system (x,y), using a power series polynomial distortion model.          */
00048 /*   The polynomial equations for the warp are:                               */
00049 /*                                                                            */
00050 /*      x  =    a1                      y  =    b1                            */
00051 /*            + a2 * x'                       + b2 * x'                       */
00052 /*            + a3 * y'                       + b3 * y'                       */
00053 /*            + a4 * x' * y'                  + b4 * x' * y'                  */
00054 /*            + a5 * x' * x'                  + b5 * x' * x'                  */
00055 /*            + a6 * y' * y'                  + b6 * y' * y'                  */
00056 /*                                                                            */
00057 /*----------------------------------------------------------------------------*/
00058 /*-Interface Information------------------------------------------------------*/
00059 void GEOMWARP (
00060 IMAGE  *in,     /*  I   Pointer to the input image.                           */
00061 short  nbnd,    /*  I   Number of bands in the output image.                  */
00062 short  nlin,    /*  I   Number of lines in the output image.                  */
00063 short  npix,    /*  I   Number of pixels/line in the output image.            */
00064 double *xc,     /*  I   Address of an array[nterms],                          */
00065                 /*      containing the x-transformation coefficients.         */
00066 double *yc,     /*  I   Address of an array[nterms],                          */
00067                 /*      containing the y-transformation coefficients.         */
00068 short  nterms,  /*  I   Number of terms in the power series polynomial.       */
00069 double option,  /*  I   Switch, determines interpolation method:              */
00070                 /*      < 0.0   -   Cubic interpolation (alpha = option).     */
00071                 /*      = 0.0   -   Bilinear interpolation.                   */
00072                 /*      > 0.0   -   Nearest neighbor interpolation.           */
00073 PIXEL  glev,    /*  I   Graylevel to fill otherwise empty pixels.             */
00074 IMAGE  **out    /*  O   Address of a pointer to the output image.             */
00075 /*----------------------------------------------------------------------------*/
00076 ) { register short i, j, k, ix, iy;
00077     char   msg[SLEN];
00078     double x, y, gl1, gl2, gl3, gl4, pinc, psum;
00079 
00080     if (TIMES) TIMING(T_GEOMWARP);
00081     if (NAMES) {
00082         MESSAGE('I',"");
00083         MESSAGE('I',"GEOMWARP");
00084         MESSAGE('I',"");
00085         sprintf(msg," Input image:                  %s",in->text);
00086         MESSAGE('I',msg);
00087         sprintf(msg," Output image size: bands:     %d",nbnd);
00088         MESSAGE('I',msg);
00089         sprintf(msg,"                    lines:     %d",nlin);
00090         MESSAGE('I',msg);
00091         sprintf(msg,"                    pixels:    %d",npix);
00092         MESSAGE('I',msg);
00093         sprintf(msg," Fill graylevel:             %12.4e",glev);
00094         MESSAGE('I',msg);
00095         sprintf(msg," Number of polynomial terms:   %d",nterms);
00096         MESSAGE('I',msg);
00097         sprintf(msg," Transformation Coefficients:  x  =   %12.4e            y  =   %12.4e",xc[0],yc[0]);
00098         MESSAGE('I',msg);
00099         sprintf(msg,"                                     +%12.4e * x'             +%12.4e * x'",xc[1],yc[1]);
00100         MESSAGE('I',msg);
00101         sprintf(msg,"                                     +%12.4e * y'             +%12.4e * y'",xc[2],yc[2]);
00102         MESSAGE('I',msg);
00103         if (nterms > 3) {
00104             sprintf(msg,"                                     +%12.4e * x' * y'        +%12.4e * x' * y'",xc[3],yc[3]);
00105             MESSAGE('I',msg);
00106         }
00107         if (nterms > 4) {
00108             sprintf(msg,"                                     +%12.4e * x' * x'        +%12.4e * x' * x'",xc[4],yc[4]);
00109             MESSAGE('I',msg);
00110         }
00111         if (nterms > 5) {
00112             sprintf(msg,"                                     +%12.4e * y' * y'        +%12.4e * y' * y'",xc[5],yc[5]);
00113             MESSAGE('I',msg);
00114         }
00115         if (option < 0.0) {
00116             sprintf(msg," Parametric cubic interpolation, alpha = %5.2f",option);
00117             MESSAGE('I',msg);
00118         } else if (option == 0.0) {
00119             MESSAGE('I'," Bilinear interpolation");
00120         } else if (option > 0.0) {
00121             MESSAGE('I'," Nearest-neighbor interpolation");
00122         }
00123         MESSAGE('I'," ...............");
00124     }
00125     
00126     /* check input */
00127     if (!CHECKIMG(in)) {
00128         MESSAGE('E'," Can't identify image.");
00129         goto the_end;
00130     } else if (MINTERMS > nterms) {
00131         sprintf(msg," Number of terms must be equal to or greater than %d.",MINTERMS);
00132         MESSAGE('E',msg);
00133         goto the_end;
00134     } else if (nterms > MAXTERMS) {
00135         sprintf(msg," Number of terms must be less than or equal to %d.",MAXTERMS);
00136         MESSAGE('E',msg);
00137         goto the_end;
00138     }
00139 
00140     /* create image of appropriate size */
00141     if (!CHECKIMG(*out)) GETMEM(nbnd,nlin,npix,out);
00142     if (!*out) goto the_end;
00143 
00144     /* initialize progress indicator */
00145     if (LINES  &&  !PROGRESS(psum=0.0)) goto the_end;
00146     pinc = 1.0/(double)nbnd/(double)nlin;
00147 
00148     /* cubic interpolation */
00149     if (option < 0.0) {
00150         for (i=0; i<nbnd; i++) {
00151             for (j=0; j<nlin; j++) {
00152                 for (k=0; k<npix; k++) {
00153                     x = y = 0.0;
00154                     switch (nterms) {
00155                     case 6:
00156                         x += xc[5]*(double)j*(double)j;
00157                         y += yc[5]*(double)j*(double)j;
00158                     case 5:
00159                         x += xc[4]*(double)k*(double)k;
00160                         y += yc[4]*(double)k*(double)k;
00161                     case 4:
00162                         x += xc[3]*(double)j*(double)k;
00163                         y += yc[3]*(double)j*(double)k;
00164                     case 3:
00165                         x += xc[0] + xc[1]*(double)k + xc[2]*(double)j;
00166                         y += yc[0] + yc[1]*(double)k + yc[2]*(double)j;
00167                     }
00168                     ix = (short)x;
00169                     iy = (short)y;
00170                     if (1 <= ix  &&  ix <= in->npix-3  &&  1 <= iy  &&  iy <= in->nlin-3) {
00171                         gl1 = CSPLINE(option,x-(double)ix,(double)in->data[i][iy-1][ix-1],(double)in->data[i][iy-1][ix],(double)in->data[i][iy-1][ix+1],(double)in->data[i][iy-1][ix+2]);
00172                         gl2 = CSPLINE(option,x-(double)ix,(double)in->data[i][iy  ][ix-1],(double)in->data[i][iy  ][ix],(double)in->data[i][iy  ][ix+1],(double)in->data[i][iy  ][ix+2]);
00173                         gl3 = CSPLINE(option,x-(double)ix,(double)in->data[i][iy+1][ix-1],(double)in->data[i][iy+1][ix],(double)in->data[i][iy+1][ix+1],(double)in->data[i][iy+1][ix+2]);
00174                         gl4 = CSPLINE(option,x-(double)ix,(double)in->data[i][iy+2][ix-1],(double)in->data[i][iy+2][ix],(double)in->data[i][iy+2][ix+1],(double)in->data[i][iy+2][ix+2]);
00175                         (*out)->data[i][j][k] = (PIXEL)CSPLINE(option,y-(double)iy,gl1,gl2,gl3,gl4);
00176                     } else if (ix == x  &&  1 <= ix  &&  ix <= in->npix-2  &&  1 <= iy  &&  iy <= in->nlin-3) {
00177                         (*out)->data[i][j][k] = (PIXEL)CSPLINE(option,y-(double)iy,(double)in->data[i][iy-1][ix],(double)in->data[i][iy][ix],(double)in->data[i][iy+1][ix],(double)in->data[i][iy+2][ix]);
00178                     } else if (1 <= ix  &&  ix <= in->npix-3  &&  iy == y  &&  1 <= iy  &&  iy <= in->nlin-2) {
00179                         (*out)->data[i][j][k] = (PIXEL)CSPLINE(option,x-(double)ix,(double)in->data[i][iy][ix-1],(double)in->data[i][iy][ix],(double)in->data[i][iy][ix+1],(double)in->data[i][iy][ix+2]);
00180                     } else if (ix == x  &&  1 <= ix  &&  ix <= in->npix-2  &&  iy == y  &&  1 <= iy  &&  iy <= in->nlin-2) {
00181                         (*out)->data[i][j][k] = in->data[i][iy][ix];
00182                     } else {
00183                         (*out)->data[i][j][k] = glev;
00184                     }
00185                 }
00186                 if (LINES  &&  !PROGRESS(psum+=pinc)) goto the_end;
00187             }
00188         }
00189 
00190     /* bilinear interpolation */
00191     } else if (option == 0.0) {
00192         for (i=0; i<nbnd; i++) {
00193             for (j=0; j<nlin; j++) {
00194                 for (k=0; k<npix; k++) {
00195                     x = y = 0.0;
00196                     switch (nterms) {
00197                     case 6:
00198                         x += xc[5]*(double)j*(double)j;
00199                         y += yc[5]*(double)j*(double)j;
00200                     case 5:
00201                         x += xc[4]*(double)k*(double)k;
00202                         y += yc[4]*(double)k*(double)k;
00203                     case 4:
00204                         x += xc[3]*(double)j*(double)k;
00205                         y += yc[3]*(double)j*(double)k;
00206                     case 3:
00207                         x += xc[0] + xc[1]*(double)k + xc[2]*(double)j;
00208                         y += yc[0] + yc[1]*(double)k + yc[2]*(double)j;
00209                     }
00210                     ix = (short)x;
00211                     iy = (short)y;
00212                     if (0 <= ix  &&  ix <= in->npix-2  &&  0 <= iy  &&  iy <= in->nlin-2) {
00213                         gl1 = (double)in->data[i][iy  ][ix] + (x-(double)ix)*(double)(in->data[i][iy  ][ix+1]-in->data[i][iy  ][ix]);
00214                         gl2 = (double)in->data[i][iy+1][ix] + (x-(double)ix)*(double)(in->data[i][iy+1][ix+1]-in->data[i][iy+1][ix]);
00215                         (*out)->data[i][j][k] = (PIXEL)(gl1 + (y-(double)iy)*(gl2-gl1));
00216                     } else if (ix == x  &&  0 <= ix  &&  ix <= in->npix-1  &&  0 <= iy  &&  iy <= in->nlin-2) {
00217                         (*out)->data[i][j][k] = (PIXEL)((double)in->data[i][iy][ix] + (y-(double)iy)*(double)(in->data[i][iy+1][ix]-in->data[i][iy][ix]));
00218                     } else if (0 <= ix  &&  ix <= in->npix-2  &&  iy == y  &&  0 <= iy  &&  iy <= in->nlin-1) {
00219                         (*out)->data[i][j][k] = (PIXEL)((double)in->data[i][iy][ix] + (x-(double)ix)*(double)(in->data[i][iy][ix+1]-in->data[i][iy][ix]));
00220                     } else if (ix == x  &&  0 <= ix  &&  ix <= in->npix-1  &&  iy == y  &&  0 <= iy  &&  iy <= in->nlin-1) {
00221                         (*out)->data[i][j][k] = in->data[i][iy][ix];
00222                     } else {
00223                         (*out)->data[i][j][k] = glev;
00224                     }
00225                 }
00226                 if (LINES  &&  !PROGRESS(psum+=pinc)) goto the_end;
00227             }
00228         }
00229 
00230     /* nearest neighbor interpolation */
00231     } else /* (option > 0.0) */ {
00232         for (i=0; i<nbnd; i++) {
00233             for (j=0; j<nlin; j++) {
00234                 for (k=0; k<npix; k++) {
00235                     x = y = 0.0;
00236                     switch (nterms) {
00237                     case 6:
00238                         x += xc[5]*(double)j*(double)j;
00239                         y += yc[5]*(double)j*(double)j;
00240                     case 5:
00241                         x += xc[4]*(double)k*(double)k;
00242                         y += yc[4]*(double)k*(double)k;
00243                     case 4:
00244                         x += xc[3]*(double)j*(double)k;
00245                         y += yc[3]*(double)j*(double)k;
00246                     case 3:
00247                         x += xc[0] + xc[1]*(double)k + xc[2]*(double)j;
00248                         y += yc[0] + yc[1]*(double)k + yc[2]*(double)j;
00249                     }
00250                     ix = rnd(x);
00251                     iy = rnd(y);
00252                     if (0 <= ix  &&  ix <= in->npix-1  &&  0 <= iy  &&  iy <= in->nlin-1) {
00253                         (*out)->data[i][j][k] = in->data[i][iy][ix];
00254                     } else {
00255                         (*out)->data[i][j][k] = glev;
00256                     }
00257                 }
00258                 if (LINES  &&  !PROGRESS(psum+=pinc)) goto the_end;
00259             }
00260         }
00261     }
00262 
00263     the_end:
00264     if (TIMES) TIMING(T_EXIT);
00265 }