Changeset 27838 for branches/tap_branches/ppStack/src/ppStackConvolve.c

Timestamp:

May 3, 2010, 8:41:49 AM (16 years ago)

Author:

eugene

Message:

updates from trunk

Location:

branches/tap_branches

Files:

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

branches/tap_branches/ppStack/src

@@ -10 +10 @@
 stamp-h1
 ppStackVersionDefinitions.h
+ppStackErrorCodes.c
+ppStackErrorCodes.h

@@ -10 +10 @@
 stamp-h1
 ppStackVersionDefinitions.h
+ppStackErrorCodes.c
+ppStackErrorCodes.h

branches/tap_branches/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
@@ -52 +59 @@
         pmFPAfileActivate(config->files, false, NULL);
         ppStackFileActivationSingle(config, PPSTACK_FILES_CONVOLVE, true, i);
+        if (options->convolve) {
+            // XXX PPSTACK.CONV.KERNEL not defined unless convolve
+            // pmFPAfileActivate(config->files, true, "PPSTACK.CONV.KERNEL");
+            pmFPAfileActivateSingle(config->files, true, "PPSTACK.CONV.KERNEL", i); // Activated file
+        }
+
         pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPSTACK.INPUT", i); // File of interest
         pmFPAview *view = ppStackFilesIterateDown(config);
@@ -69 +82 @@
         } else if (options->numCols != readout->image->numCols ||
                    options->numRows != readout->image->numRows) {
-            psError(PS_ERR_UNKNOWN, true, "Sizes of input images don't match: %dx%d vs %dx%d",
+            psError(PPSTACK_ERR_ARGUMENTS, true, "Sizes of input images don't match: %dx%d vs %dx%d",
                     readout->image->numCols, readout->image->numRows, options->numCols, options->numRows);
             psFree(rng);
@@ -79 +92 @@
         // 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);
-            options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] |= PPSTACK_MASK_MATCH;
-            psErrorClear();
-            continue;
-        }
-        options->covariances->data[i] = psMemIncrRefCounter(readout->covariance);
+            // XXX many things can cause a failure of ppStackMatch -- should some be handled differently?
+            psErrorCode error = psErrorCodeLast(); // Error code
+            switch (error) {
+                // Fatal errors
+              case PM_ERR_CONFIG:
+              case PPSTACK_ERR_CONFIG:
+              case PPSTACK_ERR_IO:
+                psError(error, false, "Unable to match image %d due to fatal error.", i);
+                return false;
+                // Non-fatal errors
+              case PM_ERR_STAMPS:
+              case PM_ERR_SMALL_AREA:
+              case PPSTACK_ERR_DATA:
+              default:
+                psErrorStackPrint(stderr, "Unable to match image %d --- ignoring.", i);
+                options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] |= PPSTACK_MASK_MATCH;
+                psErrorClear();
+                continue;
+            }
+        }
+        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 +149 @@
         pmHDU *hdu = readout->parent->parent->parent->hdu; // HDU for convolved image
         assert(hdu);
-        ppStackWriteImage(options->imageNames->data[i], hdu->header, readout->image, config);
+        if (!ppStackWriteImage(options->convImages->data[i], hdu->header, readout->image, config)) {
+            psError(PPSTACK_ERR_IO, false, "Unable to write convolved image %d", i);
+            psFree(fpaList);
+            psFree(cellList);
+            psFree(rng);
+            return false;
+        }
         psMetadata *maskHeader = psMetadataCopy(NULL, hdu->header); // Copy of header, for mask
         pmConfigMaskWriteHeader(config, maskHeader);
-        ppStackWriteImage(options->maskNames->data[i], maskHeader, readout->mask, config);
+        if (!ppStackWriteImage(options->convMasks->data[i], maskHeader, readout->mask, config)) {
+            psError(PPSTACK_ERR_IO, false, "Unable to write convolved mask %d", i);
+            psFree(fpaList);
+            psFree(cellList);
+            psFree(rng);
+            psFree(maskHeader);
+            return false;
+        }
         psFree(maskHeader);
