Changeset 27840 for branches/simtest_nebulous_branches/ppSub/src

Timestamp:

May 3, 2010, 8:50:52 AM (16 years ago)

Author:

eugene

Message:

updates from trunk

Location:

branches/simtest_nebulous_branches

Files:

• . (modified) (1 prop)
• ppSub/src (modified) (1 prop)
• ppSub/src/Makefile.am (modified) (3 diffs)
• ppSub/src/ppSub.c (modified) (4 diffs)
• ppSub/src/ppSub.h (modified) (3 diffs)
• ppSub/src/ppSubArguments.c (modified) (4 diffs)
• ppSub/src/ppSubBackground.c (modified) (2 diffs)
• ppSub/src/ppSubCamera.c (modified) (22 diffs)
• ppSub/src/ppSubConvolve.c (copied) (copied from trunk/ppSub/src/ppSubConvolve.c )
• ppSub/src/ppSubData.c (modified) (3 diffs)
• ppSub/src/ppSubDefineOutput.c (modified) (1 diff)
• ppSub/src/ppSubErrorCodes.dat (modified) (1 diff)
• ppSub/src/ppSubExit.c (copied) (copied from trunk/ppSub/src/ppSubExit.c )
• ppSub/src/ppSubFiles.c (modified) (2 diffs)
• ppSub/src/ppSubLoop.c (modified) (11 diffs)
• ppSub/src/ppSubMakePSF.c (modified) (4 diffs)
• ppSub/src/ppSubMatchPSFs.c (modified) (9 diffs)
• ppSub/src/ppSubReadoutInverse.c (modified) (2 diffs)
• ppSub/src/ppSubReadoutJpeg.c (modified) (3 diffs)
• ppSub/src/ppSubReadoutPhotometry.c (modified) (5 diffs)
• ppSub/src/ppSubReadoutStats.c (modified) (1 diff)
• ppSub/src/ppSubReadoutSubtract.c (modified) (2 diffs)
• ppSub/src/ppSubSetMasks.c (modified) (4 diffs)
• ppSub/src/ppSubThreshold.c (modified) (5 diffs)
• ppSub/src/ppSubVarianceRescale.c (modified) (1 diff)

branches/simtest_nebulous_branches/ppSub/src

@@ -13 +13 @@
 ppSubErrorCodes.c
 ppSubVersionDefinitions.h
+ppSubConvolve

@@ -13 +13 @@
 ppSubErrorCodes.c
 ppSubVersionDefinitions.h
+ppSubConvolve

branches/simtest_nebulous_branches/ppSub/src/Makefile.am

@@ -1 @@
-bin_PROGRAMS = ppSub ppSubKernel
+bin_PROGRAMS = ppSub ppSubKernel ppSubConvolve
 
 if HAVE_SVNVERSION
@@ -34 +34 @@
         ppSubData.c                     \
         ppSubErrorCodes.c               \
+        ppSubExit.c                     \
         ppSubFiles.c                    \
         ppSubLoop.c                     \
@@ -53 +54 @@
 ppSubKernel_SOURCES =           \
         ppSubKernel.c
+
+ppSubConvolve_CPPFLAGS = $(PSLIB_CFLAGS) $(PSMODULE_CFLAGS) $(PPSTATS_CFLAGS) $(PSPHOT_CFLAGS) $(PPSUB_CFLAGS)
+ppSubConvolve_LDFLAGS  = $(PSLIB_LIBS)   $(PSMODULE_LIBS)   $(PPSTATS_LIBS)   $(PSPHOT_LIBS)   $(PPSUB_LIBS)
+
+ppSubConvolve_SOURCES =         \
+        ppSubConvolve.c         \
+        ppSubExit.c             \
+        ppSubVersion.c
 
 noinst_HEADERS = \

branches/simtest_nebulous_branches/ppSub/src/ppSub.c

