IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 18162


Ignore:
Timestamp:
Jun 17, 2008, 11:24:58 AM (18 years ago)
Author:
Paul Price
Message:

Adding function psImageInterpolateVarianceFactor. This involved splitting up parts of the interpolation into separate functions, so the code could be reused. tap_psImageInterpolate2 tests all pass.

Location:
trunk/psLib/src/imageops
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/psLib/src/imageops/psImageInterpolate.c

    r18043 r18162  
    77 *  @author Paul Price, IfA
    88 *
    9  *  @version $Revision: 1.11 $ $Name: not supported by cvs2svn $
    10  *  @date $Date: 2008-06-10 02:42:41 $
     9 *  @version $Revision: 1.12 $ $Name: not supported by cvs2svn $
     10 *  @date $Date: 2008-06-17 21:24:58 $
    1111 *
    1212 *  Copyright 2004-2007 Institute for Astronomy, University of Hawaii
    … …  
    132132}
    133133
    134 // Interpolation engine using interpolation kernel
    135 static psImageInterpolateStatus interpolateKernel(double *imageValue, double *varianceValue,
    136                                                   psMaskType *maskValue, float x, float y,
    137                                                   const psImageInterpolateOptions *options)
    138 {
    139     // Parameters have been checked by psImageInterpolate()
    140 
    141     int xNum, yNum;                     // Number of interpolation kernel pixels
    142     int xCentral, yCentral;             // Central pixel of the convolution
    143     switch (options->mode) {
     134// Setup for interpolation by kernel
     135static inline void interpolateKernelSetup(int *xNum, int *yNum, // Size of interpolation kernel, returned
     136                                          int *xCentral, int *yCentral, // Central pixel of convolution
     137                                          float x, float y, // Coordinates of interest
     138                                          psImageInterpolateMode mode // Mode for interpolation
     139                                          )
     140{
     141    switch (mode) {
    144142      case PS_INTERPOLATE_BILINEAR:
    145         xNum = yNum = 2;
     143        *xNum = *yNum = 2;
    146144        // Central pixel is the pixel below the point of interest
    147         xCentral = floor(x - 0.5 + FLT_EPSILON);
    148         yCentral = floor(y - 0.5 + FLT_EPSILON);
     145        *xCentral = floor(x - 0.5 + FLT_EPSILON);
     146        *yCentral = floor(y - 0.5 + FLT_EPSILON);
    149147        break;
    150148      case PS_INTERPOLATE_BICUBE:
    151149      case PS_INTERPOLATE_GAUSS:
    152         xNum = yNum = 3;
     150        *xNum = *yNum = 3;
    153151        // Central pixel is the closest pixel to the point of interest
    154         xCentral = x;
    155         yCentral = y;
     152        *xCentral = x;
     153        *yCentral = y;
    156154        break;
    157155      case PS_INTERPOLATE_FLAT:
    … …  
    162160        psAbort("Invalid interpolation mode.");
    163161    }
    164 
    165     const psImage *image = options->image; // Image of interest
    166     int xLast = image->numCols - 1;     // Last pixel in x
    167     int yLast = image->numRows - 1;     // Last pixel in y
    168 
    169     if (xCentral - (xNum - 1) / 2 < 0 || xCentral + xNum / 2 > xLast ||
    170         yCentral - (yNum - 1) / 2 < 0 || yCentral + yNum / 2 > yLast) {
    171         // At least one pixel of the interpolation kernel is off the image
    172         if (imageValue) {
    173             *imageValue = options->badImage;
    174         }
    175         if (varianceValue) {
    176             *varianceValue = options->badVariance;
    177         }
    178         if (maskValue) {
    179             *maskValue = options->badMask;
    180         }
    181         return PS_INTERPOLATE_STATUS_OFF;
    182     }
    183     double kernel[yNum][xNum];          // Interpolation kernel for straight interpolation
    184     switch (options->mode) {
     162}
     163
     164// Generate the interpolation kernel; it should be normalised to unity
     165static inline void interpolateKernelGenerate(int xNum, int yNum, // Size of interpolation kernel
     166                                             double kernel[xNum][yNum], // Kernel, to be set
     167                                             int xCentral, int yCentral, // Central pixel of convolution
     168                                             float x, float y, // Coordinates of interest
     169                                             psImageInterpolateMode mode // Mode for interpolation
     170                                             )
     171{
     172    switch (mode) {
    185173      case PS_INTERPOLATE_BILINEAR: {
    186174          double xFrac = x - 0.5 - xCentral; // Fraction of pixel in x
    … …  
    238226        psAbort("Invalid interpolation mode.");
    239227    }
     228}
     229
     230// Interpolation engine using interpolation kernel
     231static psImageInterpolateStatus interpolateKernel(double *imageValue, double *varianceValue,
     232                                                  psMaskType *maskValue, float x, float y,
     233                                                  const psImageInterpolateOptions *options)
     234{
     235    // Parameters have been checked by psImageInterpolate()
     236
     237    int xNum, yNum;                     // Number of interpolation kernel pixels
     238    int xCentral, yCentral;             // Central pixel of the convolution
     239    interpolateKernelSetup(&xNum, &yNum, &xCentral, &yCentral, x, y, options->mode);
     240
     241    const psImage *image = options->image; // Image of interest
     242    int xLast = image->numCols - 1;     // Last pixel in x
     243    int yLast = image->numRows - 1;     // Last pixel in y
     244
     245    if (xCentral - (xNum - 1) / 2 < 0 || xCentral + xNum / 2 > xLast ||
     246        yCentral - (yNum - 1) / 2 < 0 || yCentral + yNum / 2 > yLast) {
     247        // At least one pixel of the interpolation kernel is off the image
     248        if (imageValue) {
     249            *imageValue = options->badImage;
     250        }
     251        if (varianceValue) {
     252            *varianceValue = options->badVariance;
     253        }
     254        if (maskValue) {
     255            *maskValue = options->badMask;
     256        }
     257        return PS_INTERPOLATE_STATUS_OFF;
     258    }
     259
     260    double kernel[yNum][xNum];          // Interpolation kernel for straight interpolation
     261    interpolateKernelGenerate(xNum, yNum, kernel, xCentral, yCentral, x, y, options->mode);
    240262
    241263    // Image interpolation, according to image type
    … …  
    380402}
    381403
    382 // Interpolation engine for separable interpolation kernels (either for good reasons or for practical reasons)
    383 static psImageInterpolateStatus interpolateSeparate(double *imageValue, double *varianceValue,
    384                                                     psMaskType *maskValue, float x, float y,
    385                                                     const psImageInterpolateOptions *options)
    386 {
    387     // Parameters have been checked by psImageInterpolate()
    388 
    389     int xNum, yNum;                     // Number of interpolation kernel pixels
    390     switch (options->mode) {
     404// Setup for interpolation by separable kernels
     405static inline void interpolateSeparateSetup(int *xNum, int *yNum, // Size of interpolation kernel, returned
     406                                            int *xCentral, int *yCentral, // Central pixel of convolution
     407                                            float x, float y, // Coordinates of interest
     408                                            psImageInterpolateMode mode // Mode for interpolation
     409                                            )
     410{
     411    // Central pixel is the pixel below the point of interest
     412    *xCentral = floor(x - 0.5);
     413    *yCentral = floor(y - 0.5);
     414    switch (mode) {
    391415      case PS_INTERPOLATE_LANCZOS2:
    392         xNum = yNum = 4;
     416        *xNum = *yNum = 4;
    393417        break;
    394418      case PS_INTERPOLATE_LANCZOS3:
    395         xNum = yNum = 6;
     419        *xNum = *yNum = 6;
    396420        break;
    397421      case PS_INTERPOLATE_LANCZOS4:
    398         xNum = yNum = 8;
     422        *xNum = *yNum = 8;
    399423        break;
    400424      case PS_INTERPOLATE_FLAT:
    … …  
    405429        psAbort("Invalid interpolation mode.");
    406430    }
    407 
    408     // Central pixel is the pixel below the point of interest
    409     int xCentral = floor(x - 0.5), yCentral = floor(y - 0.5); // Central pixel of the convolution
    410     const psImage *image = options->image; // Image of interest
    411     int xLast = image->numCols - 1;     // Last pixel in x
    412     int yLast = image->numRows - 1;     // Last pixel in y
    413 
    414     if (xCentral - (xNum - 1) / 2 < 0 || xCentral + xNum / 2 > xLast ||
    415         yCentral - (yNum - 1) / 2 < 0 || yCentral + yNum / 2 > yLast) {
    416         // At least one pixel of the interpolation kernel is off the image
    417         if (imageValue) {
    418             *imageValue = options->badImage;
    419         }
    420         if (varianceValue) {
    421             *varianceValue = options->badVariance;
    422         }
    423         if (maskValue) {
    424             *maskValue = options->badMask;
    425         }
    426         return PS_INTERPOLATE_STATUS_OFF;
    427     }
    428 
    429 //    bool xExact, yExact;                // Is the shift exactly on?
    430     double xKernel[xNum], yKernel[yNum];// Interpolation kernels in x and y
    431     switch (options->mode) {
     431}
     432
     433// Generate the interpolation kernels for separable case; they should be normalised to unity
     434static inline void interpolateSeparateGenerate(int xNum, int yNum, // Size of interpolation kernel
     435                                               double xKernel[xNum], double yKernel[yNum], // Kernels
     436                                               int xCentral, int yCentral, // Central pixel of convolution
     437                                               float x, float y, // Coordinates of interest
     438                                               psImageInterpolateMode mode // Mode for interpolation
     439                                               )
     440{
     441    // XXX Could put in an "exact shift" (i.e., xFrac = 0.0) version
     442    switch (mode) {
    432443      case PS_INTERPOLATE_LANCZOS2:
    433444      case PS_INTERPOLATE_LANCZOS3:
    434445      case PS_INTERPOLATE_LANCZOS4: {
    435446          double xFrac = x - xCentral - 0.5; // Fraction of pixel in x
    436 #if 0
    437           if (fabs(xFrac) < DBL_EPSILON) {
    438               xExact = true;
    439           } else {
    440 #endif
    441               lanczos(xKernel, xNum, xFrac);
    442 #if 0
    443               xExact = false;
    444           }
    445 #endif
     447          lanczos(xKernel, xNum, xFrac);
    446448          double yFrac = y - yCentral - 0.5; // Fraction of pixel in y
    447 #if 0
    448           if (fabs(yFrac) < DBL_EPSILON) {
    449               yExact = true;
    450           } else {
    451 #endif
    452               lanczos(yKernel, yNum, yFrac);
    453 #if 0
    454               yExact = false;
    455           }
    456 #endif
     449          lanczos(yKernel, yNum, yFrac);
    457450          break;
    458451      }
    … …  
    464457        psAbort("Invalid interpolation mode.");
    465458    }
     459}
     460
     461// Interpolation engine for separable interpolation kernels (either for good reasons or for practical reasons)
     462static psImageInterpolateStatus interpolateSeparate(double *imageValue, double *varianceValue,
     463                                                    psMaskType *maskValue, float x, float y,
     464                                                    const psImageInterpolateOptions *options)
     465{
     466    // Parameters have been checked by psImageInterpolate()
     467
     468    int xNum, yNum;                     // Number of interpolation kernel pixels
     469    int xCentral, yCentral; // Central pixel of the convolution
     470    interpolateSeparateSetup(&xNum, &yNum, &xCentral, &yCentral, x, y, options->mode);
     471
     472    const psImage *image = options->image; // Image of interest
     473    int xLast = image->numCols - 1;     // Last pixel in x
     474    int yLast = image->numRows - 1;     // Last pixel in y
     475
     476    if (xCentral - (xNum - 1) / 2 < 0 || xCentral + xNum / 2 > xLast ||
     477        yCentral - (yNum - 1) / 2 < 0 || yCentral + yNum / 2 > yLast) {
     478        // At least one pixel of the interpolation kernel is off the image
     479        if (imageValue) {
     480            *imageValue = options->badImage;
     481        }
     482        if (varianceValue) {
     483            *varianceValue = options->badVariance;
     484        }
     485        if (maskValue) {
     486            *maskValue = options->badMask;
     487        }
     488        return PS_INTERPOLATE_STATUS_OFF;
     489    }
     490
     491    double xKernel[xNum], yKernel[yNum]; // Interpolation kernels
     492    interpolateSeparateGenerate(xNum, yNum, xKernel, yKernel, xCentral, yCentral, x, y, options->mode);
    466493
    467494    // Image interpolation, according to image type
    … …  
    631658}
    632659
     660
     661static float varianceFactorFlat(float x, float y, const psImageInterpolateOptions *options)
     662{
     663    // There's no smearing
     664    return 1.0;
     665}
     666
     667static float varianceFactorKernel(float x, float y, const psImageInterpolateOptions *options)
     668{
     669    int xNum, yNum;                     // Number of interpolation kernel pixels
     670    int xCentral, yCentral;             // Central pixel of the convolution
     671    interpolateKernelSetup(&xNum, &yNum, &xCentral, &yCentral, x, y, options->mode);
     672
     673    const psImage *image = options->image; // Image of interest
     674    int xLast = image->numCols - 1;     // Last pixel in x
     675    int yLast = image->numRows - 1;     // Last pixel in y
     676
     677    if (xCentral - (xNum - 1) / 2 < 0 || xCentral + xNum / 2 > xLast ||
     678        yCentral - (yNum - 1) / 2 < 0 || yCentral + yNum / 2 > yLast) {
     679        // At least one pixel of the interpolation kernel is off the image
     680        return NAN;
     681    }
     682
     683    double kernel[yNum][xNum];          // Interpolation kernel for straight interpolation
     684    interpolateKernelGenerate(xNum, yNum, kernel, xCentral, yCentral, x, y, options->mode);
     685
     686    double sumKernel2 = 0.0;            // Sum of kernel squares
     687    for (int j = 0; j < yNum; j++) {
     688        for (int i = 0; i < xNum; i++) {
     689            sumKernel2 += PS_SQR(kernel[j][i]);
     690        }
     691    }
     692
     693    return sumKernel2;
     694}
     695
     696static float varianceFactorSeparate(float x, float y, const psImageInterpolateOptions *options)
     697{
     698    int xNum, yNum;                     // Number of interpolation kernel pixels
     699    int xCentral, yCentral;             // Central pixel of the convolution
     700    interpolateSeparateSetup(&xNum, &yNum, &xCentral, &yCentral, x, y, options->mode);
     701
     702    const psImage *image = options->image; // Image of interest
     703    int xLast = image->numCols - 1;     // Last pixel in x
     704    int yLast = image->numRows - 1;     // Last pixel in y
     705
     706    if (xCentral - (xNum - 1) / 2 < 0 || xCentral + xNum / 2 > xLast ||
     707        yCentral - (yNum - 1) / 2 < 0 || yCentral + yNum / 2 > yLast) {
     708        // At least one pixel of the interpolation kernel is off the image
     709        return NAN;
     710    }
     711
     712    double xKernel[xNum], yKernel[yNum]; // Interpolation kernels for separable interpolation
     713    interpolateSeparateGenerate(xNum, yNum, xKernel, yKernel, xCentral, yCentral, x, y, options->mode);
     714
     715    double ySumKernel2 = 0.0;           // Sum of kernel squared in y
     716    for (int j = 0; j < yNum; j++) {
     717        double xSumKernel2 = 0.0;       // Sum of kernel squared in x
     718        for (int i = 0; i < xNum; i++) {
     719            xSumKernel2 += PS_SQR(xKernel[i]);
     720        }
     721        ySumKernel2 += xSumKernel2 * PS_SQR(yKernel[j]);
     722    }
     723
     724    return ySumKernel2;
     725}
     726
     727float psImageInterpolateVarianceFactor(float x, float y, const psImageInterpolateOptions *options)
     728{
     729    PS_ASSERT_PTR_NON_NULL(options, PS_INTERPOLATE_STATUS_ERROR);
     730
     731    switch (options->mode) {
     732      case PS_INTERPOLATE_FLAT:
     733        return varianceFactorFlat(x, y, options);
     734      case PS_INTERPOLATE_BILINEAR:
     735      case PS_INTERPOLATE_BICUBE:
     736      case PS_INTERPOLATE_GAUSS:
     737        return varianceFactorKernel(x, y, options);
     738      case PS_INTERPOLATE_LANCZOS2:
     739      case PS_INTERPOLATE_LANCZOS3:
     740      case PS_INTERPOLATE_LANCZOS4:
     741        return varianceFactorSeparate(x, y, options);
     742      default:
     743        psError(PS_ERR_BAD_PARAMETER_VALUE, true,
     744                _("Specified interpolation method (%d) is not supported."),
     745                options->mode);
     746        return NAN;
     747    }
     748
     749    psAbort("Should never reach here.");
     750    return NAN;
     751}
     752
     753
    633754psImageInterpolateMode psImageInterpolateModeFromString(const char *name)
    634755{

  • trunk/psLib/src/imageops/psImageInterpolate.h

    r14452 r18162  
    77 * @author Paul Price, Institute for Astronomy
    88 *
    9  * @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
    10  * @date $Date: 2007-08-09 01:40:07 $
     9 * @version $Revision: 1.4 $ $Name: not supported by cvs2svn $
     10 * @date $Date: 2008-06-17 21:24:58 $
    1111 * Copyright 2004-2007 Institute for Astronomy, University of Hawaii
    1212 */
    … …  
    8484    );
    8585
     86/// Return the variance factor for the appropriate position
     87///
     88/// psImageInterpolate sets the variance appropriate for extended regions (on the scale of the interpolation
     89/// kernel), but this is not appropriate for pixel-to-pixel statistics.  This function returns the conversion
     90/// factor.
     91float psImageInterpolateVarianceFactor(float x, float y, ///< Position of interest
     92                                       const psImageInterpolateOptions *options ///< Interpolation options
     93    );
    8694
    8795#endif
Note: See TracChangeset for help on using the changeset viewer.