Changeset 34411 for branches/czw_branch/20120906/pswarp/src/pswarpLoop.c

Timestamp:

Sep 6, 2012, 1:24:39 PM (14 years ago)

Author:

watersc1

Message:

Broken background warping. Most of the code is in place, still trying to figure out why the results don't look good.

File:

• branches/czw_branch/20120906/pswarp/src/pswarpLoop.c (modified) (5 diffs)

branches/czw_branch/20120906/pswarp/src/pswarpLoop.c

@@ -179 +179 @@
         // read WCS data from the corresponding header
         pmHDU *hdu = pmFPAviewThisHDU (view, astrom->fpa);
+
+        
         if (bilevelAstrometry) {
             if (!pmAstromReadBilevelChip (chip, hdu->header)) {
@@ -196 +198 @@
         }
 
+        
         pmCell *cell;
         while ((cell = pmFPAviewNextCell (view, input->fpa, 1)) != NULL) {
@@ -227 +230 @@
 
                 pswarpTransformReadout(output, readout, config);
-
+                
                 if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
                     psError(psErrorCodeLast(), false, "Unable to write files.");
@@ -255 +258 @@
         goto DONE;
     }
-
+    
     pmCell *outCell = output->parent;   ///< Output cell
     pmChip *outChip = outCell->parent;  ///< Output chip
@@ -450 +453 @@
     return true;
 }
