Changeset 30288 for branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtraction.c

Timestamp:

Jan 17, 2011, 5:07:46 PM (16 years ago)

Author:

eugene

Message:

test chisq improvements from increasing the spatial order and DUAL vs SINGLE1,2; reject stamps on the basis of the chisq vs flux^{2 model fit (allowing for systematic errors in the fit); track the best fit with the new pmSubtractionQuality structure; split out pmSubtractionConvolveStamps from pmSubtractionCalculateEquation; replace pmSubtractionCalculateDeviations with pmSubtractionCalculateChisqAndMoments}

File:

• branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtraction.c (modified) (8 diffs)

branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtraction.c

@@ -773 +773 @@
 
     if (convolutions) {
-        // Already done
         return convolutions;
     }
@@ -787 +786 @@
 }
 
+
+bool pmSubtractionConvolveStampThread(psThreadJob *job)
+{
+    PS_ASSERT_THREAD_JOB_NON_NULL(job, false);
+
+    pmSubtractionStamp *stamp = job->args->data[0]; // List of stamps
+    pmSubtractionKernels *kernels = job->args->data[1]; // Kernels
+    int footprint = PS_SCALAR_VALUE(job->args->data[2], S32); // Stamp index
+
+    return pmSubtractionConvolveStamp(stamp, kernels, footprint);
+}
 
 bool pmSubtractionConvolveStamp (pmSubtractionStamp *stamp, pmSubtractionKernels *kernels, int footprint)
@@ -818 +828 @@
     }
 
+#ifdef TESTING
+    for (int j = 0; j < kernels->num; j++) {
+        if (stamp->convolutions1) {
+            psString convName = NULL;
+            psStringAppend(&convName, "conv1_%03d_%03d.fits", index, j);
+            psFits *fits = psFitsOpen(convName, "w");
+            psFree(convName);
+            psKernel *conv = stamp->convolutions1->data[j];
+            psFitsWriteImage(fits, NULL, conv->image, 0, NULL);
+            psFitsClose(fits);
+        }
+
+        if (stamp->convolutions2) {
+            psString convName = NULL;
+            psStringAppend(&convName, "conv2_%03d_%03d.fits", index, j);
+            psFits *fits = psFitsOpen(convName, "w");
+            psFree(convName);
+            psKernel *conv = stamp->convolutions2->data[j];
+            psFitsWriteImage(fits, NULL, conv->image, 0, NULL);
+            psFitsClose(fits);
+        }
+    }
+#endif
+
     return true;
 }
 
+bool pmSubtractionConvolveStamps(pmSubtractionStampList *stamps, pmSubtractionKernels *kernels) 
+{
+    PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false);
+    PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, false);
+
+    psTimerStart("pmSubtractionConvolveStamps");
+
+    int footprint = stamps->footprint;  // Half-size of stamps
+
+    // We iterate over each stamp and generate the convolution if needed.  We do NOT need the
+    // convolution if (a) it has already been calculated or (b) the stamp is not available for
+    // use (available = USED or CALCULATE)
+    
+    for (int i = 0; i < stamps->num; i++) {
+        pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest
+
+        bool keep = false;
+        keep |= (stamp->status == PM_SUBTRACTION_STAMP_USED);
+        keep |= (stamp->status == PM_SUBTRACTION_STAMP_CALCULATE);
+        if (!keep) continue;
+
+        bool haveConvolutions = false;
+        if (kernels->mode == PM_SUBTRACTION_MODE_1) {
+            haveConvolutions = (stamp->convolutions1 != NULL);
+        }
+        if (kernels->mode == PM_SUBTRACTION_MODE_2) {
+            haveConvolutions = (stamp->convolutions2 != NULL);
+        }
+        if (kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
+            haveConvolutions = (stamp->convolutions1 != NULL) && (stamp->convolutions2 != NULL);
+        }
+        if (haveConvolutions) {
+            continue;
+        }
+
+        if (pmSubtractionThreaded()) {
+            psThreadJob *job = psThreadJobAlloc("PSMODULES_SUBTRACTION_CONVOLVE_STAMP");
+            psArrayAdd(job->args, 1, stamp);
+            psArrayAdd(job->args, 1, kernels);
+            PS_ARRAY_ADD_SCALAR(job->args, footprint, PS_TYPE_S32);
+            if (!psThreadJobAddPending(job)) {
+                return false;
+            }
+        } else {
+            pmSubtractionConvolveStamp(stamp, kernels, footprint);
+        }
+    }
+    if (!psThreadPoolWait(true)) {
+        psError(psErrorCodeLast(), false, "Error waiting for threads.");
+        return false;
+    }
+    psLogMsg("psModules.imcombine", PS_LOG_INFO, "Convolve stamps: %f sec", psTimerClear("pmSubtractionConvolveStamps"));
+    return true;
+}
 
 int pmSubtractionRejectStamps(pmSubtractionKernels *kernels, pmSubtractionStampList *stamps,
-                              const psVector *deviations, psImage *subMask, float sigmaRej)
+                              pmSubtractionQuality *match, psImage *subMask, float sigmaRej)
 {
     PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, false);
     PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, -1);
