Changeset 27839 for branches/simmosaic_branches/pswarp/src/pswarpLoop.c

Timestamp:

May 3, 2010, 8:45:22 AM (16 years ago)

Author:

eugene

Message:

updates from trunk

Location:

branches/simmosaic_branches

Files:

• . (modified) (1 prop)
• pswarp/src/pswarpLoop.c (modified) (20 diffs)

branches/simmosaic_branches/pswarp/src/pswarpLoop.c

@@ -16 +16 @@
 
 #define WCS_NONLIN_TOL 0.001            // Non-linear tolerance for header WCS
-#define PSPHOT_FIND_PSF 1               // Use psphot's findPSF function?
 #define TESTING 0                       // Testing output?
 
 
-// XXX these are generic functions which should be moved to psModules
-// Activate a list of files
-static void fileActivation(pmConfig *config, // Configuration
-                           char **files, // Files to turn on/off
-                           bool state   // Activation state
-    )
-{
-    for (int i = 0; files[i] != NULL; i++) {
-        pmFPAfileActivate(config->files, state, files[i]);
-    }
-    return;
-}
-
-/**
- * Run down the FPA hierarchy, checking files
- */
-static bool ioChecksBefore(pmConfig *config ///< Configuration
-                           )
-{
-    pmFPAview *view = pmFPAviewAlloc(0); ///< View for checking
-    bool status = true;                 ///< Status of checks
-    status &= pmFPAfileIOChecks(config, view, PM_FPA_BEFORE);
-    view->chip = 0;
-    status &= pmFPAfileIOChecks(config, view, PM_FPA_BEFORE);
-    view->cell = 0;
-    status &= pmFPAfileIOChecks(config, view, PM_FPA_BEFORE);
-    psFree(view);
-    return status;
-}
-
-/**
- * Run up the FPA hierarchy, checking files
- */
-static bool ioChecksAfter(pmConfig *config ///< Configuration
-                          )
-{
-    pmFPAview *view = pmFPAviewAlloc(0); ///< View for checking
-    view->chip = view->cell = 0;
-    bool status = true;                 ///< Status of checks
-    status &= pmFPAfileIOChecks(config, view, PM_FPA_AFTER);
-    view->cell = -1;
-    status &= pmFPAfileIOChecks(config, view, PM_FPA_AFTER);
-    view->chip = -1;
-    status &= pmFPAfileIOChecks(config, view, PM_FPA_AFTER);
-    psFree(view);
-    return status;
-}
-
-
-/**
- * Loop over the inputs, warp them to the output skycell and then write out the output.
- */
-bool pswarpLoop(pmConfig *config)
+
+// Loop over the inputs, warp them to the output skycell and then write out the output.
+bool pswarpLoop(pmConfig *config, psMetadata *stats)
 {
     bool status;
+    bool mdok;                          // Status of MD lookup
 
     const char *skyCamera = psMetadataLookupStr(NULL, config->arguments,
-                                                "SKYCELL.CAMERA");  ///< Name of camera for skycell
+                                                "SKYCELL.CAMERA");  // Name of camera for skycell
     pmConfigCamerasCull(config, skyCamera);
     pmConfigRecipesCull(config, "PSWARP,PPSTATS,PSPHOT,PSASTRO,MASKS,JPEG");
-
 
     // load the recipe
     psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSWARP_RECIPE);
     if (!recipe) {
-        psError(PSPHOT_ERR_CONFIG, false, "missing recipe %s", PSWARP_RECIPE);
+        psError(PSWARP_ERR_CONFIG, false, "missing recipe %s", PSWARP_RECIPE);
         return false;
     }
 
     if (!pswarpSetMaskBits(config)) {
-        psError(PS_ERR_IO, false, "failed to set mask bits");
+        psError(psErrorCodeLast(), false, "failed to set mask bits");
         return NULL;
     }
@@ -112 +61 @@
 
     if (astrom->camera != input->camera) {
-        psError(PS_ERR_UNKNOWN, true, "Input camera and astrometry camera do not match.");
+        psError(PSWARP_ERR_DATA, true, "Input camera and astrometry camera do not match.");
         return false;
     }
