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

Timestamp:

Jul 10, 2007, 2:15:40 PM (19 years ago)

Author:

Paul Price

Message:

Supporting new kernel types.

File:

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

trunk/ppSub/src/ppSubReadout.c

@@ -6 +6 @@
 #include <pslib.h>
 #include <psmodules.h>
-#include <psphot.h>
 
 #include "ppSub.h"
+
+#define MASK_BAD      0x01              // Mask value for bad pixel
+#define MASK_STAMP    0x02              // Mask value for bad stamp (and bad stamp region)
 
 bool ppSubReadout(pmConfig *config, const pmFPAview *view)
@@ -14 +16 @@
     pmReadout *inRO = pmFPAfileThisReadout(config->files, view, "PPSUB.INPUT"); // Input readout
     pmReadout *refRO = pmFPAfileThisReadout(config->files, view, "PPSUB.REF"); // Reference readout
+#if 0
     pmCell *outCell = pmFPAfileThisCell(config->files, view, "PPSUB.OUTPUT"); // Output cell
     pmReadout *outRO = pmReadoutAlloc(outCell); // Output readout
+#endif
 
     psImage *input = inRO->image;       // Input image
@@ -34 +38 @@
     psVector *widths = psMetadataLookupPtr(NULL, config->arguments, "ISIS.WIDTHS"); // ISIS Gaussian widths
     psVector *orders = psMetadataLookupPtr(NULL, config->arguments, "ISIS.ORDERS"); // ISIS Polynomial orders
+    int inner = psMetadataLookupS32(NULL, config->arguments, "INNER"); // Inner radius
+    int binning = psMetadataLookupS32(NULL, config->arguments, "SPAM.BINNING"); // Binning for SPAM kernel
     psMaskType maskBad = psMetadataLookupU8(NULL, config->arguments, "MASK.BAD"); // Value to mask
     psMaskType maskBlank = psMetadataLookupU8(NULL, config->arguments, "MASK.BLANK"); // Mask for blank reg.
 
