Changeset 25027 for branches/pap/ppSub

Timestamp:

Aug 7, 2009, 4:08:25 PM (17 years ago)

Author:

Paul Price

Message:

Merging trunk (r25026) to get up-to-date on old branch.

Location:

branches/pap

Files:

• . (modified) (1 prop)
• ppSub (modified) (1 prop)
• ppSub/src/Makefile.am (modified) (2 diffs)
• ppSub/src/ppSub.h (modified) (3 diffs)
• ppSub/src/ppSubArguments.c (modified) (1 diff)
• ppSub/src/ppSubCamera.c (modified) (5 diffs)
• ppSub/src/ppSubKernel.c (modified) (1 diff)
• ppSub/src/ppSubLoop.c (modified) (2 diffs)
• ppSub/src/ppSubMatchPSFs.c (modified) (7 diffs)
• ppSub/src/ppSubReadout.c (deleted)
• ppSub/src/ppSubReadoutInverse.c (modified) (2 diffs)
• ppSub/src/ppSubReadoutPhotometry.c (modified) (2 diffs)
• ppSub/src/ppSubReadoutSubtract.c (modified) (2 diffs)
• ppSub/src/ppSubThreshold.c (copied) (copied from trunk/ppSub/src/ppSubThreshold.c )
• ppSub/src/ppSubVarianceRescale.c (copied) (copied from trunk/ppSub/src/ppSubVarianceRescale.c )

branches/pap/ppSub/src/Makefile.am

@@ -30 +30 @@
         ppSubVersion.c                  \
         ppSubBackground.c               \
+        ppSubVarianceRescale.c          \
         ppSubCamera.c                   \
         ppSubData.c                     \
@@ -44 +45 @@
         ppSubReadoutStats.c             \
         ppSubReadoutSubtract.c          \
-        ppSubSetMasks.c
+        ppSubSetMasks.c                 \
+        ppSubThreshold.c
 
 ppSubKernel_CPPFLAGS = $(PSLIB_CFLAGS) $(PSMODULE_CFLAGS) $(PPSUB_CFLAGS)

branches/pap/ppSub/src/ppSub.h

@@ -24 +24 @@
 
 #define PPSUB_RECIPE "PPSUB"            /// Name of the recipe to use
+#define WCS_TOLERANCE 0.001             // Tolerance for WCS
 
 // Output files, for activation/deactivation
@@ -77 +78 @@
     );
 
+/// Threshold low pixels in image
+bool ppSubLowThreshold(ppSubData *data  ///< Processing data
+    );
+
 /// Generate the output readout and pass the kernel info to the header
 bool ppSubDefineOutput(const char *name,///< Name of output to define
@@ -98 +103 @@
 /// Higher-order background subtraction
 bool ppSubBackground(pmConfig *config   ///< Configuration
+    );
+
+/// Perform Variance correction (rescale within a modest range)
+bool ppSubVarianceRescale(pmConfig *config   ///< Configuration
     );
 

branches/pap/ppSub/src/ppSubArguments.c

@@ -84 +84 @@
     psMetadataAddS32(arguments, PS_LIST_TAIL, "-threads", 0, "Number of threads", 0);
     psMetadataAddStr(arguments, PS_LIST_TAIL, "-dumpconfig", 0, "file to dump configuration to", NULL);
+    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);

branches/pap/ppSub/src/ppSubCamera.c

@@ -163 +163 @@
 
 
+    // 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);
+        return false;
+    }
+    if (psMetadataLookupBool(NULL, config->arguments, "-photometry")) {
+        psMetadataAddBool(recipe, PS_LIST_TAIL, "PHOTOMETRY", PS_META_REPLACE,
+                          "Perform photometry?", true);
+    }
+    if (psMetadataLookupBool(NULL, config->arguments, "-inverse")) {
+        psMetadataAddBool(recipe, PS_LIST_TAIL, "INVERSE", PS_META_REPLACE,
+                          "Generate inverse subtractions?", true);
+    }
+
+    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
     pmFPAfile *inConvImage = defineOutputFile(config, input, true, "PPSUB.INPUT.CONV", PM_FPA_FILE_IMAGE);
@@ -171 +192 @@
         return false;
     }
