Changeset 36615 for trunk/ppBackground/src/ppBackgroundStackLoop.c

Timestamp:

Mar 27, 2014, 4:21:00 PM (12 years ago)

Author:

watersc1

Message:

Working with metadata structures is far easier than pmFPAfiles. First cleanly running-to-completion version. Now to sort out the resolution issues.

File:

• trunk/ppBackground/src/ppBackgroundStackLoop.c (modified) (14 diffs)

trunk/ppBackground/src/ppBackgroundStackLoop.c

@@ -20 +20 @@
       // If we didn't load the OTA solution from an external source, we need to build one.
       if (data->fit_OTAS) {
+        printf("Fitting OTAs!\n");
         if (!ppBackgroundStackModelFitOTASolution(data)) {
           // Currently can't fail.
@@ -25 +26 @@
           return(false);
         }
-      }
+
+        if (data->OTApath) {
+          // This will write the solutions out.
+          psMetadataIterator *iter = psMetadataIteratorAlloc(data->OTA_solutions, PS_LIST_HEAD, NULL); // Iterator
+          psMetadataItem *item; // Item from iteration
+          int i = 0;
+          while ((item = psMetadataGetAndIncrement(iter))) {
+            i++;
+            if (item->type != PS_DATA_IMAGE) {
+              psString chipName = item->name;
+              psImage  *image    = item->data.V;
+              psMetadata *header = psMetadataAlloc();
+              psString solutionFileName = psStringCopy(data->OTApath);
+              psStringAppend(&solutionFileName, ".%s.fits",chipName);
+              psString resolved = pmConfigConvertFilename(solutionFileName,config,true,true);
+              psFits *solutionFits = psFitsOpen(resolved,"w");
+              if (!solutionFits) {
+                psError(2, false, "Unable to open FITS file %s to write model %s.", resolved, chipName);
+                psFitsClose(solutionFits);
+                psFree(resolved);
+                return(false);
+              }
+              if (!psFitsWriteImage(solutionFits, header, image, 0, NULL)) {
+                psError(2, false, "Unable to write FITS image %s.", resolved);
+                psFitsClose(solutionFits);
+                psFree(resolved);
+                return false;
+              }
+              if (!psFitsClose(solutionFits)) {
+                psError(2, false, "Unable to close FITS image %s.", resolved);
+                psFree(resolved);
+                return false;
+              }
+              psFree(resolved);
+              psFree(solutionFileName);
+            }
+          }
+          psFree(item);
+          psFree(iter);
+        }
+        
+      }
+      
       
       // This seems wrong, but I need a blank modelMap object, so we fit the zero-data we've stored in the calib objects