-    PS_ASSERT_VECTOR_NON_NULL(deviations, -1);
-    PS_ASSERT_VECTOR_TYPE(deviations, PS_TYPE_F32, -1);
     PS_ASSERT_IMAGE_NON_EMPTY(subMask, -1);
     PS_ASSERT_IMAGE_TYPE(subMask, PS_TYPE_IMAGE_MASK, -1);
 
-    // I used to measure the rms deviation about zero, and use that as the sigma against which to clip, but
-    // the distribution is actually something like a chi^2 or Student's t, both of which become Gaussian-like
-    // with large N.  Therefore, let's just treat this as a Gaussian distribution.
+    // Comment from PAP (r18287): I used to measure the rms deviation about zero, and use that as the
+    // sigma against which to clip, but the distribution is actually something like a chi^2 or
+    // Student's t, both of which become Gaussian-like with large N.  Therefore, let's just
+    // treat this as a Gaussian distribution.
+
+    // Comment from EAM (r29777): The residual distribution is only chisq-like if the model is
+    // a good fit to the data.  In the (likely) case that there is a systematic difference
+    // between the model and the data, the squared-residual distribution grows quadratically
+    // with increasing flux: the systematic residual flux is a constant factor times the source
+    // flux; the squared-residual is then of the form (k0 + k1*flux)^2, where k0 comes from the
+    // Gaussian distributed residual and k1*flux is the systematic residual error.
+
+    // By rejecting sources with the largest squared-residuals, the rejection biases against
+    // the brighter sources; in severe cases, this pushes the measurement to the weakest
+    // sources with the most noise.  To account for this, let's fit a 2nd order polynomial to
+    // the distribution of flux vs squared-residual, subtract that fit, and reject sources
+    // which are significantly deviant from that distribution.
 
     kernels->mean = NAN;
@@ -840 +940 @@
     kernels->numStamps = -1;
 
-    int numStamps = 0;                  // Number of used stamps
-    psVector *mask = psVectorAlloc(stamps->num, PS_TYPE_VECTOR_MASK); // Mask, for statistics
-    psVectorInit(mask, 0);
-    for (int i = 0; i < stamps->num; i++) {
-        pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest
-        if (stamp->status != PM_SUBTRACTION_STAMP_USED) {
-            mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = 0xff;
-            continue;
-        }
-        numStamps++;
-    }
-    psTrace("psModules.imcombine", 1, "Number of good stamps: %d\n", numStamps);
-
-    if (numStamps == 0) {
-        psError(PM_ERR_STAMPS, true, "No good stamps found.");
-        psFree(mask);
-        return -1;
-    }
-
-    psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV |
-                                  PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_QUARTILE); // Statistics for deviatns
-    if (!psVectorStats(stats, deviations, NULL, mask, 0xff)) {
-        psError(PM_ERR_DATA, false, "Unable to measure statistics for deviations.");
+    psTrace("psModules.imcombine", 1, "Number of good stamps: %d\n", match->nGood);
+
+    // the chisq & flux vectors are calculated by pmSubtractionCalculateChisqAndMoments
+
+    // use 3hi/3lo sigma clipping on the chisq fit
+    psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV);
+    stats->clipSigma = 5.0;
+    stats->clipIter = 2;
+    psPolynomial1D *model = psPolynomial1DAlloc (PS_POLYNOMIAL_ORD, 2);
+
+    bool result = psVectorClipFitPolynomial1D(model, stats, match->stampMask, 0xff, match->chisq, NULL, match->fluxes);
+    if (!result) {
+        psError(PM_ERR_DATA, false, "Unable to measure statistics for deviations.");
+        psFree(model);
         psFree(stats);
-        psFree(mask);
-        return -1;
-    }
-    psFree(mask);
-
-    // XXX raise an error?
+        return -1;
+    }
     if (isnan(stats->sampleMean)) {
+        psError(PM_ERR_DATA, false, "Unable to measure statistics for deviations.");
+        psFree(model);
         psFree(stats);
         return -1;
     }
 
