Changeset 29060 for branches/sc_branches/trunkTest/psModules/src/imcombine/pmSubtractionStamps.c

Timestamp:

Aug 26, 2010, 9:18:39 AM (16 years ago)

Author:

Serge CHASTEL

Message:

Merging trunk in branch

Location:

branches/sc_branches/trunkTest

Files:

• . (modified) (1 prop)
• psModules (modified) (1 prop)
• psModules/src/imcombine/pmSubtractionStamps.c (modified) (11 diffs)

branches/sc_branches/trunkTest/psModules/src/imcombine/pmSubtractionStamps.c

@@ -7 +7 @@
 #include <pslib.h>
 
+#include "pmErrorCodes.h"
+#include "pmHDU.h"
+#include "pmFPA.h"
+
 // All these includes required to get stamps out of an array of pmSources
-#include "pmMoments.h"
+#include "pmTrend2D.h"
+#include "pmResiduals.h"
+#include "pmGrowthCurve.h"
 #include "pmSpan.h"
+#include "pmFootprintSpans.h"
 #include "pmFootprint.h"
 #include "pmPeaks.h"
-#include "pmResiduals.h"
-#include "pmHDU.h"
-#include "pmFPA.h"
-#include "pmGrowthCurve.h"
-#include "pmTrend2D.h"
-#include "pmPSF.h"
+#include "pmMoments.h"
+#include "pmModelFuncs.h"
 #include "pmModel.h"
+#include "pmSourceMasks.h"
+#include "pmSourceExtendedPars.h"
+#include "pmSourceDiffStats.h"
 #include "pmSource.h"
-#include "pmErrorCodes.h"
 
 #include "pmSubtraction.h"
@@ -46 +51 @@
     psFree(list->y);
     psFree(list->flux);
-    psFree(list->window);
+    psFree(list->window1);
+    psFree(list->window2);
 }
 
@@ -225 +231 @@
     list->y = NULL;
     list->flux = NULL;
-    list->window = NULL;
     list->normFrac = normFrac;
-    list->normWindow = 0;
+    list->window1 = NULL;
+    list->window2 = NULL;
+    list->normWindow1 = 0;
+    list->normWindow2 = 0;
     list->footprint = footprint;
     list->sysErr = sysErr;
@@ -248 +256 @@
     out->y = NULL;
     out->flux = NULL;
-    out->window = psMemIncrRefCounter(in->window);
+    out->window1 = psMemIncrRefCounter(in->window1);
+    out->window2 = psMemIncrRefCounter(in->window2);
     out->footprint = in->footprint;