+    if (saveConv) {
+        inConvImage->save = true;
+        inConvMask->save = true;
+    }
     if (inVar) {
         pmFPAfile *inConvVar = defineOutputFile(config, inConvImage, false, "PPSUB.INPUT.CONV.VARIANCE",
@@ -177 +202 @@
             psError(PS_ERR_UNKNOWN, false, "Unable to define output files");
             return false;
+        }
+        if (saveConv) {
+            inConvVar->save = true;
         }
     }
@@ -188 +216 @@
         return false;
     }
+    if (saveConv) {
+        refConvImage->save = true;
+        refConvMask->save = true;
+    }
     if (refVar) {
         pmFPAfile *refConvVar = defineOutputFile(config, refConvImage, false, "PPSUB.REF.CONV.VARIANCE",
@@ -195 +227 @@
             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);
-        return false;
-    }
-    if (psMetadataLookupBool(NULL, config->arguments, "-photometry")) {
-        psMetadataAddBool(recipe, PS_LIST_TAIL, "PHOTOMETRY", PS_META_REPLACE,
-                          "Perform photometry?", true);
-    }
-    if (psMetadataLookupBool(NULL, config->arguments, "-inverse")) {
-        psMetadataAddBool(recipe, PS_LIST_TAIL, "INVERSE", PS_META_REPLACE,
-                          "Generate inverse subtractions?", true);
-    }
-
-    data->inverse = psMetadataLookupBool(NULL, recipe, "INVERSE");
-    data->photometry = psMetadataLookupBool(NULL, recipe, "PHOTOMETRY");
-
+        if (saveConv) {
+            refConvVar->save = true;
+        }
+    }
 
     // Output image

branches/pap/ppSub/src/ppSubKernel.c

@@ -29 +29 @@
     }
 
-    psTraceSetLevel("psModules.imcombine", 7);
+    (void) psTraceSetLevel("psModules.imcombine", 7);
 
     const char *inName = argv[1];       // Input file name

branches/pap/ppSub/src/ppSubLoop.c

@@ -65 +65 @@
     }
 
