IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Aug 19, 2008, 6:06:31 PM (18 years ago)
Author:
Paul Price
Message:

Changing how interpolation works with masked pixels. Previously, the pixel was defined as off the image if any part of the interpolation kernel was off the image; also, the interpolation kernel was applied to the image *without masking*, and the result of interpolation was flagged as bad/poor depending on the fractional contribution of the masked pixels. These very conservative approaches can result in significantly inflating bad pixel regions. Now, the interpolation kernel is applied, masking bad pixels, and pixels are flagged as bad/poor depending on the fractional contribution of the square of the kernel (i.e., according to its effect on the variance --- does the lack of data blow up the variance?). This also allows us to do a bit better at the edges of the image --- we no longer give up as soon as the interpolation kernel is touching the edge, but we measure the effect and return bad/poor accordingly. This made the test (tap_psImageInterpolate2.c) a bit trickier since one doesn't know whether a pixel should be bad or poor without doing the calculation, which is exactly what we're testing. Decided not to check the value, just the returned result (bad/poor) for these cases. Test passes; whether the results at the edge are correct or not, I don't know, but at least they should be masked appropriately.

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r18300 r19129  
    77 *  @author Paul Price, IfA
    88 *
    9  *  @version $Revision: 1.16 $ $Name: not supported by cvs2svn $
    10  *  @date $Date: 2008-06-24 02:04:55 $
     9 *  @version $Revision: 1.17 $ $Name: not supported by cvs2svn $
     10 *  @date $Date: 2008-08-20 04:06:31 $
    1111 *
    1212 *  Copyright 2004-2007 Institute for Astronomy, University of Hawaii
    … …  
    162162}
    163163
    164 // Generate the interpolation kernel; it should be normalised to unity
     164// Generate the interpolation kernel
     165// No need to normalise to unity
    165166static inline void interpolateKernelGenerate(int xNum, int yNum, // Size of interpolation kernel
    166167                                             double kernel[xNum][yNum], // Kernel, to be set
    … …  
    204205          double yFrac = y - yCentral - 0.5; // Fraction of pixel in y
    205206          double sigma = 0.5;           // Gaussian sigma
    206           double norm = 0.0;            // Normalisation
    207207          for (int j = 0, yPos = - (yNum - 1) / 2; j < yNum; j++, yPos++) {
    208208              for (int i = 0, xPos = - (xNum - 1) / 2; i < xNum; i++, xPos++) {
    209                   norm += kernel[j][i] = exp(-0.5 / PS_SQR(sigma) * (PS_SQR(xPos - xFrac) +
    210                                                                      PS_SQR(yPos - yFrac)));
    211               }
    212           }
    213           norm = 1.0 / norm;
    214           for (int j = 0; j < yNum; j++) {
    215               for (int i = 0; i < xNum; i++) {
    216                   kernel[j][i] *= norm;
     209                  kernel[j][i] = exp(-0.5 / PS_SQR(sigma) * (PS_SQR(xPos - xFrac) +
     210                                                             PS_SQR(yPos - yFrac)));
    217211              }
    218212          }
    … …  
    228222}
    229223
     224// Set up and check interpolation bounds
     225// This macro defines many useful values
     226#define INTERPOLATE_BOUNDS() \
     227    int xLast = image->numCols - 1, yLast = image->numRows - 1; /* Last pixels of image */ \
     228    /* Start and stop of kernel on image */ \
     229    int xStart = xCentral - (xNum - 1) / 2, xStop = xCentral + xNum / 2; \
     230    int yStart = yCentral - (yNum - 1) / 2, yStop = yCentral + yNum / 2; \
     231    if (xStart > xLast || xStop < 0 || yStart > yLast || yStop < 0) { \
     232        /* Interpolation kernel is entirely off the image */ \
     233        *imageValue = options->badImage; \
     234        if (varianceValue) { \
     235            *varianceValue = options->badVariance; \
     236        } \
     237        if (maskValue) { \
     238            *maskValue = options->badMask; \
     239        } \
     240        return PS_INTERPOLATE_STATUS_OFF; \
     241    } \
     242    int xMin, xMax, yMin, yMax; /* Minimum and maximum valid pixels on image */ \
     243    int iMin, iMax, jMin, jMax; /* Minimum and maximum valid pixels on kernel */ \
     244    bool offImage = false; /* At least one pixel of the kernel is off the image */ \
     245    if (xStart < 0) { \
     246        xMin = 0; \
     247        iMin = -xStart; \
     248        offImage = true; \
     249    } else { \
     250        xMin = xStart; \
     251        iMin = 0; \
     252    } \
     253    if (xStop > xLast) { \
     254        xMax = xLast; \
     255        iMax = xStop - xLast; \
     256        offImage = true; \
     257    } else { \
     258        xMax = xStop; \
     259        iMax = xNum; \
     260    } \
     261    if (yStart < 0) { \
     262        yMin = 0; \
     263        jMin = -yStart; \
     264        offImage = true; \
     265    } else { \
     266        yMin = yStart; \
     267        jMin = 0; \
     268    } \
     269    if (yStop > yLast) { \
     270        yMax = yLast; \
     271        jMax = yStop - yLast; \
     272        offImage = true; \
     273    } else { \
     274        yMax = yStop; \
     275        jMax = yNum; \
     276    }
     277
    230278// Interpolation engine using interpolation kernel
    231279static psImageInterpolateStatus interpolateKernel(double *imageValue, double *varianceValue,
    … …  
    240288
    241289    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     }
     290    const psImage *mask = options->mask; // Image mask
     291    const psImage *variance = options->variance; // Image variance
     292
     293    INTERPOLATE_BOUNDS();
    259294
    260295    double kernel[yNum][xNum];          // Interpolation kernel for straight interpolation
    261296    interpolateKernelGenerate(xNum, yNum, kernel, xCentral, yCentral, x, y, options->mode);
    262297
    263     // Image interpolation, according to image type
    264     #define KERNEL_IMAGE_CASE(TYPE) \
    265         case PS_TYPE_##TYPE: { \
    266             for (int j = 0, yPix = yCentral - (yNum - 1) / 2; j < yNum; j++, yPix++) { \
    267                 for (int i = 0, xPix = xCentral - (xNum - 1) / 2; i < xNum; i++, xPix++) { \
    268                     *imageValue += kernel[j][i] * image->data.TYPE[yPix][xPix]; \
    269                 } \
    270             } \
    271             break; \
    272         }
    273 
    274     // Calculate the value for the image
    275     if (imageValue) {
    276         *imageValue = 0.0;
    277         switch (image->type.type) {
    278             KERNEL_IMAGE_CASE(U8);
    279             KERNEL_IMAGE_CASE(U16);
    280             KERNEL_IMAGE_CASE(U32);
    281             KERNEL_IMAGE_CASE(U64);
    282             KERNEL_IMAGE_CASE(S8);
    283             KERNEL_IMAGE_CASE(S16);
    284             KERNEL_IMAGE_CASE(S32);
    285             KERNEL_IMAGE_CASE(S64);
    286             KERNEL_IMAGE_CASE(F32);
    287             KERNEL_IMAGE_CASE(F64);
    288           default:
    289             psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Unrecognised type for image: %x",
    290                     image->type.type);
    291             return PS_INTERPOLATE_STATUS_ERROR;
    292         }
    293     }
    294 
    295     // Check the mask value
    296     const psImage *mask = options->mask; // Image mask
    297     psImageInterpolateStatus status = PS_INTERPOLATE_STATUS_GOOD; // Status of interpolation
    298     if (maskValue) {
    299         *maskValue = 0;
    300         if (mask) {
    301             psMaskType maskVal = options->maskVal; // Mask value
    302             double badContrib = 0.0;        // Amount of kernel on bad pixels
    303             double totContrib = 0.0;        // Total amount of kernel
    304             bool badPix = false;            // Is there a bad pixel?
    305             for (int j = 0, yPix = yCentral - (yNum - 1) / 2; j < yNum; j++, yPix++) {
    306                 for (int i = 0, xPix = xCentral - (xNum - 1) / 2; i < xNum; i++, xPix++) {
    307                     if (mask->data.PS_TYPE_MASK_DATA[yPix][xPix] & maskVal) {
    308                         badContrib += fabs(kernel[j][i]);
    309                         badPix = true;
    310                     }
    311                     *maskValue |= mask->data.PS_TYPE_MASK_DATA[yPix][xPix];
    312                     totContrib += fabs(kernel[j][i]);
    313                 }
     298    float sumKernel = 0.0;              // Sum of the kernel
     299    double sumKernel2 = 0.0;            // Sum of the kernel-squared
     300    float sumBad = 0.0;                 // Sum of bad kernel-squared contributions
     301    psMaskType maskVal = options->maskVal; // Value to mask
     302    double sumImage = 0.0;              // Sum of image multiplied by kernel
     303    double sumVariance = 0.0;           // Sum of variance multiplied by kernel-squared
     304
     305    bool wantVariance = variance && varianceValue; // Does the user want the variance value?
     306    bool haveMask = mask && maskVal; // Does the user want the variance value?
     307
     308    // Add contributions in an area outside the image
     309#define INTERPOLATE_KERNEL_ADD_OFFIMAGE() { \
     310        float kernelValue = kernel[j][i]; /* Value of kernel */ \
     311        double kernelValue2 = PS_SQR(kernelValue); /* Square of kernel */ \
     312        sumBad += kernelValue2; \
     313        sumKernel2 += kernelValue2; \
     314    }
     315
     316    // Measure kernel contribution from outside image
     317    if (offImage) {
     318        // Bottom rows
     319        for (int j = 0; j < jMin; j++) {
     320            for (int i = 0; i < xNum; i++) {
     321                INTERPOLATE_KERNEL_ADD_OFFIMAGE();
    314322            }
    315 
    316             if (badPix) {
    317                 if (badContrib / totContrib >= options->poorFrac) {
    318                     *maskValue |= options->badMask;
    319                     status = PS_INTERPOLATE_STATUS_BAD;
    320                 } else {
    321                     *maskValue |= options->poorMask;
    322                     status = PS_INTERPOLATE_STATUS_POOR;
    323                 }
     323        }
     324        // Two sides
     325        for (int j = jMin; j < jMax; j++) {
     326            for (int i = 0; i < iMin; i++) {
     327                INTERPOLATE_KERNEL_ADD_OFFIMAGE();
    324328            }
    325         }
    326     }
    327 
    328     // Finally, the variance
    329     const psImage *variance = options->variance; // Image variance
    330     if (varianceValue && variance) {
    331         *varianceValue = 0.0;
    332 
    333         // Variance interpolation, according to image type
    334         #define KERNEL_VARIANCE_CASE(TYPE) \
    335             case PS_TYPE_##TYPE: { \
    336                 double sumKernel2 = 0.0; /* Sum of kernel squares */ \
    337                 for (int j = 0, yPix = yCentral - (yNum - 1) / 2; j < yNum; j++, yPix++) { \
    338                     for (int i = 0, xPix = xCentral - (xNum - 1) / 2; i < xNum; i++, xPix++) { \
    339                         double kernel2 = PS_SQR(kernel[j][i]); /* Kernel squared */ \
    340                         sumKernel2 += kernel2; \
    341                         *varianceValue += kernel2 * variance->data.TYPE[yPix][xPix]; \
    342                     } \
    343                 } \
    344                 *varianceValue /= sumKernel2; /* Normalise so that sum of kernel squares is unity */ \
    345                 break; \
     329            for (int i = iMax + 1; i < xNum; i++) {
     330                INTERPOLATE_KERNEL_ADD_OFFIMAGE();
    346331            }
    347 
    348         switch (variance->type.type) {
    349             KERNEL_VARIANCE_CASE(U8);
    350             KERNEL_VARIANCE_CASE(U16);
    351             KERNEL_VARIANCE_CASE(U32);
    352             KERNEL_VARIANCE_CASE(U64);
    353             KERNEL_VARIANCE_CASE(S8);
    354             KERNEL_VARIANCE_CASE(S16);
    355             KERNEL_VARIANCE_CASE(S32);
    356             KERNEL_VARIANCE_CASE(S64);
    357             KERNEL_VARIANCE_CASE(F32);
    358             KERNEL_VARIANCE_CASE(F64);
    359           default:
    360             psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Unrecognised type for image: %x",
    361                     image->type.type);
    362             return PS_INTERPOLATE_STATUS_ERROR;
    363         }
     332        }
     333        // Top rows
     334        for (int j = jMax + 1; j < yNum; j++) {
     335            for (int i = 0; i < xNum; i++) {
     336                INTERPOLATE_KERNEL_ADD_OFFIMAGE();
     337            }
     338        }
     339    }
     340
     341#define INTERPOLATE_KERNEL_CASE(TYPE) \
     342  case PS_TYPE_##TYPE: { \
     343      if (wantVariance) { \
     344          if (haveMask) { \
     345              /* Variance and mask */ \
     346              for (int j = jMin, yPix = yMin; j < jMax; j++, yPix++) { \
     347                  for (int i = iMin, xPix = xMin; i < iMax; i++, xPix++) { \
     348                      float kernelValue = kernel[j][i]; /* Value of kernel */ \
     349                      double kernelValue2 = PS_SQR(kernelValue); /* Square of kernel */ \
     350                      if (mask->data.PS_TYPE_MASK_DATA[yPix][xPix] & maskVal) { \
     351                          sumBad += kernelValue2; \
     352                      } else { \
     353                          sumImage += kernelValue * image->data.TYPE[yPix][xPix]; \
     354                          sumVariance += kernelValue2 * variance->data.TYPE[yPix][xPix]; \
     355                          sumKernel += kernelValue; \
     356                      } \
     357                      sumKernel2 += PS_SQR(kernelValue); \
     358                  } \
     359              } \
     360              \
     361          } else { \
     362              /* Variance, no mask */ \
     363              for (int j = jMin, yPix = yMin; j < jMax; j++, yPix++) { \
     364                  for (int i = iMin, xPix = xMin; i < iMax; i++, xPix++) { \
     365                      float kernelValue = kernel[j][i]; /* Value of kernel */ \
     366                      double kernelValue2 = PS_SQR(kernelValue); /* Square of kernel */ \
     367                      sumImage += kernelValue * image->data.TYPE[yPix][xPix]; \
     368                      sumVariance += kernelValue2 * variance->data.TYPE[yPix][xPix]; \
     369                      sumKernel += kernelValue; \
     370                      sumKernel2 += kernelValue2; \
     371                  } \
     372              } \
     373          } \
     374      } else if (haveMask) { \
     375          /* Mask, no variance */ \
     376          for (int j = jMin, yPix = yMin; j < jMax; j++, yPix++) { \
     377              for (int i = iMin, xPix = xMin; i < iMax; i++, xPix++) { \
     378                  float kernelValue = kernel[j][i]; /* Value of kernel */ \
     379                  if (mask->data.PS_TYPE_MASK_DATA[yPix][xPix] & maskVal) { \
     380                      sumBad += PS_SQR(kernelValue); \
     381                  } else { \
     382                      sumImage += kernelValue * image->data.TYPE[yPix][xPix]; \
     383                      sumKernel += kernelValue; \
     384                  } \
     385                  sumKernel2 += PS_SQR(kernelValue); \
     386              } \
     387          } \
     388      } else { \
     389          /* Neither variance nor mask */ \
     390          for (int j = jMin, yPix = yMin; j < jMax; j++, yPix++) { \
     391              for (int i = iMin, xPix = xMin; i < iMax; i++, xPix++) { \
     392                  float kernelValue = kernel[j][i]; /* Value of kernel */ \
     393                  sumImage += kernelValue * image->data.TYPE[yPix][xPix]; \
     394                  sumKernel += kernelValue; \
     395              } \
     396          } \
     397      } \
     398    } \
     399    break;
     400
     401
     402    switch (image->type.type) {
     403        INTERPOLATE_KERNEL_CASE(U8);
     404        INTERPOLATE_KERNEL_CASE(U16);
     405        INTERPOLATE_KERNEL_CASE(U32);
     406        INTERPOLATE_KERNEL_CASE(U64);
     407        INTERPOLATE_KERNEL_CASE(S8);
     408        INTERPOLATE_KERNEL_CASE(S16);
     409        INTERPOLATE_KERNEL_CASE(S32);
     410        INTERPOLATE_KERNEL_CASE(S64);
     411        INTERPOLATE_KERNEL_CASE(F32);
     412        INTERPOLATE_KERNEL_CASE(F64);
     413      default:
     414        psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Unrecognised type for image: %x",
     415                image->type.type);
     416        return PS_INTERPOLATE_STATUS_ERROR;
     417    }
     418
     419    psImageInterpolateStatus status = PS_INTERPOLATE_STATUS_ERROR; // Status of interpolation
     420    *imageValue = sumImage / sumKernel;
     421    if (wantVariance) {
     422        *varianceValue = sumVariance * PS_SQR(sumKernel) / sumKernel2;
     423    }
     424    if (sumBad == 0) {
     425        // Completely good pixel
     426        status = PS_INTERPOLATE_STATUS_GOOD;
     427    } else if (sumBad / sumKernel2 < PS_SQR(options->poorFrac)) {
     428        // Some pixels masked: poor pixel
     429        if (haveMask && maskValue) {
     430            *maskValue |= options->poorMask;
     431        }
     432        status = PS_INTERPOLATE_STATUS_POOR;
     433    } else {
     434        // Many pixels (or a few important pixels) masked: bad pixel
     435        if (haveMask && maskValue) {
     436            *maskValue |= options->badMask;
     437        }
     438        status = PS_INTERPOLATE_STATUS_BAD;
    364439    }
    365440
    … …  
    472547
    473548    const psImage *image = options->image; // Image of interest
    474     int xLast = image->numCols - 1;     // Last pixel in x
    475     int yLast = image->numRows - 1;     // Last pixel in y
    476 
    477     if (xCentral - (xNum - 1) / 2 < 0 || xCentral + xNum / 2 > xLast ||
    478         yCentral - (yNum - 1) / 2 < 0 || yCentral + yNum / 2 > yLast) {
    479         // At least one pixel of the interpolation kernel is off the image
    480         if (imageValue) {
    481             *imageValue = options->badImage;
    482         }
    483         if (varianceValue) {
    484             *varianceValue = options->badVariance;
    485         }
    486         if (maskValue) {
    487             *maskValue = options->badMask;
    488         }
    489         return PS_INTERPOLATE_STATUS_OFF;
    490     }
     549    const psImage *mask = options->mask; // Image mask
     550    const psImage *variance = options->variance; // Image variance
     551
     552    INTERPOLATE_BOUNDS();
    491553
    492554    double xKernel[xNum], yKernel[yNum]; // Interpolation kernels
    493555    interpolateSeparateGenerate(xNum, yNum, xKernel, yKernel, xCentral, yCentral, x, y, options->mode);
    494556
    495     // Image interpolation, according to image type
    496     #define SEPARATE_IMAGE_CASE(TYPE) \
    497       case PS_TYPE_##TYPE: { \
    498         for (int j = 0, yPix = yCentral - (yNum - 1) / 2; j < yNum; j++, yPix++) { \
    499             double xInterpValue = 0.0; /* Interpolation in x */ \
    500             for (int i = 0, xPix = xCentral - (xNum - 1) / 2; i < xNum; i++, xPix++) { \
    501                 xInterpValue += xKernel[i] * image->data.TYPE[yPix][xPix]; \
    502             } \
    503             *imageValue += yKernel[j] * xInterpValue; /* Interpolating in y */ \
    504         } \
    505         break; \
    506       }
    507 
    508     // Calculate the value for the image
    509     if (imageValue) {
    510         *imageValue = 0.0;
    511         switch (image->type.type) {
    512             SEPARATE_IMAGE_CASE(U8);
    513             SEPARATE_IMAGE_CASE(U16);
    514             SEPARATE_IMAGE_CASE(U32);
    515             SEPARATE_IMAGE_CASE(U64);
    516             SEPARATE_IMAGE_CASE(S8);
    517             SEPARATE_IMAGE_CASE(S16);
    518             SEPARATE_IMAGE_CASE(S32);
    519             SEPARATE_IMAGE_CASE(S64);
    520             SEPARATE_IMAGE_CASE(F32);
    521             SEPARATE_IMAGE_CASE(F64);
    522           default:
    523             psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Unrecognised type for image: %x",
    524                     image->type.type);
    525             return PS_INTERPOLATE_STATUS_ERROR;
    526         }
    527     }
    528 
    529     // Check the mask value
    530     const psImage *mask = options->mask; // Image mask
    531     psImageInterpolateStatus status = PS_INTERPOLATE_STATUS_GOOD; // Status of interpolation
    532     if (maskValue) {
    533         *maskValue = 0;
    534         if (mask) {
    535             psMaskType maskVal = options->maskVal; // Mask value
    536             double badContrib = 0.0;        // Amount of kernel on bad pixels
    537             double totContrib = 0.0;        // Total amount of kernel
    538             bool badPix = false;            // Number of bad pixels
    539             for (int j = 0, yPix = yCentral - (yNum - 1) / 2; j < yNum; j++, yPix++) {
    540                 // Interpolation in x
    541                 double xBadContrib = 0.0;   // Amount of x kernel on bad pixels
    542                 double xTotContrib = 0.0;   // Total amount of x kernel
    543                 for (int i = 0, xPix = xCentral - (xNum - 1) / 2; i < xNum; i++, xPix++) {
    544                     if (mask->data.PS_TYPE_MASK_DATA[yPix][xPix] & maskVal) {
    545                         xBadContrib += fabs(xKernel[i]);
    546                         badPix = true;
    547                     }
    548                     *maskValue |= mask->data.PS_TYPE_MASK_DATA[yPix][xPix];
    549                     xTotContrib += fabs(xKernel[i]);
    550                 }
    551                 // Interpolating in y
    552                 badContrib += fabs(yKernel[j]) * xBadContrib;
    553                 totContrib += fabs(yKernel[j]) * xTotContrib;
     557    float sumKernel = 0.0;              // Sum of the kernel
     558    double sumKernel2 = 0.0;            // Sum of the kernel-squared
     559    float sumBad = 0.0;                 // Sum of bad kernel-squared contributions
     560    psMaskType maskVal = options->maskVal; // Value to mask
     561    double sumImage = 0.0;              // Sum of image multiplied by kernel
     562    double sumVariance = 0.0;           // Sum of variance multiplied by kernel-squared
     563
     564    bool wantVariance = variance && varianceValue; // Does the user want the variance value?
     565    bool haveMask = mask && maskVal; // Does the user want the variance value?
     566
     567    // Add contributions in an area outside the image
     568#define INTERPOLATE_SEPARATE_SETUP_OFFIMAGE_COL() \
     569    float xSumBad = 0.0; \
     570    double xSumKernel2 = 0.0;
     571#define INTERPOLATE_SEPARATE_ADD_OFFIMAGE_COL() { \
     572        float kernelValue = xKernel[i]; /* Value of kernel */ \
     573        double kernelValue2 = PS_SQR(kernelValue); /* Square of kernel */ \
     574        xSumBad += kernelValue2; \
     575        xSumKernel2 += kernelValue2; \
     576    }
     577#define INTERPOLATE_SEPARATE_ADD_OFFIMAGE_ROW() { \
     578        float kernelValue = yKernel[j]; /* Value of kernel */ \
     579        double kernelValue2 = PS_SQR(kernelValue); /* Square of kernel */ \
     580        sumBad += kernelValue2 * xSumBad; \
     581        sumKernel2 += kernelValue2 * xSumKernel2; \
     582    }
     583
     584    // Measure kernel contribution from outside image
     585    if (offImage) {
     586        // Bottom rows
     587        for (int j = 0; j < jMin; j++) {
     588            INTERPOLATE_SEPARATE_SETUP_OFFIMAGE_COL();
     589            for (int i = 0; i < xNum; i++) {
     590                INTERPOLATE_SEPARATE_ADD_OFFIMAGE_COL();
    554591            }
    555 
    556             if (badPix) {
    557                 if (badContrib / totContrib >= options->poorFrac) {
    558                     *maskValue |= options->badMask;
    559                     status = PS_INTERPOLATE_STATUS_BAD;
    560                 } else {
    561                     *maskValue |= options->poorMask;
    562                     status = PS_INTERPOLATE_STATUS_POOR;
    563                 }
     592            INTERPOLATE_SEPARATE_ADD_OFFIMAGE_ROW();
     593        }
     594        // Two sides
     595        for (int j = jMin; j < jMax; j++) {
     596            INTERPOLATE_SEPARATE_SETUP_OFFIMAGE_COL();
     597            for (int i = 0; i < iMin; i++) {
     598                INTERPOLATE_SEPARATE_ADD_OFFIMAGE_COL();
    564599            }
    565         }
    566     }
    567 
    568     // Finally, the variance
    569     const psImage *variance = options->variance; // Image variance
    570     if (varianceValue && variance) {
    571         *varianceValue = 0.0;
    572 
    573         // Variance interpolation, according to image type
    574         #define SEPARATE_VARIANCE_CASE(TYPE) \
    575           case PS_TYPE_##TYPE: { \
    576             double ySumKernel2 = 0.0; /* Sum of kernel squared in y */ \
    577             for (int j = 0, yPix = yCentral - (yNum - 1) / 2; j < yNum; j++, yPix++) { \
    578                 double xSumKernel2 = 0.0; /* Sum of kernel squared in x */ \
    579                 double xInterpValue = 0.0; /* Interpolation in x */ \
    580                 for (int i = 0, xPix = xCentral - (xNum - 1) / 2; i < xNum; i++, xPix++) { \
    581                     double kernel2 = PS_SQR(xKernel[i]); /* Kernel squared */ \
    582                     xSumKernel2 += kernel2; \
    583                     xInterpValue += kernel2 * variance->data.TYPE[yPix][xPix]; \
    584                 } \
    585                 double kernel2 = PS_SQR(yKernel[j]); /* Kernel squared */ \
    586                 ySumKernel2 += xSumKernel2 * kernel2; \
    587                 *varianceValue += xInterpValue * kernel2; /* Interpolating in y */ \
    588             } \
    589             *varianceValue /= ySumKernel2; \
    590             break; \
    591           }
    592 
    593         switch (variance->type.type) {
    594             SEPARATE_VARIANCE_CASE(U8);
    595             SEPARATE_VARIANCE_CASE(U16);
    596             SEPARATE_VARIANCE_CASE(U32);
    597             SEPARATE_VARIANCE_CASE(U64);
    598             SEPARATE_VARIANCE_CASE(S8);
    599             SEPARATE_VARIANCE_CASE(S16);
    600             SEPARATE_VARIANCE_CASE(S32);
    601             SEPARATE_VARIANCE_CASE(S64);
    602             SEPARATE_VARIANCE_CASE(F32);
    603             SEPARATE_VARIANCE_CASE(F64);
    604           default:
    605             psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Unrecognised type for image: %x",
    606                     variance->type.type);
    607             return PS_INTERPOLATE_STATUS_ERROR;
    608         }
     600            for (int i = iMax + 1; i < xNum; i++) {
     601                INTERPOLATE_SEPARATE_ADD_OFFIMAGE_COL();
     602            }
     603            INTERPOLATE_SEPARATE_ADD_OFFIMAGE_ROW();
     604        }
     605        // Top rows
     606        for (int j = jMax + 1; j < yNum; j++) {
     607            INTERPOLATE_SEPARATE_SETUP_OFFIMAGE_COL();
     608            for (int i = 0; i < xNum; i++) {
     609                INTERPOLATE_SEPARATE_ADD_OFFIMAGE_COL();
     610            }
     611            INTERPOLATE_SEPARATE_ADD_OFFIMAGE_ROW();
     612        }
     613    }
     614
     615#define INTERPOLATE_SEPARATE_CASE(TYPE) \
     616  case PS_TYPE_##TYPE: { \
     617      if (wantVariance) { \
     618          if (haveMask) { \
     619              /* Variance and mask */ \
     620              for (int j = jMin, yPix = yMin; j < jMax; j++, yPix++) { \
     621                  double xSumImage = 0.0; /* Sum of image multiplied by kernel in x */ \
     622                  double xSumVariance = 0.0; /* Sum of image multiplied by kernel-squared in x */ \
     623                  float xSumKernel = 0.0; /* Sum of kernel in x */ \
     624                  double xSumKernel2 = 0.0; /* Sum of kernel-squared in x */ \
     625                  float xSumBad = 0.0; /* Sum of bad kernel-squared in x */ \
     626                  for (int i = iMin, xPix = xMin; i < iMax; i++, xPix++) { \
     627                      float kernelValue = xKernel[i]; /* Value of kernel in x */ \
     628                      double kernelValue2 = PS_SQR(kernelValue); /* Square of kernel in x */ \
     629                      if (mask->data.PS_TYPE_MASK_DATA[yPix][xPix] & maskVal) { \
     630                          xSumBad += kernelValue2; \
     631                      } else { \
     632                          xSumImage += kernelValue * image->data.TYPE[yPix][xPix]; \
     633                          xSumVariance += kernelValue2 * variance->data.TYPE[yPix][xPix]; \
     634                          xSumKernel += kernelValue; \
     635                      } \
     636                      xSumKernel2 += kernelValue2; \
     637                  } \
     638                  float kernelValue = yKernel[j]; /* Value of kernel in y */ \
     639                  double kernelValue2 = PS_SQR(kernelValue); /* Value of kernel-squared in y */ \
     640                  sumImage += kernelValue * xSumImage; \
     641                  sumVariance += kernelValue2 * xSumVariance; \
     642                  sumBad += kernelValue2 * xSumBad; \
     643                  sumKernel += kernelValue * xSumKernel; \
     644                  sumKernel2 += kernelValue2 * xSumKernel2; \
     645              } \
     646          } else { \
     647              /* Variance, no mask */ \
     648              for (int j = jMin, yPix = yMin; j < jMax; j++, yPix++) { \
     649                  double xSumImage = 0.0; /* Sum of image multiplied by kernel in x */ \
     650                  double xSumVariance = 0.0; /* Sum of image multiplied by kernel-squared in x */ \
     651                  float xSumKernel = 0.0; /* Sum of kernel in x */ \
     652                  double xSumKernel2 = 0.0; /* Sum of kernel-squared in x */ \
     653                  for (int i = iMin, xPix = xMin; i < iMax; i++, xPix++) { \
     654                      float kernelValue = xKernel[i]; /* Value of kernel */ \
     655                      double kernelValue2 = PS_SQR(kernelValue); /* Square of kernel */ \
     656                      xSumImage += kernelValue * image->data.TYPE[yPix][xPix]; \
     657                      xSumVariance += kernelValue2 * variance->data.TYPE[yPix][xPix]; \
     658                      xSumKernel += kernelValue; \
     659                      xSumKernel2 += kernelValue2; \
     660                  } \
     661                  float kernelValue = yKernel[j]; /* Value of kernel in y */ \
     662                  double kernelValue2 = PS_SQR(kernelValue); /* Value of kernel-squared in y */ \
     663                  sumImage += kernelValue * xSumImage; \
     664                  sumVariance += kernelValue2 * xSumVariance; \
     665                  sumKernel += kernelValue * xSumKernel; \
     666                  sumKernel2 += kernelValue2 * xSumKernel2; \
     667              } \
     668          } \
     669      } else if (haveMask) { \
     670          /* Mask, no variance */ \
     671          for (int j = jMin, yPix = yMin; j < jMax; j++, yPix++) { \
     672              double xSumImage = 0.0; /* Sum of image multiplied by kernel in x */ \
     673              float xSumKernel = 0.0; /* Sum of kernel in x */ \
     674              double xSumKernel2 = 0.0; /* Sum of kernel-squared in x */ \
     675              float xSumBad = 0.0; /* Sum of bad kernel-squared in x */ \
     676              for (int i = iMin, xPix = xMin; i < iMax; i++, xPix++) { \
     677                  float kernelValue = xKernel[i]; /* Value of kernel */ \
     678                  double kernelValue2 = PS_SQR(kernelValue); /* Value of kernel-squared */ \
     679                  if (mask->data.PS_TYPE_MASK_DATA[yPix][xPix] & maskVal) { \
     680                      xSumBad += kernelValue2; \
     681                  } else { \
     682                      xSumImage += kernelValue * image->data.TYPE[yPix][xPix]; \
     683                      xSumKernel += kernelValue; \
     684                  } \
     685                  xSumKernel2 += kernelValue2; \
     686              } \
     687              float kernelValue = yKernel[j]; /* Value of kernel in y */ \
     688              double kernelValue2 = PS_SQR(kernelValue); /* Value of kernel-squared in y */ \
     689              sumImage += kernelValue * xSumImage; \
     690              sumBad += kernelValue2 * xSumBad; \
     691              sumKernel += kernelValue * xSumKernel; \
     692              sumKernel2 += kernelValue2 * xSumKernel2; \
     693          } \
     694      } else {\
     695          /* Neither variance nor mask */ \
     696          for (int j = jMin, yPix = yMin; j < jMax; j++, yPix++) { \
     697              double xSumImage = 0.0; /* Sum of image multiplied by kernel in x */ \
     698              float xSumKernel = 0.0; /* Sum of kernel in x */ \
     699              for (int i = iMin, xPix = xMin; i < iMax; i++, xPix++) { \
     700                  float kernelValue = xKernel[i]; /* Value of kernel */ \
     701                  xSumImage += kernelValue * image->data.TYPE[yPix][xPix]; \
     702                  xSumKernel += kernelValue; \
     703              } \
     704              float kernelValue = yKernel[j]; /* Value of kernel in y */ \
     705              sumImage += kernelValue * xSumImage; \
     706              sumKernel += kernelValue * xSumKernel; \
     707          } \
     708      } \
     709    } \
     710    break;
     711
     712
     713    switch (image->type.type) {
     714        INTERPOLATE_SEPARATE_CASE(U8);
     715        INTERPOLATE_SEPARATE_CASE(U16);
     716        INTERPOLATE_SEPARATE_CASE(U32);
     717        INTERPOLATE_SEPARATE_CASE(U64);
     718        INTERPOLATE_SEPARATE_CASE(S8);
     719        INTERPOLATE_SEPARATE_CASE(S16);
     720        INTERPOLATE_SEPARATE_CASE(S32);
     721        INTERPOLATE_SEPARATE_CASE(S64);
     722        INTERPOLATE_SEPARATE_CASE(F32);
     723        INTERPOLATE_SEPARATE_CASE(F64);
     724      default:
     725        psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Unrecognised type for image: %x",
     726                image->type.type);
     727        return PS_INTERPOLATE_STATUS_ERROR;
     728    }
     729
     730    psImageInterpolateStatus status = PS_INTERPOLATE_STATUS_ERROR; // Status of interpolation
     731    *imageValue = sumImage / sumKernel;
     732    if (wantVariance) {
     733        *varianceValue = sumVariance * PS_SQR(sumKernel) / sumKernel2;
     734    }
     735    if (sumBad == 0) {
     736        // Completely good pixel
     737        status = PS_INTERPOLATE_STATUS_GOOD;
     738    } else if (sumBad / sumKernel2 < PS_SQR(options->poorFrac)) {
     739        // Some pixels masked: poor pixel
     740        if (haveMask && maskValue) {
     741            *maskValue |= options->poorMask;
     742        }
     743        status = PS_INTERPOLATE_STATUS_POOR;
     744    } else {
     745        // Many pixels (or a few important pixels) masked: bad pixel
     746        if (haveMask && maskValue) {
     747            *maskValue |= options->badMask;
     748        }
     749        status = PS_INTERPOLATE_STATUS_BAD;
    609750    }
    610751
    … …  
    617758                                            float x, float y, const psImageInterpolateOptions *options)
    618759{
     760    PS_ASSERT_PTR_NON_NULL(imageValue, PS_INTERPOLATE_STATUS_ERROR);
    619761    PS_ASSERT_PTR_NON_NULL(options, PS_INTERPOLATE_STATUS_ERROR);
    620762
Note: See TracChangeset for help on using the changeset viewer.