Changeset 21199

Timestamp:

Jan 28, 2009, 10:49:50 AM (17 years ago)

Author:

Paul Price

Message:

Merging pap_branch_20090108: fixed photometric normalisation. Conflicts resolved, code compiles.

Location:

trunk/ppStack/src

Files:

• ppStackArguments.c (modified) (2 diffs)
• ppStackLoop.c (modified) (23 diffs)
• ppStackMatch.c (modified) (16 diffs)
• ppStackSources.c (modified) (5 diffs)

trunk/ppStack/src/ppStackArguments.c

@@ -152 +152 @@
     psMetadataAddF32(arguments, PS_LIST_TAIL, "-threshold-mask", 0, "Threshold for mask deconvolution", NAN);
     psMetadataAddF32(arguments, PS_LIST_TAIL, "-poor-frac", 0, "Fraction of weight for poor pixels", NAN);
+    psMetadataAddF32(arguments, PS_LIST_TAIL, "-deconv-limit", 0, "Maximum deconvolution fraction limit", NAN);
     psMetadataAddF32(arguments, PS_LIST_TAIL, "-image-rej", 0,
                      "Pixel rejection fraction threshold for rejecting entire image", NAN);
@@ -236 +237 @@
     VALUE_ARG_RECIPE_FLOAT("-threshold-mask", "THRESHOLD.MASK", F32);
     VALUE_ARG_RECIPE_FLOAT("-image-rej",      "IMAGE.REJ",      F32);
+    VALUE_ARG_RECIPE_FLOAT("-deconv-limit",   "DECONV.LIMIT",   F32);
     VALUE_ARG_RECIPE_INT("-rows",             "ROWS",           S32, 0);
     VALUE_ARG_RECIPE_FLOAT("-poor-frac",      "POOR.FRACTION",  F32);

trunk/ppStack/src/ppStackLoop.c

@@ -357 +357 @@
         }
 
+        // Generate target PSF
         targetPSF = ppStackPSF(config, numCols, numRows, psfs);
         psFree(psfs);
@@ -365 +366 @@
             return false;
         }
+        psMetadataAddPtr(config->arguments, PS_LIST_TAIL, "PSF.TARGET", PS_DATA_UNKNOWN,
+                         "Target PSF for stack", targetPSF);
 
         pmChip *outChip = pmFPAviewThisChip(view, output->fpa); // Output chip
@@ -423 +426 @@
             psFree(inputMask);
             psFree(matchChi2);
+            psFree(fpaList);
+            psFree(cellList);
             return false;
         }
@@ -440 +445 @@
             psFree(inputMask);
             psFree(matchChi2);
+            psFree(fpaList);
+            psFree(cellList);
             return false;
         }
@@ -446 +453 @@
         psArray *regions = NULL, *kernels = NULL; // Regions and kernels used in subtraction
         psTimerStart("PPSTACK_MATCH");
+
         if (!ppStackMatch(readout, &regions, &kernels, &matchChi2->data.F32[i], &weightings->data.F32[i],
                           sourceLists->data[i], targetPSF, rng, config)) {
@@ -476 +484 @@
         assert(hdu);
         writeImage(imageNames->data[i],  hdu->header, readout->image, config);
-        writeImage(maskNames->data[i],   hdu->header, readout->mask, config);
+        psMetadata *maskHeader = psMetadataCopy(NULL, hdu->header); // Copy of header, for mask
+        pmConfigMaskWriteHeader(config, maskHeader);
+        writeImage(maskNames->data[i],   maskHeader, readout->mask, config);
+        psFree(maskHeader);
         writeImage(weightNames->data[i], hdu->header, readout->weight, config);
 
@@ -486 +497 @@
         cells->data[i] = psMemIncrRefCounter(inCell);
         if (!filesIterateUp(config)) {
+            psFree(fpaList);
+            psFree(cellList);
             return false;
         }
@@ -503 +516 @@
         psFree(inputMask);
         psFree(matchChi2);
+        psFree(fpaList);
+        psFree(cellList);
         return false;
     }
