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

Timestamp:

Apr 2, 2014, 1:13:06 PM (12 years ago)

Author:

watersc1

Message:

This version appears to complete and generate the expected output products. Now to test and tune the parameters.

File:

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

trunk/ppBackground/src/ppBackgroundStackLoop.c

@@ -68 +68 @@
         
       }
-      
-      
+    } // End initialization block.      
+
+
+    // Loop over the input images, and apply the models to construct the restored versions.
+    for (i = 0; i < data->stack_data->n; i++) {
+      pmFPAfile *stack = data->stack_data->data[i];
+      pmFPAview *view = pmFPAviewAlloc(0);
+
+      // PART 1:
+      // Determine the extent of the model map for this stack
+      data->x_min = 99e99; data->x_max = -99e99;
+      data->y_min = 99e99; data->y_max = -99e99;
+      // Allocate the modelMap for the region we're covering.
+      psPlane *pix = psPlaneAlloc();   // Pixel coordinates on chip
+      psPlane *tp = psPlaneAlloc();    // Focal plane coordinates
+      
+      pix->x = 0; pix->y = 0;
+      psPlaneTransformApply(tp, stack->fpa->toTPA, pix);
+      printf("%f %f -> %f %f\n",pix->x,pix->y,tp->x,tp->y);
+      if (tp->x < data->x_min) { data->x_min = tp->x; }
+      if (tp->x > data->x_max) { data->x_max = tp->x; }
+      if (tp->y < data->y_min) { data->y_min = tp->y; }
+      if (tp->y > data->y_max) { data->y_max = tp->y; }
+
+      pix->x = 6240; pix->y = 0;
+      psPlaneTransformApply(tp, stack->fpa->toTPA, pix);
+      printf("%f %f -> %f %f\n",pix->x,pix->y,tp->x,tp->y);
+      if (tp->x < data->x_min) { data->x_min = tp->x; }
+      if (tp->x > data->x_max) { data->x_max = tp->x; }
+      if (tp->y < data->y_min) { data->y_min = tp->y; }
+      if (tp->y > data->y_max) { data->y_max = tp->y; }
+
+      pix->x = 6240; pix->y = 6243;
+      psPlaneTransformApply(tp, stack->fpa->toTPA, pix);
+      printf("%f %f -> %f %f\n",pix->x,pix->y,tp->x,tp->y);
+      if (tp->x < data->x_min) { data->x_min = tp->x; }
+      if (tp->x > data->x_max) { data->x_max = tp->x; }
+      if (tp->y < data->y_min) { data->y_min = tp->y; }
+      if (tp->y > data->y_max) { data->y_max = tp->y; }
+
+      pix->x = 0; pix->y = 6243;
+      psPlaneTransformApply(tp, stack->fpa->toTPA, pix);
+      printf("%f %f -> %f %f\n",pix->x,pix->y,tp->x,tp->y);
+      if (tp->x < data->x_min) { data->x_min = tp->x; }
+      if (tp->x > data->x_max) { data->x_max = tp->x; }
+      if (tp->y < data->y_min) { data->y_min = tp->y; }
+      if (tp->y > data->y_max) { data->y_max = tp->y; }
+
+/*       data->x_min -= data->ra_min; */
+/*       data->x_max -= data->ra_min; */
+/*       data->y_min -= data->dec_min; */
+/*       data->y_max -= data->dec_min; */
+      
+      
+      psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
+      psImageBinning *binning = psImageBinningAlloc();
+      binning->nXruff = 15; // Number of samples
+      binning->nYruff = 15; 
+      //    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;
+      binning->nXfine = ceil(data->x_max - data->x_min) + 1;
+      binning->nYfine = ceil(data->y_max - data->y_min) + 1;
+      binning->nXskip = floor(data->x_min - data->ra_min) + 1;
+      binning->nYskip = floor(data->y_min - data->dec_min) + 1;
+      psImageBinningSetScale(binning,PS_IMAGE_BINNING_CENTER);
+      printf("Sizes: sky: %f %f -> %f %f :: tp: %f %f -> %f %f :: map: %d %d\n",
+             0.0,0.0,0.0,0.0,data->ra_min,data->dec_min,data->ra_max,data->dec_max,binning->nXfine,binning->nYfine);
+      printf("Sizes: corners: %f %f -> %f %f map: %d %d\n",
+             data->x_min,data->y_min,data->x_max,data->y_max,
+             binning->nXfine,binning->nYfine);
+      psImage *sizeImage = psImageAlloc(binning->nXfine,binning->nYfine,PS_TYPE_F32);
+      P_PSIMAGE_SET_COL0(sizeImage, data->x_min);
+      P_PSIMAGE_SET_ROW0(sizeImage, data->y_min);
+      data->modelMap = psImageMapAlloc(sizeImage,binning,stats);
+/*       psFree(sizeImage); */
+/*       data->modelMap = psImageMapNoImageAlloc( binning,stats); */
+/*       P_PSIMAGE_SET_COL0(data->modelMap->map, (data->x_min - binning->nXskip) / binning->nXbin); */
+/*       P_PSIMAGE_SET_ROW0(data->modelMap->map, (data->y_min - binning->nYskip) / binning->nYbin); */
+
+      // force col0/row0
+/*       data->modelMap->map->col0 = data->modelMap->binning->nXskip; */
+/*       data->modelMap->map->row0 = data->modelMap->binning->nYskip; */
+      
+      // PART 2:
+      // Solve the data into a model for this region of the sky.
+
       // 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");
@@ -83 +167 @@
         return(false);
       }
