IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Sep 5, 2007, 11:17:59 AM (19 years ago)
Author:
Paul Price
Message:

Generate kernel for the variance map as the square of the kernel for
the image. This gets completely away (finally!) from that horrible
scheme of wrapping whatever was being added in an additional function.

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r14739 r14752  
    44 *  @author GLG, MHPCC
    55 *
    6  *  @version $Revision: 1.56 $ $Name: not supported by cvs2svn $
    7  *  @date $Date: 2007-09-05 00:15:28 $
     6 *  @version $Revision: 1.57 $ $Name: not supported by cvs2svn $
     7 *  @date $Date: 2007-09-05 21:17:59 $
    88 *
    99 *  Copyright 2004-2007 Institute for Astronomy, University of Hawaii
    … …  
    7575}
    7676
    77 // Apply the appropriate weighting to the kernel value
    78 static inline float kernelWeighting(float value, // Value to apply weighting
    79                                     bool variance // Do variance weighting?
    80                                     )
    81 {
    82     return variance ? PS_SQR(value) : value;
    83 }
    84 
     77// Generate the kernel to apply to the variance from the normal kernel
     78static psKernel *varianceKernel(psKernel *normalKernel // Normal kernel
     79                                )
     80{
     81    psKernel *kernel = psKernelAlloc(normalKernel->xMin, normalKernel->xMax,
     82                                     normalKernel->yMin, normalKernel->yMax); // Kernel to return
     83
     84    // Take the square of the normal kernel
     85    for (int v = -size; v <= size; v++) {
     86        for (int u = -size; u <= size; u++) {
     87            kernel->kernel[v][u] = PS_SQR(normalKernel->kernel[v][u]);
     88        }
     89    }
     90
     91    return kernel;
     92}
    8593
    8694// Generate an image of the solved kernel
    … …  
    8896                                     const psVector *solution, // Solution to the least-squares problem
    8997                                     const pmSubtractionKernels *kernels, // Kernel basis functions
    90                                      const psImage *polyValues, // Spatial polynomial values
    91                                      bool varianceWeighting // Do variance weighting?
     98                                     const psImage *polyValues // Spatial polynomial values
    9299                                     )
    93100{
    … …  
    113120
    114121    for (int i = 0; i < numKernels; i++) {
    115         double sum = 0.0;               // Sum of multiple polynomial terms
     122        double value = 0.0;             // Value of this component, from the sum of multiple polynomial terms
    116123        for (int yOrder = 0, index = i; yOrder <= spatialOrder; yOrder++) {
    117124            for (int xOrder = 0; xOrder <= spatialOrder - yOrder; xOrder++, index += numKernels) {
    118                 sum += polyValues->data.F64[yOrder][xOrder] * solution->data.F64[index];
    119             }
    120         }
    121         float value = kernelWeighting(sum, varianceWeighting); // Value to be added to kernel
    122         float subValue = kernelWeighting(- sum, varianceWeighting); // Value to be added for subtraction
     125                value += polyValues->data.F64[yOrder][xOrder] * solution->data.F64[index];
     126            }
     127        }
    123128
    124129        switch (kernels->type) {
    … …  
    129134              if (kernels->spatialOrder > 0 && i != kernels->subIndex) {
    130135                  // The (0,0) element is subtracted from most kernels to preserve photometric scaling
    131                   kernel->kernel[0][0] += subValue;
     136                  kernel->kernel[0][0] -= value;
    132137              }
    133138              break;
    … …  
    142147
    143148              // Normalising sum of kernel component to unity
    144               float norm = kernelWeighting(1.0 / (float)((uStop - uStart + 1) * (vStop - vStart + 1)),
    145                                            varianceWeighting);
     149              float norm = 1.0 / (float)((uStop - uStart + 1) * (vStop - vStart + 1));
    146150
    147151              for (int v = vStart; v <= vStop; v++) {
    … …  
    150154                      if (kernels->spatialOrder > 0 && i != kernels->subIndex) {
    151155                          // The (0,0) element is subtracted from most kernels to preserve photometric scaling
    152                           kernel->kernel[0][0] += subValue;
     156                          kernel->kernel[0][0] -= value;
    153157                      }
    154158                  }
    … …  
    163167                  for (int v = -size; v <= size; v++) {
    164168                      for (int u = -size; u <= size; u++) {
    165                           kernel->kernel[v][u] += kernelWeighting(preCalc->kernel[v][u],
    166                                                                   varianceWeighting) * value;
     169                          kernel->kernel[v][u] += preCalc->kernel[v][u] * value;
    167170                          // Photometric scaling is built into the preCalc kernel --- no subtraction!
    168171                      }
    … …  
    176179                  if (subtract) {
    177180                      // The (0,0) element is subtracted from most kernels to preserve photometric scaling
    178                       kernel->kernel[0][0] += subValue;
     181                      kernel->kernel[0][0] -= value;
    179182                  }
    180183              }
    … …  
    186189              for (int v = -size; v <= size; v++) {
    187190                  for (int u = -size; u <= size; u++) {
    188                       kernel->kernel[v][u] += kernelWeighting(preCalc->kernel[v][u],
    189                                                               varianceWeighting) * value;
     191                      kernel->kernel[v][u] += preCalc->kernel[v][u] * value;
    190192                      // Photometric scaling is built into the preCalc kernel --- no subtraction!
    191193                  }
    … …  
    206208              for (int j = 0; j < num; j++) {
    207209                  int u = uCoords->data.S32[j], v = vCoords->data.S32[j]; // Kernel coordinates
    208                   kernel->kernel[v][u] += kernelWeighting(poly->data.F32[j], varianceWeighting) * value;
     210                  kernel->kernel[v][u] += poly->data.F32[j] * value;
    209211              }
    210212              // Photometric scaling is built into the kernel --- no subtraction!
    … …  
    939941
    940942    // The appropriate kernel
    941     psKernel *kernel = solvedKernel(NULL, solution, kernels, polyValues, false);
     943    psKernel *kernel = solvedKernel(NULL, solution, kernels, polyValues);
    942944
    943945    psFree(polyValues);
    … …  
    10911093                                           2.0 * (float)(j + size + 1 - numRows/2.0) / (float)numRows);
    10921094
    1093             kernelImage = solvedKernel(kernelImage, solution, kernels, polyValues, false);
     1095            kernelImage = solvedKernel(kernelImage, solution, kernels, polyValues);
    10941096            if (inWeight) {
    1095                 kernelWeight = solvedKernel(kernelWeight, solution, kernels, polyValues, true);
     1097                kernelWeight = varianceKernel(kernelImage);
    10961098            }
    10971099
Note: See TracChangeset for help on using the changeset viewer.