IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 30141


Ignore:
Timestamp:
Dec 21, 2010, 3:01:12 PM (15 years ago)
Author:
eugene
Message:

for psf-matched images: re-determine the psf; re-generate the models with the new psf; re-fit the fluxes; subtract all sources; as radial aperture fluxes are measured only replace the current source of interest; replace all sources before convolving to the next psf

Location:
branches/eam_branches/ipp-20101205/psphot/src
Files:
8 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ipp-20101205/psphot/src/psphot.h

    r30101 r30141  
    7575bool            psphotImageQualityReadout(pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe);
    7676
    77 bool            psphotChoosePSF (pmConfig *config, const pmFPAview *view, const char *filerule);
    78 bool            psphotChoosePSFReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe);
     77bool            psphotChoosePSF (pmConfig *config, const pmFPAview *view, const char *filerule, bool newSources);
     78bool            psphotChoosePSFReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe, bool newSources);
    7979
    8080bool            psphotGuessModels (pmConfig *config, const pmFPAview *view, const char *filerule);
    … …  
    404404bool psphotStackMatchPSFsReadout (pmConfig *config, const pmFPAview *view, psphotStackOptions *options, int index);
    405405bool psphotStackMatchPSFsPrepare (pmConfig *config, const pmFPAview *view, psphotStackOptions *options, int index);
    406 bool psphotStackMatchPSFsNext (bool *smoothAgain, pmConfig *config, const pmFPAview *view, const char *filerule, int nextSize);
    407 bool psphotStackMatchPSFsNextReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe, float currentFWHM, float targetFWHM);
     406bool psphotStackMatchPSFsNext (bool *smoothAgain, pmConfig *config, const pmFPAview *view, const char *filerule, int lastSize);
     407bool psphotStackMatchPSFsNextReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe, psVector *fwhmValues, int lastSize);
    408408
    409409// psphotStackMatchPSFsUtils
    … …  
    443443bool psphotCleanInputs (pmConfig *config, const pmFPAview *view, const char *filerule);
    444444
     445bool psphotResetModels (pmConfig *config, const pmFPAview *view, const char *filerule);
     446bool psphotResetModelsReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe);
     447
     448bool psphotRedefinePixels (pmConfig *config, const pmFPAview *view, const char *filerule);
     449bool psphotRedefinePixelsReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe);
     450
    445451#endif
     452

  • branches/eam_branches/ipp-20101205/psphot/src/psphotChoosePSF.c

    r30108 r30141  
    22
    33// generate a PSF model for inputs without PSF models already loaded
    4 bool psphotChoosePSF (pmConfig *config, const pmFPAview *view, const char *filerule)
     4bool psphotChoosePSF (pmConfig *config, const pmFPAview *view, const char *filerule, bool newSources)
    55{
    66    bool status = true;
    … …  
    1919    for (int i = 0; i < num; i++) {
    2020        if (i == chisqNum) continue; // skip chisq image
    21         if (!psphotChoosePSFReadout (config, view, filerule, i, recipe)) {
     21        if (!psphotChoosePSFReadout (config, view, filerule, i, recipe, newSources)) {
    2222            psError (PSPHOT_ERR_CONFIG, false, "failed to choose a psf model for %s entry %d", filerule, i);
    2323            return false;
    … …  
    2828
    2929// try PSF models and select best option
    30 bool psphotChoosePSFReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe) {
     30bool psphotChoosePSFReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe, bool newSources) {
    3131
    3232    bool status;
    … …  
    5050    psAssert (detections, "missing detections?");
    5151
    52     psArray *sources = detections->newSources;
     52    psArray *sources = newSources ? detections->newSources : detections->allSources;
    5353    psAssert (sources, "missing sources?");
    5454

  • branches/eam_branches/ipp-20101205/psphot/src/psphotCleanup.c

    r29548 r30141  
    2222    pmVisualCleanup ();
    2323    psLibFinalize();
    24     // fprintf (stderr, "found %d leaks at %s\n", psMemCheckLeaks (0, NULL, NULL, false), "psphot");
    25     fprintf (stderr, "found %d leaks at %s\n", psMemCheckLeaks (0, NULL, stdout, false), "psphot");
     24    fprintf (stderr, "found %d leaks at %s\n", psMemCheckLeaks (0, NULL, NULL, false), "psphot");
     25    // fprintf (stderr, "found %d leaks at %s\n", psMemCheckLeaks (0, NULL, stdout, false), "psphot");
    2626    return;
    2727}

  • branches/eam_branches/ipp-20101205/psphot/src/psphotReadoutKnownSources.c

    r29936 r30141  
    4343    }
    4444
    45     if (!psphotChoosePSF (config, view, filerule)) {
     45    if (!psphotChoosePSF (config, view, filerule, true)) {
    4646        psError(PSPHOT_ERR_PSF, false, "Failed to construct a psf model");
    4747        return psphotReadoutCleanup (config, view, filerule);

  • branches/eam_branches/ipp-20101205/psphot/src/psphotReplaceUnfit.c

    r29936 r30141  
    7575      pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
    7676    }
     77
     78    psphotVisualShowImage(readout);
    7779    psLogMsg ("psphot.replace", PS_LOG_INFO, "replaced models for %ld objects: %f sec\n", sources->n, psTimerMark ("psphot.replace"));
    7880    return true;
    … …  
    101103    return true;
    102104}
     105
     106// modify the sources to point at the corresponding pixels for the given filerule
     107bool psphotRedefinePixels (pmConfig *config, const pmFPAview *view, const char *filerule)
     108{
     109    bool status = true;
     110
     111    // select the appropriate recipe information
     112    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
     113    psAssert (recipe, "missing recipe?");
     114
     115    int num = psphotFileruleCount(config, filerule);
     116
     117    // loop over the available readouts
     118    for (int i = 0; i < num; i++) {
     119        if (!psphotRedefinePixelsReadout (config, view, filerule, i, recipe)) {
     120            psError (PSPHOT_ERR_CONFIG, false, "failed to replace all sources for %s entry %d", filerule, i);
     121            return false;
     122        }
     123    }
     124    return true;
     125}
     126
     127bool psphotRedefinePixelsReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe) {
     128
     129    bool status;
     130    pmSource *source;
     131
     132    psTimerStart ("psphot.replace");
     133
     134    // find the currently selected readout
     135    pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, index); // File of interest
     136    psAssert (file, "missing file?");
     137
     138    pmReadout *readout = pmFPAviewThisReadout(view, file->fpa);
     139    psAssert (readout, "missing readout?");
     140
     141    // XXX the sources have already been copied (merge into here?)
     142    pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
     143    psAssert (detections, "missing detections?");
     144
     145    psArray *sources = detections->allSources;
     146    psAssert (sources, "missing sources?");
     147
     148    // user-defined masks to test for good/bad pixels (build from recipe list if not yet set)
     149    psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels
     150    psAssert (maskVal, "missing mask value?");
     151
     152    for (int i = 0; i < sources->n; i++) {
     153      source = sources->data[i];
     154
     155      // sources have not yet been subtracted in this image (but this flag may be raised)
     156      source->tmpFlags &= ~PM_SOURCE_TMPF_SUBTRACTED;
     157
     158      float Xo = source->modelPSF->params->data.F32[PM_PAR_XPOS];
     159      float Yo = source->modelPSF->params->data.F32[PM_PAR_YPOS];
     160      float radius = source->modelPSF->fitRadius;
     161      if (radius < 1) {
     162          fprintf (stderr, "!");
     163      }
     164
     165      pmSourceRedefinePixels (source, readout, Xo, Yo, radius, true);
     166    }
     167    return true;
     168}
     169
     170// for now, let's store the detections on the readout->analysis for each readout
     171bool psphotResetModels (pmConfig *config, const pmFPAview *view, const char *filerule)
     172{
     173    bool status = true;
     174
     175    // select the appropriate recipe information
     176    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
     177    psAssert (recipe, "missing recipe?");
     178
     179    int num = psphotFileruleCount(config, filerule);
     180
     181    // loop over the available readouts
     182    for (int i = 0; i < num; i++) {
     183        if (!psphotResetModelsReadout (config, view, filerule, i, recipe)) {
     184            psError (PSPHOT_ERR_CONFIG, false, "failed to replace all sources for %s entry %d", filerule, i);
     185            return false;
     186        }
     187    }
     188    return true;
     189}
     190
     191bool psphotResetModelsReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe) {
     192
     193    bool status;
     194    pmSource *source;
     195
     196    psTimerStart ("psphot.replace");
     197
     198    // find the currently selected readout
     199    pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, index); // File of interest
     200    psAssert (file, "missing file?");
     201
     202    pmReadout *readout = pmFPAviewThisReadout(view, file->fpa);
     203    psAssert (readout, "missing readout?");
     204
     205    // XXX the sources have already been copied (merge into here?)
     206    pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
     207    psAssert (detections, "missing detections?");
     208
     209    psArray *sources = detections->allSources;
     210    psAssert (sources, "missing sources?");
     211
     212    // user-defined masks to test for good/bad pixels (build from recipe list if not yet set)
     213    psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels
     214    psAssert (maskVal, "missing mask value?");
     215
     216    pmPSF *psf = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF");
     217    psAssert (psf, "missing psf?");
     218
     219    for (int i = 0; i < sources->n; i++) {
     220      source = sources->data[i];
     221
     222      // sources have not yet been subtracted in this image (but this flag may be raised)
     223      source->tmpFlags &= ~PM_SOURCE_TMPF_SUBTRACTED;
     224
     225      // the guess central intensity comes from the peak:
     226      float Io = source->peak->flux;
     227      float Xo = source->modelPSF->params->data.F32[PM_PAR_XPOS];
     228      float Yo = source->modelPSF->params->data.F32[PM_PAR_YPOS];
     229      float radius = source->modelPSF->fitRadius;
     230
     231      // generate a model for this object with Io = 1.0
     232      pmModel *modelPSF = pmModelFromPSFforXY(psf, Xo, Yo, Io);
     233      if (modelPSF == NULL) {
     234          psWarning ("Failed to determine PSF model for source at (%f,%f), skipping", Xo, Yo);
     235          continue;
     236      }
     237
     238      // set the source PSF model
     239      psFree (source->modelPSF);
     240      source->modelPSF = modelPSF;
     241      source->modelPSF->fitRadius = radius;
     242
     243      pmSourceCacheModel (source, maskVal);  // ALLOC x14 (!)
     244    }
     245
     246    // psphotSaveImage(NULL, readoutIn->image, "image.in.fits");
     247    // psphotSaveImage(NULL, readoutOut->image, "image.out.sub.fits");
     248    // psphotVisualShowImage (readout);
     249
     250    psLogMsg ("psphot.replace", PS_LOG_INFO, "subtracted models for %ld objects: %f sec\n", sources->n, psTimerMark ("psphot.replace"));
     251    return true;
     252}
     253

  • branches/eam_branches/ipp-20101205/psphot/src/psphotSourceMatch.c

    r30096 r30141  
    243243            peak->footprint = pmFootprintCopyData(footprint, readout->image);
    244244
     245            // the peak does not claim ownership of the footprint (it does not free it). save a copy of this
     246            // footprint on detections->footprints so we can free it later
     247            psArrayAdd(detections->footprints, 100, peak->footprint);
     248            psFree (peak->footprint);
     249           
    245250            // create a new source
    246251            pmSource *source = pmSourceAlloc();
    … …  
    249254
    250255            // add the peak
    251             source->peak = psMemIncrRefCounter(peak);
     256            source->peak = peak;
    252257
    253258            // allocate space for moments
    … …  
    261266            psArrayAdd (detections->newSources, 100, source);
    262267            psFree (source);
    263             psFree (peak);
    264268        }
    265269        psFree (footprint);

  • branches/eam_branches/ipp-20101205/psphot/src/psphotStackMatchPSFsNext.c

    r30101 r30141  
    2626    }
    2727
    28     float targetFWHM = fwhmValues->data.F32[lastSize + 1];
    29     float currentFWHM = fwhmValues->data.F32[lastSize];
    30 
    3128    int num = psphotFileruleCount(config, filerule);
    3229
    … …  
    3633    // loop over the available readouts
    3734    for (int i = 0; i < num; i++) {
    38         if (!psphotStackMatchPSFsNextReadout (config, view, filerule, i, recipe, currentFWHM, targetFWHM)) {
    39             psError (PSPHOT_ERR_CONFIG, false, "failed to smooth image %s (%d) to target PSF %f", filerule, i, targetFWHM);
     35        if (!psphotStackMatchPSFsNextReadout (config, view, filerule, i, recipe, fwhmValues, lastSize)) {
     36            psError (PSPHOT_ERR_CONFIG, false, "failed to smooth image %s (%d) to target PSF %f", filerule, i, fwhmValues->data.F32[lastSize+1]);
    4037            psImageConvolveSetThreads(oldThreads);
    4138            return false;
    … …  
    4744}
    4845
    49 bool psphotStackMatchPSFsNextReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe, float currentFWHM, float targetFWHM) {
     46bool psphotStackMatchPSFsNextReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe, psVector *fwhmValues, int lastSize) {
    5047
    5148    bool status = false;
    … …  
    7067        minGauss = 0.5;
    7168    }
     69
     70    float targetFWHM = fwhmValues->data.F32[lastSize + 1];
     71    float currentFWHM = fwhmValues->data.F32[lastSize];
    7272
    7373    if (targetFWHM <= currentFWHM) {
    … …  
    131131    // so we save this in a vector.  if the vector is not yet defined, create it
    132132
    133     psVector *fwhmValues = psMetadataLookupVector(&status, readout->analysis, "STACK.PSF.FWHM.VALUES");
    134     psAssert(fwhmValues, "should already exist..");
    135     psVectorAppend(fwhmValues, targetFWHM);
     133    psVector *fwhmValuesOut = psMetadataLookupVector(&status, readout->analysis, "STACK.PSF.FWHM.VALUES");
     134    psAssert(fwhmValuesOut, "should already exist..");
     135    psVectorAppend(fwhmValuesOut, targetFWHM);
     136
     137    // do not generate a PSF if we already were supplied one
     138    pmPSF *psfOld = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF");
     139    if (psfOld) {
     140        // save PSF on readout->analysis
     141        char psfEntry[64];
     142        snprintf (psfEntry, 64, "PSPHOT.PSF.V%d", lastSize);
     143        if (!psMetadataAddPtr (readout->analysis, PS_LIST_TAIL, psfEntry, PS_META_REPLACE | PS_DATA_UNKNOWN, "psphot psf model", psfOld)) {
     144            psError (PSPHOT_ERR_UNKNOWN, false, "problem saving sources on readout");
     145            return false;
     146        }
     147        psMetadataRemoveKey(readout->analysis, "PSPHOT.PSF");
     148    }
    136149
    137150    return true;

  • branches/eam_branches/ipp-20101205/psphot/src/psphotStackReadout.c

    r30101 r30141  
    1010// SRC (source analysis image) : nominally CNV (optionally RAW)
    1111// OUT (psf-matched images)    : always OUT
     12
     13bool psphotStackVisualFilerule(pmConfig *config, const pmFPAview *view, const char *filerule) {
     14
     15    int num = psphotFileruleCount(config, filerule);
     16
     17    // loop over the available readouts
     18    for (int i = 0; i < num; i++) {
     19
     20        // find the currently selected readout
     21        pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, i); // File of interest
     22        psAssert (file, "missing file?");
     23
     24        pmReadout *readout = pmFPAviewThisReadout(view, file->fpa);
     25        psAssert (readout, "missing readout?");
     26
     27        psphotVisualShowImage (readout);
     28        psphotVisualShowResidualImage (readout);
     29    }
     30    return true;
     31}
    1232
    1333bool psphotStackReadout (pmConfig *config, const pmFPAview *view) {
    … …  
    92112    }
    93113
     114    if (!strcasecmp (breakPt, "TEST1")) {
     115        return psphotReadoutCleanup (config, view, STACK_SRC);
     116    }
     117
    94118    // generate the objects (object unify the sources from the different images)
    95119    // XXX this could just match the detections for the chisq image, and not bother measuring the
    … …  
    97121    psArray *objects = psphotMatchSources (config, view, STACK_SRC);
    98122
     123    if (!strcasecmp (breakPt, "TEST2")) {
     124        psFree(objects);
     125        return psphotReadoutCleanup (config, view, STACK_SRC);
     126    }
     127
    99128    // construct sources for the newly-generated sources (from other images)
    100129    if (!psphotSourceStats (config, view, STACK_SRC, false)) { // pass 1
     130        psFree(objects);
    101131        psError(PSPHOT_ERR_UNKNOWN, false, "failure to generate sources");
    102132        return psphotReadoutCleanup (config, view, STACK_SRC);
    … …  
    118148    // only run this on detections from the input images, not chisq image
    119149    if (!psphotRoughClass (config, view, STACK_SRC)) {
     150        psFree(objects);
    120151        psError (PSPHOT_ERR_UNKNOWN, false, "failed to determine rough classifications");
    121152        return psphotReadoutCleanup (config, view, STACK_SRC);
    … …  
    124155    // only run this on detections from the input images, not chisq image
    125156    if (!psphotImageQuality (config, view, STACK_SRC)) { // pass 1
     157        psFree(objects);
    126158        psError (PSPHOT_ERR_UNKNOWN, false, "failed to measure image quality");
    127159        return psphotReadoutCleanup (config, view, STACK_SRC);
    128160    }
    129161    if (!strcasecmp (breakPt, "MOMENTS")) {
     162        psFree(objects);
    130163        return psphotReadoutCleanup (config, view, STACK_SRC);
    131164    }
    … …  
    133166    // use bright stellar objects to measure PSF if we were supplied a PSF for any input file,
    134167    // this step is skipped
    135     if (!psphotChoosePSF (config, view, STACK_SRC)) { // pass 1
     168    if (!psphotChoosePSF (config, view, STACK_SRC, true)) { // pass 1
     169        psFree(objects);
    136170        psLogMsg ("psphot", 3, "failure to construct a psf model");
    137171        return psphotReadoutCleanup (config, view, STACK_SRC);
    138172    }
    139173    if (!strcasecmp (breakPt, "PSFMODEL")) {
     174        psFree(objects);
    140175        return psphotReadoutCleanup (config, view, STACK_SRC);
    141176    }
    … …  
    150185    // linear PSF fit to source peaks, subtract the models from the image (in PSF mask)
    151186    psphotFitSourcesLinearStack (config, objects, FALSE);
     187    psphotStackVisualFilerule(config, view, STACK_SRC);
    152188
    153189    // identify CRs and extended sources
    … …  
    174210    // copy the detections from SRC to OUT (for radial aperture photometry)
    175211    if (!psphotCopySources (config, view, STACK_OUT, STACK_SRC)) {
     212        psFree(objects);
    176213        psError (PSPHOT_ERR_UNKNOWN, false, "failure in peak analysis");
    177214        return psphotReadoutCleanup (config, view, STACK_SRC);
    178215    }
    179     // XXX need to do something to reassign the source pixels here
    180216
    181217    bool smoothAgain = true;
    182218    for (int nMatchedPSF = 0; smoothAgain; nMatchedPSF++) {
     219        psphotRedefinePixels (config, view, STACK_OUT);
     220
     221        psphotChoosePSF (config, view, STACK_OUT, false);
     222
     223        psphotResetModels (config, view, STACK_OUT);
     224
     225        psphotFitSourcesLinear (config, view, STACK_OUT, false);
     226
    183227        // measure circular, radial apertures (objects sorted by S/N)
    184228        psphotRadialAperturesByObject (config, objects, view, STACK_OUT, nMatchedPSF);
     229
     230        psphotReplaceAllSources (config, view, STACK_OUT);
     231
    185232        psphotStackMatchPSFsNext(&smoothAgain, config, view, STACK_OUT, nMatchedPSF);
    186233    }
Note: See TracChangeset for help on using the changeset viewer.