@@ -24 +24 @@
 int main(int argc, char *argv[])
 {
-    psExit exitValue = PS_EXIT_SUCCESS; // Exit value
+    psLibInit(NULL);
     psTimerStart("ppSub");
-    psLibInit(NULL);
+
+    pmErrorRegister();
+    ppSubErrorRegister();
+    psphotErrorRegister();
 
     ppSubData *data = NULL;             // Processing data
     pmConfig *config = pmConfigRead(&argc, argv, PPSUB_RECIPE); // Configuration
     if (!config) {
-        psErrorStackPrint(stderr, "Error reading configuration.");
-        exitValue = PS_EXIT_CONFIG_ERROR;
         goto die;
     }
@@ -39 +40 @@
 
     if (!pmModelClassInit()) {
-        psErrorStackPrint(stderr, "Error initialising model classes.\n");
-        exitValue = PS_EXIT_PROG_ERROR;
+        psError(PPSUB_ERR_PROG, false, "Unable to initialise model classes.");
         psFree(config);
         goto die;
@@ -46 +46 @@
 
     if (!psphotInit()) {
-        psErrorStackPrint(stderr, "Error initialising psphot.\n");
-        exitValue = PS_EXIT_PROG_ERROR;
+        psError(PPSUB_ERR_PROG, false, "Error initialising psphot.");
         psFree(config);
         goto die;
@@ -55 +54 @@
 
     if (!ppSubArguments(argc, argv, data)) {
-        psErrorStackPrint(stderr, "Error reading arguments.\n");
-        exitValue = PS_EXIT_CONFIG_ERROR;
+        psError(psErrorCodeLast(), false, "Error reading arguments.");
         goto die;
     }
 
     if (!ppSubCamera(data)) {
-        psErrorStackPrint(stderr, "Error setting up camera.\n");
-        exitValue = PS_EXIT_CONFIG_ERROR;
         goto die;
     }
 
     if (!ppSubLoop(data)) {
-        psErrorStackPrint(stderr, "Error performing subtraction.\n");
-        exitValue = PS_EXIT_SYS_ERROR;
         goto die;
     }
 
  die:
-    psTrace("ppSub", 1, "Finished at %f sec\n", psTimerMark("ppSub"));
-    psTimerStop();
+    {
+        psExit exitValue = ppSubExitCode(PS_EXIT_SUCCESS); // Exit code
 
-    psFree(data);
+        if (data && data->stats && data->statsFile) {
+            psString stats = psMetadataConfigFormat(data->stats); // Statistics to output
+            if (!stats || strlen(stats) == 0) {
+                psError(PPSUB_ERR_IO, false, "Unable to format statistics file");
+            } else if (fprintf(data->statsFile, "%s", stats) != strlen(stats)) {
+                psError(PPSUB_ERR_IO, true, "Unable to write statistics file");
+            }
+            psFree(stats);
+            if (fclose(data->statsFile) == EOF) {
+                psError(PPSUB_ERR_IO, true, "Unable to close statistics file");
+            }
+            data->statsFile = NULL;
+            pmConfigRunFilenameAddWrite(data->config, "STATS", data->statsName);
+            exitValue = ppSubExitCode(exitValue);
+        }
 
-    pmVisualClose(); //close plot windows, if -visual is set
-    pmModelClassCleanup();
-    pmConfigDone();
-    psLibFinalize();
+        if (config && !ppSubFilesIterateUp(config, PPSUB_FILES_ALL)) {
+            psError(psErrorCodeLast(), false, "Unable to close files.");
+            exitValue = ppSubExitCode(exitValue);
+            pmFPAfileFreeSetStrict(false);
+        }
 
-    exit(exitValue);
+        if (data) {
+            psString dump_file = psMetadataLookupStr(NULL, data->config->arguments, "-dumpconfig");
+            if (dump_file) {
+                if (!pmConfigDump(data->config, dump_file)) {
+                    psError(PPSUB_ERR_IO, false, "Unable to dump configuration.");
+                    exitValue = ppSubExitCode(exitValue);
+                }
+            }
+            psFree(data);
+        }
+
+        psTrace("ppSub", 1, "Finished at %f sec\n", psTimerMark("ppSub"));
+        psTimerStop();
+
+        pmVisualClose(); //close plot windows, if -visual is set
+        pmModelClassCleanup();
+        pmConfigDone();
+        psLibFinalize();
+
+        exitValue = ppSubExitCode(exitValue);
+        exit(exitValue);
+    }
 }

branches/simtest_nebulous_branches/ppSub/src/ppSub.h

@@ -45 +45 @@
     bool photometry;                    // Perform photometry?
     bool inverse;                       // Output inverse subtraction as well?
+    bool saveInConv;                    // Save convolved input?
+    bool saveRefConv;                   // Save convolved reference?
     psString stamps;                    // Stamps file
     pmPSF *psf;                         // Point Spread Function
@@ -152 +154 @@
     );
 
+/// Generate JPEG images
+bool ppSubResidualSampleJpeg(pmConfig *config);
+
 /// Generate inverse subtraction
 bool ppSubReadoutInverse(pmConfig *config // Configuration
@@ -160 +165 @@
 bool psMetadataCopySingle(psMetadata *target, psMetadata *source, const char *name);
 
+bool ppSubCopyPSF (pmFPAfile *output, pmFPAfile *input, pmFPAview *view);
+
+/// Return appropriate exit code
+psExit ppSubExitCode(psExit exitValue   // Current exit value
+    );
+
 ///@}
 #endif

branches/simtest_nebulous_branches/ppSub/src/ppSubArguments.c

@@ -86 +86 @@
     psMetadataAddS32(arguments, PS_LIST_TAIL, "-convolve", 0, "Image to convolve [1 or 2]", 0);
     psMetadataAddBool(arguments, PS_LIST_TAIL, "-photometry", 0, "Perform photometry?", NULL);
-    psMetadataAddBool(arguments, PS_LIST_TAIL, "-inverse", 0, "Generate inverse subtractions?", NULL);
+    psMetadataAddF32(arguments, PS_LIST_TAIL, "-zp", 0, "Zero point for photometry", NAN);
+    psMetadataAddBool(arguments, PS_LIST_TAIL, "-inverse", 0, "Generate inverse subtractions?", false);
+    psMetadataAddBool(arguments, PS_LIST_TAIL, "-save-inconv", 0, "Save input convolved images?", false);
+    psMetadataAddBool(arguments, PS_LIST_TAIL, "-save-refconv", 0, "Save reference convolved images?", false);
     psMetadataAddBool(arguments, PS_LIST_TAIL, "-visual", 0, "Show diagnostic plots", NULL);
 
@@ -140 +143 @@
     if (data->statsName && strlen(data->statsName) > 0) {
         psString resolved = pmConfigConvertFilename(data->statsName, config, true, true); // Resolved filename
+        if (!resolved) {
+            psError(psErrorCodeLast(), false, "Unable to resolve statistics file: %s", data->statsName);
+            return false;
+        }
         data->statsFile = fopen(resolved, "w");
         if (!data->statsFile) {
@@ -149 +156 @@
     }
 
+    data->saveInConv = psMetadataLookupBool(NULL, arguments, "-save-inconv");
+    data->saveRefConv = psMetadataLookupBool(NULL, arguments, "-save-refconv");
+
     if (psMetadataLookupBool(NULL, arguments, "-visual")) {
         pmVisualSetVisual(true);
@@ -156 +166 @@
     if (threads > 0) {
         if (!psThreadPoolInit(threads)) {
-            psError(PS_ERR_UNKNOWN, false, "Unable to setup %d threads", threads);
+            psError(psErrorCodeLast(), false, "Unable to setup %d threads", threads);
             return false;
         }

branches/simtest_nebulous_branches/ppSub/src/ppSubBackground.c

@@ -42 +42 @@
     if (!modelRO) {
         // Create the background model
-        if (!psphotModelBackground(config, view, "PPSUB.OUTPUT")) {
-            psError(PS_ERR_UNKNOWN, false, "Unable to model background");
+        if (!psphotModelBackgroundReadoutFileIndex(config, view, "PPSUB.OUTPUT", 0)) {
+            psError(psErrorCodeLast(), false, "Unable to model background");
             psFree(view);
             return false;
@@ -50 +50 @@
         modelRO = pmFPAfileThisReadout(config->files, view, "PSPHOT.BACKMDL");
         if (!modelRO) {
-            psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find background model");
+            psError(psErrorCodeLast(), false, "Unable to find background model");
             psFree(view);
             return false;

branches/simtest_nebulous_branches/ppSub/src/ppSubCamera.c

@@ -23 +23 @@
 
 // Define an input file
-static pmFPAfile *defineInputFile(pmConfig *config,// Configuration
+static pmFPAfile *defineInputFile(bool *success,
+                                  pmConfig *config,// Configuration
                                   pmFPAfile *bind,    // File to which to bind, or NULL
                                   char *filerule,     // Name of file rule
@@ -32 +33 @@
     bool status;
 
-    // look for the file on the RUN metadata
-    pmFPAfile *file = pmFPAfileDefineFromRun(&status, bind, config, filerule); // File to return
+    *success = false;
+
+    pmFPAfile *file = NULL;
+
+    // look for the file on the argument list
+    if (bind) {
+        file = pmFPAfileBindFromArgs(&status, bind, config, filerule, argname);
+    } else {
+        file = pmFPAfileDefineFromArgs(&status, config, filerule, argname);
+    }
+
     if (!status) {
-        psError(PS_ERR_UNKNOWN, false, "Failed to load file definition for %s", filerule);
-        return NULL;
+        psError(psErrorCodeLast(), false, "Failed to load file definition for %s", filerule);
+        return false;
     }
     if (!file) {
-        // look for the file on the argument list
-        if (bind) {
-            file = pmFPAfileBindFromArgs(&status, bind, config, filerule, argname);
-        } else {
-            file = pmFPAfileDefineFromArgs(&status, config, filerule, argname);
-        }
+        // look for the file on the RUN metadata
+        file = pmFPAfileDefineFromRun(&status, bind, config, filerule); // File to return
         if (!status) {
-            psError(PS_ERR_UNKNOWN, false, "Failed to load file definition for %s", filerule);
-            return false;
+            psError(psErrorCodeLast(), false, "Failed to load file definition for %s", filerule);
+            return NULL;
         }
     }
 
     if (!file) {
+        // no file defined
+        *success = true;
         return NULL;
     }
 
     if (file->type != fileType) {
-        psError(PS_ERR_UNKNOWN, true, "%s is not of type %s", filerule, pmFPAfileStringFromType(fileType));
-        return NULL;
-    }
-
+        psError(PPSUB_ERR_CONFIG, true, "%s is not of type %s", filerule, pmFPAfileStringFromType(fileType));
+        return NULL;
+    }
+
+    *success = true;
     return file;
 }
@@ -75 +84 @@
         pmFPAfileDefineOutput(config, template->fpa, filerule);
     if (!file) {
-        psError(PS_ERR_IO, false, _("Unable to generate output file from %s"), filerule);
+        psError(psErrorCodeLast(), false, _("Unable to generate output file from %s"), filerule);
         return NULL;
     }
     if (file->type != fileType) {
-        psError(PS_ERR_IO, true, "%s is not of type %s", filerule, pmFPAfileStringFromType(fileType));
+        psError(PPSUB_ERR_CONFIG, true, "%s is not of type %s", filerule, pmFPAfileStringFromType(fileType));
         return NULL;
     }
@@ -92 +101 @@
                                  pmFPAfile *bind,    // File to which to bind, or NULL
                                  char *filerule,     // Name of file rule
+                                 const char *argname,   // Argument name for file
                                  pmFPAfileType fileType // Type of file
     )
@@ -98 +108 @@
 
     // look for the file on the RUN metadata
-    pmFPAfile *file = pmFPAfileDefineFromRun(&status, bind, config, filerule); // File to return
+    pmFPAfile *file = pmFPAfileDefineFromRun(&status, NULL, config, filerule); // File to return
     if (!status) {
-        psError(PS_ERR_UNKNOWN, false, "Failed to load file definition for %s", filerule);
+        psError(psErrorCodeLast(), false, "Failed to load file definition for %s", filerule);
+        return NULL;
+    }
+    file = pmFPAfileBindFromArgs(&status, bind, config, filerule, argname);
+    if (!status) {
+        psError(psErrorCodeLast(), false, "Failed to load file definition for %s", filerule);
         return NULL;
     }
@@ -112 +127 @@
         file = pmFPAfileDefineOutput(config, bind ? bind->fpa : NULL, filerule);
         if (!status) {
-            psError(PS_ERR_UNKNOWN, false, "Failed to load file definition for %s", filerule);
+            psError(psErrorCodeLast(), false, "Failed to load file definition for %s", filerule);
             return false;
         }
@@ -126 +141 @@
 
     if (file->type != fileType) {
-        psError(PS_ERR_UNKNOWN, true, "%s is not of type %s", filerule, pmFPAfileStringFromType(fileType));
+        psError(PPSUB_ERR_CONFIG, true, "%s is not of type %s", filerule, pmFPAfileStringFromType(fileType));
         return NULL;
     }
@@ -136 +151 @@
 bool ppSubCamera(ppSubData *data)
 {
+    bool success = true;
+
     psAssert(data, "Require processing data");
     pmConfig *config = data->config;
@@ -141 +158 @@
 
     // Input image
-    pmFPAfile *input = defineInputFile(config, NULL, "PPSUB.INPUT", "INPUT", PM_FPA_FILE_IMAGE);
-    if (!input) {
-        psError(PS_ERR_IO, false, "Failed to build FPA from PPSUB.INPUT");
-        return false;
-    }
-    defineInputFile(config, input, "PPSUB.INPUT.MASK", "INPUT.MASK", PM_FPA_FILE_MASK);
-    pmFPAfile *inVar = defineInputFile(config, input, "PPSUB.INPUT.VARIANCE", "INPUT.VARIANCE",
-                                       PM_FPA_FILE_VARIANCE);
-    defineInputFile(config, NULL, "PPSUB.INPUT.SOURCES", "INPUT.SOURCES", PM_FPA_FILE_CMF);
+    pmFPAfile *input = defineInputFile(&success, config, NULL, "PPSUB.INPUT", "INPUT", PM_FPA_FILE_IMAGE);
+    if (!success) {
+        psError(psErrorCodeLast(), false, "Failed to build FPA from PPSUB.INPUT");
+        return false;
+    }
+
+    defineInputFile(&success, config, input, "PPSUB.INPUT.MASK", "INPUT.MASK", PM_FPA_FILE_MASK);
+    if (!success) {
+        psError(psErrorCodeLast(), false, "Failed to build FPA from PPSUB.INPUT.MASK");
+        return false;
+    }
+
+    pmFPAfile *inVar = defineInputFile(&success, config, input, "PPSUB.INPUT.VARIANCE", "INPUT.VARIANCE", PM_FPA_FILE_VARIANCE);
+    if (!success) {
+        psError(psErrorCodeLast(), false, "Failed to build FPA from PPSUB.INPUT.VARIANCE");
+        return false;
+    }
+
+    defineInputFile(&success, config, NULL, "PPSUB.INPUT.SOURCES", "INPUT.SOURCES", PM_FPA_FILE_CMF);
+    if (!success) {
+        psError(psErrorCodeLast(), false, "Failed to build FPA from PPSUB.INPUT.SOURCES");
+        return false;
+    }
 
     // Reference image
-    pmFPAfile *ref = defineInputFile(config, NULL, "PPSUB.REF", "REF", PM_FPA_FILE_IMAGE);
-    if (!ref) {
-        psError(PS_ERR_IO, false, "Failed to build FPA from PPSUB.REF");
-        return false;
-    }
-    defineInputFile(config, ref, "PPSUB.REF.MASK", "REF.MASK", PM_FPA_FILE_MASK);
-    pmFPAfile *refVar = defineInputFile(config, ref, "PPSUB.REF.VARIANCE", "REF.VARIANCE",
-                                        PM_FPA_FILE_VARIANCE);
-    defineInputFile(config, NULL, "PPSUB.REF.SOURCES", "REF.SOURCES", PM_FPA_FILE_CMF);
-
+    pmFPAfile *ref = defineInputFile(&success, config, NULL, "PPSUB.REF", "REF", PM_FPA_FILE_IMAGE);
+    if (!success) {
+        psError(psErrorCodeLast(), false, "Failed to build FPA from PPSUB.REF");
+        return false;
+    }
+
+    defineInputFile(&success, config, ref, "PPSUB.REF.MASK", "REF.MASK", PM_FPA_FILE_MASK);
+    if (!success) {
+        psError(psErrorCodeLast(), false, "Failed to build FPA from PPSUB.REF.MASK");
+        return false;
+    }
+
+    pmFPAfile *refVar = defineInputFile(&success, config, ref, "PPSUB.REF.VARIANCE", "REF.VARIANCE", PM_FPA_FILE_VARIANCE);
+    if (!success) {
+        psError(psErrorCodeLast(), false, "Failed to build FPA from PPSUB.REF.VARIANCE");
+        return false;
+    }
+
+    defineInputFile(&success, config, NULL, "PPSUB.REF.SOURCES", "REF.SOURCES", PM_FPA_FILE_CMF);
+    if (!success) {
+        psError(psErrorCodeLast(), false, "Failed to build FPA from PPSUB.REF.SOURCES");
+        return false;
+    }
 
     // Now that the camera has been determined, we can read the recipe
     psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSUB_RECIPE); // Recipe for ppSim
     if (!recipe) {
-        psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to find recipe %s", PPSUB_RECIPE);
+        psError(PPSUB_ERR_CONFIG, false, "Unable to find recipe %s", PPSUB_RECIPE);
         return false;
     }
@@ -180 +224 @@
     data->inverse = psMetadataLookupBool(NULL, recipe, "INVERSE");
     data->photometry = psMetadataLookupBool(NULL, recipe, "PHOTOMETRY");
-
-    bool mdok;                          // Status of MD lookup
-    bool saveConv = psMetadataLookupBool(&mdok, recipe, "SAVE.CONVOLVED"); // Save convolved images?
 
     // Convolved input image
@@ -189 +230 @@
                                              PM_FPA_FILE_MASK);
     if (!inConvImage || !inConvMask) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to define output files");
-        return false;
-    }
-    if (saveConv) {
-        inConvImage->save = true;
-        inConvMask->save = true;
-    }
+        psError(psErrorCodeLast(), false, "Unable to define output files");
+        return false;
+    }
+    inConvImage->save = data->saveInConv;
+    inConvMask->save = data->saveInConv;
     if (inVar) {
         pmFPAfile *inConvVar = defineOutputFile(config, inConvImage, false, "PPSUB.INPUT.CONV.VARIANCE",
                                                 PM_FPA_FILE_VARIANCE);
         if (!inConvVar) {
-            psError(PS_ERR_UNKNOWN, false, "Unable to define output files");
-            return false;
-        }
-        if (saveConv) {
-            inConvVar->save = true;
-        }
+            psError(psErrorCodeLast(), false, "Unable to define output files");
+            return false;
+        }
+        inConvVar->save = data->saveInConv;
     }
 
@@ -213 +250 @@
                                               PM_FPA_FILE_MASK);
     if (!refConvImage || !refConvMask) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to define output files");
-        return false;
-    }
-    if (saveConv) {
-        refConvImage->save = true;
-        refConvMask->save = true;
-    }
+        psError(psErrorCodeLast(), false, "Unable to define output files");
+        return false;
+    }
+    refConvImage->save = data->saveRefConv;
+    refConvMask->save = data->saveRefConv;
     if (refVar) {
         pmFPAfile *refConvVar = defineOutputFile(config, refConvImage, false, "PPSUB.REF.CONV.VARIANCE",
                                                  PM_FPA_FILE_VARIANCE);
         if (!refConvVar) {
-            psError(PS_ERR_UNKNOWN, false, "Unable to define output files");
-            return false;
-        }
-        if (saveConv) {
-            refConvVar->save = true;
-        }
+            psError(psErrorCodeLast(), false, "Unable to define output files");
+            return false;
+        }
+        refConvVar->save = data->saveRefConv;
     }
 
