Changeset 36649 for trunk/ppBackground/src/ppBackgroundStackMath.c

Timestamp:

Apr 3, 2014, 6:11:18 PM (12 years ago)

Author:

watersc1

Message:

I'm reasonably happy with the results it's making now.

File:

• trunk/ppBackground/src/ppBackgroundStackMath.c (modified) (8 diffs)

trunk/ppBackground/src/ppBackgroundStackMath.c

@@ -55 +55 @@
           }
           psMetadata *chipData = psMetadataLookupPtr(NULL, expItem->data.md, workingChip);
+          if (!chipData) { continue; }
           psImage *image      = psMetadataLookupPtr(NULL, chipData, "bkg image");
+          if (!image) { continue; }
           psVectorAppend(tmp,image->data.F32[v][u]);
         } // End loop over exposures
@@ -63 +65 @@
 
         psFree(expIter);
+        psFree(tmp);
       } // End u
     } // End v
+    // Remove the median value from this data.  We just want the tilts, not the offsets
+    psStatsInit(stats);
+    psImageStats(stats,solution,NULL,0);
+    for (v = 0; v < solution->numRows; v++) {
+      for (u = 0; u < solution->numCols; u++) {
+        solution->data.F32[v][u] -= stats->robustMedian;
+      }
+    }
+    
+    
   } // End working chip scan
 
@@ -90 +103 @@
     psMetadataItem *chipItem;
     while ((chipItem = psMetadataGetAndIncrement(chipIter))) {
+      //      const char *chipName = chipItem->name;
+      
       psImage *image = psMetadataLookupPtr(NULL, chipItem->data.md, "bkg image");
       psImage *ra    = psMetadataLookupPtr(NULL, chipItem->data.md, "bkg ra");
       psImage *dec   = psMetadataLookupPtr(NULL, chipItem->data.md, "bkg dec");
-      
+      //      psImage *camera= psMetadataLookupPtr(NULL, data->OTA_solutions, chipName);      
       psVector *obs  = psVectorAllocEmpty(image->numRows, PS_TYPE_F32);
       psVector *model= psVectorAllocEmpty(image->numRows, PS_TYPE_F32);
-      
+
+      int j = 0;
+      int used = 0;
       for (v = 0; v < image->numRows; v++) {
         for (u = 0; u < image->numCols; u++) {
           if ((ra->data.F32[v][u] < data->x_min)||(ra->data.F32[v][u] > data->x_max)||
-              (dec->data.F32[v][u] < data->y_min)||(dec->data.F32[v][u] > data->y_max)) { continue; }
-          psVectorAppend(obs,image->data.F32[v][u]);
+              (dec->data.F32[v][u] < data->y_min)||(dec->data.F32[v][u] > data->y_max)) { j++; continue; }
+          psVectorAppend(obs,image->data.F32[v][u]);// - camera->data.F32[v][u]);
           psVectorAppend(model, psImageMapEval(data->modelMap,ra->data.F32[v][u],dec->data.F32[v][u]));
-        }
-      }
-
-      psPolynomial1D *poly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD,1);
-      int status = psVectorFitPolynomial1D(poly,NULL,0,model,NULL,obs);
-      if (!status) {
-        psMetadataAddF32(chipItem->data.md,PS_LIST_TAIL,"bkg offset", PS_META_REPLACE, "background offset for this exposure/ota pair", poly->coeff[0]);
-        psMetadataAddF32(chipItem->data.md,PS_LIST_TAIL,"bkg scale", PS_META_REPLACE, "background scale for this exposure/ota pair", poly->coeff[1]);
-      }
-      psFree(poly);
+          j++;
+          used++;
+        }
+      }
+
+      if (used > 0) {
+
+        psPolynomial1D *poly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD,1);
+        int status = psVectorFitPolynomial1D(poly,NULL,0,model,NULL,obs);
+        printf("in model fit loop: %d %d %d %f %f\n",status,j,used,poly->coeff[0],poly->coeff[1]);
+        if (status && (poly->coeff[1] != 0.0)) {
+          psMetadataAddF32(chipItem->data.md,PS_LIST_TAIL,"bkg offset", PS_META_REPLACE, "background offset for this exposure/ota pair", poly->coeff[0]);
+          psMetadataAddF32(chipItem->data.md,PS_LIST_TAIL,"bkg scale", PS_META_REPLACE, "background scale for this exposure/ota pair", poly->coeff[1]);
+        }
+        psFree(poly);
+      }
     } // End OTA loop
     psFree(chipIter);
@@ -153 +176 @@
               (dec->data.F32[v][u] < data->y_min)||(dec->data.F32[v][u] > data->y_max)) { continue; }
 
-          model->data.F32[v][u] = scale * image->data.F32[v][u] - offset - camera->data.F32[v][u];
+          model->data.F32[v][u] = scale * image->data.F32[v][u] + offset - camera->data.F32[v][u];
         }
       }
@@ -168 +191 @@
 // This "averaging" is done using the psImageMapClipFit.
 bool ppBackgroundStackModelFit(ppBackgroundStackData *data) {
-  long j;
-  int u,v;
-
-  long used;
+  long j = 0;
+  int u,v;
+
+  long used = 0;
   psS16 N = psMetadataLookupS16(NULL, data->models, "N");
   psVector *X = psVectorAllocEmpty(N * 13 * 13,PS_TYPE_F32);
@@ -195 +218 @@
       psImage *ra    = psMetadataLookupPtr(NULL, chipItem->data.md, "bkg ra");
       psImage *dec   = psMetadataLookupPtr(NULL, chipItem->data.md, "bkg dec");
-      //      psImage *model = psMetadataLookupPtr(NULL, chipItem->data.md, "bkg image");
+#define DUMP_DATA 0
+#if DUMP_DATA
+      psImage *model = psMetadataLookupPtr(NULL, chipItem->data.md, "bkg image");
+
+      psF32 offset = psMetadataLookupF32(NULL,chipItem->data.md,"bkg offset");
+      psF32 scale  = psMetadataLookupF32(NULL,chipItem->data.md,"bkg scale");
+
+#endif
       for (v = 0; v < calib->numRows; v++) {
         for (u = 0; u < calib->numCols; u++) {
@@ -214 +244 @@
           used++;
           j++;
-          
-/*        printf("DATA %ld %ld %f %f %f %f\n", */
-/*               j,used, */
-/*               ra->data.F32[v][u], */
-/*               dec->data.F32[v][u], */
-/*               calib->data.F32[v][u], */
-/*               model->data.F32[v][u]); */
-                 
+#if DUMP_DATA
+          printf("DATA %d %ld %ld %f %f %f %f %f %f\n",
+                 data->model_iteration,j,used,
+                 offset,scale,
+                 ra->data.F32[v][u],
+                 dec->data.F32[v][u],
+                 calib->data.F32[v][u],
+                 model->data.F32[v][u]);
+#endif 
         }
       }
@@ -230 +261 @@
   bool fitStatus;
   bool status = psImageMapClipFit(&fitStatus,data->modelMap,stats, mask, 1, X, Y, Z, E);
-
+  data->model_iteration++;
   psFree(expIter);
   psFree(X);