Changeset 23351 for branches/cnb_branches/cnb_branch_20090301/psModules/src/imcombine/pmSubtractionMatch.c

Timestamp:

Mar 17, 2009, 12:07:42 PM (17 years ago)

Author:

beaumont

Message:

merged with head

Location:

branches/cnb_branches/cnb_branch_20090301/psModules

Files:

• . (modified) (1 prop)
• src/imcombine/pmSubtractionMatch.c (modified) (4 diffs)

branches/cnb_branches/cnb_branch_20090301/psModules/src/imcombine/pmSubtractionMatch.c

@@ -90 +90 @@
 }
 
+// Check input arguments
+static bool subtractionMatchCheck(pmReadout *conv1, pmReadout *conv2, // Convolved images
+                                  const pmReadout *ro1, const pmReadout *ro2, // Input images
+                                  int stride, // Size for convolution patches
+                                  float sysError, // Relative systematic error
+                                  psImageMaskType maskVal, // Value to mask for input
+                                  psImageMaskType maskBad, // Mask for output bad pixels
+                                  psImageMaskType maskPoor, // Mask for output poor pixels
+                                  float poorFrac, // Fraction for "poor"
+                                  float badFrac,   // Maximum fraction of bad input pixels to accept
+                                  pmSubtractionMode subMode // Mode of subtraction
+                                  )
+{
+    if (subMode != PM_SUBTRACTION_MODE_2) {
+        PM_ASSERT_READOUT_NON_NULL(conv1, false);
+        if (conv1->image) {
+            psFree(conv1->image);
+            conv1->image = NULL;
+        }
+        if (conv1->mask) {
+            psFree(conv1->mask);
+            conv1->mask = NULL;
+        }
+        if (conv1->variance) {
+            psFree(conv1->variance);
+            conv1->variance = NULL;
+        }
+    }
+    if (subMode != PM_SUBTRACTION_MODE_1) {
+        PM_ASSERT_READOUT_NON_NULL(conv2, false);
+        if (conv2->image) {
+            psFree(conv2->image);
+            conv2->image = NULL;
+        }
+        if (conv2->mask) {
+            psFree(conv2->mask);
+            conv2->mask = NULL;
+        }
+        if (conv2->variance) {
+            psFree(conv2->variance);
+            conv2->variance = NULL;
+        }
+    }
+
+    PM_ASSERT_READOUT_NON_NULL(ro1, false);
+    PM_ASSERT_READOUT_NON_NULL(ro2, false);
+    PM_ASSERT_READOUT_IMAGE(ro1, false);
+    PM_ASSERT_READOUT_IMAGE(ro2, false);
+    PS_ASSERT_IMAGES_SIZE_EQUAL(ro1->image, ro2->image, false);
+    PS_ASSERT_INT_NONNEGATIVE(stride, false);
+    if (isfinite(sysError)) {
+        PS_ASSERT_FLOAT_LARGER_THAN_OR_EQUAL(sysError, 0.0, false);
+        PS_ASSERT_FLOAT_LESS_THAN(sysError, 1.0, false);
+    }
+    // Don't care about maskVal
+    // Don't care about maskBad
+    // Don't care about maskPoor
+    PS_ASSERT_FLOAT_LARGER_THAN(poorFrac, 0.0, false);
+    PS_ASSERT_FLOAT_LESS_THAN_OR_EQUAL(poorFrac, 1.0, false);
+    if (isfinite(badFrac)) {
+        PS_ASSERT_FLOAT_LARGER_THAN(badFrac, 0.0, false);
+        PS_ASSERT_FLOAT_LESS_THAN_OR_EQUAL(badFrac, 1.0, false);
+    }
+
+    return true;
+}
+
+
+
+bool pmSubtractionMatchPrecalc(pmReadout *conv1, pmReadout *conv2, const pmReadout *ro1, const pmReadout *ro2,
+                               psMetadata *analysis, int stride, float sysError,
+                               psImageMaskType maskVal, psImageMaskType maskBad, psImageMaskType maskPoor,
+                               float poorFrac, float badFrac)
+{
+    PS_ASSERT_METADATA_NON_NULL(analysis, false);
+
+    // Extract the kernels
+    pmSubtractionMode mode = PM_SUBTRACTION_MODE_UNSURE; // Subtraction mode: which image to convolve
+    int size = 0;                       // Size of kernel
+    psList *kernelList = psListAlloc(NULL); // List of kernels
+    {
+        psMetadataIterator *iter = psMetadataIteratorAlloc(analysis, PS_LIST_HEAD,
+                                                           "^" PM_SUBTRACTION_ANALYSIS_KERNEL "$");
+        psMetadataItem *item;               // Item from iteration
+        while ((item = psMetadataGetAndIncrement(iter))) {
+            if (item->type != PS_DATA_UNKNOWN) {
+                psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Unexpected type for kernel.");
+                psFree(iter);
+                psFree(kernelList);
+                return false;
+            }
+            pmSubtractionKernels *kernel = item->data.V; // Kernel
+            PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernel, false);
+            size = PS_MAX(size, kernel->size);
+            if (mode == PM_SUBTRACTION_MODE_UNSURE) {
+                mode = kernel->mode;
+            } else if (kernel->mode != mode) {
+                // There's some confusion, so let's set the mode to dual convolution.
+                // This is only used for the subtraction mask, so it's not a big deal.
+                mode = PM_SUBTRACTION_MODE_DUAL;
+            }
+            psListAdd(kernelList, PS_LIST_TAIL, kernel);
+        }
+        psFree(iter);
+    }
+    if (psListLength(kernelList) == 0) {
+        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Unable to find kernels");
+        psFree(kernelList);
+        return false;
+    }
+    psArray *kernels = psListToArray(kernelList); // Array of kernels
+    psFree(kernelList);
+
+    // Extract the regions
+    psArray *regions = psArrayAllocEmpty(kernels->n); // Array of regions
+    {
+        psMetadataIterator *iter = psMetadataIteratorAlloc(analysis, PS_LIST_HEAD,
+                                                           "^" PM_SUBTRACTION_ANALYSIS_REGION "$");
+        psMetadataItem *item;               // Item from iteration
+        while ((item = psMetadataGetAndIncrement(iter))) {
+            if (item->type != PS_DATA_REGION) {
+                psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Unexpected type for region.");
+                psFree(iter);
+                psFree(kernels);
+                psFree(regions);
+                return false;
+            }
+            psRegion *region = item->data.V; // Region
+            psArrayAdd(regions, regions->n, region);
+        }
+        psFree(iter);
+    }
+    if (regions->n != kernels->n) {
+        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Differing number of kernels (%ld) and regions (%ld)",
+                kernels->n, regions->n);
+        psFree(regions);
+        psFree(kernels);
+        return false;
+    }
+
+    if (!subtractionMatchCheck(conv1, conv2, ro1, ro2, stride, sysError, maskVal, maskBad, maskPoor,
+                               poorFrac, badFrac, mode)) {
+        psFree(kernels);
+        psFree(regions);
+        return false;
+    }
+
+    psImage *subMask = pmSubtractionMask(ro1->mask, ro2 ? ro2->mask : NULL, maskVal, size, 0,
+                                         badFrac, useFFT); // Subtraction mask
+    if (!subMask) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to generate subtraction mask.");
+        psFree(kernels);
+        psFree(regions);
+        return false;
+    }
+
+    psMetadata *outAnalysis = psMetadataAlloc(); // Output analysis values
+
+    psTrace("psModules.imcombine", 2, "Convolving...\n");
+    for (int i = 0; i < kernels->n; i++) {
+        pmSubtractionKernels *kernel = kernels->data[i]; // Kernel of interest
+        psRegion *region = regions->data[i]; // Region of interest
+
+        if (!pmSubtractionAnalysis(outAnalysis, kernel, region, ro1->image->numCols, ro1->image->numRows)) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to generate QA data");
+            psFree(outAnalysis);
+            psFree(subMask);
+            psFree(kernels);
+            psFree(regions);
+            return false;
+        }
+
+        if (!pmSubtractionConvolve(conv1, conv2, ro1, ro2, subMask, stride, maskBad, maskPoor, poorFrac,
+                                   sysError, region, kernel, true, useFFT)) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to convolve image.");
+            psFree(outAnalysis);
+            psFree(subMask);
+            psFree(kernels);
+            psFree(regions);
+            return false;
+        }
+    }
+
+    psFree(subMask);
+    psFree(kernels);
+    psFree(regions);
+
+    if (conv1) {
+        psMetadataCopy(conv1->analysis, outAnalysis);
+    }
+    if (conv2) {
+        psMetadataCopy(conv2->analysis, outAnalysis);
+    }
+    psFree(outAnalysis);
+
+    return true;
+}
+
 
 bool pmSubtractionMatch(pmReadout *conv1, pmReadout *conv2, const pmReadout *ro1, const pmReadout *ro2,
@@ -101 +299 @@
                         psImageMaskType maskPoor, float poorFrac, float badFrac, pmSubtractionMode subMode)
 {
-    if (subMode != PM_SUBTRACTION_MODE_2) {
-        PM_ASSERT_READOUT_NON_NULL(conv1, false);
-        if (conv1->image) {
-            psFree(conv1->image);
-            conv1->image = NULL;
-        }
-        if (conv1->mask) {
-            psFree(conv1->mask);
-            conv1->mask = NULL;
-        }
-        if (conv1->variance) {
-            psFree(conv1->variance);
-            conv1->variance = NULL;
-        }
-    }
-    if (subMode != PM_SUBTRACTION_MODE_1) {
-        PM_ASSERT_READOUT_NON_NULL(conv2, false);
-        if (conv2->image) {
-            psFree(conv2->image);
-            conv2->image = NULL;
-        }
-        if (conv2->mask) {
-            psFree(conv2->mask);
-            conv2->mask = NULL;
-        }
-        if (conv2->variance) {
-            psFree(conv2->variance);
-            conv2->variance = NULL;
-        }
-    }
-
-    PM_ASSERT_READOUT_NON_NULL(ro1, false);
-    PM_ASSERT_READOUT_NON_NULL(ro2, false);
-    PM_ASSERT_READOUT_IMAGE(ro1, false);
-    PM_ASSERT_READOUT_IMAGE(ro2, false);
-    PS_ASSERT_IMAGES_SIZE_EQUAL(ro1->image, ro2->image, false);
+    if (!subtractionMatchCheck(conv1, conv2, ro1, ro2, stride, sysError, maskVal, maskBad, maskPoor,
+                               poorFrac, badFrac, subMode)) {
+        return false;
+    }
 
     PS_ASSERT_INT_NONNEGATIVE(footprint, false);
-    PS_ASSERT_INT_NONNEGATIVE(stride, false);
     // regionSize can be just about anything (except maybe negative, but it can be NAN)
     PS_ASSERT_FLOAT_LARGER_THAN(stampSpacing, 0.0, false);
@@ -168 +333 @@
     PS_ASSERT_INT_POSITIVE(iter, false);
     PS_ASSERT_FLOAT_LARGER_THAN(rej, 0.0, false);
-    if (isfinite(sysError)) {
-        PS_ASSERT_FLOAT_LARGER_THAN_OR_EQUAL(sysError, 0.0, NULL);
-        PS_ASSERT_FLOAT_LESS_THAN(sysError, 1.0, NULL);
-    }
-    // Don't care about maskVal
-    // Don't care about maskBad
-    // Don't care about maskPoor
-    PS_ASSERT_FLOAT_LARGER_THAN(poorFrac, 0.0, NULL);
-    PS_ASSERT_FLOAT_LESS_THAN_OR_EQUAL(poorFrac, 1.0, NULL);
-    if (isfinite(badFrac)) {
-        PS_ASSERT_FLOAT_LARGER_THAN(badFrac, 0.0, NULL);
-        PS_ASSERT_FLOAT_LESS_THAN_OR_EQUAL(badFrac, 1.0, NULL);
-    }
 
     // If the stamp footprint is smaller than the kernel size, then we won't get much signal in the outer
@@ -214 +366 @@
     int numCols = ro1->image->numCols, numRows = ro1->image->numRows; // Image dimensions
 
-    psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS, 0); // Random number generator
+    psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); // Random number generator
 
     memCheck("start");