Changeset 25788

Timestamp:

Oct 6, 2009, 1:19:44 PM (17 years ago)

Author:

Paul Price

Message:

Implemented unconvolved stack using rejections from convolved stack.

Location:

branches/pap/ppStack/src

Files:

• ppStackCleanup.c (modified) (1 diff)
• ppStackConvolve.c (modified) (1 diff)
• ppStackLoop.c (modified) (4 diffs)
• ppStackOptions.c (modified) (2 diffs)
• ppStackOptions.h (modified) (1 diff)
• ppStackSetup.c (modified) (3 diffs)
• ppStackThread.c (modified) (7 diffs)
• ppStackThread.h (modified) (1 diff)

branches/pap/ppStack/src/ppStackCleanup.c

@@ -74 +74 @@
 
         if (tempDelete) {
-            psString imageResolved = pmConfigConvertFilename(options->imageNames->data[i],
+            psString imageResolved = pmConfigConvertFilename(options->convImages->data[i],
                                                              config, false, false);
-            psString maskResolved = pmConfigConvertFilename(options->maskNames->data[i],
+            psString maskResolved = pmConfigConvertFilename(options->convMasks->data[i],
                                                             config, false, false);
-            psString varianceResolved = pmConfigConvertFilename(options->varianceNames->data[i],
+            psString varianceResolved = pmConfigConvertFilename(options->convVariances->data[i],
                                                                 config, false, false);
             if (unlink(imageResolved) == -1 || unlink(maskResolved) == -1 ||

branches/pap/ppStack/src/ppStackConvolve.c

@@ -114 +114 @@
         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
         {

branches/pap/ppStack/src/ppStackLoop.c

@@ -56 +56 @@
     // Start threading
     ppStackThreadInit();
-    ppStackThreadData *stack = ppStackThreadDataSetup(options, config);
+    ppStackThreadData *stack = ppStackThreadDataSetup(options, config, true);
     if (!stack) {
         psError(PS_ERR_IO, false, "Unable to initialise stack threads.");
@@ -122 +122 @@
 
 
+#if 1
+    // Unconvolved stack --- it's cheap to calculate, compared to everything else!
+    if (options->convolve) {
+        // Start threading
+        ppStackThreadData *stack = ppStackThreadDataSetup(options, config, false);
+        if (!stack) {
+            psError(PS_ERR_IO, false, "Unable to initialise stack threads.");
+            psFree(options);
+            return false;
+        }
+        psFree(options->cells); options->cells = NULL;
+        psTrace("ppStack", 2, "Final stack of convolved images....\n");
+        if (!ppStackCombineFinal(stack, options, config)) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to perform final combination.");
+            psFree(stack);
+            psFree(options);
+            return false;
+        }
+        psLogMsg("ppStack", PS_LOG_INFO, "Stage 7: Unconvolved Stack: %f sec", psTimerClear("PPSTACK_STEPS"));
+        ppStackMemDump("unconv");
+
+        psFree(stack);
+    }
+#endif
+
     // Photometry
     psTrace("ppStack", 1, "Photometering stacked image....\n");
@@ -129 +154 @@
         return false;
     }
-    psLogMsg("ppStack", PS_LOG_INFO, "Stage 7: Photometry Analysis: %f sec", psTimerClear("PPSTACK_STEPS"));
+    psLogMsg("ppStack", PS_LOG_INFO, "Stage 8: Photometry Analysis: %f sec", psTimerClear("PPSTACK_STEPS"));
     ppStackMemDump("photometry");
 
@@ -140 +165 @@
         return false;
     }
-    psLogMsg("ppStack", PS_LOG_INFO, "Stage 8: Final output: %f sec", psTimerClear("PPSTACK_STEPS"));
+    psLogMsg("ppStack", PS_LOG_INFO, "Stage 9: Final output: %f sec", psTimerClear("PPSTACK_STEPS"));
     ppStackMemDump("finish");
 

branches/pap/ppStack/src/ppStackOptions.c

@@ -12 +12 @@
         fclose(options->statsFile);
     }
-    psFree(options->imageNames);
-    psFree(options->maskNames);
-    psFree(options->varianceNames);
+    psFree(options->origImages);
+    psFree(options->convImages);
+    psFree(options->convMasks);
+    psFree(options->convVariances);
     psFree(options->psf);
     psFree(options->inputMask);
@@ -40 +41 @@
     options->stats = NULL;
     options->statsFile = NULL;
-    options->imageNames = NULL;
-    options->maskNames = NULL;
-    options->varianceNames = NULL;
+    options->origImages = NULL;
+    options->convImages = NULL;
+    options->convMasks = NULL;
+    options->convVariances = NULL;
     options->num = 0;
     options->psf = NULL;

branches/pap/ppStack/src/ppStackOptions.h

@@ -11 +11 @@
     psMetadata *stats;                  // Statistics for output
     FILE *statsFile;                    // File to which to write statistics
-    psArray *imageNames, *maskNames, *varianceNames; // Filenames for the temporary convolved images
+    psArray *origImages;                // Filenames of the original images
+    psArray *convImages, *convMasks, *convVariances; // Filenames for the temporary convolved images
     int num;                            // Number of inputs
     // Prepare

branches/pap/ppStack/src/ppStackSetup.c

@@ -46 +46 @@
     }
 
+    // Original images
+    options->origImages = psArrayAlloc(num);
+    for (int i = 0; i < num; i++) {
+        pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPSTACK.INPUT", i); // File of interest
+        options->origImages->data[i] = psStringCopy(file->name);
+    }
+
     // Generate temporary names for convolved images
     const char *tempDir = psMetadataLookupStr(NULL, recipe, "TEMP.DIR"); // Directory for temporary images
@@ -72 +79 @@
     }
 
-    options->imageNames = psArrayAlloc(num);
-    options->maskNames = psArrayAlloc(num);
-    options->varianceNames = psArrayAlloc(num);
+    options->convImages = psArrayAlloc(num);
+    options->convMasks = psArrayAlloc(num);
+    options->convVariances = psArrayAlloc(num);
     for (int i = 0; i < num; i++) {
         psString imageName = NULL, maskName = NULL, varianceName = NULL; // Names for convolved images
@@ -81 +88 @@
         psStringAppend(&varianceName, "%s/%s.%d.%s", tempDir, tempName, i, tempVariance);
         psTrace("ppStack", 5, "Temporary files: %s %s %s\n", imageName, maskName, varianceName);
-        options->imageNames->data[i] = imageName;
-        options->maskNames->data[i] = maskName;
-        options->varianceNames->data[i] = varianceName;
+        options->convImages->data[i] = imageName;
+        options->convMasks->data[i] = maskName;
+        options->convVariances->data[i] = varianceName;
     }
     psFree(outputName);

branches/pap/ppStack/src/ppStackThread.c

@@ -56 +56 @@
 }
 
