Changeset 25302

Timestamp:

Sep 8, 2009, 6:24:00 PM (17 years ago)

Author:

Paul Price

Message:

Photometer fake sources to measure magnitude biases near limiting magnitude. Had to change a few prototypes to include 'const'.

Location:

branches/pap/psphot/src

Files:

• psphot.h (modified) (3 diffs)
• psphotFake.c (modified) (8 diffs)
• psphotFitSourcesLinear.c (modified) (1 diff)
• psphotMagnitudes.c (modified) (2 diffs)
• psphotReadout.c (modified) (1 diff)

branches/pap/psphot/src/psphot.h

@@ -49 +49 @@
 bool            psphotPSFstats (pmReadout *readout, psMetadata *recipe, pmPSF *psf);
 bool            psphotMomentsStats (pmReadout *readout, psMetadata *recipe, psArray *sources);
-bool            psphotFitSourcesLinear (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, bool final);
+bool            psphotFitSourcesLinear (pmReadout *readout, psArray *sources, const psMetadata *recipe, const pmPSF *psf, bool final);
 bool            psphotReplaceUnfitSources (psArray *sources, psImageMaskType maskVal);
 
@@ -64 +64 @@
 bool            psphotGuessModel_Threaded (psThreadJob *job);
 
-bool            psphotMagnitudes (pmConfig *config, pmReadout *readout, const pmFPAview *view, psArray *sources, pmPSF *psf);
+bool            psphotMagnitudes (pmConfig *config, pmReadout *readout, const pmFPAview *view, psArray *sources, const pmPSF *psf);
 bool            psphotMagnitudes_Threaded (psThreadJob *job);
 
@@ -76 +76 @@
 bool            psphotExtendedSourceAnalysis (pmReadout *readout, psArray *sources, psMetadata *recipe);
 bool            psphotExtendedSourceFits (pmReadout *readout, psArray *sources, psMetadata *recipe);
-bool            psphotFake(pmReadout *readout, const pmPSF *psf, psMetadata *recipe, const psArray *realSources);
+bool            psphotFake(pmConfig *config, pmReadout *readout, const pmFPAview *view, const pmPSF *psf, psMetadata *recipe, const psArray *realSources);
 
 // thread-related:

branches/pap/psphot/src/psphotFake.c

@@ -21 +21 @@
     PM_ASSERT_READOUT_NON_NULL(ro, false);
     PM_ASSERT_READOUT_VARIANCE(ro, false);
+    assert(magLim);
+    assert(radius);
 
     psKernel *kernel = psImageSmoothKernel(smoothSigma, smoothNsigma); // Kernel used for smoothing
@@ -52 +54 @@
     psFree(normModel);
 
-    if (magLim) {
-        // The signal-to-noise of the smoothed image is: S/N ~ I_smooth / sqrt(variance * factor)
-
-        // since the variance factor tells us the variance in the smoothed image.  Now, the trick is working
-        // out what the intensity in the smoothed image is, and how it is related to the flux.  We are
-        // convolving Io G_* with G_PSF/2pi.w^2, where G_* is the approximately-Gaussian PSF of the star,
-        // G_PSF is the pretty-close-matching Gaussian of the convolution kernel, Io is the peak flux in the
-        // unsmoothed image, and w is the width of the Gaussian.  Now, a normalised 2D Gaussian convolved
-        // with itself has a normalisation of 1/2pi(w^2+w^2) = 1/4pi.w^2.  Therefore:
-        // 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
-        *magLim = -2.5 * log10f(fluxLim);
-        psTrace("psphot.fake", 1, "Limiting peak: %f\n", peakLim);
-        psTrace("psphot.fake", 1, "Limiting flux: %f\n", fluxLim);
-        psTrace("psphot.fake", 1, "Limiting mag: %f\n", *magLim);
-    }
-    if (radius) {
-        *radius = smoothSigma * smoothNsigma;
-    }
+    // The signal-to-noise of the smoothed image is: S/N ~ I_smooth / sqrt(variance * factor)
+    // since the variance factor tells us the variance in the smoothed image.  Now, the trick is working
+    // out what the intensity in the smoothed image is, and how it is related to the flux.  We are
+    // convolving Io G_* with G_PSF/2pi.w^2, where G_* is the approximately-Gaussian PSF of the star,
+    // G_PSF is the pretty-close-matching Gaussian of the convolution kernel, Io is the peak flux in the
+    // unsmoothed image, and w is the width of the Gaussian.  Now, a normalised 2D Gaussian convolved
+    // with itself has a normalisation of 1/2pi(w^2+w^2) = 1/4pi.w^2.  Therefore:
+    // 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
+    *magLim = -2.5 * log10f(fluxLim);
+    psTrace("psphot.fake", 1, "Limiting peak: %f\n", peakLim);
+    psTrace("psphot.fake", 1, "Limiting flux: %f\n", fluxLim);
+    psTrace("psphot.fake", 1, "Limiting mag: %f\n", *magLim);
+
+    *radius = smoothSigma * smoothNsigma;
 
     return true;
