Changeset 26649 for branches/eam_branches/20091201/ppSub/src/ppSubMatchPSFs.c

Timestamp:

Jan 20, 2010, 5:20:00 PM (16 years ago)

Author:

Paul Price

Message:

Adding scaling of kernel based on input PSFs.

File:

• branches/eam_branches/20091201/ppSub/src/ppSubMatchPSFs.c (modified) (3 diffs)

branches/eam_branches/20091201/ppSub/src/ppSubMatchPSFs.c

@@ -36 +36 @@
     psFree(subImage);
     return;
+}
+
+// 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.SOURCES")) {
+        psMetadataRemoveKey(photRO->analysis, "PSPHOT.SOURCES");
+    }
+    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");
+    psArray *inSources = psMetadataLookupPtr(NULL, inSourceRO->analysis, "PSPHOT.SOURCES"); // Input sources
+    psArray *refSources = psMetadataLookupPtr(NULL, refSourceRO->analysis, "PSPHOT.SOURCES"); // Ref sources
+
+    psFree(view);
+
+    if (!inSources || !refSources) {
+        psWarning("Unable to scale kernel, since no sources were provided.");
+        return true;
+    }
+
+
+    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 diff = sqrtf(PS_SQR(PS_MAX(inFWHM, refFWHM)) - PS_SQR(PS_MIN(inFWHM, refFWHM))); // Difference
+    float diff = PS_MAX(inFWHM, refFWHM);
+
+    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;
+    }
+
+    float scale = diff / scaleRef;      // Scaling factor
+    if (scale < scaleMin) {
+        scale = scaleMin;
+    }
+    if (scale > scaleMax) {
+        scale = scaleMax;
+    }
+
+    for (int i = 0; i < kernelWidths->n; i++) {
+        kernelWidths->data.F32[i] *= scale;
+    }
+    *kernelSize = *kernelSize * scale + 0.5;
+    *stampSize = *stampSize * scale + 0.5;
+
+    psLogMsg("ppSub", PS_LOG_INFO, "Scaled kernel parameters by %f: %d %d", scale, *kernelSize, *stampSize);
+
+    return true;
 }
 
@@ -133 +267 @@
     float kernelErr = psMetadataLookupF32(NULL, recipe, "KERNEL.ERR"); // Relative systematic error in kernel
     float sysErr = psMetadataLookupF32(NULL, recipe, "SYS.ERR"); // Relative systematic error in images
-    float skyErr = psMetadataLookupF32(NULL, recipe, "SKY.ERR"); // Relative systematic error in images
+    float skyErr = psMetadataLookupF32(NULL, recipe, "SKY.ERR"); // Additional error in sky
 
     float badFrac = psMetadataLookupF32(NULL, recipe, "BADFRAC"); // Maximum bad fraction
@@ -171 +305 @@
             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;
     }