IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 31452 for trunk


Ignore:
Timestamp:
May 5, 2011, 11:09:38 AM (15 years ago)
Author:
eugene
Message:

merge updates from eam branch 20110404

Location:
trunk/psphot
Files:
29 edited
1 copied

Legend:

Unmodified
Added
Removed

  • trunk/psphot

  • trunk/psphot/src/Makefile.am

    r31154 r31452  
    167167        psphotSourcePlots.c            \
    168168        psphotRadialPlot.c             \
     169        psphotKronMasked.c             \
    169170        psphotDeblendSatstars.c        \
    170171        psphotMosaicSubimage.c         \

  • trunk/psphot/src/psphot.h

    r31154 r31452  
    149149bool            psphotPSFstats (pmReadout *readout, pmPSF *psf);
    150150bool            psphotMomentsStats (pmReadout *readout, psArray *sources);
     151
     152bool            psphotKronMasked (pmConfig *config, const pmFPAview *view, const char *filerule);
     153bool            psphotKronMaskedReadout(pmConfig *config, psMetadata *recipe, const pmFPAview *view, pmReadout *readout, psArray *sources, pmPSF *psf);
    151154
    152155// in psphotGuessModel.c
    … …  
    310313
    311314bool psphotMakeFluxScale (psImage *image, psMetadata *recipe, pmPSF *psf);
    312 bool psphotMakeGrowthCurve (pmReadout *readout, psMetadata *recipe, pmPSF *psf);
     315bool psphotMakeGrowthCurve (pmReadout *readout, psMetadata *recipe, pmPSF *psf, psArray *sources);
    313316bool psphotDumpPSFStars (pmReadout *readout, pmPSFtry *try, float radius, psImageMaskType maskVal, psImageMaskType markVal);
    314317

  • trunk/psphot/src/psphotApResid.c

    r31154 r31452  
    11# include "psphotInternal.h"
    2 // # define DEBUG
     2# define DEBUG
    33
    44# define SKIPSTAR(MSG) { psTrace ("psphot", 3, "invalid : %s", MSG); continue; }
    … …  
    206206
    207207        // XXX make this user-configurable?
    208         if (source->errMag > 0.01) continue;
     208        if (source->psfMagErr > 0.01) continue;
    209209
    210210        // aperture residual for this source
    … …  
    222222                 source->peak->xf, source->peak->yf,
    223223                 source->modelPSF->params->data.F32[PM_PAR_XPOS], source->modelPSF->params->data.F32[PM_PAR_YPOS],
    224                  source->psfMag, source->apMag, source->errMag,
     224                 source->psfMag, source->apMag, source->psfMagErr,
    225225                 source->modelPSF->params->data.F32[PM_PAR_I0],
    226226                 source->modelPSF->params->data.F32[PM_PAR_SXX], source->modelPSF->params->data.F32[PM_PAR_SXY], source->modelPSF->params->data.F32[PM_PAR_SYY],
    … …  
    228228# endif
    229229        if (!isfinite(source->psfMag)) psAbort ("nan in psfMag");
    230         if (!isfinite(source->errMag)) psAbort ("nan in errMag");
     230        if (!isfinite(source->psfMagErr)) psAbort ("nan in psfMagErr");
    231231        if (!isfinite(source->apMag)) psAbort ("nan in apMag");
    232232        if (!isfinite(model->params->data.F32[PM_PAR_XPOS])) psAbort ("nan in xPos");
    … …  
    234234
    235235        psVectorAppend (mag, source->psfMag);
    236         psVectorAppend (dMag,source->errMag);
     236        psVectorAppend (dMag,source->psfMagErr);
    237237        psVectorAppend (apResid, dap);
    238238        psVectorAppend (xPos, model->params->data.F32[PM_PAR_XPOS]);
    … …  
    328328    psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "DAPMIFIT", PS_META_REPLACE, "ap residual scatter", psf->dApResid);
    329329    psMetadataAddS32 (readout->analysis, PS_LIST_TAIL, "NAPMIFIT", PS_META_REPLACE, "number of apresid stars", psf->nApResid);
    330     psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "APLOSS",   PS_META_REPLACE, "aperture loss (mag)", psf->growth ? psf->growth->apLoss : NAN);
     330
     331    // curve-of-growth offset
     332    if (psf->growth) {
     333        float gaussSigma = psMetadataLookupF32(&status, readout->analysis, "MOMENTS_GAUSS_SIGMA");
     334        if (!status) {
     335            gaussSigma = psMetadataLookupF32(&status, recipe, "MOMENTS_GAUSS_SIGMA");
     336        }
     337        float apScale = psMetadataLookupF32(&status, recipe, "PSF_APERTURE_SCALE");
     338        float PSF_APERTURE = (int)(apScale*gaussSigma);
     339       
     340        float offset = pmGrowthCurveCorrect (psf->growth, PSF_APERTURE);
     341        psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "APLOSS",   PS_META_REPLACE, "aperture loss (mag)", psf->growth->apLoss);
     342        psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "APREFOFF", PS_META_REPLACE, "offset to ref aperture", offset);
     343    } else {
     344        psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "APLOSS",   PS_META_REPLACE, "aperture loss (mag)", NAN);
     345        psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "APREFOFF", PS_META_REPLACE, "offset to ref aperture", NAN);
     346    }
    331347
    332348    psLogMsg ("psphot.apresid", PS_LOG_DETAIL, "aperture residual: %f +/- %f\n", psf->ApResid, psf->dApResid);

  • trunk/psphot/src/psphotBlendFit.c

    r31154 r31452  
    242242        if (source->mode &  PM_SOURCE_MODE_PAIR) continue;
    243243
     244        // do not include MOMENTS_FAILURES in the fit
     245        if (source->mode & PM_SOURCE_MODE_MOMENTS_FAILURE) continue;
     246
    244247        // limit selection to some SN limit
    245248        if (sqrt(source->peak->detValue) < FIT_SN_LIM) continue;
    … …  
    280283        if (source->mode & PM_SOURCE_MODE_EXT_LIMIT) {
    281284            if (psphotFitBlob (readout, source, newSources, psf, fitOptions, maskVal, markVal)) {
    282                 source->type = PM_SOURCE_TYPE_EXTENDED;
    283285                psTrace ("psphot", 5, "source at %7.1f, %7.1f is ext", source->peak->xf, source->peak->yf);
    284286                Next ++;
    285                 source->mode |= PM_SOURCE_MODE_NONLINEAR_FIT;
    286287                continue;
    287288            }
    … …  
    291292                psTrace ("psphot", 5, "source at %7.1f, %7.1f is psf", source->peak->xf, source->peak->yf);
    292293                Npsf ++;
    293                 source->mode |= PM_SOURCE_MODE_NONLINEAR_FIT;
    294294                continue;
    295295            }

  • trunk/psphot/src/psphotChoosePSF.c

    r31154 r31452  
    350350
    351351    // build curve-of-growth vector for this psf
    352     if (!psphotMakeGrowthCurve (readout, recipe, try->psf)) {
     352    if (!psphotMakeGrowthCurve (readout, recipe, try->psf, try->sources)) {
    353353        psError(PSPHOT_ERR_PSF, false, "Unable to construct flux scale for PSF");
    354354        psFree (models);

  • trunk/psphot/src/psphotEfficiency.c

    r30560 r31452  
    400400                source->type = PM_SOURCE_TYPE_STAR;
    401401
     402                source->modelPSF->fitRadius = sourceRadius;
     403                source->apRadius = sourceRadius;
     404
    402405                numFound++;
    403406                psArrayAdd(sources, sources->n, source);
    … …  
    468471                    source->modelPSF->params->data.F32[PM_PAR_XPOS],
    469472                    source->modelPSF->params->data.F32[PM_PAR_YPOS],
    470                     magRef, source->psfMag, source->errMag);
     473                    magRef, source->psfMag, source->psfMagErr);
    471474#endif
    472475            magDiff->data.F32[j] = source->psfMag - magRef;
    473             magErr->data.F32[j] = source->errMag;
     476            magErr->data.F32[j] = source->psfMagErr;
    474477        }
    475478        magDiff->n = sources->n;

  • trunk/psphot/src/psphotExtendedSourceFits.c

    r31154 r31452  
    333333        // this uses the footprint to judge both radius and aperture?
    334334        // XXX save the psf-based moments for output
    335         if (!pmSourceMoments (source, radius, 0.0, 0.0, maskVal)) {
     335        if (!pmSourceMoments (source, radius, 0.0, 0.0, 0.0, maskVal)) {
    336336            fprintf (stderr, "skipping (2) %f, %f\n", source->peak->xf, source->peak->yf);
    337337            // subtract the best fit from the object, leave local sky
    … …  
    379379          assert (status);
    380380
    381           // fprintf (stderr, "xfit: %f, %f : %f\n", source->peak->xf, source->peak->yf, source->peak->SN);
    382 
    383381          // limit selection to some SN limit
    384           assert (source->peak); // how can a source not have a peak?
     382          // assert (source->peak); // how can a source not have a peak?
    385383          if (sqrt(source->peak->detValue) < SNlim) {
    386384              Nfaint ++;

  • trunk/psphot/src/psphotFitSourcesLinear.c

    r31154 r31452  
    137137        if (source->mode & PM_SOURCE_MODE_CR_LIMIT) continue;
    138138
     139        // do not include MOMENTS_FAILURES in the fit
     140        if (source->mode & PM_SOURCE_MODE_MOMENTS_FAILURE) continue;
     141
    139142        // XXX count saturated stars
    140143        if (source->mode & PM_SOURCE_MODE_SATSTAR) {
    … …  
    164167
    165168        // check the integral of the model : is it large enough?
     169        // apply mask?
    166170        float modelSum = 0.0;
     171        float maskedSum = 0.0;
    167172        for (int iy = 0; iy < source->modelFlux->numRows; iy++) {
    168173            for (int ix = 0; ix < source->modelFlux->numCols; ix++) {
    169174                modelSum += source->modelFlux->data.F32[iy][ix];
     175                if (source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & maskVal) continue;
     176                maskedSum += source->modelFlux->data.F32[iy][ix];
    170177            }
    171178        }
    172179        if (modelSum < 0.5) continue; // skip sources with no model constraint (somewhat arbitrary limit)
    173         // if (modelSum < 0.01) continue; // skip sources with no model constraint (somewhat arbitrary limit)
    174180        if (modelSum < 0.8) {
    175             fprintf (stderr, "low-sig model @ %f, %f (%f sum, %f peak)\n",
    176                      source->peak->xf, source->peak->yf, modelSum, source->peak->rawFlux);
     181            fprintf (stderr, "low-sig model @ %f, %f (%f masked sum, %f sum, %f peak)\n",
     182                     source->peak->xf, source->peak->yf, maskedSum, modelSum, source->peak->rawFlux);
    177183        }
    178 
    179         pmModel *model = pmSourceGetModel (NULL, source);
     184        if (maskedSum < 0.5) {
     185            fprintf (stderr, "worrying model @ %f, %f (%f masked sum, %f sum, %f peak)\n",
     186                     source->peak->xf, source->peak->yf, maskedSum, modelSum, source->peak->rawFlux);
     187        }
     188
     189        bool isPSF = false;
     190        pmModel *model = pmSourceGetModel (&isPSF, source);
    180191
    181192        psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal));
    182193        psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, model->fitRadius, "OR", markVal);
    183194
    184         source->mode |= PM_SOURCE_MODE_LINEAR_FIT;
     195        // we call this function multiple times. for the first time, we have only PSF models for all objects
     196        // the second time has extended sources.  If we ever fit the PSF model, we should raise this bit
     197        source->mode |= PM_SOURCE_MODE_LINEAR_FIT;
     198        if (isPSF) {
     199            source->mode |= PM_SOURCE_MODE_PSFMODEL;
     200        }           
     201       
    185202        psArrayAdd (fitSources, 100, source);
    186203    }
    … …  
    260277        }
    261278    }
     279
     280# if (0)
     281    static int npass = 0;
     282    char name[128];
     283    FILE *f1 = NULL;
     284    int fd = -1;
     285
     286    snprintf (name, 128, "sparse.Aij.%02d.dat", npass);
     287    f1 = fopen (name, "w");
     288    psAssert (f1, "failed to open file\n");
     289    fd = fileno (f1);
     290    p_psVectorPrint (fd, sparse->Aij, "Aij");
     291    fclose (f1);
     292
     293    snprintf (name, 128, "sparse.Bfj.%02d.dat", npass);
     294    f1 = fopen (name, "w");
     295    psAssert (f1, "failed to open file\n");
     296    fd = fileno (f1);
     297    p_psVectorPrint (fd, sparse->Bfj, "Bfj");
     298    fclose (f1);
     299
     300    snprintf (name, 128, "sparse.Qii.%02d.dat", npass);
     301    f1 = fopen (name, "w");
     302    psAssert (f1, "failed to open file\n");
     303    fd = fileno (f1);
     304    p_psVectorPrint (fd, sparse->Qii, "Qii");
     305    fclose (f1);
     306    npass ++;
     307# endif
    262308
    263309    psSparseResort (sparse);

  • trunk/psphot/src/psphotFitSourcesLinearStack.c

    r31154 r31452  
    6969            if (!pmSourceCacheModel (source, maskVal)) continue;
    7070          }
     71
     72          // check the integral of the model : is it large enough?
     73          float modelSum = 0.0;
     74          for (int iy = 0; iy < source->modelFlux->numRows; iy++) {
     75              for (int ix = 0; ix < source->modelFlux->numCols; ix++) {
     76                  modelSum += source->modelFlux->data.F32[iy][ix];
     77              }
     78          }
     79          if (modelSum < 0.5) continue; // skip sources with no model constraint (somewhat arbitrary limit)
     80          if (modelSum < 0.8) {
     81              fprintf (stderr, "low-sig model @ %f, %f (%f sum, %f peak)\n",
     82                       source->peak->xf, source->peak->yf, modelSum, source->peak->rawFlux);
     83          }
    7184
    7285          source->mode |= PM_SOURCE_MODE_LINEAR_FIT;

  • trunk/psphot/src/psphotGuessModels.c

    r31154 r31452  
    160160        pmSource *source = sources->data[i];
    161161
    162         if (source->mode & PM_SOURCE_MODE_MOMENTS_FAILURE) {
    163             fprintf (stderr, "moment failure\n");
    164         }
    165 
    166162        // this is used to mark sources for which the model is measured. We check later that
    167163        // all are used.
    … …  
    169165
    170166        // skip non-astronomical objects (very likely defects)
     167        if (source->mode & PM_SOURCE_MODE_MOMENTS_FAILURE) continue;
    171168        if (source->type == PM_SOURCE_TYPE_DEFECT) continue;
    172169        if (source->type == PM_SOURCE_TYPE_SATURATED) continue;
    … …  
    196193        }
    197194
     195# if (0)
    198196        if (source->mode & PM_SOURCE_MODE_SATSTAR) {
    199197            fprintf (stderr, "satstar: %f,%f vs %f,%f : %c\n",
    … …  
    202200                     (useMoments ? 'T' : 'F'));
    203201        }
     202# endif
    204203
    205204        // set PSF parameters for this model (apply 2D shape model to coordinates Xo, Yo)

  • trunk/psphot/src/psphotLoadSRCTEXT.c

    r31154 r31452  
    5252
    5353            // NOTE: most of these values are irrelevant for loaded source positions
    54             source->seq       = 0;
     54            source->seq       = source->id;
    5555            PAR[PM_PAR_XPOS]  = X;
    5656            PAR[PM_PAR_YPOS]  = Y;
    … …  
    6767
    6868            source->psfMag    = 0.0;
    69             source->errMag    = 0.0;
     69            source->psfMagErr    = 0.0;
    7070            source->apMag     = 0.0;
    7171

  • trunk/psphot/src/psphotMagnitudes.c

    r31154 r31452  
    167167        pmSource *source = (pmSource *) sources->data[i];
    168168
     169        bool saveTest = false;
     170        psImage *testImage = NULL;
     171# if (0)
     172        if ((fabs(source->peak->xf-3518) < 5) && (fabs(source->peak->yf-3178) < 5)) {
     173            saveTest = true;
     174            psRegion subregion = psRegionSet (source->peak->xf - 200, source->peak->xf + 200, source->peak->yf - 200, source->peak->yf + 200);
     175            testImage = psImageSubset ((psImage *) source->pixels->parent, subregion);
     176        }
     177# endif
     178
     179        if (saveTest) {
     180            psphotSaveImage(NULL, testImage, "test.image.1.fits");
     181        }
     182
    169183        // replace object in image
    170184        if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) {
    … …  
    172186        }
    173187
     188        if (saveTest) {
     189            psphotSaveImage(NULL, testImage, "test.image.2.fits");
     190        }
     191
     192        float xPos = source->peak->xf;
     193        float yPos = source->peak->yf;
     194
     195        pmModel *model = source->modelPSF;
     196        if (model) {
     197            xPos = model->params->data.F32[PM_PAR_XPOS];
     198            yPos = model->params->data.F32[PM_PAR_YPOS];
     199        } else {
     200            bool useMoments = pmSourcePositionUseMoments(source);
     201            if (useMoments) {
     202                xPos = source->moments->Mx;
     203                yPos = source->moments->My;
     204            }
     205        }
     206
    174207        // clear the mask bit and set the circular mask pixels
    175208        psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal));
    176         psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, source->apRadius, "OR", markVal);
     209        psImageKeepCircle (source->maskObj, xPos, yPos, source->apRadius, "OR", markVal);
    177210
    178211        status = pmSourceMagnitudes (source, psf, photMode, maskVal, markVal, source->apRadius);
    179         if (status && isfinite(source->apMag)) Nap ++;
     212        if (status && isfinite(source->apFlux)) {
     213            Nap ++;
     214        } else {
     215            fprintf (stderr, "failed to measure mag for source @ %f,%f\n", source->peak->xf, source->peak->yf);
     216        }
    180217
    181218        // clear the mask bit
    … …  
    185222        pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
    186223
     224        if (saveTest) {
     225            psphotSaveImage(NULL, testImage, "test.image.3.fits");
     226        }
     227
    187228        if (backModel) {
    188229            psAssert (binning, "if backModel is defined, so should binning be");
    189             source->sky = psImageUnbinPixel(source->peak->x, source->peak->y, backModel->image, binning);
     230            source->sky = psImageUnbinPixel(xPos, yPos, backModel->image, binning);
    190231            if (isnan(source->sky) && false) {
    191232                psLogMsg ("psphot.magnitudes", PS_LOG_DETAIL, "error setting pmSource.sky");
    … …  
    197238        if (backStdev) {
    198239            psAssert (binning, "if backStdev is defined, so should binning be");
    199             source->skyErr = psImageUnbinPixel(source->peak->x, source->peak->y, backStdev->image, binning);
     240            source->skyErr = psImageUnbinPixel(xPos, yPos, backStdev->image, binning);
    200241            if (isnan(source->skyErr) && false) {
    201242                psLogMsg ("psphot.magnitudes", PS_LOG_DETAIL, "error setting pmSource.skyErr");

  • trunk/psphot/src/psphotMakeGrowthCurve.c

    r31154 r31452  
    11# include "psphotInternal.h"
    22
    3 bool psphotMakeGrowthCurve (pmReadout *readout, psMetadata *recipe, pmPSF *psf) {
     3bool psphotMakeGrowthCurve (pmReadout *readout, psMetadata *recipe, pmPSF *psf, psArray *sources) {
    44
    55    bool status;
    … …  
    2222
    2323    // measure the aperture loss as a function of radius for PSF
    24     bool IGNORE_GROWTH = psMetadataLookupBool (&status, recipe, "IGNORE_GROWTH");
    2524    float REF_RADIUS = psMetadataLookupF32 (&status, recipe, "PSF_REF_RADIUS");
    2625    float PSF_FIT_PAD   = psMetadataLookupF32 (&status, recipe, "PSF_FIT_PADDING");
     26   
     27    float gaussSigma = psMetadataLookupF32(&status, readout->analysis, "MOMENTS_GAUSS_SIGMA");
     28    if (!status) {
     29        gaussSigma = psMetadataLookupF32(&status, recipe, "MOMENTS_GAUSS_SIGMA");
     30    }
     31    float apScale = psMetadataLookupF32(&status, recipe, "PSF_APERTURE_SCALE");
     32    float PSF_APERTURE = (int)(apScale*gaussSigma);
    2733
    2834    psf->growth = pmGrowthCurveAlloc (PSF_FIT_PAD, 100.0, REF_RADIUS);
    2935
    30     if (!pmGrowthCurveGenerate (readout, psf, IGNORE_GROWTH, maskVal, markVal)) {
    31         psError(PSPHOT_ERR_APERTURE, false, "Fitting aperture corrections");
    32         psFree(psf->growth); psf->growth = NULL;
    33         return false;
     36    bool GROWTH_FROM_SOURCES = psMetadataLookupBool (&status, recipe, "GROWTH_FROM_SOURCES");
     37    if (!status) GROWTH_FROM_SOURCES = false;
     38
     39    if (GROWTH_FROM_SOURCES) {
     40        bool INTERPOLATE_AP = psMetadataLookupBool (&status, recipe, "INTERPOLATE_AP");
     41        if (!pmGrowthCurveGenerateFromSources (readout, psf, sources, INTERPOLATE_AP, maskVal, markVal)) {
     42            psError(PSPHOT_ERR_APERTURE, false, "Fitting aperture corrections");
     43            psFree(psf->growth); psf->growth = NULL;
     44            return false;
     45        }
     46
     47    } else {
     48        bool IGNORE_GROWTH = psMetadataLookupBool (&status, recipe, "IGNORE_GROWTH");
     49        if (!pmGrowthCurveGenerate (readout, psf, IGNORE_GROWTH, maskVal, markVal)) {
     50            psError(PSPHOT_ERR_APERTURE, false, "Fitting aperture corrections");
     51            psFree(psf->growth); psf->growth = NULL;
     52            return false;
     53        }
    3454    }
    3555
    3656    psLogMsg ("psphot", PS_LOG_MINUTIA, "built growth curve: %f sec\n", psTimerMark ("psphot.growth"));
    3757
     58    float offset = pmGrowthCurveCorrect (psf->growth, PSF_APERTURE);
     59    psLogMsg ("psphot", PS_LOG_DETAIL, "correction from %f to %f: %f mags\n", PSF_APERTURE, REF_RADIUS, offset);
     60
    3861    return true;
    3962}

  • trunk/psphot/src/psphotMergeSources.c

    r31154 r31452  
    137137                pmSource *source = extSourcesTXT->data[i];
    138138                source->mode |= PM_SOURCE_MODE_EXTERNAL;
     139
     140                // the supplied peak flux needs to be re-normalized
     141                source->peak->rawFlux = 1.0;
     142                source->peak->smoothFlux = 1.0;
     143                source->peak->detValue = 1.0;
    139144
    140145                // drop the loaded source modelPSF

  • trunk/psphot/src/psphotModelTest.c

    r29548 r31452  
    226226    // measure the source mags
    227227    pmSourcePhotometryModel (&fitMag, model);
    228     pmSourcePhotometryAper  (&obsMag, NULL, NULL, model, source->pixels, source->variance, source->maskObj, maskVal);
     228    pmSourcePhotometryAper  (NULL, &obsMag, NULL, NULL, model, source->pixels, source->variance, source->maskObj, maskVal);
    229229    fprintf (stderr, "ap: %f, fit: %f, apmifit: %f, nIter: %d\n", obsMag, fitMag, obsMag - fitMag, model->nIter);
    230230

  • trunk/psphot/src/psphotOutput.c

    r31154 r31452  
    7575                 source->peak->xf, source->peak->yf,
    7676                 model->params->data.F32[PM_PAR_XPOS], model->params->data.F32[PM_PAR_YPOS],
    77                  source->psfMag, source->apMag, source->errMag,
     77                 source->psfMag, source->apMag, source->psfMagErr,
    7878                 model->params->data.F32[PM_PAR_I0],
    7979                 model->params->data.F32[PM_PAR_SXX], model->params->data.F32[PM_PAR_SXY], model->params->data.F32[PM_PAR_SYY],
    … …  
    251251    psMetadataItemSupplement (&status, header, analysis, "NPSFSTAR");
    252252    psMetadataItemSupplement (&status, header, analysis, "APLOSS");
     253    psMetadataItemSupplement (&status, header, analysis, "APREFOFF");
    253254
    254255    psMetadataItemSupplement (&status, header, analysis, "FWHM_MAJ");

  • trunk/psphot/src/psphotPetrosianStats.c

    r31154 r31452  
    163163        if (refRadius->data.F32[i] > apRadius) {
    164164            if (i == 0) {
    165                 psWarning ("does this case make any sense? (refRadius[0] > apRadius)");
     165                psWarning ("does this case make any sense? (refRadius[0] %f > apRadius %f)", refRadius->data.F32[i], apRadius);
    166166                continue;
    167167            } else {
    … …  
    188188        if (!found50 && (fluxSum->data.F32[i] > flux50)) {
    189189            if (i == 0) {
    190                 psWarning ("does this case make any sense? (fluxSum[0] > flux50)");
     190                psWarning ("does this case make any sense? (fluxSum[0] %f > flux50 %f)", fluxSum->data.F32[i], flux50);
    191191                continue;
    192192            } else {
    … …  
    198198        if (!found90 && (fluxSum->data.F32[i] > flux90)) {
    199199            if (i == 0) {
    200                 psWarning ("does this case make any sense? (fluxSum[0] > flux90)");
     200                psWarning ("does this case make any sense? (fluxSum[0] %f > flux90 %f)", fluxSum->data.F32[i], flux90);
    201201                continue;
    202202            } else {

  • trunk/psphot/src/psphotRadialApertures.c

    r31154 r31452  
    9595        if (source->type == PM_SOURCE_TYPE_SATURATED) continue;
    9696        if (source->mode & PM_SOURCE_MODE_DEFECT) continue;
    97         if (source->mode & PM_SOURCE_MODE_SATSTAR) continue;
     97
     98        // XXX measure radial apertures even for saturated stars
     99        // if (source->mode & PM_SOURCE_MODE_SATSTAR) continue;
    98100
    99101        // limit selection to some SN limit
    … …  
    257259        float SBstdv = sqrt((fluxStd->data.F32[i] / nPix) - PS_SQR(SBmean));
    258260
     261        // XXX report the total flux or the mask-corrected flux?
    259262        // flux->data.F32[i]    = SBmean * Area;
     263        // fluxErr->data.F32[i] = sqrt(fluxErr->data.F32[i]) * Area / nPix;
     264
     265        fluxErr->data.F32[i] = sqrt(fluxErr->data.F32[i]);
    260266        fluxStd->data.F32[i] = SBstdv * Area;
    261         fluxErr->data.F32[i] = sqrt(fluxErr->data.F32[i]) * Area / nPix;
    262 
    263         // fill->data.F32[i] /= Area;
     267        fill->data.F32[i] /= Area;
     268
    264269        psTrace ("psphot", 5, "radial bins: %3d  %5.1f : %8.1f +/- %7.2f : %8.1f +/- %8.1f : %4.2f %6.1f\n",
    265270                 i, aperRadii->data.F32[i], flux->data.F32[i], fluxErr->data.F32[i], SBmean, SBstdv, fill->data.F32[i], Area);

  • trunk/psphot/src/psphotRadiusChecks.c

    r29936 r31452  
    184184bool psphotSetRadiusModel (pmModel *model, pmReadout *readout, pmSource *source, psImageMaskType markVal, bool deep) {
    185185
    186     psF32 *PAR = model->params->data.F32;
    187 
    188     pmMoments *moments = source->moments;
    189     if (moments == NULL) return false;
     186    pmPeak *peak = source->peak;
    190187
    191188    // set the fit radius based on the object flux limit and the model
    … …  
    199196
    200197    // redefine the pixels if needed
    201     pmSourceRedefinePixels (source, readout, PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], model->fitRadius);
    202 
    203     // set the mask to flag the excluded pixels
    204     psImageKeepCircle (source->maskObj, PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], model->fitRadius, "OR", markVal);
    205     return true;
    206 }
     198    pmSourceRedefinePixels (source, readout, peak->xf, peak->yf, model->fitRadius);
     199
     200    // set the mask to flag the excluded pixels
     201    psImageKeepCircle (source->maskObj, peak->xf, peak->yf, model->fitRadius, "OR", markVal);
     202    return true;
     203}

  • trunk/psphot/src/psphotReadout.c

    r31154 r31452  
    200200    psphotDumpChisqs (config, view, filerule);
    201201
     202    // XXX re-measure the kron mags with models subtracted
     203    psphotKronMasked(config, view, filerule);
     204
    202205    // identify CRs and extended sources (only unmeasured sources are measured)
    203206    psphotSourceSize (config, view, filerule, true); // pass 1 (detections->allSources)
    … …  
    309312pass1finish:
    310313
     314    // XXX re-measure the kron mags with models subtracted
     315    psphotKronMasked(config, view, filerule);
     316
    311317    // measure source size for the remaining sources
    312318    // NOTE: applies only to NEW (unmeasured) sources
    … …  
    345351      psErrorStackPrint(stderr, "Unable to replace sky");
    346352      psErrorClear();
    347 
    348 /*       psLogMsg("psphot", 3, "failed on psphotSkyReplace"); */
    349 /*       return(psphotReadoutCleanup(config, view, filerule)); */
    350     }
     353    }
     354
    351355    // drop the references to the image pixels held by each source
    352356    if (!psphotSourceFreePixels (config, view, filerule)) { // pass 1
    353357      psErrorStackPrint(stderr, "Unable to free source pixels");
    354358      psErrorClear();
    355 
    356 /*       psLogMsg ("psphot", 3, "failed on psphotSourceFreePixels"); */
    357 /*       return(psphotReadoutCleanup(config, view, filerule)); */
    358     }
     359    }
     360
    359361    // create the exported-metadata and free local data
    360362    return psphotReadoutCleanup(config, view, filerule);

  • trunk/psphot/src/psphotReplaceUnfit.c

    r31154 r31452  
    6767    psAssert (maskVal, "missing mask value?");
    6868
     69    // bit-mask to mark pixels not used in analysis
     70    psImageMaskType markVal = psMetadataLookupImageMask(&status, recipe, "MARK.PSPHOT");
     71    assert (markVal);
     72
     73    // maskVal is used to test for rejected pixels, and must include markVal
     74    maskVal |= markVal;
     75
    6976    for (int i = 0; i < sources->n; i++) {
    7077      source = sources->data[i];

  • trunk/psphot/src/psphotRoughClass.c

    r31154 r31452  
    207207        return false;
    208208    }
    209     psLogMsg ("psphot", 3, "psf clump  X,  Y: %f, %f : DX, DY: %f, %f : nStars %d of %d\n", psfClump.X, psfClump.Y, psfClump.dX, psfClump.dY, psfClump.nStars, psfClump.nTotal);
     209    if (!havePSF) {
     210        psLogMsg ("psphot", 3, "psf clump  X,  Y: %f, %f : DX, DY: %f, %f : nStars %d of %d\n", psfClump.X, psfClump.Y, psfClump.dX, psfClump.dY, psfClump.nStars, psfClump.nTotal);
     211    } else {
     212        psLogMsg ("psphot", 3, "psf clump  X,  Y: %f, %f : DX, DY: %f, %f : loaded from metadata\n", psfClump.X, psfClump.Y, psfClump.dX, psfClump.dY);
     213    }
    210214
    211215    // get basic parameters, or set defaults

  • trunk/psphot/src/psphotSetThreads.c

    r31154 r31452  
    2525    psFree(task);
    2626
    27     task = psThreadTaskAlloc("PSPHOT_SOURCE_STATS", 10);
     27    task = psThreadTaskAlloc("PSPHOT_SOURCE_STATS", 11);
    2828    task->function = &psphotSourceStats_Threaded;
    2929    psThreadTaskAdd(task);

  • trunk/psphot/src/psphotSourceFits.c

    r31154 r31452  
    146146    pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
    147147    source->mode |=  PM_SOURCE_MODE_BLEND_FIT;
     148    source->mode |= PM_SOURCE_MODE_NONLINEAR_FIT;
    148149    return true;
    149150}
    … …  
    186187
    187188    // build cached model and subtract
     189    source->mode |= PM_SOURCE_MODE_NONLINEAR_FIT;
    188190    pmSourceCacheModel (source, maskVal);
    189191    pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
    … …  
    233235    // XXX save the psf-based moments for output
    234236    psfMoments = *source->moments;
    235     if (!pmSourceMoments (source, radius, 0.0, 0.5, maskVal)) {
     237    if (!pmSourceMoments (source, radius, 0.0, 0.5, 0.0, maskVal)) {
    236238      *source->moments = psfMoments;
    237239      return false;
    … …  
    319321    source->type = PM_SOURCE_TYPE_EXTENDED;
    320322    source->mode |= PM_SOURCE_MODE_EXTMODEL;
     323    source->mode |= PM_SOURCE_MODE_NONLINEAR_FIT;
    321324
    322325    // build cached model and subtract
    … …  
    346349    source->modelPSF = psMemIncrRefCounter (DBL->data[0]);
    347350    source->mode     |= PM_SOURCE_MODE_PAIR;
     351    source->mode |= PM_SOURCE_MODE_NONLINEAR_FIT;
    348352
    349353    // copy most data from the primary source (modelEXT, blends stay NULL)
    … …  
    489493    // psTraceSetLevel("psLib.math.psMinimizeLMChi2", 5);
    490494    pmSourceFitModel (source, model, &options, maskVal);
    491     // fprintf (stderr, "chisq: %f, nIter: %d, radius: %f, npix: %d\n\n", model->chisqNorm, model->nIter, model->fitRadius, model->nPix);
     495    // fprintf (stderr, "chisq: %f, nIter: %d, radius: %f, npix: %d\n", model->chisqNorm, model->nIter, model->fitRadius, model->nPix);
    492496
    493497    // psTraceSetLevel("psLib.math.psMinimizeLMChi2", 0);
    … …  
    554558    // psTraceSetLevel("psLib.math.psMinimizeLMChi2", 5);
    555559    pmSourceFitPCM (pcm, source, &options, maskVal, markVal, psfSize);
    556     // fprintf (stderr, "chisq: %f, nIter: %d, radius: %f, npix: %d\n\n", model->chisqNorm, model->nIter, model->fitRadius, model->nPix);
     560    // fprintf (stderr, "chisq: %f, nIter: %d, radius: %f, npix: %d\n", model->chisqNorm, model->nIter, model->fitRadius, model->nPix);
    557561
    558562    // psTraceSetLevel("psLib.math.psMinimizeLMChi2", 0);

  • trunk/psphot/src/psphotSourceSize.c

    r31154 r31452  
    213213
    214214        psVectorAppend (ApOff, dMag);
    215         psVectorAppend (ApErr, source->errMag);
     215        psVectorAppend (ApErr, source->psfMagErr);
    216216    }
    217217    if (num == 0) {
    … …  
    520520        // set nSigmaMAG to include both systematic and poisson error terms.  we include a hard
    521521        // floor on the Ap Sys Err (to be a bit generous).  XXX put the floor in the recipe...
    522         float nSigmaMAG = (dMag - options->ApResid) / hypot(source->errMag, hypot(options->ApSysErr, 0.02));
     522        float nSigmaMAG = (dMag - options->ApResid) / hypot(source->psfMagErr, hypot(options->ApSysErr, 0.02));
    523523        source->extNsigma = nSigmaMAG;
    524524

  • trunk/psphot/src/psphotSourceStats.c

    r31154 r31452  
    9898
    9999    // generate the array of sources, define the associated pixel
     100    bool firstPass = false;
    100101    if (!detections->newSources) {
    101102        detections->newSources = psArrayAllocEmpty (peaks->n);
     103        firstPass = true;
    102104    }
    103105    sources = detections->newSources;
    … …  
    125127        if (source->mode & PM_SOURCE_MODE_MOMENTS_FAILURE) {
    126128            fprintf (stderr, "moment failure\n");
     129        }
     130
     131        if (firstPass) {
     132            source->mode2 |= PM_SOURCE_MODE2_PASS1_SRC;
    127133        }
    128134
    … …  
    159165        RADIUS = psMetadataLookupF32 (&status, recipe, "PSF_MOMENTS_RADIUS");
    160166    }
     167    float MIN_KRON_RADIUS = psMetadataLookupF32 (&status, readout->analysis, "MOMENTS_MIN_KRON");
     168    if (!status) {
     169        MIN_KRON_RADIUS = 0.75*SIGMA;
     170    }
    161171
    162172    // threaded measurement of the source magnitudes
    … …  
    186196            PS_ARRAY_ADD_SCALAR(job->args, RADIUS,  PS_TYPE_F32);
    187197            PS_ARRAY_ADD_SCALAR(job->args, SIGMA,   PS_TYPE_F32);
     198            PS_ARRAY_ADD_SCALAR(job->args, MIN_KRON_RADIUS, PS_TYPE_F32);
    188199
    189200            PS_ARRAY_ADD_SCALAR(job->args, maskVal, PS_TYPE_IMAGE_MASK);
    … …  
    215226            } else {
    216227                psScalar *scalar = NULL;
    217                 scalar = job->args->data[7];
     228                scalar = job->args->data[8];
    218229                Nmoments += scalar->data.S32;
    219                 scalar = job->args->data[8];
     230                scalar = job->args->data[9];
    220231                Nfail += scalar->data.S32;
    221                 scalar = job->args->data[9];
     232                scalar = job->args->data[10];
    222233                Nfaint += scalar->data.S32;
    223234            }
    … …  
    231242
    232243    psphotVisualShowMoments (sources);
     244
     245    // clear the mark bits
     246    // psImageMaskPixels (readout->mask, "AND", PS_NOT_IMAGE_MASK(markVal));
    233247
    234248    if (detections->allSources) {
    … …  
    373387    float RADIUS                    = PS_SCALAR_VALUE(job->args->data[3],F32);
    374388    float SIGMA                     = PS_SCALAR_VALUE(job->args->data[4],F32);
    375 
    376     psImageMaskType maskVal         = PS_SCALAR_VALUE(job->args->data[5],PS_TYPE_IMAGE_MASK_DATA);
    377     psImageMaskType markVal         = PS_SCALAR_VALUE(job->args->data[6],PS_TYPE_IMAGE_MASK_DATA);
     389    float MIN_KRON_RADIUS           = PS_SCALAR_VALUE(job->args->data[5],F32);
     390
     391    psImageMaskType maskVal         = PS_SCALAR_VALUE(job->args->data[6],PS_TYPE_IMAGE_MASK_DATA);
     392    psImageMaskType markVal         = PS_SCALAR_VALUE(job->args->data[7],PS_TYPE_IMAGE_MASK_DATA);
    378393
    379394    // if no valid pixels, massive swing or very large Mrf, likely saturated source,
    … …  
    431446        if (!(source->peak->type == PM_PEAK_SUSPECT_SATURATION)) {
    432447            // measure basic source moments (no S/N clipping on input pixels)
    433             status = pmSourceMoments (source, RADIUS, SIGMA, 0.0, maskVal);
     448            status = pmSourceMoments (source, RADIUS, SIGMA, 0.0, MIN_KRON_RADIUS, maskVal);
    434449        } else {
    435450            // For saturated stars, choose a much larger box NOTE this is slightly sleazy, but
    … …  
    447462
    448463            psTrace ("psphot", 4, "retrying moments for %d, %d\n", source->peak->x, source->peak->y);
    449             status = pmSourceMoments (source, BIG_RADIUS, BIG_SIGMA, 0.0, maskVal);
     464            status = pmSourceMoments (source, BIG_RADIUS, BIG_SIGMA, 0.0, MIN_KRON_RADIUS, maskVal);
    450465            source->mode |= PM_SOURCE_MODE_BIG_RADIUS;
    451466        }
    … …  
    456471            continue;
    457472        }
     473
     474        // XXX test of masking neighbors when measureing moments (does this fail?)
     475        // psEllipseMoments moments = {source->moments->Mxx, source->moments->Myy, source->moments->Mxy};
     476        // psEllipseAxes axes = psEllipseMomentsToAxes(moments, 20.0);
     477        // psImageMaskCircle (source->maskView, source->peak->x, source->peak->y, 3.0*axes.major, "OR", markVal);
     478
    458479        Nmoments ++;
    459480
    … …  
    492513    float sigma[NSIGMA] = {1.0, 2.0, 3.0, 4.5, 6.0, 9.0, 12.0, 18.0};
    493514    float Sout[NSIGMA];
     515    float Rmin[NSIGMA];
    494516    int   Nout[NSIGMA]; // number of stars found in clump : use this to control the number of regions measured by psphotRoughClass
    495517
    … …  
    513535
    514536            // measure basic source moments (no S/N clipping on input pixels)
    515             status = pmSourceMoments (source, 4*sigma[i], sigma[i], 0.0, maskVal);
     537            status = pmSourceMoments (source, 4*sigma[i], sigma[i], 0.0, 0.0, maskVal);
    516538        }
    517539
    … …  
    524546        pmPSFClump psfClump = pmSourcePSFClump (NULL, NULL, sources, PSF_SN_LIM, PSF_CLUMP_GRID_SCALE, MOMENTS_SX_MAX, MOMENTS_SY_MAX, MOMENTS_AR_MAX);
    525547        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]);
     548
     549        Rmin[i] = pmSourceMinKronRadius(sources, PSF_SN_LIM);
    526550
    527551#if 0
    … …  
    549573    // we are looking for sigma for which Sout = 0.65 (or some other value)
    550574    int Nstars = 0;
     575    float minKronRadius = NAN;
    551576    float Sigma = NAN;
    552577    float minS = Sout[0];
    … …  
    565590        Sigma = sigma[i] + (0.65 - Sout[i])*(sigma[i+1] - sigma[i])/(Sout[i+1] - Sout[i]);
    566591        Nstars = 0.5*(Nout[i] + Nout[i+1]);
     592        minKronRadius = Rmin[i] + (0.65 - Sout[i])*(Rmin[i+1] - Rmin[i])/(Sout[i+1] - Sout[i]);
    567593    }
    568594    psAssert (isfinite(Sigma), "did we miss a case?");
    … …  
    575601    psMetadataAddF32(analysis, PS_LIST_TAIL, "PSF_MOMENTS_RADIUS", PS_META_REPLACE, "moments limit", 4.0*Sigma);
    576602    psMetadataAddF32(analysis, PS_LIST_TAIL, "PSF_CLUMP_NSTARS", PS_META_REPLACE, "number of stars in clump", Nstars);
    577 
    578     psLogMsg ("psphot", 3, "using window function with sigma = %f\n", Sigma);
     603    psMetadataAddF32(analysis, PS_LIST_TAIL, "MOMENTS_MIN_KRON", PS_META_REPLACE, "minimum Kron Radius", minKronRadius);
     604
     605    psLogMsg ("psphot", 3, "using window function with sigma = %f (Min Kron Radius = %f, Nstars = %d)\n", Sigma, minKronRadius, Nstars);
    579606    return true;
    580607}

  • trunk/psphot/src/psphotVisual.c

    r31154 r31452  
    18671867            continue;
    18681868        }
    1869         if (source->errMag > 0.1) {
     1869        if (source->psfMagErr > 0.1) {
    18701870            xLOW->data.F32[nLOW] = source->moments->Mxx;
    18711871            yLOW->data.F32[nLOW] = source->moments->Myy;
    … …  
    25912591        x->data.F32[n] = source->psfMag;
    25922592        y->data.F32[n] = dMag;
    2593         dy->data.F32[n] = source->errMag;
     2593        dy->data.F32[n] = source->psfMagErr;
    25942594        graphdata.xmin = PS_MIN(graphdata.xmin, x->data.F32[n]);
    25952595        graphdata.xmax = PS_MAX(graphdata.xmax, x->data.F32[n]);

  • trunk/psphot/test/tap_psphot_stackphot.pro

    r31154 r31452  
    11#!/usr/bin/env mana
    22# -*-sh-*-
     3
     4$KAPA = kapa -noX
    35
    46# config for ppImage to generate chip, mask, weight
Note: See TracChangeset for help on using the changeset viewer.