IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 14802


Ignore:
Timestamp:
Sep 10, 2007, 11:31:25 AM (19 years ago)
Author:
Paul Price
Message:

Candidate stamps (or stamp regions, if no stamp list is provided) for the region are now generated in a single step before the iteration is performed.

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

Legend:

Unmodified
Added
Removed

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

    r14801 r14802  
    4747
    4848static bool getStamps(pmSubtractionStampList **stamps, // Stamps to read
    49                       const psArray *stampsData, // Stamp data from a file
    5049                      const pmReadout *reference, // Reference readout
    5150                      const pmReadout *input, // Input readout, or NULL to generate fake stamps
    … …  
    6059    )
    6160{
    62     if (stampsData && *stamps) {
    63         // We've already previously read all the stamps
    64         return true;
    65     }
    66 
    67     psImage *inImage = NULL; // Input image
    68     if (input) {
    69         inImage = input->image;
    70     }
    71 
    72     if (stampsData) {
    73         psVector *xStamp = NULL, *yStamp = NULL, *fluxStamp = NULL; // Stamp positions and fluxes
    74         if (input) {
    75             // We have x, y because the target is provided by the input image
    76             xStamp = stampsData->data[0];
    77             yStamp = stampsData->data[1];
    78         } else {
    79             // We have x, y and flux in order to generate a target
    80             xStamp = stampsData->data[0];
    81             yStamp = stampsData->data[1];
    82             fluxStamp = stampsData->data[2];
    83         }
    84 
    85         psFree(*stamps);
    86         // Apply exclusion zone if we're matching to a nominated PSF; otherwise don't care
    87         *stamps = pmSubtractionSetStamps(xStamp, yStamp, fluxStamp, reference->image, subMask,
    88                                          region, stampSpacing, input ? 0 : footprint);
    89     }
    9061    psTrace("psModules.imcombine", 3, "Finding stamps...\n");
    91     *stamps = pmSubtractionFindStamps(*stamps, reference->image, subMask, region,
    92                                       threshold, stampSpacing);
     62    *stamps = pmSubtractionStampsFind(*stamps, reference->image, subMask, region, threshold, stampSpacing);
    9363    if (!*stamps) {
    9464        psError(PS_ERR_UNKNOWN, false, "Unable to find stamps.");
    … …  
    9868    memCheck("  find stamps");
    9969
    100     if (!input && !pmSubtractionGenerateStamps(*stamps, targetWidth, footprint, size)) {
     70    if (!input && !pmSubtractionStampsGenerate(*stamps, targetWidth, footprint, size)) {
    10171        psError(PS_ERR_UNKNOWN, false, "Unable to generate target stamps.");
    10272        return false;
    10373    }
    10474
     75    memCheck("   generate stamps");
     76
    10577    psTrace("psModules.imcombine", 3, "Extracting stamps...\n");
    106     if (!pmSubtractionExtractStamps(*stamps, reference->image, inImage, weight, footprint, size)) {
     78    if (!pmSubtractionStampsExtract(*stamps, reference->image, input ? input->image : NULL,
     79                                    weight, footprint, size)) {
    10780        psError(PS_ERR_UNKNOWN, false, "Unable to extract stamps.");
    10881        return false;
    … …  
    225198    pmSubtractionKernels *kernels = NULL; // Kernel basis functions
    226199
    227     // Read stamps from file
    228     psArray *stampsData = NULL;         // Stamps data read from file
    229     if (stampsName && strlen(stampsName) > 0) {
    230         psTrace("psModules.imcombine", 3, "Reading stamps from %s...\n", stampsName);
    231         if (input) {
    232             // We have x, y because the target is provided by the input image
    233             stampsData = psVectorsReadFromFile(stampsName, "%f %f"); // Stamp positions
    234         } else {
    235             // We have x, y and flux in order to generate a target
    236             stampsData = psVectorsReadFromFile(stampsName, "%f %f %f"); // Stamp positions
    237         }
    238         if (!stampsData) {
    239             psError(PS_ERR_IO, false, "Unable to read stamps file %s", stampsName);
    240             return false;
    241         }
    242 
    243         // Correct for IRAF/FITS (unit-offset) positions to C (zero-offset) positions
    244         psVector *xStamp = stampsData->data[0]; // x positions
    245         psVector *yStamp = stampsData->data[1]; // y positions
    246         psBinaryOp(xStamp, xStamp, "-", psScalarAlloc(1.0, PS_TYPE_F32));
    247         psBinaryOp(yStamp, yStamp, "-", psScalarAlloc(1.0, PS_TYPE_F32));
    248     }
    249 
    250200    int numCols = reference->image->numCols, numRows = reference->image->numRows; // Image dimensions
    251201
    … …  
    281231            }
    282232
     233            // Read stamps from file
     234            if (stampsName && strlen(stampsName) > 0) {
     235                psTrace("psModules.imcombine", 3, "Reading stamps from %s...\n", stampsName);
     236                psVector *xStamp = NULL, *yStamp = NULL, *fluxStamp = NULL; // Stamp positions and fluxes
     237                if (input) {
     238                    // We have x, y because the target is provided by the input image
     239                    psArray *stampsData = psVectorsReadFromFile(stampsName, "%f %f"); // Stamp positions
     240                    if (!stampsData) {
     241                        psError(PS_ERR_IO, false, "Unable to read stamps file %s", stampsName);
     242                        goto ERROR;
     243                    }
     244                    xStamp = psMemIncrRefCounter(stampsData->data[0]);
     245                    yStamp = psMemIncrRefCounter(stampsData->data[1]);
     246                    psFree(stampsData);
     247                } else {
     248                    // We have x, y and flux in order to generate a target
     249                    psArray *stampsData = psVectorsReadFromFile(stampsName, "%f %f %f"); // Stamp positions
     250                    if (!stampsData) {
     251                        psError(PS_ERR_IO, false, "Unable to read stamps file %s", stampsName);
     252                        goto ERROR;
     253                    }
     254                    xStamp = psMemIncrRefCounter(stampsData->data[0]);
     255                    yStamp = psMemIncrRefCounter(stampsData->data[1]);
     256                    fluxStamp = psMemIncrRefCounter(stampsData->data[2]);
     257                    psFree(stampsData);
     258                }
     259
     260                // Correct for IRAF/FITS (unit-offset) positions to C (zero-offset) positions
     261                psBinaryOp(xStamp, xStamp, "-", psScalarAlloc(1.0, PS_TYPE_F32));
     262                psBinaryOp(yStamp, yStamp, "-", psScalarAlloc(1.0, PS_TYPE_F32));
     263
     264                stamps = pmSubtractionStampsSet(xStamp, yStamp, fluxStamp, reference->image, subMask,
     265                                                region, stampSpacing, input ? 0 : 2 * footprint);
     266            }
     267
    283268            // Define kernel basis functions
    284269            if (optimum && (type == PM_SUBTRACTION_KERNEL_ISIS || type == PM_SUBTRACTION_KERNEL_GUNK)) {
    285                 if (!getStamps(&stamps, stampsData, reference, input, subMask, weight, NULL,
     270                if (!getStamps(&stamps, reference, input, subMask, weight, NULL,
    286271                               threshold, stampSpacing, targetWidth, size, footprint)) {
    287272                    goto ERROR;
    … …  
    308293                psTrace("psModules.imcombine", 2, "Iteration %d...\n", k);
    309294
    310                 if (!getStamps(&stamps, stampsData, reference, input, subMask, weight, region,
     295                if (!getStamps(&stamps, reference, input, subMask, weight, region,
    311296                               threshold, stampSpacing, targetWidth, size, footprint)) {
    312297                    goto ERROR;
    … …  
    457442    psFree(subMask);
    458443    subMask = NULL;
    459     psFree(stampsData);
    460     stampsData = NULL;
    461444
    462445    if (!pmSubtractionBorder(convolved->image, convolved->weight, convolved->mask, size, maskBlank)) {
    … …  
    477460    psFree(stamps);
    478461    psFree(solution);
    479     psFree(stampsData);
    480462    psFree(stamps);
    481463    return false;

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

    r14801 r14802  
    55#include <stdio.h>
    66#include <pslib.h>
     7
     8// All these includes required to get stamps out of an array of pmSources
     9#include "pmMoments.h"
     10#include "pmPeaks.h"
     11#include "pmResiduals.h"
     12#include "pmHDU.h"
     13#include "pmFPA.h"
     14#include "pmGrowthCurve.h"
     15#include "pmPSF.h"
     16#include "pmModel.h"
     17#include "pmSource.h"
     18
    719
    820#include "pmSubtraction.h"
    … …  
    133145
    134146
    135 pmSubtractionStampList *pmSubtractionFindStamps(pmSubtractionStampList *stamps, const psImage *image,
     147pmSubtractionStampList *pmSubtractionStampsFind(pmSubtractionStampList *stamps, const psImage *image,
    136148                                                const psImage *subMask, const psRegion *region,
    137149                                                float threshold, float spacing)
    … …  
    170182    int numStamps = stamps->num;        // Number of stamp regions
    171183    int numFound = 0;                   // Number of stamps found
     184    int numValid = 0;                   // Number of valid regions
    172185    for (int i = 0; i < numStamps; i++) {
    173186        pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest
    … …  
    177190            continue;
    178191        }
     192        numValid++;
    179193
    180194        float xStamp = 0, yStamp = 0;   // Coordinates of stamp
    … …  
    238252    }
    239253
    240     psLogMsg("psModules.imcombine", PS_LOG_INFO, "Found %d stamps", numFound);
     254    if (numValid > 0) {
     255        psLogMsg("psModules.imcombine", PS_LOG_INFO, "Found %d stamps", numFound);
     256    }
    241257
    242258    return stamps;
    … …  
    244260
    245261
    246 pmSubtractionStampList *pmSubtractionSetStamps(const psVector *x, const psVector *y, const psVector *flux,
     262pmSubtractionStampList *pmSubtractionStampsSet(const psVector *x, const psVector *y, const psVector *flux,
    247263                                               const psImage *image, const psImage *subMask,
    248264                                               const psRegion *region, float spacing, int exclusionZone)
    249265{
    250     PS_ASSERT_VECTOR_NON_NULL(x, false);
    251     PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F32, false);
    252     PS_ASSERT_VECTOR_NON_NULL(y, false);
    253     PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F32, false);
    254     PS_ASSERT_VECTORS_SIZE_EQUAL(y, x, false);
     266    PS_ASSERT_VECTOR_NON_NULL(x, NULL);
     267    PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F32, NULL);
     268    PS_ASSERT_VECTOR_NON_NULL(y, NULL);
     269    PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F32, NULL);
     270    PS_ASSERT_VECTORS_SIZE_EQUAL(y, x, NULL);
    255271    if (flux) {
    256         PS_ASSERT_VECTOR_NON_NULL(flux, false);
    257         PS_ASSERT_VECTOR_TYPE(flux, PS_TYPE_F32, false);
    258         PS_ASSERT_VECTORS_SIZE_EQUAL(flux, x, false);
     272        PS_ASSERT_VECTOR_NON_NULL(flux, NULL);
     273        PS_ASSERT_VECTOR_TYPE(flux, PS_TYPE_F32, NULL);
     274        PS_ASSERT_VECTORS_SIZE_EQUAL(flux, x, NULL);
    259275    } else {
    260         PS_ASSERT_IMAGE_NON_NULL(image, false);
     276        PS_ASSERT_IMAGE_NON_NULL(image, NULL);
    261277    }
    262278    if (subMask) {
    263         PS_ASSERT_IMAGE_NON_NULL(subMask, false);
    264         PS_ASSERT_IMAGE_TYPE(subMask, PS_TYPE_MASK, false);
     279        PS_ASSERT_IMAGE_NON_NULL(subMask, NULL);
     280        PS_ASSERT_IMAGE_TYPE(subMask, PS_TYPE_MASK, NULL);
    265281        if (image) {
    266             PS_ASSERT_IMAGE_NON_NULL(image, false);
    267             PS_ASSERT_IMAGES_SIZE_EQUAL(image, subMask, false);
    268         }
    269     }
     282            PS_ASSERT_IMAGE_NON_NULL(image, NULL);
     283            PS_ASSERT_IMAGES_SIZE_EQUAL(image, subMask, NULL);
     284        }
     285    }
     286    PS_ASSERT_FLOAT_LARGER_THAN(spacing, 0.0, NULL);
     287    PS_ASSERT_INT_NONNEGATIVE(exclusionZone, NULL);
    270288
    271289    int numStars = x->n;                // Number of stars
    … …  
    399417
    400418
    401 bool pmSubtractionExtractStamps(pmSubtractionStampList *stamps, psImage *reference, psImage *input,
     419bool pmSubtractionStampsExtract(pmSubtractionStampList *stamps, psImage *reference, psImage *input,
    402420                                psImage *weight, int footprint, int kernelSize)
    403421{
    … …  
    467485
    468486
    469 bool pmSubtractionGenerateStamps(pmSubtractionStampList *stamps, float fwhm, int footprint, int kernelSize)
     487bool pmSubtractionStampsGenerate(pmSubtractionStampList *stamps, float fwhm, int footprint, int kernelSize)
    470488{
    471489    PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false);
    … …  
    510528    return true;
    511529}
     530
     531
     532pmSubtractionStampList *pmSubtractionStampsSetFromSources(const psArray *sources, const psImage *subMask,
     533                                                          const psRegion *region, float spacing,
     534                                                          int exclusionZone)
     535{
     536    PS_ASSERT_ARRAY_NON_NULL(sources, NULL);
     537    // Let pmSubtractionStampsSet take care of the rest of the assertions
     538
     539    int numSources = sources->n;          // Number of stars
     540
     541    psVector *x = psVectorAlloc(numSources, PS_TYPE_F32); // x coordinates
     542    psVector *y = psVectorAlloc(numSources, PS_TYPE_F32); // y coordinates
     543    psVector *flux = psVectorAlloc(numSources, PS_TYPE_F32); // Fluxes
     544
     545    for (int i = 0; i < numSources; i++) {
     546        pmSource *source = sources->data[i]; // Source of interest
     547        if (source->modelPSF) {
     548            x->data.F32[i] = source->modelPSF->params->data.F32[PM_PAR_XPOS];
     549            y->data.F32[i] = source->modelPSF->params->data.F32[PM_PAR_YPOS];
     550        } else {
     551            x->data.F32[i] = source->peak->xf;
     552            y->data.F32[i] = source->peak->yf;
     553        }
     554        flux->data.F32[i] = powf(10.0, -0.4 * source->psfMag);
     555    }
     556
     557    pmSubtractionStampList *stamps = pmSubtractionStampsSet(x, y, flux, NULL, subMask, region,
     558                                                            spacing, exclusionZone); // Stamps to return
     559    psFree(x);
     560    psFree(y);
     561    psFree(flux);
     562
     563    if (!stamps) {
     564        psError(PS_ERR_UNKNOWN, false, "Unable to set stamps from sources.");
     565    }
     566
     567    return stamps;
     568}

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

    r14801 r14802  
    6767
    6868/// Find stamps on an image
    69 pmSubtractionStampList *pmSubtractionFindStamps(pmSubtractionStampList *stamps, ///< Output stamps, or NULL
     69pmSubtractionStampList *pmSubtractionStampsFind(pmSubtractionStampList *stamps, ///< Output stamps, or NULL
    7070                                                const psImage *image, ///< Image for which to find stamps
    7171                                                const psImage *mask, ///< Mask, or NULL
    … …  
    7979/// We may optionally apply an exclusion zone around each star --- this is important when we're convolving to
    8080/// a specified PSF; less so when we're only trying to get a good subtraction.
    81 pmSubtractionStampList *pmSubtractionSetStamps(const psVector *x, ///< x coordinates for each stamp
     81pmSubtractionStampList *pmSubtractionStampsSet(const psVector *x, ///< x coordinates for each stamp
    8282                                               const psVector *y, ///< y coordinates for each stamp
    8383                                               const psVector *flux, ///< Flux for each stamp, or NULL
    … …  
    8686                                               const psRegion *region, ///< Region to search, or NULL
    8787                                               float spacing, ///< Rough spacing for stamps
    88                                                int exclusionZone // Exclusion zone around each star
     88                                               int exclusionZone ///< Exclusion zone around each star
     89    );
     90
     91///
     92pmSubtractionStampList *pmSubtractionStampsSetFromSources(
     93    const psArray *sources, ///< Sources for each stamp
     94    const psImage *subMask, ///< Mask, or NULL
     95    const psRegion *region, ///< Region to search, or NULL
     96    float spacing, ///< Rough spacing for stamps
     97    int exclusionZone ///< Exclusion zone around each star
    8998    );
    9099
    91100/// Extract stamps from the images
    92 bool pmSubtractionExtractStamps(pmSubtractionStampList *stamps, ///< Stamps
     101bool pmSubtractionStampsExtract(pmSubtractionStampList *stamps, ///< Stamps
    93102                                psImage *reference, ///< Reference image
    94103                                psImage *input, ///< Input image (or NULL)
    … …  
    99108
    100109/// Generate target for stamps based on list
    101 bool pmSubtractionGenerateStamps(pmSubtractionStampList *stamps, ///< Stamps
     110bool pmSubtractionStampsGenerate(pmSubtractionStampList *stamps, ///< Stamps
    102111                                 float fwhm, ///< FWHM for each stamp
    103112                                 int footprint, ///< Stamp footprint size
Note: See TracChangeset for help on using the changeset viewer.