Changeset 27365 for trunk/psModules/src/imcombine/pmSubtractionMatch.c

Timestamp:

Mar 19, 2010, 5:09:20 PM (16 years ago)

Author:

Paul Price

Message:

Adding program ppSubConvolve to convolve an image by the (previously calculated) convolution kernel. Made a small change to the API for initialising PSF-matching (pmSubtraction) threads. pmSubtractionMatchPrecalc now does a quick and dirty background subtraction, just like pmSubtractionMatch does. Changed a few functions to allow one of the input readouts to be NULL.

File:

• trunk/psModules/src/imcombine/pmSubtractionMatch.c (modified) (7 diffs)

trunk/psModules/src/imcombine/pmSubtractionMatch.c

@@ -132 +132 @@
     if (subMode != PM_SUBTRACTION_MODE_2) {
         PM_ASSERT_READOUT_NON_NULL(conv1, false);
+        PM_ASSERT_READOUT_NON_NULL(ro1, false);
+        PM_ASSERT_READOUT_IMAGE(ro1, false);
         if (conv1->image) {
             psFree(conv1->image);
@@ -147 +149 @@
     if (subMode != PM_SUBTRACTION_MODE_1) {
         PM_ASSERT_READOUT_NON_NULL(conv2, false);
+        PM_ASSERT_READOUT_NON_NULL(ro2, false);
+        PM_ASSERT_READOUT_IMAGE(ro2, false);
         if (conv2->image) {
             psFree(conv2->image);
@@ -161 +165 @@
     }
 
-    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 (ro1 && ro2) {
+        PS_ASSERT_IMAGES_SIZE_EQUAL(ro1->image, ro2->image, false);
+    }
     PS_ASSERT_INT_NONNEGATIVE(stride, false);
     if (isfinite(normFrac)) {
@@ -379 +381 @@
     }
 
+    int numCols, numRows;       // Size of image
+    if (ro1) {
+        numCols = ro1->image->numCols;
+        numRows = ro1->image->numRows;
+    } else if (ro2) {
+        numCols = ro2->image->numCols;
+        numRows = ro2->image->numRows;
+    } else {
+        psAbort("No input image provided.");
+    }
+
     pmSubtractionMaskInvalid(ro1, maskVal);
     pmSubtractionMaskInvalid(ro2, maskVal);
+
+    // General background subtraction, since this is done in pmSubtractionMatch
+    {
+        psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); // Random number generator
+        psStats *bg = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); // Statistics for background
+        if (ro1) {
+            psStatsInit(bg);
+            if (!psImageBackground(bg, NULL, ro1->image, ro1->mask, maskVal, rng)) {
+                psError(PM_ERR_DATA, false, "Unable to measure background statistics.");
+                psFree(bg);
+                psFree(rng);
+                return false;
+            }
+            psBinaryOp(ro1->image, ro1->image, "-", psScalarAlloc((float)bg->robustMedian, PS_TYPE_F32));
+        }
+        if (ro2) {
+            psStatsInit(bg);
+            if (!psImageBackground(bg, NULL, ro2->image, ro2->mask, maskVal, rng)) {
+                psError(PM_ERR_DATA, false, "Unable to measure background statistics.");
+                psFree(bg);
+                psFree(rng);
+                return false;
+            }
+            psBinaryOp(ro2->image, ro2->image, "-", psScalarAlloc((float)bg->robustMedian, PS_TYPE_F32));
+        }
+        psFree(bg);
+        psFree(rng);
+    }
 
     psRegion bounds = psRegionSet(NAN, NAN, NAN, NAN); // Bounds of valid pixels
@@ -396 +437 @@
     psMetadata *outHeader = psMetadataAlloc(); // Output header values
 
-    int numCols = ro1->image->numCols, numRows = ro1->image->numRows; // Size of image
     subtractionMatchAlloc(conv1, conv2, ro1, ro2, subMask, maskBad, mode, numCols, numRows);
 
@@ -404 +444 @@
         psRegion *region = regions->data[i]; // Region of interest
 
-        if (!pmSubtractionAnalysis(outAnalysis, outHeader, kernel, region,
-                                   ro1->image->numCols, ro1->image->numRows)) {
+        if (!pmSubtractionAnalysis(outAnalysis, outHeader, kernel, region, numCols, numRows)) {
             psError(psErrorCodeLast(), false, "Unable to generate QA data");
             psFree(outAnalysis);
@@ -454 +493 @@
         return false;
     }
+
+    // We need both inputs
+    PM_ASSERT_READOUT_NON_NULL(ro1, false);
+    PM_ASSERT_READOUT_NON_NULL(ro2, false);
 
     PS_ASSERT_INT_NONNEGATIVE(footprint, false);