IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Oct 29, 2008, 10:55:06 AM (18 years ago)
Author:
Paul Price
Message:

Don't rely on the input source list to have finite magnitudes (at the moment this is probably due to a mis-scaled weight map in the stacks), but use the image to set the flux of the stamp (which is only used to create the order of use). This provides some degree of flexibility, so that a subtraction can be performed even if the input is not optimal. Also, this way the multiple input methods (CMF, text list) are treated the same.

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r19234 r20465  
    262262    int numStamps = stamps->num;        // Number of stamp regions
    263263    int numFound = 0;                   // Number of stamps found
    264     int numValid = 0;                   // Number of valid regions
     264    int numSearch = 0;                  // Number of regions searched for new stamp
     265    int numGood = 0;                    // Number of good stamps in total
    265266    for (int i = 0; i < numStamps; i++) {
    266267        pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest
    267268
    268         // Only find a new stamp if we need to
    269         if (stamp->status != PM_SUBTRACTION_STAMP_REJECTED && stamp->status != PM_SUBTRACTION_STAMP_INIT) {
    270             continue;
    271         }
    272         numValid++;
    273 
    274         float xStamp = 0, yStamp = 0;   // Coordinates of stamp
    275         float fluxStamp = NAN;          // Flux of stamp
    276         bool goodStamp = false;         // Found a good stamp?
    277 
    278         // A couple different ways of finding stamps:
    279         if (stamps->x && stamps->y) {
    280             // Get the next stamp from the list
    281             psVector *xList = stamps->x->data[i], *yList = stamps->y->data[i]; // Coordinate lists
    282             psVector *fluxList = stamps->flux->data[i]; // List of stamp fluxes
    283 
    284             // Take stamp off the top of the (sorted) list
    285             if (xList->n > 0) {
    286                 int index = xList->n - 1; // Index of new stamp
    287                 xStamp = xList->data.F32[index];
    288                 yStamp = yList->data.F32[index];
    289                 fluxStamp = fluxList->data.F32[index];
    290 
    291                 // Chop off the top of the list
    292                 xList->n = index;
    293                 yList->n = index;
    294                 fluxList->n = index;
    295 
    296                 goodStamp = true;
    297             }
    298         } else {
    299             // Use a simple method of automatically finding stamps --- take the highest pixel in the subregion
    300             fluxStamp = threshold;
    301             psRegion *subRegion = stamps->regions->data[i]; // Sub-region of interest
    302             for (int y = subRegion->y0; y <= subRegion->y1; y++) {
    303                 for (int x = subRegion->x0; x <= subRegion->x1; x++) {
    304                     if (checkStampMask(x, y, subMask, mode, footprint) && image->data.F32[y][x] > fluxStamp) {
    305                         fluxStamp = image->data.F32[y][x];
    306                         xStamp = x;
    307                         yStamp = y;
    308                         goodStamp = true;
     269        switch (stamp->status) {
     270          case PM_SUBTRACTION_STAMP_NONE:
     271            continue;
     272          case PM_SUBTRACTION_STAMP_FOUND:
     273          case PM_SUBTRACTION_STAMP_CALCULATE:
     274          case PM_SUBTRACTION_STAMP_USED:
     275            numGood++;
     276            continue;
     277          case PM_SUBTRACTION_STAMP_INIT:
     278          case PM_SUBTRACTION_STAMP_REJECTED:
     279            numSearch++;
     280
     281            float xStamp = 0, yStamp = 0;   // Coordinates of stamp
     282            float fluxStamp = NAN;          // Flux of stamp
     283            bool goodStamp = false;         // Found a good stamp?
     284
     285            // A couple different ways of finding stamps:
     286            if (stamps->x && stamps->y) {
     287                // Get the next stamp from the list
     288                psVector *xList = stamps->x->data[i], *yList = stamps->y->data[i]; // Coordinate lists
     289                psVector *fluxList = stamps->flux->data[i]; // List of stamp fluxes
     290
     291                // Take stamp off the top of the (sorted) list
     292                if (xList->n > 0) {
     293                    int index = xList->n - 1; // Index of new stamp
     294                    xStamp = xList->data.F32[index];
     295                    yStamp = yList->data.F32[index];
     296                    fluxStamp = fluxList->data.F32[index];
     297
     298                    // Chop off the top of the list
     299                    xList->n = index;
     300                    yList->n = index;
     301                    fluxList->n = index;
     302
     303                    goodStamp = true;
     304                } else {
     305                    psTrace("psModules.imcombine", 9, "No sources in subregion %d", i);
     306                }
     307            } else {
     308                // Use a simple method of automatically finding stamps --- take the highest pixel in the
     309                // subregion
     310                fluxStamp = threshold;
     311                psRegion *subRegion = stamps->regions->data[i]; // Sub-region of interest
     312                for (int y = subRegion->y0; y <= subRegion->y1; y++) {
     313                    for (int x = subRegion->x0; x <= subRegion->x1; x++) {
     314                        if (checkStampMask(x, y, subMask, mode, footprint) &&
     315                            image->data.F32[y][x] > fluxStamp) {
     316                            fluxStamp = image->data.F32[y][x];
     317                            xStamp = x;
     318                            yStamp = y;
     319                            goodStamp = true;
     320                        }
    309321                    }
    310322                }
    311323            }
    312         }
    313 
    314         if (goodStamp) {
    315             stamp->x = xStamp;
    316             stamp->y = yStamp;
    317             stamp->flux = fluxStamp;
    318 
    319             // Reset the postage stamps since we're making a new stamp
    320             psFree(stamp->image1);
    321             psFree(stamp->image2);
    322             psFree(stamp->weight);
    323             psFree(stamp->convolutions1);
    324             psFree(stamp->convolutions2);
    325             stamp->image1 = stamp->image2 = stamp->weight = NULL;
    326             stamp->convolutions1 = stamp->convolutions2 = NULL;
    327 
    328             stamp->status = PM_SUBTRACTION_STAMP_FOUND;
    329             numFound++;
    330             psTrace("psModules.imcombine", 5, "Found stamp in subregion %d: %d,%d\n",
    331                     i, (int)stamp->x, (int)stamp->y);
    332         } else {
    333             stamp->status = PM_SUBTRACTION_STAMP_NONE;
    334         }
    335     }
    336 
    337     if (numValid > 0) {
     324
     325            if (goodStamp) {
     326                stamp->x = xStamp;
     327                stamp->y = yStamp;
     328                stamp->flux = fluxStamp;
     329
     330                // Reset the postage stamps since we're making a new stamp
     331                psFree(stamp->image1);
     332                psFree(stamp->image2);
     333                psFree(stamp->weight);
     334                psFree(stamp->convolutions1);
     335                psFree(stamp->convolutions2);
     336                stamp->image1 = stamp->image2 = stamp->weight = NULL;
     337                stamp->convolutions1 = stamp->convolutions2 = NULL;
     338
     339                stamp->status = PM_SUBTRACTION_STAMP_FOUND;
     340                numFound++;
     341                psTrace("psModules.imcombine", 5, "Found stamp in subregion %d: %d,%d\n",
     342                        i, (int)stamp->x, (int)stamp->y);
     343            } else {
     344                stamp->status = PM_SUBTRACTION_STAMP_NONE;
     345            }
     346        }
     347    }
     348
     349    if (numSearch > 0) {
    338350        psLogMsg("psModules.imcombine", PS_LOG_INFO, "Found %d stamps", numFound);
    339351    }
    340352
     353    if (numGood == 0 && numFound == 0) {
     354        // No good stamps
     355        psFree(stamps);
     356        return NULL;
     357    }
     358
    341359    return stamps;
    342360}
    343361
    344362
    345 pmSubtractionStampList *pmSubtractionStampsSet(const psVector *x, const psVector *y, const psVector *flux,
     363pmSubtractionStampList *pmSubtractionStampsSet(const psVector *x, const psVector *y,
    346364                                               const psImage *image, const psImage *subMask,
    347365                                               const psRegion *region, int footprint, float spacing,
    … …  
    354372    PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F32, NULL);
    355373    PS_ASSERT_VECTORS_SIZE_EQUAL(y, x, NULL);
    356     if (flux) {
    357         PS_ASSERT_VECTOR_NON_NULL(flux, NULL);
    358         PS_ASSERT_VECTOR_TYPE(flux, PS_TYPE_F32, NULL);
    359         PS_ASSERT_VECTORS_SIZE_EQUAL(flux, x, NULL);
    360     } else {
    361         PS_ASSERT_IMAGE_NON_NULL(image, NULL);
    362     }
     374    PS_ASSERT_IMAGE_NON_NULL(image, NULL);
    363375    if (subMask) {
    364376        PS_ASSERT_IMAGE_NON_NULL(subMask, NULL);
    … …  
    418430                xList->data.F32[index] = xStamp;
    419431                yList->data.F32[index] = yStamp;
    420 
    421                 if (flux) {
    422                     fluxList->data.F32[index] = flux->data.F32[i];
    423                 } else {
    424                     fluxList->data.F32[index] = image->data.F32[yPix][xPix];
    425                 }
     432                fluxList->data.F32[index] = image->data.F32[yPix][xPix];
    426433
    427434                found = true;
    … …  
    581588#endif
    582589
    583 pmSubtractionStampList *pmSubtractionStampsSetFromSources(const psArray *sources, const psImage *subMask,
    584                                                           const psRegion *region, int footprint,
    585                                                           float spacing, pmSubtractionMode mode)
     590pmSubtractionStampList *pmSubtractionStampsSetFromSources(const psArray *sources, const psImage *image,
     591                                                          const psImage *subMask, const psRegion *region,
     592                                                          int footprint, float spacing,
     593                                                          pmSubtractionMode mode)
    586594{
    587595    PS_ASSERT_ARRAY_NON_NULL(sources, NULL);
    … …  
    592600    psVector *x = psVectorAllocEmpty(numIn, PS_TYPE_F32); // x coordinates
    593601    psVector *y = psVectorAllocEmpty(numIn, PS_TYPE_F32); // y coordinates
    594     psVector *flux = psVectorAllocEmpty(numIn, PS_TYPE_F32); // Fluxes
    595602
    596603    int numOut = 0;                     // Number of sources in output list
    597604    for (int i = 0; i < numIn; i++) {
    598605        pmSource *source = sources->data[i]; // Source of interest
    599         if (!source || (source->mode & SOURCE_MASK) || !isfinite(source->psfMag) ||
    600             source->psfMag > SOURCE_FAINTEST) {
     606        if (!source || (source->mode & SOURCE_MASK) || source->psfMag > SOURCE_FAINTEST) {
    601607            continue;
    602608        }
    … …  
    608614            y->data.F32[numOut] = source->peak->yf;
    609615        }
    610         flux->data.F32[numOut] = powf(10.0, -0.4 * source->psfMag);
    611616        numOut++;
    612617    }
    613618    x->n = numOut;
    614619    y->n = numOut;
    615     flux->n = numOut;
    616 
    617     pmSubtractionStampList *stamps = pmSubtractionStampsSet(x, y, flux, NULL, subMask, region,
     620
     621    pmSubtractionStampList *stamps = pmSubtractionStampsSet(x, y, image, subMask, region,
    618622                                                            footprint, spacing, mode); // Stamps to return
    619623    psFree(x);
    620624    psFree(y);
    621     psFree(flux);
    622625
    623626    if (!stamps) {
    … …  
    647650    psBinaryOp(y, y, "-", psScalarAlloc(1.0, PS_TYPE_F32));
    648651
    649     pmSubtractionStampList *stamps = pmSubtractionStampsSet(x, y, NULL, image, subMask, region, footprint,
     652    pmSubtractionStampList *stamps = pmSubtractionStampsSet(x, y, image, subMask, region, footprint,
    650653                                                            spacing, mode);
    651654    psFree(data);
Note: See TracChangeset for help on using the changeset viewer.