Changeset 26690 for branches/eam_branches/20091201/ppStack/src/ppStackMatch.c

Timestamp:

Jan 26, 2010, 11:22:22 PM (16 years ago)

Author:

Paul Price

Message:

Renormalise image so that the normalisation provided by the convolution kernel matches that of determined from stellar photometry.

File:

• branches/eam_branches/20091201/ppStack/src/ppStackMatch.c (modified) (5 diffs)

branches/eam_branches/20091201/ppStack/src/ppStackMatch.c

@@ -18 +18 @@
 #define COVAR_FRAC 0.01                 // Truncation fraction for covariance matrix
 
-//#define TESTING                         // Enable debugging output
+#define TESTING                         // Enable debugging output
 
 #ifdef TESTING
@@ -274 +274 @@
                                                                 PM_SUBTRACTION_ANALYSIS_KERNEL);
 
-            pmSubtractionAnalysis(readout->analysis, NULL, kernels, region,
+            pmSubtractionAnalysis(conv->analysis, NULL, kernels, region,
                                   readout->image->numCols, readout->image->numRows);
 
@@ -372 +372 @@
             fake->mask = psImageCopy(NULL, readout->mask, PS_TYPE_IMAGE_MASK);
 
+#if 1
             // Add the background into the target image
             psImage *bgImage = stackBackgroundModel(readout, config); // Image of background
             psBinaryOp(fake->image, fake->image, "+", bgImage);
             psFree(bgImage);
+#endif
 
 #ifdef TESTING
@@ -456 +458 @@
                     return false;
                 }
-                                              psFree(widthsCopy);
+                psFree(widthsCopy);
             }
 
@@ -559 +561 @@
         }
 
+        // Kernel normalisation
+        {
+            double sum = 0.0;           // Sum of chi^2
+            int num = 0;                // Number of measurements of chi^2
+            psString regex = NULL;      // Regular expression
+            psStringAppend(&regex, "^%s$", PM_SUBTRACTION_ANALYSIS_NORM);
+            psMetadataIterator *iter = psMetadataIteratorAlloc(conv->analysis, PS_LIST_HEAD, regex);
+            psFree(regex);
+            psMetadataItem *item = NULL;// Item from iteration
+            while ((item = psMetadataGetAndIncrement(iter))) {
+                assert(item->type == PS_TYPE_F32);
+                float norm = item->data.F32; // Normalisation
+                sum += norm;
+                num++;
+            }
+            psFree(iter);
+            float conv = sum/num;       // Mean normalisation from convolution
+            float stars = powf(10.0, -0.4 * options->norm->data.F32[index]); // Normalisation from stars
+            float renorm =  stars / conv; // Renormalisation to apply
+            psLogMsg("ppStack", PS_LOG_INFO, "Renormalising image %d by %f (kernel: %f, stars: %f)\n",
+                     index, renorm, conv, stars);
+            psBinaryOp(readout->image, readout->image, "*", psScalarAlloc(renorm, PS_TYPE_F32));
+            psBinaryOp(readout->variance, readout->variance, "*", psScalarAlloc(PS_SQR(renorm), PS_TYPE_F32));
+        }
+
         // Reject image completely if the maximum deconvolution fraction exceeds the limit
         float deconv = psMetadataLookupF32(NULL, conv->analysis,