-    out->normWindow = in->normWindow;
+    out->normWindow1 = in->normWindow1;
+    out->normWindow2 = in->normWindow2;
 
     for (int i = 0; i < num; i++) {
@@ -638 +648 @@
     int size = stamps->footprint; // Size of postage stamps
 
-    psFree (stamps->window);
-    stamps->window = psKernelAlloc(-size, size, -size, size);
-    psImageInit(stamps->window->image, 0.0);
+    psFree (stamps->window1);
+    stamps->window1 = psKernelAlloc(-size, size, -size, size);
+    psImageInit(stamps->window1->image, 0.0);
+
+    psFree (stamps->window2);
+    stamps->window2 = psKernelAlloc(-size, size, -size, size);
+    psImageInit(stamps->window2->image, 0.0);
 
     // generate normalizations for each stamp
@@ -669 +683 @@
 
     // generate the window pixels
-    double sum = 0.0;                   // Sum inside the window
-    float maxValue = 0.0;               // Maximum value, for normalisation
+    double sum1 = 0.0;                   // Sum inside the window
+    double sum2 = 0.0;                   // Sum inside the window
+    float maxValue1 = 0.0;               // Maximum value, for normalisation
+    float maxValue2 = 0.0;               // Maximum value, for normalisation
     for (int y = -size; y <= size; y++) {
         for (int x = -size; x <= size; x++) {
@@ -691 +707 @@
             }
             float f1 = stats->robustMedian;
+
             psStatsInit (stats);
             if (!psVectorStats (stats, flux2, NULL, NULL, 0)) {
@@ -697 +714 @@
             float f2 = stats->robustMedian;
 
-            stamps->window->kernel[y][x] = f1 + f2;
+            stamps->window1->kernel[y][x] = f1;
             if (PS_SQR(x) + PS_SQR(y) <= PS_SQR(size)) {
-                sum += stamps->window->kernel[y][x];
-            }
-            maxValue = PS_MAX(maxValue, stamps->window->kernel[y][x]);
-        }
-    }
-
-    psTrace("psModules.imcombine", 3, "Window total: %f, threshold: %f\n",
-            sum, (1.0 - stamps->normFrac) * sum);
-    bool done = false;
-    for (int radius = 1; radius <= size && !done; radius++) {
-        double within = 0.0;
+                sum1 += stamps->window1->kernel[y][x];
+            }
+            maxValue1 = PS_MAX(maxValue1, stamps->window1->kernel[y][x]);
+
+            stamps->window2->kernel[y][x] = f2;
+            if (PS_SQR(x) + PS_SQR(y) <= PS_SQR(size)) {
+                sum2 += stamps->window2->kernel[y][x];
+            }
+            maxValue2 = PS_MAX(maxValue2, stamps->window2->kernel[y][x]);
+        }
+    }
+
+    psTrace("psModules.imcombine", 3, "Window total (1): %f, threshold: %f\n", sum1, (1.0 - stamps->normFrac) * sum1);
+    psTrace("psModules.imcombine", 3, "Window total (2): %f, threshold: %f\n", sum2, (1.0 - stamps->normFrac) * sum2);
+
+# if (0)
+    // this block attempts to calculate the radius based on the first radial moment
+    bool done1 = false;
+    bool done2 = false;
+    double prior1 = 0.0;
+    double prior2 = 0.0;
+    for (int y = -size; y <= size; y++) {
+        for (int x = -size; x <= size; x++) {
+            float r = hypot(x, y);
+            Sr1 += r * stamps->window1->kernel[y][x];
+            Sr2 += r * stamps->window2->kernel[y][x];
+            Sf1 += stamps->window1->kernel[y][x];
+            Sf2 += stamps->window2->kernel[y][x];
+        }
+    }
+    
+    float R1 = Sr1 / Sf1;
+    float R2 = Sr2 / Sf2;
+
+    stamps->normWindow1 = 2.5*R1;
+    stamps->normWindow1 = 2.5*R2;
+# else
+    // XXX : this block attempts to calculate the radius by looking at the curve of growth (or something vaguely equivalent).
+    // It did not do very well (though a true curve-of-growth analysis might be better...)
+    bool done1 = false;
+    bool done2 = false;
+    double prior1 = 0.0;
+    double prior2 = 0.0;
+    double delta1o = 1.0;
+    double delta2o = 1.0;
+    for (int radius = 1; radius <= size && !(done1 && done2); radius++) {
+        double within1 = 0.0;
+        double within2 = 0.0;
         for (int y = -radius; y <= radius; y++) {
             for (int x = -radius; x <= radius; x++) {
                 if (PS_SQR(x) + PS_SQR(y) <= PS_SQR(radius)) {
-                    within += stamps->window->kernel[y][x];
+                    within1 += stamps->window1->kernel[y][x];
                 }
-            }
-        }
-        psTrace("psModules.imcombine", 5, "Radius %d: %f\n", radius, within);
-        if (within > (1.0 - stamps->normFrac) * sum) {
-            stamps->normWindow = radius;
-            done = true;
-        }
-    }
-
-    psTrace("psModules.imcombine", 3, "Normalisation window radius set to %d\n", stamps->normWindow);
-    if (stamps->normWindow == 0 || stamps->normWindow >= size) {
-        psError(PM_ERR_STAMPS, true, "Unable to determine normalisation window size.");
+                if (PS_SQR(x) + PS_SQR(y) <= PS_SQR(radius)) {
+                    within2 += stamps->window2->kernel[y][x];
+                }
+            }
+        }
+        double delta1 = (within1 - prior1) / within1;
+        if (!done1 && (fabs(delta1) < stamps->normFrac)) {
+            // interpolate to the radius at which delta2 is normFrac:
+            stamps->normWindow1 = radius - (stamps->normFrac - delta1) / (delta1o - delta1);
+            fprintf (stderr, "choosing %f (%d) for 1 (%f : %f)\n", stamps->normWindow1, radius, within1, delta1);
+            done1 = true;
+        }
+        double delta2 = (within2 - prior2) / within2;
+        if (!done2 && (fabs(delta2) < stamps->normFrac)) {
+            // interpolate to the radius at which delta2 is normFrac:
+            stamps->normWindow2 = radius - (stamps->normFrac - delta2) / (delta2o - delta2);
+            fprintf (stderr, "choosing %f (%d) for 2 (%f : %f)\n", stamps->normWindow2, radius, within2, delta2);
+            done2 = true;
+        }
+        psTrace("psModules.imcombine", 5, "Radius %d: %f (%f) and %f (%f)\n", radius, within1, delta1, within2, delta2);
+
+        prior1 = within1;
+        prior2 = within2;
+        delta1o = delta1;
+        delta2o = delta2;
+
+        // if (!done1 && (within1 > (1.0 - stamps->normFrac) * sum1)) {
+        //     stamps->normWindow1 = radius;
+        //     done1 = true;
+        // }
+        // if (!done2 && (within2 > (1.0 - stamps->normFrac) * sum2)) {
+        //     stamps->normWindow2 = radius;
+        //     done2 = true;
+        // }
+
+    }
+# endif
+
+    psTrace("psModules.imcombine", 3, "Normalisation window radii set to %f and %f\n", stamps->normWindow1, stamps->normWindow2);
+    if (stamps->normWindow1 == 0 || stamps->normWindow1 >= size) {
+        psError(PM_ERR_STAMPS, true, "Unable to determine normalisation window size (1).");
+        return false;
+    }
+    if (stamps->normWindow2 == 0 || stamps->normWindow2 >= size) {
+        psError(PM_ERR_STAMPS, true, "Unable to determine normalisation window size (2).");
         return false;
     }