@@ -569 +584 @@
             numGood = 0;                    // Number of good images
             for (int i = 0; i < num; i++) {
-              if (inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PPSTACK_MASK_ALL) {
+              if (inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PPSTACK_MASK_ALL) {
                     continue;
                 }
@@ -612 +627 @@
     pmFPAview *view = NULL;             // View to readout
     psArray *inspect = NULL;            // Array of arrays of pixels to inspect
-    psVector *exptimes = psVectorAlloc(num, PS_TYPE_F32); // Exposure times for each input
     {
         int row0, col0;                 // Offset for readout
@@ -624 +638 @@
             psFree(inputMask);
             psFree(matchChi2);
-            psFree(exptimes);
             psFree(cells);
             return false;
@@ -638 +651 @@
             psFree(inputMask);
             psFree(matchChi2);
-            psFree(exptimes);
             psFree(cells);
             return false;
@@ -657 +669 @@
             psFree(view);
             psFree(outRO);
-            psFree(exptimes);
             psFree(cells);
             return false;
@@ -678 +689 @@
                 psFree(view);
                 psFree(outRO);
-                psFree(exptimes);
                 return false;
             }
@@ -705 +715 @@
                 psFree(view);
                 psFree(outRO);
-                psFree(exptimes);
                 return false;
             }
@@ -720 +729 @@
             psFree(view);
             psFree(outRO);
-            psFree(exptimes);
             return false;
         }
@@ -788 +796 @@
                 psFree(inspect);
                 psFree(rejected);
-                psFree(exptimes);
                 return false;
             }
@@ -804 +811 @@
             psFree(inspect);
             psFree(rejected);
-            psFree(exptimes);
             return false;
         }
@@ -899 +905 @@
             psFree(view);
             psFree(outRO);
-            psFree(exptimes);
             return false;
         }
@@ -928 +933 @@
                 psFree(view);
                 psFree(outRO);
-                psFree(exptimes);
                 return false;
             }
@@ -950 +954 @@
                 psFree(view);
                 psFree(outRO);
-                psFree(exptimes);
                 return false;
             }
@@ -963 +966 @@
             psFree(view);
             psFree(outRO);
-            psFree(exptimes);
             return false;
         }
@@ -1042 +1044 @@
     psFree(weightNames);
 
-    psFree(exptimes);
     psFree(inputMask);
     psFree(stack);

trunk/ppStack/src/ppStackMatch.c

@@ -6 +6 @@
 #include <pslib.h>
 #include <psmodules.h>
+#include <psphot.h>
 
 #include "ppStack.h"
@@ -47 +48 @@
 #endif
 
+// Get coordinates from a source
+static void coordsFromSource(float *x, float *y, const pmSource *source)
+{
+    assert(x && y);
+    assert(source);
+
+    if (source->modelPSF) {
+        *x = source->modelPSF->params->data.F32[PM_PAR_XPOS];
+        *y = source->modelPSF->params->data.F32[PM_PAR_YPOS];
+    } else {
+        *x = source->peak->xf;
+        *y = source->peak->yf;
+    }
+    return;
+}
+
 static psArray *stackSourcesFilter(psArray *sources, // Source list to filter
                                    int exclusion // Exclusion zone, pixels
@@ -64 +81 @@
             continue;
         }
-        x->data.F32[numGood] = source->peak->xf;
-        y->data.F32[numGood] = source->peak->yf;
+        coordsFromSource(&x->data.F32[numGood], &y->data.F32[numGood], source);
         numGood++;
     }
@@ -80 +96 @@
             continue;
         }
-        coords->data.F64[0] = source->peak->xf;
-        coords->data.F64[1] = source->peak->yf;
+        float xSource, ySource;         // Coordinates of source
+        coordsFromSource(&xSource, &ySource, source);
+
+        coords->data.F64[0] = xSource;
+        coords->data.F64[1] = ySource;
 
         long numWithin = psTreeWithin(tree, coords, exclusion); // Number within exclusion zone
@@ -101 +120 @@
 }
 
