Index: /branches/pap/psphot/src/psphot.h =================================================================== --- /branches/pap/psphot/src/psphot.h (revision 25301) +++ /branches/pap/psphot/src/psphot.h (revision 25302) @@ -49,5 +49,5 @@ bool psphotPSFstats (pmReadout *readout, psMetadata *recipe, pmPSF *psf); bool psphotMomentsStats (pmReadout *readout, psMetadata *recipe, psArray *sources); -bool psphotFitSourcesLinear (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, bool final); +bool psphotFitSourcesLinear (pmReadout *readout, psArray *sources, const psMetadata *recipe, const pmPSF *psf, bool final); bool psphotReplaceUnfitSources (psArray *sources, psImageMaskType maskVal); @@ -64,5 +64,5 @@ bool psphotGuessModel_Threaded (psThreadJob *job); -bool psphotMagnitudes (pmConfig *config, pmReadout *readout, const pmFPAview *view, psArray *sources, pmPSF *psf); +bool psphotMagnitudes (pmConfig *config, pmReadout *readout, const pmFPAview *view, psArray *sources, const pmPSF *psf); bool psphotMagnitudes_Threaded (psThreadJob *job); @@ -76,5 +76,5 @@ bool psphotExtendedSourceAnalysis (pmReadout *readout, psArray *sources, psMetadata *recipe); bool psphotExtendedSourceFits (pmReadout *readout, psArray *sources, psMetadata *recipe); -bool psphotFake(pmReadout *readout, const pmPSF *psf, psMetadata *recipe, const psArray *realSources); +bool psphotFake(pmConfig *config, pmReadout *readout, const pmFPAview *view, const pmPSF *psf, psMetadata *recipe, const psArray *realSources); // thread-related: Index: /branches/pap/psphot/src/psphotFake.c =================================================================== --- /branches/pap/psphot/src/psphotFake.c (revision 25301) +++ /branches/pap/psphot/src/psphotFake.c (revision 25302) @@ -21,4 +21,6 @@ PM_ASSERT_READOUT_NON_NULL(ro, false); PM_ASSERT_READOUT_VARIANCE(ro, false); + assert(magLim); + assert(radius); psKernel *kernel = psImageSmoothKernel(smoothSigma, smoothNsigma); // Kernel used for smoothing @@ -52,24 +54,20 @@ psFree(normModel); - if (magLim) { - // The signal-to-noise of the smoothed image is: S/N ~ I_smooth / sqrt(variance * factor) - - // since the variance factor tells us the variance in the smoothed image. Now, the trick is working - // out what the intensity in the smoothed image is, and how it is related to the flux. We are - // convolving Io G_* with G_PSF/2pi.w^2, where G_* is the approximately-Gaussian PSF of the star, - // G_PSF is the pretty-close-matching Gaussian of the convolution kernel, Io is the peak flux in the - // unsmoothed image, and w is the width of the Gaussian. Now, a normalised 2D Gaussian convolved - // with itself has a normalisation of 1/2pi(w^2+w^2) = 1/4pi.w^2. Therefore: - // I_smooth = Flux / 2*norm. - float peakLim = thresh * sqrtf(meanVar * factor); // Limiting peak value in smoothed image - float fluxLim = 2.0 * norm * peakLim; // Limiting flux in original - *magLim = -2.5 * log10f(fluxLim); - psTrace("psphot.fake", 1, "Limiting peak: %f\n", peakLim); - psTrace("psphot.fake", 1, "Limiting flux: %f\n", fluxLim); - psTrace("psphot.fake", 1, "Limiting mag: %f\n", *magLim); - } - if (radius) { - *radius = smoothSigma * smoothNsigma; - } + // The signal-to-noise of the smoothed image is: S/N ~ I_smooth / sqrt(variance * factor) + // since the variance factor tells us the variance in the smoothed image. Now, the trick is working + // out what the intensity in the smoothed image is, and how it is related to the flux. We are + // convolving Io G_* with G_PSF/2pi.w^2, where G_* is the approximately-Gaussian PSF of the star, + // G_PSF is the pretty-close-matching Gaussian of the convolution kernel, Io is the peak flux in the + // unsmoothed image, and w is the width of the Gaussian. Now, a normalised 2D Gaussian convolved + // with itself has a normalisation of 1/2pi(w^2+w^2) = 1/4pi.w^2. Therefore: + // I_smooth = Flux / 2*norm. + float peakLim = thresh * sqrtf(meanVar * factor); // Limiting peak value in smoothed image + float fluxLim = 2.0 * norm * peakLim; // Limiting flux in original + *magLim = -2.5 * log10f(fluxLim); + psTrace("psphot.fake", 1, "Limiting peak: %f\n", peakLim); + psTrace("psphot.fake", 1, "Limiting flux: %f\n", fluxLim); + psTrace("psphot.fake", 1, "Limiting mag: %f\n", *magLim); + + *radius = smoothSigma * smoothNsigma; return true; @@ -90,4 +88,6 @@ PS_ASSERT_VECTOR_NON_NULL(magOffsets, false); PS_ASSERT_VECTOR_TYPE(magOffsets, PS_TYPE_F32, false); + assert(xSrc); + assert(ySrc); int numBins = magOffsets->n; // Number of bins @@ -140,5 +140,6 @@ // *** in this section, perform the photometry for fake sources *** -bool psphotFake(pmReadout *readout, const pmPSF *psf, psMetadata *recipe, const psArray *realSources) +bool psphotFake(pmConfig *config, pmReadout *readout, const pmFPAview *view, const pmPSF *psf, + psMetadata *recipe, const psArray *realSources) { PM_ASSERT_READOUT_NON_NULL(readout, false); @@ -201,4 +202,5 @@ } + // XXX Could speed this up significantly by only convolving the central pixels of each fake source psImage *significance = psphotSignificanceImage(readout, recipe, 2, maskVal); if (!significance) { @@ -213,7 +215,10 @@ int numBins = magOffsets->n; // Number of bins int numCols = readout->image->numCols, numRows = readout->image->numRows; // Size of image - psVector *frac = psVectorAlloc(numBins, PS_TYPE_F32); // Fraction of sources in each bin + psVector *count = psVectorAlloc(numBins, PS_TYPE_S32); // Number of sources found in each bin + psArray *fakeSources = psArrayAllocEmpty(numSources * numBins); // Fake sources + psVector *fakeEnd = psVectorAllocEmpty(numSources * numBins, PS_TYPE_S32); // End index of each bin + long numFound = 0; // Number of sources found for (int i = 0; i < numBins; i++) { - int numFound = 0; // Number found + int num = 0; // Number found for (int j = 0; j < numSources; j++) { // Coordinates of interest @@ -221,10 +226,25 @@ float y = PS_MIN(yFake->data.F32[i][j] + 0.5, numRows - 1); int xPix = x, yPix = y; // Pixel coordinates - if (significance->data.F32[yPix][xPix] > thresh) { + float sig = significance->data.F32[yPix][xPix]; // Significance of pixel + if (sig > thresh) { + pmSource *source = pmSourceAlloc(); // Fake source + source->peak = pmPeakAlloc(x, y, NAN, PM_PEAK_LONE); + if (!pmSourceDefinePixels(source, readout, x, y, radius)) { + psErrorClear(); + continue; + } + source->modelPSF = pmModelFromPSFforXY(psf, x, y, 2.0 * sqrtf(sig)); + source->type = PM_SOURCE_TYPE_STAR; + + num++; + + fakeSources->data[numFound] = source; numFound++; } } - frac->data.F32[i] = (float)numFound / (float)numSources; - } + count->data.S32[i] = num; + fakeEnd->data.S32[i] = numFound; + } + fakeSources->n = numFound; psFree(xFake); @@ -232,16 +252,71 @@ psFree(significance); + if (!psphotFitSourcesLinear(readout, fakeSources, recipe, psf, true)) { + psError(PS_ERR_UNKNOWN, false, "Unable to perform linear fit on fake sources."); + psFree(fakeSources); + psFree(fakeEnd); + psFree(count); + return false; + } + if (!psphotMagnitudes(config, readout, view, fakeSources, psf)) { + psError(PS_ERR_UNKNOWN, false, "Unable to measure magnitudes of fake sources."); + psFree(fakeSources); + psFree(fakeEnd); + psFree(count); + return false; + } + + psVector *magErrMean = psVectorAlloc(numBins, PS_TYPE_F32); // Mean error in magnitude for each bin + psVector *magErrStdev = psVectorAlloc(numBins, PS_TYPE_F32); // Stdev of error in magnitude for each bin + psVector *magErr = psVectorAlloc(numSources, PS_TYPE_F32); // Magnitude errors + psVector *magMask = psVectorAlloc(numSources, PS_TYPE_VECTOR_MASK); // Mask for magnitude errors + psStats *magStats = psStatsAlloc(PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV); // Statistics + for (int i = 0, index = 0; i < numBins; i++) { + psStatsInit(magStats); + psVectorInit(magMask, 0); + + float magRef = magLim + magOffsets->data.F32[i]; // Reference magnitude + for (int j = 0; index < fakeEnd->data.S32[i]; j++, index++) { + pmSource *source = fakeSources->data[index]; // Source of interest + if (!source || !isfinite(source->psfMag)) { + magMask->data.PS_TYPE_VECTOR_MASK_DATA[j] = 0xFF; + continue; + } + magErr->data.F32[j] = source->psfMag - magRef; + } + if (!psVectorStats(magStats, magErr, NULL, magMask, 0xFF)) { + // Probably because we don't have enough sources + psErrorClear(); + } + magErrMean->data.F32[i] = magStats->sampleMean; + magErrStdev->data.F32[i] = magStats->sampleStdev; + } + psFree(magStats); + psFree(magMask); + psFree(magErr); + + psFree(fakeSources); + psFree(fakeEnd); + // XXX How do we get the results out? psMetadata *stats = psMetadataAlloc(); - psMetadataAddVector(stats, PS_LIST_TAIL, "FAKE.EFF", PS_META_REPLACE, - "Efficiency fractions", frac); + psMetadataAddS32(stats, PS_LIST_TAIL, "FAKE.NUM", PS_META_REPLACE, + "Number of sources per bin", numSources); + psMetadataAddF32(stats, PS_LIST_TAIL, "FAKE.REF", PS_META_REPLACE, + "Efficiency reference magnitude", magLim); psMetadataAddVector(stats, PS_LIST_TAIL, "FAKE.MAG", PS_META_REPLACE, "Efficiency magnitudes", magOffsets); - psMetadataAddF32(stats, PS_LIST_TAIL, "FAKE.REF", PS_META_REPLACE, - "Efficiency reference magnitude", magLim); + psMetadataAddVector(stats, PS_LIST_TAIL, "FAKE.COUNTS", PS_META_REPLACE, + "Number of sources retrieved", count); + psMetadataAddVector(stats, PS_LIST_TAIL, "FAKE.ERR", PS_META_REPLACE, + "Mean magnitude differences", magErrMean); + psMetadataAddVector(stats, PS_LIST_TAIL, "FAKE.RMS", PS_META_REPLACE, + "Stdev in magnitude differences", magErrStdev); psMetadataConfigWrite(stats, "fake.stats"); psFree(stats); - psFree(frac); + psFree(count); + psFree(magErrMean); + psFree(magErrStdev); return true; Index: /branches/pap/psphot/src/psphotFitSourcesLinear.c =================================================================== --- /branches/pap/psphot/src/psphotFitSourcesLinear.c (revision 25301) +++ /branches/pap/psphot/src/psphotFitSourcesLinear.c (revision 25302) @@ -12,5 +12,5 @@ static bool SetBorderMatrixElements (psSparseBorder *border, pmReadout *readout, psArray *sources, bool constant_weights, int SKY_FIT_ORDER, psImageMaskType markVal); -bool psphotFitSourcesLinear (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, bool final) { +bool psphotFitSourcesLinear (pmReadout *readout, psArray *sources, const psMetadata *recipe, const pmPSF *psf, bool final) { bool status; Index: /branches/pap/psphot/src/psphotMagnitudes.c =================================================================== --- /branches/pap/psphot/src/psphotMagnitudes.c (revision 25301) +++ /branches/pap/psphot/src/psphotMagnitudes.c (revision 25302) @@ -1,5 +1,5 @@ # include "psphotInternal.h" -bool psphotMagnitudes(pmConfig *config, pmReadout *readout, const pmFPAview *view, psArray *sources, pmPSF *psf) { +bool psphotMagnitudes(pmConfig *config, pmReadout *readout, const pmFPAview *view, psArray *sources, const pmPSF *psf) { bool status = false; @@ -64,5 +64,5 @@ psArrayAdd(job->args, 1, cells->data[j]); // sources - psArrayAdd(job->args, 1, psf); + psArrayAdd(job->args, 1, (pmPSF*)psf); // Casting away const psArrayAdd(job->args, 1, binning); psArrayAdd(job->args, 1, backModel); Index: /branches/pap/psphot/src/psphotReadout.c =================================================================== --- /branches/pap/psphot/src/psphotReadout.c (revision 25301) +++ /branches/pap/psphot/src/psphotReadout.c (revision 25302) @@ -224,5 +224,5 @@ psphotMagnitudes(config, readout, view, sources, psf); - if (!psphotFake(readout, psf, recipe, sources)) { + if (!psphotFake(config, readout, view, psf, recipe, sources)) { psErrorStackPrint(stderr, "Unable to determine detection efficiencies from fake sources"); psErrorClear();