Changeset 26076 for trunk/ppStack/src/ppStackConvolve.c

Timestamp:

Nov 9, 2009, 2:53:12 PM (17 years ago)

Author:

Paul Price

Message:

Merging branches/pap (unconvolved stacks, reworked combinePixels function into clearer roles, only throw out most variant pixel on each iteration, throw suspect pixels out as first rejection step) stack development. I think I've got everything, but not entirely sure, since I've already merged this branch once before (for dual convolution).

Location:

trunk/ppStack

Files:

• . (modified) (1 prop)
• src/ppStackConvolve.c (modified) (9 diffs)

trunk/ppStack/src/ppStackConvolve.c

@@ -9 +9 @@
 #include "ppStack.h"
 #include "ppStackLoop.h"
+
+//#define TESTING
+
 
 // Update the value of a concept
@@ -39 +42 @@
     options->weightings = psVectorAlloc(num, PS_TYPE_F32); // Combination weightings for images (1/noise^2)
     psVectorInit(options->weightings, 0.0);
-    options->covariances = psArrayAlloc(num); // Covariance matrices
+    options->origCovars = psArrayAlloc(num);
+    options->convCovars = psArrayAlloc(num); // Covariance matrices
+
+    psVector *renorms = psVectorAlloc(num, PS_TYPE_F32); // Renormalisation values for variances
+    psVectorInit(renorms, NAN);
 
     psList *fpaList = psListAlloc(NULL); // List of input FPAs, for concept averaging
@@ -79 +86 @@
         // Background subtraction, scaling and normalisation is performed automatically by the image matching
         psTimerStart("PPSTACK_MATCH");
+        options->origCovars->data[i] = psMemIncrRefCounter(readout->covariance);
         if (!ppStackMatch(readout, options, i, config)) {
             psErrorStackPrint(stderr, "Unable to match image %d --- ignoring.", i);
@@ -85 +93 @@
             continue;
         }
-        options->covariances->data[i] = psMemIncrRefCounter(readout->covariance);
+        options->convCovars->data[i] = psMemIncrRefCounter(readout->covariance);
+
+        float renorm = psMetadataLookupF32(NULL, readout->analysis, PM_READOUT_ANALYSIS_RENORM);
+        if (!isfinite(renorm)) {
+            renorm = 1.0;
+        }
+        renorms->data.F32[i] = renorm;
 
         if (options->stats) {
@@ -114 +128 @@
         pmHDU *hdu = readout->parent->parent->parent->hdu; // HDU for convolved image
         assert(hdu);
-        ppStackWriteImage(options->imageNames->data[i], hdu->header, readout->image, config);
+        ppStackWriteImage(options->convImages->data[i], hdu->header, readout->image, config);
         psMetadata *maskHeader = psMetadataCopy(NULL, hdu->header); // Copy of header, for mask
         pmConfigMaskWriteHeader(config, maskHeader);
-        ppStackWriteImage(options->maskNames->data[i], maskHeader, readout->mask, config);
+        ppStackWriteImage(options->convMasks->data[i], maskHeader, readout->mask, config);
         psFree(maskHeader);
-        psImageCovarianceTransfer(readout->variance, readout->covariance);
-        ppStackWriteImage(options->varianceNames->data[i], hdu->header, readout->variance, config);
+        ppStackWriteImage(options->convVariances->data[i], hdu->header, readout->variance, config);
 #ifdef TESTING
         {
@@ -129 +142 @@
             psFree(name);
         }
+        {
+            int numCols = readout->image->numCols, numRows = readout->image->numRows;
+            psImage *sn = psImageAlloc(numCols, numRows, PS_TYPE_F32);
+            for (int y = 0; y < numRows; y++) {
+                for (int x = 0; x < numCols; x++) {
+                    sn->data.F32[y][x] = readout->image->data.F32[y][x] /
+                        sqrtf(readout->variance->data.F32[y][x]);
+                }
+            }
+            psString name = NULL;
+            psStringAppend(&name, "signoise_%d.fits", i);
+            ppStackWriteImage(name, hdu->header, sn, config);
+            psFree(name);
+            psFree(sn);
+        }
 #endif
 
@@ -149 +177 @@
     psFree(rng);
 
-    psFree(options->norm); options->norm = NULL;
     psFree(options->sourceLists); options->sourceLists = NULL;
     psFree(options->psf); options->psf = NULL;
@@ -211 +238 @@
             numGood = 0;                    // Number of good images
             for (int i = 0; i < num; i++) {
-              if (options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PPSTACK_MASK_ALL) {
+                if (options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PPSTACK_MASK_ALL) {
                     continue;
                 }
@@ -233 +260 @@
     }
 
+    // Correct chi^2 for renormalisation
+    psBinaryOp(options->matchChi2, options->matchChi2, "/", renorms);
+    for (int i = 0; i < num; i++) {
+        psLogMsg("ppStack", PS_LOG_INFO, "Additional variance for image %d: %f\n",
+                 i, options->matchChi2->data.F32[i]);
+    }
+    psFree(renorms);
+
     return true;
 }