-ppStackThreadData *ppStackThreadDataSetup(const ppStackOptions *options, const pmConfig *config)
+ppStackThreadData *ppStackThreadDataSetup(const ppStackOptions *options, const pmConfig *config, bool conv)
 {
     psAssert(options, "Require options");
@@ -62 +62 @@
 
     const psArray *cells = options->cells; // Array of input cells
-    const psArray *imageNames = options->imageNames; // Names of images to read
-    const psArray *maskNames = options->maskNames; // Names of masks to read
-    const psArray *varianceNames = options->varianceNames; // Names of variance maps to read
-    const psArray *covariances = options->covariances; // Covariance matrices (already read)
+    const psArray *imageNames = conv ? options->convImages : options->origImages; // Names of images to read
+    const psArray *maskNames = conv ? options->convMasks : NULL; // Names of masks to read
+    const psArray *varianceNames = conv ? options->convVariances : NULL; // Names of variance maps to read
+    const psArray *covariances = conv ? options->covariances : NULL; // Covariance matrices (already read)
 
     PS_ASSERT_ARRAY_NON_NULL(cells, NULL);
     PS_ASSERT_ARRAYS_SIZE_EQUAL(cells, imageNames, NULL);
-    PS_ASSERT_ARRAYS_SIZE_EQUAL(cells, maskNames, NULL);
-    PS_ASSERT_ARRAYS_SIZE_EQUAL(cells, varianceNames, NULL);
-    PS_ASSERT_ARRAYS_SIZE_EQUAL(cells, covariances, NULL);
+    if (conv) {
+        PS_ASSERT_ARRAYS_SIZE_EQUAL(cells, maskNames, NULL);
+        PS_ASSERT_ARRAYS_SIZE_EQUAL(cells, varianceNames, NULL);
+        PS_ASSERT_ARRAYS_SIZE_EQUAL(cells, covariances, NULL);
+    }
 
     ppStackThreadData *stack = psAlloc(sizeof(ppStackThreadData)); // Thread data, to return
@@ -87 +89 @@
             continue;
         }
