Changeset 27838 for branches/tap_branches/psphot/src/psphotEfficiency.c

Timestamp:

May 3, 2010, 8:41:49 AM (16 years ago)

Author:

eugene

Message:

updates from trunk

Location:

branches/tap_branches

Files:

• . (modified) (1 prop)
• psphot (modified) (1 prop)
• psphot/src (modified) (1 prop)
• psphot/src/psphotEfficiency.c (modified) (15 diffs, 1 prop)

branches/tap_branches/psphot/src

@@ -19 +19 @@
 psphotMomentsStudy
 psphotPetrosianStudy
+psphotForced
+psphotMakePSF
+psphotStack

@@ -19 +19 @@
 psphotMomentsStudy
 psphotPetrosianStudy
+psphotForced
+psphotMakePSF
+psphotStack

branches/tap_branches/psphot/src/psphotEfficiency.c

@@ -4 +4 @@
                     PM_MODEL_STATUS_BADARGS | PM_MODEL_STATUS_LIMITS) // Mask to apply to models
 
-//#define TESTING
-
+# define TESTING 0
 
 // Calculate the limiting magnitude for an image
@@ -17 +16 @@
                      float *covarFactor,// Covariance factor
                      const pmReadout *ro,     // Readout of interest
-                     const pmPSF *psf,        // Point-spread function
+                     pmPSF *psf,              // Point-spread function
                      float thresh,            // Threshold for source identification
                      float smoothSigma,       // Gaussian smoothing sigma
@@ -33 +32 @@
 
     psKernel *kernel = psImageSmoothKernel(smoothSigma, smoothNsigma); // Kernel used for smoothing
-    psKernel *newCovar = psImageCovarianceCalculate(kernel, ro->covariance); // Covariance matrix
+    *covarFactor = psImageCovarianceCalculateFactor(kernel, ro->covariance);
     psFree(kernel);
-    *covarFactor = psImageCovarianceFactor(newCovar); // Variance factor
-    psFree(newCovar);
 
     psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); // Statistics for variance
@@ -128 +125 @@
     psFree(rng);
 