@@ -236 +269 @@
     pmFPAfile *outMask = defineOutputFile(config, output, false, "PPSUB.OUTPUT.MASK", PM_FPA_FILE_MASK);
     if (!output || !outMask) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to define output files");
+        psError(psErrorCodeLast(), false, "Unable to define output files");
         return false;
     }
@@ -245 +278 @@
                                              PM_FPA_FILE_VARIANCE);
         if (!outVar) {
-            psError(PS_ERR_UNKNOWN, false, "Unable to define output files");
+            psError(psErrorCodeLast(), false, "Unable to define output files");
             return false;
         }
@@ -258 +291 @@
                                               PM_FPA_FILE_MASK);
         if (!inverse || !invMask) {
-            psError(PS_ERR_UNKNOWN, false, "Unable to define output files");
+            psError(psErrorCodeLast(), false, "Unable to define output files");
             return false;
         }
@@ -267 +300 @@
                                                  PM_FPA_FILE_VARIANCE);
             if (!invVar) {
-                psError(PS_ERR_UNKNOWN, false, "Unable to define output files");
+                psError(psErrorCodeLast(), false, "Unable to define output files");
                 return false;
             }
@@ -278 +311 @@
     pmFPAfile *jpeg1 = pmFPAfileDefineOutput(config, NULL, "PPSUB.OUTPUT.JPEG1");
     if (!jpeg1) {
-        psError(PS_ERR_IO, false, _("Unable to generate output file from PPSUB.OUTPUT.JPEG1"));
+        psError(psErrorCodeLast(), false, _("Unable to generate output file from PPSUB.OUTPUT.JPEG1"));
         return false;
     }
     if (jpeg1->type != PM_FPA_FILE_JPEG) {
-        psError(PS_ERR_IO, true, "PPSUB.OUTPUT.JPEG1 is not of type JPEG");
+        psError(psErrorCodeLast(), true, "PPSUB.OUTPUT.JPEG1 is not of type JPEG");
         return false;
     }
@@ -288 +321 @@
     pmFPAfile *jpeg2 = pmFPAfileDefineOutput(config, NULL, "PPSUB.OUTPUT.JPEG2");
     if (!jpeg2) {
-        psError(PS_ERR_IO, false, _("Unable to generate output file from PPSUB.OUTPUT.JPEG2"));
+        psError(psErrorCodeLast(), false, _("Unable to generate output file from PPSUB.OUTPUT.JPEG2"));
         return false;
     }
     if (jpeg2->type != PM_FPA_FILE_JPEG) {
-        psError(PS_ERR_IO, true, "PPSUB.OUTPUT.JPEG2 is not of type JPEG");
+        psError(psErrorCodeLast(), true, "PPSUB.OUTPUT.JPEG2 is not of type JPEG");
         return false;
     }
     jpeg2->save = true;
 
+    // Output residual JPEG
+    pmFPAfile *jpeg3 = pmFPAfileDefineOutput(config, NULL, "PPSUB.OUTPUT.RESID.JPEG");
+    if (!jpeg3) {
+        psError(psErrorCodeLast(), false, _("Unable to generate output file from PPSUB.OUTPUT.RESID.JPEG"));
+        return false;
+    }
+    if (jpeg3->type != PM_FPA_FILE_JPEG) {
+        psError(psErrorCodeLast(), true, "PPSUB.OUTPUT.RESID.JPEG is not of type JPEG");
+        return false;
+    }
+    jpeg3->save = true;
+
     // Output subtraction kernel
-    pmFPAfile *kernel = defineCalcFile(config, output, "PPSUB.OUTPUT.KERNELS", PM_FPA_FILE_SUBKERNEL);
+    pmFPAfile *kernel = defineCalcFile(config, output, "PPSUB.OUTPUT.KERNELS", "KERNEL",
+                                       PM_FPA_FILE_SUBKERNEL);
     if (!kernel) {
         psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to define file PPSUB.OUTPUT.KERNELS");
@@ -305 +351 @@
 
     // psPhot input
-    if (data->photometry) {
+    if (data->photometry || 1) {
         psphotModelClassInit();        // load implementation-specific models
 
         pmFPAfile *psphot = pmFPAfileDefineFromFPA(config, output->fpa, 1, 1, "PSPHOT.INPUT");
         if (!psphot) {
-            psError(PS_ERR_IO, false, "Failed to build FPA from PSPHOT.INPUT");
+            psError(psErrorCodeLast(), false, "Failed to build FPA from PSPHOT.INPUT");
             return false;
         }
         if (psphot->type != PM_FPA_FILE_IMAGE) {
-            psError(PS_ERR_IO, true, "PSPHOT.INPUT is not of type IMAGE");
-            return false;
-        }
+            psError(psErrorCodeLast(), true, "PSPHOT.INPUT is not of type IMAGE");
+            return false;
+        }
+        // specify the number of psphot input images
+        psMetadataAddS32 (config->arguments, PS_LIST_TAIL, "PSPHOT.INPUT.NUM", PS_META_REPLACE, "number of inputs", 1);
         pmFPAfileActivate(config->files, false, "PSPHOT.INPUT");
 
-        // Internal-ish file for getting the PSF from the minuend
-        pmFPAfile *psf = pmFPAfileDefineOutputFromFile(config, psphot, "PSPHOT.PSF.LOAD");
+        // Internal file for getting the PSF from the minuend
+        pmFPAfile *psf = pmFPAfileDefineFromFPA(config, output->fpa, 1, 1, "PSPHOT.PSF.LOAD");
         if (!psf) {
-            psError(PS_ERR_IO, false, "Failed to build FPA from PSPHOT.PSF.LOAD");
+            psError(psErrorCodeLast(), false, "Failed to build FPA from PSPHOT.PSF.LOAD");
             return false;
         }
         if (psf->type != PM_FPA_FILE_PSF) {
-            psError(PS_ERR_IO, true, "PSPHOT.PSF.LOAD is not of type PSF");
+            psError(psErrorCodeLast(), true, "PSPHOT.PSF.LOAD is not of type PSF");
             return false;
         }
@@ -332 +380 @@
 
         if (!psphotDefineFiles(config, psphot)) {
-            psError(PS_ERR_UNKNOWN, false, "Unable to set up psphot files.");
+            psError(psErrorCodeLast(), false, "Unable to set up psphot files.");
             return false;
         }
@@ -343 +391 @@
                                                  PM_FPA_FILE_CMF);
         if (!outSources) {
-            psError(PS_ERR_UNKNOWN, false, "Unable to set up output source file.");
+            psError(psErrorCodeLast(), false, "Unable to set up output source file.");
             return false;
         }
@@ -352 +400 @@
                                                      PM_FPA_FILE_CMF);
             if (!invSources) {
-                psError(PS_ERR_UNKNOWN, false, "Unable to set up inverse source file.");
+                psError(psErrorCodeLast(), false, "Unable to set up inverse source file.");
                 return false;
             }

branches/simtest_nebulous_branches/ppSub/src/ppSubData.c

@@ -13 +13 @@
 static void subDataFree(ppSubData *data)
 {
-    if (data->statsFile) {
-        psString stats = psMetadataConfigFormat(data->stats); // Statistics to output
-        if (!stats || strlen(stats) == 0) {
-            psWarning("Unable to generate statistics file.");
-        } else {
-            fprintf(data->statsFile, "%s", stats);
-        }
-        psFree(stats);
-        fclose(data->statsFile);
-        pmConfigRunFilenameAddWrite(data->config, "STATS", data->statsName);
-    }
+    psAssert(!data->statsFile, "Statistics file still open.");
     psFree(data->statsName);
     psFree(data->stamps);
@@ -29 +19 @@
     psFree(data->stats);
 
-    psString dump_file = psMetadataLookupStr(NULL, data->config->arguments, "-dumpconfig");
-    if (dump_file) {
-        pmConfigDump(data->config, dump_file);
-    }
     psFree(data->config);
 
@@ -47 +33 @@
     data->photometry = false;
     data->inverse = false;
+    data->saveInConv = false;
+    data->saveRefConv = false;
     data->stamps = NULL;
     data->psf = NULL;

branches/simtest_nebulous_branches/ppSub/src/ppSubDefineOutput.c

@@ -51 +51 @@
     bool mdok;                          // Status of MD lookup
     psMetadata *analysis = inConv->analysis; // Analysis metadata with kernel information
-    pmSubtractionKernels *kernels = psMetadataLookupPtr(&mdok, analysis,
-                                                        PM_SUBTRACTION_ANALYSIS_KERNEL); // Subtraction kernel
+    pmHDU *hdu = pmHDUFromCell(inConv->parent);
+    pmSubtractionKernels *kernels = psMetadataLookupPtr(&mdok, analysis, PM_SUBTRACTION_ANALYSIS_KERNEL); // Subtraction kernel
+
     if (!kernels) {
+        hdu = pmHDUFromCell(refConv->parent);
         analysis = refConv->analysis;
         kernels = psMetadataLookupPtr(&mdok, analysis, PM_SUBTRACTION_ANALYSIS_KERNEL);
     }
     if (!kernels) {
-        psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find SUBTRACTION.KERNEL");
+        psError(PPSUB_ERR_UNKNOWN, true, "Unable to find SUBTRACTION.KERNEL");
         psFree(outRO);
         return false;
     }
-    psMetadataAddStr(outHDU->header, PS_LIST_TAIL, "PPSUB.KERNEL", 0, "Subtraction kernel",
-                     kernels->description);
+    psAssert (hdu, "unable to find HDU");
+    psMetadataAddStr(outHDU->header, PS_LIST_TAIL, "PPSUB.KERNEL", 0, "Subtraction kernel", kernels->description);
+    outHDU->header = psMetadataCopy(outHDU->header, hdu->header);
 
     // Add additional data to the header

branches/simtest_nebulous_branches/ppSub/src/ppSubErrorCodes.dat

@@ -10 +10 @@
 DATA                    Problem in data values
 NO_OVERLAP              No overlap between input and skycell
+PROG                    Programming error

branches/simtest_nebulous_branches/ppSub/src/ppSubFiles.c

@@ -23 +23 @@
 // Subtraction files
 static const char *subFiles[] = { "PPSUB.OUTPUT", "PPSUB.OUTPUT.MASK", "PPSUB.OUTPUT.VARIANCE",
-                                  "PPSUB.OUTPUT.JPEG1", "PPSUB.OUTPUT.JPEG2",
+                                  "PPSUB.OUTPUT.JPEG1", "PPSUB.OUTPUT.JPEG2", "PPSUB.OUTPUT.RESID.JPEG",
                                   NULL };
 
@@ -183 +183 @@
     }
 
+    psFree(view);
+
+    // Disable writing of data now that we've written things out
+    if (files & PPSUB_FILES_CONV) {
+        pmFPAview *view = ppSubViewReadout(); // View to readout
+        {
+            pmReadout *ro = pmFPAfileThisReadout(config->files, view, "PPSUB.INPUT.CONV");
+            if (ro) {
+                ro->data_exists = ro->parent->data_exists = ro->parent->parent->data_exists = false;
+            }
+        }
+        {
+            pmReadout *ro = pmFPAfileThisReadout(config->files, view, "PPSUB.REF.CONV");
+            if (ro) {
+                ro->data_exists = ro->parent->data_exists = ro->parent->parent->data_exists = false;
+            }
+        }
+        psFree(view);
+    }
+    if (files & PPSUB_FILES_SUB) {
+        pmFPAview *view = ppSubViewReadout(); // View to readout
+        pmReadout *ro = pmFPAfileThisReadout(config->files, view, "PPSUB.OUTPUT");
+        if (ro) {
+            ro->data_exists = ro->parent->data_exists = ro->parent->parent->data_exists = false;
+        }
+        psFree(view);
+    }
+    // Disable writing of data now that we've written things out
+    if (files & PPSUB_FILES_INV) {
+        pmFPAview *view = ppSubViewReadout(); // View to readout
+        pmReadout *ro = pmFPAfileThisReadout(config->files, view, "PPSUB.INVERSE");
+        if (ro) {
+            ro->data_exists = ro->parent->data_exists = ro->parent->parent->data_exists = false;
+        }
+        psFree(view);
+    }
+
     ppSubFilesActivate(config, PPSUB_FILES_ALL, false);
 