@@ -128 +77 @@
     psFree (view);
 
-    bool mdok;                          ///< Status of MD lookup
-    const char *statsName = psMetadataLookupStr(&mdok, config->arguments, "STATS"); ///< Filename for statistics
-    psMetadata *stats = NULL;           ///< Container for statistics
-    FILE *statsFile = NULL;             ///< File stream for statistics
-    if (mdok && statsName && strlen(statsName) > 0) {
-        psString resolved = pmConfigConvertFilename(statsName, config, true, true);
-        statsFile = fopen(resolved, "w");
-        if (!statsFile) {
-            psError(PS_ERR_IO, true, "Unable to open statistics file %s for writing.\n", resolved);
-            psFree(resolved);
-            return false;
-        }
-        psFree(resolved);
-        stats = psMetadataAlloc();
-        psMetadataAddS32(stats, PS_LIST_TAIL, "QUALITY", 0, "No problems", 0);
-    }
 
     // Turn all skycell files on to generate them, and then turn them off for the loop over the input images
     // the input, which is in a different format.
     {
-        fileActivation(config, detectorFiles, false);
-        fileActivation(config, photFiles, false);
-        fileActivation(config, independentFiles, false);
-        fileActivation(config, skycellFiles, true);
-        ioChecksBefore(config);
-        fileActivation(config, skycellFiles, false);
+        pswarpFileActivation(config, detectorFiles, false);
+        pswarpFileActivation(config, photFiles, false);
+        pswarpFileActivation(config, independentFiles, false);
+        pswarpFileActivation(config, skycellFiles, true);
+        if (!pswarpIOChecksBefore(config)) {
+            psError(psErrorCodeLast(), false, "Unable to read files.");
+            goto DONE;
+        }
+        pswarpFileActivation(config, skycellFiles, false);
     }
 
@@ -163 +99 @@
         pmChip *chip;
         pmFPAview *view = pmFPAviewAlloc(0);
-        pmFPAfileIOChecks (config, view, PM_FPA_BEFORE);
+        if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
+            psError(psErrorCodeLast(), false, "Unable to read files.");
+            goto DONE;
+        }
         while ((chip = pmFPAviewNextChip (view, input->fpa, 1)) != NULL) {
             psTrace ("pswarp", 4, "Chip %d: %x %x\n", view->chip, chip->file_exists, chip->process);
             if (!chip->process || !chip->file_exists) { continue; }
-            pmFPAfileIOChecks(config, view, PM_FPA_BEFORE);
+            if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
+                psError(psErrorCodeLast(), false, "Unable to read files.");
+                goto DONE;
+            }
             pmCell *cell;
             while ((cell = pmFPAviewNextCell (view, input->fpa, 1)) != NULL) {
                 psTrace ("pswarp", 4, "Cell %d: %x %x\n", view->cell, cell->file_exists, cell->process);
                 if (!cell->process || !cell->file_exists) { continue; }
-                pmFPAfileIOChecks (config, view, PM_FPA_BEFORE);
-                pmFPAfileIOChecks (config, view, PM_FPA_AFTER);
-            }
-            pmFPAfileIOChecks (config, view, PM_FPA_AFTER);
-        }
-        pmFPAfileIOChecks (config, view, PM_FPA_AFTER);
-        psFree(view);
-
-        fileActivation(config, detectorFiles, true);
+                if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE) ||
+                    !pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) {
+                    psError(psErrorCodeLast(), false, "Unable to read files.");
+                    goto DONE;
+                }
+            }
+            if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) {
+                psError(psErrorCodeLast(), false, "Unable to write files.");
+                goto DONE;
+            }
+        }
+        if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) {
+            psError(psErrorCodeLast(), false, "Unable to write files.");
+            goto DONE;
+        }
+        psFree(view);
+
+        pswarpFileActivation(config, detectorFiles, true);
         pmFPAfileActivate(config->files, false, "PSWARP.ASTROM");
     }