@@ -90 +88 @@
     PS_ASSERT_VECTOR_NON_NULL(magOffsets, false);
     PS_ASSERT_VECTOR_TYPE(magOffsets, PS_TYPE_F32, false);
+    assert(xSrc);
+    assert(ySrc);
 
     int numBins = magOffsets->n;                                    // Number of bins
@@ -140 +140 @@
 
 // *** in this section, perform the photometry for fake sources ***
-bool psphotFake(pmReadout *readout, const pmPSF *psf, psMetadata *recipe, const psArray *realSources)
+bool psphotFake(pmConfig *config, pmReadout *readout, const pmFPAview *view, const pmPSF *psf,
+                psMetadata *recipe, const psArray *realSources)
 {
     PM_ASSERT_READOUT_NON_NULL(readout, false);
@@ -201 +202 @@
     }
 
+    // XXX Could speed this up significantly by only convolving the central pixels of each fake source
     psImage *significance = psphotSignificanceImage(readout, recipe, 2, maskVal);
     if (!significance) {
@@ -213 +215 @@
     int numBins = magOffsets->n;                          // Number of bins
     int numCols = readout->image->numCols, numRows = readout->image->numRows; // Size of image
-    psVector *frac = psVectorAlloc(numBins, PS_TYPE_F32); // Fraction of sources in each bin
+    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
     for (int i = 0; i < numBins; i++) {
-        int numFound = 0;               // Number found
+        int num = 0;               // Number found
         for (int j = 0; j < numSources; j++) {
             // Coordinates of interest
@@ -221 +226 @@
             float y = PS_MIN(yFake->data.F32[i][j] + 0.5, numRows - 1);
             int xPix = x, yPix = y;     // Pixel coordinates
-            if (significance->data.F32[yPix][xPix] > thresh) {
+            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)) {
+                    psErrorClear();
+                    continue;
+                }
+                source->modelPSF = pmModelFromPSFforXY(psf, x, y, 2.0 * sqrtf(sig));
+                source->type = PM_SOURCE_TYPE_STAR;
+
+                num++;
+
+                fakeSources->data[numFound] = source;
                 numFound++;
             }
         }
-        frac->data.F32[i] = (float)numFound / (float)numSources;
-    }
+        count->data.S32[i] = num;
+        fakeEnd->data.S32[i] = numFound;
+    }
+    fakeSources->n = numFound;
 
     psFree(xFake);
@@ -232 +252 @@
     psFree(significance);
 
