Changeset 31713 for branches/haf_branches/ppimage2011/psvideophot/src/psvideoLoop.c

Timestamp:

Jun 24, 2011, 2:18:18 PM (15 years ago)

Author:

bills

Message:

Adapt Pavlos' code to work within the IPP framework and include it in psvideoLoop.c
Gets results similar to his original version when the star is fully within the
boundaries of the frame. His version examined pixels in the BIAS region of the chip

File:

• branches/haf_branches/ppimage2011/psvideophot/src/psvideoLoop.c (modified) (7 diffs)

branches/haf_branches/ppimage2011/psvideophot/src/psvideoLoop.c

@@ -4 +4 @@
 
 #include "psvideophot.h"
+
+// XXX: get these from a recipe
+#define DX 20
+#define DY 10
+#define SKIP 0
+#define EDGE 3
+#define THRESHOLD 100
 
 #define ESCAPE(MESSAGE) { \
@@ -12 +19 @@
 }
 
+// XXX: move these to another file 
+static double SubSum (psImage *image, double xy[2], double median);
+static void Centering(psImage *image, long edge, double *pSum, double xy[], double median);
+
 bool psvideophotLoop(pmConfig *config, psvideophotOptions *options)
 {
     psMetadata *stats = NULL;           // Statistics to output
-    float timeDetrend = 0;              // Amount of time spent in detrend
     float timePhot = 0;                 // Amount of time spent in photometry
 
     if (options->doStats) {
-        stats = psMetadataAlloc();
-        psMetadataAddS32(stats, PS_LIST_TAIL, "QUALITY", 0, "No problems", 0);
+        psMetadata *statsMD = psMetadataAlloc();
+        psMetadataAddS32(statsMD, PS_LIST_TAIL, "QUALITY", 0, "No problems", 0);
     }
 
@@ -33 +43 @@
     }
 
+
     pmConfigCamerasCull(config, NULL);
     pmConfigRecipesCull(config, "PSVIDEOPHOT,PPSTATS,PSPHOT,MASKS,PSASTRO,JPEG");
@@ -54 +65 @@
         }
 
-
-        psTimerStart(TIMER_DETREND);
+        double threshold = THRESHOLD;
+
+        psTimerStart(TIMER_PHOT);
         pmCell *cell;                   // Cell from chip
         while ((cell = pmFPAviewNextCell(view, input->fpa, 1)) != NULL) {
@@ -73 +85 @@
             }
 
-            // XXX for now, skip the video cells (cell->readouts->n > 1)
             if (cell->readouts->n <2) {
-              psWarning("Skipping Non-Video Cell for psvideophotDetrendReadout");
               continue; 
             }
-            psWarning("Found video cell\n");
+
+            psString cellName = psMetadataLookupStr(&status, cell->concepts, "CELL.NAME");
+            psWarning("Found video cell %s %d\n", cellName, view->cell);
+
+            psStats *stats = psStatsAlloc(
+                PS_STAT_SAMPLE_MEDIAN | 
+                PS_STAT_SAMPLE_STDEV |
+                PS_STAT_ROBUST_MEDIAN | 
+                PS_STAT_ROBUST_STDEV |
+                PS_STAT_MIN | PS_STAT_MAX | PS_STAT_SAMPLE_QUARTILE 
+                | PS_STAT_SAMPLE_SKEWNESS | PS_STAT_SAMPLE_KURTOSIS
+                );
 
             // process each of the readouts
+            int frame = 0;
             pmReadout *readout;         // Readout from cell
             while ((readout = pmFPAviewNextReadout (view, input->fpa, 1)) != NULL) {
-              if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
+                if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
                     ESCAPE("load failure for Readout");
-              }
-              
+                }
+                if (++frame < SKIP) {
+                    continue;
+                }
                 if (!readout->data_exists) {
                     continue;
                 }