branches/simtest_nebulous_branches/ppSub/src/ppSubLoop.c

@@ -27 +27 @@
 
     pmConfigCamerasCull(config, NULL);
-    pmConfigRecipesCull(config, "PPSUB,PPSTATS,PSPHOT,MASKS,JPEG");
+    pmConfigRecipesCull(config, "PPSUB,PPSTATS,PSPHOT,PSASTRO,MASKS,JPEG");
 
 
@@ -43 +43 @@
 
     if (!ppSubSetMasks(config)) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to set masks.");
+        psError(psErrorCodeLast(), false, "Unable to set masks.");
         return false;
     }
 
     if (!ppSubMatchPSFs(data)) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to match PSFs.");
+        psError(psErrorCodeLast(), false, "Unable to match PSFs.");
         return false;
     }
@@ -55 +55 @@
         return true;
     }
+    // generate the residual stamp grid for visualization
+    if (!ppSubResidualSampleJpeg(config)) {
+        psError(psErrorCodeLast(), false, "Unable to update.");
+        return false;
+    }
 
     psMetadataAddF32(data->stats, PS_LIST_TAIL, "TIME_MATCH", 0, "Time to match PSFs",
@@ -66 +71 @@
 
     if (!ppSubLowThreshold(data)) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to threshold images.");
+        psError(psErrorCodeLast(), false, "Unable to threshold images.");
         return false;
     }
@@ -77 +82 @@
 
     if (!ppSubDefineOutput("PPSUB.OUTPUT", config)) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to define output.");
+        psError(psErrorCodeLast(), false, "Unable to define output.");
         return false;
     }
 
     if (!data->quality && !ppSubMakePSF(data)) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to generate PSF.");
-        return false;
+        psError(psErrorCodeLast(), false, "Unable to generate PSF.");
+        return false;
+    }
+
+    // Now we've got a PSF, blow away detections in case they're confused with real output detections
+    {
+        pmFPAview *view = ppSubViewReadout(); // View to readout
+        pmReadout *out = pmFPAfileThisReadout(config->files, view, "PPSUB.OUTPUT");
+        if (psMetadataLookup(out->analysis, "PSPHOT.DETECTIONS")) {
+            psMetadataRemoveKey(out->analysis, "PSPHOT.DETECTIONS");
+        }
+        psFree(view);
     }
 
     if (!ppSubReadoutSubtract(config)) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to subtract images.");
+        psError(psErrorCodeLast(), false, "Unable to subtract images.");
         return false;
     }
@@ -99 +114 @@
     // Higher order background subtraction using psphot
     if (!ppSubBackground(config)) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to subtract background.");
+        psError(psErrorCodeLast(), false, "Unable to subtract background.");
         return false;
     }
@@ -105 +120 @@
     // Perform Variance correction (rescale within a modest range)
     if (!ppSubVarianceRescale(config)) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to rescale variance.");
+        psError(psErrorCodeLast(), false, "Unable to rescale variance.");
         return false;
     }
@@ -115 +130 @@
 
     if (!data->quality && !ppSubReadoutPhotometry("PPSUB.OUTPUT", data)) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to perform photometry.");
+        psError(psErrorCodeLast(), false, "Unable to perform photometry.");
         return false;
     }
@@ -126 +141 @@
     // Perform statistics on the cell
     if (!ppSubReadoutStats(data)) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to collect statistics");
-        return false;
-    }
-
+        psError(psErrorCodeLast(), false, "Unable to collect statistics");
+        return false;
+    }
+
+    // generate the binned image used to write the jpeg
     if (!ppSubReadoutJpeg(config)) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to update.");
+        psError(psErrorCodeLast(), false, "Unable to update.");
         return false;
     }
@@ -143 +159 @@
 
         if (data->inverse && !ppSubDefineOutput("PPSUB.INVERSE", config)) {
-            psError(PS_ERR_UNKNOWN, false, "Unable to define inverse.");
+            psError(psErrorCodeLast(), false, "Unable to define inverse.");
             return false;
         }
 
         if (!ppSubReadoutInverse(config)) {
-            psError(PS_ERR_UNKNOWN, false, "Unable to invert images.");
+            psError(psErrorCodeLast(), false, "Unable to invert images.");
             return false;
         }
@@ -164 +180 @@
 
         if (!data->quality && !ppSubReadoutPhotometry("PPSUB.INVERSE", data)) {
-            psError(PS_ERR_UNKNOWN, false, "Unable to perform photometry.");
+            psError(psErrorCodeLast(), false, "Unable to perform photometry.");
             return false;
         }

branches/simtest_nebulous_branches/ppSub/src/ppSubMakePSF.c

@@ -21 +21 @@
 
 #include "ppSub.h"
+
+psArray *ppSubSelectPSFSources(psArray *sources);
 
 bool ppSubMakePSF(ppSubData *data)
@@ -57 +59 @@
     pmFPAfile *photFile = psMetadataLookupPtr(&mdok, config->files, "PSPHOT.INPUT"); // Photometry file
     if (!pmFPACopy(photFile->fpa, minuendFile->fpa)) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to copy FPA for photometry");
+        psError(PPSUB_ERR_CONFIG, false, "Unable to copy FPA for photometry");
         psFree(view);
         return false;
     }
+
     pmReadout *photRO = pmFPAviewThisReadout(view, photFile->fpa); // Readout to photometer
-    if (psMetadataLookup(photRO->analysis, "PSPHOT.SOURCES")) {
-        psMetadataRemoveKey(photRO->analysis, "PSPHOT.SOURCES");
+
+    // we need to remove any existing PSPHOT.DETECTIONS (why not do this in psphot?)
+    if (psMetadataLookup(photRO->analysis, "PSPHOT.DETECTIONS")) {
+        psMetadataRemoveKey(photRO->analysis, "PSPHOT.DETECTIONS");
     }
+    if (psMetadataLookup(photRO->parent->parent->analysis, "PSPHOT.PSF")) {
+        psMetadataRemoveKey(photRO->parent->parent->analysis, "PSPHOT.PSF");
+    }
+
+# ifdef TESTING
+    // XXX for testing, dump these images:
+    psphotSaveImage (NULL, photRO->image, "findpsf.im.fits");
+    psphotSaveImage (NULL, photRO->variance, "findpsf.wt.fits");
+    psphotSaveImage (NULL, photRO->mask, "findpsf.mk.fits");
+# endif
 
     // Extract the loaded sources from the associated readout, and generate PSF
     // Here, we assume the image is background-subtracted
-    psArray *sources = psMetadataLookupPtr(&mdok, minuend->analysis, "PSPHOT.SOURCES");
-    if (!psphotReadoutFindPSF(config, view, sources)) {
+    pmDetections *detections = psMetadataLookupPtr(&mdok, minuend->analysis, "PSPHOT.DETECTIONS");
+    if (!detections || !detections->allSources) {
+        psError(PPSUB_ERR_CONFIG, true, "No sources from which to determine PSF.");
+        return false;
+    }
+    psArray *sources = detections->allSources;
+
+    // XXX filter sources?  limit the total number and return only brighter objects?
+    // use flags to toss totally bogus entries?
+    psArray *goodSources = ppSubSelectPSFSources (sources);
+    if (!psphotReadoutFindPSF(config, view, goodSources)) {
         // This is likely a data quality issue
         // XXX Split into multiple cases using error codes?
@@ -76 +100 @@
         ppSubDataQuality(data, PSPHOT_ERR_PSF, PPSUB_FILES_PHOT_SUB | PPSUB_FILES_PHOT_INV);
         psFree(view);
+        psFree(goodSources);
         return true;
+    }
+
+    // save the resulting PSF information on the pmFPAfile PSPHOT.PSF.LOAD
+    pmFPAfile *psfFile = psMetadataLookupPtr(&mdok, config->files, "PSPHOT.PSF.LOAD"); // PSF file
+    if (!ppSubCopyPSF(psfFile, photFile, view)) {
+        psErrorStackPrint(stderr, "PSF was not generated");
+        psWarning("PSF was not generated --- suspect bad data quality.");
+        ppSubDataQuality(data, psErrorCodeLast(), PPSUB_FILES_PHOT_SUB | PPSUB_FILES_PHOT_INV);
     }
 
@@ -82 +115 @@
     psMetadataRemoveKey(photRO->analysis, "PSPHOT.HEADER");
 
-    if (!data->quality) {
-        data->psf = psMemIncrRefCounter(psMetadataLookupPtr(NULL, photRO->parent->parent->analysis,
-                                                            "PSPHOT.PSF"));
-    }
-
+    psFree(goodSources);
     psFree(view);
 
     return true;
 }
+
+bool ppSubCopyPSF(pmFPAfile *output, pmFPAfile *input, pmFPAview *view)
+{
+    pmChip *inputChip   = pmFPAviewThisChip(view, input->fpa); // Chip with PSF info
+    pmChip *outputChip  = pmFPAviewThisChip(view, output->fpa); // Chip to store PSF info
+
+    pmReadout *inputRO  = pmFPAviewThisReadout(view, input->fpa); // Readout with PSF info
+    pmReadout *outputRO = pmFPAviewThisReadout(view, output->fpa); // Readout to store PSF info
+
+    if (!outputRO) {
+        pmCell *outputCell  = pmFPAviewThisCell(view, output->fpa);
+        outputRO = pmReadoutAlloc(outputCell);
+        outputRO->image = psMemIncrRefCounter(inputRO->image);
+    }
+
+    // Copy the PSF-related data
+    psMetadataIterator *iter = psMetadataIteratorAlloc(inputRO->analysis, PS_LIST_HEAD, "^PSF\\.CLUMP.*");
+    psMetadataItem *item;               // Item from iteration
+    while ((item = psMetadataGetAndIncrement(iter))) {
+        psMetadataAddItem(outputRO->analysis, item, PS_LIST_TAIL, PS_META_REPLACE);
+    }
+    psFree(iter);
+
+    // copy the PSF model data
+    pmPSF *psf = psMetadataLookupPtr(NULL, inputChip->analysis, "PSPHOT.PSF"); // PSF for photometry
+    if (!psf) {
+        psErrorStackPrint(stderr, "No PSF available");
+        return false;
+    }
+
+    psMetadataAddPtr(outputChip->analysis, PS_LIST_TAIL, "PSPHOT.PSF", PS_DATA_UNKNOWN | PS_META_REPLACE,
+                     "PSF from ppSubMakePSF", psf);
+
+    return true;
+}
+
+
+// XXX hardwired MAX for now
+# define MAX_NPSF 500
+
+psArray *ppSubSelectPSFSources(psArray *sources){
+
+    sources = psArraySort (sources, pmSourceSortBySN);
+
+    psArray *subset = psArrayAllocEmpty(MAX_NPSF);
+
+    int nPSF = 0;
+    for (int i = 0; (nPSF < MAX_NPSF) && (i < sources->n); i++) {
+
+        pmSource *source = sources->data[i];
+        if (!source) continue;
+
+        // skip non-astronomical objects (very likely defects)
+        if (source->type == PM_SOURCE_TYPE_DEFECT) continue;
+        if (source->type == PM_SOURCE_TYPE_SATURATED) continue;
+        if (source->mode & PM_SOURCE_MODE_SATSTAR) continue;
+
+        psArrayAdd (subset, 100, source);
+        nPSF++;
+    }
+
+    return subset;
+}

branches/simtest_nebulous_branches/ppSub/src/ppSubMatchPSFs.c

@@ -18 +18 @@
 #include <pslib.h>
 #include <psmodules.h>
+#include <psphot.h>
 
 #include "ppSub.h"
+
+#define COVAR_FRAC 0.01                 // Truncation fraction for covariance matrix
 
 // Normalise a region on an image
@@ -37 +40 @@
 }
 
