Changeset 19532 for trunk/psModules/src/imcombine/pmStack.c

Timestamp:

Sep 11, 2008, 6:12:42 PM (18 years ago)

Author:

Paul Price

Message:

Correcting stack variances for the imperfect convolution. Changed (again) the way the deviations are calculated; they're now straight chi^{2, which allows us to estimate how the variance factor needs to be changed in order to correct for it. Moving additional values into pmSubtractionKernels so that this can be calculated from a kernel that has been read in (very useful for debugging, also nice to have a record of the quality of the fit).}

File:

• trunk/psModules/src/imcombine/pmStack.c (modified) (7 diffs)

trunk/psModules/src/imcombine/pmStack.c

@@ -8 +8 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.40 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2008-09-11 20:52:16 $
+ *  @version $Revision: 1.41 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2008-09-12 04:12:42 $
  *  Copyright 2004-2007 Institute for Astronomy, University of Hawaii
  *
@@ -30 +30 @@
 #define PIXEL_LIST_BUFFER 100           // Number of entries to add to pixel list at a time
 #define PIXEL_MAP_BUFFER 2              // Number of entries to add to pixel map at a time
-//#define VARIANCE_FACTORS                // Use variance factors when calculating the variances?
+#define VARIANCE_FACTORS                // Use variance factors when calculating the variances?
 #define NUM_DIRECT_STDEV 5              // For less than this number of values, measure stdev directly
 
@@ -312 +312 @@
               // Use variance to check that the two are consistent
               float diff = pixelData->data.F32[0] - pixelData->data.F32[1];
-#if VARIANCE_FACTORS
+#ifdef VARIANCE_FACTORS
               float sigma2 = pixelVariances->data.F32[0] * varFactors->data.F32[pixelSources->data.U16[0]] +
                   pixelVariances->data.F32[1] * varFactors->data.F32[pixelSources->data.U16[1]];
@@ -346 +346 @@
                   float rej2 = PS_SQR(rej); // Rejection level squared
                   for (int i = 0; i < num; i++) {
-#if VARIANCE_FACTORS
+#ifdef VARIANCE_FACTORS
                       pixelVariances->data.F32[i] *= rej2 * varFactors->data.F32[pixelSources->data.U16[i]];
 #else
@@ -538 +538 @@
 
 /// Constructor
-pmStackData *pmStackDataAlloc(pmReadout *readout, float weight)
+pmStackData *pmStackDataAlloc(pmReadout *readout, float weight, float addVariance)
 {
     pmStackData *data = psAlloc(sizeof(pmStackData)); // Stack data, to return
@@ -547 +547 @@
     data->inspect = NULL;
     data->weight = weight;
+    data->addVariance = addVariance;
 
     return data;
@@ -602 +603 @@
         }
         varFactors->data.F32[i] = vf;
+        if (isfinite(data->addVariance)) {
+            varFactors->data.F32[i] *= data->addVariance;
+        }
         if (!haveRejects && !data->inspect) {
             data->inspect = psPixelsAllocEmpty(PIXEL_LIST_BUFFER);