-    double mean, rms;                 // Mean and RMS of deviations
-    if (numStamps < MIN_SAMPLE_STATS) {
-        mean = stats->sampleMean;
-        rms = stats->sampleStdev;
-    } else {
-        mean = stats->sampleMedian;
-        rms = 0.74 * (stats->sampleUQ - stats->sampleLQ);
-    }
-    psFree(stats);
-
-    psTrace("psModules.imcombine", 1, "Mean: %f\n", mean);
-    psTrace("psModules.imcombine", 1, "RMS deviation: %f\n", rms);
-
-    kernels->mean = mean;
-    kernels->rms = rms;
-    kernels->numStamps = numStamps;
-
-    psLogMsg("psModules.imcombine", PS_LOG_INFO, "Mean deviation from %d stamps: %lf +/- %lf",
-             numStamps, mean, rms);
-
-    if (!isfinite(sigmaRej) || sigmaRej <= 0.0) {
-        // User just wanted to calculate and record the deviation for posterity
-        return 0;
-    }
-
-    float limit = sigmaRej * rms; // Limit on maximum deviation
-    psTrace("psModules.imcombine", 1, "Deviation limit: %f\n", limit);
-
+    kernels->mean = stats->sampleMean;
+    kernels->rms = stats->sampleStdev;
+    kernels->numStamps = stats->clippedNvalues;
+
+    psLogMsg ("pmPSFtry", 4, "chisq vs flux fit: %f +/- %f\n", stats->sampleMean, stats->sampleStdev);
+    psLogMsg("psModules.imcombine", PS_LOG_INFO, "Mean deviation from %d stamps: %lf +/- %lf",  kernels->numStamps, kernels->mean, kernels->rms);
 
     psString ds9name = NULL;            // Filename for ds9 region file
@@ -914 +981 @@
     int numRejected = 0;                // Number of stamps rejected
     int numGood = 0;                    // Number of good stamps
-    double newMean = 0.0;               // New mean
     psString log = NULL;                // Log message