-#define BG_SIZE 64                      // Large box half-size in which to measure background
-
-// Generate a background model of the readout we're matching
-psImage *stackBackgroundModel(pmReadout *ro, psImageMaskType maskVal, const psArray *sources, int size)
+// Add background into the fake image
+// Based on ppSubBackground()
+static psImage *stackBackgroundModel(pmReadout *ro, // Readout for which to generate background model
+                                     const pmConfig *config // Configuration
+    )
 {
-    psImage *image = ro->image, *mask = ro->mask; // Image and mask of readout
+    psAssert(ro && ro->image, "Need readout image");
+    psAssert(config, "Need configuration");
+
+    psImage *image = ro->image;         // Image of interest
     int numCols = image->numCols, numRows = image->numRows; // Size of image
 
-    psImage *bgImage = psImageAlloc(numCols, numRows, PS_TYPE_F32); // Background image
-    psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); // Statistics for measuring background
-    psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS, 0); // Random number generator
-
-    for (int i = 0; i < sources->n; i++) {
-        pmSource *source = sources->data[i]; // Source of interest
-        if (!source) {
-            continue;
-        }
-
-        int x = source->peak->xf + 0.5, y = source->peak->yf + 0.5; // Coordinates of source
-
-        int xMin = PS_MAX(x - BG_SIZE, 0), xMax = PS_MIN(x + BG_SIZE, numCols);
-        int yMin = PS_MAX(y - BG_SIZE, 0), yMax = PS_MIN(y + BG_SIZE, numCols);
-
-        psRegion region = psRegionSet(xMin, xMax, yMin, yMax); // Region of interest
-        psImage *subImage = psImageSubset(image, region); // Subimage
-        psImage *subMask = psImageSubset(mask, region); // Subimage mask
-        if (!psImageBackground(stats, NULL, subImage, subMask, maskVal, rng)) {
-            // Can't do anything about it
-            psErrorClear();
-            continue;
-        }
-        psFree(subImage);
-        psFree(subMask);
-
-        float value = stats->robustMedian; // Value to set background to
-
-        int uMin = PS_MAX(x - size, 0), uMax = PS_MIN(x + size, numCols);
-        int vMin = PS_MAX(y - size, 0), vMax = PS_MIN(y + size, numCols);
-
-        for (int v = vMin; v < vMax; v++) {
-            for (int u = uMin; u < uMax; u++) {
-                bgImage->data.F32[v][u] = value;
-            }
-        }
-    }
-
-    psFree(stats);
-    psFree(rng);
-
-    return bgImage;
+    psMetadata *ppStackRecipe = psMetadataLookupPtr(NULL, config->recipes, PPSTACK_RECIPE);
+    psAssert(ppStackRecipe, "Need PPSTACK recipe");
+    psMetadata *psphotRecipe = psMetadataLookupPtr(NULL, config->recipes, PSPHOT_RECIPE);
+    psAssert(psphotRecipe, "Need PSPHOT recipe");
+
+    psString maskBadStr = psMetadataLookupStr(NULL, ppStackRecipe, "MASK.BAD");// Name of bits to mask for bad
+    psMaskType maskBad = pmConfigMaskGet(maskBadStr, config); // Bits to mask for bad pixels
+
+    // user-defined masks to test for good/bad pixels (build from recipe list if not yet set)
+    psMetadataAddU8(psphotRecipe, PS_LIST_TAIL, "MASK.PSPHOT", PS_META_REPLACE, "user-defined mask", maskBad);
+
+    psImage *binned = psphotBackgroundModel(ro, config); // Binned background model
+    psImageBinning *binning = psMetadataLookupPtr(NULL, psphotRecipe,
+                                                  "PSPHOT.BACKGROUND.BINNING"); // Binning for model
+    psAssert(binning, "Need binning parameters");
+    psImage *unbinned = psImageAlloc(numCols, numRows, PS_TYPE_F32); // Unbinned background model
+    if (!psImageUnbin(unbinned, binned, binning)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to unbin background model");
+        psFree(unbinned);
+        return NULL;
+    }
+
+    // XXX should these really be here?? (probably not...)
+    // pmFPAfileDropInternal(config->files, "PSPHOT.BACKMDL");
+    // pmFPAfileDropInternal(config->files, "PSPHOT.BACKMDL.STDEV");
+
+     return unbinned;
 }
