Changeset 14586

Timestamp:

Aug 21, 2007, 10:37:05 AM (19 years ago)

Author:

Paul Price

Message:

Moved all the work into psModules function pmSubtractionMatch, to
avoid code duplication (with ppStack, which also wants to do image
matching).

File:

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

trunk/ppSub/src/ppSubReadout.c

@@ -10 +10 @@
 #include "ppSub.h"
 
-//#define TESTING
 #define WCS_TOLERANCE 0.001             // Tolerance for WCS
-
-#include <unistd.h>
-static void memCheck(const char *where)
-{
-    return;
-    psMemBlock **leaks = NULL;
-    int numLeaks = psMemCheckLeaks(0, &leaks, NULL, true);
-    size_t largestSize = 0;
-    psMemId largest = 0;
-    size_t totalSize = 0;
-    for (int i = 0; i < numLeaks; i++) {
-        psMemBlock *mb = leaks[i];
-        totalSize += mb->userMemorySize;
-        if (mb->userMemorySize > largestSize) {
-            largestSize = mb->userMemorySize;
-            largest = mb->id;
-        }
-    }
-    psFree(leaks);
-    fprintf(stderr, "%s:\n", where);
-    fprintf(stderr, "    Memory in use: %zd\n", totalSize);
-    fprintf(stderr, "    Largest block: %ld\n", largest);
-    fprintf(stderr, "    sbrk(): %zd\n", (size_t)sbrk(0));
-    return;
-}
 
 bool ppSubReadout(pmConfig *config, const pmFPAview *view)
@@ -75 +49 @@
     psMaskType maskBlank = pmConfigMask(psMetadataLookupStr(NULL, config->arguments, "MASK.BLANK"),
                                         config); // Mask for blank reg.
-    const char *stampsFile = psMetadataLookupStr(NULL, config->arguments, "STAMPS"); // Filename for stamps
+    const char *stampsName = psMetadataLookupStr(NULL, config->arguments, "STAMPS"); // Filename for stamps
 
+    // Generate masks if they don't exist
     int numCols = input->numCols, numRows = input->numRows; // Image dimensions
     if (!inRO->mask) {
@@ -95 +70 @@
     }
 
