Index: /branches/pap/ippconfig/recipes/psphot.config =================================================================== --- /branches/pap/ippconfig/recipes/psphot.config (revision 25327) +++ /branches/pap/ippconfig/recipes/psphot.config (revision 25328) @@ -265,8 +265,7 @@ PSPHOT.CR.NSIGMA.SOFTEN F32 0.0025 # Softening parameter for weights -# Fake sources for detection efficiency -FAKE.NUM S32 500 # Number of fake sources per bin -@FAKE.MAG F32 -2.0 -1.0 -0.5 0.0 0.5 1.0 # Magnitude of fake sources relative to limit - +# Detection efficiency +EFF.NUM S32 500 # Number of fake sources per bin +@EFF.MAG F32 -2.0 -1.0 -0.5 0.0 0.5 1.0 # Magnitude of fake sources relative to limit # Recipe overrides for CHIP Index: /branches/pap/psphot/src/Makefile.am =================================================================== --- /branches/pap/psphot/src/Makefile.am (revision 25327) +++ /branches/pap/psphot/src/Makefile.am (revision 25328) @@ -129,5 +129,5 @@ psphotThreadTools.c \ psphotAddNoise.c \ - psphotFake.c + psphotEfficiency.c # dropped? psphotGrowthCurve.c Index: /branches/pap/psphot/src/psphot.h =================================================================== --- /branches/pap/psphot/src/psphot.h (revision 25327) +++ /branches/pap/psphot/src/psphot.h (revision 25328) @@ -76,5 +76,5 @@ bool psphotExtendedSourceAnalysis (pmReadout *readout, psArray *sources, psMetadata *recipe); bool psphotExtendedSourceFits (pmReadout *readout, psArray *sources, psMetadata *recipe); -bool psphotFake(pmConfig *config, pmReadout *readout, const pmFPAview *view, const pmPSF *psf, psMetadata *recipe, const psArray *realSources); +bool psphotEfficiency(pmConfig *config, pmReadout *readout, const pmFPAview *view, const pmPSF *psf, psMetadata *recipe, const psArray *realSources); // thread-related: Index: /branches/pap/psphot/src/psphotEfficiency.c =================================================================== --- /branches/pap/psphot/src/psphotEfficiency.c (revision 25328) +++ /branches/pap/psphot/src/psphotEfficiency.c (revision 25328) @@ -0,0 +1,438 @@ +# include "psphotInternal.h" + +#define MODEL_MASK (PM_MODEL_STATUS_NONCONVERGE | PM_MODEL_STATUS_OFFIMAGE | \ + PM_MODEL_STATUS_BADARGS | PM_MODEL_STATUS_LIMITS) // Mask to apply to models + +//#define TESTING + + +// Calculate the limiting magnitude for an image +// +// We limit ourselves to calculating the peak flux in the smoothed image, which should be close (modulo +// non-Gaussian PSF) to the limiting flux in the un-smoothed original image. +static bool effLimit(float *magLim, // Limiting magntiude, to return + int *radius, // Radius for fake sources, to return + float *minFlux, // Minimum flux for fake sources, to return + float *norm, // Normalisation of PSF (conversion: peak --> integrated flux) + const pmReadout *ro, // Readout of interest + const pmPSF *psf, // Point-spread function + float thresh, // Threshold for source identification + float smoothSigma, // Gaussian smoothing sigma + float smoothNsigma, // Smoothing limit + psImageMaskType maskVal // Value to mask + ) +{ + PM_ASSERT_READOUT_NON_NULL(ro, false); + PM_ASSERT_READOUT_VARIANCE(ro, false); + assert(magLim); + assert(radius); + assert(minFlux); + assert(norm); + + psKernel *kernel = psImageSmoothKernel(smoothSigma, smoothNsigma); // Kernel used for smoothing + psKernel *newCovar = psImageCovarianceCalculate(kernel, ro->covariance); // Covariance matrix + psFree(kernel); + float factor = psImageCovarianceFactor(newCovar); // Variance factor + psFree(newCovar); + + psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); // Statistics for variance + psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); // Random number generator + if (!psImageBackground(stats, NULL, ro->variance, ro->mask, maskVal, rng) || + !isfinite(stats->robustMedian)) { + psError(PSPHOT_ERR_DATA, false, "Unable to determine mean variance"); + psFree(stats); + psFree(rng); + return false; + } + psFree(rng); + float meanVar = stats->robustMedian; // Mean variance + psFree(stats); + + // Need to normalise out difference between Gaussian and real PSF + pmModel *normModel = pmModelFromPSFforXY(psf, 0.5 * ro->variance->numCols, + 0.5 * ro->variance->numRows, 1.0); // Model for normalisation + if (!normModel || (normModel->flags & MODEL_MASK)) { + psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to generate PSF model."); + psFree(normModel); + return false; + } + *norm = normModel->modelFlux(normModel->params); // Total flux for peak of 1.0 + psFree(normModel); + + // 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; + + *minFlux = 0.1 * sqrtf(meanVar); + + return true; +} + +/// Generate a fake image and add it in to the existing readout +static bool effGenerate(psImage **xSrc, psImage **ySrc, // Positions of sources + const pmReadout *ro, // Readout of interest + const pmPSF *psf, // Point-spread function + const psVector *magOffsets, // Magnitude offsets for fake sources + int numSources, // Number of fake sources for each bin + float refMag, // Reference magnitude + int radius, // Radius for fake sources + float minFlux // Minimum flux for fake sources + ) +{ + PM_ASSERT_READOUT_NON_NULL(ro, false); + PM_ASSERT_READOUT_IMAGE(ro, false); + 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 + int numCols = ro->image->numCols, numRows = ro->image->numRows; // Size of image + + *xSrc = psImageRecycle(*xSrc, numSources, numBins, PS_TYPE_F32); + *ySrc = psImageRecycle(*ySrc, numSources, numBins, PS_TYPE_F32); + + // Master list, for image creation + psVector *xAll = psVectorAlloc(numBins * numSources, PS_TYPE_F32); + psVector *yAll = psVectorAlloc(numBins * numSources, PS_TYPE_F32); + psVector *magAll = psVectorAlloc(numBins * numSources, PS_TYPE_F32); + + psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); // Random number generator + for (int i = 0, index = 0; i < numBins; i++) { + float mag = refMag + magOffsets->data.F32[i]; // Instrumental magnitude of sources + + for (int j = 0; j < numSources; j++, index++) { + xAll->data.F32[index] = psRandomUniform(rng) * numCols; + yAll->data.F32[index] = psRandomUniform(rng) * numRows; + magAll->data.F32[index] = mag; + } + memcpy((*xSrc)->data.F32[i], &xAll->data.F32[i * numSources], + numSources * PSELEMTYPE_SIZEOF(PS_TYPE_F32)); + memcpy((*ySrc)->data.F32[i], &yAll->data.F32[i * numSources], + numSources * PSELEMTYPE_SIZEOF(PS_TYPE_F32)); + } + psFree(rng); + + pmReadout *fakeRO = pmReadoutAlloc(NULL); // Fake readout + if (!pmReadoutFakeFromVectors(fakeRO, numCols, numRows, xAll, yAll, magAll, + NULL, NULL, psf, minFlux, radius, false, true)) { + psError(PS_ERR_UNKNOWN, false, "Unable to generate fake image"); + psFree(fakeRO); + psFree(xAll); + psFree(yAll); + psFree(magAll); + return false; + } + psFree(xAll); + psFree(yAll); + psFree(magAll); + + psBinaryOp(ro->image, ro->image, "+", fakeRO->image); + psFree(fakeRO); + + return true; +} + + +// Determine detection efficiency +bool psphotEfficiency(pmConfig *config, pmReadout *readout, const pmFPAview *view, const pmPSF *psf, + psMetadata *recipe, const psArray *realSources) +{ + PM_ASSERT_READOUT_NON_NULL(readout, false); + PM_ASSERT_READOUT_IMAGE(readout, false); + PS_ASSERT_PTR_NON_NULL(psf, false); + PS_ASSERT_METADATA_NON_NULL(recipe, false); + PS_ASSERT_ARRAY_NON_NULL(realSources, false); + + psTimerStart("psphot.fake"); + + // Collect recipe information + float fwhmMajor = psMetadataLookupF32(NULL, recipe, "FWHM_MAJ"); // PSF size in x + float fwhmMinor = psMetadataLookupF32(NULL, recipe, "FWHM_MIN"); // PSF size in y + if (!isfinite(fwhmMajor) || !isfinite(fwhmMinor) || fwhmMajor == 0.0 || fwhmMinor == 0.0) { + psError(PSPHOT_ERR_CONFIG, false, "Unable to find FWHM_MAJ and FWHM_MIN in recipe"); + return false; + } + float smoothNsigma = psMetadataLookupF32(NULL, recipe, "PEAKS_SMOOTH_NSIGMA"); // Smoothing limit + if (!isfinite(smoothNsigma)) { + psError(PSPHOT_ERR_CONFIG, false, "Unable to find PEAKS_SMOOTH_NSIGMA in recipe"); + return false; + } + float thresh = psMetadataLookupF32(NULL, recipe, "PEAKS_NSIGMA_LIMIT_2"); + if (!isfinite(thresh)) { + psError(PSPHOT_ERR_CONFIG, false, "Unable to find PEAKS_NSIGMA_LIMIT_2 in recipe"); + return false; + } + psVector *magOffsets = psMetadataLookupVector(NULL, recipe, "EFF.MAG"); // Magnitude offsets + if (!magOffsets || magOffsets->type.type != PS_TYPE_F32) { + psError(PSPHOT_ERR_CONFIG, false, "Unable to find EFF.MAG F32 vector in recipe"); + return NULL; + } + int numSources = psMetadataLookupS32(NULL, recipe, "EFF.NUM"); // Number of sources for each bin + if (numSources == 0) { + psError(PSPHOT_ERR_CONFIG, false, "Unable to find EFF.NUM in recipe"); + return NULL; + } + + float smoothSigma = 0.5*(fwhmMajor + fwhmMinor) / (2.0*sqrtf(2.0*log(2.0))); // Gaussian smoothing sigma + int numCols = readout->image->numCols, numRows = readout->image->numRows; // Size of image + + psImageMaskType maskVal = psMetadataLookupImageMask(NULL, recipe, "MASK.PSPHOT"); // Value to mask + + // remove all sources, adding noise for subtracted sources + psphotRemoveAllSources(realSources, recipe); + // psphotAddNoise(readout, realSources, recipe); + + +#if defined(TESTING) && 0 + { + psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); + for (int y = 0; y < numRows; y++) { + for (int x = 0; x < numCols; x++) { + readout->image->data.F32[y][x] = psRandomGaussian(rng); + readout->variance->data.F32[y][x] = 1.0; + readout->mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = 0; + } + } + psFree(readout->covariance); + readout->covariance = NULL; + psFree(rng); + } +#endif + + + + float magLim; // Guess at limiting magnitude + int radius; // Radius for fake sources + float minFlux; // Minimum flux for fake sources + float norm; // Normalisation of PSF + if (!effLimit(&magLim, &radius, &minFlux, &norm, readout, + psf, thresh, smoothSigma, smoothNsigma, maskVal)) { + psError(PS_ERR_UNKNOWN, false, "Unable to determine limits for image"); + return false; + } + +#ifdef TESTING + psphotSaveImage(NULL, readout->image, "orig_image.fits"); + psphotSaveImage(NULL, readout->variance, "orig_variance.fits"); + psphotSaveImage(NULL, readout->mask, "orig_mask.fits"); +#endif + + psImage *xFake = NULL, *yFake = NULL; // Coordinates of sources, each bin in a row + if (!effGenerate(&xFake, &yFake, readout, psf, magOffsets, numSources, magLim, radius, minFlux)) { + psError(PS_ERR_UNKNOWN, false, "Unable to generate fake sources"); + psFree(xFake); + psFree(yFake); + return false; + } + +#ifdef TESTING + psphotSaveImage(NULL, readout->image, "fake_image.fits"); + psphotSaveImage(NULL, readout->variance, "fake_variance.fits"); + psphotSaveImage(NULL, readout->mask, "fake_mask.fits"); +#endif + + // 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) { + psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to generate significance image"); + psFree(xFake); + psFree(yFake); + return false; + } + +#ifdef TESTING + psphotSaveImage(NULL, significance, "fake_sig.fits"); +#endif + + thresh *= thresh; // Significance image is actually significance-squared + + int numBins = magOffsets->n; // Number of bins + psVector *count = psVectorAlloc(numBins, PS_TYPE_S32); // Number of sources found in each bin + psArray *fakeSources = psArrayAlloc(numSources); // Fake sources in each bin + psArray *fakeSourcesAll = psArrayAllocEmpty(numSources * numBins); // All fake sources + for (int i = 0; i < numBins; i++) { + int numFound = 0; // Number found + + // Determine extraction size + float mag = magLim + magOffsets->data.F32[i]; // Magnitude for bin + float peak = powf(10.0, -0.4 * mag) / norm; // Peak flux + pmModel *model = pmModelFromPSFforXY(psf, numCols / 2.0, numRows / 2.0, peak); // Model for source + if (!model || (model->flags & MODEL_MASK)) { + psError(PS_ERR_UNKNOWN, true, "Unable to generate model for bin %d", i); + psFree(model); + return false; + } + float sourceRadius = PS_MAX(radius, model->modelRadius(model->params, minFlux)); // Radius for source + psFree(model); + + psArray *sources = psArrayAllocEmpty(numSources); // Sources in this bin + for (int j = 0; j < numSources; j++) { + // Coordinates of interest + float x = xFake->data.F32[i][j]; + float y = yFake->data.F32[i][j]; + int xPix = PS_MIN(x + 0.5, numCols - 1), yPix = PS_MIN(y + 0.5, numRows - 1); // Pixel coordinates + float sig = significance->data.F32[yPix][xPix]; // Significance of pixel + if (sig > thresh) { + pmSource *source = pmSourceAlloc(); // Fake source + source->peak = pmPeakAlloc(x, y, sig, PM_PEAK_LONE); + if (!pmSourceDefinePixels(source, readout, x, y, sourceRadius)) { + psErrorClear(); + continue; + } + source->peak->xf = x; + source->peak->yf = y; + source->modelPSF = pmModelFromPSFforXY(psf, x, y, 2.0 * sqrtf(sig)); + source->type = PM_SOURCE_TYPE_STAR; + + numFound++; + psArrayAdd(sources, sources->n, source); + psArrayAdd(fakeSourcesAll, fakeSourcesAll->n, source); + psFree(source); + } + } + fakeSources->data[i] = sources; + count->data.S32[i] = numFound; + } + psFree(xFake); + psFree(yFake); + psFree(significance); + + if (!psphotFitSourcesLinear(readout, fakeSourcesAll, recipe, psf, true)) { + psError(PS_ERR_UNKNOWN, false, "Unable to perform linear fit on fake sources."); + psFree(fakeSources); + psFree(count); + return false; + } + + // Disable aperture corrections; casting away const! + pmTrend2D *apTrend = psf->ApTrend; // Aperture trend + ((pmPSF*)psf)->ApTrend = NULL; + + if (!psphotMagnitudes(config, readout, view, fakeSourcesAll, psf)) { + psError(PS_ERR_UNKNOWN, false, "Unable to measure magnitudes of fake sources."); + psFree(fakeSources); + psFree(count); + ((pmPSF*)psf)->ApTrend = apTrend; // Casting away const! + return false; + } + psFree(fakeSourcesAll); + + // Re-enable aperture corrections; casting away const! + ((pmPSF*)psf)->ApTrend = apTrend; + + psVector *magDiffMean = psVectorAlloc(numBins, PS_TYPE_F32); // Mean difference in magnitude for each bin + psVector *magDiffStdev = psVectorAlloc(numBins, PS_TYPE_F32); // Stdev of diff in magnitude for each bin + psVector *magErrMean = psVectorAlloc(numBins, PS_TYPE_F32); // Mean error in magnitude for each bin + psVector *magDiff = psVectorAlloc(numSources, PS_TYPE_F32); // Magnitude differences + 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 | + PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); // Statistics + for (int i = 0; i < numBins; i++) { + psStatsInit(magStats); + psVectorInit(magMask, 0); + +#ifdef TESTING + psString name = NULL; + psStringAppend(&name, "fake_%d.dat", i); + FILE *file = fopen(name, "w"); + psFree(name); +#endif + + float magRef = magLim + magOffsets->data.F32[i]; // Reference magnitude + psArray *sources = fakeSources->data[i]; // Sources in bin + for (int j = 0; j < sources->n; j++) { + pmSource *source = sources->data[j]; // Source of interest + if (!source || !isfinite(source->psfMag)) { + magMask->data.PS_TYPE_VECTOR_MASK_DATA[j] = 0xFF; + continue; + } + +#ifdef TESTING + fprintf(file, "%f %f %f %f %f %f %f\n", source->peak->xf, source->peak->yf, + source->modelPSF->params->data.F32[PM_PAR_XPOS], + source->modelPSF->params->data.F32[PM_PAR_YPOS], + magRef, source->psfMag, source->errMag); +#endif + magDiff->data.F32[j] = source->psfMag - magRef; + magErr->data.F32[j] = source->errMag; + } + magDiff->n = sources->n; + magMask->n = sources->n; + magErr->n = sources->n; + if (!psVectorStats(magStats, magDiff, NULL, magMask, 0xFF)) { + // Probably because we don't have enough sources + psErrorClear(); + } + + if (isfinite(magStats->robustMedian) && isfinite(magStats->robustStdev)) { + magDiffMean->data.F32[i] = magStats->robustMedian; + magDiffStdev->data.F32[i] = magStats->robustStdev; + } else { + magDiffMean->data.F32[i] = magStats->sampleMean; + magDiffStdev->data.F32[i] = magStats->sampleStdev; + } + + if (!psVectorStats(magStats, magErr, NULL, magMask, 0xFF)) { + // Probably because we don't have enough sources + psErrorClear(); + } + + if (isfinite(magStats->robustMedian)) { + magErrMean->data.F32[i] = magStats->robustMedian; + } else { + magErrMean->data.F32[i] = magStats->sampleMean; + } + +#ifdef TESTING + fclose(file); +#endif + } + + psFree(magStats); + psFree(magDiff); + psFree(magMask); + psFree(magErr); + + psFree(fakeSources); + + // XXX How do we get the results out? + psMetadata *stats = psMetadataAlloc(); + 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); + psMetadataAddVector(stats, PS_LIST_TAIL, "FAKE.COUNTS", PS_META_REPLACE, + "Number of sources retrieved", count); + psMetadataAddVector(stats, PS_LIST_TAIL, "FAKE.DIFF.MEAN", PS_META_REPLACE, + "Mean magnitude differences", magDiffMean); + psMetadataAddVector(stats, PS_LIST_TAIL, "FAKE.DIFF.STDEV", PS_META_REPLACE, + "Stdev of magnitude differences", magDiffStdev); + psMetadataAddVector(stats, PS_LIST_TAIL, "FAKE.ERR.MEAN", PS_META_REPLACE, + "Mean error in magnitude differences", magErrMean); + psMetadataConfigWrite(stats, "fake.stats"); + psFree(stats); + + psFree(count); + psFree(magDiffMean); + psFree(magDiffStdev); + psFree(magErrMean); + + return true; +} Index: anches/pap/psphot/src/psphotFake.c =================================================================== --- /branches/pap/psphot/src/psphotFake.c (revision 25327) +++ (revision ) @@ -1,438 +1,0 @@ -# include "psphotInternal.h" - -#define MODEL_MASK (PM_MODEL_STATUS_NONCONVERGE | PM_MODEL_STATUS_OFFIMAGE | \ - PM_MODEL_STATUS_BADARGS | PM_MODEL_STATUS_LIMITS) // Mask to apply to models - -//#define TESTING - - -// Calculate the limiting magnitude for an image -// -// We limit ourselves to calculating the peak flux in the smoothed image, which should be close (modulo -// non-Gaussian PSF) to the limiting flux in the un-smoothed original image. -static bool fakeLimit(float *magLim, // Limiting magntiude, to return - int *radius, // Radius for fake sources, to return - float *minFlux, // Minimum flux for fake sources, to return - float *norm, // Normalisation of PSF (conversion: peak --> integrated flux) - const pmReadout *ro, // Readout of interest - const pmPSF *psf, // Point-spread function - float thresh, // Threshold for source identification - float smoothSigma, // Gaussian smoothing sigma - float smoothNsigma, // Smoothing limit - psImageMaskType maskVal // Value to mask - ) -{ - PM_ASSERT_READOUT_NON_NULL(ro, false); - PM_ASSERT_READOUT_VARIANCE(ro, false); - assert(magLim); - assert(radius); - assert(minFlux); - assert(norm); - - psKernel *kernel = psImageSmoothKernel(smoothSigma, smoothNsigma); // Kernel used for smoothing - psKernel *newCovar = psImageCovarianceCalculate(kernel, ro->covariance); // Covariance matrix - psFree(kernel); - float factor = psImageCovarianceFactor(newCovar); // Variance factor - psFree(newCovar); - - psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); // Statistics for variance - psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); // Random number generator - if (!psImageBackground(stats, NULL, ro->variance, ro->mask, maskVal, rng) || - !isfinite(stats->robustMedian)) { - psError(PSPHOT_ERR_DATA, false, "Unable to determine mean variance"); - psFree(stats); - psFree(rng); - return false; - } - psFree(rng); - float meanVar = stats->robustMedian; // Mean variance - psFree(stats); - - // Need to normalise out difference between Gaussian and real PSF - pmModel *normModel = pmModelFromPSFforXY(psf, 0.5 * ro->variance->numCols, - 0.5 * ro->variance->numRows, 1.0); // Model for normalisation - if (!normModel || (normModel->flags & MODEL_MASK)) { - psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to generate PSF model."); - psFree(normModel); - return false; - } - *norm = normModel->modelFlux(normModel->params); // Total flux for peak of 1.0 - psFree(normModel); - - // 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; - - *minFlux = 0.1 * sqrtf(meanVar); - - return true; -} - -/// Generate a fake image and add it in to the existing readout -static bool fakeGenerate(psImage **xSrc, psImage **ySrc, // Positions of sources - const pmReadout *ro, // Readout of interest - const pmPSF *psf, // Point-spread function - const psVector *magOffsets, // Magnitude offsets for fake sources - int numSources, // Number of fake sources for each bin - float refMag, // Reference magnitude - int radius, // Radius for fake sources - float minFlux // Minimum flux for fake sources - ) -{ - PM_ASSERT_READOUT_NON_NULL(ro, false); - PM_ASSERT_READOUT_IMAGE(ro, false); - 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 - int numCols = ro->image->numCols, numRows = ro->image->numRows; // Size of image - - *xSrc = psImageRecycle(*xSrc, numSources, numBins, PS_TYPE_F32); - *ySrc = psImageRecycle(*ySrc, numSources, numBins, PS_TYPE_F32); - - // Master list, for image creation - psVector *xAll = psVectorAlloc(numBins * numSources, PS_TYPE_F32); - psVector *yAll = psVectorAlloc(numBins * numSources, PS_TYPE_F32); - psVector *magAll = psVectorAlloc(numBins * numSources, PS_TYPE_F32); - - psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); // Random number generator - for (int i = 0, index = 0; i < numBins; i++) { - float mag = refMag + magOffsets->data.F32[i]; // Instrumental magnitude of sources - - for (int j = 0; j < numSources; j++, index++) { - xAll->data.F32[index] = psRandomUniform(rng) * numCols; - yAll->data.F32[index] = psRandomUniform(rng) * numRows; - magAll->data.F32[index] = mag; - } - memcpy((*xSrc)->data.F32[i], &xAll->data.F32[i * numSources], - numSources * PSELEMTYPE_SIZEOF(PS_TYPE_F32)); - memcpy((*ySrc)->data.F32[i], &yAll->data.F32[i * numSources], - numSources * PSELEMTYPE_SIZEOF(PS_TYPE_F32)); - } - psFree(rng); - - pmReadout *fakeRO = pmReadoutAlloc(NULL); // Fake readout - if (!pmReadoutFakeFromVectors(fakeRO, numCols, numRows, xAll, yAll, magAll, - NULL, NULL, psf, minFlux, radius, false, true)) { - psError(PS_ERR_UNKNOWN, false, "Unable to generate fake image"); - psFree(fakeRO); - psFree(xAll); - psFree(yAll); - psFree(magAll); - return false; - } - psFree(xAll); - psFree(yAll); - psFree(magAll); - - psBinaryOp(ro->image, ro->image, "+", fakeRO->image); - psFree(fakeRO); - - return true; -} - - -// *** in this section, perform the photometry for fake sources *** -bool psphotFake(pmConfig *config, pmReadout *readout, const pmFPAview *view, const pmPSF *psf, - psMetadata *recipe, const psArray *realSources) -{ - PM_ASSERT_READOUT_NON_NULL(readout, false); - PM_ASSERT_READOUT_IMAGE(readout, false); - PS_ASSERT_PTR_NON_NULL(psf, false); - PS_ASSERT_METADATA_NON_NULL(recipe, false); - PS_ASSERT_ARRAY_NON_NULL(realSources, false); - - psTimerStart("psphot.fake"); - - // Collect recipe information - float fwhmMajor = psMetadataLookupF32(NULL, recipe, "FWHM_MAJ"); // PSF size in x - float fwhmMinor = psMetadataLookupF32(NULL, recipe, "FWHM_MIN"); // PSF size in y - if (!isfinite(fwhmMajor) || !isfinite(fwhmMinor) || fwhmMajor == 0.0 || fwhmMinor == 0.0) { - psError(PSPHOT_ERR_CONFIG, false, "Unable to find FWHM_MAJ and FWHM_MIN in recipe"); - return false; - } - float smoothNsigma = psMetadataLookupF32(NULL, recipe, "PEAKS_SMOOTH_NSIGMA"); // Smoothing limit - if (!isfinite(smoothNsigma)) { - psError(PSPHOT_ERR_CONFIG, false, "Unable to find PEAKS_SMOOTH_NSIGMA in recipe"); - return false; - } - float thresh = psMetadataLookupF32(NULL, recipe, "PEAKS_NSIGMA_LIMIT_2"); - if (!isfinite(thresh)) { - psError(PSPHOT_ERR_CONFIG, false, "Unable to find PEAKS_NSIGMA_LIMIT_2 in recipe"); - return false; - } - psVector *magOffsets = psMetadataLookupVector(NULL, recipe, "FAKE.MAG"); // Magnitude offsets - if (!magOffsets || magOffsets->type.type != PS_TYPE_F32) { - psError(PSPHOT_ERR_CONFIG, false, "Unable to find FAKE.MAG F32 vector in recipe"); - return NULL; - } - int numSources = psMetadataLookupS32(NULL, recipe, "FAKE.NUM"); // Number of sources for each bin - if (numSources == 0) { - psError(PSPHOT_ERR_CONFIG, false, "Unable to find FAKE.NUM in recipe"); - return NULL; - } - - float smoothSigma = 0.5*(fwhmMajor + fwhmMinor) / (2.0*sqrtf(2.0*log(2.0))); // Gaussian smoothing sigma - int numCols = readout->image->numCols, numRows = readout->image->numRows; // Size of image - - psImageMaskType maskVal = psMetadataLookupImageMask(NULL, recipe, "MASK.PSPHOT"); // Value to mask - - // remove all sources, adding noise for subtracted sources - psphotRemoveAllSources(realSources, recipe); - // psphotAddNoise(readout, realSources, recipe); - - -#if defined(TESTING) && 0 - { - psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); - for (int y = 0; y < numRows; y++) { - for (int x = 0; x < numCols; x++) { - readout->image->data.F32[y][x] = psRandomGaussian(rng); - readout->variance->data.F32[y][x] = 1.0; - readout->mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = 0; - } - } - psFree(readout->covariance); - readout->covariance = NULL; - psFree(rng); - } -#endif - - - - float magLim; // Guess at limiting magnitude - int radius; // Radius for fake sources - float minFlux; // Minimum flux for fake sources - float norm; // Normalisation of PSF - if (!fakeLimit(&magLim, &radius, &minFlux, &norm, readout, - psf, thresh, smoothSigma, smoothNsigma, maskVal)) { - psError(PS_ERR_UNKNOWN, false, "Unable to determine limits for image"); - return false; - } - -#ifdef TESTING - psphotSaveImage(NULL, readout->image, "orig_image.fits"); - psphotSaveImage(NULL, readout->variance, "orig_variance.fits"); - psphotSaveImage(NULL, readout->mask, "orig_mask.fits"); -#endif - - psImage *xFake = NULL, *yFake = NULL; // Coordinates of sources, each bin in a row - if (!fakeGenerate(&xFake, &yFake, readout, psf, magOffsets, numSources, magLim, radius, minFlux)) { - psError(PS_ERR_UNKNOWN, false, "Unable to generate fake sources"); - psFree(xFake); - psFree(yFake); - return false; - } - -#ifdef TESTING - psphotSaveImage(NULL, readout->image, "fake_image.fits"); - psphotSaveImage(NULL, readout->variance, "fake_variance.fits"); - psphotSaveImage(NULL, readout->mask, "fake_mask.fits"); -#endif - - // 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) { - psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to generate significance image"); - psFree(xFake); - psFree(yFake); - return false; - } - -#ifdef TESTING - psphotSaveImage(NULL, significance, "fake_sig.fits"); -#endif - - thresh *= thresh; // Significance image is actually significance-squared - - int numBins = magOffsets->n; // Number of bins - psVector *count = psVectorAlloc(numBins, PS_TYPE_S32); // Number of sources found in each bin - psArray *fakeSources = psArrayAlloc(numSources); // Fake sources in each bin - psArray *fakeSourcesAll = psArrayAllocEmpty(numSources * numBins); // All fake sources - for (int i = 0; i < numBins; i++) { - int numFound = 0; // Number found - - // Determine extraction size - float mag = magLim + magOffsets->data.F32[i]; // Magnitude for bin - float peak = powf(10.0, -0.4 * mag) / norm; // Peak flux - pmModel *model = pmModelFromPSFforXY(psf, numCols / 2.0, numRows / 2.0, peak); // Model for source - if (!model || (model->flags & MODEL_MASK)) { - psError(PS_ERR_UNKNOWN, true, "Unable to generate model for bin %d", i); - psFree(model); - return false; - } - float sourceRadius = PS_MAX(radius, model->modelRadius(model->params, minFlux)); // Radius for source - psFree(model); - - psArray *sources = psArrayAllocEmpty(numSources); // Sources in this bin - for (int j = 0; j < numSources; j++) { - // Coordinates of interest - float x = xFake->data.F32[i][j]; - float y = yFake->data.F32[i][j]; - int xPix = PS_MIN(x + 0.5, numCols - 1), yPix = PS_MIN(y + 0.5, numRows - 1); // Pixel coordinates - float sig = significance->data.F32[yPix][xPix]; // Significance of pixel - if (sig > thresh) { - pmSource *source = pmSourceAlloc(); // Fake source - source->peak = pmPeakAlloc(x, y, sig, PM_PEAK_LONE); - if (!pmSourceDefinePixels(source, readout, x, y, sourceRadius)) { - psErrorClear(); - continue; - } - source->peak->xf = x; - source->peak->yf = y; - source->modelPSF = pmModelFromPSFforXY(psf, x, y, 2.0 * sqrtf(sig)); - source->type = PM_SOURCE_TYPE_STAR; - - numFound++; - psArrayAdd(sources, sources->n, source); - psArrayAdd(fakeSourcesAll, fakeSourcesAll->n, source); - psFree(source); - } - } - fakeSources->data[i] = sources; - count->data.S32[i] = numFound; - } - psFree(xFake); - psFree(yFake); - psFree(significance); - - if (!psphotFitSourcesLinear(readout, fakeSourcesAll, recipe, psf, true)) { - psError(PS_ERR_UNKNOWN, false, "Unable to perform linear fit on fake sources."); - psFree(fakeSources); - psFree(count); - return false; - } - - // Disable aperture corrections; casting away const! - pmTrend2D *apTrend = psf->ApTrend; // Aperture trend - ((pmPSF*)psf)->ApTrend = NULL; - - if (!psphotMagnitudes(config, readout, view, fakeSourcesAll, psf)) { - psError(PS_ERR_UNKNOWN, false, "Unable to measure magnitudes of fake sources."); - psFree(fakeSources); - psFree(count); - ((pmPSF*)psf)->ApTrend = apTrend; // Casting away const! - return false; - } - psFree(fakeSourcesAll); - - // Re-enable aperture corrections; casting away const! - ((pmPSF*)psf)->ApTrend = apTrend; - - psVector *magDiffMean = psVectorAlloc(numBins, PS_TYPE_F32); // Mean difference in magnitude for each bin - psVector *magDiffStdev = psVectorAlloc(numBins, PS_TYPE_F32); // Stdev of diff in magnitude for each bin - psVector *magErrMean = psVectorAlloc(numBins, PS_TYPE_F32); // Mean error in magnitude for each bin - psVector *magDiff = psVectorAlloc(numSources, PS_TYPE_F32); // Magnitude differences - 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 | - PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); // Statistics - for (int i = 0; i < numBins; i++) { - psStatsInit(magStats); - psVectorInit(magMask, 0); - -#ifdef TESTING - psString name = NULL; - psStringAppend(&name, "fake_%d.dat", i); - FILE *file = fopen(name, "w"); - psFree(name); -#endif - - float magRef = magLim + magOffsets->data.F32[i]; // Reference magnitude - psArray *sources = fakeSources->data[i]; // Sources in bin - for (int j = 0; j < sources->n; j++) { - pmSource *source = sources->data[j]; // Source of interest - if (!source || !isfinite(source->psfMag)) { - magMask->data.PS_TYPE_VECTOR_MASK_DATA[j] = 0xFF; - continue; - } - -#ifdef TESTING - fprintf(file, "%f %f %f %f %f %f %f\n", source->peak->xf, source->peak->yf, - source->modelPSF->params->data.F32[PM_PAR_XPOS], - source->modelPSF->params->data.F32[PM_PAR_YPOS], - magRef, source->psfMag, source->errMag); -#endif - magDiff->data.F32[j] = source->psfMag - magRef; - magErr->data.F32[j] = source->errMag; - } - magDiff->n = sources->n; - magMask->n = sources->n; - magErr->n = sources->n; - if (!psVectorStats(magStats, magDiff, NULL, magMask, 0xFF)) { - // Probably because we don't have enough sources - psErrorClear(); - } - - if (isfinite(magStats->robustMedian) && isfinite(magStats->robustStdev)) { - magDiffMean->data.F32[i] = magStats->robustMedian; - magDiffStdev->data.F32[i] = magStats->robustStdev; - } else { - magDiffMean->data.F32[i] = magStats->sampleMean; - magDiffStdev->data.F32[i] = magStats->sampleStdev; - } - - if (!psVectorStats(magStats, magErr, NULL, magMask, 0xFF)) { - // Probably because we don't have enough sources - psErrorClear(); - } - - if (isfinite(magStats->robustMedian)) { - magErrMean->data.F32[i] = magStats->robustMedian; - } else { - magErrMean->data.F32[i] = magStats->sampleMean; - } - -#ifdef TESTING - fclose(file); -#endif - } - - psFree(magStats); - psFree(magDiff); - psFree(magMask); - psFree(magErr); - - psFree(fakeSources); - - // XXX How do we get the results out? - psMetadata *stats = psMetadataAlloc(); - 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); - psMetadataAddVector(stats, PS_LIST_TAIL, "FAKE.COUNTS", PS_META_REPLACE, - "Number of sources retrieved", count); - psMetadataAddVector(stats, PS_LIST_TAIL, "FAKE.DIFF.MEAN", PS_META_REPLACE, - "Mean magnitude differences", magDiffMean); - psMetadataAddVector(stats, PS_LIST_TAIL, "FAKE.DIFF.STDEV", PS_META_REPLACE, - "Stdev of magnitude differences", magDiffStdev); - psMetadataAddVector(stats, PS_LIST_TAIL, "FAKE.ERR.MEAN", PS_META_REPLACE, - "Mean error in magnitude differences", magErrMean); - psMetadataConfigWrite(stats, "fake.stats"); - psFree(stats); - - psFree(count); - psFree(magDiffMean); - psFree(magDiffStdev); - psFree(magErrMean); - - return true; -} Index: /branches/pap/psphot/src/psphotReadout.c =================================================================== --- /branches/pap/psphot/src/psphotReadout.c (revision 25327) +++ /branches/pap/psphot/src/psphotReadout.c (revision 25328) @@ -224,5 +224,5 @@ psphotMagnitudes(config, readout, view, sources, psf); - if (!psphotFake(config, readout, view, psf, recipe, sources)) { + if (!psphotEfficiency(config, readout, view, psf, recipe, sources)) { psErrorStackPrint(stderr, "Unable to determine detection efficiencies from fake sources"); psErrorClear();