-    } // End initialization block.
-
-    // This is where a loop would likely start.
-    {
-      // Construct the offset information
-      printf("Model fit!\n");
-      if (!ppBackgroundStackDataModelFit(data)) {
-        psError(psErrorCodeLast(), false, "Error determining the exposure/OTA scaling.");
-        return(false);
-      }
-      
-      // Apply full correction
-      printf("Calib apply!\n");
-      if (!ppBackgroundStackCalibApply(data)) {
-        psError(psErrorCodeLast(), false, "Error applying the calibration models.");
-        return(false);
-      }      
-      
-      // Fit the new model
-      printf("Determining model!\n");
-      if (!ppBackgroundStackModelFit(data)) {
-        psError(psErrorCodeLast(), false, "Error determining the modelMap object.");
-        return(false);
-      }
-    } // End loop
-
-    // Loop over the input images, and apply the models to construct the restored versions.
-
-    for (i = 0; i < data->stack_data->n; i++) {
-      pmFPAfile *stack = data->stack_data->data[i];
-      pmFPAview *view = pmFPAviewAlloc(0);
-
+
+      // This is where an iterative solution loop would likely start.
+      {
+        // Construct the offset information
+        printf("Model fit!\n");
+        if (!ppBackgroundStackDataModelFit(data)) {
+          psError(psErrorCodeLast(), false, "Error determining the exposure/OTA scaling.");
+          return(false);
+        }
+        
+        // Apply full correction
+        printf("Calib apply!\n");
+        if (!ppBackgroundStackCalibApply(data)) {
+          psError(psErrorCodeLast(), false, "Error applying the calibration models.");
+          return(false);
+        }      
+        
+        // Fit the new model
+        printf("Determining model!\n");
+        if (!ppBackgroundStackModelFit(data)) {
+          psError(psErrorCodeLast(), false, "Error determining the modelMap object.");
+          return(false);
+        }
+      } // End loop
+
+      // PART 3:
+      // Define output products
       // 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");
       
@@ -134 +213 @@
       stack_model->save = true;
       stack_corr->save = true;
-
+      
       printf("I'm about to loop over the parts of this stack: %d\n",i);
       // Iterate over the images.
@@ -144 +223 @@
         return(false);
       }
-
       
       pmChip *chip;
@@ -151 +229 @@
           continue;
         }
-
+        
         if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
           psError(psErrorCodeLast(), false, "load failure for Chip");
@@ -256 +334 @@
       }
       psFree(view);
-    }      
+      psFree(data->modelMap);
+    }