Index: /branches/pap_branch_20090108/ppStack/src/ppStackMatch.c
===================================================================
--- /branches/pap_branch_20090108/ppStack/src/ppStackMatch.c	(revision 21095)
+++ /branches/pap_branch_20090108/ppStack/src/ppStackMatch.c	(revision 21096)
@@ -16,5 +16,5 @@
 #define FAINT_SOURCE_FRAC 1.0e-4         // Set minimum flux to this fraction of faintest source flux
 
-//#define TESTING                         // Enable debugging output
+#define TESTING                         // Enable debugging output
 
 #ifdef TESTING
@@ -104,4 +104,5 @@
 
 // Generate a background model of the readout we're matching
+// We only need to set the background around the sources, not everywhere.
 psImage *stackBackgroundModel(pmReadout *ro, psMaskType maskVal, const psArray *sources, int size)
 {
@@ -355,5 +356,4 @@
             psFree(bgImage);
 
-
 #ifdef TESTING
             {
@@ -580,4 +580,68 @@
     psFree(bg);
 
+#define RADIUS 10                       // Radius of photometry
+#define MIN_ERR 0.05                    // Minimum photometric error, mag
+
+    // 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) {
+            continue;
+        }
+
+        float xSrc = source->peak->xf, ySrc = source->peak->yf; // Source coordinates
+        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);
+
+#ifdef TESTING
+    {
+        psString name = NULL;
+        psStringAppend(&name, "convolved_%03d.fits", numInput);
+        psFits *fits = psFitsOpen(name, "w");
+        psFree(name);
+        psFitsWriteImage(fits, NULL, output->image, 0, NULL);
+        psFitsClose(fits);
+    }
+#endif
+
     psFree(output);