@@ -202 +153 @@
     if (bilevelAstrometry) {
         if (!pmAstromReadBilevelMosaic(input->fpa, phu->header)) {
-            psError(PS_ERR_UNKNOWN, false, "Unable to read bilevel mosaic astrometry for input FPA.");
+            psError(psErrorCodeLast(), false, "Unable to read bilevel mosaic astrometry for input FPA.");
             psFree(view);
             psFree(stats);
-            return false;
+            goto DONE;
         }
     }
@@ -212 +163 @@
 
     // files associated with the science image
-    pmFPAfileIOChecks (config, view, PM_FPA_BEFORE);
+    if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) {
+        psError(psErrorCodeLast(), false, "Unable to read files.");
+        goto DONE;
+    }
 
     pmChip *chip;
@@ -218 +172 @@
         psTrace ("pswarp", 4, "Chip %d: %x %x\n", view->chip, chip->file_exists, chip->process);
         if (!chip->process || !chip->file_exists) { continue; }
-        pmFPAfileIOChecks (config, view, PM_FPA_BEFORE);
+        if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) {
+            psError(psErrorCodeLast(), false, "Unable to read files.");
+            goto DONE;
+        }
 
         // read WCS data from the corresponding header
@@ -224 +181 @@
         if (bilevelAstrometry) {
             if (!pmAstromReadBilevelChip (chip, hdu->header)) {
-                psError(PS_ERR_UNKNOWN, false, "Unable to read bilevel chip astrometry for input FPA.");
+                psError(psErrorCodeLast(), false, "Unable to read bilevel chip astrometry for input FPA.");
                 psFree(view);
                 psFree(stats);
-                return false;
+                goto DONE;
             }
         } else {
             // we use a default FPA pixel scale of 1.0
             if (!pmAstromReadWCS (input->fpa, chip, hdu->header, 1.0)) {
-                psError(PS_ERR_UNKNOWN, false, "Unable to read WCS astrometry for input FPA.");
+                psError(psErrorCodeLast(), false, "Unable to read WCS astrometry for input FPA.");
                 psFree(view);
                 psFree(stats);
-                return false;
+                goto DONE;
             }
         }
@@ -243 +200 @@
             psTrace ("pswarp", 4, "Cell %d: %x %x\n", view->cell, cell->file_exists, cell->process);
             if (!cell->process || !cell->file_exists) { continue; }
-            pmFPAfileIOChecks (config, view, PM_FPA_BEFORE);
+            if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) {
+                psError(psErrorCodeLast(), false, "Unable to read files.");
+                goto DONE;
+            }
 
             psListAdd(cells, PS_LIST_TAIL, cell);
@@ -250 +210 @@
             pmReadout *readout;
             while ((readout = pmFPAviewNextReadout(view, input->fpa, 1)) != NULL) {
-                pmFPAfileIOChecks(config, view, PM_FPA_BEFORE);
+                if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
+                    psError(psErrorCodeLast(), false, "Unable to read files.");
+                    goto DONE;
+                }
                 if (!readout->data_exists) {
                     continue;
                 }
 
-                // Copy the sources from the astrometry carrier to the input, so they can be accessed by
+                // Copy the detections from the astrometry carrier to the input, so they can be accessed by
                 // pswarpTransformReadout
                 pmReadout *astromRO = pmFPAviewThisReadout(view, astrom->fpa); // Readout for astrometry
