Index: /trunk/psModules/src/camera/pmReadoutFake.c =================================================================== --- /trunk/psModules/src/camera/pmReadoutFake.c (revision 25229) +++ /trunk/psModules/src/camera/pmReadoutFake.c (revision 25230) @@ -50,13 +50,23 @@ } -bool pmReadoutFakeFromSources(pmReadout *readout, int numCols, int numRows, const psArray *sources, - const psVector *xOffset, const psVector *yOffset, const pmPSF *psf, - float minFlux, int radius, bool circularise, bool normalisePeak) + +bool pmReadoutFakeFromVectors(pmReadout *readout, int numCols, int numRows, + const psVector *x, const psVector *y, const psVector *mag, + const psVector *xOffset, const psVector *yOffset, + const pmPSF *psf, float minFlux, int radius, + bool circularise, bool normalisePeak) { PS_ASSERT_PTR_NON_NULL(readout, false); PS_ASSERT_INT_LARGER_THAN(numCols, 0, false); PS_ASSERT_INT_LARGER_THAN(numRows, 0, false); - PS_ASSERT_ARRAY_NON_NULL(sources, false); - + PS_ASSERT_VECTOR_NON_NULL(x, false); + PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F32, false); + PS_ASSERT_VECTOR_NON_NULL(y, false); + PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F32, false); + PS_ASSERT_VECTORS_SIZE_EQUAL(y, x, false); + PS_ASSERT_VECTOR_NON_NULL(mag, false); + PS_ASSERT_VECTOR_TYPE(mag, PS_TYPE_F32, false); + PS_ASSERT_VECTORS_SIZE_EQUAL(mag, x, false); + long numSources = x->n; // Number of sources if (xOffset || yOffset) { PS_ASSERT_VECTOR_NON_NULL(xOffset, false); @@ -64,9 +74,9 @@ PS_ASSERT_VECTORS_SIZE_EQUAL(xOffset, yOffset, false); PS_ASSERT_VECTOR_TYPE(xOffset, PS_TYPE_S32, false); - PS_ASSERT_VECTOR_TYPE_EQUAL(xOffset, yOffset, false); - if (xOffset->n != sources->n) { + PS_ASSERT_VECTOR_TYPE(yOffset, PS_TYPE_S32, false); + if (xOffset->n != numSources) { psError(PS_ERR_BAD_PARAMETER_SIZE, true, "Number of offset vectors (%ld) and sources (%ld) doesn't match", - xOffset->n, sources->n); + xOffset->n, numSources); return false; } @@ -81,26 +91,7 @@ psImageInit(readout->image, 0); - int numSources = sources->n; // Number of stars - for (int i = 0; i < numSources; i++) { - pmSource *source = sources->data[i]; // Source of interest - if (!source) { - continue; - } - if (source->mode & SOURCE_MASK) { - continue; - } - if (!isfinite(source->psfMag)) { - continue; - } - float x, y; // Coordinates of source - if (source->modelPSF) { - x = source->modelPSF->params->data.F32[PM_PAR_XPOS]; - y = source->modelPSF->params->data.F32[PM_PAR_YPOS]; - } else { - x = source->peak->xf; - y = source->peak->yf; - } - - float flux = powf(10.0, -0.4 * source->psfMag); // Flux of source + for (long i = 0; i < numSources; i++) { + float flux = powf(10.0, -0.4 * mag->data.F32[i]); // Flux of source + float xSrc = x->data.F32[i], ySrc = y->data.F32[i]; // Coordinates of source if (normalisePeak) { @@ -137,11 +128,11 @@ pmSource *fakeSource = pmSourceAlloc(); // Fake source to generate - fakeSource->peak = pmPeakAlloc(x, y, fakeModel->params->data.F32[PM_PAR_I0], PM_PEAK_LONE); + fakeSource->peak = pmPeakAlloc(xSrc, ySrc, fakeModel->params->data.F32[PM_PAR_I0], PM_PEAK_LONE); float fakeRadius = radius > 0 ? radius : PS_MAX(1.0, fakeModel->modelRadius(fakeModel->params, minFlux)); // Radius of fake source if (xOffset) { - if (!pmSourceDefinePixels(fakeSource, readout, x + xOffset->data.S32[i], - y + yOffset->data.S32[i], fakeRadius)) { + if (!pmSourceDefinePixels(fakeSource, readout, xSrc + xOffset->data.S32[i], + ySrc + yOffset->data.S32[i], fakeRadius)) { psErrorClear(); continue; @@ -153,5 +144,5 @@ } } else { - if (!pmSourceDefinePixels(fakeSource, readout, x, y, fakeRadius)) { + if (!pmSourceDefinePixels(fakeSource, readout, xSrc, ySrc, fakeRadius)) { psErrorClear(); continue; @@ -168,2 +159,53 @@ return true; } + + +bool pmReadoutFakeFromSources(pmReadout *readout, int numCols, int numRows, const psArray *sources, + const psVector *xOffset, const psVector *yOffset, const pmPSF *psf, + float minFlux, int radius, bool circularise, bool normalisePeak) +{ + PS_ASSERT_ARRAY_NON_NULL(sources, false); + + int numSources = sources->n; // Number of stars + psVector *x = psVectorAllocEmpty(numSources, PS_TYPE_F32); + psVector *y = psVectorAllocEmpty(numSources, PS_TYPE_F32); + psVector *mag = psVectorAllocEmpty(numSources, PS_TYPE_F32); + + int numGood = 0; // Number of good sources + for (int i = 0; i < numSources; i++) { + pmSource *source = sources->data[i]; // Source of interest + if (!source) { + continue; + } + if (source->mode & SOURCE_MASK) { + continue; + } + if (!isfinite(source->psfMag)) { + continue; + } + float xSrc, ySrc; // Coordinates of source + if (source->modelPSF) { + xSrc = source->modelPSF->params->data.F32[PM_PAR_XPOS]; + ySrc = source->modelPSF->params->data.F32[PM_PAR_YPOS]; + } else { + xSrc = source->peak->xf; + ySrc = source->peak->yf; + } + + x->data.F32[numGood] = xSrc; + y->data.F32[numGood] = ySrc; + mag->data.F32[numGood] = source->psfMag; + numGood++; + } + x->n = numGood; + y->n = numGood; + mag->n = numGood; + + bool status = pmReadoutFakeFromVectors(readout, numCols, numRows, x, y, mag, xOffset, yOffset, psf, + minFlux, radius, circularise, normalisePeak); + psFree(x); + psFree(y); + psFree(mag); + + return status; +} Index: /trunk/psphot/doc/efficiency.txt =================================================================== --- /trunk/psphot/doc/efficiency.txt (revision 25230) +++ /trunk/psphot/doc/efficiency.txt (revision 25230) @@ -0,0 +1,51 @@ +Strategy for efficiency analysis (a.k.a. "fake" stage) + +Want: +* Detection efficiency as a function of instrumental magnitude +* Masked fraction + +Given: +* Image +* Mask +* Variance +* PSF +* Sources on image +* Recipe: + + Magnitude bins (relative to guessed detection limit) for fake sources + + Number of fake sources for each bin + + + +Algorithm: +* Determine mean instrumental magnitude detection limit + + Have: + - Mean variance + - Smoothing size + - Covariance + - Threshold + + Calculate mean peak flux of source at threshold + + Integrate PSF to determine magnitude +* Remove all real sources from image +* Add fake sources into image +* Smooth image, variance (maybe only at positions of fake sources?) +* Count number of sources masked +* Count number of fake sources above significance level as a function of fake source magnitude + + +If the required density of fake sources is so high that they begin to +overlap, we will need to do add sources, smooth, and count +independently for each magnitude bin. This could be optimised +using the distributive property of convolutions (f*[g+h]=f*g+f*h) +if we smooth the image first, and then add in smoothed PSFs. + + +For a 5% statistical measurement of the detection efficiency per bin, we want +dN/N = 0.05 ==> N ~ 400 for Poisson statistics. +Assume a 80% masked fraction, so want 500 sources per bin. +For 5 magnitude bins, want 2500 sources +Assume worst case PSF is 15 pix FWHM --> ~ 51 pixels for 4-sigma either side +Total area is 6502500 pixels. +Chip images are 4846x4868 = 23590328 pixels = 3.6 times the worst total area. +Not too bad (especially since this is pretty much worst case, and we only +care about the central pixel, which is much smaller) +==> Let's just throw all sources on the same image. Index: /trunk/psphot/src/psphotFake.c =================================================================== --- /trunk/psphot/src/psphotFake.c (revision 25230) +++ /trunk/psphot/src/psphotFake.c (revision 25230) @@ -0,0 +1,225 @@ +# include "psphotInternal.h" + +#define TESTING + + +#ifdef TESTING +#define MIN_FLUX 0.1 // Minimum flux for faint sources +#endif + + +// 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 + float *minFlux, // Minimum flux, to return + const pmReadout *ro, // Readout of interest + float thresh, // Threshold for source identification + float xFWHM, float yFWHM, // Size of PSF + float smoothNsigma, // Smoothing limit + psImageMaskType maskVal // Value to mask + ); +{ + PM_ASSERT_READOUT_NON_NULL(ro, false); + PM_ASSERT_READOUT_VARIANCE(ro, false); + PS_ASSERT_METADATA_NON_NULL(recipe, false); + + float smoothSigma = 0.5*(FWHM_X + FWHM_Y) / (2.0*sqrtf(2.0*log(2.0))); + psKernel *kernel = psImageSmoothKernel(smoothSigma, smoothNsigma); // Kernel used for smoothing + psKernel *newCovar = psImageCovarianceCalculate(kernel, readout->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); + + if (magLim) { + *magLim = -2.5 * log10(thresh * sqrtf(meanVar * factor)); + } + if (minFlux) { + int fudge = psMetadataLookupF32(NULL, recipe, "FAKE.MINFLUX"); // Fudge factor for minimum flux + *minFlux = sqrtf(meanVar) * fudge; + } + + 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 psVector *magOffsets, // Magnitude offsets for fake sources + int numSources, // Number of fake sources for each bin + float refMag, // Reference magnitude + float minFlux // Minimum flux level for fake image + ) +{ + PM_ASSERT_READOUT_NON_NULL(ro, false); + PM_ASSERT_READOUT_IMAGE(ro, false); + PS_ASSERT_METADATA_NON_NULL(recipe, false); + PS_ASSERT_VECTOR_NON_NULL(magOffsets, false); + PS_ASSERT_VECTOR_TYPE(magOffsets, PS_TYPE_F32, false); + + int numBins = magOffset->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; i <= numBins; i++) { + psArray *data = psArrayAlloc(3); // Sources for bin + + psVector *x = psVectorAlloc(numSources, PS_TYPE_F32); + psVector *y = psVectorAlloc(numSources, PS_TYPE_F32); + + float mag = refMag + magOffset->data.F32[i]; // Instrumental magnitude of sources + + for (int j = 0; j <= numSources; j++) { + x->data.F32[j] = psRandomUniform(rng) * numCols; + y->data.F32[j] = psRandomUniform(rng) * numRows; + magAll->data.F32[i * numSources + j] = mag; + } + memcpy(xSrc->data.F32[i], x->data.F32, numSources * PSELEMTYPE_SIZEOF(PS_TYPE_F32)); + memcpy(ySrc->data.F32[i], y->data.F32, numSources * PSELEMTYPE_SIZEOF(PS_TYPE_F32)); + memcpy(&xAll->data.F32[i * numSources], x->data.F32, numSources * PSELEMTYPE_SIZEOF(PS_TYPE_F32)); + memcpy(&yAll->data.F32[i * numSources], y->data.F32, numSources * PSELEMTYPE_SIZEOF(PS_TYPE_F32)); + } + psFree(rng); + + pmReadout *fakeRO = pmReadoutAlloc(); // Fake readout + if (!pmReadoutFakeFromVectors(fakeRO, xAll, yAll, mag, NULL, NULL, psf, minFlux, NAN, false, false)) { + 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 pmPSF *psf, + const psMetadata *recipe, const psArray *realSources) +{ + PS_ASSERT_PTR_NON_NULL(config, false); + 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 xFWHM = psMetadataLookupF32(NULL, recipe, "FWHM_X"); // PSF size in x + float yFWHM = psMetadataLookupF32(NULL, recipe, "FWHM_Y"); // PSF size in y + if (!isfinite(xFWHM) || !isfinite(yFWHM)) { + psError(PSPHOT_ERR_CONFIG, false, "Unable to find FWHM_X and FWHM_Y 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 (!magOffset || magOffset->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; + } + psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); // Value to mask + + + // remove all sources, adding noise for subtracted sources + psphotRemoveAllSources(realSources, recipe); + psphotAddNoise(readout, realSources, recipe); + + float magLim; // Guess at limiting magnitude + float minFlux; // Minimum flux for fake image + if (!fakeLimit(&magLim, &minFlux, readout, thresh, xFWHM, yFWHM, smoothNsigma, maskVal)) { + psError(PS_ERR_UNKNOWN, false, "Unable to determine limits for image"); + return false; + } + + psImage *xFake = NULL, *yFake = NULL; // Coordinates of sources, each bin in a row + if (!fakeGenerate(&xFake, &yFake, readout, magOffsets, numSources, magLim, minFlux)) { + psError(PS_ERR_UNKNOWN, false, "Unable to generate fake sources"); + psFree(xFake); + psFree(yFake); + return false; + } + + 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; + } + + int numBins = magOffsets->n; // Number of bins + psVector *frac = psVectorAlloc(numBins, PS_TYPE_S32); // Fraction of sources in each bin + for (int i = 0; i < numBins; i++) { + int numFound = 0; // Number found + for (int j = 0; j < numSources; j++) { + float x = xFake->data.F32[i][j], y = yFake->data.F32[i][j]; // Coordinates of interest + int xPix = x, yPix = y; // Pixel coordinates + if (significance->data.F32[yPix][xPix] > thresh) { + numFound++; + } + } + frac->data.S32[i] = (float)numFound / (float)numSources; + } + + psFree(xFake); + psFree(yFake); + psFree(significance); + + // Putting results on recipe because that appears to be the psphot standard, but it's not a good idea + psMetadataAddVector(recipe, PS_LIST_TAIL, "FAKE.EFF", PS_META_REPLACE, + "Efficiency fractions", frac); + psMetadataAddVector(recipe, PS_LIST_TAIL, "FAKE.MAG", PS_META_REPLACE, + "Efficiency magnitudes", magOffsets); + psMetadataAddF32(recipe, PS_LIST_TAIL, "FAKE.REF", PS_META_REPLACE, + "Efficiency reference magnitude", magLim); + + psFree(frac); + + return true; +}