IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 13406


Ignore:
Timestamp:
May 16, 2007, 6:13:19 PM (19 years ago)
Author:
Paul Price
Message:

Replacing "middle man" static inline functions with macros, because
the "end man" weight functions (which are declared "static inline")
weren't disappearing in the assembly. This is a slight optimisation.

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r13398 r13406  
    44 *  @author GLG, MHPCC
    55 *
    6  *  @version $Revision: 1.10 $ $Name: not supported by cvs2svn $
    7  *  @date $Date: 2007-05-16 21:13:51 $
     6 *  @version $Revision: 1.11 $ $Name: not supported by cvs2svn $
     7 *  @date $Date: 2007-05-17 04:13:19 $
    88 *
    99 *  Copyright 2004-2007 Institute for Astronomy, University of Hawaii
    … …  
    5353}
    5454
     55// Generate an image of the solved kernel
     56// Meant to replace the following function declaration:
     57// static inline psKernel *solvedKernel(psKernel *kernel, // Kernel, to return
     58//                                      const psVector *solution, // Solution to the least-squares problem
     59//                                      const pmSubtractionKernels *kernels, // Kernel basis functions
     60//                                      const psImage *polyValues, // Spatial polynomial values
     61//                                      double (*weightFunc)(double value) // Function for weighting
     62//     );
     63// because even if the weightFunc can't be "static inline", they still appear in the assembly.
     64// TARGET: Kernel, to return (psKernel*)
     65// SOLUTION: Solution to the least-squares problem (psVector*)
     66// KERNELS: Kernel basis functions (pmSubtractionKernels*)
     67// POLYVALUES: Spatial polynomial values (psImage*)
     68// FUNC: Function for weighting (double (*weightFunc)(double value))
     69#define SOLVED_KERNEL(TARGET, SOLUTION, KERNELS, POLYVALUES, FUNC) { \
     70    int numKernels = (SOLUTION)->n - 1;   /* Number of kernel basis functions */ \
     71    assert((KERNELS)->u->n == numKernels); \
     72    assert((KERNELS)->v->n == numKernels); \
     73    assert((KERNELS)->xOrder->n == numKernels); \
     74    assert((KERNELS)->yOrder->n == numKernels); \
     75    \
     76    /* Ensure the subIndex for POIS kernels is what is expected */ \
     77    assert((KERNELS)->type != PM_SUBTRACTION_KERNEL_POIS || \
     78           ((KERNELS)->u->data.S32[(KERNELS)->subIndex] == 0 && \
     79            (KERNELS)->v->data.S32[(KERNELS)->subIndex] == 0 && \
     80            (KERNELS)->xOrder->data.S32[(KERNELS)->subIndex] == 0 && \
     81            (KERNELS)->yOrder->data.S32[(KERNELS)->subIndex] == 0)); \
     82    \
     83    int size = (KERNELS)->size;           /* Kernel half-size */ \
     84    if (!(TARGET)) { \
     85        (TARGET) = psKernelAlloc(-size, size, -size, size); \
     86    } \
     87    psImageInit((TARGET)->image, 0.0); \
     88    \
     89    for (int i = 0; i < numKernels; i++) { \
     90        int xOrder = (KERNELS)->xOrder->data.S32[i]; /* Polynomial order in x */ \
     91        int yOrder = (KERNELS)->yOrder->data.S32[i]; /* Polynomial order in y */ \
     92        float value = FUNC((POLYVALUES)->data.F64[yOrder][xOrder] * (SOLUTION)->data.F64[i]); \
     93        float subValue = FUNC(-(SOLUTION)->data.F64[i]); /* Value to subtract (add because of minus) */ \
     94        switch ((KERNELS)->type) { \
     95          case PM_SUBTRACTION_KERNEL_POIS: { \
     96              int u = (KERNELS)->u->data.S32[i]; /* Offset in x */ \
     97              int v = (KERNELS)->v->data.S32[i]; /* Offset in y */ \
     98              (TARGET)->kernel[v][u] += value; \
     99              if ((KERNELS)->spatialOrder > 0 && i != (KERNELS)->subIndex) { \
     100                  /* The (0,0) element is subtracted from most kernels to preserve photometric scaling */ \
     101                  (TARGET)->kernel[0][0] += subValue; \
     102              } \
     103              break; \
     104          } \
     105          case PM_SUBTRACTION_KERNEL_ISIS: { \
     106              psKernel *preCalc = (KERNELS)->preCalc->data[i]; /* Precalculated values */ \
     107              psKernel *subKernel = (KERNELS)->preCalc->data[(KERNELS)->subIndex]; /* Kernel to subtract */ \
     108              for (int v = -size; v <= size; v++) { \
     109                  for (int u = -size; u <= size; u++) { \
     110                      (TARGET)->kernel[v][u] += value * FUNC(preCalc->kernel[v][u]); \
     111                      /* Subtract (0,0) kernel from other kernels to preserve photometric scaling */ \
     112                      if ((KERNELS)->spatialOrder > 0 && i != (KERNELS)->subIndex) { \
     113                          (TARGET)->kernel[v][u] += subValue * FUNC(subKernel->kernel[v][u]); \
     114                      } \
     115                  } \
     116              } \
     117              break; \
     118          } \
     119          default: \
     120            psAbort("Should never get here."); \
     121        } \
     122    } \
     123}
     124
     125#if 0
    55126// Generate an image of the solved kernel
    56127static inline psKernel *solvedKernel(psKernel *kernel, // Kernel, to return
    … …  
    117188    return kernel;
    118189}
    119 
     190#endif
     191
     192// Generate the convolved pixel value
     193// Meant to replace the following function declaration:
     194// static inline double convolvePixel(const pmSubtractionKernels *kernels, // Kernel basis functions
     195//                                    int index, // Kernel basis function index
     196//                                    int x, int y, // Pixel around which to convolve
     197//                                    const psImage *image, // Image to convolve
     198//                                    const psImage *polyValues, // Spatial polynomial values
     199//                                    double (*weightFunc)(double value) // Function for weighting
     200//     );
     201// because even if the weightFunc can't be "static inline", they still appear in the assembly.
     202//
     203// TARGET: The value to 'return' (double)
     204// KERNELS: Kernel basis functions (pmSubtractionKernels*)
     205// INDEX: Kernel basis function index (int)
     206// X, Y: Pixel around which to convolve (int)
     207// IMAGE: Image to convolve (psImage*)
     208// POLYVALUES: Spatial polynomial values (psImage*)
     209// FUNC: Function for weighting (double (*weightFunc)(double value))
     210#define CONVOLVE_PIXEL(TARGET, KERNELS, INDEX, X, Y, IMAGE, POLYVALUES, FUNC) { \
     211    int xOrder = (KERNELS)->xOrder->data.S32[(INDEX)]; /* Polynomial order in x */ \
     212    int yOrder = (KERNELS)->yOrder->data.S32[(INDEX)]; /* Polynomial order in y */ \
     213    double polyValue = (POLYVALUES)->data.F64[yOrder][xOrder]; /* Value of spatial polynomial */ \
     214    \
     215    switch ((KERNELS)->type) { \
     216      case PM_SUBTRACTION_KERNEL_POIS: { \
     217          /* Convolution with a delta function is just the value specified by the offset */ \
     218          int u = (KERNELS)->u->data.S32[(INDEX)]; /* Offset in x */ \
     219          int v = (KERNELS)->v->data.S32[(INDEX)]; /* Offset in y */ \
     220          (TARGET) = FUNC(polyValue) * (IMAGE)->data.F32[(Y) + v][(X) + u]; /* Value of convolution */ \
     221          if ((KERNELS)->spatialOrder > 0 && (INDEX) != (KERNELS)->subIndex) { \
     222              /* The (0,0) element is subtracted from most kernels to preserve photometric scaling */ \
     223              (TARGET) += FUNC(-1.0) * (IMAGE)->data.F32[(Y)][(X)]; \
     224          } \
     225          break; \
     226      } \
     227      case PM_SUBTRACTION_KERNEL_ISIS: { \
     228          psKernel *kernel = (KERNELS)->preCalc->data[(INDEX)]; /* The convolution kernel */ \
     229          psKernel *subKernel = (KERNELS)->preCalc->data[(KERNELS)->subIndex]; /* Kernel to subtract */ \
     230          int size = (KERNELS)->size;     /* Kernel half-size */ \
     231          double sum = 0.0;             /* Accumulated sum from convolution */ \
     232          double sub = 0.0;             /* Accumulated sum to subtract */ \
     233          for (int v = -size; v <= size; v++) { \
     234              for (int u = -size; u <= size; u++) { \
     235                  sum += FUNC(kernel->kernel[v][u]) * (IMAGE)->data.F32[(Y) + v][(X) + u]; \
     236                  /* The (0,0) kernel is subtracted from other kernels to preserve photometric scaling */ \
     237                  if ((KERNELS)->spatialOrder > 0 && (INDEX) != (KERNELS)->subIndex) { \
     238                      sub += FUNC(subKernel->kernel[v][u]) * (IMAGE)->data.F32[(Y) + v][(X) + u]; \
     239                  } \
     240              } \
     241          } \
     242          TARGET = FUNC(polyValue) * sum + FUNC(-1.0) * sub; \
     243          break; \
     244      } \
     245      default: \
     246        psAbort("Should never get here."); \
     247    } \
     248}
     249
     250#if 0
    120251// Generate the convolved pixel value
    121252static inline double convolvePixel(const pmSubtractionKernels *kernels, // Kernel basis functions
    … …  
    165296    return NAN;
    166297}
     298#endif
    167299
    168300// Weighting function for use with convolvePixel: no weighting applied, suitable for combining image pixels
    … …  
    269401                    // Generate the convolutions
    270402                    for (int i = 0; i < numKernels; i++) {
     403#if 0
    271404                        convolutions->data.F64[i] = convolvePixel(kernels, i, x, y, reference, polyValues,
    272405                                                                  imageWeighting);
     406#else
     407                        CONVOLVE_PIXEL(convolutions->data.F64[i], kernels, i, x, y, reference, polyValues,
     408                                       imageWeighting);
     409#endif
    273410                    }
    274411
    … …  
    450587                for (int x = footprint; x <= footprint; x++) {
    451588                    for (int k = 0; k < numKernels; k++) {
     589#if 0
    452590                        convolvedStamp->data.F32[y + footprint][x + footprint] +=
    453591                            convolvePixel(kernels, k, xStamp + x, yStamp + y, refImage, polyValues,
    454592                                          imageWeighting);
     593#else
     594                        CONVOLVE_PIXEL(convolvedStamp->data.F32[y + footprint][x + footprint], kernels,
     595                                       k, xStamp + x, yStamp + y, refImage, polyValues, imageWeighting);
     596#endif
    455597                    }
    456598                }
    … …  
    532674
    533675    // The appropriate kernel
     676#if 0
    534677    psKernel *kernel = solvedKernel(NULL, solution, kernels, polyValues, imageWeighting);
     678#else
     679    psKernel *kernel = NULL;
     680    SOLVED_KERNEL(kernel, solution, kernels, polyValues, imageWeighting);
     681#endif
    535682
    536683    psFree(polyValues);
    … …  
    630777                                           2.0 * (float)(i + size - numCols/2.0) / (float)numCols,
    631778                                           2.0 * (float)(j + size - numRows/2.0) / (float)numRows);
     779#if 0
    632780            kernelImage = solvedKernel(kernelImage, solution, kernels, polyValues, imageWeighting);
    633781            if (inWeight) {
    … …  
    635783                                            varianceWeighting);
    636784            }
     785#else
     786            SOLVED_KERNEL(kernelImage, solution, kernels, polyValues, imageWeighting);
     787            if (inWeight) {
     788                SOLVED_KERNEL(kernelWeight, solution, kernels, polyValues, varianceWeighting);
     789            }
     790#endif
    637791
    638792            for (int y = j; y < j + fullSize; y++) {
Note: See TracChangeset for help on using the changeset viewer.