+
+// Loop over the inputs, warp them to the output skycell and then write out the output.
+bool pswarpLoopBackground(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
+    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(PSWARP_ERR_CONFIG, false, "missing recipe %s", PSWARP_RECIPE);
+        return false;
+    }
+
+    if (!pswarpSetMaskBits(config)) {
+        psError(psErrorCodeLast(), false, "failed to set mask bits");
+        return NULL;
+    }
+
+    // select the input data sources
+    pmFPAfile *input = psMetadataLookupPtr(NULL, config->files, "PSWARP.BKGMODEL");
+    if (!input) {
+        psError(PSWARP_ERR_CONFIG, true, "Can't find input data!\n");
+        return false;
+    }
+
+    // use the external astrometry source if supplied
+    pmFPAfile *astrom = psMetadataLookupPtr(NULL, config->files, "PSWARP.ASTROM");
+    if (!astrom) {
+        astrom = input;
+    }
+
+    if (astrom->camera != input->camera) {
+        psError(PSWARP_ERR_DATA, true, "Input camera and astrometry camera do not match.");
+        return false;
+    }
+
+    // select the output readout
+    pmFPAview *view = pmFPAviewAlloc(0);
+    view->chip = 0;
+    view->cell = 0;
+    view->readout = 0;
+    pmReadout *output = pmFPAfileThisReadout(config->files, view, "PSWARP.OUTPUT.BKGMODEL");
+    if (!output) {
+        psError(PSWARP_ERR_CONFIG, true, "Can't find output background data!\n");
+        return false;
+    }
+    psFree (view);
+
+    // 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.
+    {
+        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);
+    }
+
+    // Read the input astrometry
+    // XXX rather than use the activations here, this should just explicitly loop over the desired filerule
+    {
+        pmFPAfileActivate(config->files, true, "PSWARP.ASTROM");
+
+        pmChip *chip;
+        pmFPAview *view = pmFPAviewAlloc(0);
+        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, "AChip %d: %x %x\n", view->chip, chip->file_exists, chip->process);
+            if (!chip->process || !chip->file_exists) { continue; }
+            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, "ACell %d: %x %x\n", view->cell, cell->file_exists, cell->process);
+                if (!cell->process || !cell->file_exists) { continue; }
+                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");
+    }
+
+/*     // Turn on the source output --- we need to get rid of these so that we can measure the PSF */
+/*     pmFPAfileActivate(config->files, true, "PSWARP.OUTPUT.SOURCES"); */
+
+    // Don't care about the skycell anymore --- we've read it, and that's all we need to do.
+    pmFPAfileActivate(config->files, false, "PSWARP.SKYCELL");
+    view = pmFPAviewAlloc(0);
+
+    // find the FPA phu
+    bool bilevelAstrometry = false;
+    pmHDU *phu = pmFPAviewThisPHU(view, astrom->fpa);
+    if (phu) {
+        char *ctype = psMetadataLookupStr(NULL, phu->header, "CTYPE1");
+        if (ctype) {
+            bilevelAstrometry = !strcmp (&ctype[4], "-DIS");
+        }
+    }
+    if (bilevelAstrometry) {
+        if (!pmAstromReadBilevelMosaic(input->fpa, phu->header)) {
+            psError(psErrorCodeLast(), false, "Unable to read bilevel mosaic astrometry for input FPA.");
+            psFree(view);
+            psFree(stats);
+            goto DONE;
+        }
+    }
+
+
+
+    psList *cells = psListAlloc(NULL);  // List of cells, for concepts averaging
+
+    // files associated with the science image
+    if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) {
+        psError(psErrorCodeLast(), false, "Unable to read files.");
+        goto DONE;
+    }
+
+    pmChip *chip;
+    while ((chip = pmFPAviewNextChip (view, input->fpa, 1)) != NULL) {
+        psTrace ("pswarp", 4, "DChip %d: %x %x\n", view->chip, chip->file_exists, chip->process);
+        if (!chip->process || !chip->file_exists) { continue; }
+        if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) {
+            psError(psErrorCodeLast(), false, "Unable to read files.");
+            goto DONE;
+        }
+
+        // read WCS data from the corresponding header
+        pmHDU *hdu = pmFPAviewThisHDU (view, astrom->fpa);
+        psMetadataAddS32(config->arguments,PS_LIST_TAIL,
+                         "INTERPOLATION.MODE",PS_META_REPLACE,"", PS_INTERPOLATE_BILINEAR); ///< Mode
+        psMetadataAddS32(config->arguments,PS_LIST_TAIL,
+                         "INTERPOLATION.NUM",PS_META_REPLACE,"", 64); ///< Mode
+        
+        if (bilevelAstrometry) {
+            if (!pmAstromReadBilevelChip (chip, hdu->header)) {
+                psError(psErrorCodeLast(), false, "Unable to read bilevel chip astrometry for input FPA.");
+                psFree(view);
+                psFree(stats);
+                goto DONE;
+            }
+        } else {
+            // we use a default FPA pixel scale of 1.0
+            if (!pmAstromReadWCS (input->fpa, chip, hdu->header, 1.0)) {
+                psError(psErrorCodeLast(), false, "Unable to read WCS astrometry for input FPA.");
+                psFree(view);
+                psFree(stats);
+                goto DONE;
+            }
+        }
+#define CORRECT_INPUT_WCS 1
+#if CORRECT_INPUT_WCS
+        // Correct the input WCS
+        pmAstromWCS *WCS = pmAstromWCSfromHeader(hdu->header);
+
+        double cd1f = 1.0 * 400;
+        double cd2f = 1.0 * 400;
+
+        WCS->crpix1 = WCS->crpix1 / cd1f + 0.5;
+        WCS->crpix2 = WCS->crpix2 / cd2f - 0.5;
+
+        WCS->cdelt1 *= cd1f;
+        WCS->cdelt2 *= cd2f;
+
+        WCS->trans->x->coeff[1][0] *= cd1f;
+        WCS->trans->x->coeff[0][1] *= cd2f;
+        WCS->trans->y->coeff[1][0] *= cd1f;
+        WCS->trans->y->coeff[0][1] *= cd2f;
+
+
+        pmAstromWCStoHeader (hdu->header,WCS);
+        // End WCS work.
+#endif
+        
+        pmCell *cell;
+        while ((cell = pmFPAviewNextCell (view, input->fpa, 1)) != NULL) {
+            psTrace ("pswarp", 4, "DCell %d: %x %x\n", view->cell, cell->file_exists, cell->process);
+            if (!cell->process || !cell->file_exists) { continue; }
+            if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) {
+                psError(psErrorCodeLast(), false, "Unable to read files.");
+                goto DONE;
+            }
+
+            psListAdd(cells, PS_LIST_TAIL, cell);
+
+            // process each of the readouts
+            pmReadout *readout;
+            while ((readout = pmFPAviewNextReadout(view, input->fpa, 1)) != NULL) {
+                if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
+                    psError(psErrorCodeLast(), false, "Unable to read files.");
+                    goto DONE;
+                }
+                if (!readout->data_exists) {
+                    continue;
+                }
+
+                // 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
+                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);
+                }
+
+                for (int x = 0; x < readout->image->numCols; x++) {
+                  for (int y = 0; y < readout->image->numRows; y++) {
+                    readout->image->data.F32[y][x] = readout->image->data.F32[y][x] * (cd1f * cd2f) /
+                      (psMetadataLookupS32(&mdok,config->arguments,"BKG.XGRID") *
+                       psMetadataLookupS32(&mdok,config->arguments,"BKG.YGRID"));
+                  }
+                }
+                psMetadataAddS32(config->arguments,PS_LIST_TAIL, "INTERPOLATION.MODE", PS_META_REPLACE, "", 2);
+/*              psMetadataAddS32(config->arguments,PS_LIST_TAIL, "GRID.NX",PS_META_REPLACE,"",100); ///< Number of grid points in x */
+/*              psMetadataAddS32(config->arguments,PS_LIST_TAIL, "GRID.NY",PS_META_REPLACE,"",100); ///< Number of grid points in y */
+
+                pswarpTransformReadout(output, readout, config);
+                
+                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;
+        }
+    }
+
+    if (!output->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;
+    }
+    
+    pmCell *outCell = output->parent;   ///< Output cell
+    pmChip *outChip = outCell->parent;  ///< Output chip
+    pmFPA *outFPA = outChip->parent;    ///< Output FP
+
+    if (!pswarpPixelsLit(output, stats, config)) {
+        psError(psErrorCodeLast(), false, "Unable to calculate pixel regions.");
+        psFree(cells);
+        psFree(view);
+        goto DONE;
+    }
+/*     bool doStats = psMetadataLookupBool(&mdok,recipe,"MASK.STATS"); */
+/*     if (doStats) { */
+/*       if (!pswarpMaskStats(output, 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); */
+/*         } */
+/*     } */
+
+/*     if (!pmConceptsAverageCells(outCell, cells, NULL, NULL, false)) { */
+/*         psError(psErrorCodeLast(), false, "Unable to average cell concepts."); */
+/*         psFree(stats); */
+/*         psFree(cells); */
+/*         psFree(view); */
+/*         goto DONE; */
+/*     } */
+/*     psFree(cells); */
+
+    psRegion *trimsec = psMetadataLookupPtr(NULL, outCell->concepts, "CELL.TRIMSEC"); ///< Trim section
+    trimsec->x0 = trimsec->x1 = trimsec->y0 = trimsec->y1 = 0; ///< All pixels
+
+    if (!psMetadataCopy(outFPA->concepts, input->fpa->concepts)) {
+        psError(psErrorCodeLast(), false, "Unable to copy FPA concepts from input to output.");
+        psFree(stats);
+        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"); */
+/*     } */
+
+    pmHDU *hdu = outFPA->hdu;           ///< HDU for the output warped image
+
+    // Copy header from target
+    {
+        pmFPAview *skyView = pmFPAviewAlloc(0); ///< View into skycell
+        skyView->chip = skyView->cell = 0;
+        pmCell *cell = pmFPAfileThisCell(config->files, skyView, "PSWARP.SKYCELL"); // Skycell cell
+        psFree(skyView);
+        pmHDU *skyHDU = pmHDUFromCell(cell); ///< HDU
+        if (!skyHDU) {
+            psError(PSWARP_ERR_DATA, false, "Unable to find skycell HDU.");
+            psFree(view);
+            goto DONE;
+        }
+        hdu->header = psMetadataCopy(hdu->header, skyHDU->header);
+    }
+
+    pswarpVersionHeader(hdu->header);
+    
+    if (!pmAstromWriteWCS(hdu->header, outFPA, outChip, WCS_NONLIN_TOL)) {
+        psError(psErrorCodeLast(), false, "Unable to generate WCS header.");
+        psFree(stats);
+        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
+    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, 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); */
+/*     } */
+
+/*     // Perform statistics on the output image */
+/*     if (stats) { */
+/*         if (!ppStatsFPA(stats, output->parent->parent->parent, view, maskValue, config)) { */
+/*             psWarning("Unable to perform statistics on warped image."); */
+/*         } */
+/*     } */
+    
+
+    // 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);
+    psFree(view);
+
+ DONE:
+
+    return true;
+}