+
 
 bool ppStackMatch(pmReadout *readout, psArray **regions, psArray **kernels, float *chi2, float *weighting,
@@ -163 +172 @@
     psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSTACK_RECIPE); // ppStack recipe
     psAssert(recipe, "We've thrown an error on this before.");
+
+    float deconvLimit = psMetadataLookupF32(NULL, recipe, "DECONV.LIMIT"); // Limit on deconvolution fraction
 
     // Look up appropriate values from the ppSub recipe
@@ -195 +206 @@
         assert(outName);
         // Read convolution kernel
-        {
-            psString filename = NULL;   // Output filename
-            psStringAppend(&filename, "%s.%d.kernel", outName, numInput);
-            psString resolved = pmConfigConvertFilename(filename, config, false, false); // Resolved filename
-            psFree(filename);
-            psFits *fits = psFitsOpen(resolved, "r"); // FITS file for subtraction kernel
-            psFree(resolved);
-            if (!fits || !pmReadoutReadSubtractionKernels(output, fits)) {
-                psError(PS_ERR_IO, false, "Unable to read previously produced kernel");
-                psFitsClose(fits);
-                return false;
-            }
+        psString filename = NULL;   // Output filename
+        psStringAppend(&filename, "%s.%d.kernel", outName, numInput);
+        psString resolved = pmConfigConvertFilename(filename, config, false, false); // Resolved filename
+        psFree(filename);
+        psFits *fits = psFitsOpen(resolved, "r"); // FITS file for subtraction kernel
+        psFree(resolved);
+        if (!fits || !pmReadoutReadSubtractionKernels(output, fits)) {
+            psError(PS_ERR_IO, false, "Unable to read previously produced kernel");
             psFitsClose(fits);
-
-            // Add in variance factor
-            pmSubtractionKernels *kernels = psMetadataLookupPtr(NULL, output->analysis,
-                                                                PM_SUBTRACTION_ANALYSIS_KERNEL); // Kernels
-            float vf = pmSubtractionVarianceFactor(kernels, 0.0, 0.0, false); // Variance factor
-            psMetadataItem *vfItem = psMetadataLookup(readout->parent->concepts, "CELL.VARFACTOR");
-            if (!isfinite(vf)) {
-                vf = 1.0;
-            }
-            if (isfinite(vfItem->data.F32)) {
-                vfItem->data.F32 *= vf;
-            } else {
-                vfItem->data.F32 = vf;
-            }
+            return false;
+        }
+        psFitsClose(fits);
+
+        // Add in variance factor
+        pmSubtractionKernels *kernels = psMetadataLookupPtr(NULL, output->analysis,
+                                                            PM_SUBTRACTION_ANALYSIS_KERNEL); // Kernels
+        float vf = pmSubtractionVarianceFactor(kernels, 0.0, 0.0, false); // Variance factor
+        psMetadataItem *vfItem = psMetadataLookup(readout->parent->concepts, "CELL.VARFACTOR");
+        if (!isfinite(vf)) {
+            vf = 1.0;
+        }
+        if (isfinite(vfItem->data.F32)) {
+            vfItem->data.F32 *= vf;
+        } else {
+            vfItem->data.F32 = vf;
         }
 
@@ -244 +253 @@
         psFree(maskName);
         psFree(weightName);
+
+        psRegion *region = psMetadataLookupPtr(NULL, output->analysis,
+                                               PM_SUBTRACTION_ANALYSIS_REGION); // Convolution region
+
+        pmSubtractionAnalysis(readout->analysis, kernels, region,
+                              readout->image->numCols, readout->image->numRows);
     } else {
 #endif
@@ -285 +300 @@
             }
 
