IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 28426


Ignore:
Timestamp:
Jun 22, 2010, 3:00:02 PM (16 years ago)
Author:
Paul Price
Message:

Linear fit was ignoring pixels with any mask bit set (including the 'advisory' mask bits: SUSPECT/BURNTOOL,CONV.POOR), rather than using the set of bad mask bits. This meant that linear photometry fits resulted in a peak flux of 0, and therefore a bad PSF magnitude. This also resulted in bright sources (meeting the first S/N cut) being detected twice because they were not subtracted properly.

Location:
trunk
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • trunk/psModules/src/objects/pmSourcePhotometry.c

    r27531 r28426  
    107107        // the source peak pixel is guaranteed to be on the image, and only minimally different from the source center
    108108        double fluxScale = pmTrend2DEval (psf->FluxScale, (float)source->peak->x, (float)source->peak->y);
    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);
     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);
    114114    } else {
    115115        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]);
     116        source->psfFluxErr = source->psfFlux * (source->modelPSF->dparams->data.F32[PM_PAR_I0] / source->modelPSF->params->data.F32[PM_PAR_I0]);
    117117    }
    118118
    … …  
    120120    if (source->modelFits) {
    121121        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         }
     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, NULL, source->modelEXT);
    135         }
     133        if (source->modelEXT) {
     134            status = pmSourcePhotometryModel (&source->extMag, NULL, source->modelEXT);
     135        }
    136136    }
    137137
    … …  
    204204        }
    205205        if (mode & PM_SOURCE_PHOT_APCORR) {
    206             // XXX this should be removed -- we no longer fit for the 'sky bias'
     206            // XXX this should be removed -- we no longer fit for the 'sky bias'
    207207            rflux   = pow (10.0, 0.4*source->psfMag);
    208208            source->apMag -= PS_SQR(source->apRadius)*rflux * psf->skyBias + psf->skySat / rflux;
    … …  
    236236    flux = model->modelFlux (model->params);
    237237    if (flux > 0) {
    238         mag = -2.5*log10(flux);
     238        mag = -2.5*log10(flux);
    239239    }
    240240    if (fitMag) {
    241         *fitMag = mag;
     241        *fitMag = mag;
    242242    }
    243243    if (fitFlux) {
    244         *fitFlux = flux;
     244        *fitFlux = flux;
    245245    }
    246246
    … …  
    380380
    381381    if (source->diffStats == NULL) {
    382         source->diffStats = pmSourceDiffStatsAlloc();
     382        source->diffStats = pmSourceDiffStatsAlloc();
    383383    }
    384384
    … …  
    388388    int   nMask = 0;
    389389    int   nBad  = 0;
    390    
     390
    391391    psImage *flux     = source->pixels;
    392392    psImage *variance = source->variance;
    … …  
    394394
    395395    for (int iy = 0; iy < flux->numRows; iy++) {
    396         for (int ix = 0; ix < flux->numCols; ix++) {
     396        for (int ix = 0; ix < flux->numCols; ix++) {
    397397            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         }
     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        }
    414414    }
    415415
    416416    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);         
     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);
    420420    source->diffStats->nRatioAll  = (nGood + nMask + nBad == 0)   ? NAN : nGood / (float) (nGood + nMask + nBad);
    421421
    … …  
    628628        }
    629629    }
    630     model->nPix = Npix; 
     630    model->nPix = Npix;
    631631    model->nDOF = Npix - 1;
    632632    model->chisq = dC;
    … …  
    636636
    637637
    638 double pmSourceModelWeight(const pmSource *Mi, int term, const bool unweighted_sum, const float covarFactor)
     638double pmSourceModelWeight(const pmSource *Mi, int term, const bool unweighted_sum, const float covarFactor, psImageMaskType maskVal)
    639639{
    640640    PS_ASSERT_PTR_NON_NULL(Mi, NAN);
    … …  
    652652    for (int yi = 0; yi < Pi->numRows; yi++) {
    653653        for (int xi = 0; xi < Pi->numCols; xi++) {
    654             if (Ti->data.PS_TYPE_IMAGE_MASK_DATA[yi][xi])
     654            if (Ti->data.PS_TYPE_IMAGE_MASK_DATA[yi][xi] & maskVal)
    655655                continue;
    656656            if (!unweighted_sum) {
    … …  
    684684}
    685685
    686 double pmSourceModelDotModel (const pmSource *Mi, const pmSource *Mj, const bool unweighted_sum, const float covarFactor)
     686double pmSourceModelDotModel (const pmSource *Mi, const pmSource *Mj, const bool unweighted_sum, const float covarFactor, psImageMaskType maskVal)
    687687{
    688688    PS_ASSERT_PTR_NON_NULL(Mi, NAN);
    … …  
    727727    for (yi = yIs, yj = yJs; yi < yIe; yi++, yj++) {
    728728        for (xi = xIs, xj = xJs; xi < xIe; xi++, xj++) {
    729             if (Ti->data.PS_TYPE_IMAGE_MASK_DATA[yi][xi])
    730                 continue;
    731             if (Tj->data.PS_TYPE_IMAGE_MASK_DATA[yj][xj])
     729            if (Ti->data.PS_TYPE_IMAGE_MASK_DATA[yi][xi] & maskVal)
     730                continue;
     731            if (Tj->data.PS_TYPE_IMAGE_MASK_DATA[yj][xj] & maskVal)
    732732                continue;
    733733
    … …  
    746746}
    747747
    748 double pmSourceDataDotModel (const pmSource *Mi, const pmSource *Mj, const bool unweighted_sum, const float covarFactor)
     748double pmSourceDataDotModel (const pmSource *Mi, const pmSource *Mj, const bool unweighted_sum, const float covarFactor, psImageMaskType maskVal)
    749749{
    750750    PS_ASSERT_PTR_NON_NULL(Mi, NAN);
    … …  
    789789    for (yi = yIs, yj = yJs; yi < yIe; yi++, yj++) {
    790790        for (xi = xIs, xj = xJs; xi < xIe; xi++, xj++) {
    791             if (Ti->data.PS_TYPE_IMAGE_MASK_DATA[yi][xi])
    792                 continue;
    793             if (Tj->data.PS_TYPE_IMAGE_MASK_DATA[yj][xj])
     791            if (Ti->data.PS_TYPE_IMAGE_MASK_DATA[yi][xi] & maskVal)
     792                continue;
     793            if (Tj->data.PS_TYPE_IMAGE_MASK_DATA[yj][xj] & maskVal)
    794794                continue;
    795795

  • trunk/psModules/src/objects/pmSourcePhotometry.h

    r27531 r28426  
    4848    psImage *image,                     ///< image pixels to be used
    4949    psImage *mask,                      ///< mask of pixels to ignore
    50     psImageMaskType maskVal             ///< Value to mask
     50    psImageMaskType maskVal             ///< Value to mask
    5151);
    5252
    … …  
    5858bool pmSourceMeasureDiffStats (pmSource *source, psImageMaskType maskVal);
    5959
    60 double pmSourceDataDotModel (const pmSource *Mi, const pmSource *Mj, const bool unweighted_sum, const float covarFactor);
    61 double pmSourceModelDotModel (const pmSource *Mi, const pmSource *Mj, const bool unweighted_sum, const float covarFactor);
    62 double pmSourceModelWeight(const pmSource *Mi, int term, const bool unweighted_sum, const float covarFactor);
     60double pmSourceDataDotModel (const pmSource *Mi, const pmSource *Mj, const bool unweighted_sum, const float covarFactor, psImageMaskType maskVal);
     61double pmSourceModelDotModel (const pmSource *Mi, const pmSource *Mj, const bool unweighted_sum, const float covarFactor, psImageMaskType maskVal);
     62double pmSourceModelWeight(const pmSource *Mi, int term, const bool unweighted_sum, const float covarFactor, psImageMaskType maskVal);
    6363
    6464// retire these:

  • trunk/psphot/src/psphotFitSourcesLinear.c

    r28399 r28426  
    171171
    172172        // diagonal elements of the sparse matrix (auto-cross-product)
    173         f = pmSourceModelDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor);
     173        f = pmSourceModelDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor, maskVal);
    174174        psSparseMatrixElement (sparse, i, i, f);
    175175
    176176        // the formal error depends on the weighting scheme
    177177        if (CONSTANT_PHOTOMETRIC_WEIGHTS) {
    178             float var = pmSourceModelDotModel (SRCi, SRCi, false, covarFactor);
     178            float var = pmSourceModelDotModel (SRCi, SRCi, false, covarFactor, maskVal);
    179179            errors->data.F32[i] = 1.0 / sqrt(var);
    180180        } else {
    … …  
    184184
    185185        // find the image x model value
    186         f = pmSourceDataDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor);
     186        f = pmSourceDataDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor, maskVal);
    187187        psSparseVectorElement (sparse, i, f);
    188188
    … …  
    190190        switch (SKY_FIT_ORDER) {
    191191          case 1:
    192             f = pmSourceModelWeight (SRCi, 1, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor);
     192            f = pmSourceModelWeight (SRCi, 1, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor, maskVal);
    193193            psSparseBorderElementB (border, i, 1, f);
    194             f = pmSourceModelWeight (SRCi, 2, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor);
     194            f = pmSourceModelWeight (SRCi, 2, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor, maskVal);
    195195            psSparseBorderElementB (border, i, 2, f);
    196196
    197197          case 0:
    198             f = pmSourceModelWeight (SRCi, 0, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor);
     198            f = pmSourceModelWeight (SRCi, 0, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor, maskVal);
    199199            psSparseBorderElementB (border, i, 0, f);
    200200            break;
    … …  
    216216
    217217            // got an overlap; calculate cross-product and add to output array
    218             f = pmSourceModelDotModel (SRCi, SRCj, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor);
     218            f = pmSourceModelDotModel (SRCi, SRCj, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor, maskVal);
    219219            psSparseMatrixElement (sparse, j, i, f);
    220220        }

  • trunk/psphot/src/psphotFitSourcesLinearStack.c

    r28013 r28426  
    4343    for (int i = 0; i < objects->n; i++) {
    4444        pmPhotObj *object = objects->data[i];
    45         if (!object) continue;
    46         if (!object->sources) continue;
     45        if (!object) continue;
     46        if (!object->sources) continue;
    4747
    48         // XXX check an element of the group to see if we should use it
    49         // if (!object->flags & PM_PHOT_OBJ_BAD) continue;
     48        // XXX check an element of the group to see if we should use it
     49        // if (!object->flags & PM_PHOT_OBJ_BAD) continue;
    5050
    51         for (int j = 0; j < object->sources->n; j++) {
    52           pmSource *source = object->sources->data[j];
    53           if (!source) continue;
     51        for (int j = 0; j < object->sources->n; j++) {
     52          pmSource *source = object->sources->data[j];
     53          if (!source) continue;
    5454
    55           // turn this bit off and turn it on again if we keep this source
    56           source->mode &= ~PM_SOURCE_MODE_LINEAR_FIT;
     55          // turn this bit off and turn it on again if we keep this source
     56          source->mode &= ~PM_SOURCE_MODE_LINEAR_FIT;
    5757
    58           // generate model for sources without, or skip if we can't
    59           if (!source->modelFlux) {
     58          // generate model for sources without, or skip if we can't
     59          if (!source->modelFlux) {
    6060            if (!pmSourceCacheModel (source, maskVal)) continue;
    61           }
     61          }
    6262
    63           source->mode |= PM_SOURCE_MODE_LINEAR_FIT;
    64           psArrayAdd (fitSources, 100, source);
    65         }
     63          source->mode |= PM_SOURCE_MODE_LINEAR_FIT;
     64          psArrayAdd (fitSources, 100, source);
     65        }
    6666    }
    6767    psLogMsg ("psphot.ensemble", PS_LOG_MINUTIA, "built fitSources: %f sec (%ld objects)\n", psTimerMark ("psphot.linear"), objects->n);
    … …  
    8585
    8686        // diagonal elements of the sparse matrix (auto-cross-product)
    87         f = pmSourceModelDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS, COVAR_FACTOR);
     87        f = pmSourceModelDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS, COVAR_FACTOR, maskVal);
    8888        psSparseMatrixElement (sparse, i, i, f);
    8989
    9090        // the formal error depends on the weighting scheme
    9191        if (CONSTANT_PHOTOMETRIC_WEIGHTS) {
    92             float var = pmSourceModelDotModel (SRCi, SRCi, false, COVAR_FACTOR);
     92            float var = pmSourceModelDotModel (SRCi, SRCi, false, COVAR_FACTOR, maskVal);
    9393            errors->data.F32[i] = 1.0 / sqrt(var);
    9494        } else {
    … …  
    9797
    9898        // find the image x model value
    99         f = pmSourceDataDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS, COVAR_FACTOR);
     99        f = pmSourceDataDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS, COVAR_FACTOR, maskVal);
    100100        psSparseVectorElement (sparse, i, f);
    101101
    … …  
    104104            pmSource *SRCj = fitSources->data[j];
    105105
    106             // we only need to generate dot terms for source on the same image
    107             if (SRCj->imageID != SRCi->imageID) { continue; }
     106            // we only need to generate dot terms for source on the same image
     107            if (SRCj->imageID != SRCi->imageID) { continue; }
    108108
    109109            // skip over disjoint source images, break after last possible overlap
    … …  
    114114
    115115            // got an overlap; calculate cross-product and add to output array
    116             f = pmSourceModelDotModel (SRCi, SRCj, CONSTANT_PHOTOMETRIC_WEIGHTS, COVAR_FACTOR);
     116            f = pmSourceModelDotModel (SRCi, SRCj, CONSTANT_PHOTOMETRIC_WEIGHTS, COVAR_FACTOR, maskVal);
    117117            psSparseMatrixElement (sparse, j, i, f);
    118118        }
Note: See TracChangeset for help on using the changeset viewer.