IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Oct 16, 2008, 12:55:45 PM (18 years ago)
Author:
Paul Price
Message:

Wasn't calculating the variance correctly, which led to psphot dramatically over-estimating the significance of sources.

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r20104 r20207  
    77 *  @author Paul Price, IfA
    88 *
    9  *  @version $Revision: 1.21 $ $Name: not supported by cvs2svn $
    10  *  @date $Date: 2008-10-13 21:47:41 $
     9 *  @version $Revision: 1.22 $ $Name: not supported by cvs2svn $
     10 *  @date $Date: 2008-10-16 22:55:45 $
    1111 *
    1212 *  Copyright 2004-2007 Institute for Astronomy, University of Hawaii
    … …  
    166166// No need to normalise to unity
    167167static inline void interpolateKernelGenerate(int xNum, int yNum, // Size of interpolation kernel
    168                                              double kernel[xNum][yNum], // Kernel, to be set
     168                                             float kernel[xNum][yNum], // Kernel, to be set
    169169                                             bool *xExact, bool *yExact, // Exact shift?
    170170                                             int xCentral, int yCentral, // Central pixel of convolution
    … …  
    174174                                             )
    175175{
    176     double xFrac = x - 0.5 - xCentral; // Fraction of pixel in x
    177     double yFrac = y - 0.5 - yCentral; // Fraction of pixel in y
     176    float xFrac = x - 0.5 - xCentral; // Fraction of pixel in x
     177    float yFrac = y - 0.5 - yCentral; // Fraction of pixel in y
    178178    *xExact = (allowExact && fabs(xFrac) < DBL_EPSILON);
    179179    *yExact = (allowExact && fabs(yFrac) < DBL_EPSILON);
    … …  
    188188      }
    189189      case PS_INTERPOLATE_BICUBE: {
    190           double xFrac = x - 0.5 - xCentral; // Fraction of pixel in x
    191           double yFrac = y - 0.5 - yCentral; // Fraction of pixel in y
     190          float xFrac = x - 0.5 - xCentral; // Fraction of pixel in x
     191          float yFrac = y - 0.5 - yCentral; // Fraction of pixel in y
    192192          // Calculation variables
    193193          double xxFrac = xFrac * xFrac / 6.0;
    … …  
    208208      }
    209209      case PS_INTERPOLATE_GAUSS: {
    210           double xFrac = x - xCentral - 0.5; // Fraction of pixel in x
    211           double yFrac = y - yCentral - 0.5; // Fraction of pixel in y
    212           double sigma = 0.5;           // Gaussian sigma
     210          float xFrac = x - xCentral - 0.5; // Fraction of pixel in x
     211          float yFrac = y - yCentral - 0.5; // Fraction of pixel in y
     212          float sigma = 0.5;           // Gaussian sigma
    213213          for (int j = 0, yPos = - (yNum - 1) / 2; j < yNum; j++, yPos++) {
    214214              for (int i = 0, xPos = - (xNum - 1) / 2; i < xNum; i++, xPos++) {
    215                   kernel[j][i] = exp(-0.5 / PS_SQR(sigma) * (PS_SQR(xPos - xFrac) +
    216                                                              PS_SQR(yPos - yFrac)));
     215                  kernel[j][i] = expf(-0.5 / PS_SQR(sigma) * (PS_SQR(xPos - xFrac) +
     216                                                              PS_SQR(yPos - yFrac)));
    217217              }
    218218          }
    … …  
    326326    INTERPOLATE_BOUNDS();
    327327
    328     double kernel[yNum][xNum];          // Interpolation kernel for straight interpolation
     328    float kernel[yNum][xNum];           // Interpolation kernel for straight interpolation
    329329    bool xExact, yExact;                // Exact shifts?
    330330    interpolateKernelGenerate(xNum, yNum, kernel, &xExact, &yExact, xCentral, yCentral,
    … …  
    345345    // Add contributions in an area outside the image
    346346#define INTERPOLATE_KERNEL_ADD_OFFIMAGE() { \
    347         float kernelValue = kernel[j][i]; /* Value of kernel */ \
    348         double kernelValue2 = PS_SQR(kernelValue); /* Square of kernel */ \
    349         sumBad += kernelValue2; \
    350         sumKernel2 += kernelValue2; \
     347        sumBad += PS_SQR(kernel[j][i]); \
    351348    }
    352349
    … …  
    390387                  for (int i = iMin, xPix = xMin; i < iMax; i++, xPix++) { \
    391388                      float kernelValue = kernel[j][i]; /* Value of kernel */ \
    392                       double kernelValue2 = PS_SQR(kernelValue); /* Square of kernel */ \
     389                      float kernelValue2 = PS_SQR(kernelValue); /* Square of kernel */ \
    393390                      if (mask->data.PS_TYPE_MASK_DATA[yPix][xPix] & maskVal) { \
    394391                          sumBad += kernelValue2; \
    … …  
    397394                          sumVariance += kernelValue2 * variance->data.TYPE[yPix][xPix]; \
    398395                          sumKernel += kernelValue; \
     396                          sumKernel2 += kernelValue2; \
    399397                      } \
    400                       sumKernel2 += PS_SQR(kernelValue); \
    401398                  } \
    402399              } \
    … …  
    407404                  for (int i = iMin, xPix = xMin; i < iMax; i++, xPix++) { \
    408405                      float kernelValue = kernel[j][i]; /* Value of kernel */ \
    409                       double kernelValue2 = PS_SQR(kernelValue); /* Square of kernel */ \
     406                      float kernelValue2 = PS_SQR(kernelValue); /* Square of kernel */ \
    410407                      sumImage += kernelValue * image->data.TYPE[yPix][xPix]; \
    411408                      sumVariance += kernelValue2 * variance->data.TYPE[yPix][xPix]; \
    … …  
    420417              for (int i = iMin, xPix = xMin; i < iMax; i++, xPix++) { \
    421418                  float kernelValue = kernel[j][i]; /* Value of kernel */ \
     419                  float kernelValue2 = PS_SQR(kernelValue); /* Square of kernel */ \
    422420                  if (mask->data.PS_TYPE_MASK_DATA[yPix][xPix] & maskVal) { \
    423                       sumBad += PS_SQR(kernelValue); \
     421                      sumBad += kernelValue2; \
    424422                  } else { \
    425423                      sumImage += kernelValue * image->data.TYPE[yPix][xPix]; \
    426424                      sumKernel += kernelValue; \
     425                      sumKernel2 += kernelValue2; \
    427426                  } \
    428                   sumKernel2 += PS_SQR(kernelValue); \
    429427              } \
    430428          } \
    … …  
    463461    *imageValue = sumImage / sumKernel;
    464462    if (wantVariance) {
    465         *varianceValue = sumVariance * PS_SQR(sumKernel) / sumKernel2;
     463        *varianceValue = sumVariance / sumKernel2;
    466464    }
    467465    if (sumBad == 0) {
    468466        // Completely good pixel
    469467        status = PS_INTERPOLATE_STATUS_GOOD;
    470     } else if (sumBad / sumKernel2 < PS_SQR(options->poorFrac)) {
     468    } else if (sumBad < PS_SQR(options->poorFrac) * (sumBad + sumKernel2)) {
    471469        // Some pixels masked: poor pixel
    472470        if (haveMask && maskValue) {
    … …  
    487485
    488486// Generate Lanczos interpolation kernels
    489 static void lanczos(double values[],    // Interpolation kernel to generate
     487static void lanczos(float values[],    // Interpolation kernel to generate
    490488                    bool *exact,        // Exact shift?
    491489                    int num,            // Number of values in the kernel
    … …  
    506504    } else {
    507505        *exact = false;
    508         double norm1 = 2.0 / PS_SQR(M_PI); // Normalisation for laczos
    509         double norm2 = M_PI * 4.0 / (float)num; // Normalisation for sinc function 1
    510         double norm3 =M_PI_2 * 4.0 / (float)num; // Normalisation for sinc function 2
    511         double pos = - (num - 1)/2 - frac;  // Position of interest
     506        float norm1 = 2.0 / PS_SQR(M_PI); // Normalisation for laczos
     507        float norm2 = M_PI * 4.0 / (float)num; // Normalisation for sinc function 1
     508        float norm3 = M_PI_2 * 4.0 / (float)num; // Normalisation for sinc function 2
     509        float pos = - (num - 1)/2 - frac;  // Position of interest
    512510        for (int i = 0; i < num; i++, pos += 1.0) {
    513511            if (fabs(pos) < DBL_EPSILON) {
    514512                values[i] = 1.0;
    515513            } else {
    516                 values[i] = norm1 * sin(pos * norm2) * sin(pos * norm3) / PS_SQR(pos);
     514                values[i] = norm1 * sinf(pos * norm2) * sinf(pos * norm3) / PS_SQR(pos);
    517515            }
    518516        }
    … …  
    553551// Generate the interpolation kernels for separable case; they should be normalised to unity
    554552static inline void interpolateSeparateGenerate(int xNum, int yNum, // Size of interpolation kernel
    555                                                double xKernel[xNum], double yKernel[yNum], // Kernels
     553                                               float xKernel[xNum], float yKernel[yNum], // Kernels
    556554                                               bool *xExact, bool *yExact, // Exact shifts?
    557555                                               int xCentral, int yCentral, // Central pixel of convolution
    … …  
    566564      case PS_INTERPOLATE_LANCZOS3:
    567565      case PS_INTERPOLATE_LANCZOS4: {
    568           double xFrac = x - xCentral - 0.5; // Fraction of pixel in x
     566          float xFrac = x - xCentral - 0.5; // Fraction of pixel in x
    569567          lanczos(xKernel, xExact, xNum, xFrac, allowExact);
    570           double yFrac = y - yCentral - 0.5; // Fraction of pixel in y
     568          float yFrac = y - yCentral - 0.5; // Fraction of pixel in y
    571569          lanczos(yKernel, yExact, yNum, yFrac, allowExact);
    572570          break;
    … …  
    598596    INTERPOLATE_BOUNDS();
    599597
    600     double xKernel[xNum], yKernel[yNum]; // Interpolation kernels
     598    float xKernel[xNum], yKernel[yNum]; // Interpolation kernels
    601599    bool xExact, yExact;                // Exact shift?
    602600    interpolateSeparateGenerate(xNum, yNum, xKernel, yKernel, &xExact, &yExact,
    … …  
    605603    float sumKernel = 0.0;              // Sum of the kernel
    606604    double sumKernel2 = 0.0;            // Sum of the kernel-squared
    607     float sumBad = 0.0;                 // Sum of bad kernel-squared contributions
     605    double sumBad = 0.0;                // Sum of bad kernel-squared contributions
    608606    psMaskType maskVal = options->maskVal; // Value to mask
    609607    double sumImage = 0.0;              // Sum of image multiplied by kernel
    … …  
    617615    // Add contributions in an area outside the image
    618616#define INTERPOLATE_SEPARATE_SETUP_OFFIMAGE_COL() \
    619     float xSumBad = 0.0; \
    620     double xSumKernel2 = 0.0;
     617    float xSumBad = 0.0;
     618
    621619#define INTERPOLATE_SEPARATE_ADD_OFFIMAGE_COL() { \
    622         float kernelValue = xKernel[i]; /* Value of kernel */ \
    623         double kernelValue2 = PS_SQR(kernelValue); /* Square of kernel */ \
    624         xSumBad += kernelValue2; \
    625         xSumKernel2 += kernelValue2; \
    626     }
     620        xSumBad += PS_SQR(xKernel[i]); \
     621    }
     622
    627623#define INTERPOLATE_SEPARATE_ADD_OFFIMAGE_ROW() { \
    628         float kernelValue = yKernel[j]; /* Value of kernel */ \
    629         double kernelValue2 = PS_SQR(kernelValue); /* Square of kernel */ \
    630         sumBad += kernelValue2 * xSumBad; \
    631         sumKernel2 += kernelValue2 * xSumKernel2; \
     624        sumBad += PS_SQR(yKernel[j]) * xSumBad; \
    632625    }
    633626
    … …  
    696689                          xSumVariance += kernelValue2 * *varianceData; \
    697690                          xSumKernel += kernelValue; \
     691                          xSumKernel2 += kernelValue2; \
    698692                      } \
    699                       xSumKernel2 += kernelValue2; \
    700693                  } \
    701694                  float kernelValue = yKernel[j]; /* Value of kernel in y */ \
    … …  
    745738                      xSumImage += kernelValue * image->data.TYPE[yPix][xPix]; \
    746739                      xSumKernel += kernelValue; \
     740                      xSumKernel2 += kernelValue2; \
    747741                  } \
    748                   xSumKernel2 += kernelValue2; \
    749742              } \
    750743              float kernelValue = yKernel[j]; /* Value of kernel in y */ \
    … …  
    794787    *imageValue = sumImage / sumKernel;
    795788    if (wantVariance) {
    796         *varianceValue = sumVariance * PS_SQR(sumKernel) / sumKernel2;
     789        *varianceValue = sumVariance / sumKernel2;
    797790    }
    798791    if (sumBad == 0) {
    799792        // Completely good pixel
    800793        status = PS_INTERPOLATE_STATUS_GOOD;
    801     } else if (sumBad / sumKernel2 < PS_SQR(options->poorFrac)) {
     794    } else if (sumBad < PS_SQR(options->poorFrac) * (sumBad + sumKernel2) ) {
    802795        // Some pixels masked: poor pixel
    803796        if (haveMask && maskValue) {
    … …  
    888881    }
    889882
    890     double kernel[yNum][xNum];          // Interpolation kernel for straight interpolation
     883    float kernel[yNum][xNum];          // Interpolation kernel for straight interpolation
    891884    bool xExact, yExact;                // Exact shift?
    892885    // Pretend we're not shifting pixels (i.e., we don't care about any exact shifting) because on the off
    … …  
    896889                              x, y, options->mode, false);
    897890
    898     double sumKernel = 0.0;             // Sum of kernel
     891    float sumKernel = 0.0;             // Sum of kernel
    899892    double sumKernel2 = 0.0;            // Sum of kernel squares
    900893    for (int j = 0; j < yNum; j++) {
    … …  
    925918    }
    926919
    927     double xKernel[xNum], yKernel[yNum]; // Interpolation kernels for separable interpolation
     920    float xKernel[xNum], yKernel[yNum]; // Interpolation kernels for separable interpolation
    928921    bool xExact, yExact;                // Exact shift?
    929922    // Pretend we're not shifting pixels (i.e., we don't care about any exact shifting) because on the off
    … …  
    933926                                xCentral, yCentral, x, y, options->mode, false);
    934927
    935     double ySumKernel = 0.0;            // Sum of kernel in y
     928    float ySumKernel = 0.0;            // Sum of kernel in y
    936929    double ySumKernel2 = 0.0;           // Sum of kernel squared in y
    937930    for (int j = 0; j < yNum; j++) {
    938         double xSumKernel = 0.0;        // Sum of kernel in x
     931        float xSumKernel = 0.0;        // Sum of kernel in x
    939932        double xSumKernel2 = 0.0;       // Sum of kernel squared in x
    940933        for (int i = 0; i < xNum; i++) {
Note: See TracChangeset for help on using the changeset viewer.