IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
May 3, 2010, 8:50:52 AM (16 years ago)
Author:
eugene
Message:

updates from trunk

Location:
branches/simtest_nebulous_branches
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • branches/simtest_nebulous_branches

  • branches/simtest_nebulous_branches/psModules

  • branches/simtest_nebulous_branches/psModules/src/objects/pmSourcePhotometry.c

    r21511 r27840  
    8484
    8585    // we must have a valid model
     86    // XXX allow aperture magnitudes for sources without a model
    8687    model = pmSourceGetModel (&isPSF, source);
    8788    if (model == NULL) {
    … …  
    9091    }
    9192
     93    // XXX handle negative flux, low-significance
    9294    if (model->dparams->data.F32[PM_PAR_I0] > 0) {
    9395        SN = model->params->data.F32[PM_PAR_I0] / model->dparams->data.F32[PM_PAR_I0];
    … …  
    101103
    102104    // measure PSF model photometry
    103     if (psf->FluxScale) {
     105    // XXX TEST: do not use flux scale
     106    if (0 && psf->FluxScale) {
    104107        // the source peak pixel is guaranteed to be on the image, and only minimally different from the source center
    105108        double fluxScale = pmTrend2DEval (psf->FluxScale, (float)source->peak->x, (float)source->peak->y);
    106         if (!isfinite(fluxScale)) {
    107             // XXX objects on the edge can be slightly outside -- if we get an
    108             // error, use the full fit.
    109             psErrorClear();
    110             status = pmSourcePhotometryModel (&source->psfMag, source->modelPSF);
    111         } else {
    112             if (isfinite(fluxScale) && (fluxScale > 0.0)) {
    113                 source->psfMag = -2.5*log10(fluxScale * source->modelPSF->params->data.F32[PM_PAR_I0]);
    114             } else {
    115                 source->psfMag = NAN;
    116             }
    117         }
     109        psAssert (isfinite(fluxScale), "how can the flux scale be invalid? source at %d, %d\n", source->peak->x, source->peak->y);
     110        psAssert (fluxScale > 0.0, "how can the flux scale be negative? source at %d, %d\n", source->peak->x, source->peak->y);
     111        source->psfFlux = fluxScale * source->modelPSF->params->data.F32[PM_PAR_I0];
     112        source->psfFluxErr = fluxScale * source->modelPSF->dparams->data.F32[PM_PAR_I0];
     113        source->psfMag = -2.5*log10(source->psfFlux);
    118114    } else {
    119         status = pmSourcePhotometryModel (&source->psfMag, source->modelPSF);
    120     }
    121 
    122     // if we have a collection of model fits, one of them is a pointer to modelEXT?
    123     // XXX not guaranteed
     115        status = pmSourcePhotometryModel (&source->psfMag, &source->psfFlux, source->modelPSF);
     116        source->psfFluxErr = source->psfFlux * (source->modelPSF->dparams->data.F32[PM_PAR_I0] / source->modelPSF->params->data.F32[PM_PAR_I0]);
     117    }
     118
     119    // if we have a collection of model fits, check if one of them is a pointer to modelEXT
    124120    if (source->modelFits) {
    125       for (int i = 0; i < source->modelFits->n; i++) {
    126         pmModel *model = source->modelFits->data[i];
    127         status = pmSourcePhotometryModel (&model->mag, model);
    128       }
    129       if (source->modelEXT) {
    130         source->extMag = source->modelEXT->mag;
    131       }
     121        bool foundEXT = false;
     122        for (int i = 0; i < source->modelFits->n; i++) {
     123            pmModel *model = source->modelFits->data[i];
     124            status = pmSourcePhotometryModel (&model->mag, NULL, model);
     125            if (model == source->modelEXT) foundEXT = true;
     126        }
     127        if (foundEXT) {
     128            source->extMag = source->modelEXT->mag;
     129        } else {
     130            status = pmSourcePhotometryModel (&source->extMag, NULL, source->modelEXT);
     131        }
    132132    } else {
    133       if (source->modelEXT) {
    134         status = pmSourcePhotometryModel (&source->extMag, source->modelEXT);
    135       }
     133        if (source->modelEXT) {
     134            status = pmSourcePhotometryModel (&source->extMag, NULL, source->modelEXT);
     135        }
    136136    }
    137137
    … …  
    148148    }
    149149
     150    // measure the contribution of included pixels
     151    if (mode & PM_SOURCE_PHOT_DIFFSTATS) {
     152        pmSourceMeasureDiffStats (source, maskVal);
     153    }
     154
    150155    // measure the aperture magnitude, if (SN > AP_MIN_SN)
    151156    if (!isfinite(SN)) {
    … …  
    160165    }
    161166
    162     // replace source flux
    163     // XXX need to be certain which model and size of mask for prior subtraction
    164     // XXX full model or just analytical?
    165     // XXX use pmSourceAdd instead?
    166     if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) {
    167         pmModelAdd (source->pixels, source->maskObj, model, PM_MODEL_OP_FULL, maskVal);
    168     }
     167    // if we measure aperture magnitudes, the source must not currently be subtracted!
     168    psAssert (!(source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED), "cannot measure ap mags if source is subtracted!");
    169169
    170170    // if we are measuring aperture photometry and applying the growth correction,
    … …  
    201201    if (isfinite (source->apMag) && isPSF && psf) {
    202202        if (psf->growth && (mode & PM_SOURCE_PHOT_GROWTH)) {
    203             source->apMag += pmGrowthCurveCorrect (psf->growth, model->radiusFit);
     203            source->apMag += pmGrowthCurveCorrect (psf->growth, source->apRadius);
    204204        }
    205205        if (mode & PM_SOURCE_PHOT_APCORR) {
     206            // XXX this should be removed -- we no longer fit for the 'sky bias'
    206207            rflux   = pow (10.0, 0.4*source->psfMag);
    207             source->apMag -= PS_SQR(model->radiusFit)*rflux * psf->skyBias + psf->skySat / rflux;
     208            source->apMag -= PS_SQR(source->apRadius)*rflux * psf->skyBias + psf->skySat / rflux;
    208209        }
    209210    }
    … …  
    211212        psFree(flux);
    212213        psFree(mask);
    213     }
    214 
    215     // if source was originally subtracted, re-subtract object, leave local sky
    216     // XXX replace with pmSourceSub...
    217     if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) {
    218         pmModelSub (source->pixels, source->maskObj, model, PM_MODEL_OP_FULL, maskVal);
    219214    }
    220215
    … …  
    230225
    231226// return source model magnitude
    232 bool pmSourcePhotometryModel (float *fitMag, pmModel *model)
    233 {
    234     PS_ASSERT_PTR_NON_NULL(fitMag, false);
    235     if (model == NULL) {
    236         return false;
    237     }
    238 
    239     float fitSum = 0;
    240     *fitMag = NAN;
     227bool pmSourcePhotometryModel (float *fitMag, float *fitFlux, pmModel *model)
     228{
     229    psAssert (fitMag || fitFlux, "at least one of magnitude or flux must be requested (not NULL)");
     230    if (model == NULL) return false;
     231
     232    float mag  = NAN;
     233    float flux = NAN;
    241234
    242235    // measure fitMag
    243     fitSum = model->modelFlux (model->params);
    244     if (fitSum <= 0)
    245         return false;
    246     if (!isfinite(fitSum))
    247         return false;
    248     *fitMag = -2.5*log10(fitSum);
     236    flux = model->modelFlux (model->params);
     237    if (flux > 0) {
     238        mag = -2.5*log10(flux);
     239    }
     240    if (fitMag) {
     241        *fitMag = mag;
     242    }
     243    if (fitFlux) {
     244        *fitFlux = flux;
     245    }
     246
     247    if (flux <= 0) return false;
     248    if (!isfinite(flux)) return false;
    249249
    250250    return (true);
    … …  
    369369
    370370    *pixWeight = validSum / modelSum;
     371    return (true);
     372}
     373
     374# define FLUX_LIMIT 3.0
     375
     376// return source aperture magnitude
     377bool pmSourceMeasureDiffStats (pmSource *source, psImageMaskType maskVal)
     378{
     379    PS_ASSERT_PTR_NON_NULL(source, false);
     380
     381    if (source->diffStats == NULL) {
     382        source->diffStats = pmSourceDiffStatsAlloc();
     383    }
     384
     385    float fGood = 0.0;
     386    float fBad  = 0.0;
     387    int   nGood = 0;
     388    int   nMask = 0;
     389    int   nBad  = 0;
     390   
     391    psImage *flux     = source->pixels;
     392    psImage *variance = source->variance;
     393    psImage *mask     = source->maskObj;
     394
     395    for (int iy = 0; iy < flux->numRows; iy++) {
     396        for (int ix = 0; ix < flux->numCols; ix++) {
     397            if (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & maskVal) {
     398                nMask ++;
     399                continue;
     400            }
     401
     402            float SN = flux->data.F32[iy][ix] / sqrt(variance->data.F32[iy][ix]);
     403
     404            if (SN > +FLUX_LIMIT) {
     405                nGood ++;
     406                fGood += fabs(flux->data.F32[iy][ix]);
     407            }
     408
     409            if (SN < -FLUX_LIMIT) {
     410                nBad ++;
     411                fBad += fabs(flux->data.F32[iy][ix]);
     412            }
     413        }
     414    }
     415
     416    source->diffStats->nGood      = nGood;
     417    source->diffStats->fRatio     = (fGood + fBad         == 0.0) ? NAN : fGood / (fGood + fBad);         
     418    source->diffStats->nRatioBad  = (nGood + nBad         == 0)   ? NAN : nGood / (float) (nGood + nBad);         
     419    source->diffStats->nRatioMask = (nGood + nMask        == 0)   ? NAN : nGood / (float) (nGood + nMask);         
     420    source->diffStats->nRatioAll  = (nGood + nMask + nBad == 0)   ? NAN : nGood / (float) (nGood + nMask + nBad);
     421
    371422    return (true);
    372423}
    … …  
    558609
    559610// determine chisq, etc for linear normalization-only fit
    560 bool pmSourceChisq (pmModel *model, psImage *image, psImage *mask, psImage *variance,
    561                     psImageMaskType maskVal)
     611bool pmSourceChisq (pmModel *model, psImage *image, psImage *mask, psImage *variance, psImageMaskType maskVal, const float covarFactor)
    562612{
    563613    PS_ASSERT_PTR_NON_NULL(model, false);
    … …  
    574624            if (variance->data.F32[j][i] <= 0)
    575625                continue;
    576             dC += PS_SQR (image->data.F32[j][i]) / variance->data.F32[j][i];
     626            dC += PS_SQR (image->data.F32[j][i]) / (covarFactor * variance->data.F32[j][i]);
    577627            Npix ++;
    578628        }
    … …  
    586636
    587637
    588 double pmSourceModelWeight(const pmSource *Mi,
    589                       int term,
    590                       const bool unweighted_sum) // should the cross product be weighted?
     638double pmSourceModelWeight(const pmSource *Mi, int term, const bool unweighted_sum, const float covarFactor)
    591639{
    592640    PS_ASSERT_PTR_NON_NULL(Mi, NAN);
    … …  
    607655                continue;
    608656            if (!unweighted_sum) {
    609                 wt = Wi->data.F32[yi][xi];
     657                wt = covarFactor * Wi->data.F32[yi][xi];
    610658                if (wt == 0)
    611659                    continue;
    … …  
    636684}
    637685
    638 double pmSourceModelDotModel (const pmSource *Mi,
    639                               const pmSource *Mj,
    640                               const bool unweighted_sum) // should the cross product be weighted?
     686double pmSourceModelDotModel (const pmSource *Mi, const pmSource *Mj, const bool unweighted_sum, const float covarFactor)
    641687{
    642688    PS_ASSERT_PTR_NON_NULL(Mi, NAN);
    … …  
    690736                flux += (Pi->data.F32[yi][xi] * Pj->data.F32[yj][xj]);
    691737            } else {
    692                 wt = Wi->data.F32[yi][xi];
     738                wt = covarFactor * Wi->data.F32[yi][xi];
    693739                if (wt > 0) {
    694740                    flux += (Pi->data.F32[yi][xi] * Pj->data.F32[yj][xj]) / wt;
    … …  
    700746}
    701747
    702 double pmSourceDataDotModel (const pmSource *Mi,
    703                              const pmSource *Mj,
    704                              const bool unweighted_sum) // should the cross product be weighted?
     748double pmSourceDataDotModel (const pmSource *Mi, const pmSource *Mj, const bool unweighted_sum, const float covarFactor)
    705749{
    706750    PS_ASSERT_PTR_NON_NULL(Mi, NAN);
    … …  
    754798                flux += (Pi->data.F32[yi][xi] * Pj->data.F32[yj][xj]);
    755799            } else {
    756                 wt = Wi->data.F32[yi][xi];
     800                wt = covarFactor * Wi->data.F32[yi][xi];
    757801                if (wt > 0) {
    758802                    flux += (Pi->data.F32[yi][xi] * Pj->data.F32[yj][xj]) / wt;
    … …  
    763807    return flux;
    764808}
    765 
    766 // XXX this is test code to verify the shift is doing the right thing (seems to be)
    767 # if (0)
    768 // measure centroid of unshifted gaussian (should be 16.0,16.0)
    769         {
    770           psImage *image = source->pixels;
    771           float xo = 0.0;
    772           float yo = 0.0;
    773           float xo2 = 0.0;
    774           float yo2 = 0.0;
    775           float no = 0.0;
    776           for (int j = 0; j < image->numRows; j++)
    777           {
    778             for (int i = 0; i < image->numCols; i++) {
    779               xo += i*image->data.F32[j][i];
    780               yo += j*image->data.F32[j][i];
    781               xo2 += i*i*image->data.F32[j][i];
    782               yo2 += j*j*image->data.F32[j][i];
    783               no += image->data.F32[j][i];
    784             }
    785           }
    786           xo /= no;
    787           yo /= no;
    788           xo2 = sqrt (xo2/no - xo*xo);
    789           yo2 = sqrt (yo2/no - yo*yo);
    790           fprintf (stderr, "pre-shift centroid: %f,%f, sigma: %f,%f: flux: %f\n", xo, yo, xo2, yo2, no);
    791         }
    792 
    793 // measure centroid of unshifted gaussian (should be 16.0,16.0)
    794         {
    795           psImage *image = flux;
    796           float xo = 0.0;
    797           float yo = 0.0;
    798           float xo2 = 0.0;
    799           float yo2 = 0.0;
    800           float no = 0.0;
    801           for (int j = 0; j < image->numRows; j++)
    802           {
    803             for (int i = 0; i < image->numCols; i++) {
    804               xo += i*image->data.F32[j][i];
    805               yo += j*image->data.F32[j][i];
    806               xo2 += i*i*image->data.F32[j][i];
    807               yo2 += j*j*image->data.F32[j][i];
    808               no += image->data.F32[j][i];
    809             }
    810           }
    811           xo /= no;
    812           yo /= no;
    813           xo2 = sqrt (xo2/no - xo*xo);
    814           yo2 = sqrt (yo2/no - yo*yo);
    815           fprintf (stderr, "pre-shift centroid: %f,%f, sigma: %f,%f: flux: %f\n", xo, yo, xo2, yo2, no);
    816         }
    817 # endif
Note: See TracChangeset for help on using the changeset viewer.