+    if (!psphotFitSourcesLinear(readout, fakeSources, 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)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to measure magnitudes of fake sources.");
+        psFree(fakeSources);
+        psFree(fakeEnd);
+        psFree(count);
+        return false;
+    }
+
+    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 *magErr = psVectorAlloc(numSources, PS_TYPE_F32);   // Magnitude errors
+    psVector *magMask = psVectorAlloc(numSources, PS_TYPE_VECTOR_MASK); // Mask for magnitude errors
+    psStats *magStats = psStatsAlloc(PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV); // Statistics
+    for (int i = 0, index = 0; i < numBins; i++) {
+        psStatsInit(magStats);
+        psVectorInit(magMask, 0);
+
+        float magRef = magLim + magOffsets->data.F32[i]; // Reference magnitude
+        for (int j = 0; index < fakeEnd->data.S32[i]; j++, index++) {
+            pmSource *source = fakeSources->data[index]; // Source of interest
+            if (!source || !isfinite(source->psfMag)) {
+                magMask->data.PS_TYPE_VECTOR_MASK_DATA[j] = 0xFF;
+                continue;
+            }
+            magErr->data.F32[j] = source->psfMag - magRef;
+        }
+        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;
+    }
+    psFree(magStats);
+    psFree(magMask);
+    psFree(magErr);
+
+    psFree(fakeSources);
+    psFree(fakeEnd);
+
     // XXX How do we get the results out?
     psMetadata *stats = psMetadataAlloc();
-    psMetadataAddVector(stats, PS_LIST_TAIL, "FAKE.EFF", PS_META_REPLACE,
-                        "Efficiency fractions", frac);
+    psMetadataAddS32(stats, PS_LIST_TAIL, "FAKE.NUM", PS_META_REPLACE,
+                     "Number of sources per bin", numSources);
+    psMetadataAddF32(stats, PS_LIST_TAIL, "FAKE.REF", PS_META_REPLACE,
+                     "Efficiency reference magnitude", magLim);
     psMetadataAddVector(stats, PS_LIST_TAIL, "FAKE.MAG", PS_META_REPLACE,
                         "Efficiency magnitudes", magOffsets);
-    psMetadataAddF32(stats, PS_LIST_TAIL, "FAKE.REF", PS_META_REPLACE,
-                     "Efficiency reference magnitude", magLim);
+    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);
     psMetadataConfigWrite(stats, "fake.stats");
     psFree(stats);
 
-    psFree(frac);
+    psFree(count);
+    psFree(magErrMean);
+    psFree(magErrStdev);
 
     return true;

branches/pap/psphot/src/psphotFitSourcesLinear.c

@@ -12 +12 @@
 static bool SetBorderMatrixElements (psSparseBorder *border, pmReadout *readout, psArray *sources, bool constant_weights, int SKY_FIT_ORDER, psImageMaskType markVal);
 
-bool psphotFitSourcesLinear (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, bool final) {
+bool psphotFitSourcesLinear (pmReadout *readout, psArray *sources, const psMetadata *recipe, const pmPSF *psf, bool final) {
 
     bool status;

branches/pap/psphot/src/psphotMagnitudes.c

@@ -1 +1 @@
 # include "psphotInternal.h"
 
-bool psphotMagnitudes(pmConfig *config, pmReadout *readout, const pmFPAview *view, psArray *sources, pmPSF *psf) {
+bool psphotMagnitudes(pmConfig *config, pmReadout *readout, const pmFPAview *view, psArray *sources, const pmPSF *psf) {
 
     bool status = false;
@@ -64 +64 @@
 
             psArrayAdd(job->args, 1, cells->data[j]); // sources
-            psArrayAdd(job->args, 1, psf);
+            psArrayAdd(job->args, 1, (pmPSF*)psf);    // Casting away const
             psArrayAdd(job->args, 1, binning);
             psArrayAdd(job->args, 1, backModel);

branches/pap/psphot/src/psphotReadout.c

@@ -224 +224 @@
     psphotMagnitudes(config, readout, view, sources, psf);
 
-    if (!psphotFake(readout, psf, recipe, sources)) {
+    if (!psphotFake(config, readout, view, psf, recipe, sources)) {
         psErrorStackPrint(stderr, "Unable to determine detection efficiencies from fake sources");
         psErrorClear();