Changeset 26899 for trunk/ppSub/src/ppSubMatchPSFs.c

Timestamp:

Feb 10, 2010, 7:42:02 PM (16 years ago)

Author:

eugene

Message:

updates from eam_branches/20091201

File:

• trunk/ppSub/src/ppSubMatchPSFs.c (modified) (8 diffs)

trunk/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(PS_ERR_UNKNOWN, 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");
+
+    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");
+    // XXX assert on inSourcesRO and refSourcesRO?
+
+    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 = subImagePSF(data, inRO, inSources); // FWHM for input
+    float refFWHM = subImagePSF(data, refRO, refSources); // FWHM for reference
+
+    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(PS_ERR_UNKNOWN, false, "Unable to scale parameters.");
+        return false;
+    }
+
+    return true;
+}
+
 
 bool ppSubMatchPSFs(ppSubData *data)
@@ -71 +203 @@
     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) {
+        if (!detections->allSources) {
             psError(PS_ERR_UNKNOWN, 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
@@ -131 +270 @@
     float rej = psMetadataLookupF32(NULL, recipe, "REJ"); // Rejection threshold
     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
@@ -169 +311 @@
             break;
           default:
-            psErrorStackPrint(stderr, "Invalid value for -convolve");
+            psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Invalid value for -convolve");
             return false;
         }
+    }
+
+    if (!subScaleKernel(data, widths, &size, &footprint)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to scale kernel parameters");
+        return false;
     }
 
@@ -177 +324 @@
     if (threads > 0) {
         pmSubtractionThreadsInit(inRO, refRO);
+    }
+
+    if (inRO->covariance) {
+        psImageCovarianceTruncate(inRO->covariance, COVAR_FRAC);
+    }
+    if (refRO->covariance) {
+        psImageCovarianceTruncate(refRO->covariance, COVAR_FRAC);
     }
 
@@ -183 +337 @@
     if (kernelRO) {
         success = pmSubtractionMatchPrecalc(inConv, refConv, inRO, refRO, kernelRO->analysis,
-                                            stride, kernelErr, 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,
+                                     spacing, threshold, detections->allSources, data->stamps, type, size, order,
                                      widths, orders, inner, ringsOrder, binning, penalty, optimum,
-                                     optWidths, optOrder, optThresh, iter, rej, sysErr, kernelErr, maskVal,
-                                     maskBad, maskPoor, poorFrac, badFrac, subMode);
-    }
+                                     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);
@@ -261 +436 @@
     pmConceptsCopyFPA(refConv->parent->parent->parent, refRO->parent->parent->parent, true, true);
 
+    if (inConv->covariance) {
+        psImageCovarianceTruncate(inConv->covariance, COVAR_FRAC);
+    }
+    if (refConv->covariance) {
+        psImageCovarianceTruncate(refConv->covariance, COVAR_FRAC);
+    }
+
     if (inConv->variance) {
         psImageCovarianceTransfer(inConv->variance, inConv->covariance);