-        // Resolved names
-        psString imageResolved = pmConfigConvertFilename(imageNames->data[i], config, false, false);
-        psString maskResolved = pmConfigConvertFilename(maskNames->data[i], config, false, false);
-        psString varianceResolved = pmConfigConvertFilename(varianceNames->data[i], config, false, false);
-        stack->imageFits->data[i] = psFitsOpen(imageResolved, "r");
-        stack->maskFits->data[i] = psFitsOpen(maskResolved, "r");
-        stack->varianceFits->data[i] = psFitsOpen(varianceResolved, "r");
-        psFree(imageResolved);
-        psFree(maskResolved);
-        psFree(varianceResolved);
-        if (!stack->imageFits->data[i] || !stack->maskFits->data[i] || !stack->varianceFits->data[i]) {
-            psError(PS_ERR_UNKNOWN, false, "Unable to open convolved files %s, %s, %s",
-                    (char*)imageNames->data[i], (char*)maskNames->data[i], (char*)varianceNames->data[i]);
-            return NULL;
-        }
+
+// Open an image
+#define IMAGE_OPEN(NAME, FITS) \
+        if (NAME) { \
+            psString resolved = pmConfigConvertFilename(NAME, config, false, false); \
+            FITS = psFitsOpen(resolved, "r"); \
+            if (!(FITS)) { \
+                psError(PS_ERR_IO, false, "Unable to open file %s", (char*)(NAME)); \
+                psFree(resolved); \
+                return NULL; \
+            } \
+            psFree(resolved); \
+        }
+
+        IMAGE_OPEN(imageNames->data[i], stack->imageFits->data[i]);
+        IMAGE_OPEN(maskNames->data[i], stack->maskFits->data[i]);
+        IMAGE_OPEN(varianceNames->data[i], stack->varianceFits->data[i]);
     }
 
@@ -116 +120 @@
             }
             pmReadout *ro = pmReadoutAlloc(cell); // Readout for thread
-            ro->covariance = psMemIncrRefCounter(covariances->data[j]);
+            if (conv) {
+                ro->covariance = psMemIncrRefCounter(covariances->data[j]);
+            }
             readouts->data[j] = ro;
         }
@@ -187 +193 @@
 
 
-                int zMax = 0;
+                int zMax = 0;
                 bool keepReading = false;
-                if (pmReadoutMore(ro, imageFits, 0, &zMax, rows, config)) {
+
+                if (imageFits && pmReadoutMore(ro, imageFits, 0, &zMax, rows, config)) {
                     keepReading = true;
                     if (!pmReadoutReadChunk(ro, imageFits, 0, NULL, rows, overlap, config)) {
@@ -199 +206 @@
                 }
 
-                if (pmReadoutMoreMask(ro, maskFits, 0, &zMax, rows, config)) {
+                if (maskFits && pmReadoutMoreMask(ro, maskFits, 0, &zMax, rows, config)) {
                     keepReading = true;
                     if (!pmReadoutReadChunkMask(ro, maskFits, 0, NULL, rows, overlap, config)) {
@@ -209 +216 @@
                 }
 
-                if (pmReadoutMoreVariance(ro, varianceFits, 0, &zMax, rows, config)) {
+                if (varianceFits && pmReadoutMoreVariance(ro, varianceFits, 0, &zMax, rows, config)) {
                     keepReading = true;
                     if (!pmReadoutReadChunkVariance(ro, varianceFits, 0, NULL, rows, overlap, config)) {

branches/pap/ppStack/src/ppStackThread.h

@@ -35 +35 @@
 ppStackThreadData *ppStackThreadDataSetup(
     const ppStackOptions *options,      // Options
-    const pmConfig *config              // Configuration
+    const pmConfig *config,             // Configuration
+    bool conv                           // Use convolved products?
     );