Changeset 26649

Timestamp:

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

Author:

Paul Price

Message:

Adding scaling of kernel based on input PSFs.

Location:

branches/eam_branches/20091201/ppSub/src

Files:

• ppSubCamera.c (modified) (2 diffs)
• ppSubMakePSF.c (modified) (4 diffs)
• ppSubMatchPSFs.c (modified) (3 diffs)
• ppSubReadoutPhotometry.c (modified) (1 diff)

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

@@ -317 +317 @@
 
     // psPhot input
-    if (data->photometry) {
+    if (data->photometry || 1) {
         psphotModelClassInit();        // load implementation-specific models
 
@@ -329 +329 @@
             return false;
         }
-        // specify the number of psphot input images
-        psMetadataAddS32 (config->arguments, PS_LIST_TAIL, "PSPHOT.INPUT.NUM", PS_META_REPLACE, "number of inputs", 1);
+        // specify the number of psphot input images
+        psMetadataAddS32 (config->arguments, PS_LIST_TAIL, "PSPHOT.INPUT.NUM", PS_META_REPLACE, "number of inputs", 1);
         pmFPAfileActivate(config->files, false, "PSPHOT.INPUT");
 

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

@@ -70 +70 @@
         psMetadataRemoveKey(photRO->analysis, "PSPHOT.SOURCES");
     }
+    if (psMetadataLookup(photRO->parent->parent->analysis, "PSPHOT.PSF")) {
+        psMetadataRemoveKey(photRO->parent->parent->analysis, "PSPHOT.PSF");
+    }
 
 # ifdef TESTING
@@ -101 +104 @@
         psErrorStackPrint(stderr, "PSF was not generated");
         psWarning("PSF was not generated --- suspect bad data quality.");
-        ppSubDataQuality(data, psErrorCodeLast(), PPSUB_FILES_PHOT_SUB | PPSUB_FILES_PHOT_INV);
+        ppSubDataQuality(data, psErrorCodeLast(), PPSUB_FILES_PHOT_SUB | PPSUB_FILES_PHOT_INV);
     }
 
@@ -124 +127 @@
 
     if (!outputRO) {
-        pmCell *outputCell  = pmFPAviewThisCell(view, output->fpa); 
-        outputRO = pmReadoutAlloc(outputCell);
-        outputRO->image = psMemIncrRefCounter (inputRO->image);
+        pmCell *outputCell  = pmFPAviewThisCell(view, output->fpa);
+        outputRO = pmReadoutAlloc(outputCell);
+        outputRO->image = psMemIncrRefCounter (inputRO->image);
     }
 
@@ -135 +138 @@
     if (!nRegions) {
         psErrorStackPrint(stderr, "No PSF available");
-        return false;
+        return false;
     }
     psMetadataAddS32 (psfRegions, PS_LIST_TAIL, "PSF.CLUMP.NREGIONS",  PS_META_REPLACE, "psf clump regions", nRegions);
 
     for (int i = 0; i < nRegions; i++) {
-        char fieldName[80];
-        snprintf (fieldName, 80, "PSF.CLUMP.REGION.%03d", i);
+        char fieldName[80];
+        snprintf (fieldName, 80, "PSF.CLUMP.REGION.%03d", i);
         psMetadata *regionMD = psMetadataLookupPtr (&mdok, inputRO->analysis, fieldName);
         if (!regionMD) {
-            psWarning ("missing psf clump region metadata for entry %d", i);
-            continue;
-        }
-        psMetadataAddMetadata (psfRegions, PS_LIST_TAIL, fieldName, PS_META_REPLACE, "psf clump region", regionMD);
+            psWarning ("missing psf clump region metadata for entry %d", i);
+            continue;
+        }
+        psMetadataAddMetadata (psfRegions, PS_LIST_TAIL, fieldName, PS_META_REPLACE, "psf clump region", regionMD);
     }
 

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

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

@@ -78 +78 @@
     // erase the overlays from a previous psphot-related step
     if (pmVisualIsVisual()) {
-        psphotVisualEraseOverlays (1, "all");
+        //      psphotVisualEraseOverlays (1, "all");
     }