Changeset 34800 for trunk/pswarp/src/pswarpLoop.c

Timestamp:

Dec 11, 2012, 2:04:31 PM (13 years ago)

Author:

watersc1

Message:

Merge from my branch including background restoration and stack stage projection cell binned images.

Location:

trunk

Files:

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

trunk/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)) {
@@ -195 +197 @@
             }
         }
-
+        
         pmCell *cell;
         while ((cell = pmFPAviewNextCell (view, input->fpa, 1)) != NULL) {
@@ -227 +229 @@
 
                 pswarpTransformReadout(output, readout, config);
-
+                
                 if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
                     psError(psErrorCodeLast(), false, "Unable to write files.");
@@ -255 +257 @@
         goto DONE;
     }
-
+    
     pmCell *outCell = output->parent;   ///< Output cell
     pmChip *outChip = outCell->parent;  ///< Output chip
@@ -450 +452 @@
     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) {
+#if 0
+          // This needs to be removed because otherwise it throws an error of duplicate PSPHOT.DETECTIONS.
+            if (!chip->process || !chip->file_exists) { continue; }
+            if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
+                psError(psErrorCodeLast(), false, "Unable to read files.");
+                goto DONE;
+            }
+#endif
+            pmCell *cell;
+            while ((cell = pmFPAviewNextCell (view, input->fpa, 1)) != NULL) {
+              psTrace ("pswarp", 4, "ACell %d: %x %x %d\n", view->cell, cell->file_exists, cell->process,psErrorCodeLast());
+                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");
+    }
+
+    // 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);
+
+    //    pmAstromWCS *WCSF = pmAstromWCSfromHeader(phu->header);
+    
+    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;
+        }
+
+        // Pull information from the header of the background files so we can use it to set values.
+        pmHDU *hdu = pmFPAviewThisHDU(view,input->fpa);
+        psMetadata *header = hdu->header;
+        
+        int IMAXIS1 = psMetadataLookupS32(NULL,header,"IMNAXIS1");
+        int IMAXIS2 = psMetadataLookupS32(NULL,header,"IMNAXIS2");
+        int NAXIS1 = psMetadataLookupS32(NULL,header,"NAXIS1");
+        int NAXIS2 = psMetadataLookupS32(NULL,header,"NAXIS2");
+        char *CCDSUM = psMetadataLookupStr(NULL,header,"CCDSUM");
+        int CCDSUM1 = atoi(strtok(CCDSUM," "));
+        int CCDSUM2 = atoi(strtok(NULL," "));
+        
+        psMetadataAddF32(config->arguments,PS_LIST_TAIL,"BKG_WARP_XOFFSET", PS_META_REPLACE,
+                         "xoffset for background model data", (NAXIS1 * CCDSUM1 - IMAXIS1) / (2.0 * CCDSUM1));
+        psMetadataAddF32(config->arguments,PS_LIST_TAIL,"BKG_WARP_YOFFSET", PS_META_REPLACE,
+                         "yoffset for background model data", (NAXIS2 * CCDSUM2 - IMAXIS2) / (2.0 * CCDSUM2));
+        psTrace("pswarp",5,"%d %d %d %d %d %d %g %g %d %d",
+                psMetadataLookupS32(NULL,header,"IMNAXIS1"),
+                psMetadataLookupS32(NULL,header,"IMNAXIS2"),
+                psMetadataLookupS32(NULL,header,"NAXIS1"),
+                psMetadataLookupS32(NULL,header,"NAXIS2"),
+                CCDSUM1,CCDSUM2,
+                psMetadataLookupF32(NULL,config->arguments,"BKG_WARP_XOFFSET"),
+                psMetadataLookupF32(NULL,config->arguments,"BKG_WARP_YOFFSET"),
+                psMetadataLookupS32(NULL,config->arguments,"BKG.XGRID"),
+                psMetadataLookupS32(NULL,config->arguments,"BKG.YGRID"));
+        
+        
+        // read WCS data from the corresponding header
+        hdu = pmFPAviewThisHDU (view, astrom->fpa);
+
+        pmAstromWCS *WCS = pmAstromWCSfromHeader(hdu->header);
+
+        double cd1f = (1.0 * CCDSUM1);// * (1.0 * psMetadataLookupS32(NULL,config->arguments,"BKG.XGRID"));
+        double cd2f = (1.0 * CCDSUM2);// * (1.0 * psMetadataLookupS32(NULL,config->arguments,"BKG.YGRID"));
+
+        WCS->cdelt1 *= cd1f;
+        WCS->cdelt2 *= cd2f;
+        WCS->crpix1 = WCS->crpix1 / cd1f;
+        WCS->crpix2 = WCS->crpix2 / cd2f;
+
+        // WCS->trans->x->nX/nY
+        for (int q = 0; q <= WCS->trans->x->nX; q++) {
+          for (int r = 0; r <= WCS->trans->x->nY; r++) {
+            WCS->trans->x->coeff[q][r] *= pow(cd1f,q) * pow(cd2f,r);
+          }
+        }
+        for (int q = 0; q <= WCS->trans->y->nX; q++) {
+          for (int r = 0; r <= WCS->trans->y->nY; r++) {
+            WCS->trans->y->coeff[q][r] *= pow(cd1f,q) * pow(cd2f,r);
+          }
+        }
+        pmAstromWCStoHeader (hdu->header,WCS);
+        
+        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 (astrom->fpa, chip, hdu->header, 1.0)) {
+                psError(psErrorCodeLast(), false, "Unable to read WCS astrometry for input FPA.");
+                psFree(view);
+                psFree(stats);
+                goto DONE;
+            }
+        }
+        // Modify structure here.
+
+        
+        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;
+                }
+
+                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, "", 8);
+
+                psMetadataAddBool(config->arguments,PS_LIST_TAIL, "BACKGROUND_WARPING", PS_META_REPLACE, "", true);
+                pswarpTransformReadout(output, readout, config);
+                psMetadataAddBool(config->arguments,PS_LIST_TAIL, "BACKGROUND_WARPING", PS_META_REPLACE, "", false);
+
+                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.");
+        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; */
+/*     } */
+    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;
+    }
+
+    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);
+
+
+    // 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;
+}