Changeset 25228

Timestamp:

Sep 1, 2009, 4:30:19 PM (17 years ago)

Author:

Paul Price

Message:

First cut at fake insertion and recovery.

File:

• branches/pap/psphot/src/psphotFake.c (modified) (6 diffs)

branches/pap/psphot/src/psphotFake.c

@@ -13 +13 @@
 // 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 +32 @@
     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 +40 @@
         psFree(stats);
         psFree(rng);
-        return NAN;
+        return false;
     }
     psFree(rng);
@@ -56 +46 @@
     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 +105 @@
 
     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 +135 @@
     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;