-    memCheck("start");
-
-    // Mask for subtraction
-    psImage *subMask = pmSubtractionMask(refRO->mask, inRO->mask, maskBad, size, footprint);
-
-    memCheck("mask");
-
-#if 0
-    if (!inRO->weight) {
-        // Need to ensure input image is completely above zero, so we're not weighting by negative numbers
-        psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
-        if (!psImageBackground(stats, inRO->image, inRO->mask, maskBad, NULL)) {
-            psError(PS_ERR_UNKNOWN, false, "Unable to determine statistics of input image.");
-            return false;
-        }
-        float diff = stats->robustMedian - 5.0 * stats->robustStdev;
-        if (diff < 0.0) {
-            psLogMsg("ppSub", PS_LOG_INFO, "Adding %f to image background to get positive pixels.");
-            (void)psBinaryOp(in->image, inRO->image, "+", psScalarAlloc(diff, PS_TYPE_F32));
-        }
-    }
-#endif
-
-    // Kernel basis functions
-    pmSubtractionKernels *kernels = pmSubtractionKernelsGenerate(type, size, order, widths, orders,
-                                                                 inner, binning, ringsOrder);
-    psArray *stamps = NULL;             // Stamps for matching PSF
-    psVector *solution = NULL;          // Solution to match PSF
-
-    memCheck("kernels");
-
-    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) || regionSize == 0.0) {
-        xRegions = numCols / regionSize + 1;
-        yRegions = numRows / regionSize + 1;
-        xRegionSize = (float)numCols / (float)xRegions;
-        yRegionSize = (float)numRows / (float)yRegions;
-        region = psRegionAlloc(NAN, NAN, NAN, NAN);
+    if (!pmSubtractionMatch(outRO, refRO, inRO, footprint, regionSize, spacing, threshold, stampsName,
+                            NAN, type, size, order, widths, orders, inner, ringsOrder, binning, iter,
+                            rej, maskBad, maskBlank)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to match images.");
+        psFree(outRO);
+        return false;
     }
 
-    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");
-                if (stampsFile) {
-                    iter = 1;           // There is no iterating because we use all the stamps we have
-                    psArray *stampData = psVectorsReadFromFile(stampsFile, "%f %f"); // Stamp positions
-                    psVector *xStamp = stampData->data[0]; // x coordinates
-                    psVector *yStamp = stampData->data[1]; // y coordinates
-                    // Correct for IRAF (unit-offset) positions to C (zero-offset) positions
-                    psBinaryOp(xStamp, xStamp, "-", psScalarAlloc(1.0, PS_TYPE_F32));
-                    psBinaryOp(yStamp, yStamp, "-", psScalarAlloc(1.0, PS_TYPE_F32));
-
-                    stamps = pmSubtractionSetStamps(xStamp, yStamp, NULL, subMask, region);
-                } else {
-                    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, "Extracting stamps...\n");
-                if (!pmSubtractionExtractStamps(stamps, refRO->image, inRO->image, inRO->weight,
-                                                footprint, kernels)) {
-                    psError(PS_ERR_UNKNOWN, false, "Unable to extract stamps.");
-                    goto ERROR;
-                }
-
-                memCheck("   extract stamps");
-
-                psTrace("ppSub", 3, "Calculating equation...\n");
-                if (!pmSubtractionCalculateEquation(stamps, 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, 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
-            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");
-
-            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;
-    }
 
     // Add kernel descrption to header
-    psMetadataAddStr(outHDU->header, PS_LIST_TAIL, "PPSUB.KERNEL", 0,
-                     "Subtraction kernel", kernels->description);
-    psFree(kernels);
-
-    memCheck("convolution");
+    psMetadataAddStr(outHDU->header, PS_LIST_TAIL, "PPSUB.KERNEL", 0, "Subtraction kernel",
+                     psMetadataLookupStr(NULL, outRO->analysis, "SUBTRACTION.KERNEL.DESCRIPTION"));
 
     // Do the subtraction
     if (reverse) {
-        outRO->image = (psImage*)psBinaryOp(NULL, convImage, "-", inRO->image);
+        psBinaryOp(outRO->image, outRO->image, "-", inRO->image);
     } else {
-        outRO->image = (psImage*)psBinaryOp(NULL, inRO->image, "-", convImage);
+        psBinaryOp(outRO->image, inRO->image, "-", outRO->image);
     }
-    outRO->mask = (psImage*)psBinaryOp(NULL, convMask, "|", inRO->mask);
-    if (convWeight) {
-        outRO->weight = (psImage*)psBinaryOp(NULL, convWeight, "+", inRO->weight);
+    psBinaryOp(outRO->mask, outRO->mask, "|", inRO->mask);
+    if (outRO->weight) {
+        psBinaryOp(outRO->weight, outRO->weight, "+", inRO->weight);
     }
-
-    memCheck("subtraction");
-
-    psFree(convImage);
-    psFree(convMask);
-    psFree(convWeight);
-
-    memCheck("subtraction freed");
-
-    outRO->data_exists = true;
-    outCell->data_exists = true;
-    outCell->parent->data_exists = true;
 
 #ifdef TESTING
@@ -349 +106 @@
     if (!pmFPACopyConcepts(outCell->parent->parent, inRO->parent->parent->parent)) {
         psError(PS_ERR_UNKNOWN, false, "Unable to copy concepts from input to output.");
+        psFree(outRO);
         return false;
     }
@@ -387 +145 @@
                      "Subtraction input", inFile->filename);
 
-    memCheck("header");
-
     psFree(outRO);
 
     return true;
-
-ERROR:
-    psFree(region);
-    psFree(subMask);
-    psFree(kernels);
-    psFree(stamps);
-    psFree(solution);
-    psFree(outRO);
-    return false;
 }