IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 14480


Ignore:
Timestamp:
Aug 13, 2007, 4:14:30 PM (19 years ago)
Author:
Paul Price
Message:

Adding enhancement to extract the pixels for each stamp from the image, so that pmSubtractionCalculateEquation doesn't need to look at the input image (so I can later put in a fake target PSF). The stamps are extracted as 'kernels' so that offsets need not be performed.

Location:
trunk/psModules/src/imcombine
Files:
4 edited

Legend:

Unmodified
Added
Removed

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

    r14456 r14480  
    44 *  @author GLG, MHPCC
    55 *
    6  *  @version $Revision: 1.35 $ $Name: not supported by cvs2svn $
    7  *  @date $Date: 2007-08-09 21:20:35 $
     6 *  @version $Revision: 1.36 $ $Name: not supported by cvs2svn $
     7 *  @date $Date: 2007-08-14 02:14:30 $
    88 *
    99 *  Copyright 2004-2007 Institute for Astronomy, University of Hawaii
    … …  
    217217                                   int index, // Kernel basis function index
    218218                                   int x, int y, // Pixel around which to convolve
    219                                    const psImage *image // Image to convolve
     219                                   const psKernel *image // Image to convolve (a kernel for convenience)
    220220                                   )
    221221{
    … …  
    225225          int u = kernels->u->data.S32[index]; // Offset in x
    226226          int v = kernels->v->data.S32[index]; // Offset in y
    227           float value = image->data.F32[y + v][x + u]; // Value of convolution
     227          float value = image->kernel[y + v][x + u]; // Value of convolution
    228228          if (kernels->spatialOrder > 0 && index != kernels->subIndex) {
    229229              // The (0,0) element is subtracted from most kernels to preserve photometric scaling
    230               value -= image->data.F32[y][x];
     230              value -= image->kernel[y][x];
    231231          }
    232232          return value;
    … …  
    242242          for (int v = vStart; v <= vStop; v++) {
    243243              for (int u = uStart; u <= uStop; u++) {
    244                   sum += image->data.F32[y + v][x + u];
     244                  sum += image->kernel[y + v][x + u];
    245245              }
    246246          }
    … …  
    248248          if (kernels->spatialOrder > 0 && index != kernels->subIndex) {
    249249              // The (0,0) element is subtracted from most kernels to preserve photometric scaling
    250               sum -= image->data.F32[y][x];
     250              sum -= image->kernel[y][x];
    251251          }
    252252          return sum;
    … …  
    260260              for (int v = -size; v <= size; v++) {
    261261                  for (int u = -size; u <= size; u++) {
    262                       sum += kernel->kernel[v][u] * image->data.F32[y + v][x + u];
     262                      sum += kernel->kernel[v][u] * image->kernel[y + v][x + u];
    263263                  }
    264264              }
    … …  
    268268          int u = kernels->u->data.S32[index]; // Offset in x
    269269          int v = kernels->v->data.S32[index]; // Offset in y
    270           float value = image->data.F32[y + v][x + u]; // Value of convolution
     270          float value = image->kernel[y + v][x + u]; // Value of convolution
    271271          // The (0,0) delta function is subtracted from most kernels to preserve photometric scaling
    272272          if (kernels->spatialOrder > 0 && index != kernels->subIndex) {
    273               value -= image->data.F32[y][x];
     273              value -= image->kernel[y][x];
    274274          }
    275275          return value;
    … …  
    281281          for (int v = -size; v <= size; v++) {
    282282              for (int u = -size; u <= size; u++) {
    283                   sum += kernel->kernel[v][u] * image->data.F32[y + v][x + u];
     283                  sum += kernel->kernel[v][u] * image->kernel[y + v][x + u];
    284284                  // Photometric scaling is already built in to the precalculated kernel
    285285              }
    … …  
    289289      case PM_SUBTRACTION_KERNEL_RINGS: {
    290290          if (index == kernels->subIndex) {
    291               return image->data.F32[y][x];
     291              return image->kernel[y][x];
    292292          }
    293293          psArray *preCalc = kernels->preCalc->data[index]; // Precalculated data
    … …  
    299299          for (int j = 0; j < num; j++) {
    300300              int u = uCoords->data.S32[j], v = vCoords->data.S32[j]; // Kernel coordinates
    301               sum += image->data.F32[y + v][x + u] * poly->data.F32[j];
     301              sum += image->kernel[y + v][x + u] * poly->data.F32[j];
    302302          }
    303303          // The (0,0) kernel is subtracted from other kernels to preserve photometric scaling
    304           return sum - num * image->data.F32[y][x];
     304          return sum - num * image->kernel[y][x];
    305305      }
    306306      default:
    … …  
    436436
    437437
    438 bool pmSubtractionCalculateEquation(psArray *stamps, const psImage *reference, const psImage *input,
    439                                     const psImage *weight, const pmSubtractionKernels *kernels, int footprint)
     438bool pmSubtractionCalculateEquation(psArray *stamps, const pmSubtractionKernels *kernels, int footprint)
    440439{
    441440    PS_ASSERT_ARRAY_NON_NULL(stamps, false);
    442     PS_ASSERT_IMAGE_NON_NULL(reference, false);
    443     PS_ASSERT_IMAGE_TYPE(reference, PS_TYPE_F32, false);
    444     PS_ASSERT_IMAGE_NON_NULL(input, false);
    445     PS_ASSERT_IMAGE_TYPE(input, PS_TYPE_F32, false);
    446     PS_ASSERT_IMAGES_SIZE_EQUAL(reference, input, false);
    447441    PS_ASSERT_PTR_NON_NULL(kernels, false);
    448442    PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->v, false);
    449     if (weight) {
    450         PS_ASSERT_IMAGE_NON_NULL(weight, false);
    451         PS_ASSERT_IMAGES_SIZE_EQUAL(weight, input, false);
    452         PS_ASSERT_IMAGE_TYPE(weight, PS_TYPE_F32, false);
    453     }
    454443    PS_ASSERT_INT_NONNEGATIVE(footprint, false);
    455 
    456     // Image size
    457     int numCols = reference->numCols;
    458     int numRows = reference->numRows;
    459444
    460445    int spatialOrder = kernels->spatialOrder; // Maximum order of spatial variation
    … …  
    491476            psVectorInit(stampVector, 0.0);
    492477
    493             float xNorm = 2.0 * (float)(stamp->x - numCols/2.0) / (float)numCols; // Normalised x coord
    494             float yNorm = 2.0 * (float)(stamp->y - numRows/2.0) / (float)numRows; // Normalised y coord
    495             psImage *polyValues = spatialPolyValues(spatialOrder, xNorm, yNorm); // Spatial polynomial terms
    496 
    497             for (int y = stamp->y - footprint; y <= stamp->y + footprint; y++) {
    498                 for (int x = stamp->x - footprint; x <= stamp->x + footprint; x++) {
    499                     float invNoise2 = 1.0 / (weight ? weight->data.F32[y][x] :
    500                                              input->data.F32[y][x]); // Inverse square noise
     478            // Spatial polynomial terms
     479            psImage *polyValues = spatialPolyValues(spatialOrder, stamp->xNorm, stamp->yNorm);
     480
     481            psKernel *reference = stamp->reference; // Reference postage stamp
     482            psKernel *input = stamp->input; // Input postage stamp
     483            psKernel *weight = stamp->weight; // Weight map postage stamp
     484
     485            for (int y = - footprint; y <= footprint; y++) {
     486                for (int x = - footprint; x <= footprint; x++) {
     487                    float invNoise2 = 1.0 / weight->kernel[y][x]; // Inverse square noise
    501488
    502489                    // Generate the convolutions
    … …  
    519506                                convolutions->data.F64[j] * invNoise2;
    520507                        }
    521                         stampVector->data.F64[i] += input->data.F32[y][x] * convolutions->data.F64[i] *
     508                        stampVector->data.F64[i] += input->kernel[y][x] * convolutions->data.F64[i] *
    522509                            invNoise2;
    523510
    … …  
    527514                    // Background only terms
    528515                    stampMatrix->data.F64[bgIndex][bgIndex] += invNoise2;
    529                     stampVector->data.F64[bgIndex] += input->data.F32[y][x] * invNoise2;
     516                    stampVector->data.F64[bgIndex] += input->kernel[y][x] * invNoise2;
    530517                }
    531518            }

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

    r14455 r14480  
    66 * @author GLG, MHPCC
    77 *
    8  * @version $Revision: 1.6 $ $Name: not supported by cvs2svn $
    9  * @date $Date: 2007-08-09 20:25:52 $
     8 * @version $Revision: 1.7 $ $Name: not supported by cvs2svn $
     9 * @date $Date: 2007-08-14 02:14:30 $
    1010 * Copyright 2004-207 Institute for Astronomy, University of Hawaii
    1111 */
    … …  
    4040/// Calculate the least-squares equation to match the image quality
    4141bool pmSubtractionCalculateEquation(psArray *stamps, ///< The stamps for which to calculate the equation,
    42                                     const psImage *reference, ///< Reference image
    43                                     const psImage *input, ///< Input image
    44                                     const psImage *weight, ///< Weight image, or NULL
    4542                                    const pmSubtractionKernels *kernels, ///< Kernel parameters
    4643                                    int footprint ///< Half-size of region over which to calculate equation

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

    r14455 r14480  
    3333    stamp->x = 0;
    3434    stamp->y = 0;
     35    stamp->xNorm = 0.0;
     36    stamp->yNorm = 0.0;
    3537    stamp->matrix = NULL;
    3638    stamp->vector = NULL;
    3739    stamp->status = status;
     40
     41    stamp->reference = NULL;
     42    stamp->input = NULL;
     43    stamp->weight = NULL;
    3844
    3945    return stamp;
    … …  
    137143                stamp->x = xBest;
    138144                stamp->y = yBest;
     145                stamp->xNorm = 2.0 * (float)(xBest - numCols/2.0) / (float)numCols;
     146                stamp->yNorm = 2.0 * (float)(yBest - numRows/2.0) / (float)numRows;
     147
     148                // Reset the postage stamps since we're making a new stamp
     149                psFree(stamp->reference);
     150                psFree(stamp->input);
     151                psFree(stamp->weight);
     152                stamp->reference = stamp->input = stamp->weight = NULL;
     153
    139154                if (fluxBest > threshold) {
    140155                    stamp->status = PM_SUBTRACTION_STAMP_CALCULATE;
    … …  
    154169    return stamps;
    155170}
     171
     172
     173bool pmSubtractionExtractStamps(psArray *stamps, psImage *reference, psImage *input, psImage *weight,
     174                                int footprint, const pmSubtractionKernels *kernels)
     175{
     176    PS_ASSERT_ARRAY_NON_NULL(stamps, false);
     177    PS_ASSERT_IMAGE_NON_NULL(reference, false);
     178    PS_ASSERT_IMAGE_TYPE(reference, PS_TYPE_F32, false);
     179    if (input) {
     180        PS_ASSERT_IMAGE_NON_NULL(input, false);
     181        PS_ASSERT_IMAGES_SIZE_EQUAL(input, reference, false);
     182        PS_ASSERT_IMAGE_TYPE(input, PS_TYPE_F32, false);
     183    }
     184    if (weight) {
     185        PS_ASSERT_IMAGE_NON_NULL(weight, false);
     186        PS_ASSERT_IMAGES_SIZE_EQUAL(weight, reference, false);
     187        PS_ASSERT_IMAGE_TYPE(weight, PS_TYPE_F32, false);
     188    }
     189
     190    int numCols = reference->numCols, numRows = reference->numRows; // Size of images
     191    int size = kernels->size + footprint; // Size of postage stamps
     192
     193    if (!weight) {
     194        // Use the input as a rough approximation to the variance map, and HOPE that it's positive.
     195        weight = input;
     196    }
     197
     198    for (int i = 0; i < stamps->n; i++) {
     199        pmSubtractionStamp *stamp = stamps->data[i]; // Stamp of interest
     200
     201
     202        int x = stamp->x, y = stamp->y; // Stamp coordinates
     203        if (x < size || x > numCols - size || y < size || y > numRows - size) {
     204            psError(PS_ERR_UNKNOWN, false, "Stamp %d (%d,%d) is within the image border.\n", i, x, y);
     205            return false;
     206        }
     207
     208        psRegion region = psRegionSet(x - size, x + size + 1, y - size, y + size + 1); // Region of interest
     209
     210        psImage *refSub = psImageSubset(reference, region); // Subimage with stamp
     211        stamp->reference = psKernelAllocFromImage(refSub, size, size);
     212        psFree(refSub);                 // Drop reference
     213
     214        if (input) {
     215            psImage *inSub = psImageSubset(input, region); // Subimage with stamp
     216            stamp->input = psKernelAllocFromImage(inSub, size, size);
     217            psFree(inSub);              // Drop reference
     218        }
     219
     220        psImage *wtSub = psImageSubset(weight, region); // Subimage with stamp
     221        stamp->weight = psKernelAllocFromImage(wtSub, size, size);
     222        psFree(wtSub);                 // Drop reference
     223
     224    }
     225
     226    return true;
     227}
     228
     229

  • trunk/psModules/src/imcombine/pmSubtractionStamps.h

    r14455 r14480  
    44#include <pslib.h>
    55
     6#include "pmSubtractionKernels.h"
    67
    78/// Status of stamp
    … …  
    1718typedef struct {
    1819    int x, y;                           ///< Position
     20    float xNorm, yNorm;                 ///< Normalised position
     21    psKernel *reference;                ///< Reference image postage stamp
     22    psKernel *input;                    ///< Input image postage stamp
     23    psKernel *weight;                   ///< Weight image postage stamp
    1924    psImage *matrix;                    ///< Associated matrix
    2025    psVector *vector;                   ///< Assoicated vector
    … …  
    3136    );
    3237
     38/// Extract stamps from the images
     39bool pmSubtractionExtractStamps(psArray *stamps, ///< Stamps
     40                                psImage *reference, ///< Reference image
     41                                psImage *input, ///< Input image (or NULL)
     42                                psImage *weight, ///< Weight (variance) map
     43                                int footprint, ///< Stamp footprint size
     44                                const pmSubtractionKernels *kernels ///< Kernel basis functions (for size)
     45    );
    3346
    3447#endif
Note: See TracChangeset for help on using the changeset viewer.