Changeset 19231 for trunk/ppStack/src/ppStackMatch.c

Timestamp:

Aug 26, 2008, 4:56:07 PM (18 years ago)

Author:

Paul Price

Message:

Adding command-line option -debug-stack to use previously calculated convolved images and convolution kernel instead of convolving again.

File:

• trunk/ppStack/src/ppStackMatch.c (modified) (5 diffs)

trunk/ppStack/src/ppStackMatch.c

@@ -15 +15 @@
                      PM_SOURCE_MODE_CR_LIMIT) // Mask to apply to input sources
 
-//#define TESTING
+//#define TESTING                         // Enable debugging output
+
+
+
+#ifdef TESTING
+// Read a FITS image
+static bool readImage(psImage **target, // Target for image
+                      const char *name, // Name of FITS file
+                      const pmConfig *config // Configuration
+    )
+{
+    psString resolved = pmConfigConvertFilename(name, config, false, false); // Resolved filename
+    psFits *fits = psFitsOpen(resolved, "r");
+    psFree(resolved);
+    if (!fits) {
+        psError(PS_ERR_IO, false, "Unable to open previously produced image: %s", name);
+        return false;
+    }
+    psImage *image = psFitsReadImage(fits, psRegionSet(0,0,0,0), 0); // Image of interest
+    if (!image) {
+        psError(PS_ERR_IO, false, "Unable to read previously produced image: %s", name);
+        psFitsClose(fits);
+        return false;
+    }
+    psFitsClose(fits);
+
+    psFree(*target);
+    *target = image;
+
+    return true;
+}
+#endif
 
 bool ppStackMatch(pmReadout *readout, psArray **regions, psArray **kernels,
@@ -25 +56 @@
     assert(config);
 
-    pmReadout *output = pmReadoutAlloc(NULL); // Output readout, for holding results temporarily
-
     psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSTACK_RECIPE); // ppStack recipe
     psAssert(recipe, "We've thrown an error on this before.");
@@ -54 +83 @@
     }
 
+#ifdef TESTING
+    // Read previously produced kernel
+    static int numInput = 0;            // Index of input file
+    if (psMetadataLookupBool(NULL, config->arguments, "PPSTACK.DEBUG.STACK")) {
+        const char *outName = psMetadataLookupStr(NULL, config->arguments, "OUTPUT"); // Output root
+        assert(outName);
+        // Read convolution kernel
+        {
+            psString filename = NULL;   // Output filename
+            psStringAppend(&filename, "%s.%d.kernel", outName, numInput);
+            psString resolved = pmConfigConvertFilename(filename, config, false, false); // Resolved filename
+            psFree(filename);
+            psFits *fits = psFitsOpen(resolved, "r"); // FITS file for subtraction kernel
+            psFree(resolved);
+            if (!fits || !pmReadoutReadSubtractionKernels(readout, fits)) {
+                psError(PS_ERR_IO, false, "Unable to read previously produced kernel");
+                psFitsClose(fits);
+                numInput++;
+                return false;
+            }
+            psFitsClose(fits);
+        }
+
+        // Read image, mask, weight
+        const char *tempImage = psMetadataLookupStr(NULL, recipe, "TEMP.IMAGE"); // Suffix for image
+        const char *tempMask = psMetadataLookupStr(NULL, recipe, "TEMP.MASK"); // Suffix for mask
+        const char *tempWeight = psMetadataLookupStr(NULL, recipe, "TEMP.WEIGHT"); // Suffix for weight map
+        psString imageName = NULL, maskName = NULL, weightName = NULL; // Names for convolved images
+        psStringAppend(&imageName, "%s.%d.%s", outName, numInput, tempImage);
+        psStringAppend(&maskName, "%s.%d.%s", outName, numInput, tempMask);
+        psStringAppend(&weightName, "%s.%d.%s", outName, numInput, tempWeight);
+
+        if (!readImage(&readout->image, imageName, config) || !readImage(&readout->mask, maskName, config) ||
+            !readImage(&readout->weight, weightName, config)) {
+            psError(PS_ERR_IO, false, "Unable to read previously produced image.");
+            psFree(imageName);
+            psFree(maskName);
+            psFree(weightName);
+            numInput++;
+            return false;
+        }
+        psFree(imageName);
+        psFree(maskName);
+        psFree(weightName);
+
+        numInput++;
+        return true;
+    }
+#endif
+
+    // Normal operations here
+    pmReadout *output = pmReadoutAlloc(NULL); // Output readout, for holding results temporarily
     if (psMetadataLookupBool(&mdok, config->arguments, "HAVE.PSF")) {
         assert(psf);
@@ -116 +197 @@
         pmReadout *fake = pmReadoutAlloc(NULL); // Fake readout with target PSF
 
-        if (!pmReadoutFakeFromSources(fake, readout->image->numCols, readout->image->numRows, sources, NULL,
-                                      NULL, psf, minFlux, 0, false)) {
+        if (!pmReadoutFakeFromSources(fake, readout->image->numCols, readout->image->numRows, sources,
+                                      NULL, NULL, psf, minFlux, 0, false)) {
             psError(PS_ERR_UNKNOWN, false, "Unable to generate fake image with target PSF.");
             psFree(fake);
@@ -193 +274 @@
 
 #ifdef TESTING
-        {
-            static int num = 0;         // Image counter
-            psString filename = NULL;   // Output filename
-            psStringAppend(&filename, "stack_kernel_%d.fits", num++);
-            psFits *fits = psFitsOpen(filename, "w"); // FITS file for subtraction kernel
-            psFree(filename);
-            pmReadoutWriteSubtractionKernels(output, fits);
-            psFitsClose(fits);
-        }
-#endif
-
-    // Extract the regions and solutions used in the image matching
-    // This stops them from being freed when we iterate back up the FPA
+    {
+        psString filename = NULL;   // Output filename
+        psStringAppend(&filename, "stack_kernel_%d.fits", numInput);
+        psFits *fits = psFitsOpen(filename, "w"); // FITS file for subtraction kernel
+        psFree(filename);
+        pmReadoutWriteSubtractionKernels(output, fits);
+        psFitsClose(fits);
+    }
+#endif
+
+// Extract the regions and solutions used in the image matching
+// This stops them from being freed when we iterate back up the FPA
     *regions = psArrayAllocEmpty(ARRAY_BUFFER); // Array of regions
     {