-            // Generate a fake image to match to
-            float minFlux = INFINITY;       // Minimum flux for fake image
-            {
-                psStats *bg = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); // Statistics for bg
-                if (!psImageBackground(bg, NULL, readout->image, readout->mask, maskVal, rng)) {
-                    psWarning("Can't measure background for image.");
-                    psErrorClear();
-                } else {
-                    minFlux = bg->robustStdev;
-                }
-                psFree(bg);
-            }
-
-#if 0
-            float maxMag = -INFINITY;   // Maximum magnitude of sources
-            for (int i = 0; i < sources->n; i++) {
-                pmSource *source = sources->data[i]; // Source of interest
-                if (source->psfMag > maxMag && source->psfMag <= MAG_IGNORE &&
-                    !(source->mode & SOURCE_MASK)) {
-                    maxMag = source->psfMag;
-                }
-            }
-            minFlux = PS_MIN(FAINT_SOURCE_FRAC * powf(10.0, -0.4 * maxMag), minFlux);
-#endif
-
-
 #if 0
             // Testing the normalisation of the fake image
@@ -318 +307 @@
                 pmSource *source = pmSourceAlloc();     // Source
                 sources->data[0] = source;
-                source->peak = pmPeakAlloc(50, 50, 10000, PM_PEAK_LONE);
+                source->peak = pmPeakAlloc(500, 500, 10000, PM_PEAK_LONE);
                 source->psfMag = -13.0;
-                pmReadoutFakeFromSources(test, 100, 100, sources, NULL, NULL, psf, 0.1, 0, false, true);
+                pmReadoutFakeFromSources(test, 1000, 1000, sources, NULL, NULL, psf, 0.1, 0, false, true);
                 float sum = 0.0;
                 for (int y = 0; y < test->image->numRows; y++) {
@@ -342 +331 @@
 
             if (!pmReadoutFakeFromSources(fake, readout->image->numCols, readout->image->numRows,
-                                          stampSources, NULL, NULL, psf, minFlux, 0, false, true)) {
+                                          stampSources, NULL, NULL, psf, NAN, footprint + size,
+                                          false, true)) {
                 psError(PS_ERR_UNKNOWN, false, "Unable to generate fake image with target PSF.");
                 psFree(fake);
@@ -351 +341 @@
 
             // Add the background into the target image
-            psImage *bgImage = stackBackgroundModel(readout, maskVal, stampSources, footprint); // Image of background
+            psImage *bgImage = stackBackgroundModel(readout, config); // Image of background
             psBinaryOp(fake->image, fake->image, "+", bgImage);
             psFree(bgImage);
 
-
 #ifdef TESTING
             {
+                pmHDU *hdu = pmHDUFromCell(readout->parent);
                 psString name = NULL;
                 psStringAppend(&name, "fake_%03d.fits", numInput);
                 psFits *fits = psFitsOpen(name, "w");
                 psFree(name);
-                psFitsWriteImage(fits, NULL, fake->image, 0, NULL);
+                psFitsWriteImage(fits, hdu->header, fake->image, 0, NULL);
                 psFitsClose(fits);
             }
             {
+                pmHDU *hdu = pmHDUFromCell(readout->parent);
                 psString name = NULL;
                 psStringAppend(&name, "real_%03d.fits", numInput);
                 psFits *fits = psFitsOpen(name, "w");
                 psFree(name);
-                psFitsWriteImage(fits, NULL, readout->image, 0, NULL);
+                psFitsWriteImage(fits, hdu->header, readout->image, 0, NULL);
                 psFitsClose(fits);
             }
@@ -410 +401 @@
 #ifdef TESTING
             {
+                pmHDU *hdu = pmHDUFromCell(readout->parent);
                 psString name = NULL;
                 psStringAppend(&name, "conv_%03d.fits", numInput);
                 psFits *fits = psFitsOpen(name, "w");
                 psFree(name);
-                psFitsWriteImage(fits, NULL, output->image, 0, NULL);
+                psFitsWriteImage(fits, hdu->header, output->image, 0, NULL);
                 psFitsClose(fits);
             }
             {
+                pmHDU *hdu = pmHDUFromCell(readout->parent);
                 psString name = NULL;
                 psStringAppend(&name, "diff_%03d.fits", numInput);
@@ -423 +416 @@
                 psFree(name);
                 psBinaryOp(fake->image, output->image, "-", fake->image);
-                psFitsWriteImage(fits, NULL, fake->image, 0, NULL);
+                psFitsWriteImage(fits, hdu->header, fake->image, 0, NULL);
                 psFitsClose(fits);
             }
@@ -548 +541 @@
     }
 