-                psArray *sources = psMetadataLookupPtr(&mdok, astromRO->analysis,
-                                                       "PSPHOT.SOURCES"); // Sources from astrometry
-                if (sources) {
-                    psMetadataAddPtr(readout->analysis, PS_LIST_TAIL, "PSPHOT.SOURCES", PS_DATA_ARRAY,
-                                     "Sources from input astrometry", sources);
+                pmDetections *detections = psMetadataLookupPtr(&mdok, astromRO->analysis, "PSPHOT.DETECTIONS"); // Sources from astrometry
+                if (detections) {
+                    psMetadataAddPtr(readout->analysis, PS_LIST_TAIL, "PSPHOT.DETECTIONS", PS_DATA_ARRAY, "Sources from input astrometry", detections);
                 }
 
                 pswarpTransformReadout(output, readout, config);
 
-                pmFPAfileIOChecks(config, view, PM_FPA_AFTER);
-            }
-            pmFPAfileIOChecks(config, view, PM_FPA_AFTER);
-        }
-        pmFPAfileIOChecks(config, view, PM_FPA_AFTER);
+                if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
+                    psError(psErrorCodeLast(), false, "Unable to write files.");
+                    goto DONE;
+                }
+            }
+            if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
+                psError(psErrorCodeLast(), false, "Unable to write files.");
+                goto DONE;
+            }
+        }
+        if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
+            psError(psErrorCodeLast(), false, "Unable to write files.");
+            goto DONE;
+        }
     }
 
@@ -291 +261 @@
 
     if (!pswarpPixelsLit(output, stats, config)) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to calculate pixel regions.");
+        psError(psErrorCodeLast(), false, "Unable to calculate pixel regions.");
         psFree(cells);
         psFree(view);
-        return false;
+        goto DONE;
     }
 
@@ -314 +284 @@
     }
 
-
-    // Set variance factor
-    {
-        bool mdok;                      // Status of MD lookup
-        float varFactor = psMetadataLookupF32(&mdok, output->analysis, PSWARP_ANALYSIS_VARFACTOR);
-        long goodPix = psMetadataLookupS64(&mdok, output->analysis, PSWARP_ANALYSIS_GOODPIX);
-        varFactor /= goodPix;
-
-        psMetadataItem *vfItem = psMetadataLookup(outCell->concepts, "CELL.VARFACTOR"); ///< Item to update
-        psAssert(vfItem && vfItem->type == PS_TYPE_F32, "Concept should be as specified.");
-        if (!isfinite(vfItem->data.F32)) {
-            vfItem->data.F32 = varFactor;
-        } else {
-            vfItem->data.F32 *= varFactor;
-        }
-    }
 
     // Set covariance matrix for output
@@ -344 +298 @@
 
     if (!pmConceptsAverageCells(outCell, cells, NULL, NULL, false)) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to average cell concepts.");
+        psError(psErrorCodeLast(), false, "Unable to average cell concepts.");
         psFree(stats);
         psFree(cells);
         psFree(view);
-        return false;
+        goto DONE;
     }
     psFree(cells);
@@ -356 +310 @@
 
     if (!psMetadataCopy(outFPA->concepts, input->fpa->concepts)) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to copy FPA concepts from input to output.");
+        psError(psErrorCodeLast(), false, "Unable to copy FPA concepts from input to output.");
         psFree(stats);
         psFree(view);
-        return false;
+        goto DONE;
+    }
+
+    // Update ZP from the astrometry
+    {
+        psMetadataItem *item = psMetadataLookup(outFPA->concepts, "FPA.ZP");
+        item->data.F32 = psMetadataLookupF32(NULL, astrom->fpa->concepts, "FPA.ZP");
     }
 
@@ -372 +332 @@
         pmHDU *skyHDU = pmHDUFromCell(cell); ///< HDU
         if (!skyHDU) {
-            psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to find skycell HDU.");
+            psError(PSWARP_ERR_DATA, false, "Unable to find skycell HDU.");
             psFree(view);
-            return false;
+            goto DONE;
         }
         hdu->header = psMetadataCopy(hdu->header, skyHDU->header);
@@ -382 +342 @@
 
     if (!pmAstromWriteWCS(hdu->header, outFPA, outChip, WCS_NONLIN_TOL)) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to generate WCS header.");
+        psError(psErrorCodeLast(), false, "Unable to generate WCS header.");
         psFree(stats);
-        return false;
-    }
-
-    pmFPAfileIOChecks(config, view, PM_FPA_AFTER);
+        goto DONE;
+    }
+
+    if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
+        psError(psErrorCodeLast(), false, "Unable to write files.");
+        goto DONE;
+    }
 
     // Done with the detector side of things
