Changeset 2661 for trunk/stac/src/stacCombine.c

Timestamp:

Dec 7, 2004, 3:38:05 PM (22 years ago)

Author:

Paul Price

Message:

Working, but not getting faint CRs

File:

• trunk/stac/src/stacCombine.c (modified) (6 diffs)

trunk/stac/src/stacCombine.c

@@ -16 +16 @@
     psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV);
     (void)psVectorStats(stats, values, masks, 1);
+    float mean = stats->sampleMean;
     psFree(stats);
-    return stats->sampleMean;
+    return mean;
 #else
     // Instead, do it ourselves
@@ -45 +46 @@
 {
     psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN);
-    (void)psVectorStats(stats, values, masks, 1);
+    (void)psVectorStats(stats, values, masks, 0);
+    float median = stats->sampleMedian;
     psFree(stats);
-    return stats->sampleMedian;
+    return median;
 }
 
@@ -113 +115 @@
         }
     }
-    
+
+#ifdef TESTING
+    // chi^2 image
+    psImage *chi2 = psImageAlloc(numCols, numRows, PS_TYPE_F32);
+    static int iteration = 0;           // Number of times function has been called
+#endif
     
     for (int y = 0; y < numRows; y++) {
@@ -131 +138 @@
             }
             
-            float average = stacCombineMean(pixels, deltas, mask); // Combined value
+            float average = stacCombineMedian(pixels, deltas, mask); // Combined value
+
+#ifdef TESTING
+            // Calculate chi^2
+            chi2->data.F32[y][x] = 0.0;
+            int numGoodPix = 0;
+            for (int i = 0; i < nImages; i++) {
+                if (mask->data.U8[i]) {
+                    chi2->data.F32[y][x] += SQUARE((pixels->data.F32[i] - average) / deltas->data.F32[i]);
+                    numGoodPix++;
+                }
+            }
+            chi2->data.F32[y][x] /= (float)numGoodPix;
+#endif
             
             // Rejection iterations
             for (int rejNum = 0; rejNum < nReject; rejNum++) {
                 float max = 0.0;
-                int maxIndex = 0;
+                int maxIndex = -1;
                 for (int i = 0; i < nImages; i++) {
-                    if (mask->data.U8[i] && (ABS(pixels->data.F32[i] - average) / deltas->data.F32[i] > max)) {
-                        max = ABS(pixels->data.F32[i] - average) / deltas->data.F32[i];
+                    if (mask->data.U8[i] && ((pixels->data.F32[i] - average) / deltas->data.F32[i] > max)) {
+                        max = (pixels->data.F32[i] - average) / deltas->data.F32[i];
                         maxIndex = i;
                     }
@@ -152 +172 @@
             } // Rejection iterations
             
-            combined->data.F32[y][x] = average;
+            combined->data.F32[y][x] = stacCombineMean(pixels, deltas, mask);
         }
     } // Iterating over output pixels
@@ -160 +180 @@
     if (nReject > 0) {
         for (int i = 0; i < nImages; i++) {
-            char rejfile[MAXCHAR];              // Filename of rejection image
+            char rejfile[MAXCHAR];      // Filename of rejection image
             sprintf(rejfile,"%s.shiftrej",config->inputs->data[i]);
             psImageWriteSection((psImage*)((*rejected)->data[i]),0,0,0,NULL,0,rejfile);
         }
     }
+    // Write chi^2 image
+    iteration++;
+    char chifile[MAXCHAR];              // Filename of chi^2 image
+    sprintf(chifile,"chi2_%d.fits",iteration);
+    psImageWriteSection(chi2,0,0,0,NULL,0,chifile);
+    psFree(chi2);
 #endif