IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 32322


Ignore:
Timestamp:
Sep 5, 2011, 9:02:13 AM (15 years ago)
Author:
eugene
Message:

major re-work of psphotStackReadout to be more consistent with psphotReadout; psphotRadialApertures now assumes the 0 element is the unconvolved version; for the convolved versions, the vector of target psfs is saved as STACK.PSF.FWHM.VALUES

Location:
branches/eam_branches/ipp-20110710/psphot/src
Files:
7 edited

Legend:

Unmodified
Added
Removed

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

    r32299 r32322  
    377377    bool convolve;                      // Convolve images?
    378378    psphotStackConvolveSource convolveSource;
    379     // psArray *convImages, *convMasks, *convVariances; // Filenames for the temporary convolved images
    380 
    381     // bool matchZPs;                      // Adjust relative fluxes based on transparency analysis?
    382     // bool photometry;                    // Perform photometry?
    383     // psMetadata *stats;                  // Statistics for output
    384     // FILE *statsFile;                    // File to which to write statistics
    385     // psArray *origImages, *origMasks, *origVariances; // Filenames of the original images
    386     // psArray *origCovars;                // Original covariances matrices
    387     // int quality;                        // Bad data quality flag
    388379
    389380    // Prepare
    … …  
    392383    psVector *inputMask;                // Mask for inputs
    393384
    394     float targetSeeing;                 // Target seeing FWHM
     385    psVector *targetSeeing;             // Target seeing FWHMs
    395386    psArray *sourceLists;               // Individual lists of sources for matching
    396387    psVector *norm;                     // Normalisation for each image
    397388    psArray *psfs;
    398 
    399     // psVector *exposures;                // Exposure times
    400     // float sumExposure;                  // Sum of exposure times
    401     // float zp;                           // Zero point for output
    402     // psVector *inputMask;                // Mask for inputs
    403     // psArray *sources;                   // Matched sources
    404389
    405390    // Convolve
    … …  
    407392    psArray *regions;                   // PSF-matching regions --- required in the stacking
    408393    psVector *matchChi2;                // chi^2 for stamps from matching
    409     psVector *weightings;               // Combination weightings for images (1/noise^2)
    410     // psArray *cells;                     // Cells for convolved images --- a handle for reading again
    411     // int numCols, numRows;               // Size of image
    412     // psArray *convCovars;                // Convolved covariance matrices
    413 
    414     // Combine initial
    415     // pmReadout *outRO;                   // Output readout
    416     // pmReadout *expRO;                   // Exposure readout
    417     // psArray *inspect;                   // Array of arrays of pixels to inspect
    418 
    419     // Rejection
    420     // psArray *rejected;                  // Rejected pixels
    421394} psphotStackOptions;
    422395
    … …  
    425398bool psphotStackMatchPSFsReadout (pmConfig *config, const pmFPAview *view, psphotStackOptions *options, int index);
    426399bool psphotStackMatchPSFsPrepare (pmConfig *config, const pmFPAview *view, psphotStackOptions *options, int index);
    427 bool psphotStackMatchPSFsNext (bool *smoothAgain, pmConfig *config, const pmFPAview *view, const char *filerule, int lastSize);
    428 bool psphotStackMatchPSFsNextReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe, psVector *fwhmValues, int lastSize);
     400
     401bool psphotStackMatchPSFsNext (pmConfig *config, const pmFPAview *view, const char *filerule, int lastSize);
     402bool psphotStackMatchPSFsNextReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, int lastSize);
     403int psphotStackMatchPSFsEntries (pmConfig *config, const pmFPAview *view, const char *filerule);
    429404
    430405// psphotStackMatchPSFsUtils
    431 psVector *SetOptWidths (bool *optimum, psMetadata *recipe);
    432 pmReadout *makeFakeReadout(pmConfig *config, pmReadout *raw, psArray *sources, pmPSF *psf, psImageMaskType maskVal, int fullSize);
    433 bool rescaleData(pmReadout *readout, pmConfig *config, psphotStackOptions *options, int index);
    434 bool renormKernel(pmReadout *readout, psphotStackOptions *options, int index);
     406// psVector *SetOptWidths (bool *optimum, psMetadata *recipe);
     407// pmReadout *makeFakeReadout(pmConfig *config, pmReadout *raw, psArray *sources, pmPSF *psf, psImageMaskType maskVal, int fullSize);
     408// bool rescaleData(pmReadout *readout, pmConfig *config, psphotStackOptions *options, int index);
     409// bool renormKernel(pmReadout *readout, psphotStackOptions *options, int index);
    435410bool saveMatchData (pmReadout *readout, psphotStackOptions *options, int index);
    436411bool matchKernel(pmConfig *config, pmReadout *cnv, pmReadout *raw, psphotStackOptions *options, int index);
    437 bool dumpImageDiff(pmReadout *readoutConv, pmReadout *readoutFake, pmReadout *readoutRef, int index, char *rootname);
    438 bool dumpImage(pmReadout *readoutOut, pmReadout *readoutRef, int index, char *rootname);
    439 bool loadKernel (pmConfig *config, pmReadout *readoutCnv, psphotStackOptions *options, int index);
    440 
    441 pmPSF *psphotStackPSF(const pmConfig *config, int numCols, int numRows, const psArray *psfs, const psVector *inputMask);
     412// bool dumpImageDiff(pmReadout *readoutConv, pmReadout *readoutFake, pmReadout *readoutRef, int index, char *rootname);
     413// bool dumpImage(pmReadout *readoutOut, pmReadout *readoutRef, int index, char *rootname);
     414// bool loadKernel (pmConfig *config, pmReadout *readoutCnv, psphotStackOptions *options, int index);
     415
     416bool psphotStackRenormaliseVariance(const pmConfig *config, pmReadout *readout);
     417
     418bool psphotStackPSF(const pmConfig *config, psphotStackOptions *options);
    442419
    443420psphotStackOptions *psphotStackOptionsAlloc (int num);
    … …  
    448425bool psphotCopySourcesReadout (pmConfig *config, const pmFPAview *view, const char *ruleOut, const char *ruleSrc, int index);
    449426
    450 bool psphotRadialApertures (pmConfig *config, const pmFPAview *view, const char *filerule);
    451 bool psphotRadialAperturesReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe);
     427bool psphotRadialApertures (pmConfig *config, const pmFPAview *view, const char *filerule, int entry);
     428bool psphotRadialAperturesReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe, int entry);
    452429bool psphotRadialApertures_Threaded (psThreadJob *job);
    453430bool psphotRadialApertureSource (pmSource *source, psMetadata *recipe, psImageMaskType maskVal, const psVector *radMax, int entry);
    … …  
    489466bool psphotKronIterate_Threaded (psThreadJob *job);
    490467
     468bool psphotStackObjectsSelectForAnalysis (pmConfig *config, const pmFPAview *view, const char *filerule, psArray *objects);
     469
    491470#endif

  • branches/eam_branches/ipp-20110710/psphot/src/psphotRadialApertures.c

    r32238 r32322  
    33bool psphotRadialAperturesSortFlux (psVector *radius, psVector *pixFlux, psVector *pixVar);
    44
     5// this function measures the radial aperture fluxes for the set of readouts.  this function
     6// may be called multiple times, presumably for different versions of PSF-matched or unmatched images. 
     7
    58// for now, let's store the detections on the readout->analysis for each readout
    6 bool psphotRadialApertures (pmConfig *config, const pmFPAview *view, const char *filerule)
     9bool psphotRadialApertures (pmConfig *config, const pmFPAview *view, const char *filerule, int entry)
    710{
    811    bool status = true;
     12
     13    fprintf (stdout, "\n");
     14    psLogMsg ("psphot", PS_LOG_INFO, "--- psphot Radial Apertures ---");
    915
    1016    // select the appropriate recipe information
    … …  
    2228    // loop over the available readouts
    2329    for (int i = 0; i < num; i++) {
    24         if (!psphotRadialAperturesReadout (config, view, filerule, i, recipe)) {
     30        if (!psphotRadialAperturesReadout (config, view, filerule, i, recipe, entry)) {
    2531            psError (PSPHOT_ERR_CONFIG, false, "failed on measure extended source aperture-like parameters for %s entry %d", filerule, i);
    2632            return false;
    … …  
    3238// aperture-like measurements for extended sources
    3339// flux in simple, circular apertures
    34 bool psphotRadialAperturesReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe) {
     40// 'entry' tells us which of the matched-PSF images we are working on (0 == unmatched image, also non-stack psphot)
     41bool psphotRadialAperturesReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe, int entry) {
    3542
    3643    bool status;
    … …  
    7279    // XXX are we getting the objects out of order? does it matter?
    7380    sources = psArraySort (sources, pmSourceSortByFlux);
     81
     82    // XXX make this consistent with entry 0 == unmatched
     83    int nEntry = 1;
     84    psVector *fwhmValues = psMetadataLookupVector(&status, readout->analysis, "STACK.PSF.FWHM.VALUES");
     85    if (fwhmValues) {
     86        psAssert (entry < fwhmValues->n, "inconsistent matched-PSF entry");
     87        nEntry = fwhmValues->n;
     88    }
     89    if (entry > 0) {
     90        psLogMsg ("psphot", PS_LOG_DETAIL, "Radial Apertures for matched image %s : PSF FWHM = %f pixels\n", file->name, fwhmValues->data.F32[entry]);
     91    } else {
     92        psLogMsg ("psphot", PS_LOG_DETAIL, "Radial Apertures for unmatched image %s\n", file->name);
     93    }
    7494
    7595    // option to limit analysis to a specific region
    … …  
    98118            psArrayAdd(job->args, 1, AnalysisRegion);
    99119            psArrayAdd(job->args, 1, recipe);
     120
     121            PS_ARRAY_ADD_SCALAR(job->args, entry,  PS_TYPE_S32);
     122            PS_ARRAY_ADD_SCALAR(job->args, nEntry, PS_TYPE_S32);
    100123
    101124            PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nradial
    … …  
    115138            }
    116139            psScalar *scalar = NULL;
    117             scalar = job->args->data[4];
     140            scalar = job->args->data[6];
    118141            Nradial += scalar->data.S32;
    119142            psFree(job);
    … …  
    135158            } else {
    136159                psScalar *scalar = NULL;
    137                 scalar = job->args->data[4];
     160                scalar = job->args->data[6];
    138161                Nradial += scalar->data.S32;
    139162            }
    … …  
    159182    psMetadata *recipe      = job->args->data[3];
    160183
     184    int entry               = PS_SCALAR_VALUE(job->args->data[4],S32); // which psf-matched image are we working on? (0 == unmatched)
     185    int nEntry              = PS_SCALAR_VALUE(job->args->data[5],S32); // total number of psf-matched images + 1 unmatched
     186
    161187    // radMax stores the upper bounds of the annuli
    162188    // XXX keep the same name here as for the petrosian / elliptical apertures?
    … …  
    198224
    199225        // allocate pmSourceExtendedParameters, if not already defined
    200         if (!source->radialAper) {
    201             source->radialAper = psArrayAlloc(1);
     226        // XXX check that nPSFsizes is consistent with targets
     227        if (source->parent) {
     228            if (!source->parent->radialAper) {
     229                source->parent->radialAper = psArrayAlloc(nEntry);
     230            }
     231        } else {
     232            if (!source->radialAper) {
     233                source->radialAper = psArrayAlloc(nEntry);
     234            }
    202235        }
    203236
    … …  
    219252        pmSourceRedefinePixels (source, readout, source->peak->xf, source->peak->yf, outerRadius + 2);
    220253
    221         if (!psphotRadialApertureSource (source, recipe, maskVal, radMax, 0)) {
     254        if (!psphotRadialApertureSource (source, recipe, maskVal, radMax, entry)) {
    222255            psTrace ("psphot", 5, "failed to extract radial profile for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My);
    223256        } else {
    … …  
    233266        pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
    234267    }
    235     psScalar *scalar = job->args->data[4];
     268    psScalar *scalar = job->args->data[6];
    236269    scalar->data.S32 = Nradial;
    237270

  • branches/eam_branches/ipp-20110710/psphot/src/psphotReadout.c

    r32302 r32322  
    141141
    142142    // mark blended peaks PS_SOURCE_BLEND (detections->newSources)
     143    // XXX I've deactivated this because it was preventing galaxies close to stars from being
     144    // XXX fitted as an extended source.
    143145    if (false && !psphotBasicDeblend (config, view, filerule)) {
    144146        psError (PSPHOT_ERR_UNKNOWN, false, "failed on deblend analysis");
    … …  
    313315    psphotExtendedSourceAnalysis (config, view, filerule); // pass 1 (detections->allSources)
    314316    psphotExtendedSourceFits (config, view, filerule); // pass 1 (detections->allSources)
    315     psphotRadialApertures(config, view, filerule);
     317    psphotRadialApertures(config, view, filerule, 0);
    316318
    317319finish:

  • branches/eam_branches/ipp-20110710/psphot/src/psphotSetThreads.c

    r32299 r32322  
    5050    psFree(task);
    5151
    52     task = psThreadTaskAlloc("PSPHOT_RADIAL_APERTURES", 5);
     52    task = psThreadTaskAlloc("PSPHOT_RADIAL_APERTURES", 7);
    5353    task->function = &psphotRadialApertures_Threaded;
    5454    psThreadTaskAdd(task);

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

    r31154 r32322  
    4848    pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
    4949    psAssert (detections, "missing detections?");
    50     psAssert (detections->newSources, "new sources not defined?");
    51     psAssert (!detections->allSources, "all sources already defined?");
    52 
    53     // XXX TEST:
    54     if (detections->newSources) {
    55         psphotMatchSourcesToObjects(objects, detections->newSources, RADIUS);
    56     }
     50    psAssert (detections->allSources, "all sources not defined?");
     51
     52    psphotMatchSourcesToObjects(objects, detections->allSources, RADIUS);
    5753
    5854    return true;
    … …  
    261257            peak->assigned = true;
    262258            pmPhotObjAddSource(obj, source);
    263             psArrayAdd (detections->newSources, 100, source);
     259            psArrayAdd (detections->allSources, 100, source);
    264260            psFree (source);
    265261        }

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

    r30624 r32322  
    11# include "psphotInternal.h"
     2
     3// NOTE : element 0 of fwhmValues if the unmatched image, 
     4
     5int psphotStackMatchPSFsEntries (pmConfig *config, const pmFPAview *view, const char *filerule) {
     6
     7    int nRadialEntries = 0;
     8
     9    // find the currently selected readout
     10    pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, 0); // File of interest
     11    psAssert (file, "missing file?");
     12   
     13    pmReadout *readout = pmFPAviewThisReadout(view, file->fpa);
     14    psAssert (readout, "missing readout?");
     15   
     16    bool status = false;
     17    psVector *fwhmValues = psMetadataLookupVector(&status, readout->analysis, "STACK.PSF.FWHM.VALUES");
     18    if (!fwhmValues) {
     19        return 1;
     20    }
     21   
     22    nRadialEntries = fwhmValues->n;
     23    return nRadialEntries;
     24}
    225
    326// smooth the input image to match the next target PSF
    … …  
    528// and that the smoothing can use a 1D Gaussian kernel of width sqrt(TARGET^2 - CURRENT^2)
    629// each subsequent call
    7 bool psphotStackMatchPSFsNext(bool *smoothAgain, pmConfig *config, const pmFPAview *view, const char *filerule, int lastSize)
     30bool psphotStackMatchPSFsNext(pmConfig *config, const pmFPAview *view, const char *filerule, int lastSize)
    831{
    9     bool status = true;
    10 
    11     // select the appropriate recipe information
    12     psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
    13     psAssert (recipe, "missing recipe?");
    14 
    15     psVector *fwhmValues = psMetadataLookupVector(&status, recipe, "PSPHOT.STACK.TARGET.PSF.FWHM"); // Magnitude offsets
    16     if (!status) {
    17         // must not be a vector, only one value requested
    18         *smoothAgain = false;
    19         return true;
    20     }
    21 
    22     if (lastSize + 1 >= fwhmValues->n) {
    23         // all done with target FWHM values
    24         *smoothAgain = false;
    25         return true;
    26     }
    27 
    2832    int num = psphotFileruleCount(config, filerule);
    2933
    … …  
    3337    // loop over the available readouts
    3438    for (int i = 0; i < num; i++) {
    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]);
     39        if (!psphotStackMatchPSFsNextReadout (config, view, filerule, i, lastSize)) {
     40            psError (PSPHOT_ERR_CONFIG, false, "failed to smooth image %s (%d) to target PSF", filerule, i);
    3741            psImageConvolveSetThreads(oldThreads);
    3842            return false;
    … …  
    4448}
    4549
    46 bool psphotStackMatchPSFsNextReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe, psVector *fwhmValues, int lastSize) {
     50bool psphotStackMatchPSFsNextReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, int lastSize) {
    4751
    4852    bool status = false;
    … …  
    5862    psphotVisualShowImage(readout);
    5963
     64    // select the appropriate recipe information
     65    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
     66    psAssert (recipe, "missing recipe?");
     67
    6068    // user-defined masks to test for good/bad pixels (build from recipe list if not yet set)
    6169    psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels
    … …  
    6674        psWarning("PEAKS_MIN_GAUSS is not set in recipe; using default value");
    6775        minGauss = 0.5;
     76    }
     77
     78    psVector *fwhmValues = psMetadataLookupVector(&status, readout->analysis, "STACK.PSF.FWHM.VALUES");
     79    psAssert (fwhmValues, "need target PSFs");
     80
     81    if (lastSize + 1 >= fwhmValues->n) {
     82        return true;
    6883    }
    6984
    … …  
    128143    // psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "SIGNIFICANCE_SCALE_FACTOR", PS_META_REPLACE, "Signicance scale factor", factor);
    129144
    130     // save the output fwhm values in the readout->analysis.  we may have / will have multiple output PSF sizes,
    131     // so we save this in a vector.  if the vector is not yet defined, create it
    132 
    133     psVector *fwhmValuesOut = psMetadataLookupVector(&status, readout->analysis, "STACK.PSF.FWHM.VALUES");
    134     psAssert(fwhmValuesOut, "should already exist..");
    135     psVectorAppend(fwhmValuesOut, targetFWHM);
    136 
    137145    // do not generate a PSF if we already were supplied one
    138146    pmPSF *psfOld = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF");

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

    r31154 r32322  
    4444bool psphotStackReadout (pmConfig *config, const pmFPAview *view) {
    4545
     46    // measure the total elapsed time in psphotReadout.  dtime is the elapsed time used jointly
     47    // by the multiple threads, not the total time used by all threads.
    4648    psTimerStart ("psphotReadout");
    4749
    … …  
    5658    // optional break-point for processing
    5759    char *breakPt = psMetadataLookupStr (NULL, recipe, "BREAK_POINT");
    58     PS_ASSERT_PTR_NON_NULL (breakPt, false);
    59 
    60     // XXX or do I set OUT to be a pmFPAfile pointing at the input of interest?
     60    psAssert (breakPt, "configuration error: set BREAK_POINT");
     61
     62    // we have 3 relevant files: RAW (unconvolved), CNV (convolved stack), OUT (psf-matched stack)
     63    // select which image (RAW or CNV) is used for analysis (RAW always used for detection)
    6164    bool useRaw = psMetadataLookupBool (NULL, recipe, "PSPHOT.STACK.USE.RAW");
    6265    char *STACK_SRC = useRaw ? STACK_RAW : STACK_CNV;
    63     char *STACK_DET = STACK_RAW; // XXX optionally allow this to be CNV?
    64 
    65     // we have 3 relevant files: RAW, CNV, OUT
     66    char *STACK_DET = STACK_RAW;
    6667
    6768    // set the photcode for each image
    … …  
    7172    }
    7273
    73     // Generate the mask and weight images
    74     // XXX this should be done before we perform the convolutions
     74    // Generate the mask and weight images (if not supplied) and set mask bits
    7575    if (!psphotSetMaskAndVariance (config, view, STACK_DET)) {
    7676        return psphotReadoutCleanup (config, view, STACK_SRC);
    … …  
    8080    }
    8181
     82    // XXX I think this is not defined correctly for an array of images.
     83    // XXX I probably need to subtract the model (same model?) for both RAW and OUT.
     84    // XXX But, probably not a problem in practice since the stacks are constructed with 0.0 mean level.
     85
    8286    // generate a background model (median, smoothed image)
    83     // XXX I think this is not defined correctly for an array of images.
    84     // XXX probably need to subtract the model (same model?) for both RAW and OUT
    8587    if (!psphotModelBackground (config, view, STACK_DET)) {
    8688        return psphotReadoutCleanup (config, view, STACK_SRC);
    … …  
    9395    }
    9496
     97    // also make the chisq detection image
    9598    if (!psphotStackChisqImage(config, view, STACK_DET, STACK_SRC)) {
    9699        psError (PSPHOT_ERR_UNKNOWN, false, "failure to generate chisq image");
    … …  
    109112    }
    110113
    111     // copy the detections from DET to SRC
     114    // if DET and SRC are different images, copy the detections from DET to SRC
    112115    if (strcmp(STACK_SRC, STACK_DET)) {
    113116        if (!psphotCopySources (config, view, STACK_SRC, STACK_DET)) {
    … …  
    122125        return psphotReadoutCleanup (config, view, STACK_SRC);
    123126    }
    124 
    125     if (!strcasecmp (breakPt, "TEST1")) {
    126         return psphotReadoutCleanup (config, view, STACK_SRC);
    127     }
    128 
     127    if (!strcasecmp (breakPt, "PEAKS")) {
     128        return psphotReadoutCleanup (config, view, STACK_SRC);
     129    }
    129130    psMemDump("sourcestats");
    130131
    131     // generate the objects (object unify the sources from the different images)
     132    // classify sources based on moments, brightness
     133    // only run this on detections from the input images, not chisq image
     134    if (!psphotRoughClass (config, view, STACK_SRC)) {
     135        psError (PSPHOT_ERR_UNKNOWN, false, "failed to determine rough classifications");
     136        return psphotReadoutCleanup (config, view, STACK_SRC);
     137    }
     138
     139    // if we were not supplied a PSF model, determine the IQ stats here (detections->newSources)
     140    // only run this on detections from the input images, not chisq image
     141    if (!psphotImageQuality (config, view, STACK_SRC)) { // pass 1
     142        psError (PSPHOT_ERR_UNKNOWN, false, "failed to measure image quality");
     143        return psphotReadoutCleanup (config, view, STACK_SRC);
     144    }
     145    if (!strcasecmp (breakPt, "MOMENTS")) {
     146        return psphotReadoutCleanup (config, view, STACK_SRC);
     147    }
     148
     149    // use bright stellar objects to measure PSF
     150    if (!psphotChoosePSF (config, view, STACK_SRC, true)) { // pass 1
     151        psLogMsg ("psphot", 3, "failure to construct a psf model");
     152        return psphotReadoutCleanup (config, view, STACK_SRC);
     153    }
     154    if (!strcasecmp (breakPt, "PSFMODEL")) {
     155        return psphotReadoutCleanup (config, view, STACK_SRC);
     156    }
     157
     158    // construct an initial model for each object, set the radius to fitRadius, set circular fit mask
     159    psphotGuessModels (config, view, STACK_SRC);
     160
     161    // merge the newly selected sources into the existing list
     162    // NOTE: merge OLD and NEW
     163    psphotMergeSources (config, view, STACK_SRC);
     164
     165    // linear PSF fit to source peaks, subtract the models from the image (in PSF mask)
     166    // XXX why do this as a stack operation?
     167    // psphotFitSourcesLinearStack (config, objects, false);
     168    psphotFitSourcesLinear (config, view, STACK_SRC, false);
     169    psphotStackVisualFilerule(config, view, STACK_SRC);
     170
     171    // re-measure the kron mags with models subtracted.  this pass starts with a circular
     172    // window of size PSF_MOMENTS_RADIUS (same window used to measure the psf-scale moments)
     173    // but iterates to an appropriately larger size
     174    psphotKronIterate(config, view, STACK_SRC);
     175
     176    // identify CRs and extended sources
     177    psphotSourceSize (config, view, STACK_SRC, true);
     178
     179    // non-linear PSF and EXT fit to brighter sources
     180    // replace model flux, adjust mask as needed, fit, subtract the models (full stamp)
     181    psphotBlendFit (config, view, STACK_SRC); // pass 1 (detections->allSources)
     182
     183    // replace all sources
     184    psphotReplaceAllSources (config, view, STACK_SRC); // pass 1 (detections->allSources)
     185
     186    // linear fit to include all sources (subtract again)
     187    // NOTE : apply to ALL sources (extended + psf)
     188    psphotFitSourcesLinear (config, view, STACK_SRC, true); // pass 2 (detections->allSources)
     189
     190    // if we only do one pass, skip to extended source analysis
     191    if (!strcasecmp (breakPt, "PASS1")) goto pass1finish;
     192
     193    // NOTE: possibly re-measure background model here with objects subtracted / or masked
     194
     195    // NOTE: this block performs the 2nd pass low-significance PSF detection stage
     196    {
     197        // add noise for subtracted objects
     198        psphotAddNoise (config, view, STACK_DET); // pass 1 (detections->allSources)
     199
     200        // find fainter sources
     201        // NOTE: finds new peaks and new footprints, OLD and FULL set are saved on detections
     202        psphotFindDetections (config, view, STACK_DET, false); // pass 2 (detections->peaks, detections->footprints)
     203
     204        // remove noise for subtracted objects (ie, return to normal noise level)
     205        // NOTE: this needs to operate only on the OLD sources
     206        psphotSubNoise (config, view, STACK_DET); // pass 1 (detections->allSources)
     207
     208        // if DET and SRC are different images, copy the detections from DET to SRC
     209        if (strcmp(STACK_SRC, STACK_DET)) {
     210            // XXX how does this handle 1st vs 2nd pass sources?
     211            if (!psphotCopySources (config, view, STACK_SRC, STACK_DET)) {
     212                psError (PSPHOT_ERR_UNKNOWN, false, "failure in peak analysis");
     213                return psphotReadoutCleanup (config, view, STACK_SRC);
     214            }
     215        }
     216
     217        // define new sources based on only the new peaks & measure moments
     218        // NOTE: new sources are saved on detections->newSources
     219        psphotSourceStats (config, view, STACK_SRC, false); // pass 2 (detections->newSources)
     220
     221        // set source type
     222        // NOTE: apply only to detections->newSources
     223        if (!psphotRoughClass (config, view, STACK_SRC)) { // pass 2 (detections->newSources)
     224            psLogMsg ("psphot", 3, "failed to find a valid PSF clump for image");
     225            return psphotReadoutCleanup (config, view, STACK_SRC);
     226        }
     227
     228        // create full input models, set the radius to fitRadius, set circular fit mask
     229        // NOTE: apply only to detections->newSources
     230        psphotGuessModels (config, view, STACK_SRC); // pass 2 (detections->newSources)
     231
     232        // replace all sources so fit below applies to all at once
     233        // NOTE: apply only to OLD sources (which have been subtracted)
     234        psphotReplaceAllSources (config, view, STACK_SRC); // pass 2
     235
     236        // merge the newly selected sources into the existing list
     237        // NOTE: merge OLD and NEW
     238        // XXX check on free of sources...
     239        psphotMergeSources (config, view, STACK_SRC); // (detections->newSources + detections->allSources -> detections->allSources)
     240
     241        // NOTE: apply to ALL sources
     242        psphotFitSourcesLinear (config, view, STACK_SRC, true); // pass 3 (detections->allSources)
     243    }
     244
     245pass1finish:
     246
     247    // re-measure the kron mags with models subtracted
     248    // psphotKronMasked(config, view, STACK_SRC);
     249    psphotKronIterate(config, view, STACK_SRC);
     250
     251    // measure source size for the remaining sources
     252    // NOTE: applies only to NEW (unmeasured) sources
     253    psphotSourceSize (config, view, STACK_SRC, false); // pass 2 (detections->allSources)
     254
     255    psMemDump("psfstats");
     256
     257    // generate the objects (objects unify the sources from the different images)
    132258    // XXX this could just match the detections for the chisq image, and not bother measuring the
    133259    // source stats in that case...
    134260    psArray *objects = psphotMatchSources (config, view, STACK_SRC);
    135 
    136261    psMemDump("matchsources");
    137262
    138     if (!strcasecmp (breakPt, "TEST2")) {
    139         psFree(objects);
    140         return psphotReadoutCleanup (config, view, STACK_SRC);
    141     }
    142 
    143     // construct sources for the newly-generated sources (from other images)
    144     if (!psphotSourceStats (config, view, STACK_SRC, false)) { // pass 1
    145         psFree(objects);
    146         psError(PSPHOT_ERR_UNKNOWN, false, "failure to generate sources");
    147         return psphotReadoutCleanup (config, view, STACK_SRC);
    148     }
    149 
    150     psMemDump("sourcestats");
    151 
    152     // find blended neighbors of very saturated stars (detections->newSources)
    153     // if (!psphotDeblendSatstars (config, view)) {
    154     //     psError (PSPHOT_ERR_UNKNOWN, false, "failed on satstar deblend analysis");
    155     //     return psphotReadoutCleanup (config, view, STACK_SRC);
    156     // }
    157 
    158     // mark blended peaks PS_SOURCE_BLEND (detections->newSources)
    159     // if (!psphotBasicDeblend (config, view)) {
    160     //     psError (PSPHOT_ERR_UNKNOWN, false, "failed on deblend analysis");
    161     //     return psphotReadoutCleanup (config, view, STACK_SRC);
    162     // }
    163 
    164     // classify sources based on moments, brightness
    165     // only run this on detections from the input images, not chisq image
    166     if (!psphotRoughClass (config, view, STACK_SRC)) {
    167         psFree(objects);
    168         psError (PSPHOT_ERR_UNKNOWN, false, "failed to determine rough classifications");
    169         return psphotReadoutCleanup (config, view, STACK_SRC);
    170     }
    171     // if we were not supplied a PSF model, determine the IQ stats here (detections->newSources)
    172     // only run this on detections from the input images, not chisq image
    173     if (!psphotImageQuality (config, view, STACK_SRC)) { // pass 1
    174         psFree(objects);
    175         psError (PSPHOT_ERR_UNKNOWN, false, "failed to measure image quality");
    176         return psphotReadoutCleanup (config, view, STACK_SRC);
    177     }
    178     if (!strcasecmp (breakPt, "MOMENTS")) {
    179         psFree(objects);
    180         return psphotReadoutCleanup (config, view, STACK_SRC);
    181     }
    182 
    183     // use bright stellar objects to measure PSF if we were supplied a PSF for any input file,
    184     // this step is skipped
    185     if (!psphotChoosePSF (config, view, STACK_SRC, true)) { // pass 1
    186         psFree(objects);
    187         psLogMsg ("psphot", 3, "failure to construct a psf model");
    188         return psphotReadoutCleanup (config, view, STACK_SRC);
    189     }
    190     if (!strcasecmp (breakPt, "PSFMODEL")) {
    191         psFree(objects);
    192         return psphotReadoutCleanup (config, view, STACK_SRC);
    193     }
    194 
    195     // construct an initial model for each object, set the radius to fitRadius, set circular fit mask
    196     psphotGuessModels (config, view, STACK_SRC);
    197 
    198     // merge the newly selected sources into the existing list
    199     // NOTE: merge OLD and NEW
    200     psphotMergeSources (config, view, STACK_SRC);
    201 
    202     // linear PSF fit to source peaks, subtract the models from the image (in PSF mask)
    203     psphotFitSourcesLinearStack (config, objects, FALSE);
    204     psphotStackVisualFilerule(config, view, STACK_SRC);
    205 
    206     // identify CRs and extended sources
    207     psphotSourceSize (config, view, STACK_SRC, TRUE);
    208 
    209     // XXX do we want to do a preliminary (unconvolved) model fit here, and then
    210     // do a second detection pass? (like standard psphot)
    211 
    212     // measure aperture photometry corrections
    213     if (!psphotApResid (config, view, STACK_SRC)) {
    214         psFree (objects);
    215         psLogMsg ("psphot", 3, "failed on psphotApResid");
    216         return psphotReadoutCleanup (config, view, STACK_SRC);
    217     }
    218 
    219     psMemDump("psfstats");
    220 
    221263    psphotStackObjectsUnifyPosition (objects);
    222264
     265    psphotStackObjectsSelectForAnalysis (config, view, STACK_SRC, objects);
     266
    223267    // measure elliptical apertures, petrosians (objects sorted by S/N)
    224     psphotExtendedSourceAnalysisByObject (config, objects, view, STACK_SRC); // pass 1 (detections->allSources)
     268    // psphotExtendedSourceAnalysisByObject (config, objects, view, STACK_SRC); // pass 1 (detections->allSources)
     269    psphotExtendedSourceAnalysis (config, view, STACK_SRC); // pass 1 (detections->allSources)
    225270
    226271    // measure non-linear extended source models (exponential, deVaucouleur, Sersic) (sources sorted by S/N)
    227272    psphotExtendedSourceFits (config, view, STACK_SRC); // pass 1 (detections->allSources)
    228 
    229     // calculate source magnitudes
    230     psphotMagnitudes(config, view, STACK_SRC);
    231273
    232274    // create source children for the OUT filerule (for radial aperture photometry)
    … …  
    238280    }
    239281
     282    // measure circular, radial apertures (objects sorted by S/N)
     283    // XXX can we just use psphotRadialApertures
     284    // XXX make sure the headers are consistent with this (which PSF convolutions, ie mark 'none')
     285    // psphotRadialAperturesByObject (config, objectsRadial, view, STACK_SRC, nMatchedPSF);
     286    psphotRadialApertures (config, view, STACK_SRC, 0); // XXX entry 0 == unmatched?
    240287    psMemDump("extmeas");
    241288
    242     bool smoothAgain = true;
    243     for (int nMatchedPSF = 0; smoothAgain; nMatchedPSF++) {
     289    int nRadialEntries = psphotStackMatchPSFsEntries(config, view, STACK_OUT);
     290
     291    for (int entry = 1; entry < nRadialEntries; entry++) {
     292        // NOTE: entry 0 is the unmatched image set
    244293
    245294        // re-measure the PSF for the smoothed image (using entries in 'allSources')
    … …  
    253302
    254303        // measure circular, radial apertures (objects sorted by S/N)
    255         psphotRadialAperturesByObject (config, objectsRadial, view, STACK_OUT, nMatchedPSF);
     304        // entry 0 == unmatched? pass entry + 1?
     305        psphotRadialApertures (config, view, STACK_OUT, entry);
    256306
    257307        // replace the flux in the image so it is returned to its original state
    … …  
    259309
    260310        // smooth to the next FWHM, or set 'smoothAgain' to false if no more
    261         psphotStackMatchPSFsNext(&smoothAgain, config, view, STACK_OUT, nMatchedPSF);
     311        psphotStackMatchPSFsNext(config, view, STACK_OUT, entry);
    262312        psMemDump("matched");
    263313    }
    264314
    265     if (0 && !psphotEfficiency(config, view, STACK_OUT)) {
     315    // measure aperture photometry corrections
     316    if (!psphotApResid (config, view, STACK_SRC)) {
     317        psFree (objects);
     318        psLogMsg ("psphot", 3, "failed on psphotApResid");
     319        return psphotReadoutCleanup (config, view, STACK_SRC);
     320    }
     321
     322    // calculate source magnitudes
     323    psphotMagnitudes(config, view, STACK_SRC);
     324
     325    if (0 && !psphotEfficiency(config, view, STACK_DET)) {
    266326        psErrorStackPrint(stderr, "Unable to determine detection efficiencies from fake sources");
    267327        psErrorClear();
Note: See TracChangeset for help on using the changeset viewer.