-    fileActivation(config, detectorFiles, false);
-    fileActivation(config, independentFiles, false);
+    pswarpFileActivation(config, detectorFiles, false);
+    pswarpFileActivation(config, independentFiles, false);
 
     // We need a new PSF model for the warped frame.  It would be good to generate this analytically, but
     // that's going to be tricky.  We have a list of sources, so we use those to redetermine the PSF model.
 
-    if (psMetadataLookupBool(&mdok, config->arguments, "PSF")) {
-        fileActivation(config, photFiles, true);
-        ioChecksBefore(config);
+    if (psMetadataLookupBool(&mdok, recipe, "PSF")) {
+        pswarpFileActivation(config, photFiles, true);
+        if (!pswarpIOChecksBefore(config)) {
+            psError(psErrorCodeLast(), false, "Unable to read files.");
+            goto DONE;
+        }
 
         // supply the readout and fpa of interest to psphot
@@ -407 +373 @@
         view->chip = view->cell = view->readout = 0;
 
-        // Need to adjust the weight --- the main operation in psphotReadoutFindPSF is not determining the
-        // signficance of sources, but looking at the significance of individual pixels.
-        // We can adjust the weight directly since this is a deep copy
-        pmReadout *photRO = pmFPAviewThisReadout(view, photFile->fpa);
-        float vf = psMetadataLookupF32(NULL, photRO->parent->concepts, "CELL.VARFACTOR"); // Variance factor
-        psBinaryOp(photRO->variance, photRO->variance, "*", psScalarAlloc(vf, PS_TYPE_F32));
-
         // grab the sources of interest from the storage location (pmFPAfile PSPHOT.INPUT.CMF)
         psArray *sources = psphotLoadPSFSources (config, view);
         if (!sources) {
-            psError(PS_ERR_UNKNOWN, false, "No sources supplied to measure PSF");
-            return false;
-        }
+            psError(psErrorCodeLast(), false, "No sources supplied to measure PSF");
+            goto DONE;
+        }
+
+        pmModelClassSetLimits(PM_MODEL_LIMITS_STRICT);
 
         // measure the PSF using these sources
@@ -434 +395 @@
             psphotFilesActivate(config, false);
         }
+
+        // Ensure seeing is carried over
+        pmChip *photChip = pmFPAviewThisChip(view, photFile->fpa);                 // Chip with seeing
+        psMetadataItem *item = psMetadataLookup(outChip->concepts, "CHIP.SEEING"); // Concept with seeing
+        item->data.F32 = psMetadataLookupF32(NULL, photChip->concepts, "CHIP.SEEING");
 
 // XXX EAM : put this in a visualization function
@@ -457 +423 @@
     }
 
-    fileActivation(config, skycellFiles, true);
-
     // Perform statistics on the output image
     if (stats) {
@@ -480 +444 @@
     psFree(view);
 
-    // Ensure everything is written out, at every level
-    ioChecksAfter(config);
-
-    // Now done with the skycell side of things
-
  DONE:
-
-    // Write out summary statistics
-    if (stats) {
-        psMetadataAddF32(stats, PS_LIST_TAIL, "DT_WARP", 0, "Time for warp completion",
-                         psTimerMark("pswarp"));
-
-        const char *statsMDC = psMetadataConfigFormat(stats);
-        if (!statsMDC) {
-            psError(PS_ERR_IO, false, "Unable to get statistics MDC file.\n");
-            psFree(stats);
-            fclose(statsFile);
-            return false;
-        }
-        fprintf(statsFile, "%s", statsMDC);
-        psFree((void*)statsMDC);
-        fclose(statsFile);
-        pmConfigRunFilenameAddWrite(config, "STATS", statsName);
-
-        psFree(stats);
-    }
-
-    // Dump configuration
-    psString dump_file = psMetadataLookupStr(&status, config->arguments, "DUMP_CONFIG");
-    if (dump_file) {
-        pmConfigDump(config, dump_file);
-    }
 
     return true;