Changeset 25312

Timestamp:

Sep 9, 2009, 5:50:03 PM (17 years ago)

Author:

Paul Price

Message:

Got the magnitude measurements working properly now.

File:

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

branches/pap/psphot/src/psphotFake.c

@@ -4 +4 @@
                     PM_MODEL_STATUS_BADARGS | PM_MODEL_STATUS_LIMITS) // Mask to apply to models
 
-#define TESTING
+//#define TESTING
 
 
@@ -12 +12 @@
 // 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
+                      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
@@ -25 +27 @@
     assert(magLim);
     assert(radius);
+    assert(minFlux);
+    assert(norm);
 
     psKernel *kernel = psImageSmoothKernel(smoothSigma, smoothNsigma); // Kernel used for smoothing
@@ -53 +57 @@
         return false;
     }
-    float norm = normModel->modelFlux(normModel->params); // Total flux for peak of 1.0
+    *norm = normModel->modelFlux(normModel->params); // Total flux for peak of 1.0
     psFree(normModel);
 
@@ -65 +69 @@
     // 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
+    float fluxLim = 2.0 * *norm * peakLim; // Limiting flux in original
     *magLim = -2.5 * log10f(fluxLim);
     psTrace("psphot.fake", 1, "Limiting peak: %f\n", peakLim);
@@ -72 +76 @@
 
     *radius = smoothSigma * smoothNsigma;
+
+    *minFlux = 0.1 * sqrtf(meanVar);
 
     return true;
@@ -83 +89 @@
                          int numSources,                 // Number of fake sources for each bin
                          float refMag,                   // Reference magnitude
-                         int radius                      // Radius for fake sources
+                         int radius,                     // Radius for fake sources
+                         float minFlux                   // Minimum flux for fake sources
                          )
 {
@@ -122 +129 @@
     pmReadout *fakeRO = pmReadoutAlloc(NULL); // Fake readout
     if (!pmReadoutFakeFromVectors(fakeRO, numCols, numRows, xAll, yAll, magAll,
-                                  NULL, NULL, psf, NAN, radius, false, true)) {
+                                  NULL, NULL, psf, minFlux, radius, false, true)) {
         psError(PS_ERR_UNKNOWN, false, "Unable to generate fake image");
         psFree(fakeRO);
@@ -191 +198 @@
 
 
-
+#if defined(TESTING) && 0
     {
         psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS);
@@ -205 +212 @@
         psFree(rng);
     }
-
+#endif
 
 
@@ -211 +218 @@
     float magLim;                       // Guess at limiting magnitude
     int radius;                         // Radius for fake sources
-    if (!fakeLimit(&magLim, &radius, readout, psf, thresh, smoothSigma, smoothNsigma, maskVal)) {
+    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;
@@ -223 +233 @@
 
     psImage *xFake = NULL, *yFake = NULL; // Coordinates of sources, each bin in a row
-    if (!fakeGenerate(&xFake, &yFake, readout, psf, magOffsets, numSources, magLim, radius)) {
+    if (!fakeGenerate(&xFake, &yFake, readout, psf, magOffsets, numSources, magLim, radius, minFlux)) {
         psError(PS_ERR_UNKNOWN, false, "Unable to generate fake sources");
         psFree(xFake);
@@ -246 +256 @@
 
 #ifdef TESTING
-    psphotSaveImage(NULL, readout->mask, "fake_sig.fits");
+    psphotSaveImage(NULL, significance, "fake_sig.fits");
 #endif
 
@@ -253 +263 @@
     int numBins = magOffsets->n;                          // Number of bins
     psVector *count = psVectorAlloc(numBins, PS_TYPE_S32); // Number of sources found in each bin
-    psArray *fakeSources = psArrayAllocEmpty(numSources * numBins); // Fake sources
-    psVector *fakeEnd = psVectorAllocEmpty(numSources * numBins, PS_TYPE_S32); // End index of each bin
-    long numFound = 0;                                                           // Number of sources found
+    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 num = 0;               // Number found
-
+        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 = PS_MIN(xFake->data.F32[i][j] + 0.5, numCols - 1);
-            float y = PS_MIN(yFake->data.F32[i][j] + 0.5, numRows - 1);
-            int xPix = x, yPix = y;     // Pixel coordinates
+            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, NAN, PM_PEAK_LONE);
-                if (!pmSourceDefinePixels(source, readout, x, y, radius)) {
+                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;
 
-                num++;
-
-                fakeSources->data[numFound] = source;
                 numFound++;
+                psArrayAdd(sources, sources->n, source);
+                psArrayAdd(fakeSourcesAll, fakeSourcesAll->n, source);
+                psFree(source);
             }
         }
-        count->data.S32[i] = num;
-        fakeEnd->data.S32[i] = numFound;
-    }
-    fakeSources->n = numFound;
-
+        fakeSources->data[i] = sources;
+        count->data.S32[i] = numFound;
+    }
     psFree(xFake);
     psFree(yFake);
     psFree(significance);
 