-                psWarning ( "Cell %d: %x %x\n", view->cell, cell->file_exists, cell->process);
-            }
-            
-            
-
-        }
-        // free detrend images potentially in use: MASK, BIAS, DARK, SHUTTER, FLAT
-        timeDetrend += psTimerClear(TIMER_DETREND);
+
+                psStatsInit(stats);
+
+                if (!psImageStats(stats, readout->image, NULL, 0)) {
+                    psError(PS_ERR_UNKNOWN, false, "unable to generate image stats for cell %d readout %d\n",
+                        view->cell, view->readout);
+                    ESCAPE("failed to generate stats");
+                }
+
+                double stddev = stats->sampleStdev;
+                double median = stats->sampleMedian;
+                if (stddev > threshold) {
+
+                    double sum = 0;
+                    // XXX: use a point data type
+                    double xy[2];
+                    Centering(readout->image, EDGE, &sum, xy, median);
+                    double sum2;
+                    sum2 = SubSum(readout->image, xy, median);
+
+                    printf("%4d %9.2f %9.2f %8.2f %.1f %.1f %9.2f\n",
+                            frame, sum, stddev, median, xy[0], xy[1], sum2);
+                } else {
+                    printf("%4d %9.2f %9.2f %8.2f\n",
+                            frame, 0.0,   stddev, median);
+                }
+            }
+        }
 
         // we perform photometry on the readouts of this chip in the output
 
-        psTimerStart(TIMER_PHOT);
         timePhot += psTimerClear(TIMER_PHOT);
 
         view->cell = -1;
-        psWarning ("Got here\n");
-        //  while ((cell = pmFPAviewNextCell(view, input->fpa, 1)) != NULL) {
-        //  psWarning ( "Cell %d: %x %x\n", view->cell, cell->file_exists, cell->process);
-        //  if (!cell->process || !cell->file_exists) {
-        //       continue;
-        //   }
-        //   psWarning ("Got here: inner loop\n");
-            
-        //    if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
-        //       ESCAPE("save failure for Cell");
-        //   }
-        //}
 
         // Close chip
@@ -123 +154 @@
         }
     }
+    psFree (stats);
     
     if (psTraceGetLevel("psvideophot") >= 3) {
@@ -141 +173 @@
     }
 
-
-    psFree (stats);
-
     return true;
 }
+
+// calculate the sum of pixel values inside the window
+static double Sum(psImage *image, long edge, double median) {
+    long        x, y;
+    double      sum=0;
+    for (y=edge; y<=image->numRows-edge+1; y++) {
+        // TYPE!
+        psF32 *row = image->data.F32[y];
+        for (x=edge; x<= image->numCols - edge+1; x++) {
+            sum = sum + row[x] - median;
+        }
+    }   
+    return(sum);
+}
+
+// find the star
+static void Centering(psImage *image, long edge, double *pSum, double xy[], double median) {
+    long        x, y;
+    double      sum=0, sxf=0, syf=0;
+        
+    sum=Sum(image, edge, median);
+
+    for (y=edge; y<= image->numRows -edge+1; y++) {
+        // TYPE!!
+        psF32 *row = image->data.F32[y];
+        for (x=edge; x<=image->numCols -edge+1; x++) {
+            sxf=sxf+x*(row[x] - median);
+            syf=syf+y*(row[x] - median);
+        }
+    }
+    xy[0] = edge-1+sxf/sum;
+    xy[1] = edge-2+syf/sum;
+    *pSum = sum;
+    //printf("%f %f\n", xc, yc);
+    //return(0);
+}
+
+// calculate the sum of pixel values in an aperture around a point
+double SubSum (psImage *image, double xy[2], double median)
+{
+    double      sum=0;
+    long        x, y;
+
+    for (y = xy[1] - DY; y<= xy[1] + DY; y++) {
+        if (y < image->row0) {
+            continue;
+        }
+        if (y >= image->numRows) {
+            break;
+        }
+        psF32 *row = image->data.F32[y];
+        for (x = xy[0] - DX; x<= xy[0] + DX; x++) {
+            if (x < image->col0) {
+                continue;
+            }
+            if (x >= image->numCols) {
+                break;
+            }
+            sum += row[x] - median;
+        }
+    }   
+    return(sum);
+}