Index: trunk/psModules/src/objects/pmGrowthCurveGenerate.c =================================================================== --- trunk/psModules/src/objects/pmGrowthCurveGenerate.c (revision 29044) +++ trunk/psModules/src/objects/pmGrowthCurveGenerate.c (revision 29546) @@ -197,5 +197,5 @@ // mask the given aperture and measure the apMag psImageKeepCircle (mask, xc, yc, radius, "OR", markVal); - if (!pmSourcePhotometryAper (&apMag, model, pixels, mask, maskVal)) { + if (!pmSourcePhotometryAper (&apMag, NULL, NULL, model, pixels, NULL, mask, maskVal)) { psFree (growth); psFree (view); Index: trunk/psModules/src/objects/pmPCMdata.c =================================================================== --- trunk/psModules/src/objects/pmPCMdata.c (revision 29044) +++ trunk/psModules/src/objects/pmPCMdata.c (revision 29546) @@ -104,6 +104,6 @@ assert (source->psfImage); // XXX build if needed? - int x0 = source->peak->xf - source->psfImage->col0; - int y0 = source->peak->yf - source->psfImage->row0; + int x0 = source->peak->x - source->psfImage->col0; + int y0 = source->peak->y - source->psfImage->row0; // need to decide on the size: dynamically? statically? Index: trunk/psModules/src/objects/pmSource.c =================================================================== --- trunk/psModules/src/objects/pmSource.c (revision 29044) +++ trunk/psModules/src/objects/pmSource.c (revision 29546) @@ -132,18 +132,22 @@ // default values are NAN - source->psfMag = NAN; - source->psfFlux = NAN; - source->psfFluxErr = NAN; - source->extMag = NAN; - source->errMag = NAN; - source->apMag = NAN; - source->sky = NAN; - source->skyErr = NAN; - source->pixWeightNotBad = NAN; + source->psfMag = NAN; + source->psfFlux = NAN; + source->psfFluxErr = NAN; + source->extMag = NAN; + source->errMag = NAN; + source->apMag = NAN; + source->apMagRaw = NAN; + source->apRadius = NAN; + source->apFlux = NAN; + source->apFluxErr = NAN; + source->sky = NAN; + source->skyErr = NAN; + source->pixWeightNotBad = NAN; source->pixWeightNotPoor = NAN; - source->psfChisq = NAN; - source->crNsigma = NAN; - source->extNsigma = NAN; + source->psfChisq = NAN; + source->crNsigma = NAN; + source->extNsigma = NAN; psTrace("psModules.objects", 10, "---- end ----\n"); @@ -644,8 +648,7 @@ // The following determinations require the use of moments if (!(source->mode & noMoments)) { - // likely defect (too small to be stellar) (push out to 3 sigma) - // low S/N objects which are small are probably stellar - // XXX these limits are quite arbitrary - if (sigX < 0.05 || sigY < 0.05) { + // likely defect (bright, but too small to be stellar) + // XXX eliminate the classification? + if ((source->moments->SN > 10) && (sigX < 0.05 || sigY < 0.05)) { source->type = PM_SOURCE_TYPE_DEFECT; source->mode |= PM_SOURCE_MODE_DEFECT; Index: trunk/psModules/src/objects/pmSource.h =================================================================== --- trunk/psModules/src/objects/pmSource.h (revision 29044) +++ trunk/psModules/src/objects/pmSource.h (revision 29546) @@ -76,4 +76,6 @@ float apMagRaw; ///< raw mag in given aperture float apRadius; ///< radius for aperture magnitude + float apFlux; ///< apFlux corresponding to psfMag or extMag (depending on type) + float apFluxErr; ///< apFluxErr corresponding to psfMag or extMag (depending on type) float pixWeightNotBad; ///< PSF-weighted coverage of unmasked (not BAD) pixels @@ -234,4 +236,20 @@ ); +/** pmSourceMoments() + * + * Measure 1st moments for the given source, using the peak coordinates as the initial + * source location. The resulting moment values are applied to the source.moments + * entry. The moments are measured within the given circular radius of the source.peak + * coordinates. The return value indicates the success (TRUE) of the operation. + * + */ +bool pmSourceMomentsGetCentroid( + pmSource *source, + psF32 radius, + psF32 sigma, + psF32 minSN, + psImageMaskType maskVal, + float xGuess, float yGuess); + pmModel *pmSourceGetModel (bool *isPSF, const pmSource *source); Index: trunk/psModules/src/objects/pmSourceIO.c =================================================================== --- trunk/psModules/src/objects/pmSourceIO.c (revision 29044) +++ trunk/psModules/src/objects/pmSourceIO.c (revision 29546) @@ -1001,4 +1001,7 @@ sources = pmSourcesRead_CMF_PS1_V2 (file->fits, hdu->header); } + if (!strcmp (exttype, "PS1_V3")) { + sources = pmSourcesRead_CMF_PS1_V3 (file->fits, hdu->header); + } if (!strcmp (exttype, "PS1_SV1")) { sources = pmSourcesRead_CMF_PS1_SV1 (file->fits, hdu->header); Index: trunk/psModules/src/objects/pmSourceIO.h =================================================================== --- trunk/psModules/src/objects/pmSourceIO.h (revision 29044) +++ trunk/psModules/src/objects/pmSourceIO.h (revision 29546) @@ -75,4 +75,5 @@ psArray *pmSourcesRead_CMF_PS1_V1 (psFits *fits, psMetadata *header); psArray *pmSourcesRead_CMF_PS1_V2 (psFits *fits, psMetadata *header); +psArray *pmSourcesRead_CMF_PS1_V3 (psFits *fits, psMetadata *header); psArray *pmSourcesRead_CMF_PS1_SV1 (psFits *fits, psMetadata *header); psArray *pmSourcesRead_CMF_PS1_DV1 (psFits *fits, psMetadata *header); Index: trunk/psModules/src/objects/pmSourceIO_CMF_PS1_DV2.c =================================================================== --- trunk/psModules/src/objects/pmSourceIO_CMF_PS1_DV2.c (revision 29044) +++ trunk/psModules/src/objects/pmSourceIO_CMF_PS1_DV2.c (revision 29546) @@ -185,6 +185,11 @@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_FLUX", PS_DATA_F32, "PSF fit instrumental magnitude", source->psfFlux); psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_FLUX_SIG",PS_DATA_F32, "Sigma of PSF instrumental magnitude", source->psfFluxErr); + psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG", PS_DATA_F32, "magnitude in standard aperture", source->apMag); + psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RAW", PS_DATA_F32, "magnitude in real aperture", source->apMagRaw); psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RADIUS", PS_DATA_F32, "radius used for aperture mags", apRadius); + psMetadataAdd (row, PS_LIST_TAIL, "AP_FLUX", PS_DATA_F32, "instrumental flux in standard aperture", source->apFlux); + psMetadataAdd (row, PS_LIST_TAIL, "AP_FLUX_SIG", PS_DATA_F32, "aperture flux error", source->apFluxErr); + psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude", peakMag); psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG", PS_DATA_F32, "PSF Magnitude using supplied calibration", calMag); @@ -214,4 +219,20 @@ psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XY", PS_DATA_F32, "second moments (X*Y)", mxy); psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_YY", PS_DATA_F32, "second moments (Y*Y)", myy); + + float Mrf = source->moments ? source->moments->Mrf : NAN; + float Mrh = source->moments ? source->moments->Mrh : NAN; + float Krf = source->moments ? source->moments->KronFlux : NAN; + float dKrf = source->moments ? source->moments->KronFluxErr : NAN; + + float Kinner = source->moments ? source->moments->KronFinner : NAN; + float Kouter = source->moments ? source->moments->KronFouter : NAN; + + psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_R1", PS_DATA_F32, "first radial moment", Mrf); + psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_RH", PS_DATA_F32, "half radial moment", Mrh); + psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX", PS_DATA_F32, "Kron Flux (in 2.5 R1)", Krf); + psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_ERR", PS_DATA_F32, "Kron Flux Error", dKrf); + + psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_INNER", PS_DATA_F32, "Kron Flux (in 1.0 R1)", Kinner); + psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_OUTER", PS_DATA_F32, "Kron Flux (in 4.0 R1)", Kouter); psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NPOS", PS_DATA_S32, "nPos (n pix > 3 sigma)", diffStats.nGood); Index: trunk/psModules/src/objects/pmSourceIO_CMF_PS1_V3.c =================================================================== --- trunk/psModules/src/objects/pmSourceIO_CMF_PS1_V3.c (revision 29044) +++ trunk/psModules/src/objects/pmSourceIO_CMF_PS1_V3.c (revision 29546) @@ -66,5 +66,5 @@ psF32 xPos, yPos; psF32 xErr, yErr; - psF32 errMag, chisq, apRadius; + psF32 chisq, apRadius; psS32 nPix, nDOF; @@ -84,4 +84,14 @@ } float zeroptErr = psMetadataLookupF32 (&status2, imageHeader, "ZPT_ERR"); + + // we need a measure of the image quality (FWHM) for this image, in order to get the positional errors + float fwhmMajor = psMetadataLookupF32(&status1, readout->analysis, "FWHM_MAJ"); + if (!status1) { + fwhmMajor = psMetadataLookupF32(&status1, readout->analysis, "IQ_FW1"); + } + float fwhmMinor = psMetadataLookupF32(&status1, readout->analysis, "FWHM_MIN"); + if (!status1) { + fwhmMinor = psMetadataLookupF32(&status1, readout->analysis, "IQ_FW2"); + } // if the sequence is defined, write these in seq order; otherwise @@ -124,7 +134,6 @@ yErr = dPAR[PM_PAR_YPOS]; } else { - // in linear-fit mode, there is no error on the centroid - xErr = source->peak->dx; - yErr = source->peak->dy; + xErr = fwhmMajor * source->errMag / 2.35; + yErr = fwhmMinor * source->errMag / 2.35; } if (isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX])) { @@ -139,10 +148,21 @@ nPix = model->nPix; apRadius = source->apRadius; - errMag = model->dparams->data.F32[PM_PAR_I0] / model->params->data.F32[PM_PAR_I0]; } else { - xPos = source->peak->xf; - yPos = source->peak->yf; - xErr = source->peak->dx; - yErr = source->peak->dy; + bool useMoments = true; + useMoments = (useMoments && source->moments); // can't if there are no moments + useMoments = (useMoments && source->moments->nPixels); // can't if the moments were not measured + useMoments = (useMoments && !(source->mode && PM_SOURCE_MODE_MOMENTS_FAILURE)); // can't if the moments failed... + + if (source->moments) { + xPos = source->moments->Mx; + yPos = source->moments->My; + xErr = fwhmMajor * source->errMag / 2.35; + yErr = fwhmMinor * source->errMag / 2.35; + } else { + xPos = source->peak->xf; + yPos = source->peak->yf; + xErr = source->peak->dx; + yErr = source->peak->dy; + } axes.major = NAN; axes.minor = NAN; @@ -152,5 +172,4 @@ nPix = 0; apRadius = NAN; - errMag = NAN; } @@ -168,10 +187,10 @@ psMetadataAdd (row, PS_LIST_TAIL, "X_PSF", PS_DATA_F32, "PSF x coordinate", xPos); psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF", PS_DATA_F32, "PSF y coordinate", yPos); - 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 - 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 + psMetadataAdd (row, PS_LIST_TAIL, "X_PSF_SIG", PS_DATA_F32, "Sigma in PSF x coordinate", xErr); + psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF_SIG", PS_DATA_F32, "Sigma in PSF y coordinate", yErr); psMetadataAdd (row, PS_LIST_TAIL, "POSANGLE", PS_DATA_F32, "position angle at source (degrees)", posAngle*PS_DEG_RAD); psMetadataAdd (row, PS_LIST_TAIL, "PLTSCALE", PS_DATA_F32, "plate scale at source (arcsec/pixel)", pltScale*PS_DEG_RAD*3600.0); psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG", PS_DATA_F32, "PSF fit instrumental magnitude", source->psfMag); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG_SIG", PS_DATA_F32, "Sigma of PSF instrumental magnitude", errMag); + psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG_SIG", PS_DATA_F32, "Sigma of PSF instrumental magnitude", source->errMag); psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_FLUX", PS_DATA_F32, "PSF fit instrumental flux (counts)", source->psfFlux); psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_FLUX_SIG",PS_DATA_F32, "Sigma of PSF instrumental flux", source->psfFluxErr); @@ -179,9 +198,13 @@ psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RAW", PS_DATA_F32, "magnitude in reported aperture", source->apMagRaw); psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RADIUS", PS_DATA_F32, "radius used for aperture mags", apRadius); - psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude", peakMag); + psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG", PS_DATA_F32, "PSF Magnitude using supplied calibration", calMag); psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG_SIG", PS_DATA_F32, "measured scatter of zero point calibration", zeroptErr); + + // NOTE: RA & DEC (both double) need to be on an 8-byte boundary... psMetadataAdd (row, PS_LIST_TAIL, "RA_PSF", PS_DATA_F64, "PSF RA coordinate (degrees)", ptSky.r*PS_DEG_RAD); psMetadataAdd (row, PS_LIST_TAIL, "DEC_PSF", PS_DATA_F64, "PSF DEC coordinate (degrees)", ptSky.d*PS_DEG_RAD); + + psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude", peakMag); psMetadataAdd (row, PS_LIST_TAIL, "SKY", PS_DATA_F32, "Sky level", source->sky); psMetadataAdd (row, PS_LIST_TAIL, "SKY_SIGMA", PS_DATA_F32, "Sigma of sky level", source->skyErr); @@ -204,4 +227,18 @@ float Myy = source->moments ? source->moments->Myy : NAN; + psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XX", PS_DATA_F32, "second moments (X^2)", Mxx); + psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XY", PS_DATA_F32, "second moments (X*Y)", Mxy); + psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_YY", PS_DATA_F32, "second moments (Y*Y)", Myy); + + float M_c3 = source->moments ? 1.0*source->moments->Mxxx - 3.0*source->moments->Mxyy : NAN; + float M_s3 = source->moments ? 3.0*source->moments->Mxxy - 1.0*source->moments->Myyy : NAN; + float M_c4 = source->moments ? 1.0*source->moments->Mxxxx - 6.0*source->moments->Mxxyy + 1.0*source->moments->Myyyy : NAN; + float M_s4 = source->moments ? 4.0*source->moments->Mxxxy - 4.0*source->moments->Mxyyy : NAN; + + psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M3C", PS_DATA_F32, "third momemt cos theta", M_c3); + psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M3S", PS_DATA_F32, "third momemt sin theta", M_s3); + psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M4C", PS_DATA_F32, "fourth momemt cos theta", M_c4); + psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M4S", PS_DATA_F32, "fourth momemt sin theta", M_s4); + float Mrf = source->moments ? source->moments->Mrf : NAN; float Mrh = source->moments ? source->moments->Mrh : NAN; @@ -212,18 +249,4 @@ float Kouter = source->moments ? source->moments->KronFouter : NAN; - float M_c3 = source->moments ? 1.0*source->moments->Mxxx - 3.0*source->moments->Mxyy : NAN; - float M_s3 = source->moments ? 3.0*source->moments->Mxxy - 1.0*source->moments->Myyy : NAN; - float M_c4 = source->moments ? 1.0*source->moments->Mxxxx - 6.0*source->moments->Mxxyy + 1.0*source->moments->Myyyy : NAN; - float M_s4 = source->moments ? 4.0*source->moments->Mxxxy - 4.0*source->moments->Mxyyy : NAN; - - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XX", PS_DATA_F32, "second moments (X^2)", Mxx); - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XY", PS_DATA_F32, "second moments (X*Y)", Mxy); - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_YY", PS_DATA_F32, "second moments (Y*Y)", Myy); - - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M3C", PS_DATA_F32, "third momemt cos theta", M_c3); - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M3S", PS_DATA_F32, "third momemt sin theta", M_s3); - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M4C", PS_DATA_F32, "fourth momemt cos theta", M_c4); - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M4S", PS_DATA_F32, "fourth momemt sin theta", M_s4); - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_R1", PS_DATA_F32, "first radial moment", Mrf); psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_RH", PS_DATA_F32, "half radial moment", Mrh); @@ -236,4 +259,5 @@ psMetadataAdd (row, PS_LIST_TAIL, "FLAGS", PS_DATA_U32, "psphot analysis flags", source->mode); psMetadataAdd (row, PS_LIST_TAIL, "FLAGS2", PS_DATA_U32, "psphot analysis flags", source->mode2); + psMetadataAdd (row, PS_LIST_TAIL, "PADDING2", PS_DATA_S32, "more padding", 0); // XXX not sure how to get this : need to load Nimages with weight? Index: trunk/psModules/src/objects/pmSourceMasks.h =================================================================== --- trunk/psModules/src/objects/pmSourceMasks.h (revision 29044) +++ trunk/psModules/src/objects/pmSourceMasks.h (revision 29546) @@ -45,4 +45,5 @@ PM_SOURCE_MODE2_DIFF_WITH_SINGLE = 0x00000001, ///< diff source matched to a single positive detection PM_SOURCE_MODE2_DIFF_WITH_DOUBLE = 0x00000002, ///< diff source matched to positive detections in both images + PM_SOURCE_MODE2_MATCHED = 0x00000004, ///< diff source matched to positive detections in both images } pmSourceMode2; Index: trunk/psModules/src/objects/pmSourceMoments.c =================================================================== --- trunk/psModules/src/objects/pmSourceMoments.c (revision 29044) +++ trunk/psModules/src/objects/pmSourceMoments.c (revision 29546) @@ -64,4 +64,6 @@ void pmSourceMomentsSetVerbose(bool state){ beVerbose = state; } +// if mode & EXTERNAL or mode2 & MATCHED, do not re-calculate the centroid (use peak as centroid) + bool pmSourceMoments(pmSource *source, psF32 radius, psF32 sigma, psF32 minSN, psImageMaskType maskVal) { @@ -71,30 +73,13 @@ PS_ASSERT_FLOAT_LARGER_THAN(radius, 0.0, false); - // use sky from moments if defined, 0.0 otherwise - - // XXX this value comes from the sky model at the source center, and tends to over-estimate - // the sky in the vicinity of bright sources. we are better off assuming the model worked - // well: + // this function assumes the sky has been well-subtracted for the image psF32 sky = 0.0; + if (source->moments == NULL) { source->moments = pmMomentsAlloc(); } - // XXX if (source->moments == NULL) { - // XXX source->moments = pmMomentsAlloc(); - // XXX } else { - // XXX sky = source->moments->Sky; - // XXX } - - // First Pass: calculate the first moments (these are subtracted from the coordinates below) - // Sum = SUM (z - sky) - // X1 = SUM (x - xc)*(z - sky) - // .. etc - - psF32 peakPixel = -PS_MAX_F32; - psS32 numPixels = 0; + psF32 Sum = 0.0; psF32 Var = 0.0; - psF32 X1 = 0.0; - psF32 Y1 = 0.0; psF32 R2 = PS_SQR(radius); psF32 minSN2 = PS_SQR(minSN); @@ -109,107 +94,13 @@ // (int) so they can be used in the image index below. - int xOff = source->peak->x; - int yOff = source->peak->y; - int xPeak = source->peak->x - source->pixels->col0; // coord of peak in subimage - int yPeak = source->peak->y - source->pixels->row0; // coord of peak in subimage - - // we are guaranteed to have a valid pixel and variance at this location (right? right?) - // psF32 weightNorm = source->pixels->data.F32[yPeak][xPeak] / sqrt (source->variance->data.F32[yPeak][xPeak]); - // psAssert (isfinite(source->pixels->data.F32[yPeak][xPeak]), "peak must be on valid pixel"); - // psAssert (isfinite(source->variance->data.F32[yPeak][xPeak]), "peak must be on valid pixel"); - // psAssert (source->variance->data.F32[yPeak][xPeak] > 0, "peak must be on valid pixel"); - - for (psS32 row = 0; row < source->pixels->numRows ; row++) { - - psF32 yDiff = row - yPeak; - if (fabs(yDiff) > radius) continue; - - psF32 *vPix = source->pixels->data.F32[row]; - psF32 *vWgt = source->variance->data.F32[row]; - psImageMaskType *vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[row]; - - for (psS32 col = 0; col < source->pixels->numCols ; col++, vPix++, vWgt++) { - if (vMsk) { - if (*vMsk & maskVal) { - vMsk++; - continue; - } - vMsk++; - } - if (isnan(*vPix)) continue; - - psF32 xDiff = col - xPeak; - if (fabs(xDiff) > radius) continue; - - // radius is just a function of (xDiff, yDiff) - if (!VALID_RADIUS(xDiff, yDiff, R2)) continue; - - psF32 pDiff = *vPix - sky; - psF32 wDiff = *vWgt; - - // skip pixels below specified significance level. for a PSFs, this - // over-weights the wings of bright stars compared to those of faint stars. - // for the estimator used for extended source analysis (where the window - // function is allowed to be arbitrarily large), we need to clip to avoid - // negative second moments. - if (PS_SQR(pDiff) < minSN2*wDiff) continue; // - if ((minSN > 0.0) && (pDiff < 0)) continue; // - - // Apply a Gaussian window function. Be careful with the window function. S/N - // weighting over weights the sky for faint sources - if (sigma > 0.0) { - // XXX a lot of extra flops; can we pre-calculate? - psF32 z = (PS_SQR(xDiff) + PS_SQR(yDiff))*rsigma2; - assert (z >= 0.0); - psF32 weight = exp(-z); - - wDiff *= weight; - pDiff *= weight; - } - - Var += wDiff; - Sum += pDiff; - - psF32 xWght = xDiff * pDiff; - psF32 yWght = yDiff * pDiff; - - X1 += xWght; - Y1 += yWght; - - peakPixel = PS_MAX (*vPix, peakPixel); - numPixels++; - } - } - - // if we have less than (1/4) of the possible pixels (in circle or box), force a retry - int minPixels = PS_MIN(0.75*R2, source->pixels->numCols*source->pixels->numRows/4.0); - - // XXX EAM - the limit is a bit arbitrary. make it user defined? - if ((numPixels < minPixels) || (Sum <= 0)) { - psTrace ("psModules.objects", 3, "insufficient valid pixels (%d vs %d; %f) for source\n", numPixels, minPixels, Sum); - return (false); - } - - // calculate the first moment. - float Mx = X1/Sum; - float My = Y1/Sum; - if ((fabs(Mx) > radius) || (fabs(My) > radius)) { - psTrace ("psModules.objects", 3, "large centroid swing; invalid peak %d, %d\n", source->peak->x, source->peak->y); - return (false); - } - - 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); - - // add back offset of peak in primary image - // also offset from pixel index to pixel coordinate - // (the calculation above uses pixel index instead of coordinate) - // 0.5 PIX: moments are calculated using the pixel index and converted here to pixel coords - source->moments->Mx = Mx + xOff + 0.5; - source->moments->My = My + yOff + 0.5; - - source->moments->Sum = Sum; - source->moments->SN = Sum / sqrt(Var); - source->moments->Peak = peakPixel; - source->moments->nPixels = numPixels; + // do 2 passes : the first pass should use a somewhat smaller radius and no sigma window to + // get an unbiased (but probably noisy) centroid + if (!pmSourceMomentsGetCentroid (source, 0.75*radius, 0.0, minSN, maskVal, source->peak->xf, source->peak->yf)) { + return false; + } + // second pass applies the Gaussian window and uses the centroid from the first pass + if (!pmSourceMomentsGetCentroid (source, radius, sigma, minSN, maskVal, source->moments->Mx, source->moments->My)) { + return false; + } // Now calculate higher-order moments, using the above-calculated first moments to adjust coordinates @@ -234,8 +125,8 @@ Sum = 0.0; // the second pass may include slightly different pixels, re-determine Sum - // center of mass in subimage. Note: the calculation below uses pixel index, so we do not - // correct xCM, yCM to pixel coordinate here. - psF32 xCM = Mx + xPeak; // coord of peak in subimage - psF32 yCM = My + yPeak; // coord of peak in subimage + // center of mass in subimage. Note: the calculation below uses pixel index, so we correct + // xCM, yCM from pixel coords to pixel index here. + psF32 xCM = source->moments->Mx - 0.5 - source->pixels->col0; // coord of peak in subimage + psF32 yCM = source->moments->My - 0.5 - source->pixels->row0; // coord of peak in subimage for (psS32 row = 0; row < source->pixels->numRows ; row++) { @@ -263,6 +154,5 @@ // radius is just a function of (xDiff, yDiff) psF32 r2 = PS_SQR(xDiff) + PS_SQR(yDiff); - psF32 r = sqrt(r2); - if (r > radius) continue; + if (r2 > R2) continue; psF32 fDiff = *vPix - sky; @@ -274,13 +164,10 @@ // stars. if (PS_SQR(pDiff) < minSN2*wDiff) continue; - // if (pDiff < 0) continue; + if ((minSN > 0.0) && (pDiff < 0)) continue; // // Apply a Gaussian window function. Be careful with the window function. S/N // weighting over weights the sky for faint sources if (sigma > 0.0) { - // XXX a lot of extra flops; can we do pre-calculate? - // XXX we were re-calculating r2 (maybe the compiler caught this?) - // psF32 z = (PS_SQR(xDiff) + PS_SQR(yDiff))*rsigma2; - psF32 z = r2 * rsigma2; + psF32 z = r2 * rsigma2; assert (z >= 0.0); psF32 weight = exp(-z); @@ -292,16 +179,9 @@ Sum += pDiff; -# if (1) -# if (0) - if (fDiff < 0) continue; -# endif + // Kron Flux uses the 1st radial moment (NOT Gaussian windowed?) + psF32 r = sqrt(r2); psF32 rf = r * fDiff; psF32 rh = sqrt(r) * fDiff; psF32 rs = fDiff; -# else - psF32 rf = r * pDiff; - psF32 rh = sqrt(r) * pDiff; - psF32 rs = pDiff; -# endif psF32 x = xDiff * pDiff; @@ -363,5 +243,4 @@ // Calculate the Kron magnitude (make this block optional?) - // float radKron = 2.5*source->moments->Mrf; float radKinner = 1.0*source->moments->Mrf; float radKron = 2.5*source->moments->Mrf; @@ -397,5 +276,4 @@ // radKron is just a function of (xDiff, yDiff) psF32 r2 = PS_SQR(xDiff) + PS_SQR(yDiff); - psF32 r = sqrt(r2); psF32 pDiff = *vPix - sky; @@ -407,4 +285,5 @@ if (PS_SQR(pDiff) < minSN2*wDiff) continue; + psF32 r = sqrt(r2); if (r < radKron) { Sum += pDiff; @@ -434,6 +313,146 @@ 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", - 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); + 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); return(true); } + +bool pmSourceMomentsGetCentroid(pmSource *source, psF32 radius, psF32 sigma, psF32 minSN, psImageMaskType maskVal, float xGuess, float yGuess) { + + // First Pass: calculate the first moments (these are subtracted from the coordinates below) + // Sum = SUM (z - sky) + // X1 = SUM (x - xc)*(z - sky) + // .. etc + + psF32 sky = 0.0; + + psF32 peakPixel = -PS_MAX_F32; + psS32 numPixels = 0; + psF32 Sum = 0.0; + psF32 Var = 0.0; + psF32 X1 = 0.0; + psF32 Y1 = 0.0; + psF32 R2 = PS_SQR(radius); + psF32 minSN2 = PS_SQR(minSN); + psF32 rsigma2 = 0.5 / PS_SQR(sigma); + + float xPeak = xGuess - source->pixels->col0; // coord of peak in subimage + float yPeak = yGuess - source->pixels->row0; // coord of peak in subimage + + // we are guaranteed to have a valid pixel and variance at this location (right? right?) + // psF32 weightNorm = source->pixels->data.F32[yPeak][xPeak] / sqrt (source->variance->data.F32[yPeak][xPeak]); + // psAssert (isfinite(source->pixels->data.F32[yPeak][xPeak]), "peak must be on valid pixel"); + // psAssert (isfinite(source->variance->data.F32[yPeak][xPeak]), "peak must be on valid pixel"); + // psAssert (source->variance->data.F32[yPeak][xPeak] > 0, "peak must be on valid pixel"); + + // the moments [Sum(x*f) / Sum(f)] are calculated in pixel index values, and should + // not depend on the fractional pixel location of the source. However, the aperture + // (radius) and the Gaussian window (sigma) depend subtly on the fractional pixel + // position of the expected centroid + + for (psS32 row = 0; row < source->pixels->numRows ; row++) { + + psF32 yDiff = row + 0.5 - yPeak; + if (fabs(yDiff) > radius) continue; + + psF32 *vPix = source->pixels->data.F32[row]; + psF32 *vWgt = source->variance->data.F32[row]; + psImageMaskType *vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[row]; + + for (psS32 col = 0; col < source->pixels->numCols ; col++, vPix++, vWgt++) { + if (vMsk) { + if (*vMsk & maskVal) { + vMsk++; + continue; + } + vMsk++; + } + if (isnan(*vPix)) continue; + + psF32 xDiff = col + 0.5 - xPeak; + if (fabs(xDiff) > radius) continue; + + // radius is just a function of (xDiff, yDiff) + psF32 r2 = PS_SQR(xDiff) + PS_SQR(yDiff); + if (r2 > R2) continue; + + psF32 pDiff = *vPix - sky; + psF32 wDiff = *vWgt; + + // skip pixels below specified significance level. for a PSFs, this + // over-weights the wings of bright stars compared to those of faint stars. + // for the estimator used for extended source analysis (where the window + // function is allowed to be arbitrarily large), we need to clip to avoid + // negative second moments. + if (PS_SQR(pDiff) < minSN2*wDiff) continue; // + if ((minSN > 0.0) && (pDiff < 0)) continue; // + + // Apply a Gaussian window function. Be careful with the window function. S/N + // weighting over weights the sky for faint sources + if (sigma > 0.0) { + psF32 z = r2*rsigma2; + assert (z >= 0.0); + psF32 weight = exp(-z); + + wDiff *= weight; + pDiff *= weight; + } + + Var += wDiff; + Sum += pDiff; + + psF32 xWght = xDiff * pDiff; + psF32 yWght = yDiff * pDiff; + + X1 += xWght; + Y1 += yWght; + + peakPixel = PS_MAX (*vPix, peakPixel); + numPixels++; + } + } + + // if we have less than (1/4) of the possible pixels (in circle or box), force a retry + int minPixels = PS_MIN(0.75*R2, source->pixels->numCols*source->pixels->numRows/4.0); + + // XXX EAM - the limit is a bit arbitrary. make it user defined? + if ((numPixels < minPixels) || (Sum <= 0)) { + psTrace ("psModules.objects", 3, "insufficient valid pixels (%d vs %d; %f) for source\n", numPixels, minPixels, Sum); + return (false); + } + + // calculate the first moment. + float Mx = X1/Sum; + float My = Y1/Sum; + if ((fabs(Mx) > radius) || (fabs(My) > radius)) { + psTrace ("psModules.objects", 3, "large centroid swing; invalid peak %d, %d\n", source->peak->x, source->peak->y); + return (false); + } + + 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); + + // add back offset of peak in primary image + // also offset from pixel index to pixel coordinate + // (the calculation above uses pixel index instead of coordinate) + // 0.5 PIX: moments are calculated using the pixel index and converted here to pixel coords + + // we only update the centroid if the position is not supplied from elsewhere + bool skipCentroid = false; + skipCentroid |= (source->mode & PM_SOURCE_MODE_EXTERNAL); // skip externally supplied positions + skipCentroid |= (source->mode2 & PM_SOURCE_MODE2_MATCHED); // skip sources defined by other image positions + + if (skipCentroid) { + source->moments->Mx = source->peak->xf; + source->moments->My = source->peak->yf; + } else { + source->moments->Mx = Mx + xGuess; + source->moments->My = My + yGuess; + } + + source->moments->Sum = Sum; + source->moments->SN = Sum / sqrt(Var); + source->moments->Peak = peakPixel; + source->moments->nPixels = numPixels; + + return true; +} Index: trunk/psModules/src/objects/pmSourcePhotometry.c =================================================================== --- trunk/psModules/src/objects/pmSourcePhotometry.c (revision 29044) +++ trunk/psModules/src/objects/pmSourcePhotometry.c (revision 29546) @@ -62,18 +62,19 @@ /** - this function is used to calculate the three defined source magnitudes: - - apMag : only if S/N > AP_MIN_SN - : is optionally corrected for curve-of-growth if: - - the source is a STAR (PSF) - - the option is selected (mode & PM_SOURCE_PHOT_GROWTH) - - psfMag : all sources with non-NULL modelPSF - : is optionally corrected for aperture residual if: - - the source is a STAR (PSF) - - the option is selected (mode & PM_SOURCE_PHOT_APCORR) - - - extMag : all sources with non-NULL modelEXT + this function is used to calculate the three defined source magnitudes: + - apMag : only if S/N > AP_MIN_SN + : is optionally corrected for curve-of-growth if: + - the source is a STAR (PSF) + - the option is selected (mode & PM_SOURCE_PHOT_GROWTH) + - psfMag : all sources with non-NULL modelPSF + : is optionally corrected for aperture residual if: + - the source is a STAR (PSF) + - the option is selected (mode & PM_SOURCE_PHOT_APCORR) + + - extMag : all sources with non-NULL modelEXT **/ // XXX masked region should be (optionally) elliptical +// if mode is PM_SOURCE_PHOT_PSFONLY, we skip all other magnitudes bool pmSourceMagnitudes (pmSource *source, pmPSF *psf, pmSourcePhotometryMode mode, psImageMaskType maskVal, psImageMaskType markVal) { @@ -88,8 +89,11 @@ pmModel *model; - source->psfMag = NAN; - source->extMag = NAN; - source->errMag = NAN; - source->apMag = NAN; + source->psfMag = NAN; + source->extMag = NAN; + source->errMag = NAN; + source->apMag = NAN; + source->apMagRaw = NAN; + source->apFlux = NAN; + source->apFluxErr = NAN; // we must have a valid model @@ -114,4 +118,5 @@ // measure PSF model photometry // XXX TEST: do not use flux scale + // XXX NOTE: turn this back on? if (0 && psf->FluxScale) { // the source peak pixel is guaranteed to be on the image, and only minimally different from the source center @@ -127,4 +132,8 @@ } + if (mode == PM_SOURCE_PHOT_PSFONLY) { + return true; + } + // if we have a collection of model fits, check if one of them is a pointer to modelEXT if (source->modelFits) { @@ -148,4 +157,5 @@ // for PSFs, correct both apMag and psfMag to same system, consistent with infinite flux star in aperture RADIUS + // XXX add a flag for "ap_mag is corrected?" if ((mode & PM_SOURCE_PHOT_APCORR) && isPSF && psf && psf->ApTrend) { // the source peak pixel is guaranteed to be on the image, and only minimally different from the source center @@ -181,5 +191,7 @@ // if we are measuring aperture photometry and applying the growth correction, // we need to shift the flux in the selected pixels (but not the mask) - psImage *flux = NULL, *mask = NULL; // Star flux and mask images, to photometer + psImage *flux = NULL; + psImage *variance = NULL; + psImage *mask = NULL; // Star flux and mask images, to photometer if (mode & PM_SOURCE_PHOT_INTERP) { float dx = 0.5 - x + (int)x; @@ -198,9 +210,10 @@ } else { flux = source->pixels; + variance = source->variance; mask = source->maskObj; } // measure object aperture photometry - status = pmSourcePhotometryAper (&source->apMagRaw, model, flux, mask, maskVal); + status = pmSourcePhotometryAperSource (source, model, flux, variance, mask, maskVal); if (!status) { psTrace ("psModules.objects", 3, "fail mag : bad Ap Mag"); @@ -214,4 +227,5 @@ if (psf->growth && (mode & PM_SOURCE_PHOT_GROWTH)) { source->apMag = source->apMagRaw + pmGrowthCurveCorrect (psf->growth, source->apRadius); + // XXX correct the apFlux? } if (mode & PM_SOURCE_PHOT_APCORR) { @@ -233,7 +247,7 @@ /* -aprMag' - fitMag = flux*skySat + r^2*rflux*skyBias + ApTrend(x,y) -(aprMag - flux*skySat - r^2*rflux*skyBias) - fitMAg = ApTrend(x,y) -(aprMag - flux*skySat - r^2*rflux*skyBias) = fitMAg + ApTrend(x,y) + aprMag' - fitMag = flux*skySat + r^2*rflux*skyBias + ApTrend(x,y) + (aprMag - flux*skySat - r^2*rflux*skyBias) - fitMAg = ApTrend(x,y) + (aprMag - flux*skySat - r^2*rflux*skyBias) = fitMAg + ApTrend(x,y) */ @@ -267,5 +281,22 @@ // return source aperture magnitude -bool pmSourcePhotometryAper (float *apMag, pmModel *model, psImage *image, psImage *mask, psImageMaskType maskVal) +bool pmSourcePhotometryAperSource (pmSource *source, pmModel *model, psImage *image, psImage *variance, psImage *mask, psImageMaskType maskVal) +{ + PS_ASSERT_PTR_NON_NULL(source, false); + PS_ASSERT_PTR_NON_NULL(image, false); + PS_ASSERT_PTR_NON_NULL(mask, false); + + if (DO_SKY) { + PS_ASSERT_PTR_NON_NULL(model, false); + } + + bool status; + status = pmSourcePhotometryAper(&source->apMagRaw, &source->apFlux, &source->apFluxErr, model, image, variance, mask, maskVal); + + return status; +} + +// return source aperture magnitude +bool pmSourcePhotometryAper (float *apMag, float *apFluxOut, float *apFluxErr, pmModel *model, psImage *image, psImage *variance, psImage *mask, psImageMaskType maskVal) { PS_ASSERT_PTR_NON_NULL(apMag, false); @@ -277,34 +308,34 @@ } - float apSum = 0; float sky = 0; - *apMag = NAN; + float apFlux = 0; + float apFluxVar = 0; if (DO_SKY) { sky = model->params->data.F32[PM_PAR_SKY]; - } else { - sky = 0; } psF32 **imData = image->data.F32; psImageMaskType **mkData = mask->data.PS_TYPE_IMAGE_MASK_DATA; - int nAperPix = 0; - - // measure apMag + psF32 **varData = (variance) ? variance->data.F32 : image->data.F32; // if variance is not supplied, assume gain of 1.0, no read noise + + // measure apFlux and apFluxVar, save apMag if not NAN + // XXX note that these fluxes/mags are uncorrected for masked pixels + // XXX raise a bit if the aperture has a masked pixel (not marked)? for (int iy = 0; iy < image->numRows; iy++) { for (int ix = 0; ix < image->numCols; ix++) { - if (mkData[iy][ix] & maskVal) - continue; - apSum += imData[iy][ix] - sky; - nAperPix ++; - // fprintf (stderr, "aper: %d %d %f %f %f\n", ix, iy, sky, imData[iy][ix], apSum); - } - } - if (apSum <= 0) { + if (mkData[iy][ix] & maskVal) continue; + apFlux += imData[iy][ix] - sky; + apFluxVar += varData[iy][ix]; + } + } + if (apFluxOut) *apFluxOut = apFlux; + if (apFluxErr) *apFluxErr = sqrt(apFluxVar); + + if (apFlux <= 0) { *apMag = NAN; - return true; - } - - *apMag = -2.5*log10(apSum); + } else { + *apMag = -2.5*log10(apFlux); + } return true; } @@ -419,7 +450,13 @@ } + // NOTE: until 2010.10.01, these measurements included a 3sigma-per-pixel significance + // this followed what we understood as the definition given to us + // by Armin, but it always seemed a poor idea -- a faint source is unlikely to have any 3sigma pixels. + // changed to remove the per-pixel filter. + for (int iy = 0; iy < flux->numRows; iy++) { for (int ix = 0; ix < flux->numCols; ix++) { // only count up the stats in the unmarked region (ie, the aperture) + // skip the marked pixels; these are not relevant if (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & markVal) { continue; @@ -430,14 +467,12 @@ } - float SN = flux->data.F32[iy][ix] / sqrt(variance->data.F32[iy][ix]); - - if (SN > +FLUX_LIMIT) { + float value = flux->data.F32[iy][ix]; + + if (value > 0.0) { nGood ++; - fGood += fabs(flux->data.F32[iy][ix]); - } - - if (SN < -FLUX_LIMIT) { + fGood += fabs(value); + } else { nBad ++; - fBad += fabs(flux->data.F32[iy][ix]); + fBad += fabs(value); } } @@ -613,14 +648,14 @@ switch (term) { - case 0: + case 0: factor = 1; break; - case 1: + case 1: factor = xi + Pi->col0; break; - case 2: + case 2: factor = yi + Pi->row0; break; - default: + default: psAbort("invalid term for pmSourceWeight"); } @@ -691,14 +726,14 @@ switch (term) { - case 0: + case 0: factor = 1; break; - case 1: + case 1: factor = xi + Pi->col0; break; - case 2: + case 2: factor = yi + Pi->row0; break; - default: + default: psAbort("invalid term for pmSourceWeight"); } Index: trunk/psModules/src/objects/pmSourcePhotometry.h =================================================================== --- trunk/psModules/src/objects/pmSourcePhotometry.h (revision 29044) +++ trunk/psModules/src/objects/pmSourcePhotometry.h (revision 29546) @@ -35,4 +35,5 @@ PM_SOURCE_PHOT_INTERP = 0x0008, PM_SOURCE_PHOT_DIFFSTATS = 0x0010, + PM_SOURCE_PHOT_PSFONLY = 0x0020, } pmSourcePhotometryMode; @@ -44,7 +45,19 @@ bool pmSourcePhotometryAper( - float *apMag, ///< aperture flux magnitude + float *apMag, + float *apFluxOut, + float *apFluxErr, pmModel *model, ///< model used for photometry psImage *image, ///< image pixels to be used + psImage *variance, ///< variance pixels to be used + psImage *mask, ///< mask of pixels to ignore + psImageMaskType maskVal ///< Value to mask +); + +bool pmSourcePhotometryAperSource( + pmSource *source, ///< aperture flux magnitude + pmModel *model, ///< model used for photometry + psImage *image, ///< image pixels to be used + psImage *variance, ///< variance pixels to be used psImage *mask, ///< mask of pixels to ignore psImageMaskType maskVal ///< Value to mask