-    if (!psphotFitSourcesLinear(readout, fakeSources, recipe, psf, true)) {
+    if (!psphotFitSourcesLinear(readout, fakeSourcesAll, recipe, psf, true)) {
         psError(PS_ERR_UNKNOWN, false, "Unable to perform linear fit on fake sources.");
         psFree(fakeSources);
-        psFree(fakeEnd);
         psFree(count);
         return false;
     }
-    if (!psphotMagnitudes(config, readout, view, fakeSources, psf)) {
+
+    // 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(fakeEnd);
         psFree(count);
-        return false;
-    }
-
+        ((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 *magErrStdev = psVectorAlloc(numBins, PS_TYPE_F32); // Stdev of error in magnitude for each bin
-    psVector *magErrErrMean = 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 *magErrErr = psVectorAlloc(numSources, PS_TYPE_F32); // Error in magnitude error
     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
-    psVector *robustMedian = psVectorAlloc(numBins, PS_TYPE_F32);
-    psVector *robustStdev = psVectorAlloc(numBins, PS_TYPE_F32);
-    for (int i = 0, index = 0; i < numBins; i++) {
+    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
-        int j;                          // Index
-        for (j = 0; index < fakeEnd->data.S32[i]; j++, index++) {
-            pmSource *source = fakeSources->data[index]; // Source of interest
+        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;
@@ -331 +363 @@
             }
 
-            fprintf(file, "%f %f %f %f %f\n", source->peak->xf, source->peak->yf,
+#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);
-            magErr->data.F32[j] = source->psfMag - magRef;
-            magErrErr->data.F32[j] = source->errMag;
-        }
-        magErr->n = j;
-        magMask->n = j;
+#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();
         }
-        magErrMean->data.F32[i] = magStats->sampleMean;
-        magErrStdev->data.F32[i] = magStats->sampleStdev;
-        robustMedian->data.F32[i] = magStats->robustMedian;
-        robustStdev->data.F32[i] = magStats->robustStdev;
-        if (!psVectorStats(magStats, magErrErr, NULL, magMask, 0xFF)) {
-            // Probably because we don't have enough sources
-            psErrorClear();
-        }
-        magErrErrMean->data.F32[i] = magStats->sampleMean;
+
+        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);
-    psFree(fakeEnd);
 
     // XXX How do we get the results out?
@@ -372 +421 @@
     psMetadataAddVector(stats, PS_LIST_TAIL, "FAKE.COUNTS", PS_META_REPLACE,
                         "Number of sources retrieved", count);
-    psMetadataAddVector(stats, PS_LIST_TAIL, "FAKE.ERR", PS_META_REPLACE,
-                        "Mean magnitude differences", magErrMean);
-    psMetadataAddVector(stats, PS_LIST_TAIL, "FAKE.RMS", PS_META_REPLACE,
-                        "Stdev in magnitude differences", magErrStdev);
-    psMetadataAddVector(stats, PS_LIST_TAIL, "FAKE.RM", PS_META_REPLACE,
-                        "Robust median magnitude differences", robustMedian);
-    psMetadataAddVector(stats, PS_LIST_TAIL, "FAKE.RS", PS_META_REPLACE,
-                        "Robust stdev in magnitude differences", robustStdev);
-    psMetadataAddVector(stats, PS_LIST_TAIL, "FAKE.MEANERR", PS_META_REPLACE,
-                        "Mean error in magnitude differences", magErrErrMean);
+    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);
-    psFree(magErrStdev);
-
-    psFree(robustMedian);
-    psFree(robustStdev);
 
     return true;