+// Measure the PSF for an image
+static float subImagePSF(ppSubData *data, // Processing data
+                         const pmReadout *ro, // Readout for which to measure PSF
+                         psArray *sources     // Sources with positions at which to measure PSF
+    )
+{
+    psAssert(data, "Require processing data");
+    pmConfig *config = data->config;    // Configuration
+    psAssert(config, "Require configuration");
+
+    psAssert(ro, "Need readout");
+    psAssert(sources, "Need sources.");
+
+    pmFPAfile *photFile = psMetadataLookupPtr(NULL, config->files, "PSPHOT.INPUT"); // Photometry file
+    psAssert(photFile, "Need photometry file.");
+    if (!pmFPACopy(photFile->fpa, ro->parent->parent->parent)) {
+        psError(PPSUB_ERR_CONFIG, false, "Unable to copy FPA for photometry");
+        return false;
+    }
+
+    pmFPAview *view = ppSubViewReadout(); // View to readout
+    pmReadout *photRO = pmFPAviewThisReadout(view, photFile->fpa); // Readout to photometer
+
+    if (psMetadataLookup(photRO->analysis, "PSPHOT.DETECTIONS")) {
+        psMetadataRemoveKey(photRO->analysis, "PSPHOT.DETECTIONS");
+    }
+    if (psMetadataLookup(photRO->parent->parent->analysis, "PSPHOT.PSF")) {
+        psMetadataRemoveKey(photRO->parent->parent->analysis, "PSPHOT.PSF");
+    }
+
+    // Extract the loaded sources from the associated readout, and generate PSF
+    // Here, we assume the image is background-subtracted
+    if (!psphotReadoutFindPSF(config, view, sources)) {
+        psErrorStackPrint(stderr, "Unable to determine PSF");
+        psWarning("Unable to determine PSF.");
+        psFree(view);
+        return true;
+    }
+
+    psFree(view);
+
+    psMetadata *header = psMetadataLookupMetadata(NULL, photRO->analysis, "PSPHOT.HEADER");
+    psAssert(header, "Require header.");
+    float fwhm = psMetadataLookupF32(NULL, header, "FWHM_MAJ");
+
+    if (psMetadataLookup(photRO->analysis, "PSPHOT.DETECTIONS")) {
+        psMetadataRemoveKey(photRO->analysis, "PSPHOT.DETECTIONS");
+    }
+    if (psMetadataLookup(photRO->parent->parent->analysis, "PSPHOT.PSF")) {
+        psMetadataRemoveKey(photRO->parent->parent->analysis, "PSPHOT.PSF");
+    }
+
+    return fwhm;
+}
+
+// Scale the kernel parameters according to the PSFs
+static bool subScaleKernel(ppSubData *data, // Processing data
+                           psVector *kernelWidths, // Widths for kernel
+                           int *kernelSize,        // Size of kernel
+                           int *stampSize          // Size of stamps (footprint)
+    )
+{
+    psAssert(data, "Require processing data");
+    pmConfig *config = data->config;    // Configuration
+    psAssert(config, "Require configuration");
+
+    // Nothing to do if pre-calculated kernel exists
+    pmFPAview *view = ppSubViewReadout(); // View to readout
+    pmReadout *kernelRO = pmFPAfileThisReadout(config->files, view, "PPSUB.OUTPUT.KERNELS"); // RO with kernel
+    if (kernelRO) {
+        psFree(view);
+        return true;
+    }
+
+    // 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, "SCALE")) {
+        // No scaling requested
+        psFree(view);
+        return true;
+    }
+
+    // Input images
+    pmReadout *inRO = pmFPAfileThisReadout(config->files, view, "PPSUB.INPUT"); // Input readout
+    pmReadout *refRO = pmFPAfileThisReadout(config->files, view, "PPSUB.REF"); // Reference readout
+
+    // Input sources
+    pmReadout *inSourceRO = pmFPAfileThisReadout(config->files, view, "PPSUB.INPUT.SOURCES");
+    pmReadout *refSourceRO = pmFPAfileThisReadout(config->files, view, "PPSUB.REF.SOURCES");
+    if (!inSourceRO || !refSourceRO) {
+        psWarning("Unable to scale kernel, since no sources were provided.");
+        return true;
+    }
+
+    pmDetections *inDetections = psMetadataLookupPtr(NULL, inSourceRO->analysis, "PSPHOT.DETECTIONS"); // Input sources
+    pmDetections *refDetections = psMetadataLookupPtr(NULL, refSourceRO->analysis, "PSPHOT.DETECTIONS"); // Ref sources
+
+    psFree(view);
+
+    if (!inDetections || !refDetections) {
+        psWarning("Unable to scale kernel, since no sources were provided.");
+        return true;
+    }
+
+    psArray *inSources = inDetections->allSources;
+    psAssert (inSources, "missing inSources?");
+
+    psArray *refSources = refDetections->allSources;
+    psAssert (refSources, "missing refSources?");
+
+    float inFWHM = psMetadataLookupF32(NULL, inRO->parent->parent->concepts, "CHIP.SEEING"); // FWHM for input
+    if (!isfinite(inFWHM) || inFWHM == 0.0) {
+        inFWHM = subImagePSF(data, inRO, inSources);
+    }
+    float refFWHM = psMetadataLookupF32(NULL, refRO->parent->parent->concepts, "CHIP.SEEING"); // FWHM for ref
+    if (!isfinite(refFWHM) || refFWHM == 0.0) {
+        refFWHM = subImagePSF(data, refRO, refSources);
+    }
+    psLogMsg("ppSub", PS_LOG_INFO, "Input FWHM: %f\nReference FWHM: %f\n", inFWHM, refFWHM);
+    if (!isfinite(inFWHM) || !isfinite(refFWHM)) {
+        psWarning("Unable to scale kernel, since unable to measure PSFs.");
+        return true;
+    }
+
+    float scaleRef = psMetadataLookupF32(NULL, recipe, "SCALE.REF"); // Reference for scaling
+    float scaleMin = psMetadataLookupF32(NULL, recipe, "SCALE.MIN"); // Minimum for scaling
+    float scaleMax = psMetadataLookupF32(NULL, recipe, "SCALE.MAX"); // Maximum for scaling
+    if (!isfinite(scaleRef) || !isfinite(scaleMin) || !isfinite(scaleMax)) {
+        psError(PPSUB_ERR_ARGUMENTS, false,
+                "Scale parameters (SCALE.REF=%f, SCALE.MIN=%f, SCALE.MAX=%f) not set in recipe.",
+                scaleRef, scaleMin, scaleMax);
+        return false;
+    }
+
+    if (!pmSubtractionParamsScale(kernelSize, stampSize, kernelWidths, inFWHM, refFWHM,
+                                  scaleRef, scaleMin, scaleMax)) {
+        psError(PPSUB_ERR_DATA, false, "Unable to scale parameters.");
+        return false;
+    }
+
+    return true;
+}
+
 
 bool ppSubMatchPSFs(ppSubData *data)
@@ -71 +218 @@
     pmReadout *kernelRO = pmFPAfileThisReadout(config->files, view, "PPSUB.OUTPUT.KERNELS"); // RO with kernel
 
-    psFree(view);
-
     // Sources in image, used for stamps: these must be loaded from previous analysis stages
     pmReadout *inSourceRO = pmFPAfileThisReadout(config->files, view, "PPSUB.INPUT.SOURCES");
     pmReadout *refSourceRO = pmFPAfileThisReadout(config->files, view, "PPSUB.REF.SOURCES");
-    psArray *inSources = inSourceRO ? psMetadataLookupPtr(&mdok, inSourceRO->analysis, "PSPHOT.SOURCES") :
-        NULL; // Source list from input image
-    psArray *refSources = refSourceRO ? psMetadataLookupPtr(&mdok, refSourceRO->analysis, "PSPHOT.SOURCES") :
-        NULL ; // Source list from reference image
-
-    psArray *sources = NULL;            // Merged list of sources
-    if (inSources && refSources) {
+
+    pmDetections *inDetections  = inSourceRO  ? psMetadataLookupPtr(&mdok, inSourceRO->analysis,  "PSPHOT.DETECTIONS") : NULL; // Input sources
+    pmDetections *refDetections = refSourceRO ? psMetadataLookupPtr(&mdok, refSourceRO->analysis, "PSPHOT.DETECTIONS") : NULL; // Ref sources
+
+    psFree(view);
+
+    pmDetections *detections = NULL;    // Merged detection set
+    if (inDetections && refDetections) {
+        psArray *inSources  = inDetections->allSources;
+        psArray *refSources = refDetections->allSources;
+
+        psAssert (inSources, "missing in sources?");
+        psAssert (refSources, "missing ref sources?");
+
+        detections = pmDetectionsAlloc();
         float radius = psMetadataLookupF32(NULL, recipe, "SOURCE.RADIUS"); // Matching radius
         psArray *lists = psArrayAlloc(2); // Source lists
         lists->data[0] = psMemIncrRefCounter(inSources);
         lists->data[1] = psMemIncrRefCounter(refSources);
-        sources = pmSourceMatchMerge(lists, radius);
+        detections->allSources = pmSourceMatchMerge(lists, radius);
         psFree(lists);
-        if (!sources) {
-            psError(PS_ERR_UNKNOWN, false, "Unable to merge source lists");
+        if (!detections->allSources) {
+            psError(PPSUB_ERR_DATA, false, "Unable to merge source lists");
+            psFree(detections);
             return false;
         }
-    } else if (inSources) {
-        sources = psMemIncrRefCounter(inSources);
-    } else if (refSources) {
-        sources = psMemIncrRefCounter(refSources);
-    }
-
-    psMetadataAddArray(inConv->analysis, PS_LIST_TAIL, "PSPHOT.SOURCES", PS_META_REPLACE,
-                       "Merged source list", sources);
-    psMetadataAddArray(refConv->analysis, PS_LIST_TAIL, "PSPHOT.SOURCES", PS_META_REPLACE,
-                       "Merged source list", sources);
-    psFree(sources);                    // Drop reference
+    }
+    if (!detections && inDetections) {
+        detections = psMemIncrRefCounter(inDetections);
+    }
+    if (!detections && refDetections) {
+        detections = psMemIncrRefCounter(refDetections);
+    }
+
+    psMetadataAddPtr(inConv->analysis,  PS_LIST_TAIL, "PSPHOT.DETECTIONS", PS_META_REPLACE | PS_DATA_UNKNOWN, "Merged source list", detections);
+    psMetadataAddPtr(refConv->analysis, PS_LIST_TAIL, "PSPHOT.DETECTIONS", PS_META_REPLACE | PS_DATA_UNKNOWN, "Merged source list", detections);
+    psFree(detections);                    // Drop reference
 
     int footprint = psMetadataLookupS32(NULL, recipe, "STAMP.FOOTPRINT"); // Stamp half-size
@@ -115 +269 @@
     pmSubtractionKernelsType type = pmSubtractionKernelsTypeFromString(typeStr); // Type of kernel
     if (type == PM_SUBTRACTION_KERNEL_NONE) {
-        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Unrecognised kernel type: %s", typeStr);
+        psError(PPSUB_ERR_ARGUMENTS, true, "Unrecognised kernel type: %s", typeStr);
         return false;
     }
@@ -130 +284 @@
     int iter = psMetadataLookupS32(NULL, recipe, "ITER"); // Rejection iterations
     float rej = psMetadataLookupF32(NULL, recipe, "REJ"); // Rejection threshold
-    float sys = psMetadataLookupF32(NULL, recipe, "SYS"); // Relative systematic error in kernel
+    float kernelErr = psMetadataLookupF32(NULL, recipe, "KERNEL.ERR"); // Relative systematic error in kernel
+    float normFrac = psMetadataLookupF32(NULL, recipe, "NORM.FRAC"); // Fraction of window for normalisn windw
+    float sysErr = psMetadataLookupF32(NULL, recipe, "SYS.ERR"); // Relative systematic error in images
+    float skyErr = psMetadataLookupF32(NULL, recipe, "SKY.ERR"); // Additional error in sky
+    float covarFrac = psMetadataLookupF32(NULL, recipe, "COVAR.FRAC"); // Fraction for covariance calculation
 
     float badFrac = psMetadataLookupF32(NULL, recipe, "BADFRAC"); // Maximum bad fraction
@@ -168 +326 @@
             break;
           default:
-            psErrorStackPrint(stderr, "Invalid value for -convolve");
+            psError(PPSUB_ERR_ARGUMENTS, false, "Invalid value for -convolve");
             return false;
         }
     }
 
