Changeset 14457 for trunk/ppSub/src/ppSubReadout.c

Timestamp:

Aug 9, 2007, 2:11:07 PM (19 years ago)

Author:

Paul Price

Message:

Adding support for doing subtraction in multiple regions with independent solutions.

File:

• trunk/ppSub/src/ppSubReadout.c (modified) (5 diffs)

trunk/ppSub/src/ppSubReadout.c

@@ -10 +10 @@
 #include "ppSub.h"
 
-#define TESTING
+//#define TESTING
 #define WCS_TOLERANCE 0.001             // Tolerance for WCS
 
@@ -52 +52 @@
     int size = psMetadataLookupS32(NULL, config->arguments, "KERNEL.SIZE"); // Kernel half-size
     int order = psMetadataLookupS32(NULL, config->arguments, "SPATIAL.ORDER"); // Spatial polynomial order
+    float regionSize = psMetadataLookupF32(NULL, config->arguments, "REGION.SIZE"); // Size of iso-kernel regs
     float spacing = psMetadataLookupF32(NULL, config->arguments, "STAMP.SPACING"); // Typical stamp spacing
     int footprint = psMetadataLookupS32(NULL, config->arguments, "STAMP.FOOTPRINT"); // Stamp half-size
@@ -119 +120 @@
     memCheck("kernels");
 
