Index: /branches/ipp-magic-v0/ppSub/src/ppSubReadoutPhotometry.c
===================================================================
--- /branches/ipp-magic-v0/ppSub/src/ppSubReadoutPhotometry.c	(revision 21375)
+++ /branches/ipp-magic-v0/ppSub/src/ppSubReadoutPhotometry.c	(revision 21375)
@@ -0,0 +1,159 @@
+/** @file ppSubReadoutPhotometry.c
+ *
+ *  @brief
+ *
+ *  @ingroup ppSub
+ *
+ *  @author IfA
+ *  @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2009-02-06 03:20:02 $
+ *  Copyright 2009 Institute for Astronomy, University of Hawaii
+ */
+
+#include "ppSub.h"
+
+/**
+ * Photometry stage 2: find and measure sources on the subtracted image
+ */
+bool ppSubReadoutPhotometry (pmConfig *config, psMetadata *stats, const pmFPAview *view) {
+
+    bool mdok = false;
+
+    // Look up recipe values
+    psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSUB_RECIPE); // Recipe for ppSim
+    psAssert(recipe, "We checked this earlier, so it should be here.");
+
+    if (!psMetadataLookupBool(NULL, recipe, "PHOTOMETRY")) {
+        pmFPAfileDropInternal(config->files, "PSPHOT.BACKMDL");
+        pmFPAfileDropInternal(config->files, "PSPHOT.BACKMDL.STDEV");
+        return true;
+    }
+
+    // The PSF (measured in ppSubMakePSF) is stored on the chip->analysis of PSPHOT.INPUT
+    // In order to use an incoming PSF, it must be stored on the chip->analysis of PSPHOT.PSF.LOAD
+    pmChip *psfInputChip = pmFPAfileThisChip(config->files, view, "PSPHOT.INPUT");
+    psAssert (psfInputChip, "should have been generated for ppSubMakePSF");
+
+    pmChip *psfLoadChip = pmFPAfileThisChip(config->files, view, "PSPHOT.PSF.LOAD");
+    psAssert (psfLoadChip, "PSPHOT.PSF.LOAD should have been defined in ppSubCamera");
+
+    pmPSF *psf = psMetadataLookupPtr(NULL, psfInputChip->analysis, "PSPHOT.PSF");
+    if (!psf) {
+        psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find PSF from psphot");
+        return false;
+    }
+    psMetadataAddPtr(psfLoadChip->analysis, PS_LIST_TAIL, "PSPHOT.PSF", PS_DATA_UNKNOWN | PS_META_REPLACE, "PSF from matched addition", psf);
+
+    // psphotReadoutMinimal performs the photometry analysis on PSPHOT.INPUT; we need to move
+    // around the pointers so PSPHOT.INPUT corresponds to the output image; previously, it was
+    // equivalent to the minuend image.
+    pmFPAfile *outputFile = psMetadataLookupPtr(&mdok, config->files, "PPSUB.OUTPUT");
+    pmReadout *outRO = pmFPAviewThisReadout(view, outputFile->fpa); // Readout with the sources
+
+    // XXX possibly rename this to PPSUB.RESID?
+    pmFPAfile *photFile = psMetadataLookupPtr(&mdok, config->files, "PSPHOT.RESID");
+    pmReadout *photRO = pmFPAviewThisReadout(view, photFile->fpa); // Readout with the sources
+    if (!photRO) {
+        pmCell *cell = pmFPAfileThisCell(config->files, view, "PSPHOT.RESID"); // Output cell
+        photRO = pmReadoutAlloc(cell); // Output readout: subtraction
+    }
+    photRO->image = psImageCopy(photRO->image, outRO->image, PS_TYPE_F32);
+    if (outRO->variance) {
+        photRO->variance = psImageCopy(photRO->variance, outRO->variance, PS_TYPE_F32);
+    } else {
+        psFree (photRO->variance);
+        photRO->variance = NULL;
+    }
+    if (outRO->mask) {
+        photRO->mask = psImageCopy(photRO->mask, outRO->mask, PS_TYPE_IMAGE_MASK);
+    } else {
+        psFree (photRO->mask);
+        photRO->mask = NULL;
+    }
+
+    // pmFPAfile *photFile = outputFile;
+    psMetadataAddPtr (config->files, PS_LIST_TAIL, "PSPHOT.INPUT", PS_DATA_UNKNOWN | PS_META_REPLACE, "psphot input : view on another pmFPAfile", photFile);
+
+    // old-style variance renormalization
+    if (!ppSubReadoutRenormPhot (config, recipe, photRO)) {
+        psError(PS_ERR_UNKNOWN, false, "failure in renormalization");
+        return false;
+    }
+
+    if (!psphotReadoutMinimal(config, view)) {
+        psWarning("Unable to perform photometry on subtracted image.");
+        psErrorStackPrint(stderr, "Error stack from photometry:");
+        psErrorClear();
+    }
+
+    if (stats) {
+        psArray *sources = psMetadataLookupPtr(NULL, photRO->analysis, "PSPHOT.SOURCES"); // Sources
+        psMetadataAddS32(stats, PS_LIST_TAIL, "NUM_SOURCES", 0, "Number of sources detected",
+                         sources ? sources->n : 0);
+        psMetadataAddF32(stats, PS_LIST_TAIL, "TIME_PHOT", 0, "Time to do photometry",
+                         psTimerClear("PPSUB_PHOT"));
+    }
+
+    pmFPAfileDropInternal(config->files, "PSPHOT.BACKMDL");
+    pmFPAfileDropInternal(config->files, "PSPHOT.BACKMDL.STDEV");
+
+    return true;
+}
+
+    // XXX not sure that this is still needed (only if psphotReadoutMinimal measures the background)
+    // if (!pmFPAfileDropInternal(config->files, "PSPHOT.BACKMDL") ||
+    //  !pmFPAfileDropInternal (config->files, "PSPHOT.BACKMDL.STDEV") ||
+    //  !pmFPAfileDropInternal (config->files, "PSPHOT.BACKGND")) {
+    //  psError(PS_ERR_UNKNOWN, false, "Unable to drop PSPHOT internal files.");
+    //  return false;
+    // }
+
+    // pmFPAfileActivate(config->files, false, "PSPHOT.INPUT");
+    // pmFPAfileActivate(config->files, false, "PSPHOT.LOAD.PSF");
+
+
+#ifdef TESTING
+    // Record data about sources: not everything gets into the output CMF files
+    {
+        pmReadout *photRO = pmFPAviewThisReadout(view, photFile->fpa); // Readout with the sources
+        psArray *sources = psMetadataLookupPtr(NULL, photRO->analysis, "PSPHOT.SOURCES"); // Sources
+        FILE *sourceFile = fopen("sources.dat", "w"); // File for sources
+        fprintf(sourceFile,
+                "# x y mag mag_err psf_chisq cr_nsigma ext_nsigma psf_qf flags m_x m_y m_xx m_xy m_yy\n");
+        for (int i = 0; i < sources->n; i++) {
+            pmSource *source = sources->data[i];
+            if (!source) {
+                continue;
+            }
+
+            float x, y;             // Position of source
+            float chi2;             // chi^2 for source
+            if (source->modelPSF) {
+                x = source->modelPSF->params->data.F32[PM_PAR_XPOS];
+                y = source->modelPSF->params->data.F32[PM_PAR_YPOS];
+                chi2 = source->modelPSF->chisq;
+            } else if (source->peak) {
+                x = source->peak->xf;
+                y = source->peak->yf;
+                chi2 = NAN;
+            } else {
+                psWarning("No position available for source.");
+                continue;
+            }
+
+            float xMoment = NAN, yMoment = NAN, xxMoment = NAN, xyMoment = NAN, yyMoment = NAN;
+            if (source->moments) {
+                xMoment = source->moments->Mx;
+                yMoment = source->moments->My;
+                xxMoment = source->moments->Mxx;
+                xyMoment = source->moments->Mxy;
+                yyMoment = source->moments->Myy;
+            }
+
+            fprintf(sourceFile, "%f %f %f %f %f %f %f %f %d %f %f %f %f %f\n",
+                    x, y, source->psfMag, source->errMag, chi2, source->crNsigma, source->extNsigma,
+                    source->pixWeight, source->mode, xMoment, yMoment, xxMoment, xyMoment, yyMoment);
+        }
+        fclose(sourceFile);
+    }
+#endif
Index: /branches/ipp-magic-v0/ppSub/src/ppSubReadoutRenorm.c
===================================================================
--- /branches/ipp-magic-v0/ppSub/src/ppSubReadoutRenorm.c	(revision 21375)
+++ /branches/ipp-magic-v0/ppSub/src/ppSubReadoutRenorm.c	(revision 21375)
@@ -0,0 +1,65 @@
+/** @file ppSubReadoutRenorm.c
+ *
+ *  @brief
+ *
+ *  @ingroup ppSub
+ *
+ *  @author IfA
+ *  @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2009-02-06 03:20:02 $
+ *  Copyright 2009 Institute for Astronomy, University of Hawaii
+ */
+
+#include "ppSub.h"
+
+bool ppSubReadoutRenormPixels (pmConfig *config, psMetadata *recipe, pmReadout *readout) {
+
+    bool mdok = false;
+
+    if (!psMetadataLookupBool(&mdok, recipe, "RENORM")) return true;
+
+    // Statistics for renormalisation
+    psStatsOptions renormMean = psStatsOptionFromString(psMetadataLookupStr(&mdok, recipe, "RENORM.MEAN"));
+    psStatsOptions renormStdev = psStatsOptionFromString(psMetadataLookupStr(&mdok, recipe, "RENORM.STDEV"));
+    if (renormMean == PS_STAT_NONE || renormStdev == PS_STAT_NONE) {
+        psError(PS_ERR_BAD_PARAMETER_VALUE, false,
+                "Unable to parse renormalisation statistics from recipe.");
+        return false;
+    }
+    psImageMaskType maskValue = pmConfigMaskGet("BLANK", config); // Bits to mask
+    if (!pmReadoutVarianceRenormPixels(readout, maskValue, renormMean, renormStdev, NULL)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to renormalise variances.");
+        return false;
+    }
+    return true;
+}
+
+/**
+ * old-style variance renormalization
+ */
+bool ppSubReadoutRenormPhot (pmConfig *config, psMetadata *recipe, pmReadout *readout) {
+
+    bool mdok = false;
+
+    if (!psMetadataLookupBool(&mdok, recipe, "RENORM")) return true;
+
+    // Statistics for renormalisation
+    psStatsOptions renormMean = psStatsOptionFromString(psMetadataLookupStr(&mdok, recipe, "RENORM.MEAN"));
+    psStatsOptions renormStdev = psStatsOptionFromString(psMetadataLookupStr(&mdok, recipe, "RENORM.STDEV"));
+    if (renormMean == PS_STAT_NONE || renormStdev == PS_STAT_NONE) {
+        psError(PS_ERR_BAD_PARAMETER_VALUE, false,
+                "Unable to parse renormalisation statistics from recipe.");
+        return false;
+    }
+
+    // other renorm parameters
+    int renormNum = psMetadataLookupS32(&mdok, recipe, "RENORM.NUM"); // Number of samples for renormalisation
+    float renormWidth = psMetadataLookupS32(&mdok, recipe, "RENORM.WIDTH"); // Width of Gaussian phot
+
+    psImageMaskType maskValue = pmConfigMaskGet("BLANK", config); // Bits to mask
+    if (!pmReadoutVarianceRenormPhot(readout, maskValue, renormNum, renormWidth, renormMean, renormStdev, NULL)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to renormalise variances.");
+        return false;
+    }
+    return true;
+}
Index: /branches/ipp-magic-v0/ppSub/src/ppSubVarianceFactors.c
===================================================================
--- /branches/ipp-magic-v0/ppSub/src/ppSubVarianceFactors.c	(revision 21375)
+++ /branches/ipp-magic-v0/ppSub/src/ppSubVarianceFactors.c	(revision 21375)
@@ -0,0 +1,117 @@
+/** @file ppSubVarianceFactors.c
+ *
+ *  @brief
+ *
+ *  @ingroup ppSub
+ *
+ *  @author IfA
+ *  @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2009-02-06 03:20:02 $
+ *  Copyright 2009 Institute for Astronomy, University of Hawaii
+ */
+
+#include "ppSub.h"
+
+/**
+ * Calculate the variance factor for the output image based on the input images
+ */
+bool ppSubVarianceFactors (pmConfig *config, psMetadata *stats, const pmFPAview *view) {
+
+    // Look up recipe values
+    psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSUB_RECIPE); // Recipe for ppSim
+    psAssert(recipe, "We checked this earlier, so it should be here.");
+
+    // convolved input images
+    pmCell *inCell = pmFPAfileThisCell(config->files, view, "PPSUB.INPUT.CONV"); // Input cell
+    pmCell *refCell = pmFPAfileThisCell(config->files, view, "PPSUB.REF.CONV"); // Reference cell
+
+    pmReadout *inConv = pmFPAfileThisReadout(config->files, view, "PPSUB.INPUT.CONV"); // Input cell
+    pmReadout *refConv = pmFPAfileThisReadout(config->files, view, "PPSUB.REF.CONV"); // Reference cell
+
+    pmReadout *outRO = pmFPAfileThisReadout(config->files, view, "PPSUB.OUTPUT"); // Output readout
+
+    // Update variance factors
+    // It's not possible to do this perfectly, since we're combining different images:
+    // vf_sum sigma_sum^2 = vf_1 * sigma_1^2 + vf_2 * sigma_2^2
+    // It's not possible to write vf_sum as a function of vf_1 and vf_2 with no reference to the sigmas.
+    // Instead, we're going to cheat.
+
+    bool mdok;                      // Status of MD lookup
+    pmSubtractionKernels *kernels = psMetadataLookupPtr(&mdok, outRO->analysis, PM_SUBTRACTION_ANALYSIS_KERNEL); // Kernels
+    if (!mdok) {
+        psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find subtraction kernels.");
+        return false;
+    }
+
+    float vfIn = psMetadataLookupF32(NULL, inCell->concepts, "CELL.VARFACTOR"); // Variance factor for input
+    if (!isfinite(vfIn)) {
+        vfIn = 1.0;
+    }
+    float vfRef = psMetadataLookupF32(NULL, refCell->concepts, "CELL.VARFACTOR"); // Variance factor for ref
+    if (!isfinite(vfRef)) {
+        vfRef = 1.0;
+    }
+
+    if (kernels->mode == PM_SUBTRACTION_MODE_1) {
+        vfIn *= pmSubtractionVarianceFactor(kernels, 0.0, 0.0, false);
+    }
+    if (kernels->mode == PM_SUBTRACTION_MODE_2) {
+        vfRef *= pmSubtractionVarianceFactor(kernels, 0.0, 0.0, false);
+    }
+    if (kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
+        vfIn *= pmSubtractionVarianceFactor(kernels, 0.0, 0.0, false);
+        vfRef *= pmSubtractionVarianceFactor(kernels, 0.0, 0.0, true);
+    }
+
+    psStats *vfStats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); // Statistics
+    psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS, 0); // Random number generator
+
+    // mask these values when examining the background
+    psImageMaskType maskVal = pmConfigMaskGet("MASK.VALUE", config); // Bits to mask going in to pmSubtractionMatch
+    psImageMaskType maskBad = pmConfigMaskGet("BLANK", config); // Bits to mask for bad pixels
+
+    // measure the mean of the convolved input variance image
+    if (!psImageBackground(vfStats, NULL, inConv->variance, inConv->mask, maskVal | maskBad, rng)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to measure mean variance for convolved input");
+        psFree(vfStats);
+        psFree(rng);
+        return false;
+    }
+    float inMeanVar = vfStats->robustMedian; // Mean variance of input
+
+    // measure the mean of the convolved reference variance image
+    if (!psImageBackground(vfStats, NULL, refConv->variance, refConv->mask, maskVal | maskBad, rng)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to measure mean variance for convolved reference");
+        psFree(vfStats);
+        psFree(rng);
+        return false;
+    }
+    float refMeanVar = vfStats->robustMedian; // Mean variance of reference
+
+    psFree(rng);
+    psFree(vfStats);
+
+    float vf = (vfIn * inMeanVar + vfRef * refMeanVar) / (inMeanVar + refMeanVar); // Resulting var factor
+    psMetadataItem *vfItem = psMetadataLookup(outRO->parent->concepts, "CELL.VARFACTOR");
+    vfItem->data.F32 = vf;
+
+    // Statistics on the matching
+    if (stats) {
+        psMetadataCopySingle(stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_MODE);
+        psMetadataCopySingle(stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_STAMPS);
+        psMetadataCopySingle(stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_DEV_MEAN);
+        psMetadataCopySingle(stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_DEV_RMS);
+        psMetadataCopySingle(stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_NORM);
+        psMetadataCopySingle(stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_BGDIFF);
+        psMetadataCopySingle(stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_MX);
+        psMetadataCopySingle(stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_MY);
+        psMetadataCopySingle(stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_MXX);
+        psMetadataCopySingle(stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_MXY);
+        psMetadataCopySingle(stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_MYY);
+
+        psMetadataAddF32(stats, PS_LIST_TAIL, "TIME_MATCH", 0, "Time to match PSFs",
+                         psTimerClear("PPSUB_MATCH"));
+    }
+
+    return true;
+}