IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 14319


Ignore:
Timestamp:
Jul 19, 2007, 11:42:12 AM (19 years ago)
Author:
Paul Price
Message:

CONVOLVE_PIXEL is not needed, since it never needs to do weighting.
Therefore, we can happily use the function version without worrying
about performance.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/psModules/src/imcombine/pmSubtraction.c

    r14318 r14319  
    44 *  @author GLG, MHPCC
    55 *
    6  *  @version $Revision: 1.28 $ $Name: not supported by cvs2svn $
    7  *  @date $Date: 2007-07-19 21:26:12 $
     6 *  @version $Revision: 1.29 $ $Name: not supported by cvs2svn $
     7 *  @date $Date: 2007-07-19 21:42:12 $
    88 *
    99 *  Copyright 2004-2007 Institute for Astronomy, University of Hawaii
    … …  
    274274#endif
    275275
    276 #ifndef USE_FUNCTIONS_INSTEAD_OF_MACROS
    277 // Generate the convolved pixel value
    278 // Meant to replace the following function declaration:
    279 // static inline double convolvePixel(const pmSubtractionKernels *kernels, // Kernel basis functions
    280 //                                    int index, // Kernel basis function index
    281 //                                    int x, int y, // Pixel around which to convolve
    282 //                                    const psImage *image, // Image to convolve
    283 //                                    const psImage *polyValues, // Spatial polynomial values
    284 //                                    double (*weightFunc)(double value) // Function for weighting
    285 //     );
    286 // because even though weightFunc is "static inline", it still appears in the assembly.
    287 //
    288 // TARGET: The value to 'return' (double)
    289 // KERNELS: Kernel basis functions (pmSubtractionKernels*)
    290 // INDEX: Kernel basis function index (int)
    291 // X, Y: Pixel around which to convolve (int)
    292 // IMAGE: Image to convolve (psImage*)
    293 // POLYVALUES: Spatial polynomial values (psImage*)
    294 // FUNC: Function for weighting (double (*weightFunc)(double value))
    295 #define CONVOLVE_PIXEL(TARGET, KERNELS, INDEX, X, Y, IMAGE, POLYVALUES, FUNC) { \
    296     int xOrder = (KERNELS)->xOrder->data.S32[(INDEX)]; /* Polynomial order in x */ \
    297     int yOrder = (KERNELS)->yOrder->data.S32[(INDEX)]; /* Polynomial order in y */ \
    298     double polyValue = (POLYVALUES)->data.F64[yOrder][xOrder]; /* Value of spatial polynomial */ \
    299     \
    300     switch ((KERNELS)->type) { \
    301       case PM_SUBTRACTION_KERNEL_POIS: { \
    302           /* Convolution with a delta function is just the value specified by the offset */ \
    303           int u = (KERNELS)->u->data.S32[(INDEX)]; /* Offset in x */ \
    304           int v = (KERNELS)->v->data.S32[(INDEX)]; /* Offset in y */ \
    305           (TARGET) = FUNC(polyValue) * (IMAGE)->data.F32[(Y) + v][(X) + u]; /* Value of convolution */ \
    306           if ((KERNELS)->spatialOrder > 0 && (INDEX) != (KERNELS)->subIndex) { \
    307               /* The (0,0) element is subtracted from most kernels to preserve photometric scaling */ \
    308               (TARGET) += FUNC(-1.0) * (IMAGE)->data.F32[(Y)][(X)]; \
    309           } \
    310           break; \
    311       } \
    312       /* Method for SPAM and FRIES is the same */ \
    313       case PM_SUBTRACTION_KERNEL_SPAM: \
    314       case PM_SUBTRACTION_KERNEL_FRIES: { \
    315           int uStart = (KERNELS)->u->data.S32[(INDEX)]; \
    316           int uStop = (KERNELS)->uStop->data.S32[(INDEX)]; \
    317           int vStart = (KERNELS)->v->data.S32[(INDEX)]; \
    318           int vStop = (KERNELS)->vStop->data.S32[(INDEX)]; \
    319           double sum = 0.0; \
    320           for (int v = vStart; v <= vStop; v++) { \
    321               for (int u = uStart; u <= uStop; u++) { \
    322                   sum += FUNC(polyValue) * (IMAGE)->data.F32[(Y) + v][(X) + u]; \
    323               } \
    324           } \
    325           if ((KERNELS)->spatialOrder > 0 && (INDEX) != (KERNELS)->subIndex) { \
    326               /* The (0,0) element is subtracted from most kernels to preserve photometric scaling */ \
    327               sum += FUNC(-1.0) * (IMAGE)->data.F32[(Y)][(X)]; \
    328           } \
    329           (TARGET) = sum; \
    330           break; \
    331       } \
    332       case PM_SUBTRACTION_KERNEL_GUNK: { \
    333           double value;                 /* The value to return */ \
    334           if ((INDEX) < (KERNELS)->inner) { \
    335               /* Using pre-calculated function */ \
    336               psKernel *kernel = (KERNELS)->preCalc->data[(INDEX)]; /* The convolution kernel */ \
    337               int size = (KERNELS)->size;     /* Kernel half-size */ \
    338               double sum = 0.0;             /* Accumulated sum from convolution */ \
    339               for (int v = -size; v <= size; v++) { \
    340                   for (int u = -size; u <= size; u++) { \
    341                       sum += FUNC(kernel->kernel[v][u]) * (IMAGE)->data.F32[(Y) + v][(X) + u]; \
    342                   } \
    343               } \
    344               value = FUNC(polyValue) * sum; \
    345           } else { \
    346               /* Using delta function */ \
    347               int u = (KERNELS)->u->data.S32[(INDEX)]; /* Offset in x */ \
    348               int v = (KERNELS)->v->data.S32[(INDEX)]; /* Offset in y */ \
    349               value = FUNC(polyValue) * (IMAGE)->data.F32[(Y) + v][(X) + u]; /* Value of convolution */ \
    350           } \
    351           /* The (0,0) delta function is subtracted from most kernels to preserve photometric scaling */ \
    352           if ((KERNELS)->spatialOrder > 0 && (INDEX) != (KERNELS)->subIndex) { \
    353               value += FUNC(-1.0) * (IMAGE)->data.F32[(Y)][(X)]; \
    354           } \
    355           (TARGET) = value; \
    356       } \
    357       case PM_SUBTRACTION_KERNEL_ISIS: { \
    358           psKernel *kernel = (KERNELS)->preCalc->data[(INDEX)]; /* The convolution kernel */ \
    359           psKernel *subKernel = (KERNELS)->preCalc->data[(KERNELS)->subIndex]; /* Kernel to subtract */ \
    360           int size = (KERNELS)->size;     /* Kernel half-size */ \
    361           double sum = 0.0;             /* Accumulated sum from convolution */ \
    362           double sub = 0.0;             /* Accumulated sum to subtract */ \
    363           for (int v = -size, yLocation = (Y) - size; v <= size; v++, yLocation++) { \
    364               for (int u = -size, xLocation = (X) - size; u <= size; u++, xLocation++) { \
    365                   sum += FUNC(kernel->kernel[v][u]) * (IMAGE)->data.F32[yLocation][xLocation]; \
    366                   /* The (0,0) kernel is subtracted from other kernels to preserve photometric scaling */ \
    367                   if ((KERNELS)->spatialOrder > 0 && (INDEX) != (KERNELS)->subIndex) { \
    368                       sub += FUNC(subKernel->kernel[v][u]) * (IMAGE)->data.F32[yLocation][xLocation]; \
    369                   } \
    370               } \
    371           } \
    372           TARGET = FUNC(polyValue) * sum + FUNC(-1.0) * sub; \
    373           break; \
    374       } \
    375       default: \
    376         psAbort("Should never get here."); \
    377     } \
    378 }
    379 #else
    380276// Generate the convolved pixel value
    381277static inline double convolvePixel(const pmSubtractionKernels *kernels, // Kernel basis functions
    … …  
    383279                                   int x, int y, // Pixel around which to convolve
    384280                                   const psImage *image, // Image to convolve
    385                                    const psImage *polyValues, // Spatial polynomial values
    386                                    double (*weightFunc)(double value) // Function for weighting
    387     )
     281                                   const psImage *polyValues // Spatial polynomial values
     282                                   )
    388283{
    389284    int xOrder = kernels->xOrder->data.S32[index]; // Polynomial order in x
    … …  
    396291          int u = kernels->u->data.S32[index]; // Offset in x
    397292          int v = kernels->v->data.S32[index]; // Offset in y
    398           double value = weightFunc(polyValue) * image->data.F32[y + v][x + u]; // Value of convolution
     293          double value = polyValue * image->data.F32[y + v][x + u]; // Value of convolution
    399294          if (kernels->spatialOrder > 0 && index != kernels->subIndex) {
    400295              // The (0,0) element is subtracted from most kernels to preserve photometric scaling
    401               value += weightFunc(-1.0) * image->data.F32[y][x];
     296              value -= image->data.F32[y][x];
    402297          }
    403298          return value;
    404299      }
    405       /* Method for SPAM and FRIES is the same */
     300        // Method for SPAM and FRIES is the same
    406301      case PM_SUBTRACTION_KERNEL_SPAM:
    407302      case PM_SUBTRACTION_KERNEL_FRIES: {
    … …  
    413308          for (int v = vStart; v <= vStop; v++) {
    414309              for (int u = uStart; u <= uStop; u++) {
    415                   sum += weightFunc(polyValue) * image->data.F32[y + v][x + u];
     310                  sum += polyValue * image->data.F32[y + v][x + u];
    416311              }
    417312          }
    418313          if (kernels->spatialOrder > 0 && index != kernels->subIndex) {
    419               /* The (0,0) element is subtracted from most kernels to preserve photometric scaling */
    420               sum += weightFunc(-1.0) * image->data.F32[y][x];
     314              // The (0,0) element is subtracted from most kernels to preserve photometric scaling
     315              sum -= image->data.F32[y][x];
    421316          }
    422317          return sum;
    … …  
    431326              for (int v = -size; v <= size; v++) {
    432327                  for (int u = -size; u <= size; u++) {
    433                       sum += weightFunc(kernel->kernel[v][u]) * image->data.F32[y + v][x + u];
     328                      sum += kernel->kernel[v][u] * image->data.F32[y + v][x + u];
    434329                  }
    435330              }
    436               value = weightFunc(polyValue) * sum;
     331              value = polyValue * sum;
    437332          } else {
    438333              // Using delta function
    439334              int u = kernels->u->data.S32[index]; // Offset in x
    440335              int v = kernels->v->data.S32[index]; // Offset in y
    441               value = weightFunc(polyValue) * image->data.F32[y + v][x + u]; // Value of convolution
     336              value = polyValue * image->data.F32[y + v][x + u]; // Value of convolution
    442337          }
    443           /* The (0,0) delta function is subtracted from most kernels to preserve photometric scaling */
     338          // The (0,0) delta function is subtracted from most kernels to preserve photometric scaling
    444339          if (kernels->spatialOrder > 0 && index != kernels->subIndex) {
    445               value += weightFunc(-1.0) * image->data.F32[y][x];
     340              value -= image->data.F32[y][x];
    446341          }
    447342          return value;
    … …  
    455350          for (int v = -size; v <= size; v++) {
    456351              for (int u = -size; u <= size; u++) {
    457                   sum += weightFunc(kernel->kernel[v][u]) * image->data.F32[y + v][x + u];
     352                  sum += kernel->kernel[v][u] * image->data.F32[y + v][x + u];
    458353                  // The (0,0) kernel is subtracted from other kernels to preserve photometric scaling
    459354                  if (kernels->spatialOrder > 0 && index != kernels->subIndex) {
    460                       sub += weightFunc(subKernel->kernel[v][u]) * image->data.F32[y + v][x + u];
     355                      sub += subKernel->kernel[v][u] * image->data.F32[y + v][x + u];
    461356                  }
    462357              }
    463358          }
    464           return weightFunc(polyValue) * sum + weightFunc(-1.0) * sub;
     359          return polyValue * sum - sub;
    465360      }
    466361      default:
    … …  
    469364    return NAN;
    470365}
    471 #endif
    472 
    473 // Weighting function for use with convolvePixel: no weighting applied, suitable for combining image pixels
     366
     367// Weighting function for use with solvedKernel: no weighting applied, suitable for combining image pixels
    474368static inline double imageWeighting(double value)
    475369{
    … …  
    477371}
    478372
    479 // Weighting function for use with convolvePixel: weighting suitable for combining variances
     373// Weighting function for use with solvedKernel: weighting suitable for combining variances
    480374static inline double varianceWeighting(double value)
    481375{
    … …  
    674568                for (int x = stamp->x - footprint; x <= stamp->x + footprint; x++) {
    675569                    float invNoise2 = 1.0 / (weight ? weight->data.F32[y][x] :
    676                                              input->data.F32[y][x]); // Inverse square noise
     570                                             reference->data.F32[y][x]); // Inverse square noise
    677571
    678572                    // Generate the convolutions
    679573                    for (int i = 0; i < numKernels; i++) {
    680 #ifdef USE_FUNCTIONS_INSTEAD_OF_MACROS
    681                         convolutions->data.F64[i] = convolvePixel(kernels, i, x, y, reference, polyValues,
    682                                                                   imageWeighting);
    683 #else
    684                         CONVOLVE_PIXEL(convolutions->data.F64[i], kernels, i, x, y, reference, polyValues,
    685                                        imageWeighting);
    686 #endif
     574                        convolutions->data.F64[i] = convolvePixel(kernels, i, x, y, reference, polyValues);
    687575                    }
    688576
Note: See TracChangeset for help on using the changeset viewer.