+    if (!ppSubLowThreshold(data)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to threshold images.");
+        return false;
+    }
+
     // Set up subtraction files
     if (!ppSubFilesIterateDown(config, PPSUB_FILES_SUB)) {
@@ -95 +100 @@
     if (!ppSubBackground(config)) {
         psError(PS_ERR_UNKNOWN, false, "Unable to subtract background.");
+        return false;
+    }
+
+    // Perform Variance correction (rescale within a modest range)
+    if (!ppSubVarianceRescale(config)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to rescale variance.");
         return false;
     }

branches/pap/ppSub/src/ppSubMatchPSFs.c

@@ -18 +18 @@
 #include <pslib.h>
 #include <psmodules.h>
-#include <psphot.h>
 
 #include "ppSub.h"
+
+// Normalise a region on an image
+static void normaliseRegion(psImage *image, // Image to normalise
+                            const psRegion *region, // Region of image to normalise
+                            float norm  // Normalisation
+                            )
+{
+    if (!image) {
+        return;
+    }
+    psAssert(region, "Expect region");
+    psImage *subImage = psImageSubset(image, *region); // Sub-image
+    psBinaryOp(subImage, subImage, "*", psScalarAlloc(norm, PS_TYPE_F32));
+    psFree(subImage);
+    return;
+}
+
 
 bool ppSubMatchPSFs(ppSubData *data)
@@ -60 +76 @@
     pmReadout *inSourceRO = pmFPAfileThisReadout(config->files, view, "PPSUB.INPUT.SOURCES");
     pmReadout *refSourceRO = pmFPAfileThisReadout(config->files, view, "PPSUB.REF.SOURCES");
-    psArray *inSources = psMetadataLookupPtr(&mdok, inSourceRO->analysis, "PSPHOT.SOURCES"); // Source list
-    psArray *refSources = psMetadataLookupPtr(&mdok, refSourceRO->analysis, "PSPHOT.SOURCES"); // Source list
+    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
@@ -92 +110 @@
     float spacing = psMetadataLookupF32(NULL, recipe, "STAMP.SPACING"); // Typical stamp spacing
     float threshold = psMetadataLookupF32(NULL, recipe, "STAMP.THRESHOLD"); // Threshold for stmps
-    const char *stampsName = psMetadataLookupStr(&mdok, config->arguments, "STAMPS"); // Filename for stamps
 
     const char *typeStr = psMetadataLookupStr(NULL, recipe, "KERNEL.TYPE"); // Kernel type
@@ -135 +152 @@
 
     bool dual = psMetadataLookupBool(&mdok, recipe, "DUAL"); // Dual convolution?
-    pmSubtractionMode subMode = dual ? PM_SUBTRACTION_MODE_DUAL : PM_SUBTRACTION_MODE_UNSURE; // Subtracn mode
+    pmSubtractionMode subMode;          // Subtraction mode
+    if (dual) {
+        subMode = PM_SUBTRACTION_MODE_DUAL;
+    } else {
+        int convolve = psMetadataLookupS32(NULL, config->arguments, "-convolve"); // Image number to convolve
+        switch (convolve) {
+          case 0:
+            subMode = PM_SUBTRACTION_MODE_UNSURE;
+            break;
+          case 1:
+            subMode = PM_SUBTRACTION_MODE_1;
+            break;
+          case 2:
+            subMode = PM_SUBTRACTION_MODE_2;
+            break;
+          default:
+            psErrorStackPrint(stderr, "Invalid value for -convolve");
+            return false;
+        }
+    }
 
     int threads = psMetadataLookupS32(NULL, config->arguments, "-threads"); // Number of threads
@@ -149 +185 @@
     } else {
         success = pmSubtractionMatch(inConv, refConv, inRO, refRO, footprint, stride, regionSize,
-                                     spacing, threshold, sources, stampsName, type, size, order,
+                                     spacing, threshold, sources, data->stamps, type, size, order,
                                      widths, orders, inner, ringsOrder, binning, penalty, optimum,
                                      optWidths, optOrder, optThresh, iter, rej, sys, maskVal,
@@ -165 +201 @@
             ppSubDataQuality(data, error, PPSUB_FILES_ALL);
             return true;
+        } else if (error == PM_ERR_SMALL_AREA) {
+            psErrorStackPrint(stderr, "Insufficient area for PSF matching");
+            psWarning("Insufficient area for PSF matching --- suspect bad data quality.");
+            ppSubDataQuality(data, error, PPSUB_FILES_ALL);
+            return true;
         } else {
             psError(PS_ERR_UNKNOWN, false, "Unable to match images.");
@@ -171 +212 @@
     }
 
+    // Need to be careful with the normalisation
+    // We will normalise everything to the normalisation of the *input* image
+    {
+        // Since the entries are MULTI, we have to retrieve them differently
+        psMetadataIterator *regIter = psMetadataIteratorAlloc(inConv->analysis, PS_LIST_HEAD,
+                                                              "^" PM_SUBTRACTION_ANALYSIS_REGION "$");
+        psMetadataIterator *modeIter = psMetadataIteratorAlloc(inConv->analysis, PS_LIST_HEAD,
+                                                              "^" PM_SUBTRACTION_ANALYSIS_MODE "$");
+        psMetadataIterator *normIter = psMetadataIteratorAlloc(inConv->analysis, PS_LIST_HEAD,
+                                                              "^" PM_SUBTRACTION_ANALYSIS_NORM "$");
+        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
+            psMetadataItem *modeItem = psMetadataGetAndIncrement(modeIter); // Item with mode
+            psAssert(modeItem && modeItem->type == PS_TYPE_S32, "Expect subtraction mode");
+            pmSubtractionMode mode = modeItem->data.S32; // Subtraction mode
+            psMetadataItem *normItem = psMetadataGetAndIncrement(normIter); // Item with normalisation
+            psAssert(normItem && normItem->type == PS_TYPE_F32 && isfinite(normItem->data.F32),
+                     "Expect normalisation");
+            float norm = normItem->data.F32; // Normalisation
+
+            switch (mode) {
+              case PM_SUBTRACTION_MODE_1: // Convolved the input to match template
+              case PM_SUBTRACTION_MODE_DUAL: // Convolved both; template should have flux conserved
+                psLogMsg("ppSub", PS_LOG_INFO, "Correcting image for normalisation of %f\n", norm);
+                normaliseRegion(inConv->image, region, 1.0 / norm);
+                normaliseRegion(refConv->image, region, 1.0 / norm);
+                normaliseRegion(inConv->variance, region, 1.0 / PS_SQR(norm));
+                normaliseRegion(refConv->variance, region, 1.0 / PS_SQR(norm));
+                break;
+              case PM_SUBTRACTION_MODE_2:       // Convolved the template to match input
+                // We're already happy!
+                psLogMsg("ppSub", PS_LOG_INFO, "Image normalisation is correct\n");
+                break;
+              default:
+                psAbort("Invalid subtraction mode: %x", mode);
+            }
+        }
+        psFree(regIter);
+        psFree(modeIter);
+        psFree(normIter);
+    }
+
+
     pmConceptsCopyFPA(inConv->parent->parent->parent, inRO->parent->parent->parent, true, true);
     pmConceptsCopyFPA(refConv->parent->parent->parent, refRO->parent->parent->parent, true, true);