+    if (!subScaleKernel(data, widths, &size, &footprint)) {
+        psError(PPSUB_ERR_DATA, false, "Unable to scale kernel parameters");
+        return false;
+    }
+
     int threads = psMetadataLookupS32(NULL, config->arguments, "-threads"); // Number of threads
     if (threads > 0) {
-        pmSubtractionThreadsInit(inRO, refRO);
-    }
+        pmSubtractionThreadsInit();
+    }
+
+    if (inRO->covariance) {
+        psImageCovarianceTruncate(inRO->covariance, COVAR_FRAC);
+    }
+    if (refRO->covariance) {
+        psImageCovarianceTruncate(refRO->covariance, COVAR_FRAC);
+    }
+
+    // Data doesn't exist until it's created in pmSubtraction...
+    inConv->data_exists = inConv->parent->data_exists = inConv->parent->parent->data_exists = false;
+    refConv->data_exists = refConv->parent->data_exists = refConv->parent->parent->data_exists = false;
 
     // Match the PSFs
@@ -182 +356 @@
     if (kernelRO) {
         success = pmSubtractionMatchPrecalc(inConv, refConv, inRO, refRO, kernelRO->analysis,
-                                            stride, sys, maskVal, maskBad, maskPoor, poorFrac, badFrac);
+                                            stride, kernelErr, covarFrac, maskVal, maskBad, maskPoor,
+                                            poorFrac, badFrac);
     } else {
         success = pmSubtractionMatch(inConv, refConv, inRO, refRO, footprint, stride, regionSize,
-                                     spacing, threshold, sources, data->stamps, type, size, order,
-                                     widths, orders, inner, ringsOrder, binning, penalty, optimum,
-                                     optWidths, optOrder, optThresh, iter, rej, sys, maskVal,
-                                     maskBad, maskPoor, poorFrac, badFrac, subMode);
-    }
+                                     spacing, threshold, detections ? detections->allSources : NULL,
+                                     data->stamps, type, size, order, widths, orders, inner, ringsOrder,
+                                     binning, penalty, optimum, optWidths, optOrder, optThresh, iter, rej,
+                                     normFrac, sysErr, skyErr, kernelErr, covarFrac,
+                                     maskVal, maskBad, maskPoor, poorFrac, badFrac, subMode);
+    }
+
+# ifdef TESTING
+    // XXX for testing
+    psphotSaveImage (NULL, refRO->image,    "refRO.im.fits");
+    psphotSaveImage (NULL, refRO->variance, "refRO.wt.fits");
+    psphotSaveImage (NULL, refRO->mask,     "refRO.mk.fits");
+
+    psphotSaveImage (NULL, inRO->image,    "inRO.im.fits");
+    psphotSaveImage (NULL, inRO->variance, "inRO.wt.fits");
+    psphotSaveImage (NULL, inRO->mask,     "inRO.mk.fits");
+
+    psphotSaveImage (NULL, inConv->image,    "inConv.im.fits");
+    psphotSaveImage (NULL, inConv->variance, "inConv.wt.fits");
+    psphotSaveImage (NULL, inConv->mask,     "inConv.mk.fits");
+
+    psphotSaveImage (NULL, refConv->image,    "refConv.im.fits");
+    psphotSaveImage (NULL, refConv->variance, "refConv.wt.fits");
+    psphotSaveImage (NULL, refConv->mask,     "refConv.mk.fits");
+# endif
 
     psFree(optWidths);
-    pmSubtractionThreadsFinalize(inRO, refRO);
+    pmSubtractionThreadsFinalize();
 
     if (!success) {
@@ -207 +402 @@
             return true;
         } else {
-            psError(PS_ERR_UNKNOWN, false, "Unable to match images.");
+            psError(PPSUB_ERR_DATA, false, "Unable to match images.");
             return false;
         }
@@ -260 +455 @@
     pmConceptsCopyFPA(refConv->parent->parent->parent, refRO->parent->parent->parent, true, true);
 
-    psImageCovarianceTransfer(inConv->variance, inConv->covariance);
-    psImageCovarianceTransfer(refConv->variance, refConv->covariance);
+    if (inConv->covariance) {
+        psImageCovarianceTruncate(inConv->covariance, COVAR_FRAC);
+    }
+    if (refConv->covariance) {
+        psImageCovarianceTruncate(refConv->covariance, COVAR_FRAC);
+    }
+
+    if (inConv->variance) {
+        psImageCovarianceTransfer(inConv->variance, inConv->covariance);
+    }
+    if (refConv->variance) {
+        psImageCovarianceTransfer(refConv->variance, refConv->covariance);
+    }
 
     return true;

branches/simtest_nebulous_branches/ppSub/src/ppSubReadoutInverse.c

@@ -20 +20 @@
     invRO->variance = psMemIncrRefCounter(outRO->variance);
     invRO->covariance = psMemIncrRefCounter(outRO->covariance);
+    invRO->analysis = psMetadataCopy(invRO->analysis, outRO->analysis);
 
     invRO->data_exists = invRO->parent->data_exists = invRO->parent->parent->data_exists = true;
@@ -34 +35 @@
     pmHDU *invHDU = invFPA->hdu;          // Inverse HDU
     if (!pmAstromWriteWCS(invHDU->header, outFPA, outChip, WCS_TOLERANCE)) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to write WCS astrometry to PPSUB.INVERSE.");
+        psError(psErrorCodeLast(), false, "Unable to write WCS astrometry to PPSUB.INVERSE.");
         return false;
     }
     // Read from newly written astrometry so that it exists in the "inverse" FPA (for sources)
     if (!pmAstromReadWCS(invFPA, invChip, invHDU->header, 1.0)) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to read WCS astrometry.");
+        psError(psErrorCodeLast(), false, "Unable to read WCS astrometry.");
         return false;
     }

branches/simtest_nebulous_branches/ppSub/src/ppSubReadoutJpeg.c

@@ -33 +33 @@
     pmReadout *ro1 = pmReadoutAlloc(cell1), *ro2 = pmReadoutAlloc(cell2); // Binned readouts
     if (!pmReadoutRebin(ro1, outRO, maskBad, bin1, bin1)) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to bin output (1st binning)");
+        psError(PPSUB_ERR_DATA, false, "Unable to bin output (1st binning)");
         psFree(ro1);
         psFree(ro2);
@@ -39 +39 @@
     }
     if (!pmReadoutRebin(ro2, ro1, 0, bin2, bin2)) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to bin output (2nd binning)");
+        psError(PPSUB_ERR_DATA, false, "Unable to bin output (2nd binning)");
         psFree(ro1);
         psFree(ro2);
@@ -49 +49 @@
     return true;
 }
+
+bool ppSubResidualSampleJpeg(pmConfig *config)
+{
+    return true;
+
+    pmFPAview *view = ppSubViewReadout(); // View to readout
+
+    // we save sample difference stamps on the convolved image readout->analysis metadata
+    pmReadout *inConv = pmFPAfileThisReadout(config->files, view, "PPSUB.INPUT.CONV"); // Input convolved
+
+    psArray *samples = psArrayAllocEmpty(16); // Array of sample stamp images
+
+    // we may have multiple entries with the same name; pull them off into a separate array:
+    psString regex = NULL;          // Regular expression
+    psStringAppend(&regex, "^%s$", "SUBTRACTION.SAMPLE.STAMP.SET");
+    psMetadataIterator *iter = psMetadataIteratorAlloc(inConv->analysis, PS_LIST_HEAD, regex); // Iterator
+    psFree(regex);
+
+    psMetadataItem *item = NULL;// Item from iteration
+    while ((item = psMetadataGetAndIncrement(iter))) {
+        assert(item->type == PS_DATA_ARRAY);
+        psArray *sampleStamps = item->data.V;
+        samples = psArrayAdd(samples, 16, sampleStamps);
+    }
+    psFree(iter);
+    psAssert (samples, "no sample stamps?");
+    psAssert (samples->n, "no sample stamps?");
+
+    // get the kernel sizes
+    psArray *kernels = samples->data[0];
+    psAssert (kernels, "no valid kernel?");
+    psAssert (kernels->n, "no valid kernel?");
+
+    psImage *kernel = kernels->data[0];
+    psAssert (kernel, "missing valid kernel?");
+
+    int DX = kernel->numCols;
+    int DY = kernel->numRows;
+
+    // each array contains up to 9 sample stamps.  generate an image mosaicking these together in 3x3 blocks.
+    int innerBorder = 1;
+    int outerBorder = 2;
+    int NXblock = sqrt(samples->n);
+    int NYblock = samples->n / NXblock;
+    if (samples->n % NXblock) NYblock ++;
+
+    int NXpix = NXblock * (3 * (DX + innerBorder) + outerBorder);
+    int NYpix = NYblock * (3 * (DY + innerBorder) + outerBorder);
+
+    // output cell
+    pmCell *cell = pmFPAfileThisCell(config->files, view, "PPSUB.OUTPUT.RESID.JPEG");
+    pmReadout *readout = pmReadoutAlloc(cell);
+    readout->image = psImageAlloc(NXpix, NYpix, PS_TYPE_F32);
+
+    cell->data_exists = true;
+    cell->parent->data_exists = true;
+
+    psImageInit (readout->image, 0.0);
+
+    for (int i = 0; i < samples->n; i++) {
+
+        int xBlock = i % NXblock;
+        int yBlock = i / NYblock;
+
+        psArray *kernels = samples->data[i];
+
+        for (int j = 0; j < kernels->n; j++) {
+
+            psImage *kernel = kernels->data[j];
+
+            int xSub = j % 3;
+            int ySub = j / 3;
+
+            int xPix = xBlock * (3 * (DX + innerBorder) + outerBorder) + xSub * (DX + innerBorder);
+            int yPix = yBlock * (3 * (DX + innerBorder) + outerBorder) + ySub * (DY + innerBorder);
+
+            for (int y = 0; y < kernel->numRows; y++) {
+                for (int x = 0; x < kernel->numCols; x++) {
+                    readout->image->data.F32[y + yPix][x + xPix] = kernel->data.F32[y][x];
+                }
+            }
+        }
+    }
+
+    {
+        psFits *fits = psFitsOpen ("resid.stamps.fits", "w");
+        psFitsWriteImage (fits, NULL, readout->image, 0, NULL);
+        psFitsClose (fits);
+    }
+
+    return true;
+}

branches/simtest_nebulous_branches/ppSub/src/ppSubReadoutPhotometry.c

@@ -24 +24 @@
 bool ppSubReadoutPhotometry(const char *name, ppSubData *data)
 {
+    bool mdok = false;
+
     psAssert(data, "Require processing data");
     pmConfig *config = data->config;    // Configuration
@@ -42 +44 @@
     }
 
-    // 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
+    // select the view of interest
     pmFPAview *view = ppSubViewReadout(); // View to readout
