IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 29546


Ignore:
Timestamp:
Oct 25, 2010, 2:54:22 PM (16 years ago)
Author:
eugene
Message:

merge changes from eam_branches/ipp-20100823

Location:
trunk/psModules/src/objects
Files:
12 edited

Legend:

Unmodified
Added
Removed

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

    r29004 r29546  
    197197        // mask the given aperture and measure the apMag
    198198        psImageKeepCircle (mask, xc, yc, radius, "OR", markVal);
    199         if (!pmSourcePhotometryAper (&apMag, model, pixels, mask, maskVal)) {
     199        if (!pmSourcePhotometryAper (&apMag, NULL, NULL, model, pixels, NULL, mask, maskVal)) {
    200200            psFree (growth);
    201201            psFree (view);

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

    r29004 r29546  
    104104    assert (source->psfImage); // XXX build if needed?
    105105
    106     int x0 = source->peak->xf - source->psfImage->col0;
    107     int y0 = source->peak->yf - source->psfImage->row0;
     106    int x0 = source->peak->x - source->psfImage->col0;
     107    int y0 = source->peak->y - source->psfImage->row0;
    108108
    109109    // need to decide on the size: dynamically? statically?

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

    r29004 r29546  
    132132
    133133    // default values are NAN
    134     source->psfMag     = NAN;
    135     source->psfFlux    = NAN;
    136     source->psfFluxErr = NAN;
    137     source->extMag = NAN;
    138     source->errMag = NAN;
    139     source->apMag  = NAN;
    140     source->sky    = NAN;
    141     source->skyErr = NAN;
    142     source->pixWeightNotBad = NAN;
     134    source->psfMag           = NAN;
     135    source->psfFlux          = NAN;
     136    source->psfFluxErr       = NAN;
     137    source->extMag           = NAN;   
     138    source->errMag           = NAN;
     139    source->apMag            = NAN;
     140    source->apMagRaw         = NAN;
     141    source->apRadius         = NAN;
     142    source->apFlux           = NAN;
     143    source->apFluxErr        = NAN;
     144    source->sky              = NAN;
     145    source->skyErr           = NAN;   
     146    source->pixWeightNotBad  = NAN;
    143147    source->pixWeightNotPoor = NAN;
    144148
    145     source->psfChisq = NAN;
    146     source->crNsigma = NAN;
    147     source->extNsigma = NAN;
     149    source->psfChisq         = NAN;
     150    source->crNsigma         = NAN;
     151    source->extNsigma        = NAN;
    148152
    149153    psTrace("psModules.objects", 10, "---- end ----\n");
    … …  
    644648        // The following determinations require the use of moments
    645649        if (!(source->mode & noMoments)) {
    646             // likely defect (too small to be stellar) (push out to 3 sigma)
    647             // low S/N objects which are small are probably stellar
    648             // XXX these limits are quite arbitrary
    649             if (sigX < 0.05 || sigY < 0.05) {
     650            // likely defect (bright, but too small to be stellar)
     651            // XXX eliminate the classification?
     652            if ((source->moments->SN > 10) && (sigX < 0.05 || sigY < 0.05)) {
    650653                source->type = PM_SOURCE_TYPE_DEFECT;
    651654                source->mode |= PM_SOURCE_MODE_DEFECT;

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

    r29004 r29546  
    7676    float apMagRaw;                     ///< raw mag in given aperture
    7777    float apRadius;                     ///< radius for aperture magnitude
     78    float apFlux;                       ///< apFlux corresponding to psfMag or extMag (depending on type)
     79    float apFluxErr;                    ///< apFluxErr corresponding to psfMag or extMag (depending on type)
    7880
    7981    float pixWeightNotBad;              ///< PSF-weighted coverage of unmasked (not BAD) pixels
    … …  
    234236);
    235237
     238/** pmSourceMoments()
     239 *
     240 * Measure 1st moments for the given source, using the peak coordinates as the initial
     241 * source location. The resulting moment values are applied to the source.moments
     242 * entry. The moments are measured within the given circular radius of the source.peak
     243 * coordinates.  The return value indicates the success (TRUE) of the operation.
     244 *
     245 */
     246bool pmSourceMomentsGetCentroid(
     247  pmSource *source,
     248  psF32 radius,
     249  psF32 sigma,
     250  psF32 minSN,
     251  psImageMaskType maskVal,
     252  float xGuess, float yGuess);
     253
    236254pmModel *pmSourceGetModel (bool *isPSF, const pmSource *source);
    237255

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

    r29004 r29546  
    10011001                sources = pmSourcesRead_CMF_PS1_V2 (file->fits, hdu->header);
    10021002            }
     1003            if (!strcmp (exttype, "PS1_V3")) {
     1004                sources = pmSourcesRead_CMF_PS1_V3 (file->fits, hdu->header);
     1005            }
    10031006            if (!strcmp (exttype, "PS1_SV1")) {
    10041007                sources = pmSourcesRead_CMF_PS1_SV1 (file->fits, hdu->header);

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

    r29004 r29546  
    7575psArray *pmSourcesRead_CMF_PS1_V1 (psFits *fits, psMetadata *header);
    7676psArray *pmSourcesRead_CMF_PS1_V2 (psFits *fits, psMetadata *header);
     77psArray *pmSourcesRead_CMF_PS1_V3 (psFits *fits, psMetadata *header);
    7778psArray *pmSourcesRead_CMF_PS1_SV1 (psFits *fits, psMetadata *header);
    7879psArray *pmSourcesRead_CMF_PS1_DV1 (psFits *fits, psMetadata *header);

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

    r29004 r29546  
    185185        psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_FLUX",    PS_DATA_F32, "PSF fit instrumental magnitude",             source->psfFlux);
    186186        psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_FLUX_SIG",PS_DATA_F32, "Sigma of PSF instrumental magnitude",        source->psfFluxErr);
     187
    187188        psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG",           PS_DATA_F32, "magnitude in standard aperture",             source->apMag);
     189        psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RAW",       PS_DATA_F32, "magnitude in real aperture",                 source->apMagRaw);
    188190        psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RADIUS",    PS_DATA_F32, "radius used for aperture mags",              apRadius);
     191        psMetadataAdd (row, PS_LIST_TAIL, "AP_FLUX",          PS_DATA_F32, "instrumental flux in standard aperture",     source->apFlux);
     192        psMetadataAdd (row, PS_LIST_TAIL, "AP_FLUX_SIG",      PS_DATA_F32, "aperture flux error",                        source->apFluxErr);
     193
    189194        psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude",           peakMag);
    190195        psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG",      PS_DATA_F32, "PSF Magnitude using supplied calibration",   calMag);
    … …  
    214219        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XY",       PS_DATA_F32, "second moments (X*Y)",                      mxy);
    215220        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_YY",       PS_DATA_F32, "second moments (Y*Y)",                      myy);
     221
     222        float Mrf  = source->moments ? source->moments->Mrf : NAN;
     223        float Mrh  = source->moments ? source->moments->Mrh : NAN;
     224        float Krf  = source->moments ? source->moments->KronFlux : NAN;
     225        float dKrf = source->moments ? source->moments->KronFluxErr : NAN;
     226
     227        float Kinner = source->moments ? source->moments->KronFinner : NAN;
     228        float Kouter = source->moments ? source->moments->KronFouter : NAN;
     229
     230        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_R1",       PS_DATA_F32, "first radial moment",                       Mrf);
     231        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_RH",       PS_DATA_F32, "half radial moment",                        Mrh);
     232        psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX",        PS_DATA_F32, "Kron Flux (in 2.5 R1)",                     Krf);
     233        psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_ERR",    PS_DATA_F32, "Kron Flux Error",                          dKrf);
     234
     235        psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_INNER",  PS_DATA_F32, "Kron Flux (in 1.0 R1)",                     Kinner);
     236        psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_OUTER",  PS_DATA_F32, "Kron Flux (in 4.0 R1)",                     Kouter);
    216237
    217238        psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NPOS",        PS_DATA_S32, "nPos (n pix > 3 sigma)",                    diffStats.nGood);

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

    r29012 r29546  
    6666    psF32 xPos, yPos;
    6767    psF32 xErr, yErr;
    68     psF32 errMag, chisq, apRadius;
     68    psF32 chisq, apRadius;
    6969    psS32 nPix, nDOF;
    7070
    … …  
    8484    }
    8585    float zeroptErr = psMetadataLookupF32 (&status2, imageHeader, "ZPT_ERR");
     86
     87    // we need a measure of the image quality (FWHM) for this image, in order to get the positional errors
     88    float fwhmMajor = psMetadataLookupF32(&status1, readout->analysis, "FWHM_MAJ");
     89    if (!status1) {
     90        fwhmMajor = psMetadataLookupF32(&status1, readout->analysis, "IQ_FW1");
     91    }
     92    float fwhmMinor = psMetadataLookupF32(&status1, readout->analysis, "FWHM_MIN");
     93    if (!status1) {
     94        fwhmMinor = psMetadataLookupF32(&status1, readout->analysis, "IQ_FW2");
     95    }
    8696
    8797    // if the sequence is defined, write these in seq order; otherwise
    … …  
    124134                yErr = dPAR[PM_PAR_YPOS];
    125135            } else {
    126                 // in linear-fit mode, there is no error on the centroid
    127                 xErr = source->peak->dx;
    128                 yErr = source->peak->dy;
     136                xErr = fwhmMajor * source->errMag / 2.35;
     137                yErr = fwhmMinor * source->errMag / 2.35;
    129138            }
    130139            if (isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX])) {
    … …  
    139148            nPix = model->nPix;
    140149            apRadius = source->apRadius;
    141             errMag = model->dparams->data.F32[PM_PAR_I0] / model->params->data.F32[PM_PAR_I0];
    142150        } else {
    143             xPos = source->peak->xf;
    144             yPos = source->peak->yf;
    145             xErr = source->peak->dx;
    146             yErr = source->peak->dy;
     151            bool useMoments = true;
     152            useMoments = (useMoments && source->moments);          // can't if there are no moments
     153            useMoments = (useMoments && source->moments->nPixels); // can't if the moments were not measured
     154            useMoments = (useMoments && !(source->mode && PM_SOURCE_MODE_MOMENTS_FAILURE)); // can't if the moments failed...
     155
     156            if (source->moments) {
     157                xPos = source->moments->Mx;
     158                yPos = source->moments->My;
     159                xErr = fwhmMajor * source->errMag / 2.35;
     160                yErr = fwhmMinor * source->errMag / 2.35;
     161            } else {
     162                xPos = source->peak->xf;
     163                yPos = source->peak->yf;
     164                xErr = source->peak->dx;
     165                yErr = source->peak->dy;
     166            }
    147167            axes.major = NAN;
    148168            axes.minor = NAN;
    … …  
    152172            nPix = 0;
    153173            apRadius = NAN;
    154             errMag = NAN;
    155174        }
    156175
    … …  
    168187        psMetadataAdd (row, PS_LIST_TAIL, "X_PSF",            PS_DATA_F32, "PSF x coordinate",                           xPos);
    169188        psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF",            PS_DATA_F32, "PSF y coordinate",                           yPos);
    170         psMetadataAdd (row, PS_LIST_TAIL, "X_PSF_SIG",        PS_DATA_F32, "Sigma in PSF x coordinate",                  xErr); // XXX this is only measured for non-linear fits
    171         psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF_SIG",        PS_DATA_F32, "Sigma in PSF y coordinate",                  yErr); // XXX this is only measured for non-linear fits
     189        psMetadataAdd (row, PS_LIST_TAIL, "X_PSF_SIG",        PS_DATA_F32, "Sigma in PSF x coordinate",                  xErr);
     190        psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF_SIG",        PS_DATA_F32, "Sigma in PSF y coordinate",                  yErr);
    172191        psMetadataAdd (row, PS_LIST_TAIL, "POSANGLE",         PS_DATA_F32, "position angle at source (degrees)",         posAngle*PS_DEG_RAD);
    173192        psMetadataAdd (row, PS_LIST_TAIL, "PLTSCALE",         PS_DATA_F32, "plate scale at source (arcsec/pixel)",       pltScale*PS_DEG_RAD*3600.0);
    174193        psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG",     PS_DATA_F32, "PSF fit instrumental magnitude",             source->psfMag);
    175         psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG_SIG", PS_DATA_F32, "Sigma of PSF instrumental magnitude",        errMag);
     194        psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG_SIG", PS_DATA_F32, "Sigma of PSF instrumental magnitude",        source->errMag);
    176195        psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_FLUX",    PS_DATA_F32, "PSF fit instrumental flux (counts)",         source->psfFlux);
    177196        psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_FLUX_SIG",PS_DATA_F32, "Sigma of PSF instrumental flux",             source->psfFluxErr);
    … …  
    179198        psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RAW",       PS_DATA_F32, "magnitude in reported aperture",             source->apMagRaw);
    180199        psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RADIUS",    PS_DATA_F32, "radius used for aperture mags",              apRadius);
    181         psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude",           peakMag);
     200
    182201        psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG",      PS_DATA_F32, "PSF Magnitude using supplied calibration",   calMag);
    183202        psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG_SIG",  PS_DATA_F32, "measured scatter of zero point calibration", zeroptErr);
     203       
     204        // NOTE: RA & DEC (both double) need to be on an 8-byte boundary...
    184205        psMetadataAdd (row, PS_LIST_TAIL, "RA_PSF",           PS_DATA_F64, "PSF RA coordinate (degrees)",                ptSky.r*PS_DEG_RAD);
    185206        psMetadataAdd (row, PS_LIST_TAIL, "DEC_PSF",          PS_DATA_F64, "PSF DEC coordinate (degrees)",               ptSky.d*PS_DEG_RAD);
     207
     208        psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude",           peakMag);
    186209        psMetadataAdd (row, PS_LIST_TAIL, "SKY",              PS_DATA_F32, "Sky level",                                  source->sky);
    187210        psMetadataAdd (row, PS_LIST_TAIL, "SKY_SIGMA",        PS_DATA_F32, "Sigma of sky level",                         source->skyErr);
    … …  
    204227        float Myy = source->moments ? source->moments->Myy : NAN;
    205228
     229        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XX",       PS_DATA_F32, "second moments (X^2)",                      Mxx);
     230        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XY",       PS_DATA_F32, "second moments (X*Y)",                      Mxy);
     231        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_YY",       PS_DATA_F32, "second moments (Y*Y)",                      Myy);
     232
     233        float M_c3 = source->moments ? 1.0*source->moments->Mxxx - 3.0*source->moments->Mxyy : NAN;
     234        float M_s3 = source->moments ? 3.0*source->moments->Mxxy - 1.0*source->moments->Myyy : NAN;
     235        float M_c4 = source->moments ? 1.0*source->moments->Mxxxx - 6.0*source->moments->Mxxyy + 1.0*source->moments->Myyyy : NAN;
     236        float M_s4 = source->moments ? 4.0*source->moments->Mxxxy - 4.0*source->moments->Mxyyy : NAN;
     237
     238        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M3C",      PS_DATA_F32, "third momemt cos theta",                    M_c3);
     239        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M3S",      PS_DATA_F32, "third momemt sin theta",                    M_s3);
     240        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M4C",      PS_DATA_F32, "fourth momemt cos theta",                   M_c4);
     241        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M4S",      PS_DATA_F32, "fourth momemt sin theta",                   M_s4);
     242
    206243        float Mrf  = source->moments ? source->moments->Mrf : NAN;
    207244        float Mrh  = source->moments ? source->moments->Mrh : NAN;
    … …  
    212249        float Kouter = source->moments ? source->moments->KronFouter : NAN;
    213250
    214         float M_c3 = source->moments ? 1.0*source->moments->Mxxx - 3.0*source->moments->Mxyy : NAN;
    215         float M_s3 = source->moments ? 3.0*source->moments->Mxxy - 1.0*source->moments->Myyy : NAN;
    216         float M_c4 = source->moments ? 1.0*source->moments->Mxxxx - 6.0*source->moments->Mxxyy + 1.0*source->moments->Myyyy : NAN;
    217         float M_s4 = source->moments ? 4.0*source->moments->Mxxxy - 4.0*source->moments->Mxyyy : NAN;
    218 
    219         psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XX",       PS_DATA_F32, "second moments (X^2)",                      Mxx);
    220         psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XY",       PS_DATA_F32, "second moments (X*Y)",                      Mxy);
    221         psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_YY",       PS_DATA_F32, "second moments (Y*Y)",                      Myy);
    222 
    223         psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M3C",      PS_DATA_F32, "third momemt cos theta",                    M_c3);
    224         psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M3S",      PS_DATA_F32, "third momemt sin theta",                    M_s3);
    225         psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M4C",      PS_DATA_F32, "fourth momemt cos theta",                   M_c4);
    226         psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M4S",      PS_DATA_F32, "fourth momemt sin theta",                   M_s4);
    227 
    228251        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_R1",       PS_DATA_F32, "first radial moment",                       Mrf);
    229252        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_RH",       PS_DATA_F32, "half radial moment",                        Mrh);
    … …  
    236259        psMetadataAdd (row, PS_LIST_TAIL, "FLAGS",            PS_DATA_U32, "psphot analysis flags",                     source->mode);
    237260        psMetadataAdd (row, PS_LIST_TAIL, "FLAGS2",           PS_DATA_U32, "psphot analysis flags",                     source->mode2);
     261        psMetadataAdd (row, PS_LIST_TAIL, "PADDING2",         PS_DATA_S32, "more padding", 0);
    238262
    239263        // XXX not sure how to get this : need to load Nimages with weight?

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

    r29004 r29546  
    4545    PM_SOURCE_MODE2_DIFF_WITH_SINGLE = 0x00000001, ///< diff source matched to a single positive detection
    4646    PM_SOURCE_MODE2_DIFF_WITH_DOUBLE = 0x00000002, ///< diff source matched to positive detections in both images
     47    PM_SOURCE_MODE2_MATCHED          = 0x00000004, ///< diff source matched to positive detections in both images
    4748} pmSourceMode2;
    4849

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

    r29044 r29546  
    6464void pmSourceMomentsSetVerbose(bool state){ beVerbose = state; }
    6565
     66// if mode & EXTERNAL or mode2 & MATCHED, do not re-calculate the centroid (use peak as centroid)
     67
    6668bool pmSourceMoments(pmSource *source, psF32 radius, psF32 sigma, psF32 minSN, psImageMaskType maskVal)
    6769{
    … …  
    7173    PS_ASSERT_FLOAT_LARGER_THAN(radius, 0.0, false);
    7274
    73     // use sky from moments if defined, 0.0 otherwise
    74 
    75     // XXX this value comes from the sky model at the source center, and tends to over-estimate
    76     // the sky in the vicinity of bright sources.  we are better off assuming the model worked
    77     // well:
     75    // this function assumes the sky has been well-subtracted for the image
    7876    psF32 sky = 0.0;
     77
    7978    if (source->moments == NULL) {
    8079      source->moments = pmMomentsAlloc();
    8180    }
    82     // XXX if (source->moments == NULL) {
    83     // XXX     source->moments = pmMomentsAlloc();
    84     // XXX } else {
    85     // XXX     sky = source->moments->Sky;
    86     // XXX }
    87 
    88     // First Pass: calculate the first moments (these are subtracted from the coordinates below)
    89     // Sum = SUM (z - sky)
    90     // X1  = SUM (x - xc)*(z - sky)
    91     // .. etc
    92 
    93     psF32 peakPixel = -PS_MAX_F32;
    94     psS32 numPixels = 0;
     81
    9582    psF32 Sum = 0.0;
    9683    psF32 Var = 0.0;
    97     psF32 X1 = 0.0;
    98     psF32 Y1 = 0.0;
    9984    psF32 R2 = PS_SQR(radius);
    10085    psF32 minSN2 = PS_SQR(minSN);
    … …  
    10994    // (int) so they can be used in the image index below.
    11095
    111     int xOff  = source->peak->x;
    112     int yOff  = source->peak->y;
    113     int xPeak = source->peak->x - source->pixels->col0; // coord of peak in subimage
    114     int yPeak = source->peak->y - source->pixels->row0; // coord of peak in subimage
    115 
    116     // we are guaranteed to have a valid pixel and variance at this location (right? right?)
    117     // psF32 weightNorm = source->pixels->data.F32[yPeak][xPeak] / sqrt (source->variance->data.F32[yPeak][xPeak]);
    118     // psAssert (isfinite(source->pixels->data.F32[yPeak][xPeak]), "peak must be on valid pixel");
    119     // psAssert (isfinite(source->variance->data.F32[yPeak][xPeak]), "peak must be on valid pixel");
    120     // psAssert (source->variance->data.F32[yPeak][xPeak] > 0, "peak must be on valid pixel");
    121 
    122     for (psS32 row = 0; row < source->pixels->numRows ; row++) {
    123 
    124         psF32 yDiff = row - yPeak;
    125         if (fabs(yDiff) > radius) continue;
    126 
    127         psF32 *vPix = source->pixels->data.F32[row];
    128         psF32 *vWgt = source->variance->data.F32[row];
    129         psImageMaskType *vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[row];
    130 
    131         for (psS32 col = 0; col < source->pixels->numCols ; col++, vPix++, vWgt++) {
    132             if (vMsk) {
    133                 if (*vMsk & maskVal) {
    134                     vMsk++;
    135                     continue;
    136                 }
    137                 vMsk++;
    138             }
    139             if (isnan(*vPix)) continue;
    140 
    141             psF32 xDiff = col - xPeak;
    142             if (fabs(xDiff) > radius) continue;
    143 
    144             // radius is just a function of (xDiff, yDiff)
    145             if (!VALID_RADIUS(xDiff, yDiff, R2)) continue;
    146 
    147             psF32 pDiff = *vPix - sky;
    148             psF32 wDiff = *vWgt;
    149 
    150             // skip pixels below specified significance level.  for a PSFs, this
    151             // over-weights the wings of bright stars compared to those of faint stars.
    152             // for the estimator used for extended source analysis (where the window
    153             // function is allowed to be arbitrarily large), we need to clip to avoid
    154             // negative second moments.
    155             if (PS_SQR(pDiff) < minSN2*wDiff) continue; //
    156             if ((minSN > 0.0) && (pDiff < 0)) continue; //
    157 
    158             // Apply a Gaussian window function.  Be careful with the window function.  S/N
    159             // weighting over weights the sky for faint sources
    160             if (sigma > 0.0) {
    161                 // XXX a lot of extra flops; can we pre-calculate?
    162                 psF32 z  = (PS_SQR(xDiff) + PS_SQR(yDiff))*rsigma2;
    163                 assert (z >= 0.0);
    164                 psF32 weight  = exp(-z);
    165 
    166                 wDiff *= weight;
    167                 pDiff *= weight;
    168             }
    169 
    170             Var += wDiff;
    171             Sum += pDiff;
    172 
    173             psF32 xWght = xDiff * pDiff;
    174             psF32 yWght = yDiff * pDiff;
    175 
    176             X1  += xWght;
    177             Y1  += yWght;
    178 
    179             peakPixel = PS_MAX (*vPix, peakPixel);
    180             numPixels++;
    181         }
    182     }
    183 
    184     // if we have less than (1/4) of the possible pixels (in circle or box), force a retry
    185     int minPixels = PS_MIN(0.75*R2, source->pixels->numCols*source->pixels->numRows/4.0);
    186 
    187     // XXX EAM - the limit is a bit arbitrary.  make it user defined?
    188     if ((numPixels < minPixels) || (Sum <= 0)) {
    189         psTrace ("psModules.objects", 3, "insufficient valid pixels (%d vs %d; %f) for source\n", numPixels, minPixels, Sum);
    190         return (false);
    191     }
    192 
    193     // calculate the first moment.
    194     float Mx = X1/Sum;
    195     float My = Y1/Sum;
    196     if ((fabs(Mx) > radius) || (fabs(My) > radius)) {
    197         psTrace ("psModules.objects", 3, "large centroid swing; invalid peak %d, %d\n", source->peak->x, source->peak->y);
    198         return (false);
    199     }
    200 
    201     psTrace ("psModules.objects", 5, "sky: %f  Mx: %f  My: %f  Sum: %f  X1: %f  Y1: %f  Npix: %d\n", sky, Mx, My, Sum, X1, Y1, numPixels);
    202 
    203     // add back offset of peak in primary image
    204     // also offset from pixel index to pixel coordinate
    205     // (the calculation above uses pixel index instead of coordinate)
    206     // 0.5 PIX: moments are calculated using the pixel index and converted here to pixel coords
    207     source->moments->Mx = Mx + xOff + 0.5;
    208     source->moments->My = My + yOff + 0.5;
    209 
    210     source->moments->Sum = Sum;
    211     source->moments->SN  = Sum / sqrt(Var);
    212     source->moments->Peak = peakPixel;
    213     source->moments->nPixels = numPixels;
     96    // do 2 passes : the first pass should use a somewhat smaller radius and no sigma window to
     97    // get an unbiased (but probably noisy) centroid
     98    if (!pmSourceMomentsGetCentroid (source, 0.75*radius, 0.0, minSN, maskVal, source->peak->xf, source->peak->yf)) {
     99        return false;
     100    }
     101    // second pass applies the Gaussian window and uses the centroid from the first pass
     102    if (!pmSourceMomentsGetCentroid (source, radius, sigma, minSN, maskVal, source->moments->Mx, source->moments->My)) {
     103        return false;
     104    }
    214105
    215106    // Now calculate higher-order moments, using the above-calculated first moments to adjust coordinates
    … …  
    234125    Sum = 0.0;  // the second pass may include slightly different pixels, re-determine Sum
    235126
    236     // center of mass in subimage.  Note: the calculation below uses pixel index, so we do not
    237     // correct xCM, yCM to pixel coordinate here.
    238     psF32 xCM = Mx + xPeak; // coord of peak in subimage
    239     psF32 yCM = My + yPeak; // coord of peak in subimage
     127    // center of mass in subimage.  Note: the calculation below uses pixel index, so we correct
     128    // xCM, yCM from pixel coords to pixel index here.
     129    psF32 xCM = source->moments->Mx - 0.5 - source->pixels->col0; // coord of peak in subimage
     130    psF32 yCM = source->moments->My - 0.5 - source->pixels->row0; // coord of peak in subimage
    240131
    241132    for (psS32 row = 0; row < source->pixels->numRows ; row++) {
    … …  
    263154            // radius is just a function of (xDiff, yDiff)
    264155            psF32 r2  = PS_SQR(xDiff) + PS_SQR(yDiff);
    265             psF32 r  = sqrt(r2);
    266             if (r > radius) continue;
     156            if (r2 > R2) continue;
    267157
    268158            psF32 fDiff = *vPix - sky;
    … …  
    274164            // stars.
    275165            if (PS_SQR(pDiff) < minSN2*wDiff) continue;
    276             // if (pDiff < 0) continue;
     166            if ((minSN > 0.0) && (pDiff < 0)) continue; //
    277167
    278168            // Apply a Gaussian window function.  Be careful with the window function.  S/N
    279169            // weighting over weights the sky for faint sources
    280170            if (sigma > 0.0) {
    281                 // XXX a lot of extra flops; can we do pre-calculate?
    282                 // XXX we were re-calculating r2 (maybe the compiler caught this?)
    283                 // psF32 z  = (PS_SQR(xDiff) + PS_SQR(yDiff))*rsigma2;
    284                 psF32 z  = r2 * rsigma2;
     171                psF32 z = r2 * rsigma2;
    285172                assert (z >= 0.0);
    286173                psF32 weight  = exp(-z);
    … …  
    292179            Sum += pDiff;
    293180
    294 # if (1)
    295 # if (0)
    296             if (fDiff < 0) continue;
    297 # endif
     181            // Kron Flux uses the 1st radial moment (NOT Gaussian windowed?)
     182            psF32 r = sqrt(r2);
    298183            psF32 rf = r * fDiff;
    299184            psF32 rh = sqrt(r) * fDiff;
    300185            psF32 rs = fDiff;
    301 # else
    302             psF32 rf = r * pDiff;
    303             psF32 rh = sqrt(r) * pDiff;
    304             psF32 rs = pDiff;
    305 # endif
    306186
    307187            psF32 x = xDiff * pDiff;
    … …  
    363243
    364244    // Calculate the Kron magnitude (make this block optional?)
    365     // float radKron = 2.5*source->moments->Mrf;
    366245    float radKinner = 1.0*source->moments->Mrf;
    367246    float radKron   = 2.5*source->moments->Mrf;
    … …  
    397276            // radKron is just a function of (xDiff, yDiff)
    398277            psF32 r2  = PS_SQR(xDiff) + PS_SQR(yDiff);
    399             psF32 r  = sqrt(r2);
    400278
    401279            psF32 pDiff = *vPix - sky;
    … …  
    407285            if (PS_SQR(pDiff) < minSN2*wDiff) continue;
    408286
     287            psF32 r  = sqrt(r2);
    409288            if (r < radKron) {
    410289                Sum += pDiff;
    … …  
    434313
    435314    psTrace ("psModules.objects", 3, "peak %f %f (%f = %f) Mx: %f  My: %f  Sum: %f  Mxx: %f  Mxy: %f  Myy: %f  sky: %f  Npix: %d\n",
    436              source->peak->xf, source->peak->yf, source->peak->flux, source->peak->SN, source->moments->Mx,   source->moments->My, Sum, source->moments->Mxx,   source->moments->Mxy,   source->moments->Myy, sky, numPixels);
     315             source->peak->xf, source->peak->yf, source->peak->flux, source->peak->SN, source->moments->Mx,   source->moments->My, Sum, source->moments->Mxx,   source->moments->Mxy,   source->moments->Myy, sky, source->moments->nPixels);
    437316
    438317    return(true);
    439318}
     319
     320bool pmSourceMomentsGetCentroid(pmSource *source, psF32 radius, psF32 sigma, psF32 minSN, psImageMaskType maskVal, float xGuess, float yGuess) {
     321
     322    // First Pass: calculate the first moments (these are subtracted from the coordinates below)
     323    // Sum = SUM (z - sky)
     324    // X1  = SUM (x - xc)*(z - sky)
     325    // .. etc
     326
     327    psF32 sky = 0.0;
     328
     329    psF32 peakPixel = -PS_MAX_F32;
     330    psS32 numPixels = 0;
     331    psF32 Sum = 0.0;
     332    psF32 Var = 0.0;
     333    psF32 X1 = 0.0;
     334    psF32 Y1 = 0.0;
     335    psF32 R2 = PS_SQR(radius);
     336    psF32 minSN2 = PS_SQR(minSN);
     337    psF32 rsigma2 = 0.5 / PS_SQR(sigma);
     338
     339    float xPeak = xGuess - source->pixels->col0; // coord of peak in subimage
     340    float yPeak = yGuess - source->pixels->row0; // coord of peak in subimage
     341
     342    // we are guaranteed to have a valid pixel and variance at this location (right? right?)
     343    // psF32 weightNorm = source->pixels->data.F32[yPeak][xPeak] / sqrt (source->variance->data.F32[yPeak][xPeak]);
     344    // psAssert (isfinite(source->pixels->data.F32[yPeak][xPeak]), "peak must be on valid pixel");
     345    // psAssert (isfinite(source->variance->data.F32[yPeak][xPeak]), "peak must be on valid pixel");
     346    // psAssert (source->variance->data.F32[yPeak][xPeak] > 0, "peak must be on valid pixel");
     347
     348    // the moments [Sum(x*f) / Sum(f)] are calculated in pixel index values, and should
     349    // not depend on the fractional pixel location of the source.  However, the aperture
     350    // (radius) and the Gaussian window (sigma) depend subtly on the fractional pixel
     351    // position of the expected centroid
     352
     353    for (psS32 row = 0; row < source->pixels->numRows ; row++) {
     354
     355        psF32 yDiff = row + 0.5 - yPeak;
     356        if (fabs(yDiff) > radius) continue;
     357
     358        psF32 *vPix = source->pixels->data.F32[row];
     359        psF32 *vWgt = source->variance->data.F32[row];
     360        psImageMaskType *vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[row];
     361
     362        for (psS32 col = 0; col < source->pixels->numCols ; col++, vPix++, vWgt++) {
     363            if (vMsk) {
     364                if (*vMsk & maskVal) {
     365                    vMsk++;
     366                    continue;
     367                }
     368                vMsk++;
     369            }
     370            if (isnan(*vPix)) continue;
     371
     372            psF32 xDiff = col + 0.5 - xPeak;
     373            if (fabs(xDiff) > radius) continue;
     374
     375            // radius is just a function of (xDiff, yDiff)
     376            psF32 r2  = PS_SQR(xDiff) + PS_SQR(yDiff);
     377            if (r2 > R2) continue;
     378
     379            psF32 pDiff = *vPix - sky;
     380            psF32 wDiff = *vWgt;
     381
     382            // skip pixels below specified significance level.  for a PSFs, this
     383            // over-weights the wings of bright stars compared to those of faint stars.
     384            // for the estimator used for extended source analysis (where the window
     385            // function is allowed to be arbitrarily large), we need to clip to avoid
     386            // negative second moments.
     387            if (PS_SQR(pDiff) < minSN2*wDiff) continue; //
     388            if ((minSN > 0.0) && (pDiff < 0)) continue; //
     389
     390            // Apply a Gaussian window function.  Be careful with the window function.  S/N
     391            // weighting over weights the sky for faint sources
     392            if (sigma > 0.0) {
     393                psF32 z  = r2*rsigma2;
     394                assert (z >= 0.0);
     395                psF32 weight  = exp(-z);
     396
     397                wDiff *= weight;
     398                pDiff *= weight;
     399            }
     400
     401            Var += wDiff;
     402            Sum += pDiff;
     403
     404            psF32 xWght = xDiff * pDiff;
     405            psF32 yWght = yDiff * pDiff;
     406
     407            X1  += xWght;
     408            Y1  += yWght;
     409
     410            peakPixel = PS_MAX (*vPix, peakPixel);
     411            numPixels++;
     412        }
     413    }
     414
     415    // if we have less than (1/4) of the possible pixels (in circle or box), force a retry
     416    int minPixels = PS_MIN(0.75*R2, source->pixels->numCols*source->pixels->numRows/4.0);
     417
     418    // XXX EAM - the limit is a bit arbitrary.  make it user defined?
     419    if ((numPixels < minPixels) || (Sum <= 0)) {
     420        psTrace ("psModules.objects", 3, "insufficient valid pixels (%d vs %d; %f) for source\n", numPixels, minPixels, Sum);
     421        return (false);
     422    }
     423
     424    // calculate the first moment.
     425    float Mx = X1/Sum;
     426    float My = Y1/Sum;
     427    if ((fabs(Mx) > radius) || (fabs(My) > radius)) {
     428        psTrace ("psModules.objects", 3, "large centroid swing; invalid peak %d, %d\n", source->peak->x, source->peak->y);
     429        return (false);
     430    }
     431
     432    psTrace ("psModules.objects", 5, "sky: %f  Mx: %f  My: %f  Sum: %f  X1: %f  Y1: %f  Npix: %d\n", sky, Mx, My, Sum, X1, Y1, numPixels);
     433
     434    // add back offset of peak in primary image
     435    // also offset from pixel index to pixel coordinate
     436    // (the calculation above uses pixel index instead of coordinate)
     437    // 0.5 PIX: moments are calculated using the pixel index and converted here to pixel coords
     438
     439    // we only update the centroid if the position is not supplied from elsewhere
     440    bool skipCentroid = false;
     441    skipCentroid |= (source->mode  & PM_SOURCE_MODE_EXTERNAL); // skip externally supplied positions
     442    skipCentroid |= (source->mode2 & PM_SOURCE_MODE2_MATCHED); // skip sources defined by other image positions
     443
     444    if (skipCentroid) {
     445        source->moments->Mx = source->peak->xf;
     446        source->moments->My = source->peak->yf;
     447    } else {
     448        source->moments->Mx = Mx + xGuess;
     449        source->moments->My = My + yGuess;
     450    }
     451
     452    source->moments->Sum = Sum;
     453    source->moments->SN  = Sum / sqrt(Var);
     454    source->moments->Peak = peakPixel;
     455    source->moments->nPixels = numPixels;
     456
     457    return true;
     458}

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

    r29004 r29546  
    6262
    6363/**
    64     this function is used to calculate the three defined source magnitudes:
    65     - apMag  : only if S/N > AP_MIN_SN
    66              : is optionally corrected for curve-of-growth if:
    67         - the source is a STAR (PSF)
    68         - the option is selected (mode & PM_SOURCE_PHOT_GROWTH)
    69     - psfMag : all sources with non-NULL modelPSF
    70              : is optionally corrected for aperture residual if:
    71         - the source is a STAR (PSF)
    72         - the option is selected (mode & PM_SOURCE_PHOT_APCORR)
    73 
    74     - extMag : all sources with non-NULL modelEXT
     64   this function is used to calculate the three defined source magnitudes:
     65   - apMag  : only if S/N > AP_MIN_SN
     66   : is optionally corrected for curve-of-growth if:
     67   - the source is a STAR (PSF)
     68   - the option is selected (mode & PM_SOURCE_PHOT_GROWTH)
     69   - psfMag : all sources with non-NULL modelPSF
     70   : is optionally corrected for aperture residual if:
     71   - the source is a STAR (PSF)
     72   - the option is selected (mode & PM_SOURCE_PHOT_APCORR)
     73
     74   - extMag : all sources with non-NULL modelEXT
    7575**/
    7676
    7777// XXX masked region should be (optionally) elliptical
     78// if mode is PM_SOURCE_PHOT_PSFONLY, we skip all other magnitudes
    7879bool pmSourceMagnitudes (pmSource *source, pmPSF *psf, pmSourcePhotometryMode mode, psImageMaskType maskVal, psImageMaskType markVal)
    7980{
    … …  
    8889    pmModel *model;
    8990
    90     source->psfMag = NAN;
    91     source->extMag = NAN;
    92     source->errMag = NAN;
    93     source->apMag  = NAN;
     91    source->psfMag    = NAN;
     92    source->extMag    = NAN;
     93    source->errMag    = NAN;
     94    source->apMag     = NAN;
     95    source->apMagRaw  = NAN;
     96    source->apFlux    = NAN;
     97    source->apFluxErr = NAN;
    9498
    9599    // we must have a valid model
    … …  
    114118    // measure PSF model photometry
    115119    // XXX TEST: do not use flux scale
     120    // XXX NOTE: turn this back on?
    116121    if (0 && psf->FluxScale) {
    117122        // the source peak pixel is guaranteed to be on the image, and only minimally different from the source center
    … …  
    127132    }
    128133
     134    if (mode == PM_SOURCE_PHOT_PSFONLY) {
     135        return true;
     136    }
     137
    129138    // if we have a collection of model fits, check if one of them is a pointer to modelEXT
    130139    if (source->modelFits) {
    … …  
    148157
    149158    // for PSFs, correct both apMag and psfMag to same system, consistent with infinite flux star in aperture RADIUS
     159    // XXX add a flag for "ap_mag is corrected?"
    150160    if ((mode & PM_SOURCE_PHOT_APCORR) && isPSF && psf && psf->ApTrend) {
    151161        // the source peak pixel is guaranteed to be on the image, and only minimally different from the source center
    … …  
    181191    // if we are measuring aperture photometry and applying the growth correction,
    182192    // we need to shift the flux in the selected pixels (but not the mask)
    183     psImage *flux = NULL, *mask = NULL; // Star flux and mask images, to photometer
     193    psImage *flux = NULL;
     194    psImage *variance = NULL;
     195    psImage *mask = NULL; // Star flux and mask images, to photometer
    184196    if (mode & PM_SOURCE_PHOT_INTERP) {
    185197        float dx = 0.5 - x + (int)x;
    … …  
    198210    } else {
    199211        flux = source->pixels;
     212        variance = source->variance;
    200213        mask = source->maskObj;
    201214    }
    202215
    203216    // measure object aperture photometry
    204     status = pmSourcePhotometryAper  (&source->apMagRaw, model, flux, mask, maskVal);
     217    status = pmSourcePhotometryAperSource (source, model, flux, variance, mask, maskVal);
    205218    if (!status) {
    206219        psTrace ("psModules.objects", 3, "fail mag : bad Ap Mag");
    … …  
    214227        if (psf->growth && (mode & PM_SOURCE_PHOT_GROWTH)) {
    215228            source->apMag = source->apMagRaw + pmGrowthCurveCorrect (psf->growth, source->apRadius);
     229            // XXX correct the apFlux?
    216230        }
    217231        if (mode & PM_SOURCE_PHOT_APCORR) {
    … …  
    233247
    234248/*
    235 aprMag' - fitMag = flux*skySat + r^2*rflux*skyBias + ApTrend(x,y)
    236 (aprMag - flux*skySat - r^2*rflux*skyBias) - fitMAg = ApTrend(x,y)
    237 (aprMag - flux*skySat - r^2*rflux*skyBias) = fitMAg + ApTrend(x,y)
     249  aprMag' - fitMag = flux*skySat + r^2*rflux*skyBias + ApTrend(x,y)
     250  (aprMag - flux*skySat - r^2*rflux*skyBias) - fitMAg = ApTrend(x,y)
     251  (aprMag - flux*skySat - r^2*rflux*skyBias) = fitMAg + ApTrend(x,y)
    238252
    239253*/
    … …  
    267281
    268282// return source aperture magnitude
    269 bool pmSourcePhotometryAper (float *apMag, pmModel *model, psImage *image, psImage *mask, psImageMaskType maskVal)
     283bool pmSourcePhotometryAperSource (pmSource *source, pmModel *model, psImage *image, psImage *variance, psImage *mask, psImageMaskType maskVal)
     284{
     285    PS_ASSERT_PTR_NON_NULL(source, false);
     286    PS_ASSERT_PTR_NON_NULL(image, false);
     287    PS_ASSERT_PTR_NON_NULL(mask, false);
     288
     289    if (DO_SKY) {
     290        PS_ASSERT_PTR_NON_NULL(model, false);
     291    }
     292
     293    bool status;
     294    status = pmSourcePhotometryAper(&source->apMagRaw, &source->apFlux, &source->apFluxErr, model, image, variance, mask, maskVal);
     295
     296    return status;
     297}
     298
     299// return source aperture magnitude
     300bool pmSourcePhotometryAper (float *apMag, float *apFluxOut, float *apFluxErr, pmModel *model, psImage *image, psImage *variance, psImage *mask, psImageMaskType maskVal)
    270301{
    271302    PS_ASSERT_PTR_NON_NULL(apMag, false);
    … …  
    277308    }
    278309
    279     float apSum = 0;
    280310    float sky = 0;
    281     *apMag = NAN;
     311    float apFlux = 0;
     312    float apFluxVar = 0;
    282313
    283314    if (DO_SKY) {
    284315        sky = model->params->data.F32[PM_PAR_SKY];
    285     } else {
    286         sky = 0;
    287316    }
    288317
    289318    psF32 **imData = image->data.F32;
    290319    psImageMaskType **mkData = mask->data.PS_TYPE_IMAGE_MASK_DATA;
    291     int nAperPix = 0;
    292 
    293     // measure apMag
     320    psF32 **varData = (variance) ? variance->data.F32 : image->data.F32; // if variance is not supplied, assume gain of 1.0, no read noise
     321
     322    // measure apFlux and apFluxVar, save apMag if not NAN
     323    // XXX note that these fluxes/mags are uncorrected for masked pixels
     324    // XXX raise a bit if the aperture has a masked pixel (not marked)?
    294325    for (int iy = 0; iy < image->numRows; iy++) {
    295326        for (int ix = 0; ix < image->numCols; ix++) {
    296             if (mkData[iy][ix] & maskVal)
    297                 continue;
    298             apSum += imData[iy][ix] - sky;
    299             nAperPix ++;
    300             // fprintf (stderr, "aper: %d %d  %f  %f  %f\n", ix, iy, sky, imData[iy][ix], apSum);
    301         }
    302     }
    303     if (apSum <= 0) {
     327            if (mkData[iy][ix] & maskVal) continue;
     328            apFlux += imData[iy][ix] - sky;
     329            apFluxVar += varData[iy][ix];
     330        }
     331    }
     332    if (apFluxOut) *apFluxOut = apFlux;
     333    if (apFluxErr) *apFluxErr = sqrt(apFluxVar);
     334
     335    if (apFlux <= 0) {
    304336        *apMag = NAN;
    305         return true;
    306     }
    307 
    308     *apMag = -2.5*log10(apSum);
     337    } else {
     338        *apMag = -2.5*log10(apFlux);
     339    }
    309340    return true;
    310341}
    … …  
    419450    }
    420451
     452    // NOTE: until 2010.10.01, these measurements included a 3sigma-per-pixel significance
     453    // this followed what we understood as the definition given to us
     454    // by Armin, but it always seemed a poor idea -- a faint source is unlikely to have any 3sigma pixels.
     455    // changed to remove the per-pixel filter.
     456
    421457    for (int iy = 0; iy < flux->numRows; iy++) {
    422458        for (int ix = 0; ix < flux->numCols; ix++) {
    423459            // only count up the stats in the unmarked region (ie, the aperture)
     460            // skip the marked pixels; these are not relevant
    424461            if (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & markVal) {
    425462                continue;
    … …  
    430467            }
    431468
    432             float SN = flux->data.F32[iy][ix] / sqrt(variance->data.F32[iy][ix]);
    433 
    434             if (SN > +FLUX_LIMIT) {
     469            float value = flux->data.F32[iy][ix];
     470
     471            if (value > 0.0) {
    435472                nGood ++;
    436                 fGood += fabs(flux->data.F32[iy][ix]);
    437             }
    438 
    439             if (SN < -FLUX_LIMIT) {
     473                fGood += fabs(value);
     474            } else {
    440475                nBad ++;
    441                 fBad += fabs(flux->data.F32[iy][ix]);
     476                fBad += fabs(value);
    442477            }
    443478        }
    … …  
    613648
    614649            switch (term) {
    615             case 0:
     650              case 0:
    616651                factor = 1;
    617652                break;
    618             case 1:
     653              case 1:
    619654                factor = xi + Pi->col0;
    620655                break;
    621             case 2:
     656              case 2:
    622657                factor = yi + Pi->row0;
    623658                break;
    624             default:
     659              default:
    625660                psAbort("invalid term for pmSourceWeight");
    626661            }
    … …  
    691726
    692727            switch (term) {
    693             case 0:
     728              case 0:
    694729                factor = 1;
    695730                break;
    696             case 1:
     731              case 1:
    697732                factor = xi + Pi->col0;
    698733                break;
    699             case 2:
     734              case 2:
    700735                factor = yi + Pi->row0;
    701736                break;
    702             default:
     737              default:
    703738                psAbort("invalid term for pmSourceWeight");
    704739            }

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

    r29004 r29546  
    3535    PM_SOURCE_PHOT_INTERP    = 0x0008,
    3636    PM_SOURCE_PHOT_DIFFSTATS = 0x0010,
     37    PM_SOURCE_PHOT_PSFONLY   = 0x0020,
    3738} pmSourcePhotometryMode;
    3839
    … …  
    4445
    4546bool pmSourcePhotometryAper(
    46     float   *apMag,                     ///< aperture flux magnitude
     47    float *apMag,
     48    float *apFluxOut,
     49    float *apFluxErr,
    4750    pmModel *model,                     ///< model used for photometry
    4851    psImage *image,                     ///< image pixels to be used
     52    psImage *variance,                  ///< variance pixels to be used
     53    psImage *mask,                      ///< mask of pixels to ignore
     54    psImageMaskType maskVal             ///< Value to mask
     55);
     56
     57bool pmSourcePhotometryAperSource(
     58    pmSource *source,                   ///< aperture flux magnitude
     59    pmModel *model,                     ///< model used for photometry
     60    psImage *image,                     ///< image pixels to be used
     61    psImage *variance,                  ///< variance pixels to be used
    4962    psImage *mask,                      ///< mask of pixels to ignore
    5063    psImageMaskType maskVal             ///< Value to mask
Note: See TracChangeset for help on using the changeset viewer.