-    psStringAppend(&log, "Rejecting stamps, mean = %f, threshold = %f\n", mean, limit);
+
+    // save DS9 region files for the stamps and mark for rejection and replacement
     for (int i = 0; i < stamps->num; i++) {
         pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest
-        if (stamp->status == PM_SUBTRACTION_STAMP_USED) {
+        if (stamp->status  != PM_SUBTRACTION_STAMP_USED) { continue; }
+        if (match->stampMask->data.PS_TYPE_VECTOR_MASK_DATA[i]) {
             // Should we reject stars with low deviation?  Well, if this is really a Gaussian-like
             // distribution and they're low, then we have the right to ask why.  Isn't it suspicious that
@@ -926 +994 @@
             // subtract well, in which case very few (if any) stars will be legitimately rejected for being
             // low.
-            if (fabsf(deviations->data.F32[i] - mean) > limit) {
-                // Mask out the stamp in the image so you it's not found again
-                psTrace("psModules.imcombine", 3, "Rejecting stamp %d (%d,%d)\n", i,
-                        (int)(stamp->x - 0.5), (int)(stamp->y - 0.5));
-                psStringAppend(&log, "Stamp %d (%d,%d): %f\n", i,
-                               (int)(stamp->x - 0.5), (int)(stamp->y - 0.5),
-                               fabsf(deviations->data.F32[i] - mean));
-                numRejected++;
-                for (int y = stamp->y - footprint; y <= stamp->y + footprint; y++) {
-                    for (int x = stamp->x - footprint; x <= stamp->x + footprint; x++) {
-                        subMask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] |= PM_SUBTRACTION_MASK_REJ;
-                    }
-                }
-                pmSubtractionStampPrint(ds9, stamp->x, stamp->y, stamps->footprint, "red");
-
-                // Set stamp for replacement
-                stamp->x = 0;
-                stamp->y = 0;
-                stamp->xNorm = NAN;
-                stamp->yNorm = NAN;
-                stamp->status = PM_SUBTRACTION_STAMP_REJECTED;
-                // Recalculate convolutions
-                psFree(stamp->convolutions1);
-                psFree(stamp->convolutions2);
-                stamp->convolutions1 = stamp->convolutions2 = NULL;
-                psFree(stamp->image1);
-                psFree(stamp->image2);
-                psFree(stamp->weight);
-                stamp->image1 = stamp->image2 = stamp->weight = NULL;
-                psFree(stamp->matrix);
-                stamp->matrix = NULL;
-                psFree(stamp->vector);
-                stamp->vector = NULL;
-            } else {
-                numGood++;
-                newMean += deviations->data.F32[i];
-                pmSubtractionStampPrint(ds9, stamp->x, stamp->y, stamps->footprint, "green");
-            }
-        }
-    }
-    newMean /= numGood;
+            psTrace("psModules.imcombine", 3, "Rejecting stamp %d (%d,%d)\n", i,
+                    (int)(stamp->x - 0.5), (int)(stamp->y - 0.5));
+            psStringAppend(&log, "Stamp %d (%d,%d): %f : %f : %f\n", 
+                           i, (int)(stamp->x - 0.5), (int)(stamp->y - 0.5),
+                           match->chisq->data.F32[i], match->fluxes->data.F32[i], match->chisq->data.F32[i] - psPolynomial1DEval(model, match->fluxes->data.F32[i])); 
+            numRejected++;
+            for (int y = stamp->y - footprint; y <= stamp->y + footprint; y++) {
+                for (int x = stamp->x - footprint; x <= stamp->x + footprint; x++) {
+                    subMask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] |= PM_SUBTRACTION_MASK_REJ;
+                }
+            }
+            pmSubtractionStampPrint(ds9, stamp->x, stamp->y, stamps->footprint, "red");
+
+            // Set stamp for replacement
+            stamp->x = 0;
+            stamp->y = 0;
+            stamp->xNorm = NAN;
+            stamp->yNorm = NAN;
+            stamp->status = PM_SUBTRACTION_STAMP_REJECTED;
+            // Recalculate convolutions
+            psFree(stamp->convolutions1);
+            psFree(stamp->convolutions2);
+            stamp->convolutions1 = stamp->convolutions2 = NULL;
+            psFree(stamp->image1);
+            psFree(stamp->image2);
+            psFree(stamp->weight);
+            stamp->image1 = stamp->image2 = stamp->weight = NULL;
+            psFree(stamp->matrix);
+            stamp->matrix = NULL;
+            psFree(stamp->vector);
+            stamp->vector = NULL;
+        } else {
+            numGood++;
+            pmSubtractionStampPrint(ds9, stamp->x, stamp->y, stamps->footprint, "green");
+        }
+    }
 
     if (numRejected == 0) {
@@ -978 +1041 @@
     }
 
+    psFree(model);
+    psFree(stats);
+
     if (numRejected > 0) {
-        psLogMsg("psModules.imcombine", PS_LOG_INFO,
-                 "%d good stamps; %d rejected.\nMean deviation: %lf --> %lf\n",
-                 numGood, numRejected, mean, newMean);
+        psLogMsg("psModules.imcombine", PS_LOG_INFO, "%d good stamps; %d rejected.\n", numGood, numRejected);
     } else {
-        psLogMsg("psModules.imcombine", PS_LOG_INFO,
-                 "%d good stamps; 0 rejected.\nMean deviation: %lf\n",
-                 numGood, mean);
+        psLogMsg("psModules.imcombine", PS_LOG_INFO, "%d good stamps; 0 rejected.\n", numGood);
     }
 
@@ -1479 +1541 @@
   return true;
 }
+
+static void pmSubtractionQualityFree(pmSubtractionQuality *quality) {
+
+    psFree (quality->fluxes);
+    psFree (quality->chisq);
+    psFree (quality->moments);
+    psFree (quality->stampMask);
+}    
+
+pmSubtractionQuality *pmSubtractionQualityAlloc() {
+
+    pmSubtractionQuality *quality = psAlloc(sizeof(pmSubtractionQuality)); // Stamp list to return
+    psMemSetDeallocator(quality, (psFreeFunc)pmSubtractionQualityFree);
+
+    quality->fluxes = NULL;
+    quality->chisq = NULL;
+    quality->moments = NULL;
+    quality->stampMask = NULL;
+
+    quality->score = NAN;
+    quality->mode = PM_SUBTRACTION_MODE_ERR;
+    quality->spatialOrder = -1;
+    quality->nGood = 0;
+    
+    return quality;
+}