+      printf("Determining blank modelMap!\n");
       if (!ppBackgroundStackModelFit(data)) {
         psError(psErrorCodeLast(), false, "Error determining the blank modelMap object.");
@@ -34 +78 @@
       
       // Apply OTA solution
+      printf("Calib apply!\n");
       if (!ppBackgroundStackCalibApply(data)) {
         psError(psErrorCodeLast(), false, "Error applying the calibration models.");
@@ -43 +88 @@
     {
       // Construct the offset information
+      printf("Model fit!\n");
       if (!ppBackgroundStackDataModelFit(data)) {
         psError(psErrorCodeLast(), false, "Error determining the exposure/OTA scaling.");
@@ -49 +95 @@
       
       // Apply full correction
+      printf("Calib apply!\n");
       if (!ppBackgroundStackCalibApply(data)) {
         psError(psErrorCodeLast(), false, "Error applying the calibration models.");
@@ -55 +102 @@
       
       // Fit the new model
+      printf("Determining model!\n");
       if (!ppBackgroundStackModelFit(data)) {
         psError(psErrorCodeLast(), false, "Error determining the modelMap object.");
@@ -66 +114 @@
       pmFPAfile *stack = data->stack_data->data[i];
       pmFPAview *view = pmFPAviewAlloc(0);
-      // Define output image
-      pmFPAfile *stack_model = pmFPAfileDefineOutput(config,stack->fpa,"PPBACKGROUND.STACK.MODEL");
-
+
+      // Define output image.  Why is this always so hard to do?
+      pmFPA *tmp_fpa1,*tmp_fpa2;
+      tmp_fpa1 = pmFPAConstruct(config->camera,config->cameraName);
+      tmp_fpa2 = pmFPAConstruct(config->camera,config->cameraName);
+
+      pmFPAfile *stack_model = pmFPAfileDefineOutput(config,tmp_fpa1,"PPBACKGROUND.STACK.MODEL");
+      
       if (!stack_model) {
         psError(psErrorCodeLast(), false, "Unable to generate output model");
@@ -74 +127 @@
       }
       
-      pmFPAfile *stack_corr  = pmFPAfileDefineOutput(config,stack->fpa,"PPBACKGROUND.STACK.OUTPUT");
+      pmFPAfile *stack_corr  = pmFPAfileDefineOutput(config,tmp_fpa2,"PPBACKGROUND.STACK.OUTPUT");
       if (!stack_corr) {
         psError(psErrorCodeLast(), false, "Unable to generate output result");
@@ -82 +135 @@
       stack_corr->save = true;
 
+      printf("I'm about to loop over the parts of this stack: %d\n",i);
       // Iterate over the images.
+      pmFPAfileActivate(config->files,true,"PPBACKGROUND.STACK");
+      pmFPAfileActivate(config->files,true,"PPBACKGROUND.STACK.MODEL");
+      pmFPAfileActivate(config->files,true,"PPBACKGROUND.STACK.OUTPUT");
+      if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
+        psError(psErrorCodeLast(), false, "load failure for Chip");
+        return(false);
+      }
+
+      
       pmChip *chip;
       while ((chip = pmFPAviewNextChip(view, stack->fpa, 1))) {
@@ -93 +156 @@
           return(false);
         }
-
+        printf("  I'm in a chip\n");
         pmCell *cell;
         
@@ -109 +172 @@
             continue;
           }
+          printf("    I'm in a cell\n");
+
 
           // process each of the readouts
@@ -120 +185 @@
               continue;
             }
-
-            pmReadout *model = pmFPAfileThisReadout(config->files,view,"PPBACKGROUND.STACK.MODEL");
-            pmReadout *corr  = pmFPAfileThisReadout(config->files,view,"PPBACKGROUND.STACK.OUTPUT");
-
+            printf("      I'm in a readout\n");
+
+            // Futz with things to get an acceptable output product.
+            pmCell *model_cell = pmFPAviewThisCell(view,stack_model->fpa);
+            pmCell *corr_cell  = pmFPAviewThisCell(view,stack_corr->fpa);
+            
+            pmReadout *model = pmReadoutAlloc(model_cell);
+            pmReadout *corr  = pmReadoutAlloc(corr_cell);
+            model->data_exists = true;
+            corr->data_exists = true;
+            model->parent->data_exists = true;
+            corr->parent->data_exists = true;
+            model->parent->parent->data_exists = true;
+            corr->parent->parent->data_exists = true;
+            model->image = psImageAlloc(readout->image->numCols,readout->image->numRows,PS_TYPE_F32);
+            corr->image = psImageAlloc(readout->image->numCols,readout->image->numRows,PS_TYPE_F32);
+
+            model_cell->concepts = psMemIncrRefCounter(cell->concepts);
+            model_cell->conceptsRead = cell->conceptsRead;
+            corr_cell->concepts = psMemIncrRefCounter(cell->concepts);
+            corr_cell->conceptsRead = cell->conceptsRead;
+            
             psPlane *pix = psPlaneAlloc();   // Pixel coordinates on chip
             psPlane *fp = psPlaneAlloc();    // Focal plane coordinates
@@ -130 +213 @@
             int x,y;
             for (y = 0; y < readout->image->numRows; y++) {
+              pix->y = y;
               for (x = 0; x < readout->image->numCols; x++) {
+                pix->x = x;
                 // Calculate model for each pixel of output
-                psPlaneTransformApply(fp, chip->toFPA, pix);
-                psPlaneTransformApply(tp, stack->fpa->toTPA, fp);
+                //              psPlaneTransformApply(fp, chip->toFPA, pix);
+                psPlaneTransformApply(tp, stack->fpa->toTPA, pix);
                 
                 model->image->data.F32[y][x] = psImageMapEval(data->modelMap,tp->x,tp->y);
                 corr->image->data.F32[y][x] = readout->image->data.F32[y][x] + model->image->data.F32[y][x];
+                
               }
             }
@@ -143 +229 @@
             psFree(tp);
           } // Close readout
-
+          printf("    I'm done with that readout\n");
           // Close output image
-          // Close cells (XXX shouldn't pmFPAfileClose iterate down as needed?)
-          view->cell = -1;
-          while ((cell = pmFPAviewNextCell(view, stack->fpa, 1)) != NULL) {
-            if (!cell->process || !cell->file_exists) {
-              continue;
-            }
-            if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
-                psError(psErrorCodeLast(), false, "save failure for Cell");
-                return(false);
-            }
-          }
         } // Close Cell
-        
+        // Close cells (XXX shouldn't pmFPAfileClose iterate down as needed?)
+        view->cell = -1;
+        while ((cell = pmFPAviewNextCell(view, stack->fpa, 1)) != NULL) {
+          if (!cell->process || !cell->file_exists) {
+            continue;
+          }
+          if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
+            psError(psErrorCodeLast(), false, "save failure for Cell");
+            return(false);
+          }
+        }
+
         // Close chip
         if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {