Changeset 35512 for branches/eam_branches/ipp-20130419/pswarp/src/pswarpLoop.c

Timestamp:

May 5, 2013, 2:20:46 PM (13 years ago)

Author:

eugene

Message:

stats is basically working now

File:

• branches/eam_branches/ipp-20130419/pswarp/src/pswarpLoop.c (modified) (9 diffs)

branches/eam_branches/ipp-20130419/pswarp/src/pswarpLoop.c

@@ -36 +36 @@
     if (!pswarpSetMaskBits(config)) {
         psError(psErrorCodeLast(), false, "failed to set mask bits");
-        return NULL;
+        return false;
     }
 
@@ -83 +83 @@
     // XXX pmFPAfileActivate(config->files, true, "PSWARP.OUTPUT.SOURCES");
 
-    // Done with the detector side of things
+    // the loops below load the input data
     pswarpFileActivation(config, detectorFiles, true);
-    pswarpFileActivation(config, independentFiles, true);
+
+    // pswarpFileActivation(config, independentFiles, true);
 
     pmFPAview *view = pmFPAviewAlloc(0);
@@ -95 +96 @@
     if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) {
         psError(psErrorCodeLast(), false, "Unable to read files.");
-        goto DONE;
+        goto FAIL;
     }
 
@@ -107 +108 @@
         if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) {
             psError(psErrorCodeLast(), false, "Unable to read files.");
-            goto DONE;
+            goto FAIL;
         }
 
@@ -116 +117 @@
             if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) {
                 psError(psErrorCodeLast(), false, "Unable to read files.");
-                goto DONE;
+                goto FAIL;
             }
 
@@ -124 +125 @@
                 if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
                     psError(psErrorCodeLast(), false, "Unable to read files.");
-                    goto DONE;
+                    goto FAIL;
                 }
                 if (!readout->data_exists) {
@@ -145 +146 @@
                 if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
                     psError(psErrorCodeLast(), false, "Unable to write files.");
-                    goto DONE;
+                    goto FAIL;
                 }
             }
             if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
                 psError(psErrorCodeLast(), false, "Unable to write files.");
-                goto DONE;
+                goto FAIL;
             }
         }
         if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
             psError(psErrorCodeLast(), false, "Unable to write files.");
-            goto DONE;
+            goto FAIL;
         }
     }
     if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
         psError(psErrorCodeLast(), false, "Unable to write files.");
-        goto DONE;
+        goto FAIL;
     }
 
@@ -167 +168 @@
     pswarpFileActivation(config, independentFiles, false);
 
-    // XXXX turn these on later .. pswarpUpdateStatistics (output);
+    if (!pswarpUpdateStatistics (output->fpa, stats, input->fpa, astrom->fpa, config)) {
+      psError(psErrorCodeLast(), false, "problem generating statistics.");
+      goto FAIL;
+    }
+
     // XXXX turn these on later .. pswarpMakePSF (output);
-
     // XXX have not yet written the output to disk
-
     psFree(view);
-
- DONE:
     return true;
+
+ FAIL:
+    psFree (view);
+    return false;
 }
 
-# if (0)
 // once the output fpa elements have been built, loop over the fpa and generate stats
 // for each readout
-bool pswarpUpdateStatistics (pmFPA *output, pmFPA *input, pmConfig *config)  {
-
-    // output mask bits
-    psImageMaskType maskValue = psMetadataLookupImageMask(&status, recipe, "MASK.OUTPUT");
-    psAssert (status, "MASK.OUTPUT was not defined");
+bool pswarpTransformToTarget (pmFPA *output, pmReadout *input, pmConfig *config)  {
 
     pmFPAview *view = pmFPAviewAlloc(0);
@@ -195 +195 @@
 
         pmCell *cell;
-        while ((cell = pmFPAviewNextCell (view, output->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; }
-
-            // process each of the readouts
-            pmReadout *readout;
-            while ((readout = pmFPAviewNextReadout(view, output->fpa, 1)) != NULL) {
-                if (!readout->data_exists) {
-                    psWarning("No overlap between input and skycell.");
-                    if (stats) {
-                        psMetadataAddS32(stats, PS_LIST_TAIL, "QUALITY", PS_META_REPLACE,
-                                         "No overlap between input and skycell", PSWARP_ERR_NO_OVERLAP);
-                    }
-                    psphotFilesActivate(config, false);
-                    psFree(cells);
-                    psFree(view);
-                    goto DONE;
-                }
-    
-                if (!pswarpPixelsLit(readout, stats, config)) {
-                    psError(psErrorCodeLast(), false, "Unable to calculate pixel regions.");
-                    psFree(view);
-                    goto DONE;
-                }
-                bool doStats = psMetadataLookupBool(&mdok,recipe,"MASK.STATS");
-                if (doStats) {
-                    if (!pswarpMaskStats(readout, stats, config)) {
-                        psError(psErrorCodeLast(), false, "Unable to calculate mask stats.");
-                        psFree(cells);
-                        psFree(view);
-                        goto DONE;
-                    }
-                }
-
-                // Set covariance matrix for output
-                {
-                    psList *covariances = psMetadataLookupPtr(&mdok, output->analysis,
-                                                              PSWARP_ANALYSIS_COVARIANCES); // Covariance matrices
-                    psAssert(covariances, "Should be there");
-                    psArray *covars = psListToArray(covariances); // Array of covariance matrices
-                    psKernel *covar = psImageCovarianceAverage(covars);
-                    psFree(covars);
-                    psMetadataRemoveKey(output->analysis, PSWARP_ANALYSIS_COVARIANCES);
-
-                    // Correct covariance matrix scale for the mean (square root of the) Jacobian
-                    double jacobian = psMetadataLookupF64(NULL, output->analysis, PSWARP_ANALYSIS_JACOBIAN); // Jacobian
-                    int goodPixels = psMetadataLookupS32(NULL, output->analysis, PSWARP_ANALYSIS_GOODPIX);   // Good pixels
-                    jacobian /= goodPixels;
-                    output->covariance = psImageCovarianceScale(covar, jacobian);
-                    psFree(covar);
-
-                    if (output->variance) {
-                        psImageCovarianceTransfer(output->variance, output->covariance);
-                    }
-                }
-            }
-
-            // XXX need to get the input Cells for this output cell (save on the analysis somewhere
-            // XXX skip this for now
-# if (0)
-            psArray *inputCells = psMetadataLookupPtr (NULL, cell->analysis, "something");
-            if (!pmConceptsAverageCells(cell, inputCells, NULL, NULL, false)) {
-                psError(psErrorCodeLast(), false, "Unable to average cell concepts.");
-                psFree(view);
-                goto DONE;
-            }
-# endif
-            // XXX Is this too ad-hoc?
-            psRegion *trimsec = psMetadataLookupPtr(NULL, cell->concepts, "CELL.TRIMSEC"); ///< Trim section
-            trimsec->x0 = trimsec->x1 = trimsec->y0 = trimsec->y1 = 0; ///< All pixels
-
-        }
-
-        // Perform statistics on the output image
-        if (stats) {
-            // XXX this should be done at each level...
-            if (!ppStatsFPA(stats, output->parent->parent->parent, view, maskValue, config)) {
-                psWarning("Unable to perform statistics on warped image.");
-            }
-        }
-    }
-    if (!psMetadataCopy(output->concepts, input->concepts)) {
-        psError(psErrorCodeLast(), false, "Unable to copy FPA concepts from input to output.");
-        psFree(view);
-        goto DONE;
-    }
-
-    // Update ZP from the astrometry
-    {
-        psMetadataItem *item = psMetadataLookup(outFPA->concepts, "FPA.ZP");
-        item->data.F32 = psMetadataLookupF32(NULL, astrom->fpa->concepts, "FPA.ZP");
-    }
-
-    // need to update the chip and fpa level astrometry appropriate, see
-    // code in psastro
-    
-# if (0)
-    // Add MD5 information for readout
-    const char *chipName = psMetadataLookupStr(NULL, output->parent->parent->concepts, "CHIP.NAME");
-    const char *cellName = psMetadataLookupStr(NULL, output->parent->concepts, "CELL.NAME");
-    psString headerName = NULL; ///< Header name for MD5
-    psStringAppend(&headerName, "MD5_%s_%s_%d", chipName, cellName, view->readout);
-    psVector *md5 = psImageMD5(output->image); ///< md5 hash
-    psString md5string = psMD5toString(md5); ///< String
-    psFree(md5);
-    psMetadataAddStr(hdu->header, PS_LIST_TAIL, headerName, PS_META_REPLACE,
-                     "Image MD5", md5string);
-    psFree(md5string);
-    psFree(headerName);
-# endif
-
-# if (0)
-    if (!pmAstromWriteWCS(hdu->header, outFPA, outChip, WCS_NONLIN_TOL)) {
-        psError(psErrorCodeLast(), false, "Unable to generate WCS header.");
-        goto DONE;
-    }
-# endif
-    return true;
-}
-# endif
-
-# if (0)
-// 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.
-bool pswarpMakePSF (pmConfig *config) {
-
-    // load the recipe
-    psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSWARP_RECIPE);
-    if (!recipe) {
-        psError(PSWARP_ERR_CONFIG, false, "missing recipe %s", PSWARP_RECIPE);
-        return false;
-    }
-
-    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
-        pmFPAfile *photFile = psMetadataLookupPtr(NULL, config->files, "PSPHOT.INPUT");
-        pmFPACopy(photFile->fpa, outFPA);
-
-        pmFPAview *view = pmFPAviewAlloc(0); ///< View into skycell
-        view->chip = view->cell = view->readout = 0;
-
-        // grab the sources of interest from the storage location (pmFPAfile PSPHOT.INPUT.CMF)
-        psArray *sources = psphotLoadPSFSources (config, view);
-        if (!sources) {
-            psError(psErrorCodeLast(), false, "No sources supplied to measure PSF");
-            goto DONE;
-        }
-
-        pmModelClassSetLimits(PM_MODEL_LIMITS_STRICT);
-
-        // measure the PSF using these sources
-        if (!psphotReadoutFindPSF(config, view, "PSPHOT.INPUT", sources)) {
-            // This is likely a data quality issue
-            // XXX Split into multiple cases using error codes?
-            psErrorStackPrint(stderr, "Unable to determine PSF");
-            psWarning("Unable to determine PSF --- suspect bad data quality.");
-            if (stats && psMetadataLookupS32(NULL, stats, "QUALITY") == 0) {
-                psMetadataAddS32(stats, PS_LIST_TAIL, "QUALITY", PS_META_REPLACE,
-                                 "Unable to determine PSF", psErrorCodeLast());
-            }
-            psErrorClear();
-            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
-#if (TESTING)
-        {
-            #define PSF_SIZE 20         ///< Half-size of PSF
-            #define PSF_FLUX 10000      ///< Central flux for PSF
-            pmChip *photChip = pmFPAviewThisChip(view, photFile->fpa);
-            pmPSF *psf = psMetadataLookupPtr(NULL, photChip->analysis, "PSPHOT.PSF");
-            psImage *image = psImageAlloc(2 * PSF_SIZE + 1, 2 * PSF_SIZE + 1, PS_TYPE_F32);
-            psImageInit(image, 0);
-            pmModel *model = pmModelFromPSFforXY(psf, PSF_SIZE, PSF_SIZE, PSF_FLUX);
-            pmModelAdd(image, NULL, model, PM_MODEL_OP_FULL, 0);
-            psFree(model);
-            psFits *fits = psFitsOpen("psf.fits", "w");
-            psFitsWriteImage(fits, NULL, image, 0, NULL);
-            psFitsClose(fits);
-            psFree(image);
-        }
-#endif
-
-        psFree(view);
-    }
-    return true;
-}
-# endif
-
-// once the output fpa elements have been built, loop over the fpa and generate stats
-// for each readout
-bool pswarpTransformToTarget (pmFPA *output, pmReadout *input, pmConfig *config)  {
-
-    pmFPAview *view = pmFPAviewAlloc(0);
-    
-    pmChip *chip;
-    while ((chip = pmFPAviewNextChip (view, output, 1)) != NULL) {
-        psTrace ("pswarp", 4, "Chip %d: %x %x\n", view->chip, chip->file_exists, chip->process);
-        if (!chip->process || !chip->file_exists) { continue; }
-
-        pmCell *cell;
         while ((cell = pmFPAviewNextCell (view, output, 1)) != NULL) {
             psTrace ("pswarp", 4, "Cell %d: %x %x\n", view->cell, cell->file_exists, cell->process);