Changeset 20953 for trunk/psModules/src/objects/pmSourceMatch.c

Timestamp:

Dec 10, 2008, 5:24:56 PM (17 years ago)

Author:

Paul Price

Message:

Added photometric determination and measurement clipping for relphot. Seems to work.

File:

• trunk/psModules/src/objects/pmSourceMatch.c (modified) (8 diffs)

trunk/psModules/src/objects/pmSourceMatch.c

@@ -63 +63 @@
         pmSource *source = sources->data[i]; // Source of interest
         if (!source || (source->mode & SOURCE_MASK) || !isfinite(source->psfMag) ||
-            source->psfMag > SOURCE_FAINTEST) {
+            !isfinite(source->errMag) || source->psfMag > SOURCE_FAINTEST) {
             continue;
         }
@@ -122 +122 @@
     match->image = psVectorAllocEmpty(num, PS_TYPE_U32);
     match->index = psVectorAllocEmpty(num, PS_TYPE_U32);
+    match->mask = psVectorAllocEmpty(num, PS_TYPE_MASK);
 
     return match;
@@ -142 +143 @@
     match->image->data.S32[num] = image;
     match->index->data.S32[num] = index;
+    match->index->data.PS_TYPE_MASK_DATA[num] = 0;
     match->num++;
 
-    match->mag->n = match->magErr->n = match->image->n = match->index->n = match->num;
+    match->mag->n = match->magErr->n = match->image->n = match->index->n = match->mask->n = match->num;
 }
 
@@ -287 +289 @@
                                 psVector *stars, // Star magnitudes
                                 const psArray *matches, // Array of matches
