IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 14456


Ignore:
Timestamp:
Aug 9, 2007, 11:20:35 AM (19 years ago)
Author:
Paul Price
Message:

Changing the system by which the kernel weighting is performed.
Previously, was passing function pointers (the functions were defined
'static inline' to attempt to avoid the function call overhead). This
is annoying for a couple of reasons:

  1. Because we're passing function pointers, rather than calling the function directly, the optimisation isn't performed, and you get the function call overhead. You can get around this by having a macro version of the function to which you're passing the function pointer, but that's a maintenance nightmare.
  1. Using the function version (as opposed to the macro version) produces relocation errors when linking with gcc-4.1.1-r3 (​http://gcc.gnu.org/bugzilla/show_bug.cgi?id=30153).

In light of 1 and 2, I'm doing away with it at last.

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r14455 r14456  
    44 *  @author GLG, MHPCC
    55 *
    6  *  @version $Revision: 1.34 $ $Name: not supported by cvs2svn $
    7  *  @date $Date: 2007-08-09 20:25:52 $
     6 *  @version $Revision: 1.35 $ $Name: not supported by cvs2svn $
     7 *  @date $Date: 2007-08-09 21:20:35 $
    88 *
    99 *  Copyright 2004-2007 Institute for Astronomy, University of Hawaii
    … …  
    2525#include "pmSubtraction.h"
    2626
    27 #define USE_FUNCTIONS_INSTEAD_OF_MACROS          // Can I pass the address of a static inline function?
    28 
    2927//////////////////////////////////////////////////////////////////////////////////////////////////////////////
    3028// Private (file-static) functions
    … …  
    7169}
    7270
    73 #ifndef USE_FUNCTIONS_INSTEAD_OF_MACROS
    74 // Generate an image of the solved kernel
    75 // Meant to replace the following function declaration:
    76 // static inline psKernel *solvedKernel(psKernel *kernel, // Kernel, to return
    77 //                                      const psVector *solution, // Solution to the least-squares problem
    78 //                                      const pmSubtractionKernels *kernels, // Kernel basis functions
    79 //                                      const psImage *polyValues, // Spatial polynomial values
    80 //                                      double (*weightFunc)(double value) // Function for weighting
    81 //     );
    82 // because even though the weightFunc is "static inline", it still appears in the assembly.
    83 // TARGET: Kernel, to return (psKernel*)
    84 // SOLUTION: Solution to the least-squares problem (psVector*)
    85 // KERNELS: Kernel basis functions (pmSubtractionKernels*)
    86 // POLYVALUES: Spatial polynomial values (psImage*)
    87 // FUNC: Function for weighting (double (*weightFunc)(double value))
    88 #define SOLVED_KERNEL(TARGET, SOLUTION, KERNELS, POLYVALUES, FUNC) { \
    89     int numKernels = (SOLUTION)->n - 1;   /* Number of kernel basis functions */ \
    90     assert((KERNELS)->u->n == numKernels); \
    91     assert((KERNELS)->v->n == numKernels); \
    92     assert((KERNELS)->xOrder->n == numKernels); \
    93     assert((KERNELS)->yOrder->n == numKernels); \
    94     \
    95     /* Ensure the subIndex for POIS kernels is what is expected */ \
    96     assert(((KERNELS)->type != PM_SUBTRACTION_KERNEL_POIS && \
    97            (KERNELS)->type != PM_SUBTRACTION_KERNEL_SPAM && \
    98            (KERNELS)->type != PM_SUBTRACTION_KERNEL_FRIES) || \
    99            ((KERNELS)->u->data.S32[(KERNELS)->subIndex] == 0 && \
    100             (KERNELS)->v->data.S32[(KERNELS)->subIndex] == 0 && \
    101             (KERNELS)->xOrder->data.S32[(KERNELS)->subIndex] == 0 && \
    102             (KERNELS)->yOrder->data.S32[(KERNELS)->subIndex] == 0)); \
    103     \
    104     int size = (KERNELS)->size;           /* Kernel half-size */ \
    105     if (!(TARGET)) { \
    106         (TARGET) = psKernelAlloc(-size, size, -size, size); \
    107     } \
    108     psImageInit((TARGET)->image, 0.0); \
    109     \
    110     for (int i = 0; i < numKernels; i++) { \
    111         int xOrder = (KERNELS)->xOrder->data.S32[i]; /* Polynomial order in x */ \
    112         int yOrder = (KERNELS)->yOrder->data.S32[i]; /* Polynomial order in y */ \
    113         float value = FUNC((POLYVALUES)->data.F64[yOrder][xOrder] * (SOLUTION)->data.F64[i]); \
    114         float subValue = FUNC(-(SOLUTION)->data.F64[i]); /* Value to subtract (add because of minus) */ \
    115         switch ((KERNELS)->type) { \
    116           case PM_SUBTRACTION_KERNEL_POIS: { \
    117               int u = (KERNELS)->u->data.S32[i]; /* Offset in x */ \
    118               int v = (KERNELS)->v->data.S32[i]; /* Offset in y */ \
    119               (TARGET)->kernel[v][u] += value; \
    120               if ((KERNELS)->spatialOrder > 0 && i != (KERNELS)->subIndex) { \
    121                   /* The (0,0) element is subtracted from most kernels to preserve photometric scaling */ \
    122                   (TARGET)->kernel[0][0] += subValue; \
    123               } \
    124               break; \
    125           } \
    126           /* SPAM and FRIES use the same method */ \
    127           case PM_SUBTRACTION_KERNEL_SPAM: \
    128           case PM_SUBTRACTION_KERNEL_FRIES: { \
    129             int uStart = (KERNELS)->u->data.S32[i]; \
    130             int uStop = (KERNELS)->uStop->data.S32[i]; \
    131             int vStart = (KERNELS)->v->data.S32[i]; \
    132             int vStop = (KERNELS)->vStop->data.S32[i]; \
    133             /* Normalising sum of kernel component to unity */ \
    134             value /= FUNC((uStop - uStart) * (vStop - vStart)); \
    135             for (int v = vStart; v <= vStop; v++) { \
    136                 for (int u = uStart; u <= uStop; u++) { \
    137                     (TARGET)->kernel[v][u] += value; \
    138                 } \
    139             } \
    140             if ((KERNELS)->spatialOrder > 0 && i != (KERNELS)->subIndex) { \
    141                 /* The (0,0) element is subtracted from most kernels to preserve photometric scaling */ \
    142                 (TARGET)->kernel[0][0] += subValue; \
    143             } \
    144             break; \
    145           } \
    146           case PM_SUBTRACTION_KERNEL_GUNK: { \
    147               if (i < (KERNELS)->inner) { \
    148                   /* Using pre-calculated function */ \
    149                   psKernel *preCalc = (KERNELS)->preCalc->data[i]; /* Precalculated values */ \
    150                   for (int v = -size; v <= size; v++) { \
    151                       for (int u = -size; u <= size; u++) { \
    152                           (TARGET)->kernel[v][u] += value * FUNC(preCalc->kernel[v][u]); \
    153                       } \
    154                   } \
    155               } else { \
    156                   /* Using delta function */ \
    157                   int u = (KERNELS)->u->data.S32[i]; /* Offset in x */ \
    158                   int v = (KERNELS)->v->data.S32[i]; /* Offset in y */ \
    159                   (TARGET)->kernel[v][u] += value; \
    160               } \
    161               /* The (0,0) kernel is subtracted from other kernels to preserve photometric scaling */ \
    162               if ((KERNELS)->spatialOrder > 0 && i != (KERNELS)->subIndex) { \
    163                   (TARGET)->kernel[0][0] += subValue; \
    164               } \
    165               break; \
    166           } \
    167           case PM_SUBTRACTION_KERNEL_ISIS: { \
    168               psKernel *preCalc = (KERNELS)->preCalc->data[i]; /* Precalculated values */ \
    169               for (int v = -size; v <= size; v++) { \
    170                   for (int u = -size; u <= size; u++) { \
    171                       (TARGET)->kernel[v][u] += value * FUNC(preCalc->kernel[v][u]); \
    172                       /* Photometric scaling is already built in to the precalculated kernel */ \
    173                   } \
    174               } \
    175               break; \
    176           } \
    177           case PM_SUBTRACTION_KERNEL_RINGS: { \
    178               if (i == (KERNELS)->subIndex) { \
    179                   (TARGET)->kernel[0][0] += value; \
    180                   break; \
    181               } \
    182               psArray *preCalc = (KERNELS)->preCalc->data[i]; /* Precalculated data */ \
    183               psVector *uCoords = preCalc->data[0]; /* u coordinates */ \
    184               psVector *vCoords = preCalc->data[1]; /* v coordinates */ \
    185               psVector *poly = preCalc->data[2]; /* Polynomial values */ \
    186               int num = uCoords->n;     /* Number of pixels */ \
    187               value /= weightFunc(num); /* Normalising sum of kernel component to unity */ \
    188               for (int j = 0; j < num; j++) { \
    189                   int u = uCoords->data.S32[j], v = vCoords->data.S32[j]; /* Kernel coordinates */ \
    190                   (TARGET)->kernel[v][u] += value * FUNC(poly->data.F32[j]); \
    191               } \
    192               /* The (0,0) kernel is subtracted from other kernels to preserve photometric scaling */ \
    193               if (kernels->spatialOrder > 0) { \
    194                   kernel->kernel[0][0] += subValue; \
    195               } \
    196               break; \
    197           } \
    198           default: \
    199             psAbort("Should never get here."); \
    200         } \
    201     } \
    202 }
    203 #else
     71// Apply the appropriate weighting to the kernel value
     72static inline float kernelWeighting(float value, // Value to apply weighting
     73                                    bool variance // Do variance weighting?
     74                                    )
     75{
     76    return variance ? PS_SQR(value) : value;
     77}
     78
     79
    20480// Generate an image of the solved kernel
    20581static inline psKernel *solvedKernel(psKernel *kernel, // Kernel, to return
    … …  
    20783                                     const pmSubtractionKernels *kernels, // Kernel basis functions
    20884                                     const psImage *polyValues, // Spatial polynomial values
    209                                      double (*weightFunc)(double value) // Function for weighting
    210     )
     85                                     bool varianceWeighting // Do variance weighting?
     86                                     )
    21187{
    21288    int numKernels = kernels->num;      // Number of kernel basis functions
    … …  
    231107
    232108    for (int i = 0; i < numKernels; i++) {
    233         double value = 0.0;              // Value of kernel coefficient
     109        double sum = 0.0;               // Sum of multiple polynomial terms
    234110        for (int yOrder = 0, index = i; yOrder <= spatialOrder; yOrder++) {
    235111            for (int xOrder = 0; xOrder <= spatialOrder - yOrder; xOrder++, index += numKernels) {
    236                 value += polyValues->data.F64[yOrder][xOrder] * solution->data.F64[index];
    237             }
    238         }
     112                sum += polyValues->data.F64[yOrder][xOrder] * solution->data.F64[index];
     113            }
     114        }
     115        float value = kernelWeighting(sum, varianceWeighting); // Value to be added to kernel
     116        float subValue = kernelWeighting(- sum, varianceWeighting); // Value to be added for subtraction
    239117
    240118        switch (kernels->type) {
    … …  
    242120              int u = kernels->u->data.S32[i]; // Offset in x
    243121              int v = kernels->v->data.S32[i]; // Offset in y
    244               kernel->kernel[v][u] += weightFunc(value);
     122              kernel->kernel[v][u] += value;
    245123              if (kernels->spatialOrder > 0 && i != kernels->subIndex) {
    246124                  // The (0,0) element is subtracted from most kernels to preserve photometric scaling
    247                   kernel->kernel[0][0] += weightFunc(- value);
     125                  kernel->kernel[0][0] += subValue;
    248126              }
    249127              break;
    … …  
    258136
    259137              // Normalising sum of kernel component to unity
    260               float norm = 1.0 / weightFunc((uStop - uStart) * (vStop - vStart));
     138              float norm = kernelWeighting(1.0 / (float)((uStop - uStart) * (vStop - vStart)),
     139                                           varianceWeighting);
    261140
    262141              for (int v = vStart; v <= vStop; v++) {
    263142                  for (int u = uStart; u <= uStop; u++) {
    264                       kernel->kernel[v][u] += norm * weightFunc(value);
     143                      kernel->kernel[v][u] += norm * value;
    265144                      if (kernels->spatialOrder > 0 && i != kernels->subIndex) {
    266145                          // The (0,0) element is subtracted from most kernels to preserve photometric scaling
    267                           kernel->kernel[0][0] += weightFunc(- value);
     146                          kernel->kernel[0][0] += subValue;
    268147                      }
    269148                  }
    … …  
    278157                  for (int v = -size; v <= size; v++) {
    279158                      for (int u = -size; u <= size; u++) {
    280                           kernel->kernel[v][u] += weightFunc(preCalc->kernel[v][u] * value);
     159                          kernel->kernel[v][u] += kernelWeighting(preCalc->kernel[v][u],
     160                                                                  varianceWeighting) * value;
    281161                          // Photometric scaling is built into the preCalc kernel --- no subtraction!
    282162                      }
    … …  
    287167                  int u = kernels->u->data.S32[i]; // Offset in x
    288168                  int v = kernels->v->data.S32[i]; // Offset in y
    289                   kernel->kernel[v][u] += weightFunc(value);
     169                  kernel->kernel[v][u] += value;
    290170                  if (subtract) {
    291171                      // The (0,0) element is subtracted from most kernels to preserve photometric scaling
    292                       kernel->kernel[0][0] += weightFunc(- value);
     172                      kernel->kernel[0][0] += subValue;
    293173                  }
    294174              }
    … …  
    300180              for (int v = -size; v <= size; v++) {
    301181                  for (int u = -size; u <= size; u++) {
    302                       kernel->kernel[v][u] += weightFunc(preCalc->kernel[v][u] * value);
     182                      kernel->kernel[v][u] += kernelWeighting(preCalc->kernel[v][u],
     183                                                              varianceWeighting) * value;
    303184                      // Photometric scaling is built into the preCalc kernel --- no subtraction!
    304185                  }
    … …  
    308189          case PM_SUBTRACTION_KERNEL_RINGS: {
    309190              if (i == kernels->subIndex) {
    310                   kernel->kernel[0][0] += weightFunc(value);
     191                  kernel->kernel[0][0] += value;
    311192                  break;
    312193              }
    … …  
    319200              for (int j = 0; j < num; j++) {
    320201                  int u = uCoords->data.S32[j], v = vCoords->data.S32[j]; // Kernel coordinates
    321                   kernel->kernel[v][u] += weightFunc(value * poly->data.F32[j]);
     202                  kernel->kernel[v][u] += kernelWeighting(poly->data.F32[j], varianceWeighting) * value;
    322203              }
    323               kernel->kernel[0][0] += weightFunc(- value * num);
     204              kernel->kernel[0][0] += kernelWeighting(num, varianceWeighting) * subValue;
    324205              break;
    325206          }
    … …  
    331212    return kernel;
    332213}
    333 #endif
    334214
    335215// Generate the convolved pixel value
    … …  
    429309    return NAN;
    430310}
    431 
    432 // Weighting function for use with solvedKernel: no weighting applied, suitable for combining image pixels
    433 static inline double imageWeighting(double value)
    434 {
    435     return value;
    436 }
    437 
    438 // Weighting function for use with solvedKernel: weighting suitable for combining variances
    439 static inline double varianceWeighting(double value)
    440 {
    441     return PS_SQR(value);
    442 }
    443 
    444311
    445312// Mark a pixel as blank in the image, mask and weight
    … …  
    830697            psImage *polyValues = spatialPolyValues(kernels->spatialOrder, xNorm, yNorm); // Polynomial terms
    831698
    832 #ifdef USE_FUNCTIONS_INSTEAD_OF_MACROS
    833             kernelImage = solvedKernel(kernelImage, solution, kernels, polyValues, imageWeighting);
    834 #else
    835             SOLVED_KERNEL(kernelImage, solution, kernels, polyValues, imageWeighting);
    836 #endif
     699            kernelImage = solvedKernel(kernelImage, solution, kernels, polyValues, false);
    837700            for (int y = yStamp - footprint, j = 0; y <= yStamp + footprint; y++, j++) {
    838701                for (int x = xStamp - footprint, i = 0; x <= xStamp + footprint; x++, i++) {
    … …  
    923786
    924787    // The appropriate kernel
    925 #ifdef USE_FUNCTIONS_INSTEAD_OF_MACROS
    926     psKernel *kernel = solvedKernel(NULL, solution, kernels, polyValues, imageWeighting);
    927 #else
    928     psKernel *kernel = NULL;
    929     SOLVED_KERNEL(kernel, solution, kernels, polyValues, imageWeighting);
    930 #endif
     788    psKernel *kernel = solvedKernel(NULL, solution, kernels, polyValues, false);
    931789
    932790    psFree(polyValues);
    … …  
    1077935                                           2.0 * (float)(i + size + 1 - numCols/2.0) / (float)numCols,
    1078936                                           2.0 * (float)(j + size + 1 - numRows/2.0) / (float)numRows);
    1079 #ifdef USE_FUNCTIONS_INSTEAD_OF_MACROS
    1080             kernelImage = solvedKernel(kernelImage, solution, kernels, polyValues, imageWeighting);
     937            kernelImage = solvedKernel(kernelImage, solution, kernels, polyValues, false);
    1081938            if (inWeight) {
    1082                 kernelWeight = solvedKernel(kernelWeight, solution, kernels, polyValues,
    1083                                             varianceWeighting);
    1084             }
    1085 #else
    1086             SOLVED_KERNEL(kernelImage, solution, kernels, polyValues, imageWeighting);
    1087             if (inWeight) {
    1088                 SOLVED_KERNEL(kernelWeight, solution, kernels, polyValues, varianceWeighting);
    1089             }
    1090 #endif
     939                kernelWeight = solvedKernel(kernelWeight, solution, kernels, polyValues, true);
     940            }
    1091941
    1092942            for (int y = j; y < PS_MIN(j + fullSize, yMax); y++) {
Note: See TracChangeset for help on using the changeset viewer.