+    // Reject image completely if the maximum deconvolution fraction exceeds the limit
+    float deconv = psMetadataLookupF32(NULL, output->analysis,
+                                       PM_SUBTRACTION_ANALYSIS_DECONV_MAX); // Maximum deconvolution fraction
+    if (deconv > deconvLimit) {
+        psWarning("Maximum deconvolution fraction (%f) exceeds limit (%f) --- rejecting\n",
+                  deconv, deconvLimit);
+        psFree(output);
+        return NULL;
+    }
+
     // Renormalise the variances if desired
     if (renorm) {
@@ -580 +583 @@
     psFree(bg);
 
+
+#if 0
+#define RADIUS 10                       // Radius of photometry
+#define MIN_ERR 0.05                    // Minimum photometric error, mag
+#define MAX_MAG -13                     // Maximum magnitude for source
+
+    // Ensure the normalisation is correct
+    // XXX Ideally, would like to do proper PSF-fit photometry, but will settle for aperture photometry
+    int numSources = sources->n;        // Number of sources
+    psVector *ratio = psVectorAlloc(numSources, PS_TYPE_F32); // Flux ratios for sources
+    psVector *ratioMask = psVectorAlloc(numSources, PS_TYPE_MASK); // Mask for flux ratios
+    psVectorInit(ratioMask, 0xFF);
+    psImage *image = readout->image;    // Image of interest
+    psImage *mask = readout->mask;      // Mask of interest
+    int numCols = image->numCols, numRows = image->numRows; // Size of image
+    for (int i = 0; i < numSources; i++) {
+        pmSource *source = sources->data[i]; // Source of interest
+        if (!source || source->mode & SOURCE_MASK || !isfinite(source->psfMag) || !isfinite(source->errMag) ||
+            source->errMag > MIN_ERR || source->psfMag > MAX_MAG) {
+            continue;
+        }
+
+        float xSrc, ySrc;              // Source coordinates
+        coordsFromSource(&xSrc, &ySrc, source);
+        int xMin = PS_MAX(0, xSrc - RADIUS), xMax = PS_MIN(numCols - 1, xSrc + RADIUS); // Bounds in x
+        int yMin = PS_MAX(0, ySrc - RADIUS), yMax = PS_MIN(numRows - 1, ySrc + RADIUS); // Bounds in y
+        int numPix = 0;                 // Number of pixels
+        float sum = 0.0;                // Sum of pixels
+        for (int y = yMin; y <= yMax; y++) {
+            for (int x = xMin; x <= xMax; x++) {
+                if (mask->data.PS_TYPE_MASK_DATA[y][x] & maskBad) {
+                    continue;
+                }
+                sum += image->data.F32[y][x];
+                numPix++;
+            }
+        }
+        if (sum >= 0 && numPix > 0) {
+            float mag = -2.5 * log10(sum * M_PI * PS_SQR(RADIUS) / numPix); // Instrumental magnitude
+            ratio->data.F32[i] = mag - source->psfMag;
+            ratioMask->data.PS_TYPE_MASK_DATA[i] = 0;
+        }
+    }
+
+    psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); // Statistics
+    if (!psVectorStats(stats, ratio, NULL, ratioMask, 0xFF)) {
+        psWarning("Unable to measure normalisation --- assuming correct.");
+    } else {
+        psLogMsg("ppStack", PS_LOG_INFO, "Renormalising image by %f (+/- %f) mag\n",
+                 stats->robustMedian, stats->robustStdev);
+        float norm = powf(10.0, -0.4 * stats->robustMedian); // Normalisation to apply
+        psBinaryOp(image, image, "*", psScalarAlloc(norm, PS_TYPE_F32));
+    }
+    psFree(stats);
+    psFree(ratio);
+    psFree(ratioMask);
+#endif
+
+#ifdef TESTING
+    {
+        pmHDU *hdu = pmHDUFromCell(readout->parent);
+        psString name = NULL;
+        psStringAppend(&name, "convolved_%03d.fits", numInput);
+        psFits *fits = psFitsOpen(name, "w");
+        psFree(name);
+        psFitsWriteImage(fits, hdu->header, output->image, 0, NULL);
+        psFitsClose(fits);
+    }
+#endif
+
     psFree(output);
 