-    pmChip *psfInputChip = pmFPAfileThisChip(config->files, view, "PSPHOT.INPUT"); // Chip with PSF
-    psAssert (psfInputChip, "should have been generated for ppSubMakePSF");
-    pmChip *psfLoadChip = pmFPAfileThisChip(config->files, view, "PSPHOT.PSF.LOAD"); // Chip to have PSF
-    psAssert (psfLoadChip, "PSPHOT.PSF.LOAD should have been defined in ppSubCamera");
-    pmPSF *psf = psMetadataLookupPtr(NULL, psfInputChip->analysis, "PSPHOT.PSF"); // PSF for photometry
-    if (!psf) {
-        psErrorStackPrint(stderr, "No PSF available");
-        psWarning("No PSF available --- suspect bad data quality.");
-        ppSubDataQuality(data, psErrorCodeLast(), PPSUB_FILES_PHOT_SUB | PPSUB_FILES_PHOT_INV);
-        return true;
-    }
-    psMetadataAddPtr(psfLoadChip->analysis, PS_LIST_TAIL, "PSPHOT.PSF", PS_DATA_UNKNOWN | PS_META_REPLACE,
-                     "PSF from ppSubMakePSF", psf);
-
-    bool mdok = false;
 
     // 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.
-    pmReadout *inRO = pmFPAfileThisReadout(config->files, view, name); // Readout with image and sources
-    if (psMetadataLookup(inRO->analysis, "PSPHOT.SOURCES")) {
-        psMetadataRemoveKey(inRO->analysis, "PSPHOT.SOURCES");
-    }
-
+    // around the pointers so PSPHOT.INPUT corresponds to the output image of interest (on one
+    // pass this is the subtraction image, in another it is negative of the subtraction
     pmFPAfile *photFile = psMetadataLookupPtr(&mdok, config->files, "PSPHOT.INPUT"); // Photometry file
     pmFPAfile *inFile = psMetadataLookupPtr(&mdok, config->files, name); // Input file
     if (!pmFPACopy(photFile->fpa, inFile->fpa)) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to copy FPA for photometry");
+        psError(PPSUB_ERR_CONFIG, false, "Unable to copy FPA for photometry");
         psFree(view);
         return false;
     }
+
+    // drop references to PSPHOT.DETECTIONS on both of these  (why is this needed for both??)
     pmReadout *photRO = pmFPAviewThisReadout(view, photFile->fpa); // Readout to photometer
-    if (psMetadataLookup(photRO->analysis, "PSPHOT.SOURCES")) {
-        psMetadataRemoveKey(photRO->analysis, "PSPHOT.SOURCES");
+    if (psMetadataLookup(photRO->analysis, "PSPHOT.DETECTIONS")) {
+        psMetadataRemoveKey(photRO->analysis, "PSPHOT.DETECTIONS");
+    }
+    pmReadout *inRO = pmFPAfileThisReadout(config->files, view, name); // Readout with image and sources
+    if (psMetadataLookup(inRO->analysis, "PSPHOT.DETECTIONS")) {
+        psMetadataRemoveKey(inRO->analysis, "PSPHOT.DETECTIONS");
     }
 
-    psMetadataAddPtr(photRO->parent->parent->analysis, PS_LIST_TAIL, "PSPHOT.PSF",
-                     PS_META_REPLACE | PS_DATA_UNKNOWN, "Point-spread function", data->psf);
+    // grab the PSF information from the pmFPAfile PSPHOT.PSF.LOAD
+    pmFPAfile *psfFile = psMetadataLookupPtr(&mdok, config->files, "PSPHOT.PSF.LOAD"); // PSF file
+    ppSubCopyPSF (photFile, psfFile, view);
 
-    psFree(photRO->analysis);
-    photRO->analysis = psMetadataAlloc();
+    // psphotSaveImage (photFile->fpa->hdu->header, photRO->image, "findsrc.im.fits");
+    // psphotSaveImage (photFile->fpa->hdu->header, photRO->variance, "findsrc.wt.fits");
+    // psphotSaveImage (photFile->fpa->hdu->header, photRO->mask, "findsrc.mk.fits");
+
+    // erase the overlays from a previous psphot-related step
+    if (pmVisualIsVisual()) {
+        //      psphotVisualEraseOverlays (1, "all");
+    }
 
     if (!psphotReadoutMinimal(config, view)) {
@@ -96 +90 @@
 
     // If no sources were found, there's no error,  but we want to trigger 'bad quality'
-    psArray *sources = psMetadataLookupPtr(NULL, photRO->analysis, "PSPHOT.SOURCES"); // Sources
-    if (!sources) {
+    pmDetections *detections = psMetadataLookupPtr(NULL, photRO->analysis, "PSPHOT.DETECTIONS"); // Sources
+    if (!detections) {
         ppSubDataQuality(data, PSPHOT_ERR_DATA, PPSUB_FILES_PHOT_SUB | PPSUB_FILES_PHOT_INV);
     }
+    psArray *sources = detections->allSources;
+    psAssert (sources, "missing sources?");
 
     if (data->stats) {
@@ -114 +110 @@
 
     if (!data->quality) {
-        if (!psMetadataCopySingle(inRO->analysis, photRO->analysis, "PSPHOT.SOURCES")) {
-            psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to copy PSPHOT.SOURCES");
+        if (!psMetadataCopySingle(inRO->analysis, photRO->analysis, "PSPHOT.DETECTIONS")) {
+            psError(PPSUB_ERR_PROG, false, "Unable to copy PSPHOT.DETECTIONS");
             return false;
         }
         if (!psMetadataCopySingle(inRO->analysis, photRO->analysis, "PSPHOT.HEADER")) {
-            psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to copy PSPHOT.HEADER");
+            psError(PPSUB_ERR_PROG, false, "Unable to copy PSPHOT.HEADER");
             return false;
+        }
+        if (!psMetadataCopySingle(inRO->analysis, photRO->analysis, PM_DETEFF_ANALYSIS)) {
+            psError(PPSUB_ERR_PROG, false, "Unable to copy Detection Efficiency");
+            return false;
+        }
+
+        // Ensure photometry information is put in the header
+        pmHDU *hdu = pmHDUFromReadout(inRO); // HDU for readout
+        if (hdu) {
+            psMetadata *photHeader = psMetadataLookupMetadata(NULL, inRO->analysis, "PSPHOT.HEADER"); // Header
+            hdu->header = psMetadataCopy(hdu->header, photHeader);
         }
     }
@@ -130 +137 @@
 }
 
-
-
-
-
 #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
+        pmDetections *detections = psMetadataLookupPtr(NULL, photRO->analysis, "PSPHOT.DETECTIONS"); // Sources
+        psArray *sources = detections->allSources;
         FILE *sourceFile = fopen("sources.dat", "w"); // File for sources
         fprintf(sourceFile,

branches/simtest_nebulous_branches/ppSub/src/ppSubReadoutStats.c

@@ -23 +23 @@
     pmFPAfile *output = psMetadataLookupPtr(NULL, config->files, "PPSUB.OUTPUT"); // Output file
     if (!output) {
-        psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find file PPSUB.OUTPUT.\n");
+        psError(PPSUB_ERR_PROG, true, "Unable to find file PPSUB.OUTPUT.\n");
         return false;
     }

branches/simtest_nebulous_branches/ppSub/src/ppSubReadoutSubtract.c

@@ -71 +71 @@
     pmFPA *outFPA = outFile->fpa;       // Output FPA
     if (!pmConceptsCopyFPA(outFPA, inFPA, true, true)) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to copy concepts from input to output.");
+        psError(PPSUB_ERR_CONFIG, false, "Unable to copy concepts from input to output.");
         psFree(outRO);
         psFree(view);
@@ -84 +84 @@
     psFree(view);
     if (!outHDU || !inHDU) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to find HDU at FPA level to copy astrometry.");
+        psError(PPSUB_ERR_PROG, true, "Unable to find HDU at FPA level to copy astrometry.");
         return false;
     }
     if (!pmAstromReadWCS(outFPA, outChip, inHDU->header, 1.0)) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to read WCS astrometry from input FPA.");
+        psError(psErrorCodeLast(), false, "Unable to read WCS astrometry from input FPA.");
         return false;
     }
     if (!pmAstromWriteWCS(outHDU->header, outFPA, outChip, WCS_TOLERANCE)) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to write WCS astrometry to output FPA.");
+        psError(psErrorCodeLast(), false, "Unable to write WCS astrometry to output FPA.");
         return false;
     }

branches/simtest_nebulous_branches/ppSub/src/ppSubSetMasks.c

@@ -28 +28 @@
     psImageMaskType maskValue, markValue; // Mask values
     if (!pmConfigMaskSetBits(&maskValue, &markValue, config)) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to determine mask value.");
+        psError(PPSUB_ERR_CONFIG, false, "Unable to determine mask value.");
         return false;
     }
 
     // Set the mask bits needed by psphot (in psphot recipe)
-    psphotSetMaskRecipe (config, maskValue, markValue);
+    psphotSetMaskRecipe(config, maskValue, markValue);
 
     // Look up recipe values
@@ -79 +79 @@
     // Mask the NAN values
     if (!pmReadoutMaskNonfinite(inRO, satValue)) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to mask non-finite pixels in input.");
+        psError(PPSUB_ERR_DATA, false, "Unable to mask non-finite pixels in input.");
         return false;
     }
     if (!pmReadoutMaskNonfinite(refRO, satValue)) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to mask non-finite pixels in reference.");
+        psError(PPSUB_ERR_DATA, false, "Unable to mask non-finite pixels in reference.");
         return false;
     }
@@ -94 +94 @@
     psImageInterpolateMode interpMode = psImageInterpolateModeFromString(interpModeStr); // Interp
     if (interpMode == PS_INTERPOLATE_NONE) {
-        psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Unknown interpolation mode: %s", interpModeStr);
+        psError(PPSUB_ERR_CONFIG, false, "Unknown interpolation mode: %s", interpModeStr);
         return false;
     }
@@ -104 +104 @@
     // Interpolate over bad pixels, so the bad pixels don't explode
     if (!pmReadoutInterpolateBadPixels(inRO, maskVal, interpMode, poorFrac, maskPoor, maskBad)) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to interpolate bad pixels for input image.");
+        psError(PPSUB_ERR_DATA, false, "Unable to interpolate bad pixels for input image.");
         return false;
     }
     if (!pmReadoutInterpolateBadPixels(refRO, maskVal, interpMode, poorFrac, maskPoor, maskBad)) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to interpolate bad pixels for reference image.");
+        psError(PPSUB_ERR_DATA, false, "Unable to interpolate bad pixels for reference image.");
         return false;
     }

branches/simtest_nebulous_branches/ppSub/src/ppSubThreshold.c

@@ -27 +27 @@
     psImageMaskType maskIgnore,         // Ignore pixels with this mask
     psImageMaskType maskThresh,         // Give pixels this mask if below threshold
+    psRegion *region,                   // Region of interest
     const char *description             // Description of image
     )
@@ -37 +38 @@
     psImage *image = ro->image;         // Image
     psImage *mask = ro->mask;           // Mask
-    int numCols = ro->image->numCols, numRows = ro->image->numRows; // Size of image
+
+    psImage *subImage = psImageSubset(image, *region); // Image with region of interest
+    psImage *subMask = psImageSubset(mask, *region);  // Maks with region of interest
 
     psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); // Statistics
     psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS);                               // Random number generator