-    int numCols = input->numCols, numRows = input->numRows; // Image dimensions
-    if (!inRO->mask) {
-        inRO->mask = psImageAlloc(numCols, numRows, PS_TYPE_MASK);
-        psImageInit(inRO->mask, 0);
+    if (!inRO->mask && !pmReadoutGenerateMask(inRO)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to generate mask for input image");
+        return false;
     }
-    if (!refRO->mask) {
-        refRO->mask = psImageAlloc(numCols, numRows, PS_TYPE_MASK);
-        psImageInit(refRO->mask, 0);
+    if (!inRO->weight && !pmReadoutGenerateWeight(inRO, true)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to generate weight map for input image");
+        return false;
+    }
+    if (!refRO->mask && !pmReadoutGenerateMask(refRO)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to generate mask for reference image");
+        return false;
+    }
+    if (!refRO->weight && !pmReadoutGenerateWeight(refRO, true)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to generate weight map for reference image");
+        return false;
     }
 
-    // Mask for subtraction
-    psImage *subMask = pmSubtractionMask(inRO->mask, refRO->mask, maskBad, size, footprint);
+    // Worried about the masks for bad pixels and bad stamps colliding, so make our own mask
+    int numCols = input->numCols, numRows = input->numRows; // Image dimensions
+    psImage *stampMask = psImageAlloc(numCols, numRows, PS_TYPE_MASK); // Mask to use for stamps
+    for (int y = 0; y < numRows; y++) {
+        for (int x = 0; x < numCols; x++) {
+            stampMask->data.PS_TYPE_MASK_DATA[y][x] =
+                (refRO->mask->data.PS_TYPE_MASK_DATA[y][x] & maskBad) ? MASK_BAD : 0;
+        }
+    }
 
 #if 0
@@ -66 +86 @@
 #endif
 
-    pmSubtractionKernels *kernels = pmSubtractionKernelsGenerate(type, size, order,
-                                                                 widths, orders); // Kernel basis functions
+    pmSubtractionKernels *kernels = pmSubtractionKernelsGenerate(type, size, order, widths, orders,
+                                                                 inner, binning); // Kernel basis functions
     psArray *stamps = NULL;             // Stamps for matching PSF
     psVector *solution = NULL;          // Solution to match PSF
@@ -73 +93 @@
     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);
+        stamps = pmSubtractionFindStamps(stamps, refRO->image, stampMask, MASK_BAD, MASK_STAMP,
+                                         threshold, spacing, size + footprint);
         if (!stamps) {
             psError(PS_ERR_UNKNOWN, false, "Unable to find stamps on reference image.");
@@ -91 +112 @@
         }
 
-        numRejected = pmSubtractionRejectStamps(stamps, refRO->image, inRO->image, subMask,
+        numRejected = pmSubtractionRejectStamps(stamps, refRO->image, inRO->image, stampMask, MASK_STAMP,
                                                 solution, footprint, rej, kernels);
         if (numRejected < 0) {
@@ -99 +120 @@
         psLogMsg("ppSub", PS_LOG_INFO, "%d stamps rejected on iteration %d.", numRejected, i);
     }
-
-    if (numRejected > 0) {
-        solution = pmSubtractionSolveEquation(solution, stamps);
-        if (!solution) {
-            psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation.");
-            goto ERROR;
-        }
-    }
+    psFree(stampMask);
 
     psImage *convImage = NULL, *convWeight = NULL, *convMask = NULL; // Convolved images
-    if (!pmSubtractionConvolve(&convImage, &convWeight, &convMask, refRO->image, refRO->weight, subMask,
-                               maskBlank, solution, kernels)) {
+    if (!pmSubtractionConvolve(&convImage, &convWeight, &convMask,
+                               refRO->image, refRO->weight, refRO->mask,
+                               MASK_BAD, maskBlank, solution, kernels)) {
         psError(PS_ERR_UNKNOWN, false, "Unable to convolve reference image.");
         goto ERROR;
     }
-    psFree(subMask);
 
     // Do the subtraction
     if (reverse) {
-        outRO->image = (psImage*)psBinaryOp(NULL, convImage, "-", inRO->image);
+        (void)psBinaryOp(inRO->image, convImage, "-", inRO->image);
     } else {
-        outRO->image = (psImage*)psBinaryOp(NULL, inRO->image, "-", convImage);
+        (void)psBinaryOp(inRO->image, inRO->image, "-", convImage);
     }
-    outRO->mask = (psImage*)psBinaryOp(NULL, convMask, "|", inRO->mask);
-    if (convWeight) {
-        outRO->weight = (psImage*)psBinaryOp(NULL, convWeight, "+", inRO->weight);
-    }
+    (void)psBinaryOp(inRO->mask, convMask, "|", inRO->mask);
+    (void)psBinaryOp(inRO->weight, convWeight, "+", inRO->weight);
 
     psFree(convImage);
     psFree(convMask);
     psFree(convWeight);
-
-    outRO->data_exists = true;
-    outCell->data_exists = true;
-    outCell->parent->data_exists = true;
-
-    if (!pmFPACopyConcepts(outCell->parent->parent, inRO->parent->parent->parent)) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to copy concepts from input to output.");
-        return false;
-    }
 
 #if 0
@@ -173 +176 @@
 #endif
 
-    psFree(kernels);
-    psFree(stamps);
-    psFree(solution);
 
-    if (psMetadataLookupBool(NULL, config->arguments, "PHOTOMETRY")) {
-        pmFPAfile *photFile = psMetadataLookupPtr(NULL, config->files, "PSPHOT.INPUT");
-        pmFPACopy(photFile->fpa, inRO->parent->parent->parent);
+        psFree(kernels);
+        psFree(stamps);
+        psFree(solution);
+        return true;
 
-        psphotReadout(config, view);
-
-        pmFPAfileActivate(config->files, false, "PSPHOT.INPUT");
-    }
-
-    return true;
-
-ERROR:
-    psFree(kernels);
-    psFree(stamps);
-    psFree(solution);
-    return false;
+    ERROR:
+        psFree(kernels);
+        psFree(stamps);
+        psFree(solution);
+        return false;
 }