IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Apr 9, 2007, 5:59:13 PM (19 years ago)
Author:
Paul Price
Message:

Fixing lanczos and gauss interpolations.

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r12756 r12776  
    77 *  @author Paul Price, IfA
    88 *
    9  *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
    10  *  @date $Date: 2007-04-05 23:57:42 $
     9 *  @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
     10 *  @date $Date: 2007-04-10 03:59:13 $
    1111 *
    1212 *  Copyright 2004-2007 Institute for Astronomy, University of Hawaii
    … …  
    132132
    133133    int xNum, yNum;                     // Number of interpolation kernel pixels
     134    int xCentral, yCentral;             // Central pixel of the convolution
    134135    switch (options->mode) {
    135136      case PS_INTERPOLATE_BILINEAR:
    136137        xNum = yNum = 2;
     138        // Central pixel is the pixel below the point of interest
     139        xCentral = floor(x - 0.5 + FLT_EPSILON);
     140        yCentral = floor(y - 0.5 + FLT_EPSILON);
    137141        break;
    138142      case PS_INTERPOLATE_BICUBE:
    139143      case PS_INTERPOLATE_GAUSS:
    140144        xNum = yNum = 3;
     145        // Central pixel is the closest pixel to the point of interest
     146        xCentral = x;
     147        yCentral = y;
    141148        break;
    142149      case PS_INTERPOLATE_FLAT:
    … …  
    149156
    150157    const psImage *image = options->image; // Image of interest
    151     int xFloor = floor(x - 0.5 + FLT_EPSILON); // Pixel below point of interest in x
    152     int yFloor = floor(y - 0.5 + FLT_EPSILON); // Pixel below point of interest in y
    153158    int xLast = image->numCols - 1;     // Last pixel in x
    154159    int yLast = image->numRows - 1;     // Last pixel in y
    155160
    156     if (xFloor - (xNum - 1) / 2 < 0 || xFloor + xNum / 2 > xLast ||
    157         yFloor - (yNum - 1) / 2 < 0 || yFloor + yNum / 2 > yLast) {
     161    if (xCentral - (xNum - 1) / 2 < 0 || xCentral + xNum / 2 > xLast ||
     162        yCentral - (yNum - 1) / 2 < 0 || yCentral + yNum / 2 > yLast) {
    158163        // At least one pixel of the interpolation kernel is off the image
    159164        if (imageValue) {
    … …  
    171176    switch (options->mode) {
    172177      case PS_INTERPOLATE_BILINEAR: {
    173           double xFrac = x - 0.5 - xFloor; // Fraction of pixel in x
    174           double yFrac = y - 0.5 - yFloor; // Fraction of pixel in y
     178          double xFrac = x - 0.5 - xCentral; // Fraction of pixel in x
     179          double yFrac = y - 0.5 - yCentral; // Fraction of pixel in y
    175180          kernel[0][0] = (1.0 - xFrac) * (1.0 - yFrac);
    176181          kernel[0][1] = xFrac * (1.0 - yFrac);
    … …  
    180185      }
    181186      case PS_INTERPOLATE_BICUBE: {
    182           double xFrac = x - 0.5 - xFloor; // Fraction of pixel in x
    183           double yFrac = y - 0.5 - yFloor; // Fraction of pixel in y
     187          double xFrac = x - 0.5 - xCentral; // Fraction of pixel in x
     188          double yFrac = y - 0.5 - yCentral; // Fraction of pixel in y
    184189          // Calculation variables
    185190          double xxFrac = xFrac * xFrac / 6.0;
    … …  
    200205      }
    201206      case PS_INTERPOLATE_GAUSS: {
    202           double xFrac = x - 0.5 - xFloor; // Fraction of pixel in x
    203           double yFrac = y - 0.5 - yFloor; // Fraction of pixel in y
    204           double xGaussFrac = 2.0 * erf((double)xNum / 4.0) - 1.0; // Fraction of Gaussian in x
    205           double yGaussFrac = 2.0 * erf((double)yNum / 4.0) - 1.0; // Fraction of Gaussian in x
    206           double norm = 1.0 / (double)xNum / (double)yNum *
    207               (1.0 - xGaussFrac) / xGaussFrac * (1.0 - yGaussFrac) / yGaussFrac; // Normalisation
     207          double xFrac = x - xCentral - 0.5; // Fraction of pixel in x
     208          double yFrac = y - yCentral - 0.5; // Fraction of pixel in y
     209          double sigma = 0.5;           // Gaussian sigma
     210          double norm = 0.0;            // Normalisation
    208211          for (int j = 0, yPos = - (yNum - 1) / 2; j < yNum; j++, yPos++) {
    209212              for (int i = 0, xPos = - (xNum - 1) / 2; i < xNum; i++, xPos++) {
    210                   kernel[j][i] = norm * exp(-0.5 * (PS_SQR(xPos - xFrac) + PS_SQR(yPos - yFrac)));
     213                  norm += kernel[j][i] = exp(-0.5 / PS_SQR(sigma) * (PS_SQR(xPos - xFrac) +
     214                                                                     PS_SQR(yPos - yFrac)));
     215              }
     216          }
     217          norm = 1.0 / norm;
     218          for (int j = 0; j < yNum; j++) {
     219              for (int i = 0; i < xNum; i++) {
     220                  kernel[j][i] *= norm;
    211221              }
    212222          }
    … …  
    224234    #define KERNEL_IMAGE_CASE(TYPE) \
    225235        case PS_TYPE_##TYPE: { \
    226             for (int j = 0, yPix = yFloor - (yNum - 1) / 2; j < yNum; j++, yPix++) { \
    227                 for (int i = 0, xPix = xFloor - (xNum - 1) / 2; i < xNum; i++, xPix++) { \
     236            for (int j = 0, yPix = yCentral - (yNum - 1) / 2; j < yNum; j++, yPix++) { \
     237                for (int i = 0, xPix = xCentral - (xNum - 1) / 2; i < xNum; i++, xPix++) { \
    228238                    value += values[j][i] = kernel[j][i] * image->data.TYPE[yPix][xPix]; \
    229239                } \
    … …  
    262272        int badPix = 0;                 // Number of bad pixels
    263273        *maskValue = 0;
    264         for (int j = 0, yPix = yFloor - (yNum - 1) / 2; j < yNum; j++, yPix++) {
    265             for (int i = 0, xPix = xFloor - (xNum - 1) / 2; i < xNum; i++, xPix++) {
     274        for (int j = 0, yPix = yCentral - (yNum - 1) / 2; j < yNum; j++, yPix++) {
     275            for (int i = 0, xPix = xCentral - (xNum - 1) / 2; i < xNum; i++, xPix++) {
    266276                if (mask->data.PS_TYPE_MASK_DATA[yPix][xPix] & maskVal) {
    267277                    badValue += values[j][i];
    … …  
    289299        #define KERNEL_VARIANCE_CASE(TYPE) \
    290300            case PS_TYPE_##TYPE: { \
    291                 for (int j = 0, yPix = yFloor - (yNum - 1) / 2; j < yNum; j++, yPix++) { \
    292                     for (int i = 0, xPix = xFloor - (xNum - 1) / 2; i < xNum; i++, xPix++) { \
     301                for (int j = 0, yPix = yCentral - (yNum - 1) / 2; j < yNum; j++, yPix++) { \
     302                    for (int i = 0, xPix = xCentral - (xNum - 1) / 2; i < xNum; i++, xPix++) { \
    293303                        *varianceValue += PS_SQR(kernel[j][i]) * variance->data.TYPE[yPix][xPix]; \
    294304                    } \
    … …  
    337347        }
    338348    } else {
     349        double norm3 = 0.0;             // Normalisation
    339350        double norm1 = 2.0 / PS_SQR(M_PI); // Normalisation for laczos
    340         double norm2 = 2.0 / (num / 2);  // Normalisation for sinc functions
    341         double pos = frac - (num - 1)/2;  // Position of interest
     351        double norm2 = 4.0 / (float)num; // 2.0 / (num / 2);  // Normalisation for sinc functions
     352        double pos = - (num - 1)/2 - frac;  // Position of interest
    342353        for (int i = 0; i < num; i++, pos += 1.0) {
    343             values[i] = norm1 * sin(M_PI * pos * norm2) * sin(M_PI_2 * pos * norm2) / PS_SQR(pos);
     354            norm3 += values[i] = norm1 * sin(M_PI * pos * norm2) * sin(M_PI_2 * pos * norm2) / PS_SQR(pos);
     355        }
     356        norm3 = 1.0 / norm3;
     357        for (int i = 0; i < num; i++, pos += 1.0) {
     358            values[i] *= norm3;
    344359        }
    345360    }
    … …  
    373388    }
    374389
     390    // Central pixel is the pixel below the point of interest
     391    int xCentral = floor(x - 0.5), yCentral = floor(y - 0.5); // Central pixel of the convolution
    375392    const psImage *image = options->image; // Image of interest
    376     int xFloor = floor(x - 0.5 + FLT_EPSILON); // Pixel below point of interest in x
    377     int yFloor = floor(y - 0.5 + FLT_EPSILON); // Pixel below point of interest in y
    378393    int xLast = image->numCols - 1;     // Last pixel in x
    379394    int yLast = image->numRows - 1;     // Last pixel in y
    380395
    381     if (xFloor - (xNum - 1) / 2 < 0 || xFloor + xNum / 2 > xLast ||
    382         yFloor - (yNum - 1) / 2 < 0 || yFloor + yNum / 2 > yLast) {
     396    if (xCentral - (xNum - 1) / 2 < 0 || xCentral + xNum / 2 > xLast ||
     397        yCentral - (yNum - 1) / 2 < 0 || yCentral + yNum / 2 > yLast) {
    383398        // At least one pixel of the interpolation kernel is off the image
    384399        if (imageValue) {
    … …  
    400415      case PS_INTERPOLATE_LANCZOS3:
    401416      case PS_INTERPOLATE_LANCZOS4: {
    402           double xFrac = x - 0.5 - xFloor; // Fraction of pixel in x
     417          double xFrac = x - xCentral - 0.5; // Fraction of pixel in x
    403418#if 0
    404419          if (fabs(xFrac) < DBL_EPSILON) {
    … …  
    411426          }
    412427#endif
    413           double yFrac = y - 0.5 - yFloor; // Fraction of pixel in y
     428          double yFrac = y - yCentral - 0.5; // Fraction of pixel in y
    414429#if 0
    415430          if (fabs(yFrac) < DBL_EPSILON) {
    … …  
    435450    #define SEPARATE_IMAGE_CASE(TYPE) \
    436451      case PS_TYPE_##TYPE: { \
    437         for (int j = 0, yPix = yFloor - (yNum - 1) / 2; j < yNum; j++, yPix++) { \
     452        for (int j = 0, yPix = yCentral - (yNum - 1) / 2; j < yNum; j++, yPix++) { \
    438453            double xInterpValue = 0.0; /* Interpolation in x */ \
    439             for (int i = 0, xPix = xFloor - (xNum - 1) / 2; i < xNum; i++, xPix++) { \
     454            for (int i = 0, xPix = xCentral - (xNum - 1) / 2; i < xNum; i++, xPix++) { \
    440455                xInterpValue += values[j][i] = xKernel[i] * image->data.TYPE[yPix][xPix]; \
    441456            } \
    … …  
    475490        int badPix = 0;                 // Number of bad pixels
    476491        *maskValue = 0;
    477         for (int j = 0, yPix = yFloor - (yNum - 1) / 2; j < yNum; j++, yPix++) {
     492        for (int j = 0, yPix = yCentral - (yNum - 1) / 2; j < yNum; j++, yPix++) {
    478493            // Interpolation in x
    479494            double xInterpValue = 0.0;
    480             for (int i = 0, xPix = xFloor - (xNum - 1) / 2; i < xNum; i++, xPix++) {
     495            for (int i = 0, xPix = xCentral - (xNum - 1) / 2; i < xNum; i++, xPix++) {
    481496                if (mask->data.PS_TYPE_MASK_DATA[yPix][xPix] & maskVal) {
    482497                    xInterpValue += values[j][i];
    … …  
    506521        #define SEPARATE_VARIANCE_CASE(TYPE) \
    507522          case PS_TYPE_##TYPE: { \
    508             for (int j = 0, yPix = yFloor - (yNum - 1) / 2; j < yNum; j++, yPix++) { \
     523            for (int j = 0, yPix = yCentral - (yNum - 1) / 2; j < yNum; j++, yPix++) { \
    509524                double xInterpValue = 0.0; /* Interpolation in x */ \
    510                 for (int i = 0, xPix = xFloor - (xNum - 1) / 2; i < xNum; i++, xPix++) { \
     525                for (int i = 0, xPix = xCentral - (xNum - 1) / 2; i < xNum; i++, xPix++) { \
    511526                    xInterpValue += PS_SQR(xKernel[i]) * variance->data.TYPE[yPix][xPix]; \
    512527                } \
Note: See TracChangeset for help on using the changeset viewer.