-                                const psVector *zp // Zero points for each image (including airmass term)
+                                const psVector *zp, // Zero points for each image (including airmass term)
+                                const psVector *photo // Photometric image?
                                 )
 {
+    psAssert(zp && zp->type.type == PS_TYPE_F32, "Need zero points");
+    psAssert(matches, "Need list of matches");
+
     int numImages = zp->n;              // Number of images
     int numStars = matches->n;          // Number of stars
@@ -297 +303 @@
     psAssert(stars && stars->type.type == PS_TYPE_F32, "Need star magnitudes");
     psAssert(stars->n == numStars, "Not enough stars: %ld\n", stars->n);
-    psAssert(matches, "Need list of matches");
-    psAssert(zp && zp->type.type == PS_TYPE_F32, "Need zero points");
     psAssert(zp->n == numImages, "Not enough ZPs: %ld", zp->n);
-
-    psVector *transErr = psVectorAlloc(numImages, PS_TYPE_F32); // Errors in the transparency
-    psVector *starsErr = psVectorAlloc(numStars, PS_TYPE_F32); // Errors in the star magnitudes
+    psAssert(!photo || photo->type.type == PS_TYPE_U8, "Photometric determination is wrong type");
+    psAssert(!photo || photo->n == numImages, "Not enough photometric determinations: %ld", photo->n);
 
     // Solve the star magnitudes
     psVectorInit(stars, 0.0);
-    psVectorInit(starsErr, 0.0);
+    int numGoodStars = 0;               // Number of stars with good measurements
+    for (int i = 0; i < numStars; i++) {
+        pmSourceMatch *match = matches->data[i]; // Matched stars
+        int numMeasurements = 0;        // Number of unmasked measurements for star
+        double star = 0.0, starErr = 0.0; // Accumulators for star
+        for (int j = 0; j < match->num; j++) {
+            if (match->mask->data.PS_TYPE_MASK_DATA[j]) {
+                continue;
+            }
+            numMeasurements++;
+            int index = match->image->data.U32[j]; // Image index
+            float mag = match->mag->data.F32[j]; // Measured magnitude
+            double magErr = match->magErr->data.F32[j]; // Error in measured magnitude
+            double invErr2 = 1.0 / PS_SQR(magErr); // Inverse square error
+            float cal = zp->data.F32[index]; // Calibration to apply to image
+            if (!photo || !photo->data.U8[index]) {
+                cal -= trans->data.F32[index];
+            }
+            star += (mag + cal) * invErr2;
+            starErr += invErr2;
+        }
+        if (numMeasurements > 1) {
+            // It's only a good star (contributing to the chi^2) if there's more than 1 measurement
+            numGoodStars++;
+        }
+        stars->data.F32[i] = star / starErr;
+    }
+
+    // Solve for the transparencies
+    // We solve even for the "photometric" images since they may jump out of that status upon iteration
+    psVector *accum = psVectorAlloc(numImages, PS_TYPE_F64); // Transparency accumulator
+    psVector *accumErr = psVectorAlloc(numImages, PS_TYPE_F64); // Transparency accumulator
+    psVectorInit(accum, 0.0);
+    psVectorInit(accumErr, 0.0);
     for (int i = 0; i < numStars; i++) {
         pmSourceMatch *match = matches->data[i]; // Matched stars
         for (int j = 0; j < match->num; j++) {
+            if (match->mask->data.PS_TYPE_MASK_DATA[j]) {
+                continue;
+            }
+            int index = match->image->data.U32[j]; // Image index
+            float mag = match->mag->data.F32[j]; // Measured magnitude
+            double magErr = match->magErr->data.F32[j]; // Error in measured magnitude
+            double invErr2 = 1.0 / PS_SQR(magErr); // Inverse square error
+            float cal = zp->data.F32[index]; // Calibration to apply to image
+            accum->data.F64[index] += (mag + cal - stars->data.F32[i]) * invErr2;
+            accumErr->data.F64[index] += invErr2;
+        }
+    }
+    for (int i = 0; i < numImages; i++) {
+        trans->data.F32[i] = accum->data.F64[i] / accumErr->data.F64[i];
+
+        psTrace("psModules.objects", 3, "Transparency for image %d: %f\n", i, trans->data.F32[i]);
+    }
+    psFree(accum);
+    psFree(accumErr);
+
+    // Once more through to evaluate chi^2
+    float chi2 = 0.0;                   // chi^2 for iteration
+    int dof = 0;                        // Degrees of freedom
+    for (int i = 0; i < numStars; i++) {
+        pmSourceMatch *match = matches->data[i]; // Matched stars
+        for (int j = 0; j < match->num; j++) {
+            if (match->mask->data.PS_TYPE_MASK_DATA[j]) {
+                continue;
+            }
             int index = match->image->data.U32[j]; // Image index
             float mag = match->mag->data.F32[j]; // Measured magnitude
             float magErr = match->magErr->data.F32[j]; // Error in measured magnitude
-            float invErr2 = 1.0 / PS_SQR(magErr); // Inverse square error
-            stars->data.F32[i] += (mag + zp->data.F32[index] - trans->data.F32[index]) * invErr2;
-            starsErr->data.F32[i] += invErr2;
-        }
-    }
-    for (int i = 0; i < numStars; i++) {
-        float inverse = 1.0 / starsErr->data.F32[i]; // Inverse error
-        stars->data.F32[i] *= inverse;
-        starsErr->data.F32[i] = sqrtf(inverse);
-    }
-
-    // Solve for the transparencies
-    psVectorInit(trans, 0.0);
-    psVectorInit(transErr, 0.0);
+            float cal = zp->data.F32[index]; // Calibration to apply to image
+            if (!photo || !photo->data.U8[index]) {
+                cal -= trans->data.F32[index];
+            }
+            float dev = (mag + cal - stars->data.F32[index]) / magErr; // Deviation from fit
+            if (isfinite(dev)) {
+                chi2 += PS_SQR(dev);
+                dof++;
+            }
+        }
+    }
+    dof -= numGoodStars + numImages;
+    chi2 /= dof;
+
+    return chi2;
+}
+
+// Determine which images are photometric, based on estimated transparencies
+// Returns number of photometric images, or -1 on error
+static int sourceMatchRelphotPhotometric(psVector *photo, // Photometric determination
+                                         const psVector *trans, // Estimated transparencies
+                                         int transIter, // Iterations for transparency
+                                         float transClip, // Clipping level for transparency
+                                         float photoLevel // Level below which we declare photometric
+                                         )
+{
+    psAssert(photo && photo->type.type == PS_TYPE_U8, "Need photometric determination");
+    psAssert(trans && trans->type.type == PS_TYPE_F32, "Need transparencies");
+
+    int numImages = photo->n;              // Number of images
+
+    psAssert(trans->n == numImages, "Not enough transparencies: %ld", trans->n);
+    psAssert(transIter >= 0, "Iterations for transparency must be non-negative: %d", transIter);
+    psAssert(transClip > 0, "Clipping level for transparency must be positive: %f", transClip);
+    psAssert(photoLevel > 0, "Photometric level must be positive: %f", photoLevel);
+
+    psStats *stats = psStatsAlloc(PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV); // Statistics
+    stats->clipIter = transIter;
+    stats->clipSigma = transClip;
+
+    if (!psVectorStats(stats, trans, NULL, NULL, 0)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to perform statistics on transparencies.");
+        psFree(stats);
+        return -1;
+    }
+
+    float thresh = stats->clippedMean + photoLevel * stats->clippedStdev; // Threshold for clouds
+    psFree(stats);
+
+    int numPhoto = 0;                   // Number of photometric images
+    for (int i = 0; i < numImages; i++) {
+        if (trans->data.F32[i] < thresh) {
+            photo->data.U8[i] = 0xFF;
+            numPhoto++;
+        } else {
+            photo->data.U8[i] = 0;
+        }
+    }
+
+    return numPhoto;
+}
+
+
+// Reject star measurements
+// Returns the fraction of measurements that were rejected
+static float sourceMatchRelphotReject(psVector *trans, // Transparencies
+                                      psVector *stars, // Star magnitudes
+                                      const psArray *matches, // Array of matches
+                                      const psVector *zp, // Zero points for each image
+                                      const psVector *photo, // Photometric image?
+                                      float starClip, // Clipping for stars
+                                      float sysErr // Systematic error
+                                )
+{
+    psAssert(zp && zp->type.type == PS_TYPE_F32, "Need zero points");
+    psAssert(matches, "Need list of matches");
+
+    int numImages = zp->n;              // Number of images
+    int numStars = matches->n;          // Number of stars
+
+    psAssert(trans && trans->type.type == PS_TYPE_F32, "Need transparencies");
+    psAssert(trans->n == numImages, "Not enough transparencies: %ld\n", trans->n);
+    psAssert(stars && stars->type.type == PS_TYPE_F32, "Need star magnitudes");
+    psAssert(stars->n == numStars, "Not enough stars: %ld\n", stars->n);
+    psAssert(zp->n == numImages, "Not enough ZPs: %ld", zp->n);
+    psAssert(!photo || photo->type.type == PS_TYPE_U8, "Photometric determination is wrong type");
+    psAssert(!photo || photo->n == numImages, "Not enough photometric determinations: %ld", photo->n);
+
+    starClip = PS_SQR(starClip);
+    sysErr = PS_SQR(sysErr);
+
+    int numRejected = 0;                // Number rejected
+    int numMeasurements = 0;            // Number of measurements
     for (int i = 0; i < numStars; i++) {
         pmSourceMatch *match = matches->data[i]; // Matched stars
         for (int j = 0; j < match->num; j++) {
+            if (match->mask->data.PS_TYPE_MASK_DATA[j]) {
+                continue;
+            }
+            numMeasurements++;
             int index = match->image->data.U32[j]; // Image index
             float mag = match->mag->data.F32[j]; // Measured magnitude
             float magErr = match->magErr->data.F32[j]; // Error in measured magnitude
-            float invErr2 = 1.0 / PS_SQR(magErr); // Inverse square error
-            trans->data.F32[index] += (mag - stars->data.F32[i]) * invErr2;
-            transErr->data.F32[index] += invErr2;
-        }
-    }
-    for (int i = 0; i < numImages; i++) {
-        float inverse = 1.0 / transErr->data.F32[i]; // Inverse error
-        trans->data.F32[i] *= inverse;
-        transErr->data.F32[i] = sqrtf(inverse);
-
-        psTrace("psModules.objects", 3, "Transparency for image %d: %f +/- %f\n",
-                i, trans->data.F32[i], transErr->data.F32[i]);
-    }
-
-    // Once more through to evaluate chi^2
-    float chi2 = 0.0;                   // chi^2 for iteration
-    for (int i = 0; i < numStars; i++) {
-        pmSourceMatch *match = matches->data[i]; // Matched stars
-        for (int j = 0; j < match->num; j++) {
-            int index = match->image->data.U32[j]; // Image index
-            float mag = match->mag->data.F32[j]; // Measured magnitude
-            float magErr = match->magErr->data.F32[j]; // Error in measured magnitude
-            chi2 += PS_SQR(mag - stars->data.F32[index] - trans->data.F32[index] + zp->data.F32[index]) /
-                PS_SQR(magErr);
-        }
-    }
-    chi2 /= numStars + numImages;
-
-    psFree(transErr);
-    psFree(starsErr);
-
-    return chi2;
+            float cal = zp->data.F32[index]; // Calibration to apply to image
+            if (!photo || !photo->data.U8[index]) {
+                cal -= trans->data.F32[index];
+            }
+            float dev = mag + cal - stars->data.F32[i]; // Deviation
+            if (PS_SQR(dev) > starClip * (PS_SQR(magErr) + sysErr)) {
+                numRejected++;
+                match->mask->data.PS_TYPE_MASK_DATA[j] = 0xFF;
+            }
+        }
+    }
+
+    return (float)numRejected / (float)numMeasurements;
 }
 