branches/pap/ppSub/src/ppSubReadoutInverse.c

@@ -7 +7 @@
 
 #include "ppSub.h"
-
 
 bool ppSubReadoutInverse(pmConfig *config)
@@ -24 +23 @@
     invRO->data_exists = invRO->parent->data_exists = invRO->parent->parent->data_exists = true;
 
-    pmConceptsCopyFPA(invRO->parent->parent->parent, outRO->parent->parent->parent, true, true);
+    // Get concepts from reference
+    pmFPAfile *refFile = psMetadataLookupPtr(NULL, config->files, "PPSUB.REF.CONV"); // File with concepts
+    pmFPA *invFPA = invRO->parent->parent->parent; // Inverse FPA
+    pmConceptsCopyFPA(invFPA, refFile->fpa, true, true);
+
+    // Get astrometry from (forward) subtraction
+    pmChip *outChip = outRO->parent->parent;       // Output chip
+    pmFPA *outFPA = outChip->parent;               // Output FPA
+    pmChip *invChip = invRO->parent->parent; // Inverse chip
+    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.");
+        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.");
+        return false;
+    }
 
     return true;

branches/pap/ppSub/src/ppSubReadoutPhotometry.c

@@ -65 +65 @@
     // 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");
+    }
 
     pmFPAfile *photFile = psMetadataLookupPtr(&mdok, config->files, "PSPHOT.INPUT"); // Photometry file
@@ -92 +95 @@
     }
 
+    // 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) {
+        ppSubDataQuality(data, PSPHOT_ERR_DATA, PPSUB_FILES_PHOT_SUB | PPSUB_FILES_PHOT_INV);
+    }
+
     if (data->stats) {
-        psArray *sources = psMetadataLookupPtr(NULL, photRO->analysis, "PSPHOT.SOURCES"); // Sources
         bool mdok;
         int numSources = psMetadataLookupS32(&mdok, data->stats, "NUM_SOURCES"); // Number of sources

branches/pap/ppSub/src/ppSubReadoutSubtract.c

@@ -20 +20 @@
 
 #include "ppSub.h"
-
-#define WCS_TOLERANCE 0.001             // Tolerance for WCS
 
 bool ppSubReadoutSubtract(pmConfig *config)
@@ -46 +44 @@
         subtrahend = pmFPAfileThisReadout(config->files, view, "PPSUB.REF.CONV");
     }
-
-// XXX this should be replaced by SAVE options to the inConv and refConf pmFPAfiles
-#ifdef TESTING
-    {
-        pmReadoutMaskApply(minuend, maskVal);
-        psFits *fits = psFitsOpen("minuend.fits", "w");
-        psFitsWriteImage(fits, NULL, minuend->image, 0, NULL);
-        psFitsClose(fits);
-    }
-    {
-        pmReadoutMaskApply(subtrahend, maskVal);
-        psFits *fits = psFitsOpen("subtrahend.fits", "w");
-        psFitsWriteImage(fits, NULL, subtrahend->image, 0, NULL);
-            psFitsClose(fits);
-    }
-#endif
 
     // Do the actual subtraction