-        psImageCovarianceTransfer(readout->variance, readout->covariance);
-        ppStackWriteImage(options->varianceNames->data[i], hdu->header, readout->variance, config);
+        if (!ppStackWriteImage(options->convVariances->data[i], hdu->header, readout->variance, config)) {
+            psError(PPSTACK_ERR_IO, false, "Unable to write convolved variance %d", i);
+            psFree(fpaList);
+            psFree(cellList);
+            psFree(rng);
+            return false;
+        }
 #ifdef TESTING
         {
@@ -129 +182 @@
             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
 
@@ -136 +204 @@
         psListAdd(fpaList, PS_LIST_TAIL, inCell->parent->parent);
 
+        // Correct ZP
+        if (options->matchZPs) {
+            // I think I need to take off the exposure time because we're going to set the new exposure time
+            psMetadataItem *zpItem = psMetadataLookup(inCell->parent->parent->concepts, "FPA.ZP");
+            zpItem->data.F32 += options->norm->data.F32[i] + 2.5*log10(options->sumExposure);
+
+            psMetadataItem *expItem = psMetadataLookup(inCell->parent->parent->concepts, "FPA.EXPOSURE");
+            expItem->data.F32 = options->sumExposure;
+
+            expItem = psMetadataLookup(inCell->concepts, "CELL.EXPOSURE");
+            expItem->data.F32 = options->sumExposure;
+        }
+
         options->cells->data[i] = psMemIncrRefCounter(inCell);
         if (!ppStackFilesIterateUp(config)) {
@@ -149 +230 @@
     psFree(rng);
 
-    psFree(options->norm); options->norm = NULL;
     psFree(options->sourceLists); options->sourceLists = NULL;
     psFree(options->psf); options->psf = NULL;
@@ -155 +235 @@
 
     if (numGood == 0) {
-        psError(PS_ERR_UNKNOWN, false, "No good images to combine.");
+        psError(PPSTACK_ERR_REJECTED, false, "No good images to combine.");
         psFree(fpaList);
         psFree(cellList);
@@ -167 +247 @@
         pmFPAview view;                 // View for output
         view.chip = view.cell = view.readout = 0;
+
         pmCell *outCell = pmFPAfileThisCell(config->files, &view, "PPSTACK.OUTPUT"); // Output cell
         pmFPA *outFPA = outCell->parent->parent; // Output FPA
+
+        pmCell *unconvCell = pmFPAfileThisCell(config->files, &view, "PPSTACK.UNCONV"); // Unconvolved cell
+        pmFPA *unconvFPA = unconvCell->parent->parent;                                  // Unconvolved FPA
+
         pmConceptsAverageFPAs(outFPA, fpaList);
         pmConceptsAverageCells(outCell, cellList, NULL, NULL, true);
+
+        pmConceptsAverageFPAs(unconvFPA, fpaList);
+        pmConceptsAverageCells(unconvCell, cellList, NULL, NULL, true);
+
         psFree(fpaList);
         psFree(cellList);
@@ -177 +266 @@
         UPDATE_CONCEPT(outCell, "CELL.EXPOSURE", options->sumExposure);
         UPDATE_CONCEPT(outCell, "CELL.DARKTIME", NAN);
+        UPDATE_CONCEPT(outFPA,  "FPA.ZP",        options->zp);
+
+        UPDATE_CONCEPT(unconvFPA,  "FPA.EXPOSURE",  options->sumExposure);
+        UPDATE_CONCEPT(unconvCell, "CELL.EXPOSURE", options->sumExposure);
+        UPDATE_CONCEPT(unconvCell, "CELL.DARKTIME", NAN);
+        UPDATE_CONCEPT(unconvFPA,  "FPA.ZP",        options->zp);
     }
 
@@ -191 +286 @@
         assert(values->n == numGood);
         if (!psVectorSortInPlace(values)) {
-            psError(PS_ERR_UNKNOWN, false, "Unable to sort vector.");
+            psError(PPSTACK_ERR_PROG, false, "Unable to sort vector.");
             psFree(values);
             return false;
@@ -211 +306 @@
             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;
                 }
@@ -229 +324 @@
 
     if (numGood == 0) {
-        psError(PS_ERR_UNKNOWN, false, "No good images to combine.");
+        psError(PPSTACK_ERR_REJECTED, false, "No good images to combine.");
         return false;
     }
+
+    // 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;

@@ -10 +10 @@
 stamp-h1
 ppStackVersionDefinitions.h
+ppStackErrorCodes.c
+ppStackErrorCodes.h