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

Timestamp:

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

Author:

eugene

Message:

updates from eam_branches/20091201

File:

• trunk/ppSub/src/ppSubMakePSF.c (modified) (4 diffs)

trunk/ppSub/src/ppSubMakePSF.c

@@ -22 +22 @@
 #include "ppSub.h"
 
+psArray *ppSubSelectPSFSources(psArray *sources);
+
 bool ppSubMakePSF(ppSubData *data)
 {
@@ -61 +63 @@
         return false;
     }
+
     pmReadout *photRO = pmFPAviewThisReadout(view, photFile->fpa); // Readout to photometer
-    if (psMetadataLookup(photRO->analysis, "PSPHOT.SOURCES")) {
-        psMetadataRemoveKey(photRO->analysis, "PSPHOT.SOURCES");
-    }
-
-    // psphotSaveImage (photFile->fpa->hdu->header, photRO->image, "findpsf.im.fits");
-    // psphotSaveImage (photFile->fpa->hdu->header, photRO->variance, "findpsf.wt.fits");
-    // psphotSaveImage (photFile->fpa->hdu->header, photRO->mask, "findpsf.mk.fits");
-
-    // XXX for testing, can i dump a complete copy of the psphot sources here?  
-    // I actually only need the X,Y coordinates to test this.
+
+    // we need to remove any existing PSPHOT.DETECTIONS (why not do this in psphot?)
+    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");
+    }
+
+# ifdef TESTING
+    // XXX for testing, dump these images:
+    psphotSaveImage (NULL, photRO->image, "findpsf.im.fits");
+    psphotSaveImage (NULL, photRO->variance, "findpsf.wt.fits");
+    psphotSaveImage (NULL, photRO->mask, "findpsf.mk.fits");
+# endif
 
     // Extract the loaded sources from the associated readout, and generate PSF
     // Here, we assume the image is background-subtracted
-    psArray *sources = psMetadataLookupPtr(&mdok, minuend->analysis, "PSPHOT.SOURCES");
-    if (!psphotReadoutFindPSF(config, view, sources)) {
+    pmDetections *detections = psMetadataLookupPtr(&mdok, minuend->analysis, "PSPHOT.DETECTIONS");
+    psAssert(detections, "missing detections?");
+    psArray *sources = detections->allSources;
+    psAssert(sources, "missing sources?");
+
+    // XXX filter sources?  limit the total number and return only brighter objects?
+    // use flags to toss totally bogus entries?
+    psArray *goodSources = ppSubSelectPSFSources (sources);
+    if (!psphotReadoutFindPSF(config, view, goodSources)) {
         // This is likely a data quality issue
         // XXX Split into multiple cases using error codes?
@@ -83 +98 @@
         ppSubDataQuality(data, PSPHOT_ERR_PSF, PPSUB_FILES_PHOT_SUB | PPSUB_FILES_PHOT_INV);
         psFree(view);
+        psFree(goodSources);
         return true;
+    }
+
+    // save the resulting PSF information on the pmFPAfile PSPHOT.PSF.LOAD
+    pmFPAfile *psfFile = psMetadataLookupPtr(&mdok, config->files, "PSPHOT.PSF.LOAD"); // PSF file
+    if (!ppSubCopyPSF (psfFile, photFile, view)) {
+        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);
     }
 
@@ -89 +113 @@
     psMetadataRemoveKey(photRO->analysis, "PSPHOT.HEADER");
 
-    if (!data->quality) {
-        data->psf = psMemIncrRefCounter(psMetadataLookupPtr(NULL, photRO->parent->parent->analysis,
-                                                            "PSPHOT.PSF"));
-    }
-
+    psFree(goodSources);
     psFree(view);
 
     return true;
 }
+
+bool ppSubCopyPSF (pmFPAfile *output, pmFPAfile *input, pmFPAview *view) {
+
+    bool mdok = false;
+
+    pmChip *inputChip   = pmFPAviewThisChip(view, input->fpa); // Chip with PSF info
+    pmChip *outputChip  = pmFPAviewThisChip(view, output->fpa); // Chip to store PSF info
+
+    pmReadout *inputRO  = pmFPAviewThisReadout(view, input->fpa); // Readout with PSF info
+    pmReadout *outputRO = pmFPAviewThisReadout(view, output->fpa); // Readout to store PSF info
+
+    if (!outputRO) {
+        pmCell *outputCell  = pmFPAviewThisCell(view, output->fpa);
+        outputRO = pmReadoutAlloc(outputCell);
+        outputRO->image = psMemIncrRefCounter (inputRO->image);
+    }
+
+    // copy the PSF-related data to PSPHOT.PSF.LOAD for safe-keeping
+    psMetadata *psfRegions = psMetadataAlloc();
+
+    int nRegions = psMetadataLookupS32 (&mdok, inputRO->analysis, "PSF.CLUMP.NREGIONS");
+    if (!nRegions) {
+        psErrorStackPrint(stderr, "No PSF available");
+        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);
+        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);
+    }
+
+    // XXX why am I replacing the entire analysis MD?
+    psFree(outputRO->analysis);
+    outputRO->analysis = psfRegions;
+
+    // copy the PSF model data
+    pmPSF *psf = psMetadataLookupPtr(NULL, inputChip->analysis, "PSPHOT.PSF"); // PSF for photometry
+    if (!psf) {
+        psErrorStackPrint(stderr, "No PSF available");
+        return false;
+    }
+
+    psMetadataAddPtr(outputChip->analysis, PS_LIST_TAIL, "PSPHOT.PSF", PS_DATA_UNKNOWN | PS_META_REPLACE, "PSF from ppSubMakePSF", psf);
+
+    return true;
+}
+
+
+// XXX hardwired MAX for now
+# define MAX_NPSF 500
+
+psArray *ppSubSelectPSFSources(psArray *sources){
+
+    sources = psArraySort (sources, pmSourceSortBySN);
+
+    psArray *subset = psArrayAllocEmpty(MAX_NPSF);
+
+    int nPSF = 0;
+    for (int i = 0; (nPSF < MAX_NPSF) && (i < sources->n); i++) {
+
+        pmSource *source = sources->data[i];
+        if (!source) continue;
+
+        // skip non-astronomical objects (very likely defects)
+        if (source->type == PM_SOURCE_TYPE_DEFECT) continue;
+        if (source->type == PM_SOURCE_TYPE_SATURATED) continue;
+        if (source->mode & PM_SOURCE_MODE_SATSTAR) continue;
+
+        psArrayAdd (subset, 100, source);
+        nPSF++;
+    }
+
+    return subset;
+}