Changeset 36580 for trunk/ppBackground/src/ppBackgroundStackCamera.c

Timestamp:

Mar 7, 2014, 6:08:39 PM (12 years ago)

Author:

watersc1

Message:

ppBackgroundStack compiles. I know it's not useful yet, because I haven't implemented the 'apply' and 'save the outputs' parts.

File:

• trunk/ppBackground/src/ppBackgroundStackCamera.c (modified) (9 diffs)

trunk/ppBackground/src/ppBackgroundStackCamera.c

@@ -31 +31 @@
 {
     bool status;                        // Status of file definition
-
+    pmConfig *config = data->config;    // Because I'm reusing code.
+    int u,v;
+    // FIX Figure out what stacks we have to deal with.
+    psString stackName = data->stacks->data[0];
+    fileArguments("IMAGE", stackName, "Input image", data->config);
+    pmFPAfile *stack = pmFPAfileDefineFromArgs(&status, data->config, "PPBACKGROUND.STACK",
+                                               "IMAGE");
+    if (!status || !stack) {
+      psError(psErrorCodeLast(), false, "Failed to build file from PPBACKGROUND.STACK");
+      return false;
+    }
+
+    // Read over the input background models.    
     psMetadataIterator *iter = psMetadataIteratorAlloc(data->contents, PS_LIST_HEAD, NULL); // Iterator
     psMetadataItem *item; // Item from iteration
-    int i = 0;
+    //    int i = 0;
     while ((item = psMetadataGetAndIncrement(iter))) {
       if (item->type != PS_DATA_METADATA) {
-        psError(PPBACKGROUND_ERR_ARGUMENTS, true
+        psError(PPBACKGROUND_ERR_ARGUMENTS, true,
                 "Component %s of the input metadata is not of type METADATA", item->name);
         psFree(iter);
@@ -47 +59 @@
       psString smfFileName = psMetadataLookupStr(NULL, input, "astrom");
 
-      // FIX Figure out how to calculate the smf name
-      // psString cam_file = pmFPAfileNameFromRule("PSASTRO.OUTPUT");
-      // pmFPAfile *smfFile = pmFPAfileBindFromArgs(&status, NULL, config);
-      // pmFPAfile *smfFile = pmFPAfileDefineInput(config,
+      // Read the smf file from this item
       fileArguments("astrom",smfFileName,"",config);
       pmFPAfile *smfFile = pmFPAfileDefineFromArgs(&status, config, "PSWARP.ASTROM","astrom");
 
       // Read the SMF data
-
       pmFPAview *view = pmFPAviewAlloc(0); // Pointer into FPA hierarchy
       if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
@@ -95 +103 @@
         }
         // read WCS data from the corresponding header
-        pmHDU *hdu = pmFPAviewThisHDU (view, astromFile->fpa);
+        pmHDU *hdu = pmFPAviewThisHDU (view, smfFile->fpa);
         if (bilevelAstrometry) {
           if (!pmAstromReadBilevelChip (chip, hdu->header)) {
@@ -108 +116 @@
           // we use a default FPA pixel scale of 1.0
           psWarning("Reading WCS astrometry for chip %s.", chipName);
-          if (!pmAstromReadWCS(astromFile->fpa, chip, hdu->header, 1.0)) {
+          if (!pmAstromReadWCS(smfFile->fpa, chip, hdu->header, 1.0)) {
             psError(psErrorCodeLast(), false, "Unable to read WCS astrometry for input FPA.");
             psFree(view);
@@ -116 +124 @@
 
         // FIX Load model data for this chip
-        psString modelFileName = psMetadataLookupStr(NULL, input->models, chipName);
+        psString modelFileName = psMetadataLookupStr(NULL, input, chipName);
         psFits *modelFits = psFitsOpen(modelFileName,"r");
-        psImage *image = psFitsReadImage(modelFits,NULL,0);
+        int nC, nR;
+        psElemType *type;
+        psRegion region;
+        psFitsImageSize(&nC,&nR,type,modelFits,region);
+        psImage *image = psFitsReadImage(modelFits,region,0);
         psMetadata *header = psFitsReadHeader(NULL,modelFits);
 
         // Allocate the data structures for this chip
-        psImage *raim = psImageAlloc(image->cols,image->rows,PS_TYPE_F32);
-        psImage *decim = psImageAlloc(image->cols,image->rows,PS_TYPE_F32);
-        psImage *model = psImageAlloc(image->cols,image->rows,PS_TYPE_F32);
+        psImage *raim = psImageAlloc(image->numCols,image->numRows,PS_TYPE_F32);
+        psImage *decim = psImageAlloc(image->numCols,image->numRows,PS_TYPE_F32);
+        psImage *model = psImageAlloc(image->numCols,image->numRows,PS_TYPE_F32);
 
         // Read header and construct original positions
@@ -142 +154 @@
         psSphere *sky = psSphereAlloc(); // Sky coordinates
         
-        for (v = 0; v < image->nrows; v++) {
+        for (v = 0; v < image->numRows; v++) {
           pix->y = (v - yoffset) * ybin;
-          for (u = 0; u < image->ncols; u++) {
+          for (u = 0; u < image->numCols; u++) {
             pix->x = (u - xoffset) * xbin;
             psPlaneTransformApply(fp, chip->toFPA, pix);
             psPlaneTransformApply(tp, smfFile->fpa->toTPA, fp);
-            psDeproject(sky, tp, astromFile->fpa->toSky);
-
-            // FIX project sky to stack projection cell grid.
-            if (!data->radians) {
-              sky->r *= 180.0 / M_PI;
-              sky->d *= 180.0 / M_PI;
-            }
-            raim->data.F32[v][u] = sky->r;
-            decim->data.F32[v][u] = sky->d;
+            psDeproject(sky, tp, smfFile->fpa->toSky);
+
+            psProject(tp,sky,stack->fpa->toSky);
+
+            raim->data.F32[v][u] = tp->x;
+            decim->data.F32[v][u] = tp->y;
             model->data.F32[v][u] = 0.0;
 
             // Check the bounds so we'll know how large of an area to model in the map
-            if (sky->r < data->ra_min) { data->ra_min = sky->r; }
-            else if (sky->r > data->ra_max) { data->ra_max = sky->r; }
-            if (sky->d < data->dec_min) { data->dec_min = sky->d; }
-            else if (sky->d > data->dec_max) { data->dec_max = sky->d; }
+            if (tp->x < data->ra_min) { data->ra_min = tp->x; }
+            else if (tp->x > data->ra_max) { data->ra_max = tp->x; }
+            if (tp->y < data->dec_min) { data->dec_min = tp->y; }
+            else if (tp->y > data->dec_max) { data->dec_max = tp->y; }
           }
         }
@@ -187 +196 @@
           // FIX this should try to find the imagefile. 
           if (data->fit_OTAS) { // We are fitting OTAs, so allocate a new image
-            psImage *solution = psImageAlloc(image->cols,image->rows,PS_TYPE_F32);
+            psImage *solution = psImageAlloc(image->numCols,image->numRows,PS_TYPE_F32);
             psMetadataAddPtr(data->OTA_solutions,PS_LIST_TAIL, chipName, PS_DATA_UNKNOWN | PS_META_REPLACE, "OTA solution element", solution);
           }
@@ -194 +203 @@
             psStringAppend(&solutionFileName, ".%s.fits",chipName);
             psFits *solutionFits = psFitsOpen(solutionFileName,"r");
-            psImage *solution = psFitsReadImage(solutionFits,NULL,0);
+            psImage *solution = psFitsReadImage(solutionFits,region,0);
+            psMetadataAddPtr(data->OTA_solutions,PS_LIST_TAIL, chipName, PS_DATA_UNKNOWN | PS_META_REPLACE, "OTA solution element", solution);
           }
         }
@@ -215 +225 @@
     }
 
+    // Allocate the modelMap for the region we're covering.
+
+    psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
+    psImageBinning *binning = psImageBinningAlloc();
+    binning->nXruff = 100; // Number of samples
+    binning->nYruff = 100; 
+    binning->nXfine = ceil(data->ra_max - data->ra_min) + 1; // This is the range we're looking at
+    binning->nYfine = ceil(data->dec_max - data->dec_min) + 1;
+
+    data->modelMap = psImageMapNoImageAlloc( binning,stats);
+    
     return true;
 }