Changeset 25027 for branches/pap/ppSub/src/ppSubMatchPSFs.c

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/ppSubMatchPSFs.c (modified) (7 diffs)

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);