-    int numRejected = -1;               // Number of rejected stamps in each iteration
-    for (int i = 0; i < iter && numRejected != 0; i++) {
-        stamps = pmSubtractionFindStamps(stamps, refRO->image, subMask, threshold, spacing);
-        if (!stamps) {
-            psError(PS_ERR_UNKNOWN, false, "Unable to find stamps on reference image.");
-            goto ERROR;
-        }
-
-        memCheck("  find stamps");
-
-        if (!pmSubtractionCalculateEquation(stamps, refRO->image, inRO->image, inRO->weight,
-                                            kernels, footprint)) {
-            psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation.");
-            goto ERROR;
-         }
-
-        memCheck("  calculate equation");
-
-        solution = pmSubtractionSolveEquation(solution, stamps);
-        if (!solution) {
-            psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation.");
-            goto ERROR;
-        }
-
-        memCheck("  solve equation");
-
-        numRejected = pmSubtractionRejectStamps(stamps, refRO->image, inRO->image, subMask,
-                                                solution, footprint, rej, kernels);
-        if (numRejected < 0) {
-            psError(PS_ERR_UNKNOWN, false, "Unable to reject stamps.");
-            goto ERROR;
-        }
-        psLogMsg("ppSub", PS_LOG_INFO, "%d stamps rejected on iteration %d.", numRejected, i);
-
-        memCheck("  reject stamps");
-    }
-
-    if (numRejected > 0) {
-        solution = pmSubtractionSolveEquation(solution, stamps);
-        if (!solution) {
-            psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation.");
-            goto ERROR;
-        }
-    }
-    psFree(stamps);
-    stamps = NULL;
-
-    memCheck("solution");
+    int xRegions = 1, yRegions = 1;     // Number of iso-kernel regions
+    float xRegionSize = 0, yRegionSize = 0; // Size of iso-kernel regions
+    psRegion *region = NULL;            // Iso-kernel region
+    if (isfinite(regionSize)) {
+        xRegions = numCols / regionSize + 1;
+        yRegions = numRows / regionSize + 1;
+        xRegionSize = (float)numCols / (float)xRegions;
+        yRegionSize = (float)numRows / (float)yRegions;
+        region = psRegionAlloc(NAN, NAN, NAN, NAN);
+    }
+
+    psImage *convImage = NULL, *convWeight = NULL, *convMask = NULL; // Convolved images
+
+    // Iterate over iso-kernel regions
+    for (int j = 0; j < yRegions; j++) {
+        for (int i = 0; i < xRegions; i++) {
+            psTrace("ppSub", 1, "Subtracting region %d of %d...\n",
+                    j * xRegions + i + 1, xRegions * yRegions);
+            if (region) {
+                *region = psRegionSet((int)(i * xRegionSize), (int)((i + 1) * xRegionSize),
+                                      (int)(j * yRegionSize), (int)((j + 1) * yRegionSize));
+                psString string = psRegionToString(*region);
+                psTrace("ppSub", 3, "Iso-kernel region: %s out of %d,%d\n", string, numCols, numRows);
+                psFree(string);
+            }
+
+            int numRejected = -1;               // Number of rejected stamps in each iteration
+            for (int k = 0; k < iter && numRejected != 0; k++) {
+                psTrace("ppSub", 2, "Iteration %d...\n", k);
+                psTrace("ppSub", 3, "Finding stamps...\n");
+                stamps = pmSubtractionFindStamps(stamps, refRO->image, subMask, region, threshold, spacing);
+                if (!stamps) {
+                    psError(PS_ERR_UNKNOWN, false, "Unable to find stamps on reference image.");
+                    goto ERROR;
+                }
+
+                memCheck("  find stamps");
+
+                psTrace("ppSub", 3, "Calculating equation...\n");
+                if (!pmSubtractionCalculateEquation(stamps, refRO->image, inRO->image, inRO->weight,
+                                                    kernels, footprint)) {
+                    psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation.");
+                    goto ERROR;
+                }
+
+                memCheck("  calculate equation");
+
+                psTrace("ppSub", 3, "Solving equation...\n");
+                solution = pmSubtractionSolveEquation(solution, stamps);
+                if (!solution) {
+                    psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation.");
+                    goto ERROR;
+                }
+
+                memCheck("  solve equation");
+
+                psTrace("ppSub", 3, "Rejecting stamps...\n");
+                numRejected = pmSubtractionRejectStamps(stamps, refRO->image, inRO->image, subMask,
+                                                        solution, footprint, rej, kernels);
+                if (numRejected < 0) {
+                    psError(PS_ERR_UNKNOWN, false, "Unable to reject stamps.");
+                    goto ERROR;
+                }
+                psLogMsg("ppSub", PS_LOG_INFO, "%d stamps rejected on iteration %d.", numRejected, k);
+
+                memCheck("  reject stamps");
+            }
+
+            if (numRejected > 0) {
+                psTrace("ppSub", 3, "Solving equation...\n");
+                solution = pmSubtractionSolveEquation(solution, stamps);
+                if (!solution) {
+                    psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation.");
+                    goto ERROR;
+                }
+            }
+            psFree(stamps);
+            stamps = NULL;
+
+            memCheck("solution");
 
 #ifdef TESTING
-    {
-        // Generate image with convolution kernels
-        int fullSize = 2 * size + 1 + 1;    // Full size of kernel
-        psImage *convKernels = psImageAlloc(5 * fullSize - 1, 5 * fullSize - 1, PS_TYPE_F32);
-        psImageInit(convKernels, NAN);
-        for (int j = -2; j <= 2; j++) {
-            for (int i = -2; i <= 2; i++) {
-                psImage *kernel = pmSubtractionKernelImage(solution, kernels, (float)i / 2.0,
-                                                           (float)j / 2.0); // Image of the kernel
-                if (!kernel) {
-                    psError(PS_ERR_UNKNOWN, false, "Unable to generate kernel image.");
-                    psFree(convKernels);
-                    goto ERROR;
-                }
-
-                if (psImageOverlaySection(convKernels, kernel, (i + 2) * fullSize,
-                                          (j + 2) * fullSize, "=") == 0) {
-                    psError(PS_ERR_UNKNOWN, false, "Unable to overlay kernel image.");
-                    psFree(kernel);
-                    psFree(convKernels);
-                    goto ERROR;
-                }
-                psFree(kernel);
-            }
-        }
-
-        // XXX What do we do with this image?
-
-        psFits *kernelFile = psFitsOpen("kernel.fits", "w");
-        (void)psFitsWriteImage(kernelFile, NULL, convKernels, 0, NULL);
-        psFitsClose(kernelFile);
-
-        psFree(convKernels);
-    }
-
-    {
-        // Generate images of the kernel components
-        psMetadata *header = psMetadataAlloc(); // Header
-        for (int i = 0; i < solution->n; i++) {
-            psString name = NULL;       // Header keyword
-            psStringAppend(&name, "SOLN%04d", i);
-            psMetadataAddF64(header, PS_LIST_TAIL, name, 0, NULL, solution->data.F64[i]);
-            psFree(name);
-        }
-        psArray *kernelImages = pmSubtractionKernelSolutions(solution, kernels, 0.0, 0.0);
-        psFits *kernelFile = psFitsOpen("kernels.fits", "w");
-        (void)psFitsWriteImageCube(kernelFile, header, kernelImages, NULL);
-        psFitsClose(kernelFile);
-        psFree(kernelImages);
-        psFree(header);
-    }
+            psTrace("ppSub", 2, "Generating diagnostics...\n");
+            {
+                // Generate image with convolution kernels
+                int fullSize = 2 * size + 1 + 1;    // Full size of kernel
+                psImage *convKernels = psImageAlloc(5 * fullSize - 1, 5 * fullSize - 1, PS_TYPE_F32);
+                psImageInit(convKernels, NAN);
+                for (int j = -2; j <= 2; j++) {
+                    for (int i = -2; i <= 2; i++) {
+                        psImage *kernel = pmSubtractionKernelImage(solution, kernels, (float)i / 2.0,
+                                                                   (float)j / 2.0); // Image of the kernel
+                        if (!kernel) {
+                            psError(PS_ERR_UNKNOWN, false, "Unable to generate kernel image.");
+                            psFree(convKernels);
+                            goto ERROR;
+                        }
+
+                        if (psImageOverlaySection(convKernels, kernel, (i + 2) * fullSize,
+                                                  (j + 2) * fullSize, "=") == 0) {
+                            psError(PS_ERR_UNKNOWN, false, "Unable to overlay kernel image.");
+                            psFree(kernel);
+                            psFree(convKernels);
+                            goto ERROR;
+                        }
+                        psFree(kernel);
+                    }
+                }
+
+                // XXX What do we do with this image?
+
+                psFits *kernelFile = psFitsOpen("kernel.fits", "w");
+                (void)psFitsWriteImage(kernelFile, NULL, convKernels, 0, NULL);
+                psFitsClose(kernelFile);
+
+                psFree(convKernels);
+            }
+
+            {
+                // Generate images of the kernel components
+                psMetadata *header = psMetadataAlloc(); // Header
+                for (int i = 0; i < solution->n; i++) {
+                    psString name = NULL;       // Header keyword
+                    psStringAppend(&name, "SOLN%04d", i);
+                    psMetadataAddF64(header, PS_LIST_TAIL, name, 0, NULL, solution->data.F64[i]);
+                    psFree(name);
+                }
+                psArray *kernelImages = pmSubtractionKernelSolutions(solution, kernels, 0.0, 0.0);
+                psFits *kernelFile = psFitsOpen("kernels.fits", "w");
+                (void)psFitsWriteImageCube(kernelFile, header, kernelImages, NULL);
+                psFitsClose(kernelFile);
+                psFree(kernelImages);
+                psFree(header);
+            }
 #endif
 
-    memCheck("diag outputs");
-
-    psImage *convImage = NULL, *convWeight = NULL, *convMask = NULL; // Convolved images
-    if (!pmSubtractionConvolve(&convImage, &convWeight, &convMask, refRO->image, refRO->weight, subMask,
-                               maskBlank, solution, kernels)) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to convolve reference image.");
-        goto ERROR;
-    }
+            memCheck("diag outputs");
+
+            psTrace("ppSub", 2, "Convolving...\n");
+            if (!pmSubtractionConvolve(&convImage, &convWeight, &convMask, refRO->image, refRO->weight,
+                                       subMask, maskBlank, region, solution, kernels)) {
+                psError(PS_ERR_UNKNOWN, false, "Unable to convolve reference image.");
+                goto ERROR;
+            }
+
+            psFree(solution);
+            solution = NULL;
+        }
+    }
+
+    psFree(region);
+    region = NULL;
     psFree(subMask);
     subMask = NULL;
 
+    if (!pmSubtractionBorder(convImage, convWeight, convMask, kernels, maskBlank)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to set border region of convolved image.");
+        goto ERROR;
+    }
+
     psFree(kernels);
-    psFree(solution);
 
     memCheck("convolution");
@@ -265 +311 @@
         for (int x = 0; x < outRO->image->numCols; x++) {
             if (isnan(outRO->image->data.F32[y][x]) && !(outRO->mask->data.U8[y][x] & maskBlank)) {
-                printf("%d %d --> %d\n", x, y, outRO->mask->data.U8[y][x]);
+                printf("Unmasked NAN at %d %d --> %d\n", x, y, outRO->mask->data.U8[y][x]);
             }
         }
@@ -315 +361 @@
 
 ERROR:
+    psFree(region);
     psFree(subMask);
     psFree(kernels);