Changeset 15516 for branches/eam_branch_20071023/psModules/src/imcombine/pmSubtractionMatch.c

Timestamp:

Nov 8, 2007, 12:59:04 PM (19 years ago)

Author:

eugene

Message:

accepting updates from HEAD

File:

• branches/eam_branch_20071023/psModules/src/imcombine/pmSubtractionMatch.c (modified) (17 diffs)

branches/eam_branch_20071023/psModules/src/imcombine/pmSubtractionMatch.c

@@ -50 +50 @@
 
 static bool getStamps(pmSubtractionStampList **stamps, // Stamps to read
-                      const pmReadout *reference, // Reference readout
-                      const pmReadout *input, // Input readout, or NULL to generate fake stamps
+                      const pmReadout *ro1, // Readout 1
+                      const pmReadout *ro2, // Readout 2
                       const psImage *subMask, // Mask for subtraction, or NULL
                       psImage *weight,  // Weight map
@@ -57 +57 @@
                       float threshold,  // Threshold for stamp finding
                       float stampSpacing, // Spacing between stamps
-                      float targetWidth,// Target width for fake stamps
                       int size,         // Kernel half-size
-                      int footprint     // Convolution footprint for stamps
+                      int footprint,     // Convolution footprint for stamps
+                      pmSubtractionMode mode // Mode for subtraction
     )
 {
     psTrace("psModules.imcombine", 3, "Finding stamps...\n");
-    *stamps = pmSubtractionStampsFind(*stamps, reference->image, subMask, region, threshold, stampSpacing);
+    *stamps = pmSubtractionStampsFind(*stamps, ro1->image, subMask, region, threshold, stampSpacing, mode);
     if (!*stamps) {
         psError(PS_ERR_UNKNOWN, false, "Unable to find stamps.");
@@ -71 +71 @@
     memCheck("  find stamps");
 
-    if (!input && !pmSubtractionStampsGenerate(*stamps, targetWidth, footprint, size)) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to generate target stamps.");
-        return false;
-    }
-
-    memCheck("   generate stamps");
-
     psTrace("psModules.imcombine", 3, "Extracting stamps...\n");
-    if (!pmSubtractionStampsExtract(*stamps, reference->image, input ? input->image : NULL,
-                                    weight, footprint, size)) {
+    if (!pmSubtractionStampsExtract(*stamps, ro1->image, ro2 ? ro2->image : NULL, weight, footprint, size)) {
         psError(PS_ERR_UNKNOWN, false, "Unable to extract stamps.");
         return false;
@@ -94 +86 @@
 
 
-bool pmSubtractionMatch(pmReadout *convolved, const pmReadout *reference, const pmReadout *input,
+bool pmSubtractionMatch(pmReadout *convolved, const pmReadout *ro1, const pmReadout *ro2,
                         int footprint, float regionSize, float stampSpacing, float threshold,
-                        const psArray *sources, const char *stampsName, float targetWidth,
+                        const psArray *sources, const char *stampsName,
                         pmSubtractionKernelsType type, int size, int spatialOrder,
                         const psVector *isisWidths, const psVector *isisOrders,
                         int inner, int ringsOrder, int binning, bool optimum, const psVector *optFWHMs,
                         int optOrder, float optThreshold, int iter, float rej, psMaskType maskBad,
-                        psMaskType maskBlank, float badFrac)
+                        psMaskType maskBlank, float badFrac, pmSubtractionMode mode)
 {
     PS_ASSERT_PTR_NON_NULL(convolved, false);
-    PS_ASSERT_PTR_NON_NULL(reference, false);
-    PS_ASSERT_IMAGE_NON_NULL(reference->image, false);
-    PS_ASSERT_IMAGE_TYPE(reference->image, PS_TYPE_F32, false);
-    if (reference->mask) {
-        PS_ASSERT_IMAGE_NON_NULL(reference->mask, false);
-        PS_ASSERT_IMAGE_TYPE(reference->mask, PS_TYPE_MASK, false);
-        PS_ASSERT_IMAGES_SIZE_EQUAL(reference->mask, reference->image, false);
-    }
-    if (reference->weight) {
-        PS_ASSERT_IMAGE_NON_NULL(reference->weight, false);
-        PS_ASSERT_IMAGE_TYPE(reference->weight, PS_TYPE_F32, false);
-        PS_ASSERT_IMAGES_SIZE_EQUAL(reference->weight, reference->image, false);
-    }
-    if (input) {
-        PS_ASSERT_IMAGE_NON_NULL(input->image, false);
-        PS_ASSERT_IMAGE_TYPE(input->image, PS_TYPE_F32, false);
-        PS_ASSERT_IMAGES_SIZE_EQUAL(input->image, reference->image, false);
-        if (input->mask) {
-            PS_ASSERT_IMAGE_NON_NULL(input->mask, false);
-            PS_ASSERT_IMAGE_TYPE(input->mask, PS_TYPE_MASK, false);
-            PS_ASSERT_IMAGES_SIZE_EQUAL(input->mask, reference->image, false);
-        }
-        if (input->weight) {
-            PS_ASSERT_IMAGE_NON_NULL(input->weight, false);
-            PS_ASSERT_IMAGE_TYPE(input->weight, PS_TYPE_F32, false);
-            PS_ASSERT_IMAGES_SIZE_EQUAL(input->weight, reference->image, false);
+    PS_ASSERT_PTR_NON_NULL(ro1, false);
+    PS_ASSERT_IMAGE_NON_NULL(ro1->image, false);
+    PS_ASSERT_IMAGE_TYPE(ro1->image, PS_TYPE_F32, false);
+    if (ro1->mask) {
+        PS_ASSERT_IMAGE_NON_NULL(ro1->mask, false);
+        PS_ASSERT_IMAGE_TYPE(ro1->mask, PS_TYPE_MASK, false);
+        PS_ASSERT_IMAGES_SIZE_EQUAL(ro1->mask, ro1->image, false);
+    }
+    if (ro1->weight) {
+        PS_ASSERT_IMAGE_NON_NULL(ro1->weight, false);
+        PS_ASSERT_IMAGE_TYPE(ro1->weight, PS_TYPE_F32, false);
+        PS_ASSERT_IMAGES_SIZE_EQUAL(ro1->weight, ro1->image, false);
+    }
+    if (ro2) {
+        PS_ASSERT_IMAGE_NON_NULL(ro2->image, false);
+        PS_ASSERT_IMAGE_TYPE(ro2->image, PS_TYPE_F32, false);
+        PS_ASSERT_IMAGES_SIZE_EQUAL(ro2->image, ro1->image, false);
+        if (ro2->mask) {
+            PS_ASSERT_IMAGE_NON_NULL(ro2->mask, false);
+            PS_ASSERT_IMAGE_TYPE(ro2->mask, PS_TYPE_MASK, false);
+            PS_ASSERT_IMAGES_SIZE_EQUAL(ro2->mask, ro1->image, false);
+        }
+        if (ro2->weight) {
+            PS_ASSERT_IMAGE_NON_NULL(ro2->weight, false);
+            PS_ASSERT_IMAGE_TYPE(ro2->weight, PS_TYPE_F32, false);
+            PS_ASSERT_IMAGES_SIZE_EQUAL(ro2->weight, ro1->image, false);
         }
     } else if (!stampsName && !sources) {
@@ -144 +136 @@
     }
     // stampsName may be anything
-    // targetWidth can be just about anything (except maybe negative, but it can be NAN)
     // We'll check kernel type when we allocate the kernels
     PS_ASSERT_INT_POSITIVE(size, false);
@@ -196 +187 @@
     }
 
-    psImage *inImage = NULL, *inMask = NULL; // Input image, mask, weight
-    if (input) {
-        inImage = input->image;
-        inMask = input->mask;
-    }
-
     // Where does our weight map come from?
     psImage *weight = NULL;             // Weight image to use
-    if (input && input->weight) {
-        weight = input->weight;
-    } else if (reference->weight) {
-        weight = reference->weight;
-    } else if (input) {
-        weight = input->image;
+    if (ro1->weight && ro2 && ro2->weight) {
+        weight = (psImage*)psBinaryOp(NULL, ro1->weight, "+", ro2->weight);
+    } else if (ro1->weight) {
+        weight = psMemIncrRefCounter(ro1->weight);
+    } else if (ro2) {
+        if (ro2->weight) {
+            weight = psMemIncrRefCounter(ro2->weight);
+        } else {
+            weight = (psImage*)psBinaryOp(NULL, ro1->image, "+", ro2->image);
+        }
     } else {
-        weight = reference->image;
+        weight = psMemIncrRefCounter(ro1->image);
     }
 
@@ -221 +210 @@
     psVector *solution = NULL;          // Solution to match PSF
     pmSubtractionKernels *kernels = NULL; // Kernel basis functions
-
-    int numCols = reference->image->numCols, numRows = reference->image->numRows; // Image dimensions
+    int numCols = ro1->image->numCols, numRows = ro1->image->numRows; // Image dimensions
 
     memCheck("start");
 
-    subMask = pmSubtractionMask(reference->mask, inMask, maskBad, size, footprint, badFrac, useFFT);
+    subMask = pmSubtractionMask(ro1->mask, ro2 ? ro2->mask : NULL, maskBad, size, footprint,
+                                badFrac, useFFT);
     if (!subMask) {
         psError(PS_ERR_UNKNOWN, false, "Unable to generate subtraction mask.");
+        psFree(weight);
         return false;
     }
@@ -260 +250 @@
 
             if (sources) {
-                stamps = pmSubtractionStampsSetFromSources(sources, subMask, region, stampSpacing,
-                                                           input ? 0 : 2 * footprint);
+                stamps = pmSubtractionStampsSetFromSources(sources, subMask, region, stampSpacing, mode);
             } else if (stampsName && strlen(stampsName) > 0) {
-                // Read stamps from file
-                psTrace("psModules.imcombine", 3, "Reading stamps from %s...\n", stampsName);
-                const char *stampFormat = input ? "%f %f" : "%f %f %f"; // Format for reading stamp file
-                psArray *stampsData = stampsData = psVectorsReadFromFile(stampsName, stampFormat);
-                psVector *xStamp = NULL, *yStamp = NULL, *fluxStamp = NULL; // Stamp positions and fluxes
-                if (!stampsData) {
-                    psError(PS_ERR_IO, false, "Unable to read stamps file %s", stampsName);
-                    goto ERROR;
-                }
-                xStamp = stampsData->data[0];
-                yStamp = stampsData->data[1];
-                if (!input) {
-                    fluxStamp = stampsData->data[2];
-                }
-
-                // Correct for IRAF/FITS (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 = pmSubtractionStampsSet(xStamp, yStamp, fluxStamp, reference->image, subMask,
-                                                region, stampSpacing, input ? 0 : footprint + size);
-                psFree(stampsData);
+                stamps = pmSubtractionStampsSetFromFile(stampsName, subMask, region, stampSpacing, mode);
+            }
+
+            // We get the stamps here; we will also attempt to get stamps at the first iteration, but it
+            // doesn't matter.
+            if (!getStamps(&stamps, ro1, ro2, subMask, weight, NULL, threshold, stampSpacing,
+                           size, footprint, mode)) {
+                goto MATCH_ERROR;
+            }
+
+            if (mode == PM_SUBTRACTION_MODE_UNSURE || mode == PM_SUBTRACTION_MODE_TARGET) {
+                pmSubtractionMode newMode = pmSubtractionOrder(stamps, footprint); // Subtraction mode
+                switch (newMode) {
+                  case PM_SUBTRACTION_MODE_1:
+                    psLogMsg("psModules.imcombine", PS_LOG_INFO, "Convolving image 1 to match image 2.");
+                    break;
+                  case PM_SUBTRACTION_MODE_2:
+                    if (mode == PM_SUBTRACTION_MODE_TARGET) {
+                        psError(PS_ERR_BAD_PARAMETER_VALUE, true,
+                                "Input PSF is larger than target PSF --- can't match image.");
+                        goto MATCH_ERROR;
+                    }
+                    psLogMsg("psModules.imcombine", PS_LOG_INFO, "Convolving image 2 to match image 1.");
+                    break;
+                  default:
+                    psError(PS_ERR_UNKNOWN, false, "Unable to determine subtraction order.");
+                    goto MATCH_ERROR;
+                }
+                mode = newMode;
             }
 
             // Define kernel basis functions
             if (optimum && (type == PM_SUBTRACTION_KERNEL_ISIS || type == PM_SUBTRACTION_KERNEL_GUNK)) {
-                if (!getStamps(&stamps, reference, input, subMask, weight, NULL,
-                               threshold, stampSpacing, targetWidth, size, footprint)) {
-                    goto ERROR;
-                }
                 kernels = pmSubtractionKernelsOptimumISIS(type, size, inner, spatialOrder, optFWHMs, optOrder,
-                                                          stamps, footprint, optThreshold);
+                                                          stamps, footprint, optThreshold, mode);
                 if (!kernels) {
                     psErrorClear();
@@ -314 +306 @@
                 psLogMsg("psModules.imcombine", PS_LOG_INFO, "Iteration %d.", k);
 
-                if (!getStamps(&stamps, reference, input, subMask, weight, region,
-                               threshold, stampSpacing, targetWidth, size, footprint)) {
-                    goto ERROR;
+                if (!getStamps(&stamps, ro1, ro2, subMask, weight, region, threshold, stampSpacing,
+                               size, footprint, mode)) {
+                    goto MATCH_ERROR;
                 }
 
                 psTrace("psModules.imcombine", 3, "Calculating equation...\n");
-                if (!pmSubtractionCalculateEquation(stamps, kernels, footprint)) {
+                if (!pmSubtractionCalculateEquation(stamps, kernels, footprint, mode)) {
                     psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation.");
-                    goto ERROR;
+                    goto MATCH_ERROR;
                 }
 
@@ -331 +323 @@
                 if (!solution) {
                     psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation.");
-                    goto ERROR;
+                    goto MATCH_ERROR;
                 }
 
                 memCheck("  solve equation");
 
+                psVector *deviations = pmSubtractionCalculateDeviations(stamps, solution, footprint, kernels,
+                                                                        mode); // Deviations for each stamp
+                if (!deviations) {
+                    psError(PS_ERR_UNKNOWN, false, "Unable to calculate deviations.");
+                    goto MATCH_ERROR;
+                }
+
+                memCheck("   calculate deviations");
+
                 psTrace("psModules.imcombine", 3, "Rejecting stamps...\n");
-                numRejected = pmSubtractionRejectStamps(stamps, subMask, solution, footprint, rej, kernels);
+                numRejected = pmSubtractionRejectStamps(stamps, deviations, subMask, rej, footprint);
                 if (numRejected < 0) {
                     psError(PS_ERR_UNKNOWN, false, "Unable to reject stamps.");
-                    goto ERROR;
-                }
+                    psFree(deviations);
+                    goto MATCH_ERROR;
+                }
+                psFree(deviations);
+
                 memCheck("  reject stamps");
             }
@@ -350 +354 @@
                 if (!solution) {
                     psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation.");
-                    goto ERROR;
-                }
-                (void)pmSubtractionRejectStamps(stamps, subMask, solution, footprint, NAN, kernels);
+                    goto MATCH_ERROR;
+                }
+                psVector *deviations = pmSubtractionCalculateDeviations(stamps, solution, footprint, kernels,
+                                                                        mode); // Deviations for each stamp
+                if (!deviations) {
+                    psError(PS_ERR_UNKNOWN, false, "Unable to calculate deviations.");
+                    goto MATCH_ERROR;
+                }
+                pmSubtractionRejectStamps(stamps, deviations, subMask, NAN, footprint);
+                psFree(deviations);
             }
             psFree(stamps);
@@ -372 +383 @@
                             psError(PS_ERR_UNKNOWN, false, "Unable to generate kernel image.");
                             psFree(convKernels);
-                            goto ERROR;
+                            goto MATCH_ERROR;
                         }
 
@@ -380 +391 @@
                             psFree(kernel);
                             psFree(convKernels);
-                            goto ERROR;
+                            goto MATCH_ERROR;
                         }
                         psFree(kernel);
@@ -416 +427 @@
 
             psTrace("psModules.imcombine", 2, "Convolving...\n");
-            if (!pmSubtractionConvolve(&convolved->image, &convolved->weight, &convolved->mask,
-                                       reference->image, reference->weight, subMask, maskBlank, region,
-                                       solution, kernels, useFFT)) {
-                psError(PS_ERR_UNKNOWN, false, "Unable to convolve reference image.");
-                goto ERROR;
+            if (!pmSubtractionConvolve(convolved, ro1, ro2, subMask, maskBlank, region,
+                                       solution, kernels, mode, useFFT)) {
+                psError(PS_ERR_UNKNOWN, false, "Unable to convolve image.");
+                goto MATCH_ERROR;
             }
             psFree(kernels);
@@ -461 +471 @@
     psFree(subMask);
     subMask = NULL;
+    psFree(weight);
+    weight = NULL;
 
     if (!pmSubtractionBorder(convolved->image, convolved->weight, convolved->mask, size, maskBlank)) {
         psError(PS_ERR_UNKNOWN, false, "Unable to set border of convolved image.");
-        goto ERROR;
+        goto MATCH_ERROR;
     }
 
@@ -472 +484 @@
     return true;
 
-ERROR:
+MATCH_ERROR:
     psFree(region);
     psFree(regionString);
@@ -479 +491 @@
     psFree(stamps);
     psFree(solution);
+    psFree(weight);
     return false;
 }
+
+// Calculate the second order moments for an image
+static float subtractionOrderMoment(const psKernel *kernel, // Image for which to measure moments
+                                    int radius   // Maximum radius
+                                    )
+{
+    assert(kernel && kernel->kernel);
+
+    int xMin = PS_MAX(kernel->xMin, -radius), xMax = PS_MIN(kernel->xMax, radius); // Bounds in x
+    int yMin = PS_MAX(kernel->yMin, -radius), yMax = PS_MIN(kernel->yMax, radius); // Bounds in y
+
+    float xCentroid = 0.0, yCentroid = 0.0; // Centroid (first moment)
+    float sum = 0.0;       // Sum (zero-th moment)
+    for (int y = yMin; y <= yMax; y++) {
+        for (int x = xMin; x <= xMax; x++) {
+            xCentroid += kernel->kernel[y][x] * x;
+            yCentroid += kernel->kernel[y][x] * y;
+            sum += kernel->kernel[y][x];
+        }
+    }
+    xCentroid /= sum;
+    yCentroid /= sum;
+
+    float eta20 = 0.0, eta02 = 0.0;     // Second moments
+    for (int y = yMin; y <= yMax; y++) {
+        float yDiff = y - yCentroid;
+        for (int x = xMin; x <= xMax; x++) {
+            float xDiff = x - xCentroid;
+            eta20 += PS_SQR(xDiff) * kernel->kernel[y][x];
+            eta02 += PS_SQR(yDiff) * kernel->kernel[y][x];
+        }
+    }
+
+    // Normalise to calculate the scale-invariant
+    float sum2 = PS_SQR(sum);
+    eta20 /= sum2;
+    eta02 /= sum2;
+    // eta11 /= sum2;
+
+    return eta20 + eta02;
+}
+
+#if 0
+// Calculate the deviations for a particular subtraction order
+static psVector *subtractionOrderDeviation(float *sumKernel, // Sum of the kernel
+                                           pmSubtractionStampList *stamps, // Stamps to convolve
+                                           const pmSubtractionKernels *kernels, // Kernel basis functions
+                                           int footprint, // Stamp footprint
+                                           pmSubtractionMode mode // Mode of subtraction
+                                           )
+{
+    assert(stamps);
+    assert(footprint >= 0);
+    assert(mode == PM_SUBTRACTION_MODE_1 || mode == PM_SUBTRACTION_MODE_2);
+
+    if (!pmSubtractionCalculateEquation(stamps, kernels, footprint, mode)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation.");
+        return NULL;
+    }
+
+    psVector *solution = pmSubtractionSolveEquation(NULL, stamps);
+    if (!solution) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation.");
+        return NULL;
+    }
+
+    if (sumKernel) {
+        float sum = 0.0;                // Sum of the kernel
+        psImage *image = pmSubtractionKernelImage(solution, kernels, 0.0, 0.0); // Image of kernel
+        for (int y = 0; y < image->numRows; y++) {
+            for (int x = 0; x < image->numCols; x++) {
+                sum += image->data.F32[y][x];
+            }
+        }
+        psFree(image);
+        *sumKernel = sum;
+    }
+
+    psVector *deviations = pmSubtractionCalculateDeviations(stamps, solution, footprint, kernels, mode);
+    psFree(solution);
+    if (!deviations) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to calculate deviations.");
+        return NULL;
+    }
+
+    return deviations;
+}
+#endif
+
+pmSubtractionMode pmSubtractionOrder(pmSubtractionStampList *stamps, int radius)
+{
+    PS_ASSERT_INT_POSITIVE(radius, PM_SUBTRACTION_MODE_ERR);
+
+    psVector *mask = psVectorAlloc(stamps->num, PS_TYPE_MASK); // Mask for stamps
+    psVector *moments = psVectorAlloc(stamps->num, PS_TYPE_F32); // Moments
+    for (int i = 0; i < stamps->num; i++) {
+        pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest
+        if (stamp->status != PM_SUBTRACTION_STAMP_CALCULATE && stamp->status != PM_SUBTRACTION_STAMP_USED) {
+            mask->data.PS_TYPE_MASK_DATA[i] = 0xff;
+            continue;
+        }
+        mask->data.PS_TYPE_MASK_DATA[i] = 0;
+        moments->data.F32[i] = subtractionOrderMoment(stamp->image1, radius) /
+            subtractionOrderMoment(stamp->image2, radius);
+    }
+
+    psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN);
+    if (!psVectorStats(stats, moments, NULL, mask, 0xff)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to calculate statistics for moments ratio.");
+        psFree(mask);
+        psFree(moments);
+        psFree(stats);
+        return PM_SUBTRACTION_MODE_ERR;
+    }
+    psFree(moments);
+    psFree(mask);
+
+    psLogMsg("psModules.imcombine", PS_LOG_INFO, "Median ratio of second moments: %lf", stats->robustMedian);
+    pmSubtractionMode mode = (stats->robustMedian <= 1.0 ? PM_SUBTRACTION_MODE_1 : PM_SUBTRACTION_MODE_2);
+    psFree(stats);
+
+    return mode;
+}
+
+
+