@@ -372 +501 @@
                                int maxIter, // Maximum number of iterations
                                float tol, // Relative tolerance for convergence
-                               float cloudClip, // Clipping for clouds
-                               float starClip // Clipping for stars
+                               float rejLimit, // Limit on rejection between iterations
+                               int transIter, // Clipping iterations for transparency
+                               float transClip, // Clipping level for transparency
+                               float photoLevel, // Level at which we declare image is photometric
+                               float starClip, // Clipping for stars
+                               float sysErr // Systematic error in measurements
                                )
 {
@@ -379 +512 @@
     PS_ASSERT_VECTOR_NON_NULL(zp, NULL);
     PS_ASSERT_VECTOR_TYPE(zp, PS_TYPE_F32, NULL);
-    PS_ASSERT_FLOAT_LARGER_THAN(cloudClip, 0.0, NULL);
+    PS_ASSERT_FLOAT_LARGER_THAN(transClip, 0.0, NULL);
 
     int numImages = zp->n;              // Number of images
@@ -385 +518 @@
     psVector *trans = psVectorAlloc(numImages, PS_TYPE_F32); // Transparencies for each image, magnitudes
     psVectorInit(trans, 0.0);
+    psVector *photo = psVectorAlloc(numImages, PS_TYPE_U8); // Photometric determination for each image
+    psVectorInit(photo, 0);
     psVector *stars = psVectorAlloc(numStars, PS_TYPE_F32); // Magnitudes for each star
 
-    float chi2 = sourceMatchRelphotIterate(trans, stars, matches, zp); // chi^2 for solution
+    float chi2 = sourceMatchRelphotIterate(trans, stars, matches, zp, photo); // chi^2 for solution
     psTrace("psModules.objects", 1, "Initial: chi^2 = %f\n", chi2);
     float lastChi2 = INFINITY;          // chi^2 on last iteration
-    for (int i = 0; i < maxIter && (lastChi2 - chi2) / chi2 > tol; i++) {
+    float fracRej = INFINITY;        // Fraction of measurements rejected
+    for (int i = 0; i < maxIter && (fabsf(lastChi2 - chi2) > tol * chi2 || fracRej > rejLimit); i++) {
         lastChi2 = chi2;
-        chi2 = sourceMatchRelphotIterate(trans, stars, matches, zp); // chi^2 for solution
-        psTrace("psModules.objects", 1, "iter = %d: chi^2 = %f\n", i, chi2);
+
+        // Identify photometric nights
+        int numPhoto = sourceMatchRelphotPhotometric(photo, trans, transIter, transClip,
+                                                     photoLevel); // Number of photometric images
+        if (numPhoto < 0) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to perform photometric determination");
+            psFree(trans);
+            psFree(photo);
+            psFree(stars);
+            return NULL;
+        }
+        psTrace("psModules.objects", 3, "Determined %d/%d are photometric", numPhoto, numImages);
+
+        fracRej = sourceMatchRelphotReject(trans, stars, matches, zp, photo, starClip, sysErr);
+        psTrace("psModules.objects", 3, "%f%% of measurements rejected", fracRej * 100);
+
+        chi2 = sourceMatchRelphotIterate(trans, stars, matches, zp, photo); // chi^2 for solution
+        psTrace("psModules.objects", 1, "iter = %d: chi^2 = %f rejected = %f\n", i, chi2, fracRej);
+    }
+
+    psFree(photo);
+    psFree(stars);
+
+    if (fabsf(lastChi2 - chi2) > tol * chi2 || fracRej > rejLimit) {
+        psWarning("Unable to converge to relphot solution (%f,%f)", (lastChi2 - chi2) / chi2, fracRej);
     }