IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 14513


Ignore:
Timestamp:
Aug 15, 2007, 1:17:56 PM (19 years ago)
Author:
Paul Price
Message:

Adding function to set stamps based on a list of x,y positions (as opposed to find a set of stamps automatically). Adding another function to generate target stamps so that we can convolve to the target PSF.

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

Legend:

Unmodified
Added
Removed
  • trunk/psModules/src/imcombine/pmSubtraction.c

    r14480 r14513  
    44 *  @author GLG, MHPCC
    55 *
    6  *  @version $Revision: 1.36 $ $Name: not supported by cvs2svn $
    7  *  @date $Date: 2007-08-14 02:14:30 $
     6 *  @version $Revision: 1.37 $ $Name: not supported by cvs2svn $
     7 *  @date $Date: 2007-08-15 23:17:56 $
    88 *
    99 *  Copyright 2004-2007 Institute for Astronomy, University of Hawaii
     
    542542                stamp->status = PM_SUBTRACTION_STAMP_REJECTED;
    543543                psTrace("psModules.imcombine", 3, "Rejecting stamp %d (%d,%d) because of bad equation\n",
    544                         i, stamp->x, stamp->y);
     544                        i, (int)(stamp->x + 0.5), (int)(stamp->y + 0.5));
    545545            } else {
    546546                stamp->status = PM_SUBTRACTION_STAMP_USED;
     
    718718            deviations->data.F32[i] = sqrt(stats->sampleMean / 2.0);
    719719            psTrace("psModules.imcombine", 1, "Deviation for stamp %d (%d,%d): %f\n",
    720                     i, stamp->x, stamp->y, deviations->data.F32[i]);
     720                    i, (int)(stamp->x + 0.5), (int)(stamp->y + 0.5), deviations->data.F32[i]);
    721721            totalSquareDev += PS_SQR(deviations->data.F32[i]);
    722722            numStamps++;
     
    736736        if (stamp->status == PM_SUBTRACTION_STAMP_USED && deviations->data.F32[i] > limit) {
    737737            // Mask out the stamp in the image so you it's not found again
    738             psTrace("psModules.imcombine", 3, "Rejecting stamp %d (%d,%d)\n", i, stamp->x, stamp->y);
     738            psTrace("psModules.imcombine", 3, "Rejecting stamp %d (%d,%d)\n", i,
     739                    (int)(stamp->x + 0.5), (int)(stamp->y + 0.5));
    739740            numRejected++;
    740741            for (int y = stamp->y - footprint; y <= stamp->y + footprint; y++) {
  • trunk/psModules/src/imcombine/pmSubtractionStamps.c

    r14480 r14513  
    99#include "pmSubtractionStamps.h"
    1010
     11
     12#define PM_SUBTRACTION_MASK_BAD_STAMP (PM_SUBTRACTION_MASK_BORDER | PM_SUBTRACTION_MASK_FOOTPRINT | \
     13                                       PM_SUBTRACTION_MASK_REJ) // Mask indicates a stamp would be bad
     14
     15
    1116//////////////////////////////////////////////////////////////////////////////////////////////////////////////
    1217// Private (file-static) functions
     
    3136    psMemSetDeallocator(stamp, (psFreeFunc)subtractionStampFree);
    3237
    33     stamp->x = 0;
    34     stamp->y = 0;
    35     stamp->xNorm = 0.0;
    36     stamp->yNorm = 0.0;
     38    stamp->x = NAN;
     39    stamp->y = NAN;
     40    stamp->xNorm = NAN;
     41    stamp->yNorm = NAN;
    3742    stamp->matrix = NULL;
    3843    stamp->vector = NULL;
     
    131136                    for (int x = xStart; x <= xStop ; x++) {
    132137                        if ((!subMask || !(subMask->data.PS_TYPE_MASK_DATA[y][x] &
    133                                            (PM_SUBTRACTION_MASK_BORDER | PM_SUBTRACTION_MASK_FOOTPRINT |
    134                                             PM_SUBTRACTION_MASK_REJ))) &&
     138                                           PM_SUBTRACTION_MASK_BAD_STAMP)) &&
    135139                            image->data.F32[y][x] > fluxBest) {
    136140                            fluxBest = image->data.F32[y][x];
     
    143147                stamp->x = xBest;
    144148                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;
    147149
    148150                // Reset the postage stamps since we're making a new stamp
     
    156158                    numFound++;
    157159                    psTrace("psModules.imcombine", 5, "Found stamp in region (%d,%d): %d,%d\n",
    158                             i, j, stamp->x, stamp->y);
     160                            i, j, (int)stamp->x, (int)stamp->y);
    159161                } else {
    160162                    stamp->status = PM_SUBTRACTION_STAMP_NONE;
     
    166168    psLogMsg("psModules.imcombine", PS_LOG_INFO, "Found %d stamps in %dx%d regions",
    167169             numFound, xNumStamps, yNumStamps);
     170
     171    return stamps;
     172}
     173
     174
     175psArray *pmSubtractionSetStamps(const psVector *x, const psVector *y,
     176                                const psImage *subMask, const psRegion *region)
     177{
     178    PS_ASSERT_VECTOR_NON_NULL(x, false);
     179    PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F32, false);
     180    PS_ASSERT_VECTOR_NON_NULL(y, false);
     181    PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F32, false);
     182    PS_ASSERT_VECTORS_SIZE_EQUAL(y, x, false);
     183    if (subMask) {
     184        PS_ASSERT_IMAGE_NON_NULL(subMask, false);
     185        PS_ASSERT_IMAGE_TYPE(subMask, PS_TYPE_MASK, false);
     186    }
     187
     188    int num = x->n;                     // Number of stamps
     189    psArray *stamps = psArrayAllocEmpty(num); // Stamps, to return
     190
     191    for (int i = 0; i < num; i++) {
     192        int xStamp = x->data.F32[i] + 0.5, yStamp = y->data.F32[i] + 0.5; // Pixel coords of stamp
     193        if (region && (xStamp < region->x0 || xStamp > region->x1 ||
     194                       yStamp < region->y0 || yStamp > region->y1)) {
     195            // Not interested at this time
     196            continue;
     197        }
     198
     199        if (subMask && subMask->data.PS_TYPE_MASK_DATA[yStamp][xStamp] & PM_SUBTRACTION_MASK_BAD_STAMP) {
     200            // It's bad
     201            continue;
     202        }
     203
     204        pmSubtractionStamp *stamp = pmSubtractionStampAlloc(PM_SUBTRACTION_STAMP_CALCULATE); // Stamp
     205        stamp->x = x->data.F32[i];
     206        stamp->y = y->data.F32[i];
     207        stamps->data[stamps->n++] = stamp;
     208    }
    168209
    169210    return stamps;
     
    198239    for (int i = 0; i < stamps->n; i++) {
    199240        pmSubtractionStamp *stamp = stamps->data[i]; // Stamp of interest
    200 
    201 
    202         int x = stamp->x, y = stamp->y; // Stamp coordinates
     241        if (!stamp || stamp->status != PM_SUBTRACTION_STAMP_CALCULATE) {
     242            continue;
     243        }
     244
     245        if (isnan(stamp->xNorm)) {
     246            stamp->xNorm = 2.0 * (stamp->x - (float)numCols/2.0) / (float)numCols;
     247        }
     248        if (isnan(stamp->yNorm)) {
     249            stamp->yNorm = 2.0 * (stamp->y - (float)numRows/2.0) / (float)numRows;
     250        }
     251
     252        int x = stamp->x + 0.5, y = stamp->y + 0.5; // Stamp coordinates
    203253        if (x < size || x > numCols - size || y < size || y > numRows - size) {
    204254            psError(PS_ERR_UNKNOWN, false, "Stamp %d (%d,%d) is within the image border.\n", i, x, y);
     
    228278
    229279
     280psArray *pmSubtractionGenerateStamps(const psVector *x, const psVector *y, const psVector *flux,
     281                                     float sigma, int footprint, const pmSubtractionKernels *kernels)
     282{
     283    PS_ASSERT_VECTOR_NON_NULL(x, false);
     284    PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F32, false);
     285    PS_ASSERT_VECTOR_NON_NULL(y, false);
     286    PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F32, false);
     287    PS_ASSERT_VECTORS_SIZE_EQUAL(y, x, false);
     288    PS_ASSERT_VECTOR_NON_NULL(flux, false);
     289    PS_ASSERT_VECTOR_TYPE(flux, PS_TYPE_F32, false);
     290    PS_ASSERT_VECTORS_SIZE_EQUAL(flux, y, false);
     291    PS_ASSERT_FLOAT_LARGER_THAN(sigma, 0.0, false);
     292
     293    int size = kernels->size + footprint; // Size of postage stamps
     294    int num = x->n;                     // Number of stamps
     295    psArray *stamps = psArrayAlloc(num); // Stamps, to return
     296
     297    for (int i = 0; i < num; i++) {
     298        pmSubtractionStamp *stamp = pmSubtractionStampAlloc(PM_SUBTRACTION_STAMP_CALCULATE); // Stamp
     299        stamps->data[i] = stamp;
     300
     301        stamp->x = x->data.F32[i];
     302        stamp->y = y->data.F32[i];
     303
     304        float xStamp = x->data.F32[i] - (int)(stamp->x + 0.5); // x coordinate of star in stamp frame
     305        float yStamp = y->data.F32[i] - (int)(stamp->y + 0.5); // y coordinate of star in stamp frame
     306        float fluxStamp = flux->data.F32[i]; // Flux of star
     307
     308        stamp->input = psKernelAlloc(-size, size, -size, size);
     309        psKernel *input = stamp->input; // Target stamp
     310
     311        // Put in a Gaussian, just for fun!
     312        for (int y = -size; y <= size; y++) {
     313            for (int x = -size; x <= size; x++) {
     314                input->kernel[y][x] = fluxStamp / sigma * 0.5 * M_2_SQRTPI * M_SQRT1_2 *
     315                    expf(0.5 * (PS_SQR(x + xStamp) + PS_SQR(y + yStamp)) / PS_SQR(sigma));
     316            }
     317        }
     318
     319    }
     320
     321    return stamps;
     322}
  • trunk/psModules/src/imcombine/pmSubtractionStamps.h

    r14480 r14513  
    1717/// A stamp for image subtraction
    1818typedef struct {
    19     int x, y;                           ///< Position
     19    float x, y;                         ///< Position
    2020    float xNorm, yNorm;                 ///< Normalised position
    2121    psKernel *reference;                ///< Reference image postage stamp
     
    3030psArray *pmSubtractionFindStamps(psArray *stamps, ///< Output stamps, or NULL
    3131                                 const psImage *image, ///< Image for which to find stamps
    32                                  const psImage *mask, ///< Mask
    33                                  const psRegion *region, ///< Region in which to search for stamps
     32                                 const psImage *mask, ///< Mask, or NULL
     33                                 const psRegion *region, ///< Region in which to search for stamps, or NULL
    3434                                 float threshold, ///< Threshold for stamps in the image
    3535                                 float spacing ///< Rough spacing for stamps
     36    );
     37
     38/// Set stamps based on a list of x,y
     39psArray *pmSubtractionSetStamps(const psVector *x, ///< x coordinates for each stamp
     40                                const psVector *y, ///< y coordinates for each stamp
     41                                const psImage *mask, ///< Mask, or NULL
     42                                const psRegion *region ///< Region in which to search for stamps, or NULL
    3643    );
    3744
     
    4552    );
    4653
     54/// Generate target for stamps based on list
     55psArray *pmSubtractionGenerateStamps(const psVector *x, ///< x coordinates for each stamp
     56                                     const psVector *y, ///< y coordinates for each stamp
     57                                     const psVector *flux, ///< Flux for each stamp
     58                                     float sigma, ///< Gaussian width for each stamp
     59                                     int footprint, ///< Stamp footprint size
     60                                     const pmSubtractionKernels *kernels ///< Kernel basis functions
     61    );
     62
    4763#endif
Note: See TracChangeset for help on using the changeset viewer.