+    bool oldThreads = pmReadoutFakeThreads(true); // Old threading status
+
     pmReadout *fakeRO = pmReadoutAlloc(NULL); // Fake readout
     if (!pmReadoutFakeFromVectors(fakeRO, numCols, numRows, xAll, yAll, magAll,
@@ -142 +141 @@
     psFree(yAll);
 
+    pmReadoutFakeThreads(oldThreads);
+
     psBinaryOp(ro->image, ro->image, "+", fakeRO->image);
     psFree(fakeRO);
@@ -148 +149 @@
 }
 
+bool psphotEfficiency (pmConfig *config, const pmFPAview *view)
+{
+    bool status = true;
+
+    // select the appropriate recipe information
+    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
+    psAssert (recipe, "missing recipe?");
+
+    int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM");
+    psAssert (status, "programming error: must define PSPHOT.INPUT.NUM");
+
+    // skip the chisq image (optionally?)
+    int chisqNum = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.CHISQ.NUM");
+    if (!status) chisqNum = -1;
+
+    // loop over the available readouts
+    for (int i = 0; i < num; i++) {
+        if (i == chisqNum) continue; // skip chisq image
+        if (!psphotEfficiencyReadout (config, view, "PSPHOT.INPUT", i, recipe)) {
+            psError (PSPHOT_ERR_CONFIG, false, "failed to measure detection efficiency for PSPHOT.INPUT entry %d", i);
+            return false;
+        }
+    }
+    return true;
+}
 
 // Determine detection efficiency
-bool psphotEfficiency(pmConfig *config, pmReadout *readout, const pmFPAview *view, const pmPSF *psf,
-                psMetadata *recipe, const psArray *realSources)
+bool psphotEfficiencyReadout(pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe)
 {
+    bool status = true;
+
+    psTimerStart("psphot.fake");
+
+    // find the currently selected readout
+    pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest
+    psAssert (file, "missing file?");
+
+    pmReadout *readout = pmFPAviewThisReadout(view, file->fpa);
+    psAssert (readout, "missing readout?");
+
+    pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
+    psAssert (detections, "missing detections?");
+
+    psArray *realSources = detections->allSources;
+    psAssert (realSources, "missing sources?");
+
+    // XXX do we need to skip this step if we do not have a psf?
+    pmPSF *psf = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF");
+    psAssert (psf, "missing psf?");
+
     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
+    float smoothNsigma = psMetadataLookupF32(&status, recipe, "PEAKS_SMOOTH_NSIGMA"); // Smoothing limit
+    psAssert (status && isfinite(smoothNsigma), "Unable to find PEAKS_SMOOTH_NSIGMA in recipe (or invalid value)");
+
+    float thresh = psMetadataLookupF32(&status, recipe, "PEAKS_NSIGMA_LIMIT_2");
+    psAssert (status && isfinite(thresh), "Unable to find PEAKS_NSIGMA_LIMIT_2 in recipe (or invalid value)");
+
+    psVector *magOffsets = psMetadataLookupVector(&status, recipe, "EFF.MAG"); // Magnitude offsets
+    psAssert (status, "Unable to find EFF.MAG F32 vector in recipe");
+    psAssert (magOffsets->type.type == PS_TYPE_F32, "Unable to find EFF.MAG F32 vector in recipe");
+
+    int numSources = psMetadataLookupS32(&status, recipe, "EFF.NUM"); // Number of sources for each bin
+    psAssert (status && (numSources > 0), "Unable to find EFF.NUM in recipe (or invalid value)");
+
+    float minGauss = psMetadataLookupF32(&status, recipe, "PEAKS_MIN_GAUSS"); // Minimum valid fraction of kernel
+    psAssert (status && isfinite(minGauss), "PEAKS_MIN_GAUSS is not set in recipe (or invalid)");
+
+    // find the PSF size information (why is this not part of the psf structure?)
+    float fwhmMajor = psMetadataLookupF32(NULL, readout->analysis, "FWHM_MAJ"); // PSF size in x
+    float fwhmMinor = psMetadataLookupF32(NULL, readout->analysis, "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 minGauss = psMetadataLookupF32(NULL, recipe, "PEAKS_MIN_GAUSS"); // Minimum valid fraction of kernel
-    if (!isfinite(minGauss)) {
-        psWarning("PEAKS_MIN_GAUSS is not set in recipe; using default value");
-        minGauss = 0.5;
+        psError(PSPHOT_ERR_CONFIG, false, "Unable to find FWHM_MAJ and FWHM_MIN in readout->analysis");
+        return false;
     }
 
@@ -202 +237 @@
     // remove all sources, adding noise for subtracted sources
     psphotRemoveAllSources(realSources, recipe);
-    //    psphotAddNoise(readout, realSources, recipe);
-
-
-#if defined(TESTING) && 0
+
+#if TESTING
     {
         psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS);
@@ -234 +267 @@
     }
 
-#ifdef TESTING
+#if TESTING
     psphotSaveImage(NULL, readout->image, "orig_image.fits");
     psphotSaveImage(NULL, readout->variance, "orig_variance.fits");
@@ -249 +282 @@
     }
 
-#ifdef TESTING
+#if TESTING
     psphotSaveImage(NULL, readout->image, "fake_image.fits");
     psphotSaveImage(NULL, readout->variance, "fake_variance.fits");
@@ -364 +397 @@
     psFree(significance);
 
-    if (!psphotFitSourcesLinear(readout, fakeSourcesAll, recipe, psf, true)) {
+    if (!psphotFitSourcesLinearReadout(recipe, readout, fakeSourcesAll, psf, true)) {
         psError(PS_ERR_UNKNOWN, false, "Unable to perform linear fit on fake sources.");
         psFree(fakeSources);
@@ -371 +404 @@
     }
 
-    // Disable aperture corrections; casting away const!
+    // Disable aperture corrections (save current value)
     pmTrend2D *apTrend = psf->ApTrend;  // Aperture trend
-    ((pmPSF*)psf)->ApTrend = NULL;
-
-    if (!psphotMagnitudes(config, readout, view, fakeSourcesAll, psf)) {
+    psf->ApTrend = NULL;
+
+    if (!psphotMagnitudesReadout(config, recipe, view, readout, 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!
+        psf->ApTrend = apTrend; // Casting away const!
         return false;
     }
     psFree(fakeSourcesAll);
 
-    // Re-enable aperture corrections; casting away const!
-    ((pmPSF*)psf)->ApTrend = apTrend;
+    // Replace aperture corrections
+    psf->ApTrend = apTrend;
 
     psVector *magDiffMean = psVectorAlloc(numBins, PS_TYPE_F32); // Mean difference in magnitude for each bin
@@ -399 +432 @@
         psVectorInit(magMask, 0);
 
-#ifdef TESTING
+#if TESTING
         psString name = NULL;
         psStringAppend(&name, "fake_%d.dat", i);
@@ -415 +448 @@
             }
 
-#ifdef TESTING
+#if 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],
@@ -451 +484 @@
         }
 
-#ifdef TESTING
+#if TESTING
         fclose(file);
 #endif
@@ -476 +509 @@
     psLogMsg("psphot", PS_LOG_INFO, "Detection efficiency: %lf sec\n", psTimerClear("psphot.fake"));
 
-
     return true;
 }

@@ -19 +19 @@
 psphotMomentsStudy
 psphotPetrosianStudy
+psphotForced
+psphotMakePSF
+psphotStack