Changeset 24557 for branches/eam_branches/20090522/ppSub

Timestamp:

Jun 25, 2009, 2:00:56 PM (17 years ago)

Author:

eugene

Message:

merging changes from head

Location:

branches/eam_branches/20090522

Files:

• . (modified) (1 prop)
• ppSub (modified) (1 prop)
• ppSub/src/ppSubCamera.c (modified) (5 diffs)
• ppSub/src/ppSubMatchPSFs.c (modified) (3 diffs)
• ppSub/src/ppSubReadoutInverse.c (modified) (2 diffs)
• ppSub/src/ppSubReadoutPhotometry.c (modified) (1 diff)

branches/eam_branches/20090522/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/eam_branches/20090522/ppSub/src/ppSubMatchPSFs.c

@@ -21 +21 @@
 
 #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)
@@ -166 +183 @@
             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.");
@@ -172 +194 @@
     }
 
+    // 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/eam_branches/20090522/ppSub/src/ppSubReadoutInverse.c

@@ -23 +23 @@
     invRO->data_exists = invRO->parent->data_exists = invRO->parent->parent->data_exists = 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
-    pmFPA *invFPA = invRO->parent->parent->parent; // Inverse FPA
-
-    pmConceptsCopyFPA(invFPA, outFPA, true, true);
-
-#if 0
-    // Copy astrometry over
-    pmHDU *outHDU = outFPA->hdu;          // Output HDU
     pmChip *invChip = invRO->parent->parent; // Inverse chip
-    psFree(view);
-    if (!outHDU || !inHDU) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to find HDU at FPA level to copy astrometry.");
-        return false;
-    }
-    if (!pmAstromReadWCS(invFPA, invChip, outHDU->header, 1.0)) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to read WCS astrometry from PPSUB.OUTPUT.");
-        return false;
-    }
-#endif
     pmHDU *invHDU = invFPA->hdu;          // Inverse HDU
     if (!pmAstromWriteWCS(invHDU->header, outFPA, outChip, WCS_TOLERANCE)) {
@@ -48 +37 @@
         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/eam_branches/20090522/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