@@ -733 +820 @@
     for (int y = -size; y <= size; y++) {
         for (int x = -size; x <= size; x++) {
-            stamps->window->kernel[y][x] /= maxValue;
+            stamps->window1->kernel[y][x] /= maxValue1;
+        }
+    }
+    // re-normalize so chisquare values are sensible
+    for (int y = -size; y <= size; y++) {
+        for (int x = -size; x <= size; x++) {
+            stamps->window2->kernel[y][x] /= maxValue2;
         }
     }
@@ -739 +832 @@
 #if 0
     {
-        psFits *fits = psFitsOpen ("window.fits", "w");
-        psFitsWriteImage (fits, NULL, stamps->window->image, 0, NULL);
+        psFits *fits = NULL;
+        fits = psFitsOpen ("window1.fits", "w");
+        psFitsWriteImage (fits, NULL, stamps->window1->image, 0, NULL);
+        psFitsClose (fits);
+        fits = psFitsOpen ("window2.fits", "w");
+        psFitsWriteImage (fits, NULL, stamps->window2->image, 0, NULL);
         psFitsClose (fits);
     }
@@ -747 +844 @@
     psFree (stats);
     psFree (flux1);
-    psFree(flux2);
+    psFree (flux2);
     psFree (norm1);
     psFree (norm2);
     return true;
