Index: /branches/pap/psphot/src/psphotFake.c =================================================================== --- /branches/pap/psphot/src/psphotFake.c (revision 25227) +++ /branches/pap/psphot/src/psphotFake.c (revision 25228) @@ -13,26 +13,18 @@ // 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 float fakeLimit(const pmReadout *ro, // Readout of interest - const psMetadata *recipe // psphot recipe - ); +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, NAN); - PM_ASSERT_READOUT_VARIANCE(ro, NAN); - PS_ASSERT_METADATA_NON_NULL(recipe, NAN); - - float xFWHM = psMetadataLookupF32(NULL, recipe, "FWHM_X"); - float yFWHM = psMetadataLookupF32(NULL, recipe, "FWHM_Y"); - if (!isfinite(xFWHM) || !isfinite(yFWHM)) { - psError(PSPHOT_ERR_CONFIG, false, "Unable to find FWHM_X and FWHM_Y"); - return NAN; - } + 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))); - float smoothNsigma = psMetadataLookupF32(NULL, recipe, "PEAKS_SMOOTH_NSIGMA"); - if (!isfinite(smoothSigma) || !isfinite(smoothNsigma)) { - psError(PSPHOT_ERR_CONFIG, false, "Unable to set smoothing parameters"); - return NAN; - } - psKernel *kernel = psImageSmoothKernel(smoothSigma, smoothNsigma); // Kernel used for smoothing psKernel *newCovar = psImageCovarianceCalculate(kernel, readout->covariance); // Covariance matrix @@ -40,6 +32,4 @@ float factor = psImageCovarianceFactor(newCovar); // Variance factor psFree(newCovar); - - psImageMaskType maskVal = psMetadataLookupImageMask(NULL, recipe, "MASK.PSPHOT"); // Value to mask psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); // Statistics for variance @@ -50,5 +40,5 @@ psFree(stats); psFree(rng); - return NAN; + return false; } psFree(rng); @@ -56,40 +46,32 @@ psFree(stats); - float thresh = psMetadataLookupF32(NULL, recipe, "PEAKS_NSIGMA_LIMIT_2"); - if (!isfinite(thresh)) { - psError(PSPHOT_ERR_CONFIG, false, "Unable to find threshold"); - return NAN; - } - - return -2.5 * log10(thresh * sqrtf(meanVar * factor)); + 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 psMetadata *recipe, // psphot recipe - float refMag // Reference magnitude + 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, NULL); - PM_ASSERT_READOUT_IMAGE(ro, NULL); - PS_ASSERT_METADATA_NON_NULL(recipe, NULL); - - psMetadata *fake = psMetadataLookupMetadata(NULL, recipe, "FAKE"); // Fake information - if (!fake) { - psError(PSPHOT_ERR_CONFIG, false, "Unable to find FAKE information in recipe"); - return NULL; - } - - 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; - } + 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 numSources = psMetadataLookupS32(NULL, recipe, "FAKE.NUM"); // Number of sources for each bin - int numCols = ro->image->numCols, numRows = ro->image->numRows; // Size of image - - *xSrc = psImageRecycle(*xSrc, numSources, numBins, PS_TYPE_F32); @@ -123,6 +105,21 @@ pmReadout *fakeRO = pmReadoutAlloc(); // Fake readout - pmReadoutFakeFromVectors(fakeRO, xAll, yAll, mag, - + 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; +} @@ -138,93 +135,90 @@ PS_ASSERT_ARRAY_NON_NULL(realSources, false); - psTimerStart("psphotFake"); - - // Generate array of fake sources and define the pixels with pmSourceDefinePixels - psArray *fakeOrig = fakeGenerate(config, readout, recipe); - - psArray *fakeMeasured = psArrayAlloc(fakeOrig->n); // Measured values for the sources in fakeOrig - for (int i = 0; i < fakeOrig->n; i++) { - fakeMeasured->data[i] = pmSourceCopy(fakeOrig->data[i]); - } - - int numCols = readout->image->numCols, numRows = readout->image->numRows; // Size of image - - pmReadout *fakeRO = pmReadoutAlloc(NULL); // Readout with fake sources - - if (!pmReadoutFakeFromSources(fakeRO, numCols, numRows, fakeOrig, NULL, NULL, - psf, MIN_FLUX, 0, false, false)) { - psError(psErrorCodeLast(), false, "Unable to generate fake readout"); - psFree(fakeRO); - psFree(fakeMeasured); - psFree(fakeOrig); - return false; - } - - -bool pmReadoutFakeFromSources(pmReadout *readout, ///< Output readout, or NULL - int numCols, int numRows, ///< Dimension of image - const psArray *sources, ///< Array of pmSource - const psVector *xOffset, ///< x offsets for sources (source -> img), or NULL - const psVector *yOffset, ///< y offsets for sources (source -> img), or NULL - const pmPSF *psf, ///< PSF for sources - float minFlux, ///< Minimum flux to bother about; for setting source radius - int radius, ///< Fixed radius for sources - bool circularise, ///< Circularise PSF model? - bool normalise ///< Normalise the peak value? - ); - - - - - // remove all sources + 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); - - // add noise for real (subtracted) objects - psphotAddNoise (readout, realSources, recipe); - - // XXX fake sources should measure peak->x,y, force sources should not - psImageMaskType maskVal = 0xff; - psImage *significance = psphotSignificanceImage(readout, recipe, 1, maskVal); - ppSimDetections (significance, recipe, fakeSources); - psFree (significance); - - // remove noise for subtracted objects (ie, return to normal noise level) - psphotSubNoise (readout, realSources, recipe); - - // replace all sources - psphotReplaceAllSources (realSources, recipe); - - // construct an initial model for each object - psphotGuessModels (config, readout, realSources, psf); - psphotGuessModels (config, readout, fakeSources, psf); - - // linear fit to real + fake sources - psArray *sources = ppSimMergeSources (realSources, fakeSources); - psphotFitSourcesLinear (readout, sources, recipe, psf, FALSE); // XXX option to NOT subtract - psphotReplaceAllSources (sources, recipe); - psFree (sources); // only frees the merged references - - // calculate source magnitudes (for which set??) - psphotMagnitudes(config, readout, view, fakeSources, psf); - - // drop the references to the image pixels held by each source - psphotSourceFreePixels (realSources); - psphotSourceFreePixels (fakeSources); - - // create the exported-metadata and free local data - // XXX this places the sources on readout->analysis as PSPHOT.SOURCES. modify? - // (or don't supply the sources, and do this with a different function) - psphotReadoutCleanup(config, readout, recipe, NULL, psf, NULL); - - pmCell *fakeCell = pmFPAfileThisCell (config->files, view, "PPSIM.FAKE.SOURCES"); psAssert (fakeCell, "no cell?"); - pmChip *fakeChip = pmFPAfileThisChip (config->files, view, "PPSIM.FAKE.SOURCES"); psAssert (fakeChip, "no chip?"); - pmReadout *fakeReadout = pmFPAfileThisReadout (config->files, view, "PPSIM.FAKE.SOURCES"); - if (!fakeReadout) { - fakeReadout = pmReadoutAlloc (fakeCell); - psFree (fakeReadout); // there is a copy on 'cell' as well - } - psAssert (fakeReadout, "no fakeReadout?"); - pmChipSetDataStatus (fakeChip, true); - psMetadataAddArray (fakeReadout->analysis, PS_LIST_TAIL, "PSPHOT.SOURCES", PS_META_REPLACE, "fake photometry ", fakeSources); + 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;