-    if (!psImageBackground(stats, NULL, image, mask, maskIgnore, rng)) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to determine threshold.");
+    if (!psImageBackground(stats, NULL, subImage, subMask, maskIgnore, rng)) {
+        psError(PPSUB_ERR_DATA, false, "Unable to determine threshold.");
         psFree(rng);
         psFree(stats);
@@ -49 +52 @@
     psFree(rng);
 
+    psFree(subImage);
+    psFree(subMask);
+
     float threshold = stats->robustMedian - thresh * stats->robustStdev; // Threshold below which to clip
     psFree(stats);
     psLogMsg("ppSub", PS_LOG_INFO, "Masking pixels below %f in %s", threshold, description);
 
-    for (int y = 0; y < numRows; y++) {
-        for (int x = 0; x < numCols; x++) {
+    for (int y = region->y0; y < region->y1; y++) {
+        for (int x = region->x0; x < region->x1; x++) {
             if (mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskIgnore) {
                 continue;
@@ -91 +97 @@
     pmReadout *in = pmFPAfileThisReadout(config->files, view, "PPSUB.INPUT.CONV"); // Input image
     if (!in) {
-        psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to find readout.");
-        return false;
-    }
-    if (!lowThreshold(in, thresh, maskVal, maskThresh, "input convolved image")) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to threshold input image.");
+        psError(PPSUB_ERR_UNKNOWN, false, "Unable to find readout.");
         return false;
     }
@@ -101 +103 @@
     pmReadout *ref = pmFPAfileThisReadout(config->files, view, "PPSUB.REF.CONV"); // Reference image
     if (!ref) {
-        psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to find readout.");
+        psError(PPSUB_ERR_UNKNOWN, false, "Unable to find readout.");
         return false;
     }
-    if (!lowThreshold(ref, thresh, maskVal, maskThresh, "reference convolved image")) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to threshold input image.");
-        return false;
+
+    psMetadataIterator *regIter = psMetadataIteratorAlloc(in->analysis, PS_LIST_HEAD,
+                                                          "^" PM_SUBTRACTION_ANALYSIS_REGION "$");
+    psMetadataItem *regItem;        // Item with region
+    while ((regItem = psMetadataGetAndIncrement(regIter))) {
+        psAssert(regItem->type == PS_DATA_REGION && regItem->data.V, "Expect region type");
+        psRegion *region = regItem->data.V; // Region of interest
+        if (!lowThreshold(in, thresh, maskVal, maskThresh, region, "input convolved image")) {
+            psError(psErrorCodeLast(), false, "Unable to threshold input image.");
+            return false;
+        }
+        if (!lowThreshold(ref, thresh, maskVal, maskThresh, region, "reference convolved image")) {
+            psError(psErrorCodeLast(), false, "Unable to threshold input image.");
+            return false;
+        }
     }
+    psFree(regIter);
 
     psFree(view);

branches/simtest_nebulous_branches/ppSub/src/ppSubVarianceRescale.c

@@ -28 +28 @@
     bool mdok; // Status of metadata lookups
 
-    psMetadata *ppSubRecipe = psMetadataLookupPtr(NULL, config->recipes, PPSUB_RECIPE); // Recipe for ppSub
-    psAssert(ppSubRecipe, "Need PPSUB recipe");
+    psMetadata *recipe = psMetadataLookupPtr(NULL, config->recipes, PPSUB_RECIPE); // Recipe for ppSub
+    psAssert(recipe, "Need PPSUB recipe");
 
-    if (!psMetadataLookupBool(&mdok, ppSubRecipe, "VARIANCE.RESCALE")) return true; 
+    if (!psMetadataLookupBool(&mdok, recipe, "RENORM")) return true;
+
+    int num = psMetadataLookupS32(&mdok, recipe, "RENORM.NUM");
+    if (!mdok) {
+        psError(PPSUB_ERR_ARGUMENTS, true, "RENORM.NUM is not set in the recipe");
+        return false;
+    }
+    float minValid = psMetadataLookupF32(&mdok, recipe, "RENORM.MIN");
+    if (!mdok) {
+        psError(PPSUB_ERR_ARGUMENTS, true, "RENORM.MIN is not set in the recipe");
+        return false;
+    }
+    float maxValid = psMetadataLookupF32(&mdok, recipe, "RENORM.MAX");
+    if (!mdok) {
+        psError(PPSUB_ERR_ARGUMENTS, true, "RENORM.MAX is not set in the recipe");
+        return false;
+    }
 
     psImageMaskType maskBad = pmConfigMaskGet("BLANK", config); // Bits to mask
 
     pmFPAview *view = ppSubViewReadout(); // View to readout
-    pmReadout *outRO = pmFPAfileThisReadout(config->files, view, "PPSUB.OUTPUT"); // Output image
+    pmReadout *readout = pmFPAfileThisReadout(config->files, view, "PPSUB.OUTPUT"); // Output image
 
-    psImage *image    = outRO->image;    // Image of interest
-    psImage *variance = outRO->variance; // Variance image
-    psImage *mask     = outRO->mask;     // Mask of interest
-
-    psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS);
-
-    psLogMsg("psLib.psImageBackground", PS_LOG_DETAIL, "Generating the signal-to-noise image");
-
-    // generate an image representing the signal/noise:
-    psImage *SN = psImageCopy (NULL, outRO->image, PS_TYPE_F32); 
-
-    int nx = image->numCols; // Size of image
-    int ny = image->numRows; // Size of image
-
-    for (int y = 0; y < ny; y++) {
-        for (int x = 0; x < nx; x++) {
-            if (mask && mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskBad) {
-                SN->data.F32[y][x] = 0.0;
-                continue;
-            } 
-            if (!isfinite(image->data.F32[y][x])) {
-                SN->data.F32[y][x] = 0.0;
-                continue;
-            } 
-            if (!isfinite(variance->data.F32[y][x]) || (variance->data.F32[y][x] < 0.0)) {
-                SN->data.F32[y][x] = 0.0;
-                continue;
-            } 
-            SN->data.F32[y][x] = image->data.F32[y][x] / sqrt(variance->data.F32[y][x]);
-        }
+    if (!readout->variance) {
+        // Nothing to renormalise
+        psWarning("Renormalisation of the variance requested, but no variance provided.");
+        return true;
     }
 
-    // these should not be chosen by the user: only a ROBUST version makes sense
-    psStats *stats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
-
-    int Nsubset = psMetadataLookupS32(&mdok, ppSubRecipe, "VARIANCE.RESCALE.NSAMPLE");
-    psAssert (mdok, "missing VARIANCE.RESCALE.NSAMPLE in ppSub recipe");
-    
-    const int Npixels = nx*ny; // Total number of pixels
-    Nsubset = PS_MIN(Nsubset, Npixels); // Number of pixels in subset
-    psLogMsg("psLib.psImageBackground", PS_LOG_DETAIL, "Searching for %d pixels in S/N image", Nsubset);
-
-    psVector *values = psVectorAllocEmpty(Nsubset, PS_TYPE_F32); // Vector containing subsample
-
-    // Minimum and maximum values
-    float min = +PS_MAX_F32;
-    float max = -PS_MAX_F32;
-
-    // select a subset of the image pixels to measure the stats
-    long n = 0;                         // Number of actual pixels in subset
-    if (Nsubset >= Npixels) {
-        // if we have an image smaller than Nsubset, just loop over the image pixels
-        for (int iy = 0; iy < ny; iy++) {
-            for (int ix = 0; ix < nx; ix++) {
-                if (!isfinite(image->data.F32[iy][ix]) || (mask && mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & maskBad)) {
-                    continue;
-                }
-
-                float value = SN->data.F32[iy][ix];
-                min = PS_MIN(value, min);
-                max = PS_MAX(value, max);
-                values->data.F32[n] = value;
-                n++;
-            }
-        }
-    } else {
-        for (long i = 0; i < Nsubset; i++) {
-            double frnd = psRandomUniform(rng);
-            int pixel = Npixels * frnd;
-            int ix = pixel % nx;
-            int iy = pixel / nx;
-
-            if (!isfinite(image->data.F32[iy][ix]) || (mask && mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & maskBad)) {
-                continue;
-            }
-
-            float value = SN->data.F32[iy][ix];
-            min = PS_MIN(value, min);
-            max = PS_MAX(value, max);
-            values->data.F32[n] = value;
-            n++;
-        }
-    }
-    if (n < 0.01*Nsubset) {
-        psLogMsg("psLib.psImageBackground", PS_LOG_DETAIL, "Unable to measure image background: too few data points (%ld)", n);
-        psFree(values);
-        return false;
-    }
-    values->n = n;
-    psFree (SN);
-
-    psLogMsg("psLib.psImageBackground", PS_LOG_DETAIL, "Using for %ld pixels in S/N image", values->n);
-
-    if (!psVectorStats (stats, values, NULL, NULL, 0)) {
-        if (psTraceGetLevel("psLib.imageops") >= 5) {
-            FILE *f = fopen ("vector.dat", "w");
-            int fd = fileno(f);
-            p_psVectorPrint (fd, values, "values");
-            fclose (f);
-        }
-        psError(PS_ERR_UNKNOWN, false, "Unable to measure statistics for image background "
-                "(%dx%d, (row0,col0) = (%d,%d)",
-                image->numRows, image->numCols, image->row0, image->col0);
-        psFree(values);
-        return false;
-    }
-    if (psTraceGetLevel("psLib.imageops") >= 6) {
-        FILE *f = fopen ("vector.dat", "w");
-        int fd = fileno(f);
-        p_psVectorPrint (fd, values, "values");
-        fclose (f);
-    }
-    psFree (values);
-    psLogMsg ("ppSub", PS_LOG_INFO, "background stdev %f\n", stats->robustStdev);
-
-    float minValid = psMetadataLookupF32(&mdok, ppSubRecipe, "VARIANCE.RESCALE.MIN.VALID");
-    psAssert (mdok, "missing VARIANCE.RESCALE.MIN.VALID in ppSub recipe");
-
-    float maxValid = psMetadataLookupF32(&mdok, ppSubRecipe, "VARIANCE.RESCALE.MAX.VALID");
-    psAssert (mdok, "missing VARIANCE.RESCALE.MAX.VALID in ppSub recipe");
-
-    psLogMsg("psLib.psImageBackground", PS_LOG_DETAIL, "Stdev of S/N image: %f (mean: %f)", stats->robustStdev, stats->robustMedian);
-
-    // valid range of correction; 0.66 to 1.5 (skip if in range 0.98 to 1.02)
-    if ((stats->robustStdev < minValid) || (stats->robustStdev > maxValid)) {
-        psWarning ("background stdev (%f) is out of allowed range; not correcting\n", stats->robustStdev);
-        // XXX set quality to poor
-        psFree (stats);
-        return true;
-    }
-
-    // valid range of correction; 0.66 to 1.5 (skip if in range 0.98 to 1.02)
-    // if ((stats->robustStdev > 0.98) && (stats->robustStdev > 1.02)) {
-    //  psLogMsg ("ppSub", PS_LOG_INFO, "background stdev (%f) is close to 1.0; do not modify variance\n", stats->robustStdev);
-    //  // XXX set quality to poor
-    //  psFree (stats);
-    //  return true;
-    // }
-
-    float varianceFactor = PS_SQR(stats->robustStdev);
-    psLogMsg ("ppSub", PS_LOG_INFO, "background stdev is not 1.0: multiplying variance by %f\n", varianceFactor);
-    for (int y = 0; y < ny; y++) {
-        for (int x = 0; x < nx; x++) {
-            if (mask && mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskBad) {
-                continue;
-            } 
-            if (!isfinite(image->data.F32[y][x])) {
-                continue;
-            } 
-            if (!isfinite(variance->data.F32[y][x]) || (variance->data.F32[y][x] < 0.0)) {
-                continue;
-            } 
-            variance->data.F32[y][x] *= varianceFactor;
-        }
-    }
-
-    pmFPA *outFPA = outRO->parent->parent->parent;
-    pmHDU *outHDU = outFPA->hdu;        // Output HDU
-
-    char history[80];
-    snprintf (history, 80, "Rescaled variance by %6.4f (stdev by %6.4f)", varianceFactor, stats->robustStdev);
-    psMetadataAddStr(outHDU->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, NULL, history);
-
-    psFree (stats);
-    return true;
+    return pmReadoutVarianceRenormalise(readout, maskBad, num, minValid, maxValid);
 }

@@ -13 +13 @@
 ppSubErrorCodes.c
 ppSubVersionDefinitions.h
+ppSubConvolve