+}
+
+static pthread_mutex_t getPenaltiesMutex = PTHREAD_MUTEX_INITIALIZER;
+
+// kernels->penalty is an overall scaling factor (user-supplied)
+bool pmSubtractionKernelPenaltiesStamp(pmSubtractionStamp *stamp, pmSubtractionKernels *kernels)
+{
+    // we only need the penalties if we are doing dual convolution
+    if (kernels->mode != PM_SUBTRACTION_MODE_DUAL) return true;
+
+    // we only calculate the penalties once.
+    if (kernels->havePenalties) return true;
+
+    // in a threaded context, only one thread can calculate the penalties.  attempt to grab a
+    // mutex before continuing
+    pthread_mutex_lock(&getPenaltiesMutex);
+
+    // did someone else already get the mutex and do this?
+    if (kernels->havePenalties) {
+        pthread_mutex_unlock(&getPenaltiesMutex);
+        return true;
+    }
+
+    for (int i = 0; i < kernels->num; i++) {
+        pmSubtractionKernelPenalties(stamp, kernels, i);
+    }
+
+    kernels->havePenalties = true;
+    pthread_mutex_unlock(&getPenaltiesMutex);
+    return true;
+}
+
+# define EMPIRICAL 0
+
+// kernels->penalty is an overall scaling factor (user-supplied)
+bool pmSubtractionKernelPenalties(pmSubtractionStamp *stamp, pmSubtractionKernels *kernels, int index)
+{
+    float penalty1, penalty2;
+    float fwhm1, fwhm2;
+
+    // XXX this is annoyingly hack-ish
+    pmSubtractionGetFWHMs(&fwhm1, &fwhm2);
+
+    bool zeroNull = false;
+    int uOrder = kernels->u->data.S32[index];
+    int vOrder = kernels->v->data.S32[index];
+    if (uOrder % 2 == 0 && vOrder % 2 == 0) zeroNull = true;
+
+    if (EMPIRICAL) {
+        psKernel *convolution1 = stamp->convolutions1->data[index];
+        penalty1 = pmSubtractionKernelPenaltySingle(convolution1, zeroNull);
+
+        psKernel *convolution2 = stamp->convolutions2->data[index];
+        penalty2 = pmSubtractionKernelPenaltySingle(convolution2, zeroNull);
+    } else {
+        pmSubtractionKernelPreCalc *kernel = kernels->preCalc->data[index];
+        float M2 = pmSubtractionKernelPenaltySingle(kernel->kernel, zeroNull);
+
+        penalty1 = M2 + PS_SQR(fwhm1 / 2.35); // rescale the unconvolved second-moment by the image second moment 
+        penalty2 = M2 + PS_SQR(fwhm2 / 2.35); // rescale the unconvolved second-moment by the image second moment 
+    }
+    kernels->penalties1->data.F32[index] = kernels->penalty * penalty1;
+    psAssert (isfinite(kernels->penalties1->data.F32[index]), "invalid penalty");
+
+    kernels->penalties2->data.F32[index] = kernels->penalty * penalty2;
+    psAssert (isfinite(kernels->penalties2->data.F32[index]), "invalid penalty");
+
+    fprintf(stderr, "penalty1: %f, penalty2: %f\n", penalty1, penalty2);
+
+    return true;
+}
+
+float pmSubtractionKernelPenaltySingle(psKernel *kernel, bool zeroNull)
+{
+    // Calculate moments
+    double penalty = 0.0;                   // Moment, for penalty
+    double sum = 0.0, sum2 = 0.0;           // Sum of kernel component
+    float min = INFINITY, max = -INFINITY;  // Minimum and maximum kernel value
+    for (int v = kernel->yMin; v <= kernel->yMax; v++) {
+        for (int u = kernel->xMin; u <= kernel->xMax; u++) {
+            double value = kernel->kernel[v][u];
+            if (false && zeroNull && (u == 0) && (v == 0)) {
+                value += 1.0;
+            }
+            double value2 = PS_SQR(value);
+            sum += value;
+            sum2 += value2;
+            penalty += value2 * PS_SQR((PS_SQR(u) + PS_SQR(v)));
+            min = PS_MIN(value, min);
+            max = PS_MAX(value, max);
+        }
+    }
+    penalty *= 1.0 / sum2;
+
+    if (1) {
+        fprintf(stderr, "min: %lf, max: %lf, moment: %lf, flux^2: %lf\n", min, max, penalty, sum2);
+        // psTrace("psModules.imcombine", 7, "Kernel %d: %f %d %d %f\n", index, fwhm, uOrder, vOrder, penalty);
+    }
+
+    return penalty;
 }