IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 32348


Ignore:
Timestamp:
Sep 6, 2011, 1:32:31 PM (15 years ago)
Author:
eugene
Message:
  • clarified the output to show the stages and their timing more clearly
  • psphotStackReadout now uses the same functions as psphotReadout (deprecated FitLinearStack, ExtendedAnalysisByObject, RadialAnalysisByObject, etc)
  • psphotStack now does 2nd (5 sigma) detection pass and radial photometry of the un-matched image (as well as the matched images)
  • for speed, fitting uses a smaller subset of pixels than subtraction: the model is extended when the subtraction is performed to avoid a discontinuity
  • extended source analysis (petrosians) is threaded
  • iterative kron radius / magnitude analysis with down-weighted neighbors replaces the masked neighbor kron analysis. the result of this analysis is used by fitting guesses and is the value written to the PSF and XFIT tables.
  • the target matched PSFs are determined in advance and the 0 element of this array is always the unmatched image (with target PSF of NAN)
  • dQ, the difference between the upper and lower quartile, is measured (but not yet used to adjust the sky level. sky should be FITTED_MEAN - 2.5*dQ
  • fake sources inserted for the efficiency analysis are removed at the end of the analysis
  • radial apertures is now threaded
  • edge case petrosian mags are now handled without aborts (may still have a problem with the measured error of the petrosian radius)
  • new sersic model guess based on the psf size, the moments, and the kron flux. the psf - kron mag is used to guess the index, then the index is used to convert measured moments into model guess parameters. this uses empirical relationships between the quantities generated from fake images.
  • define the fit region and window region based on the 1st radial moments
  • use the kron S/N to determine if analysis is performed for extended sources, not the psf S/N
  • remove the old 'basic deblend' which did not really do any deblending and is superceeded by the deblending in KronIterate and the footprint analysis
Location:
trunk/psphot
Files:
50 edited
23 copied

Legend:

Unmodified
Added
Removed

  • trunk/psphot

  • trunk/psphot/src

    • Property svn:ignore

      •  

        old new  
        2222psphotMakePSF
        2323psphotStack
         24psphotModelTest

  • trunk/psphot/src/Makefile.am

    r31452 r32348  
    2525libpsphot_la_LDFLAGS = $(PSPHOT_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS)
    2626
    27 bin_PROGRAMS = psphot psphotForced psphotMakePSF psphotStack
     27bin_PROGRAMS = psphot psphotForced psphotMakePSF psphotStack psphotModelTest
    2828# bin_PROGRAMS = psphotPetrosianStudy psphotTest psphotMomentsStudy
    2929
    … …  
    4343psphotStack_LDFLAGS = $(PSPHOT_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS)
    4444psphotStack_LDADD = libpsphot.la
     45
     46psphotModelTest_CFLAGS = $(PSPHOT_CFLAGS) $(PSMODULE_CFLAGS) $(PSLIB_CFLAGS)
     47psphotModelTest_LDFLAGS = $(PSPHOT_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS)
     48psphotModelTest_LDADD = libpsphot.la
    4549
    4650# psphotMomentsStudy_CFLAGS = $(PSPHOT_CFLAGS) $(PSMODULE_CFLAGS) $(PSLIB_CFLAGS)
    … …  
    102106        psphotCleanup.c
    103107
     108# a psphot-variant that simply generates the PSF model
     109psphotModelTest_SOURCES = \
     110        psphotModelTest.c            \
     111        psphotModelTestArguments.c   \
     112        psphotParseCamera.c        \
     113        psphotImageLoop.c   \
     114        psphotMosaicChip.c         \
     115        psphotCleanup.c
     116
    104117# psphotTest_SOURCES = \
    105118#         psphotTest.c
    … …  
    117130        psphotVisual.c                 \
    118131        psphotCullPeaks.c              \
     132        psphotFootprintSaddles.c       \
    119133        psphotVersion.c                \
    120134        psphotModelGroupInit.c         \
    … …  
    129143        psphotMakePSFReadout.c         \
    130144        psphotModelBackground.c        \
     145        psphotModelTestReadout.c       \
    131146        psphotMaskBackground.c         \
    132147        psphotSubtractBackground.c     \
    … …  
    157172        psphotExtendedSourceAnalysisByObject.c \
    158173        psphotExtendedSourceFits.c     \
     174        psphotSersicModelClass.c       \
    159175        psphotKernelFromPSF.c          \
    160176        psphotFitSet.c                 \
    … …  
    168184        psphotRadialPlot.c             \
    169185        psphotKronMasked.c             \
     186        psphotKronIterate.c            \
    170187        psphotDeblendSatstars.c        \
    171188        psphotMosaicSubimage.c         \
    … …  
    191208        psphotEfficiency.c
    192209
    193 # XXX need to fix this for the new apis
    194 #       psphotModelTest.c             
    195 
    196210# re-instate these
    197211#       psphotIsophotal.c              \

  • trunk/psphot/src/psphot.c

    r31673 r32348  
    2727    }
    2828
    29     psLogMsg ("psphot", 3, "complete psphot run: %f sec\n", psTimerMark ("complete"));
     29    psLogMsg ("psphot", PS_LOG_WARN, "complete psphot run: %f sec\n", psTimerMark ("complete"));
    3030
    3131    psErrorCode exit_status = psphotGetExitStatus();

  • trunk/psphot/src/psphot.h

    r31452 r32348  
    1919    PSPHOT_FORCED,
    2020    PSPHOT_MAKE_PSF,
     21    PSPHOT_MODEL_TEST,
    2122} psphotImageLoopMode;
    2223
    … …  
    117118bool            psphotExtendedSourceAnalysis (pmConfig *config, const pmFPAview *view, const char *filerule);
    118119bool            psphotExtendedSourceAnalysisReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe);
     120bool            psphotExtendedSourceAnalysis_Threaded (psThreadJob *job);
    119121
    120122bool            psphotExtendedSourceFits (pmConfig *config, const pmFPAview *view, const char *filerule);
    … …  
    327329bool psphotMakePSFReadout(pmConfig *config, const pmFPAview *view, const char *filerule);
    328330
     331pmConfig *psphotModelTestArguments(int argc, char **argv);
     332bool psphotModelTestReadout(pmConfig *config, const pmFPAview *view, const char *filerule);
     333
    329334int psphotFileruleCount(const pmConfig *config, const char *filerule);
    330335bool psphotFileruleCountSet(const pmConfig *config, const char *filerule, int num);
    … …  
    372377    bool convolve;                      // Convolve images?
    373378    psphotStackConvolveSource convolveSource;
    374     // psArray *convImages, *convMasks, *convVariances; // Filenames for the temporary convolved images
    375 
    376     // bool matchZPs;                      // Adjust relative fluxes based on transparency analysis?
    377     // bool photometry;                    // Perform photometry?
    378     // psMetadata *stats;                  // Statistics for output
    379     // FILE *statsFile;                    // File to which to write statistics
    380     // psArray *origImages, *origMasks, *origVariances; // Filenames of the original images
    381     // psArray *origCovars;                // Original covariances matrices
    382     // int quality;                        // Bad data quality flag
    383379
    384380    // Prepare
    … …  
    387383    psVector *inputMask;                // Mask for inputs
    388384
    389     float targetSeeing;                 // Target seeing FWHM
     385    psVector *targetSeeing;             // Target seeing FWHMs
    390386    psArray *sourceLists;               // Individual lists of sources for matching
    391387    psVector *norm;                     // Normalisation for each image
    392388    psArray *psfs;
    393 
    394     // psVector *exposures;                // Exposure times
    395     // float sumExposure;                  // Sum of exposure times
    396     // float zp;                           // Zero point for output
    397     // psVector *inputMask;                // Mask for inputs
    398     // psArray *sources;                   // Matched sources
    399389
    400390    // Convolve
    … …  
    402392    psArray *regions;                   // PSF-matching regions --- required in the stacking
    403393    psVector *matchChi2;                // chi^2 for stamps from matching
    404     psVector *weightings;               // Combination weightings for images (1/noise^2)
    405     // psArray *cells;                     // Cells for convolved images --- a handle for reading again
    406     // int numCols, numRows;               // Size of image
    407     // psArray *convCovars;                // Convolved covariance matrices
    408 
    409     // Combine initial
    410     // pmReadout *outRO;                   // Output readout
    411     // pmReadout *expRO;                   // Exposure readout
    412     // psArray *inspect;                   // Array of arrays of pixels to inspect
    413 
    414     // Rejection
    415     // psArray *rejected;                  // Rejected pixels
    416394} psphotStackOptions;
    417395
    … …  
    420398bool psphotStackMatchPSFsReadout (pmConfig *config, const pmFPAview *view, psphotStackOptions *options, int index);
    421399bool psphotStackMatchPSFsPrepare (pmConfig *config, const pmFPAview *view, psphotStackOptions *options, int index);
    422 bool psphotStackMatchPSFsNext (bool *smoothAgain, pmConfig *config, const pmFPAview *view, const char *filerule, int lastSize);
    423 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);
    424404
    425405// psphotStackMatchPSFsUtils
    426 psVector *SetOptWidths (bool *optimum, psMetadata *recipe);
    427 pmReadout *makeFakeReadout(pmConfig *config, pmReadout *raw, psArray *sources, pmPSF *psf, psImageMaskType maskVal, int fullSize);
    428 bool rescaleData(pmReadout *readout, pmConfig *config, psphotStackOptions *options, int index);
    429 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);
    430410bool saveMatchData (pmReadout *readout, psphotStackOptions *options, int index);
    431411bool matchKernel(pmConfig *config, pmReadout *cnv, pmReadout *raw, psphotStackOptions *options, int index);
    432 bool dumpImageDiff(pmReadout *readoutConv, pmReadout *readoutFake, pmReadout *readoutRef, int index, char *rootname);
    433 bool dumpImage(pmReadout *readoutOut, pmReadout *readoutRef, int index, char *rootname);
    434 bool loadKernel (pmConfig *config, pmReadout *readoutCnv, psphotStackOptions *options, int index);
    435 
    436 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);
    437419
    438420psphotStackOptions *psphotStackOptionsAlloc (int num);
    … …  
    443425bool psphotCopySourcesReadout (pmConfig *config, const pmFPAview *view, const char *ruleOut, const char *ruleSrc, int index);
    444426
    445 bool psphotRadialApertures (pmConfig *config, const pmFPAview *view, const char *filerule);
    446 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);
     429bool psphotRadialApertures_Threaded (psThreadJob *job);
    447430bool psphotRadialApertureSource (pmSource *source, psMetadata *recipe, psImageMaskType maskVal, const psVector *radMax, int entry);
    448431
    … …  
    469452psArray *psphotSourceChildrenByObject (pmConfig *config, const pmFPAview *view, const char *filerule, psArray *objectsSrc);
    470453
     454bool psphotSersicModelClassGuessPCM (pmPCMdata *pcm, pmSource *source);
     455void psphotSersicModelClassInit ();
     456void psphotSersicModelClassCleanup ();
     457
     458bool psphotSetRadiusMoments (float *fitRadius, float *windowRadius, pmReadout *readout, pmSource *source, psImageMaskType markVal);
     459bool psphotSetRadiusMomentsExact (float *fitRadius, float *windowRadius, pmReadout *readout, pmSource *source, psImageMaskType markVal);
     460
     461bool psphotFootprintSaddles(pmReadout *readout, psArray *footprints);
     462bool psphotMaskFootprint (pmReadout *readout, pmSource *source, psImageMaskType markVal);
     463
     464bool psphotKronIterate (pmConfig *config, const pmFPAview *view, const char *filerule);
     465bool psphotKronIterateReadout(pmConfig *config, psMetadata *recipe, const pmFPAview *view, pmReadout *readout, psArray *sources, pmPSF *psf);
     466bool psphotKronIterate_Threaded (psThreadJob *job);
     467
     468bool psphotStackObjectsSelectForAnalysis (pmConfig *config, const pmFPAview *view, const char *filerule, psArray *objects);
     469
    471470#endif

  • trunk/psphot/src/psphotAddNoise.c

    r29936 r32348  
    8585    }
    8686    if (add) {
    87         psLogMsg ("psphot.noise", PS_LOG_INFO, "add noise for %ld objects: %f sec\n", sources->n, psTimerMark ("psphot.noise"));
     87        psLogMsg ("psphot.noise", PS_LOG_WARN, "add noise for %ld objects: %f sec\n", sources->n, psTimerMark ("psphot.noise"));
    8888    } else {
    89         psLogMsg ("psphot.noise", PS_LOG_INFO, "sub noise for %ld objects: %f sec\n", sources->n, psTimerMark ("psphot.noise"));
     89        psLogMsg ("psphot.noise", PS_LOG_WARN, "sub noise for %ld objects: %f sec\n", sources->n, psTimerMark ("psphot.noise"));
    9090    }
    9191

  • trunk/psphot/src/psphotApResid.c

    r31452 r32348  
    11# include "psphotInternal.h"
    2 # define DEBUG
     2// # define DEBUG
    33
    44# define SKIPSTAR(MSG) { psTrace ("psphot", 3, "invalid : %s", MSG); continue; }
    … …  
    99{
    1010    bool status = true;
     11
     12    fprintf (stdout, "\n");
     13    psLogMsg ("psphot", PS_LOG_INFO, "--- psphot Aperture Residuals ---");
    1114
    1215    // select the appropriate recipe information
    … …  
    347350
    348351    psLogMsg ("psphot.apresid", PS_LOG_DETAIL, "aperture residual: %f +/- %f\n", psf->ApResid, psf->dApResid);
    349     psLogMsg ("psphot.apresid", PS_LOG_INFO, "measure full-frame aperture residuals for %d sources: %f sec\n", Npsf, psTimerMark ("psphot.apresid"));
     352    psLogMsg ("psphot.apresid", PS_LOG_WARN, "measure full-frame aperture residuals for %d sources: %f sec\n", Npsf, psTimerMark ("psphot.apresid"));
    350353
    351354    psFree (xPos);

  • trunk/psphot/src/psphotBlendFit.c

    r31452 r32348  
    55{
    66    bool status = true;
     7
     8    fprintf (stdout, "\n");
     9    psLogMsg ("psphot", PS_LOG_INFO, "--- psphot Fit Sources (Non-Linear) ---");
    710
    811    // select the appropriate recipe information
    … …  
    186189    psFree (fitOptions);
    187190
    188     psLogMsg ("psphot.psphotBlendFit", PS_LOG_INFO, "fit models: %f sec for %d objects (%d psf, %d ext, %d failed, %ld skipped)\n", psTimerMark ("psphot.fit.nonlinear"), Nfit, Npsf, Next, Nfail, sources->n - Nfit);
     191    psLogMsg ("psphot.psphotBlendFit", PS_LOG_WARN, "fit models: %f sec for %d objects (%d psf, %d ext, %d failed, %ld skipped)\n", psTimerMark ("psphot.fit.nonlinear"), Nfit, Npsf, Next, Nfail, sources->n - Nfit);
    189192
    190193    psphotVisualShowResidualImage (readout, false);
    … …  
    233236        pmSource *source = sources->data[i];
    234237
     238# define TEST_X -420.0
     239# define TEST_Y 300.0
     240   
     241        if ((fabs(source->peak->xf - TEST_X) < 5) && (fabs(source->peak->yf - TEST_Y) < 5)) {
     242            fprintf (stderr, "test galaxy\n");
     243        }
     244
     245# undef TEST_X
     246# undef TEST_Y
     247
    235248        // skip non-astronomical objects (very likely defects)
    236249        if (source->mode &  PM_SOURCE_MODE_BLEND) continue;
    … …  
    246259
    247260        // limit selection to some SN limit
    248         if (sqrt(source->peak->detValue) < FIT_SN_LIM) continue;
    249 
     261        if (source->mode & PM_SOURCE_MODE_EXT_LIMIT) {
     262            if (source->moments->KronFlux < FIT_SN_LIM * source->moments->KronFluxErr) continue;
     263        } else {
     264            if (sqrt(source->peak->detValue) < FIT_SN_LIM) continue;
     265        }
    250266        // exclude sources outside optional analysis region
    251267        if (source->peak->xf < AnalysisRegion.x0) continue;
    … …  
    267283        }
    268284
    269         // replace object in image
     285        // replace object in image & remove excess noise
    270286        if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) {
    271287            pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
    … …  
    305321        Nfail ++;
    306322
    307         // re-subtract the object, leave local sky
     323        // re-subtract the object, leave local sky, re-bump noise
    308324        pmSourceCacheModel (source, maskVal);
    309325        pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);

  • trunk/psphot/src/psphotChoosePSF.c

    r31673 r32348  
    55{
    66    bool status = true;
     7
     8    fprintf (stdout, "\n");
     9    psLogMsg ("psphot", PS_LOG_INFO, "--- psphot Choose PSF ---");
    710
    811    // select the appropriate recipe information
    … …  
    376379
    377380    char *modelName = pmModelClassGetName (psf->type);
    378     psLogMsg ("psphot.pspsf", PS_LOG_INFO, "select psf model: %f sec\n", psTimerMark ("psphot.choose.psf"));
     381    psLogMsg ("psphot.pspsf", PS_LOG_WARN, "select psf model: %f sec\n", psTimerMark ("psphot.choose.psf"));
    379382    psLogMsg ("psphot.pspsf", PS_LOG_INFO, "psf model %s, ApResid: %f +/- %f\n", modelName, psf->ApResid, psf->dApResid);
    380383
    … …  
    443446            psVectorAppend (fwhmMinor, FWHM_MINOR);
    444447
    445             if (modelPSF->params->n >= 7) {
     448            if (modelPSF->params->n > 7) {
    446449              psVectorAppend (psfExtra1, modelPSF->params->data.F32[7]);
    447450            }
    448             if (modelPSF->params->n >= 8) {
     451            if (modelPSF->params->n > 8) {
    449452              psVectorAppend (psfExtra2, modelPSF->params->data.F32[8]);
    450453            }
    … …  
    517520    }
    518521
    519     // psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "ANGLE",    PS_META_REPLACE, "PSF angle",           axes.theta);
     522    psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "ANGLE",    PS_META_REPLACE, "PSF angle",           axes.theta);
    520523    psMetadataAddS32 (readout->analysis, PS_LIST_TAIL, "NPSFSTAR", PS_META_REPLACE, "Number of stars used to make PSF", psf->nPSFstars);
    521524

  • trunk/psphot/src/psphotEfficiency.c

    r31452 r32348  
    9999
    100100/// Generate a fake image and add it in to the existing readout
    101 static bool effGenerate(psImage **xSrc, psImage **ySrc, // Positions of sources
     101static pmReadout *effGenerate(psImage **xSrc, psImage **ySrc, // Positions of sources
    102102                        const pmReadout *ro,            // Readout of interest
    103103                        const pmPSF *psf,               // Point-spread function
    … …  
    152152        psFree(xAll);
    153153        psFree(yAll);
    154         return false;
     154        return NULL;
    155155    }
    156156    psFree(magAll);
    … …  
    161161
    162162    psBinaryOp(ro->image, ro->image, "+", fakeRO->image);
    163     psFree(fakeRO);
    164 
    165     return true;
     163
     164    // return the readout so we can subtract it later
     165    return fakeRO;
    166166}
    167167
    … …  
    169169{
    170170    bool status = true;
     171
     172    fprintf (stdout, "\n");
     173    psLogMsg ("psphot", PS_LOG_INFO, "--- psphot Efficiency ---");
    171174
    172175    // select the appropriate recipe information
    … …  
    290293
    291294    psImage *xFake = NULL, *yFake = NULL; // Coordinates of sources, each bin in a row
    292     if (!effGenerate(&xFake, &yFake, readout, psf, magOffsets,
    293                      numSources, magLim, radius, minFlux)) {
     295    pmReadout *fakeRO = effGenerate(&xFake, &yFake, readout, psf, magOffsets, numSources, magLim, radius, minFlux);
     296    if (!fakeRO) {
    294297        psError(PS_ERR_UNKNOWN, false, "Unable to generate fake sources");
    295298        psFree(xFake);
    … …  
    416419    psFree(significance);
    417420
     421    // psphotFitSourcesLinearReadout subtracts the model fits
    418422    if (!psphotFitSourcesLinearReadout(recipe, readout, fakeSourcesAll, psf, true)) {
    419423        psError(PS_ERR_UNKNOWN, false, "Unable to perform linear fit on fake sources.");
    … …  
    427431    psf->ApTrend = NULL;
    428432
     433    // measure the magnitudes and fluxes for the sources
    429434    if (!psphotMagnitudesReadout(config, recipe, view, readout, fakeSourcesAll, psf)) {
    430435        psError(PS_ERR_UNKNOWN, false, "Unable to measure magnitudes of fake sources.");
    … …  
    433438        psf->ApTrend = apTrend; // Casting away const!
    434439        return false;
     440    }
     441
     442    // replace the subtracted model fits
     443    for (int i = 0; i < fakeSourcesAll->n; i++) {
     444        pmSource *source = fakeSourcesAll->data[i];
     445
     446        // replace other sources?
     447        if (!(source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED)) continue;
     448       
     449        pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
    435450    }
    436451    psFree(fakeSourcesAll);
    … …  
    515530    psFree(fakeSources);
    516531
     532    // subtract the faked sources from the original image
     533    psBinaryOp(readout->image, readout->image, "-", fakeRO->image);
     534    psFree(fakeRO);
     535
    517536    pmDetEff *de = pmDetEffAlloc(magLim, numSources, numBins); // Detection efficiency
    518537    de->magOffsets = psVectorCopy(NULL, magOffsets, PS_TYPE_F32);
    … …  
    529548    psFree(de);
    530549
    531     psLogMsg("psphot", PS_LOG_INFO, "Detection efficiency: %lf sec\n", psTimerClear("psphot.fake"));
     550    psLogMsg("psphot", PS_LOG_WARN, "Detection efficiency: %lf sec\n", psTimerClear("psphot.fake"));
    532551
    533552    return true;

  • trunk/psphot/src/psphotEllipticalContour.c

    r29004 r32348  
    127127psF32 psphotEllipticalContourFunc (psVector *deriv, const psVector *params, const psVector *coord) {
    128128
    129     static int pass = 0;
    130 
    131129    psF32 *par = params->data.F32;
    132130
    … …  
    145143
    146144    // value is X
    147     // if (coord->data.F32[1] == 0) {
    148     if (pass == 0) {
    149         pass = 1;
    150 
     145    if (coord->data.F32[1] < 0.5) {
    151146        float value = par[PAR_RMIN]*cs_alpha*r;
    152147
    … …  
    161156
    162157    // value is Y
    163     // if (coord->data.F32[1] == 1) {
    164     if (pass == 1) {
    165         pass = 0;
    166 
     158    if (coord->data.F32[1] > 0.5) {
    167159        float value = par[PAR_RMIN]*sn_alpha*r;
    168160

  • trunk/psphot/src/psphotExtendedSourceAnalysis.c

    r31154 r32348  
    22
    33// measure the elliptical radial profile and use this to measure the petrosian parameters for the sources
    4 // XXX this function needs to be threaded
    54
    65// for now, let's store the detections on the readout->analysis for each readout
    … …  
    87{
    98    bool status = true;
     9
     10    fprintf (stdout, "\n");
     11    psLogMsg ("psphot", PS_LOG_INFO, "--- psphot Extended Source Analysis (Petrosians) ---");
    1012
    1113    // select the appropriate recipe information
    … …  
    5961    }
    6062
    61     // user-defined masks to test for good/bad pixels (build from recipe list if not yet set)
    62     psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels
    63     assert (maskVal);
     63    // determine the number of allowed threads
     64    int nThreads = psMetadataLookupS32(&status, config->arguments, "NTHREADS"); // Number of threads
     65    if (!status) {
     66        nThreads = 0;
     67    }
    6468
    6569    // get the sky noise from the background analysis; if this is missing, get the user-supplied value
    … …  
    7074    }
    7175
     76    // source analysis is done in S/N order (brightest first)
     77    sources = psArraySort (sources, pmSourceSortByFlux);
     78
     79    // option to limit analysis to a specific region
     80    char *region = psMetadataLookupStr (&status, recipe, "ANALYSIS_REGION");
     81    psRegion *AnalysisRegion = psRegionAlloc(0,0,0,0);
     82    *AnalysisRegion = psRegionForImage (readout->image, psRegionFromString (region));
     83    if (psRegionIsNaN (*AnalysisRegion)) psAbort("analysis region mis-defined");
     84
     85    // choose Cx, Cy (see psphotThreadTools.c for overview of the concepts)
     86    int Cx = 1, Cy = 1;
     87    psphotChooseCellSizes (&Cx, &Cy, readout, nThreads);
     88
     89    psArray *cellGroups = psphotAssignSources (Cx, Cy, sources);
     90
     91    for (int i = 0; i < cellGroups->n; i++) {
     92
     93        psArray *cells = cellGroups->data[i];
     94
     95        for (int j = 0; j < cells->n; j++) {
     96
     97            // allocate a job -- if threads are not defined, this just runs the job
     98            psThreadJob *job = psThreadJobAlloc ("PSPHOT_EXTENDED_ANALYSIS");
     99
     100            psArrayAdd(job->args, 1, readout);
     101            psArrayAdd(job->args, 1, cells->data[j]); // sources
     102            psArrayAdd(job->args, 1, AnalysisRegion);
     103            psArrayAdd(job->args, 1, recipe);
     104
     105            PS_ARRAY_ADD_SCALAR(job->args, skynoise, PS_TYPE_F32);
     106
     107            PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Next
     108            PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Npetro
     109            PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nannuli
     110
     111// set this to 0 to run without threading
     112# if (1)           
     113            if (!psThreadJobAddPending(job)) {
     114                psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
     115                psFree(AnalysisRegion);
     116                return false;
     117            }
     118# else
     119            if (!psphotExtendedSourceAnalysis_Threaded(job)) {
     120                psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
     121                psFree(AnalysisRegion);
     122                return false;
     123            }
     124            psScalar *scalar = NULL;
     125            scalar = job->args->data[5];
     126            Next += scalar->data.S32;
     127            scalar = job->args->data[6];
     128            Npetro += scalar->data.S32;
     129            scalar = job->args->data[7];
     130            Nannuli += scalar->data.S32;
     131            psFree(job);
     132# endif
     133        }
     134
     135        // wait for the threads to finish and manage results
     136        if (!psThreadPoolWait (false)) {
     137            psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
     138            psFree(AnalysisRegion);
     139            return false;
     140        }
     141
     142        // we have only supplied one type of job, so we can assume the types here
     143        psThreadJob *job = NULL;
     144        while ((job = psThreadJobGetDone()) != NULL) {
     145            if (job->args->n < 1) {
     146                fprintf (stderr, "error with job\n");
     147            } else {
     148                psScalar *scalar = NULL;
     149                scalar = job->args->data[5];
     150                Next += scalar->data.S32;
     151                scalar = job->args->data[6];
     152                Npetro += scalar->data.S32;
     153                scalar = job->args->data[7];
     154                Nannuli += scalar->data.S32;
     155            }
     156            psFree(job);
     157        }
     158    }
     159    psFree (cellGroups);
     160    psFree(AnalysisRegion);
     161
     162    psLogMsg ("psphot", PS_LOG_WARN, "extended source analysis: %f sec for %d objects\n", psTimerMark ("psphot.extended"), Next);
     163    psLogMsg ("psphot", PS_LOG_INFO, "  %d petrosian\n", Npetro);
     164    psLogMsg ("psphot", PS_LOG_INFO, "  %d annuli\n", Nannuli);
     165
     166    psphotVisualShowResidualImage (readout, false);
     167
     168    bool doPetrosian    = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_PETROSIAN");
     169    if (doPetrosian) {
     170        psphotVisualShowPetrosians (sources);
     171    }
     172
     173    return true;
     174}
     175
     176bool psphotExtendedSourceAnalysis_Threaded (psThreadJob *job) {
     177
     178    bool status;
     179
     180    int Next = 0;
     181    int Npetro = 0;
     182    int Nannuli = 0;
     183
     184    // arguments: readout, sources, models, region, psfSize, maskVal, markVal
     185    pmReadout *readout      = job->args->data[0];
     186    psArray *sources        = job->args->data[1];
     187    psRegion *region        = job->args->data[2];
     188    psMetadata *recipe      = job->args->data[3];
     189
     190    float skynoise          = PS_SCALAR_VALUE(job->args->data[4],F32);
     191
    72192    // S/N limit to perform full non-linear fits
    73193    float SN_LIM = psMetadataLookupF32 (&status, recipe, "EXTENDED_SOURCE_SN_LIM");
    … …  
    78198    bool doPetroStars   = psMetadataLookupBool (&status, recipe, "PETROSIAN_FOR_STARS");
    79199
    80     // source analysis is done in S/N order (brightest first)
    81     sources = psArraySort (sources, pmSourceSortByFlux);
    82 
    83     // option to limit analysis to a specific region
    84     char *region = psMetadataLookupStr (&status, recipe, "ANALYSIS_REGION");
    85     psRegion AnalysisRegion = psRegionForImage (readout->image, psRegionFromString (region));
    86     if (psRegionIsNaN (AnalysisRegion)) psAbort("analysis region mis-defined");
     200    // user-defined masks to test for good/bad pixels (build from recipe list if not yet set)
     201    psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels
     202    assert (maskVal);
    87203
    88204    // choose the sources of interest
    … …  
    90206
    91207        pmSource *source = sources->data[i];
     208
     209        // if we have checked the source validity on the basis of the object set, then
     210        // we either skip these tests below or we skip the source completely
     211        if (source->tmpFlags & PM_SOURCE_TMPF_STACK_SKIP) continue;
     212        if (source->tmpFlags & PM_SOURCE_TMPF_STACK_KEEP) goto keepSource;
    92213
    93214        // skip PSF-like and non-astronomical objects
    … …  
    104225
    105226        // limit selection to some SN limit
    106         assert (source->peak); // how can a source not have a peak?
    107         if (sqrt(source->peak->detValue) < SN_LIM) continue;
    108 
    109         // limit selection by analysis region
    110         if (source->peak->x < AnalysisRegion.x0) continue;
    111         if (source->peak->y < AnalysisRegion.y0) continue;
    112         if (source->peak->x > AnalysisRegion.x1) continue;
    113         if (source->peak->y > AnalysisRegion.y1) continue;
     227        // assert (source->peak); // how can a source not have a peak?
     228        // limit selection to some SN limit
     229        bool skipSource = false;
     230        if (source->mode & PM_SOURCE_MODE_EXT_LIMIT) {
     231            skipSource = (source->moments->KronFlux < SN_LIM * source->moments->KronFluxErr);
     232        } else {
     233            skipSource = (sqrt(source->peak->detValue) < SN_LIM);
     234        }
     235        if (skipSource) continue;
     236
     237        // limit selection by analysis region (this automatically apply
     238        if (source->peak->x < region->x0) continue;
     239        if (source->peak->y < region->y0) continue;
     240        if (source->peak->x > region->x1) continue;
     241        if (source->peak->y > region->y1) continue;
     242
     243    keepSource:
    114244
    115245        // replace object in image
    … …  
    159289    }
    160290
    161     psLogMsg ("psphot", PS_LOG_INFO, "extended source analysis: %f sec for %d objects\n", psTimerMark ("psphot.extended"), Next);
    162     psLogMsg ("psphot", PS_LOG_INFO, "  %d petrosian\n", Npetro);
    163     psLogMsg ("psphot", PS_LOG_INFO, "  %d annuli\n", Nannuli);
    164 
    165     psphotVisualShowResidualImage (readout, false);
    166 
    167     if (doPetrosian) {
    168         psphotVisualShowPetrosians (sources);
    169     }
    170 
    171     // fprintf (stderr, "xsrc : tried %ld objects\n", sources->n);
     291    psScalar *scalar = NULL;
     292
     293    // change the value of a scalar on the array (wrap this and put it in psArray.h)
     294    scalar = job->args->data[5];
     295    scalar->data.S32 = Next;
     296
     297    scalar = job->args->data[6];
     298    scalar->data.S32 = Npetro;
     299
     300    scalar = job->args->data[7];
     301    scalar->data.S32 = Nannuli;
    172302
    173303    return true;

  • trunk/psphot/src/psphotExtendedSourceFits.c

    r31673 r32348  
    66    bool status = true;
    77
     8    fprintf (stdout, "\n");
     9    psLogMsg ("psphot", PS_LOG_INFO, "--- psphot Extended Source Fits ---");
     10
    811    // select the appropriate recipe information
    912    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
    … …  
    1215    // perform full extended source non-linear fits?
    1316    if (!psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_FITS")) {
    14         psLogMsg ("psphot", PS_LOG_INFO, "skipping extended source measurements\n");
     17        psLogMsg ("psphot", PS_LOG_INFO, "skipping extended source fits\n");
    1518        return true;
    1619    }
    … …  
    4548    int Nfail = 0;
    4649
     50    psphotSersicModelClassInit();
     51
    4752    psTimerStart ("psphot.extended");
    4853
    … …  
    103108
    104109      if (item->type != PS_DATA_METADATA) {
     110        // XXX we could cull the bad entries or build a validated model folder
    105111        psAbort ("Invalid type for EXTENDED_SOURCE_MODEL entry %s, not a metadata folder", item->name);
    106         // XXX we could cull the bad entries or build a validated model folder
    107112      }
    108113
    … …  
    247252    psFree(AnalysisRegion);
    248253
    249     psLogMsg ("psphot", PS_LOG_INFO, "extended source model fits: %f sec for %d objects\n", psTimerMark ("psphot.extended"), Next);
     254    psphotSersicModelClassCleanup();
     255
     256    psphotVisualShowResidualImage (readout, false);
     257
     258    psLogMsg ("psphot", PS_LOG_WARN, "extended source model fits: %f sec for %d objects\n", psTimerMark ("psphot.extended"), Next);
    250259    psLogMsg ("psphot", PS_LOG_INFO, "  %d convolved models (%d passed)\n", Nconvolve, NconvolvePass);
    251260    psLogMsg ("psphot", PS_LOG_INFO, "  %d plain models (%d passed)\n", Nplain, NplainPass);
    … …  
    265274    int Nplain = 0;
    266275    int NplainPass = 0;
    267     bool savePics = false;
    268     float radius;
     276    float fitRadius, windowRadius;
    269277    psScalar *scalar = NULL;
    270     pmMoments psfMoments;
    271278
    272279    // arguments: readout, sources, models, region, psfSize, maskVal, markVal
    … …  
    275282    psMetadata *models      = job->args->data[2];
    276283    psRegion *region        = job->args->data[3];
    277     int psfSize             = PS_SCALAR_VALUE(job->args->data[4],PS_TYPE_IMAGE_MASK_DATA);
     284    int psfSize             = PS_SCALAR_VALUE(job->args->data[4],S32);
    278285    psImageMaskType maskVal = PS_SCALAR_VALUE(job->args->data[5],PS_TYPE_IMAGE_MASK_DATA);
    279286    psImageMaskType markVal = PS_SCALAR_VALUE(job->args->data[6],PS_TYPE_IMAGE_MASK_DATA);
    280 
    281     // pthread_t tid = pthread_self();     // Thread identifier
    282287
    283288    // Define source fitting parameters for extended source fits
    … …  
    299304
    300305        // skip PSF-like and non-astronomical objects
    301         if (source->type == PM_SOURCE_TYPE_STAR) continue;
     306        if (!(source->mode & PM_SOURCE_MODE_EXT_LIMIT)) continue;
    302307        if (source->type == PM_SOURCE_TYPE_DEFECT) continue;
    303308        if (source->type == PM_SOURCE_TYPE_SATURATED) continue;
    … …  
    309314        if (source->peak->y > region->y1) continue;
    310315
    311         // if model is NULL, we don't have a starting guess
    312         // XXX use the parameters guessed from moments
    313         // if (source->modelEXT == NULL) continue;
    314 
    315         // fprintf (stderr, "fit %d,%d in thread %d\n", source->peak->x, source->peak->y, (int) tid);
    316 
    317316        // replace object in image
    318317        if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) {
    … …  
    321320        Next ++;
    322321
    323         // set the radius based on the footprint (also sets the mask pixels)
    324         if (!psphotSetRadiusFootprint(&radius, readout, source, markVal, 1.0)) {
    325             fprintf (stderr, "skipping (1) %f, %f\n", source->peak->xf, source->peak->yf);
    326             pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
    327             // XXX raise an error of some kind
    328             continue;
    329         }
    330 
    331         // XXX note that this changes the source moments that are published...
    332         // recalculate the source moments using the larger extended-source moments radius
    333         // at this stage, skip Gaussian windowing, and do not clip pixels by S/N
    334         // this uses the footprint to judge both radius and aperture?
    335         // XXX save the psf-based moments for output
    336         psfMoments = *source->moments;
    337         if (!pmSourceMoments (source, radius, 0.0, 0.0, 0.0, maskVal)) {
    338             // subtract the best fit from the object, leave local sky
    339             fprintf (stderr, "skipping (2) %f, %f\n", source->peak->xf, source->peak->yf);
    340             *source->moments = psfMoments;
    341             pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
    342             // XXX raise an error flag of some kind
    343             continue;
    344         }
    345 
    346         // save the modelFlux here in case we need to subtract it (for failure)
    347         psImage *modelFluxStart = psMemIncrRefCounter (source->modelFlux);
    348         if (!modelFluxStart) {
    349             pmSourceCacheModel (source, maskVal);
    350             modelFluxStart = psMemIncrRefCounter (source->modelFlux);
    351             if (!modelFluxStart) {
    352                 fprintf (stderr, "skipping (3) %f, %f\n", source->peak->xf, source->peak->yf);
    353                 *source->moments = psfMoments;
    354                 pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
    355                 // XXX raise an error of some kind?
    356                 continue;
    357             }
    358         }
    359        
    360         if (savePics) {
    361           psphotSaveImage (NULL, readout->image, "image.xp.fits");
    362         }
    363 
    364         // array to store the pointers to the model flux images while the models are being fitted
    365         psArray *modelFluxes = psArrayAllocEmpty (models->list->n);
     322        // set the fit radius based on the first radial moment (also sets the mask pixels)
     323        psphotSetRadiusMomentsExact(&fitRadius, &windowRadius, readout, source, markVal);
     324
     325        // UPDATE : we have changed the moments calculation.  There is now an iteration within
     326        // psphotKronMasked to determine moments appropriate for a larger object.  The values
     327        // Mrf, KronFlux, and KronFluxErr are calculated for the iterated radius.  The other
     328        // values are left at the psf-based values.
    366329
    367330        // allocate the array to store the model fits
    … …  
    380343
    381344          // check the SNlim and skip model if source is too faint
    382           float SNlim = psMetadataLookupF32 (&status, model, "SNLIM_VALUE");
     345          float FIT_SN_LIM = psMetadataLookupF32 (&status, model, "SNLIM_VALUE");
    383346          assert (status);
    384347
    385348          // limit selection to some SN limit
    386349          // assert (source->peak); // how can a source not have a peak?
    387           if (sqrt(source->peak->detValue) < SNlim) {
     350          // limit selection to some SN limit
     351          bool skipSource = false;
     352          if (source->mode & PM_SOURCE_MODE_EXT_LIMIT) {
     353              skipSource = (source->moments->KronFlux < FIT_SN_LIM * source->moments->KronFluxErr);
     354          } else {
     355              skipSource = (sqrt(source->peak->detValue) < FIT_SN_LIM);
     356          }
     357          if (skipSource) {
    388358              Nfaint ++;
    389359              continue;
    … …  
    397367          bool convolved = psMetadataLookupBool (&status, model, "PSF_CONVOLVED_VALUE");
    398368          assert (status);
    399 
    400           // XXX psTraceSetLevel ("psLib.math.psMinimizeLMChi2", 6);
    401           // XXX psTraceSetLevel ("psphot.psphotModelWithPSF_LMM", 6);
    402369
    403370          // fit the model as convolved or not
    … …  
    410377                  continue;
    411378              }
    412               psTrace ("psphot", 4, "fit psf-conv model for %f, %f : %s chisq = %f\n", source->moments->Mx, source->moments->My, pmModelClassGetName (modelFit->type), modelFit->chisq);
     379              psTrace ("psphot", 4, "fit psf-conv model for %f, %f : %s chisq = %f (npix: %d, niter: %d)\n",
     380                       source->moments->Mx, source->moments->My, pmModelClassGetName (modelFit->type), modelFit->chisq, modelFit->nPix, modelFit->nIter);
    413381              Nconvolve ++;
    414382              if (!(modelFit->flags & (PM_MODEL_STATUS_BADARGS | PM_MODEL_STATUS_NONCONVERGE | PM_MODEL_STATUS_OFFIMAGE))) {
    … …  
    425393                  continue;
    426394              }
    427               pmSourceCacheModel (source, maskVal);
    428               psTrace ("psphot", 4, "fit plain model for %f, %f : %s chisq = %f\n", source->moments->Mx, source->moments->My, pmModelClassGetName (modelFit->type), modelFit->chisq);
     395              psTrace ("psphot", 4, "fit plain model for %f, %f : %s chisq = %f (npix: %d, niter: %d)\n",
     396                       source->moments->Mx, source->moments->My, pmModelClassGetName (modelFit->type), modelFit->chisq, modelFit->nPix, modelFit->nIter);
    429397              Nplain ++;
    430398              if (!(modelFit->flags & (PM_MODEL_STATUS_BADARGS | PM_MODEL_STATUS_NONCONVERGE | PM_MODEL_STATUS_OFFIMAGE))) {
    … …  
    434402          }
    435403
    436           // save each of the model flux images and store the best
    437           psArrayAdd (modelFluxes, 4, source->modelFlux);
     404          // XXX deprecate?
     405          // XXX pmSourceExtFitPars *extFitPars = pmSourceExtFitParsAlloc();
     406          // XXX extFitPars->Mxx = source->moments->Mxx;
     407          // XXX extFitPars->Mxy = source->moments->Mxy;
     408          // XXX extFitPars->Myy = source->moments->Myy;
     409          // XXX extFitPars->Mrf = source->moments->Mrf;
     410          // XXX extFitPars->Mrh = source->moments->Mrh;
     411          // XXX extFitPars->peakMag = (source->peak->rawFlux > 0) ? -2.5*log10(source->peak->rawFlux) : NAN;
     412         
     413          // save kron mag, but assign apMag & psfMag on output (not yet calculated)
     414          // XXX extFitPars->krMag = -2.5*log10(source->moments->KronFlux);
     415
     416          // psArrayAdd (source->extFitPars, 4, extFitPars);
     417          // psFree (extFitPars);
    438418
    439419          // test for fit quality / result
    440           modelFit->fitRadius = radius;
     420          modelFit->fitRadius = fitRadius;
    441421          psArrayAdd (source->modelFits, 4, modelFit);
    442422
    … …  
    463443        }
    464444
     445        // clear the circular mask
     446        psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal));
     447
    465448        if (minModel == -1) {
    466           // re-subtract the object, leave local sky
     449          // no valid extended fit; re-subtract the object, leave local sky
    467450          psTrace ("psphot", 5, "failed to fit extended source model to object at %f, %f", source->moments->Mx, source->moments->My);
    468451
    469           // replace original model, subtract it
    470           psFree (source->modelFlux);
    471           source->modelFlux = modelFluxStart;
    472 
    473           *source->moments = psfMoments;
     452          // ensure the modelEXT is cached
     453          pmSourceCacheModel (source, maskVal);
     454
    474455          pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
    475456
    476           psFree (modelFluxes);
    477 
    478           if (savePics) {
    479               psphotSaveImage (NULL, readout->image, "image.xp.fits");
    480               char key[10];
    481               fprintf (stdout, "continue? ");
    482               if (!fgets (key, 8, stdin)) {
    483                   psWarning("Couldn't read anything.");
    484               } else if (key[0] == 'n') {
    485                   savePics = false;
    486               }
    487           }
    488457          continue;
    489458        }
    … …  
    493462        source->modelEXT = psMemIncrRefCounter (source->modelFits->data[minModel]);
    494463
    495         // save the modelFlux for the best model
    496         psFree (source->modelFlux);
    497         source->modelFlux = psMemIncrRefCounter (modelFluxes->data[minModel]);
    498 
    499         // replace the original moments
    500         *source->moments = psfMoments;
     464        // adjust the window so the subtraction covers the faint wings
     465        psphotSetRadiusMoments(&fitRadius, &windowRadius, readout, source, markVal);
     466
     467        // cache the model flux
     468        if (source->modelEXT->isPCM) {
     469            pmPCMCacheModel (source, maskVal, psfSize);
     470        } else {
     471            pmSourceCacheModel (source, maskVal);
     472        }
    501473
    502474        // subtract the best fit from the object, leave local sky
    503475        pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
    504476
    505         // the initial model flux is no longer needed
    506         psFree (modelFluxStart);
    507         psFree (modelFluxes);
    508 
    509477        psTrace ("psphot", 4, "best ext model for %f, %f : %s chisq = %f\n", source->moments->Mx, source->moments->My, pmModelClassGetName (source->modelEXT->type), source->modelEXT->chisq);
    510478        psTrace ("psphot", 5, "extended source model for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My);
    511 
    512         if (savePics) {
    513           psphotSaveImage (NULL, readout->image, "image.xm.fits");
    514           char key[10];
    515           fprintf (stdout, "continue? ");
    516           if (!fgets (key, 8, stdin)) {
    517               psWarning("Couldn't read anything.");
    518           } else if (key[0] == 'n') {
    519               savePics = false;
    520           }
    521         }
    522479    }
    523480    psFree (fitOptions);
    524 
    525     // fprintf (stderr, "xfit : tried %ld objects\n", sources->n);
    526481
    527482    // change the value of a scalar on the array (wrap this and put it in psArray.h)

  • trunk/psphot/src/psphotFindDetections.c

    r31154 r32348  
    77{
    88    bool status = true;
     9
     10    fprintf (stdout, "\n");
     11    psLogMsg ("psphot", PS_LOG_INFO, "--- psphot Find Detections ---");
    912
    1013    // select the appropriate recipe information

  • trunk/psphot/src/psphotFindFootprints.c

    r31154 r32348  
    5858    psphotCullPeaks(readout, significance, recipe, detections->footprints);
    5959    detections->peaks = pmFootprintArrayToPeaks(detections->footprints);
    60     psLogMsg ("psphot", PS_LOG_INFO, "%ld peaks, %ld total footprints: %f sec\n", detections->peaks->n, detections->footprints->n, psTimerMark ("psphot.footprints"));
     60
     61    // psphotFootprintSaddles (readout, detections->footprints);
     62    //
     63    // int nSaddle = 0;
     64    // for (int i = 0; i < detections->peaks->n; i++) {
     65    //  pmPeak *peak = detections->peaks->data[i];
     66    // 
     67    //  if (peak->saddlePoints) nSaddle += peak->saddlePoints->n;
     68    // }
     69    // fprintf (stderr, "%d saddle points for %ld peaks\n", nSaddle, detections->peaks->n);
     70
     71    psLogMsg ("psphot", PS_LOG_WARN, "%ld peaks, %ld total footprints: %f sec\n", detections->peaks->n, detections->footprints->n, psTimerMark ("psphot.footprints"));
    6172
    6273    return detections;

  • trunk/psphot/src/psphotFindPeaks.c

    r31154 r32348  
    6868        pmPeaksWriteText (peaks, output);
    6969    }
    70     psLogMsg ("psphot", PS_LOG_INFO, "%ld peaks: %f sec\n", peaks->n, psTimerMark ("psphot.peaks"));
     70    psLogMsg ("psphot", PS_LOG_WARN, "%ld peaks: %f sec\n", peaks->n, psTimerMark ("psphot.peaks"));
    7171
    7272    return peaks;

  • trunk/psphot/src/psphotFitSourcesLinear.c

    r31452 r32348  
    1717    bool status = true;
    1818
     19    fprintf (stdout, "\n");
     20    psLogMsg ("psphot", PS_LOG_INFO, "--- psphot Fit Source (Linear) ---");
     21
    1922    // select the appropriate recipe information
    2023    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
    … …  
    2326    int num = psphotFileruleCount(config, filerule);
    2427
     28    // skip the chisq image (optionally?)
     29    int chisqNum = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.CHISQ.NUM");
     30    if (!status) chisqNum = -1;
     31
    2532    // loop over the available readouts
    2633    for (int i = 0; i < num; i++) {
     34        if (i == chisqNum) continue; // skip chisq image
    2735
    2836        // find the currently selected readout
    … …  
    190198        pmModel *model = pmSourceGetModel (&isPSF, source);
    191199
     200        // clear the 'mark' pixels and remask on the fit aperture
    192201        psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal));
    193202        psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, model->fitRadius, "OR", markVal);
    … …  
    349358        model->dparams->data.F32[PM_PAR_I0] = errors->data.F32[i];
    350359
    351         // subtract object
     360        // clear the 'mark' pixels so the subtraction covers the full window
     361        psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal));
     362
     363        // subtract object & add noise
    352364        pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
    353365    }
    … …  
    374386    psFree (border);
    375387
    376     psLogMsg ("psphot.ensemble", PS_LOG_INFO, "measure ensemble of PSFs: %f sec\n", psTimerMark ("psphot.linear"));
     388    psLogMsg ("psphot.ensemble", PS_LOG_WARN, "measure ensemble of PSFs: %f sec\n", psTimerMark ("psphot.linear"));
    377389
    378390    psphotVisualPlotChisq (sources);

  • trunk/psphot/src/psphotFitSourcesLinearStack.c

    r31452 r32348  
    6464          // turn this bit off and turn it on again if we keep this source
    6565          source->mode &= ~PM_SOURCE_MODE_LINEAR_FIT;
     66
     67          // skip non-astronomical objects (very likely defects)
     68          if (source->type == PM_SOURCE_TYPE_DEFECT) continue;
     69          if (source->type == PM_SOURCE_TYPE_SATURATED) continue;
     70
     71          // do not include CRs in the full ensemble fit
     72          if (source->mode & PM_SOURCE_MODE_CR_LIMIT) continue;
     73
     74          // do not include MOMENTS_FAILURES in the fit
     75          if (source->mode & PM_SOURCE_MODE_MOMENTS_FAILURE) continue;
     76
     77          if (final) {
     78              if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) continue;
     79          } else {
     80              // if (source->mode & PM_SOURCE_MODE_BLEND) continue;
     81          }
    6682
    6783          // generate model for sources without, or skip if we can't
    … …  
    8399          }
    84100
    85           source->mode |= PM_SOURCE_MODE_LINEAR_FIT;
     101          bool isPSF = false;
     102          pmModel *model = pmSourceGetModel (&isPSF, source);
     103
     104          // clear the 'mark' pixels and remask on the fit aperture
     105          psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal));
     106          psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, model->fitRadius, "OR", markVal);
     107
     108          // we call this function multiple times. for the first time, we have only PSF models for all objects
     109          // the second time has extended sources.  If we ever fit the PSF model, we should raise this bit
     110          source->mode |= PM_SOURCE_MODE_LINEAR_FIT;
     111          if (isPSF) {
     112              source->mode |= PM_SOURCE_MODE_PSFMODEL;
     113          }         
    86114          psArrayAdd (fitSources, 100, source);
    87115        }
    … …  
    166194        model->params->data.F32[PM_PAR_I0] = norm->data.F32[i];
    167195        model->dparams->data.F32[PM_PAR_I0] = errors->data.F32[i];
     196
     197        // clear the 'mark' pixels so the subtraction covers the full window
     198        psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal));
    168199
    169200        // subtract object
    … …  
    190221    psFree (errors);
    191222
    192     psLogMsg ("psphot.ensemble", PS_LOG_INFO, "measure ensemble of PSFs: %f sec\n", psTimerMark ("psphot.linear"));
     223    psLogMsg ("psphot.ensemble", PS_LOG_WARN, "measure ensemble of PSFs: %f sec\n", psTimerMark ("psphot.linear"));
    193224    return true;
    194225}

  • trunk/psphot/src/psphotGuessModels.c

    r31452 r32348  
    141141    psFree (cellGroups);
    142142
    143     psLogMsg ("psphot.models", 4, "built models for %ld objects: %f sec\n", sources->n, psTimerMark ("psphot.models"));
     143    psLogMsg ("psphot.models", PS_LOG_WARN, "built models for %ld objects: %f sec\n", sources->n, psTimerMark ("psphot.models"));
    144144    return true;
    145145}

  • trunk/psphot/src/psphotImageLoop.c

    r31154 r32348  
    129129                    }
    130130                    break;
     131                  case PSPHOT_MODEL_TEST:
     132                    if (!psphotModelTestReadout (config, view, "PSPHOT.INPUT")) {
     133                        psError(psErrorCodeLast(), false, "failure in psphotReadout for chip %d, cell %d, readout %d\n", view->chip, view->cell, view->readout);
     134                        psFree (view);
     135                        return false;
     136                    }
     137                    break;
    131138                }
    132139            }

  • trunk/psphot/src/psphotImageQuality.c

    r30624 r32348  
    289289#endif
    290290
    291     psLogMsg ("psphot", PS_LOG_INFO, "Image Quality Stats from %ld psf stars : FWHM (major, minor) [pixels]: %f, %f\n",
     291    psLogMsg ("psphot", PS_LOG_WARN, "Image Quality Stats from %ld psf stars : FWHM (major, minor) [pixels]: %f, %f\n",
    292292              M2->n, fwhm_major, fwhm_minor);
    293293

  • trunk/psphot/src/psphotKronMasked.c

    r31673 r32348  
    11# include "psphotInternal.h"
    2 
     2# ifndef ROUND
     3# define ROUND(X) ((int) ((X) + 0.5*SIGN(X)))
     4# endif
     5
     6bool psphotKronMask (pmSource *source, psImage *kronMask, psImageMaskType markVal);
    37bool psphotKronMag (pmSource *source, float radius, float minKronRadius, psImageMaskType maskVal);
    48
    … …  
    4549}
    4650
     51int psphotKapaChannel (int channel);
     52bool psphotVisualShowMask (int kapaFD, psImage *inImage, const char *name, int channel);
     53bool psphotVisualRangeImage (int kapaFD, psImage *inImage, const char *name, int channel, float min, float max);
     54
    4755bool psphotKronMaskedReadout(pmConfig *config, psMetadata *recipe, const pmFPAview *view, pmReadout *readout, psArray *sources, pmPSF *psf) {
    4856
    … …  
    7280    }
    7381
     82    float EXT_FIT_MAX_RADIUS = psMetadataLookupF32 (&status, recipe, "EXT_FIT_MAX_RADIUS");
     83
    7484    // bit-masks to test for good/bad pixels
    7585    psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT");
    … …  
    8393    maskVal |= markVal;
    8494
     95    // psphotSaveImage (NULL, readout->mask, "kron.unmasked.fits");
     96
     97    pmSourcePhotometryMode photMode = PM_SOURCE_PHOT_PSFONLY;
     98
     99    // XXX tmp visualization
     100    // int kapa = psphotKapaChannel (1);
     101
    85102    // generate the mask image: increment counter for every source overlapping the pixel
    86     psImage *kronMask = psImageAlloc (readout->image->numCols, readout->image->numRows, PS_TYPE_U8);
     103    psImage *kronMask = psImageAlloc (readout->image->numCols, readout->image->numRows, PS_TYPE_S32);
    87104    psImageInit (kronMask, 0);
    88     int Nx = kronMask->numCols;
    89     int Ny = kronMask->numRows;
    90105    for (int i = 0; i < sources->n; i++) {
    91106
    … …  
    96111        if (!source->moments) continue;
    97112
    98         // replace object in image
    99         if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) {
    100             pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
    101         }
    102 
    103         psphotKronMag (source, RADIUS, MIN_KRON_RADIUS, maskVal);
    104 
    105         // re-subtract the object, leave local sky
    106         pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
    107 
    108         continue;
    109 
    110         // XXX skip this code
    111         // generate the pixel masks
    112         // int Xo = source->moments->Mx;
    113         // int Yo = source->moments->My;
    114         float dX = source->moments->Mx - source->peak->xf;
    115         float dY = source->moments->My - source->peak->yf;
    116         float dR = hypot(dX, dY);
    117        
    118         float Xo = (dR < 2.0) ? source->moments->Mx : source->peak->xf;
    119         float Yo = (dR < 2.0) ? source->moments->My : source->peak->yf;
    120        
    121         int Kr = 2.5*source->moments->Mrf;
    122         int Kr2 = Kr*Kr;
    123        
    124         for (int iy = Yo - Kr; iy < Yo + Kr + 1; iy++) {
    125             if (iy < 0) continue;
    126             if (iy >= Ny) continue;
    127             for (int ix = Xo - Kr; ix < Xo + Kr + 1; ix++) {
    128                 if (ix < 0) continue;
    129                 if (ix >= Nx) continue;
    130 
    131                 if (PS_SQR(ix - Xo) + PS_SQR(iy - Yo) > Kr2) continue;
    132                 if (kronMask->data.U8[iy][ix] < 0xff) {
    133                     kronMask->data.U8[iy][ix] ++;
    134                 }
    135             }
    136         }
    137     }           
     113        // replace object in image
     114        if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) {
     115            pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
     116        }
     117
     118        // save the PSF-based values
     119        // source->moments->KronFluxPSF    = source->moments->KronFlux;
     120        // source->moments->KronFluxPSFErr = source->moments->KronFluxErr;
     121        // source->moments->KronRadiusPSF  = source->moments->Mrf;
     122
     123        // iterate to the window radius
     124        float windowRadius = RADIUS;
     125        for (int j = 0; j < 4; j++) {
     126            // XXX use some S/N criterion to limit?
     127            // if (source->moments->KronFlux / source->moments->KronFluxErr > 10.0) {
     128            // windowRadius = PS_MAX (RADIUS, 10*source->moments->Mrf);
     129            // }
     130
     131            // re-allocate image, weight, mask arrays for each peak with box big enough to fit BIG_RADIUS
     132            pmSourceRedefinePixels (source, readout, source->peak->x, source->peak->y, windowRadius + 2);
     133
     134            // mask the pixels not contained by the footprint
     135            // if (psphotMaskFootprint (readout, source, markVal)) {
     136            //  // psphotVisualShowMask (kapa, source->maskObj, "mask", 2);
     137            //  // psphotVisualRangeImage (kapa, source->pixels, "image", 1, -50, 300);
     138            //  // fprintf (stderr, "masked\n");
     139            // }           
     140
     141            // mask the pixels associated with near neighbors
     142            if (psphotKronMask (source, kronMask, markVal)) {
     143                // psphotVisualShowMask (kapa, source->maskObj, "mask", 2);
     144                // psphotVisualRangeImage (kapa, source->pixels, "image", 1, -50, 300);
     145                // fprintf (stderr, "masked\n");
     146            }       
     147            // this function populates moments->Mrf,KronFlux,KronFluxErr
     148            // XXX what about KronInner, KronOuter, etc?
     149            psphotKronMag (source, windowRadius, MIN_KRON_RADIUS, maskVal);
     150            windowRadius = PS_MIN(PS_MAX(RADIUS, 4.0*source->moments->Mrf), EXT_FIT_MAX_RADIUS);
     151            psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal));
     152        }
     153        pmSourceMagnitudes (source, psf, photMode, maskVal, markVal, source->apRadius);
     154        float kmag = -2.5*log10(source->moments->KronFlux);
     155        if (source->psfMag - kmag > 0.25) {
     156            // psphotVisualShowMask (kapa, source->maskObj, "mask", 1);
     157            // psphotVisualRangeImage (kapa, source->pixels, "image", 0, -50, 300);
     158            // fprintf (stderr, "psf: %f, kron: %f, dmag: %f\n", source->psfMag, kmag, source->psfMag - kmag);
     159            // fprintf (stderr, "continue\n");
     160        }
     161        // re-subtract the object, leave local sky
     162        pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
     163    }
     164
     165    // psphotSaveImage (NULL, readout->mask, "kron.masked.fits");
    138166    // psphotSaveImage (NULL, kronMask, "kronmask.fits");
     167
    139168    psLogMsg ("psphot.kron", PS_LOG_DETAIL, "measure masked kron magnitudes : %f sec for %ld objects\n", psTimerMark ("psphot.kron"), sources->n);
    140169
    … …  
    142171    return true;
    143172}
     173
     174# define WEIGHTED 0
    144175
    145176bool psphotKronMag (pmSource *source, float radius, float minKronRadius, psImageMaskType maskVal) {
    … …  
    150181    PS_ASSERT_FLOAT_LARGER_THAN(radius, 0.0, false);
    151182
    152     psF32 Sum = 0.0;
    153     psF32 Var = 0.0;
    154183    psF32 R2 = PS_SQR(radius);
     184
     185# if (WEIGHTED)
     186    float rsigma2 = 0.5 / PS_SQR(radius / 4.0); // use a Gaussian window with sigma = R_window / 2
     187# endif
    155188
    156189    // a note about coordinates: coordinates of objects throughout psphot refer to the primary
    … …  
    168201    psF32 RS = 0.0;
    169202
    170 # if (1)
     203    // the peak position is less accurate but less subject to extreme deviations
    171204    float dX = source->moments->Mx - source->peak->xf;
    172205    float dY = source->moments->My - source->peak->yf;
    … …  
    174207    float Xo = (dR < 2.0) ? source->moments->Mx : source->peak->xf;
    175208    float Yo = (dR < 2.0) ? source->moments->My : source->peak->yf;
    176 # else
    177     float Xo = source->moments->Mx;
    178     float Yo = source->moments->My;
    179 # endif
    180209
    181210    // center of mass in subimage.  Note: the calculation below uses pixel index, so we correct
    … …  
    210239            if (r2 > R2) continue;
    211240
     241# if (WEIGHTED)
     242            float z = r2 * rsigma2;
     243            assert (z >= 0.0);
     244            float weight  = exp(-z);
     245            psF32 pDiff = *vPix * weight;
     246# else
    212247            psF32 pDiff = *vPix;
    213 
    214             Sum += pDiff;
    215 
    216             // Kron Flux uses the 1st radial moment (NOT Gaussian windowed)
     248# endif
     249
     250            // Kron Flux uses the 1st radial moment (maybe Gaussian windowed?)
    217251            psF32 rf = pDiff * sqrt(r2);
    218252            psF32 rs = pDiff;
    … …  
    234268
    235269    int nKronPix = 0;
    236     Sum = Var = 0.0;
     270    float Sum = 0.0;
     271    float Var = 0.0;
    237272
    238273    for (psS32 row = 0; row < source->pixels->numRows ; row++) {
    … …  
    263298            if (r2 > radKron2) continue;
    264299
     300# if (WEIGHTED)
     301            float z = r2 * rsigma2;
     302            assert (z >= 0.0);
     303            float weight  = exp(-z);
     304            psF32 pDiff = *vPix * weight;
     305            psF32 wDiff = *vWgt * weight;
     306# else
    265307            psF32 pDiff = *vPix;
    266308            psF32 wDiff = *vWgt;
     309# endif
    267310
    268311            Sum += pDiff;
    … …  
    272315    }
    273316
    274     source->moments->Mrh = Mrf;
    275 
    276     // XXX for a test, save the old values here:
    277     source->moments->KronFouter    = source->moments->KronFlux;
    278     source->moments->KronCoreErr = source->moments->KronFluxErr;
    279 
    280     source->moments->KronFlux    = Sum       * M_PI * PS_SQR(radKron) / nKronPix;
    281     source->moments->KronFluxErr = sqrt(Var) * M_PI * PS_SQR(radKron) / nKronPix;
     317    source->moments->Mrf = Mrf;
     318    source->moments->KronFlux    = Sum;
     319    source->moments->KronFluxErr = sqrt(Var);
     320    // source->moments->KronFlux    = Sum       * M_PI * PS_SQR(radKron) / nKronPix;
     321    // source->moments->KronFluxErr = sqrt(Var) * M_PI * PS_SQR(radKron) / nKronPix;
    282322
    283323    return true;
    284324}
     325
     326enum {Y_P, Y_M, X_P, X_M};
     327
     328bool psphotKronMask (pmSource *source, psImage *kronMask, psImageMaskType markVal) {
     329
     330    int X, Y, e, dXs, dYs;
     331
     332    if (!source->peak) return false;
     333    if (!source->peak->saddlePoints) return false;
     334    if (!source->peak->saddlePoints->n) return false;
     335
     336    // return true;
     337
     338    int xMax = source->maskObj->numCols;
     339    int yMax = source->maskObj->numRows;
     340    // int xOff = source->maskObj->col0;
     341    // int yOff = source->maskObj->row0;
     342
     343    int Xo = ROUND(source->peak->xf) - source->maskObj->col0;
     344    int Yo = ROUND(source->peak->yf) - source->maskObj->row0;
     345
     346    for (int i = 0; i < source->peak->saddlePoints->n; i++) {
     347
     348        psVector *saddlePoint = source->peak->saddlePoints->data[i];
     349        int Xs = ROUND(saddlePoint->data.S32[0]) - source->maskObj->col0;
     350        int Ys = ROUND(saddlePoint->data.S32[1]) - source->maskObj->row0;
     351       
     352        // fprintf (stderr, "mask %d,%d @ %d,%d (%d,%d)\n", Xo, Yo, Xs, Ys, saddlePoint->data.S32[0], saddlePoint->data.S32[1]);
     353
     354        // We want to mask the pixels between the edge of the image and the line perpendicular
     355        // to the connecting line.  Call dX,dY = (Xs-Xo,Ys-Yo).  There are 4 cases:
     356        // +y : dY > 0, fabs(dY) > fabs(dX)
     357        // -y : dY < 0, fabs(dY) > fabs(dX)
     358        // +x : dX > 0, fabs(dX) > fabs(dY)
     359        // -x : dX < 0, fabs(dX) > fabs(dY)
     360
     361        int dX = Xs - Xo;
     362        int dY = Ys - Yo;
     363
     364        // +y : dY > 0, fabs(dY) > fabs(dX)
     365        // -y : dY < 0, fabs(dY) > fabs(dX)
     366        if (fabs(dY) > fabs(dX)) {
     367            // points to right of saddle
     368            Y = Ys;
     369            e = 0;
     370            dXs = (dY > 0) ? +dY : -dY;  // this forces dXs > 0
     371            dYs = (dY > 0) ? -dX : +dX;
     372            for (X = Xs; X < xMax; X++) {
     373                if (dY > 0) {
     374                    for (int Ym = Y; Ym < yMax; Ym++) {
     375                        source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[Ym][X] |= markVal;
     376                        // kronMask->data.S32[Ym+yOff][X+xOff] = source->id;
     377                    }
     378                } else {
     379                    for (int Ym = Y; Ym >= 0; Ym--) {
     380                        source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[Ym][X] |= markVal;
     381                        // kronMask->data.S32[Ym+yOff][X+xOff] = source->id;
     382                    }
     383                }
     384                e += dYs;
     385                int e2 = 2 * e;
     386                if (e2 > dXs) {
     387                    Y ++;
     388                    e -= dXs;
     389                }
     390                if (e2 < -dXs) {
     391                    Y --;
     392                    e += dXs;
     393                }
     394                if (Y >= yMax) break;
     395                if (Y < 0) break;
     396            }
     397            // points to left of saddle
     398            Y = Ys;
     399            e = 0;
     400            for (X = Xs; X >= 0; X--) {
     401                if (dY > 0) {
     402                    for (int Ym = Y; Ym < yMax; Ym++) {
     403                        source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[Ym][X] |= markVal;
     404                        // kronMask->data.S32[Ym+yOff][X+xOff] = source->id;
     405                    }
     406                } else {
     407                    for (int Ym = Y; Ym >= 0; Ym--) {
     408                        source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[Ym][X] |= markVal;
     409                        // kronMask->data.S32[Ym+yOff][X+xOff] = source->id;
     410                    }
     411                }
     412                e -= dYs;
     413                int e2 = 2 * e;
     414                if (e2 > dXs) {
     415                    Y ++;
     416                    e -= dXs;
     417                }
     418                if (e2 < -dXs) {
     419                    Y --;
     420                    e += dXs;
     421                }
     422                if (Y >= yMax) break;
     423                if (Y < 0) break;
     424            }
     425        } else {
     426            // points to right of saddle
     427            X = Xs;
     428            e = 0;
     429            dXs = (dX > 0) ? -dY : +dY;
     430            dYs = (dX > 0) ? +dX : -dX;  // this forces dYs > 0
     431            for (Y = Ys; Y < yMax; Y++) {
     432                if (dX > 0) {
     433                    for (int Xm = X; Xm < xMax; Xm++) {
     434                        source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[Y][Xm] |= markVal;
     435                        // kronMask->data.S32[Y+yOff][Xm+xOff] = source->id;
     436                    }
     437                } else {
     438                    for (int Xm = X; Xm >= 0; Xm--) {
     439                        source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[Y][Xm] |= markVal;
     440                        // kronMask->data.S32[Y+yOff][Xm+xOff] = source->id;
     441                    }
     442                }
     443                e += dXs;
     444                int e2 = 2 * e;
     445                if (e2 > dYs) {
     446                    X ++;
     447                    e -= dYs;
     448                }
     449                if (e2 < -dYs) {
     450                    X --;
     451                    e += dYs;
     452                }
     453                if (X >= yMax) break;
     454                if (X < 0) break;
     455            }
     456            // points to left of saddle
     457            X = Xs;
     458            e = 0;
     459            for (Y = Ys; Y >= 0; Y--) {
     460                if (dX > 0) {
     461                    for (int Xm = X; Xm < xMax; Xm++) {
     462                        source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[Y][Xm] |= markVal;
     463                        // kronMask->data.S32[Y+yOff][Xm+xOff] = source->id;
     464                    }
     465                } else {
     466                    for (int Xm = X; Xm >= 0; Xm--) {
     467                        source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[Y][Xm] |= markVal;
     468                        // kronMask->data.S32[Y+yOff][Xm+xOff] = source->id;
     469                    }
     470                }
     471                e -= dXs;
     472                int e2 = 2 * e;
     473                if (e2 > dYs) {
     474                    X ++;
     475                    e -= dYs;
     476                }
     477                if (e2 < -dYs) {
     478                    X --;
     479                    e += dYs;
     480                }
     481                if (X >= yMax) break;
     482                if (X < 0) break;
     483            }
     484        }
     485    }
     486    return true;
     487}

  • trunk/psphot/src/psphotLoadSRCTEXT.c

    r31452 r32348  
    7373            dPAR[PM_PAR_I0]   = 0.0;
    7474
    75             pmPSF_AxesToModel (PAR, axes);
     75            pmPSF_AxesToModel (PAR, axes, modelType);
    7676
    7777            float peakFlux    = 1.0;

  • trunk/psphot/src/psphotMagnitudes.c

    r31673 r32348  
    44{
    55    bool status = true;
     6
     7    fprintf (stdout, "\n");
     8    psLogMsg ("psphot", PS_LOG_INFO, "--- psphot Magnitudes ---");
    69
    710    // select the appropriate recipe information
    … …  
    146149    psFree (cellGroups);
    147150
    148     psLogMsg ("psphot.magnitudes", PS_LOG_DETAIL, "measure magnitudes : %f sec for %ld objects (%d with apertures)\n", psTimerMark ("psphot.mags"), sources->n, Nap);
     151    psLogMsg ("psphot.magnitudes", PS_LOG_WARN, "measure magnitudes : %f sec for %ld objects (%d with apertures)\n", psTimerMark ("psphot.mags"), sources->n, Nap);
    149152    return true;
    150153}

  • trunk/psphot/src/psphotModelBackground.c

    r31154 r32348  
    145145    psMetadataAddPtr(analysis, PS_LIST_TAIL, "PSPHOT.BACKGROUND.BINNING", PS_DATA_UNKNOWN | PS_META_REPLACE, "Background binning", binning);
    146146
     147    psVector *dQ = psVectorAllocEmpty (100, PS_TYPE_F32);
    147148    psF32 **modelData = model->data.F32;
    148149    psF32 **modelStdevData = modelStdev->data.F32;
    … …  
    216217                }
    217218                modelStdevData[iy][ix] = psStatsGetValue(stats, statsOptionWidth);
     219                // fprintf (stderr, "dQ : %f - %f - %f = %f\n", stats->robustLQ, stats->robustMedian, stats->robustUQ, stats->robustUQ + stats->robustLQ - 2*stats->robustMedian);
     220                psVectorAppend (dQ, stats->robustUQ + stats->robustLQ - 2*stats->robustMedian);
    218221
    219222                // supply sample to plotting routing
    … …  
    322325    }
    323326
    324     psLogMsg ("psphot", PS_LOG_INFO, "built background image: %f sec\n", psTimerMark ("psphot.background"));
    325 
    326     psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "SKY_MEAN", PS_META_REPLACE, "sky mean", Value);
     327    psLogMsg ("psphot", PS_LOG_WARN, "built background image: %f sec\n", psTimerMark ("psphot.background"));
     328
     329    psStats *statsDQ = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN);
     330    psVectorStats (statsDQ, dQ, NULL, NULL, 0);
     331
     332    psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "SKY_MEAN",  PS_META_REPLACE, "sky mean", Value);
    327333    psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "SKY_STDEV", PS_META_REPLACE, "sky stdev", ValueStdev);
    328     psLogMsg ("psphot", PS_LOG_INFO, "image sky : mean %f stdev %f", Value, ValueStdev);
     334    psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "SKY_DQ",    PS_META_REPLACE, "sky quartile slope", statsDQ->sampleMedian);
     335    psLogMsg ("psphot", PS_LOG_INFO, "image sky : mean %f, stdev %f, dQ %f", Value, ValueStdev, statsDQ->sampleMedian);
    329336
    330337    // measure image and background stats and save for later output
    … …  
    334341                                      PS_STAT_MAX);
    335342    psImageStats (statsBck, model, NULL, 0);
    336     psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "MSKY_MN", PS_META_REPLACE, "sky model mean",          statsBck->sampleMean);
     343    psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "MSKY_MN",  PS_META_REPLACE, "sky model mean",          statsBck->sampleMean);
    337344    psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "MSKY_SIG", PS_META_REPLACE, "sky model stdev",         statsBck->sampleStdev);
     345    psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "MSKY_DEV", PS_META_REPLACE, "sky stdev",               ValueStdev);
     346    psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "MSKY_DQ",  PS_META_REPLACE, "sky quartile slope",      statsDQ->sampleMedian);
    338347    psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "MSKY_MAX", PS_META_REPLACE, "sky model maximum value", statsBck->max);
    339348    psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "MSKY_MIN", PS_META_REPLACE, "sky model minimum value", statsBck->min);
    … …  
    343352              statsBck->min, statsBck->sampleMean, statsBck->max, statsBck->sampleStdev);
    344353
     354    psFree(statsDQ);
     355    psFree(dQ);
     356
    345357    psFree(stats);
    346358    psFree(statsBck);
    … …  
    401413    int num = psphotFileruleCount(config, filerule);
    402414
     415    fprintf (stdout, "\n");
     416    psLogMsg ("psphot", PS_LOG_INFO, "--- psphot Model Background ---");
     417
    403418    // loop over the available readouts
    404419    for (int i = 0; i < num; i++) {

  • trunk/psphot/src/psphotModelTest.c

    r31452 r32348  
    1 # include "psphotInternal.h"
    2 # define PM_SOURCE_FIT_PSF_X_EXT PM_SOURCE_FIT_PSF_AND_SKY
     1# include "psphotStandAlone.h"
    32
    4 // XXX add more test information?
    5 bool psphotModelTest (pmConfig *config, const pmFPAview *view, const char *filerule, psMetadata *recipe) {
     3int main (int argc, char **argv) {
    64
    7     bool status;
    8     int modelType = -1;
    9     float obsMag, fitMag, value;
    10     char name[64];
    11     pmPSF *psf = NULL;
    12     pmSourceFitMode fitMode;
     5    psMemInit();
     6    psTimerStart ("complete");
     7    pmErrorRegister();                  // register psModule's error codes/messages
     8    psphotInit();
    139
    14     // bit-masks to test for good/bad pixels
    15     psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT");
    16     assert (maskVal);
     10    // load command-line arguments, options, and system config data
     11    pmConfig *config = psphotModelTestArguments (argc, argv);
     12    assert(config);
    1713
    18     // bit-mask to mark pixels not used in analysis
    19     psImageMaskType markVal = psMetadataLookupImageMask(&status, recipe, "MARK.PSPHOT");
    20     assert (markVal);
     14    psphotVersionPrint();
    2115
    22     // maskVal is used to test for rejected pixels, and must include markVal
    23     maskVal |= markVal;
    24 
    25     // run model fitting tests on a single source?
    26     if (!psMetadataLookupBool (&status, recipe, "TEST_FIT")) return false;
    27 
    28     psTimerStart ("modelTest");
    29 
    30     // find the currently selected readout
    31     pmReadout  *readout = pmFPAfileThisReadout (config->files, view, "PSPHOT.INPUT");
    32     PS_ASSERT_PTR_NON_NULL (readout, false);
    33 
    34     // use poissonian errors or local-sky errors
    35     bool POISSON_ERRORS = psMetadataLookupBool (&status, recipe, filerule);
    36     if (!status) POISSON_ERRORS = true;
    37     pmSourceFitModelInit (15, 0.1, 1.0, POISSON_ERRORS);
    38 
    39     // find the various fitting parameters (try test values first)
    40     float INNER = psMetadataLookupF32 (&status, recipe, "TEST_FIT_INNER_RADIUS");
    41     if (!status || !isfinite(INNER)) {
    42         INNER = psMetadataLookupF32 (&status, recipe, "SKY_INNER_RADIUS");
    43     }
    44     float OUTER = psMetadataLookupF32 (&status, recipe, "TEST_FIT_OUTER_RADIUS");
    45     if (!status || !isfinite(OUTER)) {
    46         OUTER = psMetadataLookupF32 (&status, recipe, "SKY_OUTER_RADIUS");
    47     }
    48     float RADIUS = psMetadataLookupF32 (&status, recipe, "TEST_FIT_RADIUS");
    49     if (!status || !isfinite(RADIUS)) {
    50         RADIUS = psMetadataLookupF32 (&status, recipe, "PSF_FIT_RADIUS");
    51     }
    52     float mRADIUS = psMetadataLookupF32 (&status, recipe, "TEST_MOMENTS_RADIUS");
    53     if (!status || !isfinite(mRADIUS)) {
    54         mRADIUS = psMetadataLookupF32 (&status, recipe, "PSF_MOMENTS_RADIUS");
     16    // load input data (config and images (signal, noise, mask)
     17    if (!psphotParseCamera (config)) {
     18        psErrorStackPrint(stderr, "Error setting up the camera\n");
     19        exit (psphotGetExitStatus());
    5520    }
    5621
    57     // define the source of interest
    58     float xObj     = psMetadataLookupF32 (&status, recipe, "TEST_FIT_X");
    59     float yObj     = psMetadataLookupF32 (&status, recipe, "TEST_FIT_Y");
    60     if (!isfinite(xObj) || !isfinite(yObj)) psAbort ("object position is not defined");
    61 
    62     // what fitting mode to use?
    63     fitMode = PM_SOURCE_FIT_EXT;
    64     char *fitModeWord = psMetadataLookupStr (&status, recipe, "TEST_FIT_MODE");
    65     if (fitModeWord && !strcasecmp (fitModeWord, "PSF")) fitMode = PM_SOURCE_FIT_PSF;
    66     if (fitModeWord && !strcasecmp (fitModeWord, "CONV")) fitMode = PM_SOURCE_FIT_PSF_X_EXT;
    67     if (fitModeWord && !strcasecmp (fitModeWord, "DEFAULT")) fitMode = PM_SOURCE_FIT_EXT;
    68 
    69     // construct the source structures
    70     pmSource *source = pmSourceAlloc();
    71     source->peak = pmPeakAlloc (xObj, yObj, 0, 0);
    72     pmSourceDefinePixels (source, readout, xObj, yObj, OUTER);
    73 
    74     // in fitMode, psf sets the model type
    75     if (fitMode == PM_SOURCE_FIT_PSF) {
    76         psf = psphotLoadPSF (config, view, recipe);
    77         if (!psf) psAbort("PSF_INPUT_FILE not supplied");
    78         modelType = psf->type;
    79         source->type = PM_SOURCE_TYPE_STAR;
    80     }
    81     if (fitMode == PM_SOURCE_FIT_EXT) {
    82         // find the model: supplied by user or first in the PSF_MODEL list
    83         char *modelName  = psMetadataLookupStr (&status, recipe, "TEST_FIT_MODEL");
    84         if (!status || !strcasecmp (modelName, "DEFAULT")) {
    85             // get the list pointers for the PSF_MODEL entries
    86 
    87             psList *list = NULL;
    88             psMetadataItem *mdi = psMetadataLookup (recipe, "PSF_MODEL");
    89             if (mdi == NULL) psAbort("missing PSF_MODEL selection");
    90             if (mdi->type == PS_DATA_STRING) {
    91                 list = psListAlloc(NULL);
    92                 psListAdd (list, PS_LIST_HEAD, mdi);
    93             } else {
    94                 if (mdi->type != PS_DATA_METADATA_MULTI) psAbort("missing PSF_MODEL selection");
    95                 list = psMemIncrRefCounter(mdi->data.list);
    96             }
    97 
    98             // take the first list element
    99             psMetadataItem *item = psListGet (list, PS_LIST_HEAD);
    100             modelName = item->data.V;
    101         }
    102         modelType = pmModelClassGetType (modelName);
    103         if (modelType < 0) psAbort("unknown model %s", modelName);
    104         source->type = PM_SOURCE_TYPE_EXTENDED;
    105     }
    106     if (fitMode == PM_SOURCE_FIT_PSF_X_EXT) {
    107         // we need to load BOTH a psf and an ext model
    108         psf = psphotLoadPSF (config, view, recipe);
    109         if (!psf) psAbort("PSF_INPUT_FILE not supplied");
    110 
    111         // find the model: supplied by user or first in the PSF_MODEL list
    112         char *modelName  = psMetadataLookupStr (&status, recipe, "TEST_FIT_MODEL");
    113         if (!status || !strcasecmp (modelName, "DEFAULT")) {
    114             // get the list pointers for the PSF_MODEL entries
    115 
    116             psList *list = NULL;
    117             psMetadataItem *mdi = psMetadataLookup (recipe, "PSF_MODEL");
    118             if (mdi == NULL) psAbort("missing PSF_MODEL selection");
    119             if (mdi->type == PS_DATA_STRING) {
    120                 list = psListAlloc(NULL);
    121                 psListAdd (list, PS_LIST_HEAD, mdi);
    122             } else {
    123                 if (mdi->type != PS_DATA_METADATA_MULTI) psAbort("missing PSF_MODEL selection");
    124                 list = psMemIncrRefCounter(mdi->data.list);
    125             }
    126 
    127             // take the first list element
    128             psMetadataItem *item = psListGet (list, PS_LIST_HEAD);
    129             modelName = item->data.V;
    130         }
    131         modelType = pmModelClassGetType (modelName);
    132         if (modelType < 0) psAbort("unknown model %s", modelName);
    133         source->type = PM_SOURCE_TYPE_EXTENDED;
     22    // call psphot for each readout
     23    if (!psphotImageLoop (config, PSPHOT_MODEL_TEST)) {
     24        psErrorStackPrint(stderr, "Error in the psphot image loop\n");
     25        exit (psphotGetExitStatus());
    13426    }
    13527
    136     // find the local sky
    137     status = pmSourceLocalSky (source, PS_STAT_SAMPLE_MEDIAN, INNER, maskVal, markVal);
    138     if (!status) psAbort("pmSourceLocalSky error");
     28    psLogMsg ("psphot", 3, "complete psphot run: %f sec\n", psTimerMark ("complete"));
    13929
    140     // get the source moments
    141     status = pmSourceMoments (source, mRADIUS, 0.0, 1.0, maskVal);
    142     if (!status) psAbort("psSourceMoments error");
    143     source->peak->value = source->moments->Peak;
     30    psErrorCode exit_status = psphotGetExitStatus();
     31    psphotCleanup (config);
     32    exit (exit_status);
     33}
    14434
    145     fprintf (stderr, "sum: %f @ (%f, %f)\n", source->moments->Sum, source->moments->Mx, source->moments->My);
    146     fprintf (stderr, "moments: %f, %f - %f\n", source->moments->Mxx, source->moments->Myy, source->moments->Mxy);
    147 
    148     psEllipseMoments moments;
    149     moments.x2 = source->moments->Mxx;
    150     moments.y2 = source->moments->Myy;
    151     moments.xy = source->moments->Mxy;
    152     psEllipseAxes axes = psEllipseMomentsToAxes (moments, 20.0);
    153 
    154     fprintf (stderr, "axes: %f @ (%f, %f)\n", axes.theta*180/M_PI, axes.major, axes.minor);
    155 
    156     // get the initial model parameter guess
    157     pmModel *model = pmSourceModelGuess (source, modelType);
    158     source->modelEXT = model;
    159 
    160     // if any parameters are defined by the user, take those values
    161     int nParams = pmModelClassParameterCount (modelType);
    162     psF32 *params = model->params->data.F32;
    163     params[PM_PAR_XPOS] = xObj; // XXX use the user-supplied value,
    164     params[PM_PAR_YPOS] = yObj; // XXX or use the centroid
    165     for (int i = 0; i < nParams; i++) {
    166         if (i == PM_PAR_XPOS) continue;
    167         if (i == PM_PAR_YPOS) continue;
    168 
    169         sprintf (name, "TEST_FIT_PAR%d", i);
    170         value = psMetadataLookupF32 (&status, recipe, name);
    171         if (status && isfinite (value)) {
    172             params[i] = value;
    173         }
    174     }
    175 
    176     float area = params[4]*params[5];
    177     fprintf (stderr, "peak: %f @ (%f, %f)\n", source->moments->Sum*area, (double)source->peak->x, (double)source->peak->y);
    178 
    179     // for PSF fitting, set the shape parameters based on the PSF & source position
    180     if (fitMode == PM_SOURCE_FIT_PSF) {
    181         source->modelPSF = pmModelFromPSF (model, psf);
    182         psFree (model);
    183         model = source->modelPSF;
    184         params = model->params->data.F32;
    185     }
    186 
    187     // list model input shape
    188     psEllipseShape shape;
    189     shape.sx  = 1.4 / model->params->data.F32[4];
    190     shape.sy  = 1.4 / model->params->data.F32[5];
    191     shape.sxy = model->params->data.F32[6];
    192     axes = psEllipseShapeToAxes (shape, 20.0);
    193 
    194     fprintf (stderr, "guess: %f @ (%f, %f)\n", axes.theta*180/M_PI, axes.major, axes.minor);
    195 
    196     fprintf (stderr, "input parameters: \n");
    197     for (int i = 0; i < nParams; i++) {
    198         fprintf (stderr, "%d : %f\n", i, params[i]);
    199     }
    200 
    201     // define the pixels used for the fit
    202     psImageKeepCircle (source->maskObj, xObj, yObj, RADIUS, "OR", markVal);
    203     psphotSaveImage (NULL, source->maskObj, "mask1.fits");
    204 
    205     char *fitset = psMetadataLookupStr (&status, recipe, "TEST_FIT_SET");
    206     if (status) {
    207         status = psphotFitSet (source, model, fitset, fitMode, maskVal);
    208         exit (0);
    209     }
    210 
    211     if (fitMode == PM_SOURCE_FIT_PSF_X_EXT) {
    212         // build the psf for the object
    213         source->modelPSF = pmModelFromPSF (model, psf);
    214         source->modelEXT = model;
    215 
    216         // what fraction of the PSF is used? (radius in pixels : 2 -> 5x5 box)
    217         int psfSize = psMetadataLookupS32 (&status, recipe, "PCM_BOX_SIZE");
    218         assert (status);
    219 
    220         model = psphotPSFConvModel (readout, source, modelType, maskVal, markVal, psfSize);
    221         params = model->params->data.F32;
    222     } else {
    223         status = pmSourceFitModel (source, model, fitMode, maskVal);
    224     }
    225 
    226     // measure the source mags
    227     pmSourcePhotometryModel (&fitMag, model);
    228     pmSourcePhotometryAper  (NULL, &obsMag, NULL, NULL, model, source->pixels, source->variance, source->maskObj, maskVal);
    229     fprintf (stderr, "ap: %f, fit: %f, apmifit: %f, nIter: %d\n", obsMag, fitMag, obsMag - fitMag, model->nIter);
    230 
    231     // write out positive object
    232     psphotSaveImage (NULL, source->pixels, "object.fits");
    233 
    234     // subtract object, leave local sky
    235     // pmModelSub (source->pixels, source->maskObj, model, PM_MODEL_OP_FULL, maskVal);
    236     pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
    237 
    238     fprintf (stderr, "output parameters: \n");
    239     for (int i = 0; i < nParams; i++) {
    240         fprintf (stderr, "%d : %f\n", i, params[i]);
    241     }
    242 
    243     // write out
    244     psphotSaveImage (NULL, source->pixels, "resid.fits");
    245     psphotSaveImage (NULL, source->maskObj, "mask.fits");
    246 
    247     psLogMsg ("psphot", PS_LOG_INFO, "model test : %f sec\n", psTimerMark ("modelTest"));
    248 
    249     exit (0);
    250 }
     35// all functions which return to this level must raise one of the top-level error codes if they
     36// exit with an error.  these error codes are used to specify the program exit status

  • trunk/psphot/src/psphotOutput.c

    r31452 r32348  
    325325    psMetadataItemSupplement (&status, header, analysis, "MSKY_NY");
    326326
     327    psMetadataItemSupplement (&status, header, analysis, "MSKY_DEV");
     328    psMetadataItemSupplement (&status, header, analysis, "MSKY_DQ");
     329
    327330    psMetadataItemSupplement (&status, header, analysis, "DETEFF.MAGREF");
    328331

  • trunk/psphot/src/psphotPetrosianProfile.c

    r30624 r32348  
    3838    // convert the isophotal radius vs angle measurements to an elliptical contour
    3939    if (!psphotEllipticalContour (source, petrosian)) {
    40         psLogMsg ("psphot", 3, "failed to measure elliptical contour");
     40        // psLogMsg ("psphot", 3, "failed to measure elliptical contour");
    4141        psFree (petrosian);
    4242        return false;

  • trunk/psphot/src/psphotPetrosianStats.c

    r31452 r32348  
    143143    // if we failed to reach the PETROSIAN_RATIO, use the lowest significant ratio instead (flag this!)
    144144    if (!anyPetro) {
    145         // interpolate Rvec between i-1 and i to PETROSIAN_RATIO to get flux (Fvec) and radius (rvec)
    146         if (lowestSignificantRadius == 0) {
    147             // assume Fmax @ R = 0.0
    148             petRadius    = InterpolateValues     (1.0, 0.0, petRatio->data.F32[lowestSignificantRadius], refRadius->data.F32[lowestSignificantRadius], PETROSIAN_RATIO);
    149             petRadiusErr = InterpolateValuesErrX (1.0, 0.0, petRatio->data.F32[lowestSignificantRadius], refRadius->data.F32[lowestSignificantRadius], PETROSIAN_RATIO, 0.0, petRatioErr->data.F32[lowestSignificantRadius]);
    150 
    151         } else {
    152             int n0 = lowestSignificantRadius-1;
    153             int n1 = lowestSignificantRadius;
    154             petRadius    = InterpolateValues     (petRatio->data.F32[n0], refRadius->data.F32[n0], petRatio->data.F32[n1], refRadius->data.F32[n1], PETROSIAN_RATIO);
    155             petRadiusErr = InterpolateValuesErrX (petRatio->data.F32[n0], refRadius->data.F32[n0], petRatio->data.F32[n1], refRadius->data.F32[n1], PETROSIAN_RATIO, petRatioErr->data.F32[n0], petRatioErr->data.F32[n1]);
     145        petRadius = refRadius->data.F32[lowestSignificantRadius];
     146        petRadiusErr = NAN;
     147        if (!isfinite(petRadius)) {
     148            fprintf (stderr, "nan pet radius\n");
    156149        }
    157150    }
    … …  
    170163                petArea    = InterpolateValues     (refRadius->data.F32[i-1], areaSum->data.F32[i-1], refRadius->data.F32[i], areaSum->data.F32[i], apRadius);
    171164                petApix    = InterpolateValues     (refRadius->data.F32[i-1], apixSum->data.F32[i-1], refRadius->data.F32[i], apixSum->data.F32[i], apRadius);
     165                if (!isfinite(petFlux)) {
     166                    fprintf (stderr, "nan pet flux\n");
     167                }
    172168                break;
    173169            }
    … …  
    188184        if (!found50 && (fluxSum->data.F32[i] > flux50)) {
    189185            if (i == 0) {
    190                 psWarning ("does this case make any sense? (fluxSum[0] %f > flux50 %f)", fluxSum->data.F32[i], flux50);
    191                 continue;
     186                R50    = InterpolateValues     (fluxSum->data.F32[i], refRadius->data.F32[i], fluxSum->data.F32[i+1], refRadius->data.F32[i+1], flux50);
     187                R50err = InterpolateValuesErrX (fluxSum->data.F32[i], refRadius->data.F32[i], fluxSum->data.F32[i+1], refRadius->data.F32[i+1], flux50, sqrt(fluxSumErr2->data.F32[i]), sqrt(fluxSumErr2->data.F32[i+1]));
     188                found50 = true;
    192189            } else {
    193190                R50    = InterpolateValues     (fluxSum->data.F32[i-1], refRadius->data.F32[i-1], fluxSum->data.F32[i], refRadius->data.F32[i], flux50);
    … …  
    198195        if (!found90 && (fluxSum->data.F32[i] > flux90)) {
    199196            if (i == 0) {
    200                 psWarning ("does this case make any sense? (fluxSum[0] %f > flux90 %f)", fluxSum->data.F32[i], flux90);
    201                 continue;
     197                R90    = InterpolateValues     (fluxSum->data.F32[i], refRadius->data.F32[i], fluxSum->data.F32[i+1], refRadius->data.F32[i+1], flux90);
     198                R90err = InterpolateValuesErrX (fluxSum->data.F32[i], refRadius->data.F32[i], fluxSum->data.F32[i+1], refRadius->data.F32[i+1], flux90, sqrt(fluxSumErr2->data.F32[i]), sqrt(fluxSumErr2->data.F32[i+1]));
     199                found90 = true;
    202200            } else {
    203201                R90    = InterpolateValues     (fluxSum->data.F32[i-1], refRadius->data.F32[i-1], fluxSum->data.F32[i], refRadius->data.F32[i], flux90);

  • trunk/psphot/src/psphotRadialApertures.c

    r31452 r32348  
    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;
    … …  
    6269        return true;
    6370    }
     71
     72    // determine the number of allowed threads
     73    int nThreads = psMetadataLookupS32(&status, config->arguments, "NTHREADS"); // Number of threads
     74    if (!status) {
     75        nThreads = 0;
     76    }
     77
     78    // source analysis is done in S/N order (brightest first)
     79    // XXX are we getting the objects out of order? does it matter?
     80    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    }
     94
     95    // option to limit analysis to a specific region
     96    char *region = psMetadataLookupStr (&status, recipe, "ANALYSIS_REGION");
     97    psRegion *AnalysisRegion = psRegionAlloc(0,0,0,0);
     98    *AnalysisRegion = psRegionForImage (readout->image, psRegionFromString (region));
     99    if (psRegionIsNaN (*AnalysisRegion)) psAbort("analysis region mis-defined");
     100
     101    // choose Cx, Cy (see psphotThreadTools.c for overview of the concepts)
     102    int Cx = 1, Cy = 1;
     103    psphotChooseCellSizes (&Cx, &Cy, readout, nThreads);
     104
     105    psArray *cellGroups = psphotAssignSources (Cx, Cy, sources);
     106
     107    for (int i = 0; i < cellGroups->n; i++) {
     108
     109        psArray *cells = cellGroups->data[i];
     110
     111        for (int j = 0; j < cells->n; j++) {
     112
     113            // allocate a job -- if threads are not defined, this just runs the job
     114            psThreadJob *job = psThreadJobAlloc ("PSPHOT_RADIAL_APERTURES");
     115
     116            psArrayAdd(job->args, 1, readout);
     117            psArrayAdd(job->args, 1, cells->data[j]); // sources
     118            psArrayAdd(job->args, 1, AnalysisRegion);
     119            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);
     123
     124            PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nradial
     125
     126// set this to 0 to run without threading
     127# if (1)           
     128            if (!psThreadJobAddPending(job)) {
     129                psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
     130                psFree(AnalysisRegion);
     131                return false;
     132            }
     133# else
     134            if (!psphotRadialApertures_Threaded(job)) {
     135                psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
     136                psFree(AnalysisRegion);
     137                return false;
     138            }
     139            psScalar *scalar = NULL;
     140            scalar = job->args->data[6];
     141            Nradial += scalar->data.S32;
     142            psFree(job);
     143# endif
     144        }
     145
     146        // wait for the threads to finish and manage results
     147        if (!psThreadPoolWait (false)) {
     148            psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
     149            psFree(AnalysisRegion);
     150            return false;
     151        }
     152
     153        // we have only supplied one type of job, so we can assume the types here
     154        psThreadJob *job = NULL;
     155        while ((job = psThreadJobGetDone()) != NULL) {
     156            if (job->args->n < 1) {
     157                fprintf (stderr, "error with job\n");
     158            } else {
     159                psScalar *scalar = NULL;
     160                scalar = job->args->data[6];
     161                Nradial += scalar->data.S32;
     162            }
     163            psFree(job);
     164        }
     165    }
     166    psFree (cellGroups);
     167    psFree(AnalysisRegion);
     168
     169    psLogMsg ("psphot", PS_LOG_WARN, "radial source apertures: %f sec for %d objects\n", psTimerMark ("psphot.radial"), Nradial);
     170    return true;
     171}
     172 
     173bool psphotRadialApertures_Threaded (psThreadJob *job) {
     174
     175    bool status;
     176    int Nradial = 0;
     177
     178    // arguments: readout, sources, models, region, psfSize, maskVal, markVal
     179    pmReadout *readout      = job->args->data[0];
     180    psArray *sources        = job->args->data[1];
     181    psRegion *region        = job->args->data[2];
     182    psMetadata *recipe      = job->args->data[3];
     183
     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
    64186
    65187    // radMax stores the upper bounds of the annuli
    … …  
    68190    psAssert (radMax, "annular bins (RADIAL.ANNULAR.BINS.UPPER) are not defined in the recipe");
    69191    psAssert (radMax->n, "no valid annular bins (RADIAL.ANNULAR.BINS.UPPER) are define");
    70     float outerRadius = radMax->data.F32[radMax->n - 1];
    71192
    72193    // user-defined masks to test for good/bad pixels (build from recipe list if not yet set)
    … …  
    77198    float SN_LIM = psMetadataLookupF32 (&status, recipe, "RADIAL_APERTURES_SN_LIM");
    78199
    79     // source analysis is done in S/N order (brightest first)
    80     // XXX are we getting the objects out of order? does it matter?
    81     sources = psArraySort (sources, pmSourceSortByFlux);
    82 
    83     // option to limit analysis to a specific region
    84     char *region = psMetadataLookupStr (&status, recipe, "ANALYSIS_REGION");
    85     psRegion AnalysisRegion = psRegionForImage (readout->image, psRegionFromString (region));
    86     if (psRegionIsNaN (AnalysisRegion)) psAbort("analysis region mis-defined");
     200    float outerRadius = radMax->data.F32[radMax->n - 1];
    87201
    88202    // choose the sources of interest
    … …  
    104218
    105219        // limit selection by analysis region
    106         if (source->peak->x < AnalysisRegion.x0) continue;
    107         if (source->peak->y < AnalysisRegion.y0) continue;
    108         if (source->peak->x > AnalysisRegion.x1) continue;
    109         if (source->peak->y > AnalysisRegion.y1) continue;
     220        if (source->peak->x < region->x0) continue;
     221        if (source->peak->y < region->y0) continue;
     222        if (source->peak->x > region->x1) continue;
     223        if (source->peak->y > region->y1) continue;
    110224
    111225        // allocate pmSourceExtendedParameters, if not already defined
    112         if (!source->radialAper) {
    113             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            }
    114235        }
    115236
    … …  
    131252        pmSourceRedefinePixels (source, readout, source->peak->xf, source->peak->yf, outerRadius + 2);
    132253
    133         if (!psphotRadialApertureSource (source, recipe, maskVal, radMax, 0)) {
     254        if (!psphotRadialApertureSource (source, recipe, maskVal, radMax, entry)) {
    134255            psTrace ("psphot", 5, "failed to extract radial profile for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My);
    135256        } else {
    … …  
    145266        pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
    146267    }
    147 
    148     psLogMsg ("psphot", PS_LOG_INFO, "radial source apertures: %f sec for %d objects\n", psTimerMark ("psphot.radial"), Nradial);
     268    psScalar *scalar = job->args->data[6];
     269    scalar->data.S32 = Nradial;
     270
    149271    return true;
    150272}

  • trunk/psphot/src/psphotRadialAperturesByObject.c

    r31154 r32348  
    162162            // re-subtract the object, leave local sky
    163163            pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
    164 
    165             // psLogMsg("psphot", PS_LOG_INFO, "radial apertures for %d", index);
    166             // psphotVisualShowImage(readout);
    167164        }
    168165    }

  • trunk/psphot/src/psphotRadialProfile.c

    r31154 r32348  
    3939    // convert the isophotal radius vs angle measurements to an elliptical contour
    4040    if (!psphotEllipticalContour (source)) {
    41         psLogMsg ("psphot", 3, "failed to measure elliptical contour");
     41        // psLogMsg ("psphot", 3, "failed to measure elliptical contour");
    4242        return false;
    4343    }

  • trunk/psphot/src/psphotRadiusChecks.c

    r31452 r32348  
    142142}
    143143
     144# define MIN_WINDOW 5.0
     145# define SCALE1 5.0
     146# define SCALE2 12.0
     147
     148// call this function whenever you (re)-define the EXT model
     149// XXX this function does not shrink the window
     150bool psphotSetRadiusMoments (float *fitRadius, float *windowRadius, pmReadout *readout, pmSource *source, psImageMaskType markVal) {
     151
     152    psAssert (source, "source not defined??");
     153    psAssert (source->moments, "moments not defined??");
     154
     155    *fitRadius = SCALE1 * source->moments->Mrf;
     156    *fitRadius = PS_MIN (PS_MAX(*fitRadius, MIN_WINDOW), EXT_FIT_MAX_RADIUS);
     157
     158    *windowRadius = SCALE2 * source->moments->Mrf;
     159    *windowRadius = PS_MIN (PS_MAX(*windowRadius, 2.5*MIN_WINDOW), 2.5*EXT_FIT_MAX_RADIUS);
     160
     161    // redefine the pixels if needed
     162    pmSourceRedefinePixels (source, readout, source->peak->xf, source->peak->yf, *windowRadius);
     163
     164    // set the mask to flag the excluded pixels
     165    psImageKeepCircle (source->maskObj, source->peak->xf, source->peak->yf, *fitRadius, "OR", markVal);
     166
     167    return true;
     168}
     169# undef SCALE1
     170# undef SCALE2
     171# undef MIN_WINDOW
     172
     173# define MIN_WINDOW 5.0
     174# define SCALE1 5.0
     175# define PAD_WINDOW 3.0
     176
     177// call this function whenever you (re)-define the EXT model
     178// XXX alternate function to set exactly the desired window size
     179bool psphotSetRadiusMomentsExact (float *fitRadius, float *windowRadius, pmReadout *readout, pmSource *source, psImageMaskType markVal) {
     180
     181    psRegion newRegion;
     182
     183    psAssert (source, "source not defined??");
     184    psAssert (source->moments, "moments not defined??");
     185
     186    *fitRadius = SCALE1 * source->moments->Mrf;
     187    *fitRadius = PS_MIN (PS_MAX(*fitRadius, MIN_WINDOW), EXT_FIT_MAX_RADIUS);
     188
     189    *windowRadius = *fitRadius + PAD_WINDOW;
     190
     191    // check to see if new region is completely contained within old region
     192    newRegion = psRegionForSquare (source->peak->xf, source->peak->yf, *windowRadius);
     193    newRegion = psRegionForImage (readout->image, newRegion);
     194
     195    // redefine the pixels to match
     196    pmSourceRedefinePixelsByRegion (source, readout, newRegion);
     197
     198    // set the mask to flag the excluded pixels
     199    psImageKeepCircle (source->maskObj, source->peak->xf, source->peak->yf, *fitRadius, "OR", markVal);
     200
     201    return true;
     202}
     203
    144204// call this function whenever you (re)-define the EXT model
    145205bool psphotSetRadiusFootprint (float *radius, pmReadout *readout, pmSource *source, psImageMaskType markVal, float factor) {
    … …  
    181241}
    182242
     243// call this function whenever you (re)-define the EXT model
     244bool psphotMaskFootprint (pmReadout *readout, pmSource *source, psImageMaskType markVal) {
     245
     246    psAssert (source, "source not defined??");
     247    psAssert (source->peak, "peak not defined??");
     248
     249    pmPeak *peak = source->peak;
     250
     251    // set the radius based on the footprint:
     252    if (!peak->footprint) return false;
     253    pmFootprint *footprint = peak->footprint;
     254    if (!footprint->spans) return false;
     255    if (footprint->spans->n < 1) return false;
     256
     257    int Xo = source->maskObj->col0;
     258    int Yo = source->maskObj->row0;
     259
     260    // mark all pixels
     261    for (int j = 0; j < source->maskObj->numRows; j++) {
     262        for (int i = 0; i < source->maskObj->numCols; i++) {
     263            source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[j][i] |= markVal;
     264        }
     265    }
     266
     267    psImageMaskType clearVal = PS_NOT_IMAGE_MASK(markVal);
     268
     269    for (int j = 0; j < footprint->spans->n; j++) {
     270        pmSpan *span = footprint->spans->data[j];
     271
     272        // mask the rows before and after each span
     273        int minX = span->x0 - Xo - 2;
     274        int maxX = span->x1 - Xo + 2;
     275        int myY = span->y - Yo;
     276
     277        // unmark pixels inside the footprint
     278        for (int i = minX; i <= maxX; i++) {
     279            source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[myY][i] &= clearVal;
     280        }
     281    }
     282    return true;
     283}
     284
    183285// alternative EXT radius based on model guess (for use without footprints)
    184286bool psphotSetRadiusModel (pmModel *model, pmReadout *readout, pmSource *source, psImageMaskType markVal, bool deep) {

  • trunk/psphot/src/psphotReadout.c

    r31452 r32348  
    102102        return psphotReadoutCleanup (config, view, filerule);
    103103    }
    104 
    105 # if (0)
    106     // XXX test to mask outliers, not very helpful
    107     // mask the high values in the image (with MARK)
    108     if (!psphotMaskBackground (config, view, filerule)) {
    109         return psphotReadoutCleanup (config, view, filerule);
    110     }
    111 
    112     // generate a background model (median, smoothed image)
    113     if (!psphotModelBackground (config, view, filerule)) {
    114         return psphotReadoutCleanup (config, view, filerule);
    115     }
    116 # endif
    117 
    118104    if (!psphotSubtractBackground (config, view, filerule)) {
    119105        return psphotReadoutCleanup (config, view, filerule);
    … …  
    137123    }
    138124
    139     // construct sources and measure basic stats (saved on detections->newSources)
     125    // construct sources and measure moments and other basic stats (saved on detections->newSources)
     126    // all sources use the auto-scaled window appropriate to a PSF, except for the saturated
     127    // stars : these use a larger window (3x the basic window)
    140128    if (!psphotSourceStats (config, view, filerule, true)) { // pass 1
    141129        psError(PSPHOT_ERR_UNKNOWN, false, "failure to generate sources");
    … …  
    153141
    154142    // mark blended peaks PS_SOURCE_BLEND (detections->newSources)
    155     if (!psphotBasicDeblend (config, view, filerule)) {
     143    // XXX I've deactivated this because it was preventing galaxies close to stars from being
     144    // XXX fitted as an extended source.
     145    if (false && !psphotBasicDeblend (config, view, filerule)) {
    156146        psError (PSPHOT_ERR_UNKNOWN, false, "failed on deblend analysis");
    157147        return psphotReadoutCleanup (config, view, filerule);
    … …  
    200190    psphotDumpChisqs (config, view, filerule);
    201191
    202     // XXX re-measure the kron mags with models subtracted
    203     psphotKronMasked(config, view, filerule);
     192    // re-measure the kron mags with models subtracted.  this pass uses a circular window of size PSF_MOMENTS_RADIUS (same window used to measure the psf-scale moments)
     193   
     194    // but this is chosen above to be appropriate for the PSF objects (not galaxies)
     195    // psphotKronMasked(config, view, filerule);
     196    psphotKronIterate(config, view, filerule);
    204197
    205198    // identify CRs and extended sources (only unmeasured sources are measured)
    … …  
    240233        psphotSubNoise (config, view, filerule); // pass 1 (detections->allSources)
    241234
    242         // define new sources based on only the new peaks
     235        // define new sources based on only the new peaks & measure moments
    243236        // NOTE: new sources are saved on detections->newSources
    244237        psphotSourceStats (config, view, filerule, false); // pass 2 (detections->newSources)
    … …  
    312305pass1finish:
    313306
    314     // XXX re-measure the kron mags with models subtracted
    315     psphotKronMasked(config, view, filerule);
     307    // re-measure the kron mags with models subtracted
     308    // psphotKronMasked(config, view, filerule);
     309    psphotKronIterate(config, view, filerule);
    316310
    317311    // measure source size for the remaining sources
    … …  
    321315    psphotExtendedSourceAnalysis (config, view, filerule); // pass 1 (detections->allSources)
    322316    psphotExtendedSourceFits (config, view, filerule); // pass 1 (detections->allSources)
    323     psphotRadialApertures(config, view, filerule);
     317    psphotRadialApertures(config, view, filerule, 0);
    324318
    325319finish:
    … …  
    359353    }
    360354
     355    psLogMsg ("psphot.readout", PS_LOG_WARN, "complete psphot readout : %f sec\n", psTimerMark ("psphotReadout"));
     356
    361357    // create the exported-metadata and free local data
    362358    return psphotReadoutCleanup(config, view, filerule);

  • trunk/psphot/src/psphotRoughClass.c

    r31673 r32348  
    1111{
    1212    bool status = true;
     13
     14    fprintf (stdout, "\n");
     15    psLogMsg ("psphot", PS_LOG_INFO, "--- psphot Rough Class ---");
    1316
    1417    // select the appropriate recipe information
    … …  
    128131    psphotDumpMoments (recipe, sources);
    129132
    130     psLogMsg ("psphot.roughclass", PS_LOG_INFO, "rough classification: %f sec\n", psTimerMark ("psphot.rough"));
     133    psLogMsg ("psphot.roughclass", PS_LOG_WARN, "rough classification: %f sec\n", psTimerMark ("psphot.rough"));
    131134
    132135    psphotVisualPlotMoments (recipe, readout->analysis, sources);

  • trunk/psphot/src/psphotSetThreads.c

    r31452 r32348  
    3030    psFree(task);
    3131
     32    task = psThreadTaskAlloc("PSPHOT_KRON_ITERATE", 7);
     33    task->function = &psphotKronIterate_Threaded;
     34    psThreadTaskAdd(task);
     35    psFree(task);
     36
    3237    task = psThreadTaskAlloc("PSPHOT_BLEND_FIT", 10);
    3338    task->function = &psphotBlendFit_Threaded;
    … …  
    4045    psFree(task);
    4146
     47    task = psThreadTaskAlloc("PSPHOT_EXTENDED_ANALYSIS", 8);
     48    task->function = &psphotExtendedSourceAnalysis_Threaded;
     49    psThreadTaskAdd(task);
     50    psFree(task);
     51
     52    task = psThreadTaskAlloc("PSPHOT_RADIAL_APERTURES", 7);
     53    task->function = &psphotRadialApertures_Threaded;
     54    psThreadTaskAdd(task);
     55    psFree(task);
     56
    4257    return true;
    4358}

  • trunk/psphot/src/psphotSourceFits.c

    r31452 r32348  
    211211bool psphotFitBlob (pmReadout *readout, pmSource *source, psArray *newSources, pmPSF *psf, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal) {
    212212
    213     float radius;
     213    float fitRadius, windowRadius;
    214214    bool okEXT, okDBL;
    215215    pmModel *ONE = NULL;
    … …  
    217217    pmModel *EXT = NULL;
    218218    psArray *DBL = NULL;
    219     pmMoments psfMoments;
     219    // pmMoments psfMoments;
    220220
    221221    // skip the source if we don't think it is extended
    … …  
    225225    if (source->type == PM_SOURCE_TYPE_SATURATED) return false;
    226226
     227# define TEST_X -420.0
     228# define TEST_Y 300.0
     229   
     230    if ((fabs(source->peak->xf - TEST_X) < 5) && (fabs(source->peak->yf - TEST_Y) < 5)) {
     231        fprintf (stderr, "test galaxy\n");
     232    }
     233
     234# undef TEST_X
     235# undef TEST_Y
     236
    227237    // set the radius based on the footprint (also sets the mask pixels)
    228     if (!psphotSetRadiusFootprint(&radius, readout, source, markVal, 1.0)) return false;
    229 
    230     // XXX note that this changes the source moments that are published...
    231     // XXX all published moments should use the same measurement
    232     // recalculate the source moments using the larger extended-source moments radius
    233     // at this stage, skip Gaussian windowing, and do not clip pixels by S/N
    234     // this uses the footprint to judge both radius and aperture?
    235     // XXX save the psf-based moments for output
    236     psfMoments = *source->moments;
    237     if (!pmSourceMoments (source, radius, 0.0, 0.5, 0.0, maskVal)) {
    238       *source->moments = psfMoments;
    239       return false;
    240     }
     238    if (!psphotSetRadiusMoments(&fitRadius, &windowRadius, readout, source, markVal)) return false;
     239    // fprintf (stderr, "rad: %6.1f %6.1f  | %5.2f %5.2f %5.2f  ", source->peak->xf, source->peak->yf, source->moments->Mrf, fitRadius, windowRadius);
    241240
    242241    psTrace ("psphot", 5, "trying blob...\n");
    … …  
    263262        ONE = DBL->data[0];
    264263        if (ONE) {
    265             if (!isfinite(ONE->params->data.F32[PM_PAR_I0])) psAbort("nan in fit");
     264            psAssert (isfinite(ONE->params->data.F32[PM_PAR_I0]), "nan in fit");
    266265            chiDBL = ONE->chisqNorm; // save chisq for double-star/galaxy comparison
    267             ONE->fitRadius = radius;
     266            ONE->fitRadius = fitRadius;
    268267        }
    269268
    … …  
    271270        ONE = DBL->data[1];
    272271        if (ONE) {
    273             if (!isfinite(ONE->params->data.F32[PM_PAR_I0])) psAbort("nan in fit");
    274             ONE->fitRadius = radius;
     272            psAssert (isfinite(ONE->params->data.F32[PM_PAR_I0]), "nan in fit");
     273            ONE->fitRadius = fitRadius;
    275274        }
    276275    }
    … …  
    284283        okEXT = psphotEvalEXT (tmpSrc, EXT);
    285284        chiEXT = EXT ? EXT->chisqNorm : NAN;
     285        EXT->fitRadius = fitRadius;
    286286    }
    287287
    … …  
    294294
    295295    if (okEXT && okDBL) {
    296         psTrace ("psphot", 5, "blob chisq: %f vs %f\n", chiEXT, chiDBL);
    297296        // XXX EAM : a bogus bias: need to examine this better
    298297        if (3*chiEXT > chiDBL) goto keepDBL;
    … …  
    303302    if (!okEXT && okDBL) goto keepDBL;
    304303
     304    psTrace ("psphot", 4, "both failed: blob chisq: %f vs %f for %f,%f\n", chiEXT, chiDBL, source->peak->xf, source->peak->yf);
     305
    305306    // both models failed; reject them both
    306307    // XXX -- change type flags to psf in this case, and make sure we subtract it?
    307308    // reset the psf moments
    308     *source->moments = psfMoments;
     309    // XXX *source->moments = psfMoments;
    309310
    310311    psFree (EXT);
    … …  
    313314
    314315keepEXT:
     316    psTrace ("psphot", 4, "goto EXT : blob chisq: %f vs %f for %f,%f\n", chiEXT, chiDBL, source->peak->xf, source->peak->yf);
    315317    // sub EXT
    316318    psFree (DBL);
    … …  
    338340
    339341    // reset the psf moments
    340     *source->moments = psfMoments;
     342    // XXX *source->moments = psfMoments;
    341343    return true;
    342344
    343345keepDBL:
     346    psTrace ("psphot", 4, "goto DBL : blob chisq: %f vs %f for %f,%f\n", chiEXT, chiDBL, source->peak->xf, source->peak->yf);
    344347    // sub DLB
    345348    psFree (EXT);
    … …  
    372375            psLogMsg ("psphot", 1, "PAR %d : %f +/- %f\n", i, source->modelPSF->params->data.F32[i], source->modelPSF->dparams->data.F32[i]);
    373376        }
    374         psphotVisualShowResidualImage (readout, false);
    375377    }
    376378# endif
    377 
    378     // reset the (original) psf moments
    379     *source->moments = psfMoments;
    380     *newSrc->moments = psfMoments;
    381379
    382380    psArrayAdd (newSources, 100, newSrc);
    … …  
    386384
    387385escape:
    388     // reset the psf moments
    389     *source->moments = psfMoments;
    390386    psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal));
    391387    psFree (tmpSrc);
    … …  
    476472    // for sersic models, use a grid search to choose an index, then float the params there
    477473    if (modelType == pmModelClassGetType("PS_MODEL_SERSIC")) {
    478         // for the test fits, use a somewhat smaller radius
    479         psphotSetRadiusFootprint(&model->fitRadius, readout, source, markVal, 0.5);
    480474        psphotFitSersicIndex (model, readout, source, fitOptions, maskVal, markVal);
    481     }
    482 
    483     if (!psphotSetRadiusModel (model, readout, source, markVal, true)) {
    484         psphotSetRadiusFootprint(&model->fitRadius, readout, source, markVal, 1.0);
    485475    }
    486476
    … …  
    494484    pmSourceFitModel (source, model, &options, maskVal);
    495485    // fprintf (stderr, "chisq: %f, nIter: %d, radius: %f, npix: %d\n", model->chisqNorm, model->nIter, model->fitRadius, model->nPix);
    496 
    497486    // psTraceSetLevel("psLib.math.psMinimizeLMChi2", 0);
     487
    498488    return (model);
    499489}
     490
     491# define TIMING 0
    500492
    501493pmModel *psphotFitPCM (pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, pmModelType modelType, psImageMaskType maskVal, psImageMaskType markVal, int psfSize) {
    … …  
    517509        return NULL;
    518510    }
     511
     512# define TEST_X -540.0
     513# define TEST_Y 540.0
     514   
     515    if ((fabs(source->peak->xf - TEST_X) < 5) && (fabs(source->peak->yf - TEST_Y) < 5)) {
     516        psTraceSetLevel("psModules.objects.pmPCM_MinimizeChisq", 5);
     517    }
     518
     519    float t1, t2, t4, t5;
     520    if (TIMING) { psTimerStart ("psphotFitPCM"); }
    519521
    520522    pmPCMdata *pcm = pmPCMinit (source, &options, model, maskVal, psfSize);
    … …  
    524526        return model;
    525527    }
    526 
    527     // use the source moments, etc to guess basic model parameters
    528     if (!pmSourceModelGuessPCM (pcm, source, maskVal, markVal)) {
    529         psFree (pcm);
    530         model->flags |= PM_MODEL_STATUS_BADARGS;
    531         return model;
    532     }
    533 
    534     // for sersic models, use a grid search to choose an index, then float the params there
     528    if (TIMING) { t1 = psTimerMark ("psphotFitPCM"); }
     529
     530    // get the guess for sersic models
    535531    if (modelType == pmModelClassGetType("PS_MODEL_SERSIC")) {
    536         // for the test fits, use a somewhat smaller radius
    537         psphotSetRadiusFootprint(&model->fitRadius, readout, source, markVal, 0.5);
    538 
    539         if (!psphotFitSersicIndexPCM (pcm, readout, source, fitOptions, maskVal, markVal, psfSize)) {
     532        // use the source moments, etc to guess basic model parameters
     533        if (!psphotSersicModelClassGuessPCM (pcm, source)) {
    540534            psFree (pcm);
    541535            model->flags |= PM_MODEL_STATUS_BADARGS;
    542536            return model;
    543537        }
    544     }
    545 
    546     if (!psphotSetRadiusModel (model, readout, source, markVal, true)) {
    547         psphotSetRadiusFootprint(&model->fitRadius, readout, source, markVal, 1.0);
    548     }
     538    } else {
     539        // use the source moments, etc to guess basic model parameters
     540        if (!pmSourceModelGuessPCM (pcm, source, maskVal, markVal)) {
     541            psFree (pcm);
     542            model->flags |= PM_MODEL_STATUS_BADARGS;
     543            return model;
     544        }
     545    }
     546
     547    if (TIMING) { t2 = psTimerMark ("psphotFitPCM"); }
    549548
    550549    if (modelType == pmModelClassGetType("PS_MODEL_SERSIC")) {
    … …  
    555554    // update the pcm elements if we have changed the circumstance (options.mode or source->pixels)
    556555    pmPCMupdate(pcm, source, &options, model);
     556    if (TIMING) { t4 = psTimerMark ("psphotFitPCM"); }
    557557
    558558    // psTraceSetLevel("psLib.math.psMinimizeLMChi2", 5);
    559559    pmSourceFitPCM (pcm, source, &options, maskVal, markVal, psfSize);
    560     // fprintf (stderr, "chisq: %f, nIter: %d, radius: %f, npix: %d\n", model->chisqNorm, model->nIter, model->fitRadius, model->nPix);
     560    if (TIMING) { t5 = psTimerMark ("psphotFitPCM"); }
     561
     562    if (TIMING) {
     563        int nPixBig = source->pixels->numCols * source->pixels->numRows;
     564        fprintf (stderr, "psphotFitPCM : nIter: %2d, radius: %6.1f, npix: %5d of %5d, t1: %6.4f, t2: %6.4f, t4: %6.4f, t5: %6.4f\n", model->nIter, model->fitRadius, model->nPix, nPixBig, t1, t2, t4, t5);
     565    }
     566
     567    if ((fabs(source->peak->xf - TEST_X) < 5) && (fabs(source->peak->yf - TEST_Y) < 5)) {
     568        psTraceSetLevel("psModules.objects.pmPCM_MinimizeChisq", 0);
     569    }
    561570
    562571    // psTraceSetLevel("psLib.math.psMinimizeLMChi2", 0);
    … …  
    566575}
    567576
     577# undef TEST_X
     578# undef TEST_Y
     579
    568580// note that these should be 1/2n of the standard sersic index
    569 float indexGuess[] = {0.5, 0.33, 0.25, 0.167, 0.125, 0.083};
    570 # define N_INDEX_GUESS 6
     581// float indexGuess[] = {0.8, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 5.5, 6.0};
     582// float indexGuess[] = {0.5, 0.33, 0.25, 0.167, 0.125, 0.083};
     583float indexGuess[] = {1.0, 2.0, 3.0, 4.0};
     584# define N_INDEX_GUESS 4
    571585
    572586// A sersic model is very sensitive to the index.  attempt to find the index first by grid search in just the index
    … …  
    586600    float chiSquare[N_INDEX_GUESS];
    587601
     602# define TEST_X -540.0
     603# define TEST_Y 540.0
     604   
     605    if ((fabs(source->peak->xf - TEST_X) < 5) && (fabs(source->peak->yf - TEST_Y) < 5)) {
     606        psTraceSetLevel("psModules.objects.pmPCM_MinimizeChisq", 5);
     607    }
     608
    588609    for (int i = 0; i < N_INDEX_GUESS; i++) {
    589         model->params->data.F32[PM_PAR_7] = indexGuess[i];
     610        model->params->data.F32[PM_PAR_7] = 0.5/indexGuess[i];
    590611
    591612        if (!model->modelGuess(model, source)) {
    … …  
    594615        }
    595616
    596         // each time we change the model guess, we need to adjust the radius
    597         // XXX this did not work : we do not need such a large radius -- just uses moments-based radius
    598         if (false && !psphotSetRadiusModel (model, readout, source, markVal, false)) {
    599             psphotSetRadiusFootprint(&model->fitRadius, readout, source, markVal, 0.5);
    600         }
    601        
    602617        pmSourceFitModel (source, model, &options, maskVal);
    603         // fprintf (stderr, "chisq: %f, nIter: %d, radius: %f, npix: %d\n", model->chisqNorm, model->nIter, model->fitRadius, model->nPix);
     618        // fprintf (stderr, "index: %f, chisq: %f, nIter: %d, radius: %f, npix: %d\n", indexGuess[i], model->chisqNorm, model->nIter, model->fitRadius, model->nPix);
    604619
    605620        chiSquare[i] = model->chisqNorm;
    … …  
    616631    assert (iMin >= 0);
    617632
     633    if ((fabs(source->peak->xf - TEST_X) < 5) && (fabs(source->peak->yf - TEST_Y) < 5)) {
     634        psTraceSetLevel("psModules.objects.pmPCM_MinimizeChisq", 0);
     635    }
     636
    618637    model->flags = PM_MODEL_STATUS_NONE; // do not attempt to handle failures here, let the next iteration deal with it
    619     model->params->data.F32[PM_PAR_7] = indexGuess[iMin];
     638    model->params->data.F32[PM_PAR_7] = 0.5/indexGuess[iMin];
    620639    model->modelGuess(model, source);
    621 
    622     // each time we change the model guess, we need to adjust the radius
    623     // if (!psphotSetRadiusModel (model, readout, source, markVal, true)) {
    624     //  psphotSetRadiusFootprint(&model->fitRadius, readout, source, markVal);
    625     // }
    626640
    627641    return true;
    … …  
    648662    float xMin = NAN;
    649663    float chiSquare[N_INDEX_GUESS];
     664
     665    if ((fabs(source->peak->xf - TEST_X) < 5) && (fabs(source->peak->yf - TEST_Y) < 5)) {
     666        psTraceSetLevel("psModules.objects.pmPCM_MinimizeChisq", 5);
     667    }
    650668
    651669    for (int i = 0; i < N_INDEX_GUESS; i++) {
    … …  
    657675        }
    658676
    659         // each time we change the model guess, we need to adjust the radius
    660         // XXX this did not work : we do not need such a large radius -- just uses moments-based radius
    661         if (false && !psphotSetRadiusModel (model, readout, source, markVal, false)) {
    662             psphotSetRadiusFootprint(&model->fitRadius, readout, source, markVal, 0.5);
    663         }
    664        
    665677# if (0)
     678        // this block is to test the relative speed of straight and PCM fits
    666679        pmSourceFitModel (source, model, &options, maskVal);
    667680# else
    … …  
    669682        pmSourceFitPCM (pcm, source, &options, maskVal, markVal, psfSize);
    670683# endif
     684        fprintf (stderr, "index: %f, chisq: %f, nIter: %d, radius: %f, npix: %d\n", indexGuess[i], model->chisqNorm, model->nIter, model->fitRadius, model->nPix);
    671685        // fprintf (stderr, "chisq: %f, nIter: %d, radius: %f, npix: %d\n", model->chisqNorm, model->nIter, model->fitRadius, model->nPix);
    672686
    … …  
    684698    assert (iMin >= 0);
    685699   
     700    if ((fabs(source->peak->xf - TEST_X) < 5) && (fabs(source->peak->yf - TEST_Y) < 5)) {
     701        psTraceSetLevel("psModules.objects.pmPCM_MinimizeChisq", 0);
     702    }
     703
    686704    model->flags = PM_MODEL_STATUS_NONE; // do not attempt to handle failures here, let the next iteration deal with it
    687705    model->params->data.F32[PM_PAR_7] = indexGuess[iMin];

  • trunk/psphot/src/psphotSourceMatch.c

    r31154 r32348  
    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        }

  • trunk/psphot/src/psphotSourceSize.c

    r31452 r32348  
    4545    bool status = true;
    4646
     47    fprintf (stdout, "\n");
     48    psLogMsg ("psphot", PS_LOG_INFO, "--- psphot Source Size ---");
     49
    4750    // select the appropriate recipe information
    4851    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
    … …  
    209212        pmSourceMagnitudes (source, psf, photMode, maskVal, markVal, source->apRadius);
    210213
    211         float kMag = -2.5*log10(source->moments->KronFlux);
     214        float kMag = -2.5*log10(source->moments->KronFluxPSF);
    212215        float dMag = source->psfMag - kMag;
    213216
    … …  
    334337        psF32 Mxy = source->moments->Mxy;
    335338
    336         float KronMag = -2.5*log10(source->moments->KronFlux);
     339        float KronMag = -2.5*log10(source->moments->KronFluxPSF);
    337340
    338341        float Mminor = 0.5*(Mxx + Myy) - 0.5*sqrt(PS_SQR(Mxx - Myy) + 4.0*PS_SQR(Mxy));
    … …  
    509512        pmSourceSub (source, PM_MODEL_OP_FULL, options->maskVal);
    510513
    511         float kMag = -2.5*log10(source->moments->KronFlux);
     514        float kMag = -2.5*log10(source->moments->KronFluxPSF);
    512515        float dMag = source->psfMag - kMag;
    513516

  • trunk/psphot/src/psphotSourceStats.c

    r31673 r32348  
    99{
    1010    bool status = true;
     11
     12    fprintf (stdout, "\n");
     13    psLogMsg ("psphot", PS_LOG_INFO, "--- psphot Source Stats ---");
    1114
    1215    // select the appropriate recipe information
    … …  
    239242    psFree (cellGroups);
    240243
    241     psLogMsg ("psphot", PS_LOG_INFO, "%ld sources, %d moments, %d faint, %d failed: %f sec\n", sources->n, Nmoments, Nfaint, Nfail, psTimerMark ("psphot.stats"));
     244    psLogMsg ("psphot", PS_LOG_WARN, "%ld sources, %d moments, %d faint, %d failed: %f sec\n", sources->n, Nmoments, Nfaint, Nfail, psTimerMark ("psphot.stats"));
    242245
    243246    psphotVisualShowMoments (sources);
    … …  
    554557        // determine the PSF parameters from the source moment values
    555558        pmPSFClump psfClump = pmSourcePSFClump (NULL, NULL, sources, PSF_SN_LIM, PSF_CLUMP_GRID_SCALE, MOMENTS_SX_MAX, MOMENTS_SY_MAX, MOMENTS_SX_MIN, MOMENTS_SY_MIN, MOMENTS_AR_MAX);
    556         psLogMsg ("psphot", 3, "radius %.1f, nStars: %d of %d in clump, nSigma: %5.2f, X,  Y: %f, %f (%f, %f)\n", sigma[i], psfClump.nStars, psfClump.nTotal, psfClump.nSigma, psfClump.X, psfClump.Y, sqrt(psfClump.X) / sigma[i], sqrt(psfClump.Y) / sigma[i]);
     559        psLogMsg ("psphot", 3, "sigma guess (pix) %.1f, nStars: %d of %d in clump, nSigma: %5.2f, X,  Y: %f, %f (%f, %f)\n", sigma[i], psfClump.nStars, psfClump.nTotal, psfClump.nSigma, psfClump.X, psfClump.Y, sqrt(psfClump.X) / sigma[i], sqrt(psfClump.Y) / sigma[i]);
    557560
    558561        Rmin[i] = pmSourceMinKronRadius(sources, PSF_SN_LIM);

  • trunk/psphot/src/psphotStackImageLoop.c

    r31154 r32348  
    5252                psMemDump("load");
    5353
    54                 // PSF matching
     54                // Generate the 1st PSF-matched image set (larger target PSFs are generated by smoothing this image)
    5555                if (!psphotStackMatchPSFs (config, view)) {
    5656                    psError(psErrorCodeLast(), false, "failure in psphotStackMatchPSFs for chip %d, cell %d, readout %d\n", view->chip, view->cell, view->readout);

  • trunk/psphot/src/psphotStackMatchPSFs.c

    r31154 r32348  
    2525    }
    2626
    27     // loop over the available readouts (ignore chisq image)
     27    // loop over the available readouts
    2828    for (int i = 0; i < num; i++) {
     29        // set up the PSF-matching parameters describing the input images
    2930        if (!psphotStackMatchPSFsPrepare (config, view, options, i)) {
    3031            psError (PSPHOT_ERR_CONFIG, false, "failed to set PSF matching options for entry %d", i);
    … …  
    3334    }
    3435
    35     // Generate target PSF
     36    // XXX convolve == false might not be valid at the moment
    3637    if (options->convolve) {
    37         options->psf = psphotStackPSF(config, options->numCols, options->numRows, options->psfs, options->inputMask);
    38         if (!options->psf) {
     38        // Determine the 1st target PSF (either AUTO or defined by PSPHOT.STACK.TARGET.PSF.FWHM)
     39        // NOTE: this also set the full list of target FWHMs (options->targetSeeing)
     40        if (!psphotStackPSF(config, options)) {
    3941            psError(psErrorCodeLast(), false, "Unable to determine output PSF.");
    4042            return false;
    4143        }
    42         psMetadataAddPtr(config->arguments, PS_LIST_TAIL, "PSF.TARGET", PS_DATA_UNKNOWN, "Target PSF for stack", options->psf);
    43         options->targetSeeing = pmPSFtoFWHM(options->psf, 0.5 * options->numCols, 0.5 * options->numRows); // FWHM for target
    44         psLogMsg("psphotStack", PS_LOG_INFO, "Target seeing FWHM: %f\n", options->targetSeeing);
    45 
    46         // XXX is this needed to supply the psf to psphot??
    47         // pmChip *outChip = pmFPAfileThisChip(config->files, view, "PPSTACK.TARGET.PSF"); // Output chip
    48         // psMetadataAddPtr(outChip->analysis, PS_LIST_TAIL, "PSPHOT.PSF", PS_DATA_UNKNOWN, "Target PSF", options->psf);
    49         // outChip->data_exists = true;
    5044    }
    5145
    … …  
    6357
    6458// convolve the image to match desired PSF
     59// XXX is this code consistent with 'convolve' = false?
    6560bool psphotStackMatchPSFsReadout (pmConfig *config, const pmFPAview *view, psphotStackOptions *options, int index) {
    6661
    … …  
    109104        saveMatchData(readoutOut, options, index);
    110105    }
    111     rescaleData(readoutOut, config, options, index);
     106
     107    // renormalize the stack variances to have sigma = 1.0
     108    if (!psphotStackRenormaliseVariance(config, readoutOut)) {
     109        psError(psErrorCodeLast(), false, "Unable to renormalise variance.");
     110        return false;
     111    }
    112112
    113113    // save the output fwhm values in the readout->analysis.  we may have / will have multiple output PSF sizes,
    114114    // so we save this in a vector.  if the vector is not yet defined, create it
    115     bool mdok = false;
    116     psVector *fwhmValues = psMetadataLookupVector(&mdok, readoutOut->analysis, "STACK.PSF.FWHM.VALUES");
    117     if (!fwhmValues) {
    118         fwhmValues = psVectorAllocEmpty(10, PS_TYPE_F32);
    119         psMetadataAddVector(readoutOut->analysis, PS_LIST_TAIL, "STACK.PSF.FWHM.VALUES", PS_META_REPLACE, "PSF sizes", fwhmValues);
    120         psFree(fwhmValues); // drops the extra copy
     115    // NOTE: fwhmValues as defined here has 1 + nMatched PSF : 0 == unmatched
     116    psVector *fwhmValues = psVectorAllocEmpty(10, PS_TYPE_F32);
     117    psVectorAppend(fwhmValues, NAN); // XXX this corresponds to the unmatched image set
     118    for (int i = 0; i < options->targetSeeing->n; i++) {
     119        psVectorAppend(fwhmValues, options->targetSeeing->data.F32[i]);
    121120    }
    122     psVectorAppend(fwhmValues, options->targetSeeing);
     121    psMetadataAddVector(readoutSrc->analysis, PS_LIST_TAIL, "STACK.PSF.FWHM.VALUES", PS_META_REPLACE, "PSF sizes", fwhmValues);
     122    psMetadataAddVector(readoutOut->analysis, PS_LIST_TAIL, "STACK.PSF.FWHM.VALUES", PS_META_REPLACE, "PSF sizes", fwhmValues);
     123    psFree(fwhmValues); // drops the extra copy
    123124
    124125    return true;

  • trunk/psphot/src/psphotStackMatchPSFsNext.c

    r30624 r32348  
    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");

  • trunk/psphot/src/psphotStackMatchPSFsPrepare.c

    r28013 r32348  
    44bool determineSeeing (pmPSF *psf, psphotStackOptions *options, int index);
    55
    6 // set up the stacking parameters
     6// set up the PSF-matching parameters describing the input images
    77bool psphotStackMatchPSFsPrepare (pmConfig *config, const pmFPAview *view, psphotStackOptions *options, int index) {
    88

  • trunk/psphot/src/psphotStackMatchPSFsUtils.c

    r31154 r32348  
    22# define ARRAY_BUFFER 16                 // Number to add to array at a time
    33
    4 // XXX better name
    5 bool readImage(psImage **target, // Target for image
    6                const char *name, // Name of FITS file
    7                const pmConfig *config // Configuration
    8     )
    9 {
    10     psString resolved = pmConfigConvertFilename(name, config, false, false); // Resolved filename
    11     psFits *fits = psFitsOpen(resolved, "r");
    12     psFree(resolved);
    13     if (!fits) {
    14         psError(PSPHOT_ERR_IO, false, "Unable to open previously produced image: %s", name);
    15         return false;
    16     }
    17     psImage *image = psFitsReadImage(fits, psRegionSet(0,0,0,0), 0); // Image of interest
    18     if (!image) {
    19         psError(PSPHOT_ERR_IO, false, "Unable to read previously produced image: %s", name);
    20         psFitsClose(fits);
    21         return false;
    22     }
    23     psFitsClose(fits);
    24 
    25     psFree(*target);
    26     *target = image;
    27 
    28     return true;
    29 }
     4psVector *SetOptWidths (bool *optimum, psMetadata *recipe);
     5pmReadout *makeFakeReadout(pmConfig *config, pmReadout *raw, psArray *sources, pmPSF *psf, psImageMaskType maskVal, int fullSize);
     6bool saveMatchData (pmReadout *readout, psphotStackOptions *options, int index);
     7bool matchKernel(pmConfig *config, pmReadout *cnv, pmReadout *raw, psphotStackOptions *options, int index);
     8bool dumpImageDiff(pmReadout *readoutConv, pmReadout *readoutFake, pmReadout *readoutRef, int index, char *rootname);
     9bool dumpImage(pmReadout *readoutOut, pmReadout *readoutRef, int index, char *rootname);
    3010
    3111// Get coordinates from a source
    … …  
    148128
    149129// Renormalise a readout's variance map
    150 bool stackRenormaliseReadout(const pmConfig *config, // Configuration
    151                              pmReadout *readout      // Readout to renormalise
     130bool psphotStackRenormaliseVariance(const pmConfig *config, // Configuration
     131                              pmReadout *readout      // Readout to renormalise
    152132    )
    153133{
    … …  
    180160    return pmReadoutVarianceRenormalise(readout, maskBad, num, minValid, maxValid);
    181161}
    182 
    183 // This is a hack to use the temporary convolved images and kernel generated previously.
    184 // This makes the 'matching' operation much faster, allowing debugging of the stack process easier.
    185 // It implicitly assumes the output root name is the same between invocations.
    186 
    187 # if (0)
    188 bool loadKernel (pmConfig *config, pmReadout *readoutCnv, psphotStackOptions *options, int index) {
    189 
    190     // Read the convolution kernel from the saved file
    191     pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPSTACK.CONV.KERNEL", index);
    192     psAssert(file, "Require file");
    193 
    194     pmFPAview *view = pmFPAviewAlloc(0); // View to readout of interest
    195     view->chip = view->cell = view->readout = 0;
    196     psString filename = pmFPAfileNameFromRule(file->filerule, file, view); // Filename of interest
    197 
    198     // Read convolution kernel data
    199     psString resolved = pmConfigConvertFilename(filename, config, false, false); // Resolved filename
    200     psFree(filename);
    201     psFits *fits = psFitsOpen(resolved, "r"); // FITS file for subtraction kernel
    202     psFree(resolved);
    203     if (!fits || !pmReadoutReadSubtractionKernels(readoutCnv, fits)) {
    204         psError(PSPHOT_ERR_IO, false, "Unable to read previously produced kernel");
    205         psFitsClose(fits);
    206         return false;
    207     }
    208     psFitsClose(fits);
    209 
    210     // read the convolved pixels (image, mask, variance) -- names are pre-defined
    211     if (!readImage(&readoutCnv->image,    options->convImages->data[index],    config) ||
    212         !readImage(&readoutCnv->mask,     options->convMasks->data[index],     config) ||
    213         !readImage(&readoutCnv->variance, options->convVariances->data[index], config)) {
    214         psError(PSPHOT_ERR_IO, false, "Unable to read previously produced image.");
    215         return false;
    216     }
    217 
    218     // XXX ??? not sure what is happening here -- consult Paul
    219     psRegion *region = psMetadataLookupPtr(NULL, readoutCnv->analysis, PM_SUBTRACTION_ANALYSIS_REGION); // Convolution region
    220     pmSubtractionKernels *kernels = psMetadataLookupPtr(NULL, readoutCnv->analysis, PM_SUBTRACTION_ANALYSIS_KERNEL);
    221 
    222     pmSubtractionAnalysis(readoutCnv->analysis, NULL, kernels, region, readoutCnv->image->numCols, readoutCnv->image->numRows);
    223 
    224     psKernel *kernel = pmSubtractionKernel(kernels, 0.0, 0.0, false); // Convolution kernel
    225 
    226     // update the covariance matrix
    227     // XXX why is this needed if we have correctly read the saved data?
    228     bool oldThreads = psImageCovarianceSetThreads(true);              // Old thread setting
    229     psKernel *covar = psImageCovarianceCalculate(kernel, readoutCnv->covariance); // Covariance matrix
    230     psImageCovarianceSetThreads(oldThreads);
    231     psFree(readoutCnv->covariance);
    232     readoutCnv->covariance = covar;
    233     psFree(kernel);
    234     return true;
    235 }
    236 # endif
    237162
    238163bool dumpImage(pmReadout *readoutOut, pmReadout *readoutRef, int index, char *rootname) {
    … …  
    362287        widthsCopy = psVectorCopy(NULL, widths, PS_TYPE_F32); // Copy of kernel widths
    363288
    364         // we need to register the FWHM values for use downstream
    365         pmSubtractionSetFWHMs(options->targetSeeing, options->inputSeeing->data.F32[index]);
     289        // we need to register the FWHM values for use by pmSubtraction code
     290        pmSubtractionSetFWHMs(options->targetSeeing->data.F32[0], options->inputSeeing->data.F32[index]);
    366291
    367292        pmSubtractionParamScaleOptions(scale, scaleRef, scaleMin, scaleMax);
    … …  
    469394}
    470395
    471 // Kernel normalisation for convolved readout
    472 bool renormKernel(pmReadout *readout, psphotStackOptions *options, int index) {
    473 
    474     double sum = 0.0;           // Sum of chi^2
    475     int num = 0;                // Number of measurements of chi^2
    476     psString regex = NULL;      // Regular expression
    477     psStringAppend(&regex, "^%s$", PM_SUBTRACTION_ANALYSIS_NORM);
    478     psMetadataIterator *iter = psMetadataIteratorAlloc(readout->analysis, PS_LIST_HEAD, regex);
    479     psFree(regex);
    480     psMetadataItem *item = NULL;// Item from iteration
    481     while ((item = psMetadataGetAndIncrement(iter))) {
    482         assert(item->type == PS_TYPE_F32);
    483         float norm = item->data.F32; // Normalisation
    484         sum += norm;
    485         num++;
    486     }
    487     psFree(iter);
    488     float conv = sum/num;       // Mean normalisation from convolution
    489     float stars = powf(10.0, -0.4 * options->norm->data.F32[index]); // Normalisation from stars
    490     float renorm =  stars / conv; // Renormalisation to apply
    491     psLogMsg("psphotStack", PS_LOG_INFO, "Renormalising image %d by %f (kernel: %f, stars: %f)\n", index, renorm, conv, stars);
    492 
    493     psBinaryOp(readout->image, readout->image, "*", psScalarAlloc(renorm, PS_TYPE_F32));
    494     psBinaryOp(readout->variance, readout->variance, "*", psScalarAlloc(PS_SQR(renorm), PS_TYPE_F32));
    495     return true;
    496 }
    497 
    498 // adjust scaling for readout (remove background, ..., determine weighting)
    499 bool rescaleData(pmReadout *readout, pmConfig *config, psphotStackOptions *options, int index) {
    500 
    501     psMetadata *stackRecipe = psMetadataLookupMetadata(NULL, config->recipes, "PPSTACK"); // ppStack recipe
    502     psAssert(stackRecipe, "We've thrown an error on this before.");
    503 
    504     // Look up appropriate values from the ppSub recipe
    505     psMetadata *subRecipe = psMetadataLookupMetadata(NULL, config->recipes, "PPSUB"); // PPSUB recipe
    506     psAssert(subRecipe, "recipe missing");
    507 
    508     psString maskValStr = psMetadataLookupStr(NULL, subRecipe, "MASK.VAL"); // Name of bits to mask going in
    509     psString maskBadStr = psMetadataLookupStr(NULL, stackRecipe, "MASK.BAD"); // Name of bits to mask for bad
    510 
    511     psImageMaskType maskVal = pmConfigMaskGet(maskValStr, config); // Bits to mask going in to pmSubtractionMatch
    512     psImageMaskType maskBad = pmConfigMaskGet(maskBadStr, config); // Bits to mask for bad pixels
    513 
    514     // Ensure the background value is zero
    515     psStats *bg = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); // Statistics for background
    516     psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); // Random number generator
    517 
    518     // XXX why is this in config->arguments and not recipe?
    519     if (!psMetadataLookupBool(NULL, config->arguments, "PPSTACK.SKIP.BG.SUB")) {
    520         if (!psImageBackground(bg, NULL, readout->image, readout->mask, maskVal | maskBad, rng)) {
    521             psAbort("Can't measure background for image.");
    522             // XXX we used to clear error: why is this acceptable? psErrorClear();
    523         }
    524 
    525         float value = psStatsGetValue(bg, PS_STAT_ROBUST_MEDIAN);
    526         float stdev = psStatsGetValue(bg, PS_STAT_ROBUST_STDEV);
    527 
    528         psLogMsg("psphotStack", PS_LOG_INFO, "Correcting convolved image background by %lf (+/- %lf)", value, stdev);
    529         psBinaryOp(readout->image, readout->image, "-", psScalarAlloc(value, PS_TYPE_F32));
    530     }
    531 
    532     if (!stackRenormaliseReadout(config, readout)) {
    533         psFree(rng);
    534         psFree(bg);
    535         return false;
    536     }
    537 
    538     // Measure the variance level for the weighting
    539     if (psMetadataLookupBool(NULL, stackRecipe, "WEIGHTS")) {
    540         if (!psImageBackground(bg, NULL, readout->variance, readout->mask, maskVal | maskBad, rng)) {
    541             psError(PSPHOT_ERR_DATA, false, "Can't measure mean variance for image.");
    542             psFree(rng);
    543             psFree(bg);
    544             return false;
    545         }
    546         options->weightings->data.F32[index] = 1.0 / (psStatsGetValue(bg, PS_STAT_ROBUST_MEDIAN) * psImageCovarianceFactor(readout->covariance));
    547     } else {
    548         options->weightings->data.F32[index] = 1.0;
    549     }
    550     psLogMsg("psphotStack", PS_LOG_INFO, "Weighting for image %d is %f\n", index, options->weightings->data.F32[index]);
    551 
    552     psFree(rng);
    553     psFree(bg);
    554     return true;
    555 }
    556 
    557396# define NOISE_FRACTION 0.01             // Set minimum flux to this fraction of noise
    558397# define SOURCE_MASK (PM_SOURCE_MODE_FAIL | PM_SOURCE_MODE_DEFECT | PM_SOURCE_MODE_SATURATED | PM_SOURCE_MODE_CR_LIMIT | PM_SOURCE_MODE_EXT_LIMIT) // Mask to apply to input sources

  • trunk/psphot/src/psphotStackObjects.c

    r28013 r32348  
    4747    return true;
    4848}
     49
     50// mark good vs bad objects
     51bool psphotStackObjectsSelectForAnalysis (pmConfig *config, const pmFPAview *view, const char *filerule, psArray *objects) {
     52
     53    bool status = false;
     54
     55    // find the currently selected readout
     56    pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, 0); // use the 0-index image to represent the image area
     57    psAssert (file, "missing file?");
     58
     59    pmReadout *readout = pmFPAviewThisReadout(view, file->fpa);
     60    psAssert (readout, "missing readout?");
     61
     62    // select the appropriate recipe information
     63    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
     64    psAssert (recipe, "missing recipe?");
     65
     66    // option to limit analysis to a specific region
     67    char *region = psMetadataLookupStr (&status, recipe, "ANALYSIS_REGION");
     68    psRegion AnalysisRegion = psRegionForImage (readout->image, psRegionFromString (region));
     69    if (psRegionIsNaN (AnalysisRegion)) psAbort("analysis region mis-defined");
     70
     71    // S/N limit to perform full non-linear fits
     72    float SN_LIM = psMetadataLookupF32 (&status, recipe, "EXTENDED_SOURCE_SN_LIM");
     73
     74    bool doPetroStars   = psMetadataLookupBool (&status, recipe, "PETROSIAN_FOR_STARS");
     75
     76    for (int i = 0; i < objects->n; i++) {
     77        pmPhotObj *object = objects->data[i];
     78        if (!object) continue;
     79        if (!object->sources) continue;
     80
     81        // we check each source for an object and keep the object if any source is valid
     82
     83        bool keepObject = false;
     84        for (int j = 0; j < object->sources->n; j++) {
     85
     86            pmSource *source = object->sources->data[j];
     87            if (!source) continue;
     88            if (!source->peak) continue;
     89
     90            // skip PSF-like and non-astronomical objects
     91            if (source->type == PM_SOURCE_TYPE_DEFECT) continue;
     92            if (source->type == PM_SOURCE_TYPE_SATURATED) continue;
     93            if (source->mode & PM_SOURCE_MODE_DEFECT) continue;
     94            if (source->mode & PM_SOURCE_MODE_SATSTAR) continue;
     95           
     96        // optionally allow non-extended objects to get petrosians as well
     97        if (!doPetroStars) {
     98            if (!(source->mode & PM_SOURCE_MODE_EXT_LIMIT)) continue;
     99            if (source->type == PM_SOURCE_TYPE_STAR) continue;
     100        }
     101
     102            // limit selection to some SN limit
     103            // assert (source->peak); // how can a source not have a peak?
     104            // limit selection to some SN limit
     105            bool skipSource = false;
     106            if (source->mode & PM_SOURCE_MODE_EXT_LIMIT) {
     107                skipSource = (source->moments->KronFlux < SN_LIM * source->moments->KronFluxErr);
     108            } else {
     109                skipSource = (sqrt(source->peak->detValue) < SN_LIM);
     110            }
     111            if (skipSource) continue;
     112
     113            // limit selection by analysis region (this automatically apply
     114            if (source->peak->x < AnalysisRegion.x0) continue;
     115            if (source->peak->y < AnalysisRegion.y0) continue;
     116            if (source->peak->x > AnalysisRegion.x1) continue;
     117            if (source->peak->y > AnalysisRegion.y1) continue;
     118           
     119            keepObject = true;
     120        }
     121
     122        for (int j = 0; j < object->sources->n; j++) {
     123            pmSource *source = object->sources->data[j];
     124            if (!source) continue;
     125            if (!source->peak) continue;
     126
     127            // we have to set a bit in either case to tell psphotExtendedSourceAnalysis to
     128            // avoid the single-detection tests
     129
     130            if (keepObject) {
     131                source->tmpFlags |=  PM_SOURCE_TMPF_STACK_KEEP;
     132                source->tmpFlags &= ~PM_SOURCE_TMPF_STACK_SKIP;
     133            } else {
     134                source->tmpFlags |=  PM_SOURCE_TMPF_STACK_SKIP;
     135                source->tmpFlags &= ~PM_SOURCE_TMPF_STACK_KEEP;
     136            }       
     137        }
     138    }
     139
     140    psLogMsg ("psphot", PS_LOG_INFO, "marked good vs bad objects\n");
     141    return true;
     142}

  • trunk/psphot/src/psphotStackOptions.c

    r28160 r32348  
    1919    psFree (options->norm);
    2020    psFree (options->matchChi2);
    21     psFree (options->weightings);
     21    psFree (options->targetSeeing);
    2222
    2323    return;
    … …  
    3636    options->convolve = false;
    3737    options->convolveSource = PSPHOT_CNV_SRC_NONE;
    38     options->targetSeeing = NAN;
     38    options->targetSeeing = NULL;
    3939
    4040    options->psfs        = psArrayAlloc(num);
    … …  
    4747    options->norm        = psVectorAlloc(num, PS_TYPE_F32);
    4848    options->matchChi2   = psVectorAlloc(num, PS_TYPE_F32); // chi^2 for stamps when matching
    49     options->weightings  = psVectorAlloc(num, PS_TYPE_F32); // Combination weightings for images (1/noise^2)
    5049
    5150    psVectorInit(options->inputMask,   0);
    … …  
    5352    psVectorInit(options->norm,        NAN);
    5453    psVectorInit(options->matchChi2,   NAN);
    55     psVectorInit(options->weightings,  NAN);
    5654
    5755    return options;

  • trunk/psphot/src/psphotStackPSF.c

    r31154 r32348  
    11# include "psphotInternal.h"
    22
    3 pmPSF *psphotStackPSF(const pmConfig *config, int numCols, int numRows, const psArray *psfs, const psVector *inputMask)
    4 {
     3// determine the 1st target PSF (either AUTO or defined by PSPHOT.STACK.TARGET.PSF.FWHM)
     4bool psphotStackPSF(const pmConfig *config, psphotStackOptions *options)  {
     5
    56    bool mdok = false;
    6     pmPSF *psf = NULL;
     7
     8    int numCols = options->numCols;
     9    int numRows = options->numRows;
     10    psArray *psfs = options->psfs;
     11    psVector *inputMask = options->inputMask;
    712
    813    // Get the recipe values
    … …  
    1722       
    1823    char *psfModel = psMetadataLookupStr(NULL, recipe, "PSF.MODEL"); // Model for PSF
     24
     25    options->targetSeeing = psVectorAllocEmpty(4, PS_TYPE_F32);
    1926
    2027    if (autoPSF) {
    … …  
    3845
    3946        // Solve for the target PSF
    40         psf = pmPSFEnvelope(numCols, numRows, psfs, psfInstances, psfRadius, psfModel, psfOrder, psfOrder, maskVal);
    41         if (!psf) {
     47        options->psf = pmPSFEnvelope(numCols, numRows, psfs, psfInstances, psfRadius, psfModel, psfOrder, psfOrder, maskVal);
     48        if (!options->psf) {
    4249            psError(PSPHOT_ERR_PSF, false, "Unable to determine output PSF.");
    43             return NULL;
     50            return false;
    4451        }
    4552
    46     } else {
    47 
    48         // externally-defined PSF
    49         // XXX need to test for compatibility of target with inputs
    50 
    51         float targetFWHM = psMetadataLookupF32 (&mdok, psphotRecipe, "PSPHOT.STACK.TARGET.PSF.FWHM");
    52         if (!mdok) {
    53             psVector *fwhmValues = psMetadataLookupVector(&mdok, psphotRecipe, "PSPHOT.STACK.TARGET.PSF.FWHM"); // Magnitude offsets
    54             psAssert (mdok, "missing psphot recipe value PSPHOT.STACK.TARGET.PSF.FWHM");
    55             targetFWHM = fwhmValues->data.F32[0];
    56         }
    57 
    58         // measured scale factors (fwhm = Sxx * 2.35 * scaleFactor / sqrt(2.0))
    59         // GAUSS  : 1.000
    60         // PGAUSS : 1.006
    61         // QGAUSS : 1.151
    62         // RGAUSS : 0.883
    63         // PS1_V1 : 1.134
    64        
    65         float scaleFactor = NAN;
    66         if (!strcmp(psfModel, "PS_MODEL_GAUSS")) {
    67             scaleFactor = 1.000;
    68         }
    69         if (!strcmp(psfModel, "PS_MODEL_PGAUSS")) {
    70             scaleFactor = 1.0006;
    71         }
    72         if (!strcmp(psfModel, "PS_MODEL_QGAUSS")) {
    73             scaleFactor = 1.151;
    74         }
    75         if (!strcmp(psfModel, "PS_MODEL_RGAUSS")) {
    76             scaleFactor = 0.883;
    77         }
    78         if (!strcmp(psfModel, "PS_MODEL_PS1_V1")) {
    79             scaleFactor = 1.134;
    80         }
    81         psAssert (isfinite(scaleFactor), "invalid model for PSF");
    82 
    83         float Sxx = sqrt(2.0)*targetFWHM / 2.35 / scaleFactor;
    84 
    85         // XXX probably should make the model type (and par 7) optional from recipe
    86         // psf = pmPSFBuildSimple(psfModel, Sxx, Sxx, 0.0, 1.0);
    87         psf = pmPSFBuildSimple(psfModel, Sxx, Sxx, 0.0, 0.2);
    88         if (!psf) {
    89             psError(PSPHOT_ERR_PSF, false, "Unable to build dummy PSF.");
    90             return NULL;
    91         }
     53        psMetadataAddPtr(config->arguments, PS_LIST_TAIL, "PSF.TARGET", PS_DATA_UNKNOWN, "Target PSF for stack", options->psf);
     54        float targetSeeing = pmPSFtoFWHM(options->psf, 0.5 * options->numCols, 0.5 * options->numRows); // FWHM for target
     55        psVectorAppend(options->targetSeeing, targetSeeing);
     56        psLogMsg("psphotStack", PS_LOG_INFO, "Target seeing FWHM (auto-scaled): %f\n", targetSeeing);
     57        return true;
    9258    }
    9359
    94     return psf;
     60    // externally-defined PSF
     61    // XXX need to test for compatibility of target with inputs
     62
     63    // is a single target FWHM specified, or a set of values?  set up the vector options->targetSeeing and the local 1st value
     64    float targetSeeing = psMetadataLookupF32 (&mdok, psphotRecipe, "PSPHOT.STACK.TARGET.PSF.FWHM");
     65    if (!mdok) {
     66        psVector *fwhmValues = psMetadataLookupVector(&mdok, psphotRecipe, "PSPHOT.STACK.TARGET.PSF.FWHM"); // Magnitude offsets
     67        psAssert (mdok, "missing psphot recipe value PSPHOT.STACK.TARGET.PSF.FWHM");
     68        for (int i = 0; i < fwhmValues->n; i++) {
     69            psVectorAppend(options->targetSeeing, fwhmValues->data.F32[i]);
     70        }           
     71        targetSeeing = fwhmValues->data.F32[0];
     72    } else {
     73        psVectorAppend(options->targetSeeing, targetSeeing);
     74    }
     75
     76    // measured scale factors (fwhm = Sxx * 2.35 * scaleFactor / sqrt(2.0))
     77    // GAUSS  : 1.000
     78    // PGAUSS : 1.006
     79    // QGAUSS : 1.151
     80    // RGAUSS : 0.883
     81    // PS1_V1 : 1.134
     82       
     83    float scaleFactor = NAN;
     84    if (!strcmp(psfModel, "PS_MODEL_GAUSS")) {
     85        scaleFactor = 1.000;
     86    }
     87    if (!strcmp(psfModel, "PS_MODEL_PGAUSS")) {
     88        scaleFactor = 1.0006;
     89    }
     90    if (!strcmp(psfModel, "PS_MODEL_QGAUSS")) {
     91        scaleFactor = 1.151;
     92    }
     93    if (!strcmp(psfModel, "PS_MODEL_RGAUSS")) {
     94        scaleFactor = 0.883;
     95    }
     96    if (!strcmp(psfModel, "PS_MODEL_PS1_V1")) {
     97        scaleFactor = 1.134;
     98    }
     99    psAssert (isfinite(scaleFactor), "invalid model for PSF");
     100
     101    float Sxx = sqrt(2.0)*targetSeeing / 2.35 / scaleFactor;
     102
     103    // XXX probably should make the model type (and par 7) optional from recipe
     104    // psf = pmPSFBuildSimple(psfModel, Sxx, Sxx, 0.0, 1.0);
     105    options->psf = pmPSFBuildSimple(psfModel, Sxx, Sxx, 0.0, 0.2);
     106    if (!options->psf) {
     107        psError(PSPHOT_ERR_PSF, false, "Unable to build dummy PSF.");
     108        return false;
     109    }
     110
     111    psMetadataAddPtr(config->arguments, PS_LIST_TAIL, "PSF.TARGET", PS_DATA_UNKNOWN, "Target PSF for stack", options->psf);
     112    psLogMsg("psphotStack", PS_LOG_INFO, "Target seeing FWHM (1 of %ld): %f\n", options->targetSeeing->n, targetSeeing);
     113    return true;
    95114}

  • trunk/psphot/src/psphotStackReadout.c

    r31154 r32348  
    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();

  • trunk/psphot/src/psphotSubtractBackground.c

    r31154 r32348  
    103103        psphotSaveImage (NULL, image, name);
    104104    }
    105     psLogMsg ("psphot", PS_LOG_INFO, "subtracted background model: %f sec\n", psTimerMark ("psphot.background"));
     105    psLogMsg ("psphot", PS_LOG_WARN, "subtracted background model: %f sec\n", psTimerMark ("psphot.background"));
    106106
    107107    // the pmReadout selected in this function are all view on entries in config->files
Note: See TracChangeset for help on using the changeset viewer.