trunk/ppStack/src/ppStackSources.c

@@ -9 +9 @@
 //#define TESTING                         // Enable debugging output
 
+// Size of fake image; set by hand because it's trouble to get it from other places
+#define FAKE_COLS 4861
+#define FAKE_ROWS 4913
 
 float ppStackSourcesTransparency(const psArray *sourceLists, const pmFPAview *view, const pmConfig *config)
@@ -60 +63 @@
         pmCell *cell = pmFPAviewThisCell(view, file->fpa); // Cell of interest
 
+#ifdef TESTING
+        pmReadout *fake = pmReadoutAlloc(NULL); // Fake readout
+        pmPSF *psf = psMetadataLookupPtr(NULL, config->arguments, "PSF.TARGET"); // PSF for fake image
+        pmReadoutFakeFromSources(fake, FAKE_COLS, FAKE_ROWS, sourceLists->data[i], NULL, NULL, psf, 5, 0, false, true);
+        psString name = NULL;
+        psStringAppend(&name, "start_%03d.fits", i);
+        psFits *fits = psFitsOpen(name, "w");
+        psFree(name);
+        psFitsWriteImage(fits, NULL, fake->image, 0, NULL);
+        psFitsClose(fits);
+        psFree(fake);
+#endif
+
+
         float exptime = psMetadataLookupF32(NULL, cell->concepts, "CELL.EXPOSURE"); // Exposure time
         float airmass = psMetadataLookupF32(NULL, file->fpa->concepts, "FPA.AIRMASS"); // Airmass
@@ -105 +122 @@
         FILE *outMatches = fopen("source_match.dat", "w"); // Output matches
         psVector *mag = psVectorAlloc(num, PS_TYPE_F32); // Magnitudes for each star
+        psVector *magErr = psVectorAlloc(num, PS_TYPE_F32); // Errors for each star
         for (int i = 0; i < matches->n; i++) {
             pmSourceMatch *match = matches->data[i]; // Match of interest
             psVectorInit(mag, NAN);
+            psVectorInit(magErr, NAN);
             for (int j = 0; j < match->num; j++) {
                 if (match->mask->data.PS_TYPE_VECTOR_MASK_DATA[j]) {
@@ -114 +133 @@
                 int index = match->image->data.U32[j]; // Image index
                 mag->data.F32[index] = match->mag->data.F32[j] - zp->data.F32[index] - trans->data.F32[index];
+                magErr->data.F32[index] = match->magErr->data.F32[j];
             }
             for (int j = 0; j < num; j++) {
-                fprintf(outMatches, "%f ", mag->data.F32[j]);
+                fprintf(outMatches, "%f (%f) ", mag->data.F32[j], magErr->data.F32[j]);
             }
             fprintf(outMatches, "\n");
         }
         psFree(mag);
+        psFree(magErr);
         fclose(outMatches);
     }
@@ -162 +183 @@
         psVector *trans = pmSourceMatchRelphot(matches, zp, iter, tol, starLimit, transIter, transRej,
                                                transThresh, starRej, starSys);
+        for (int i = 0; i < num; i++) {
+            fprintf(stderr, "Transparency of image %d: %f\n", i, trans->data.F32[i]);
+        }
         psFree(trans);
     }