Changeset 23189

Timestamp:

Mar 5, 2009, 9:18:09 AM (17 years ago)

Author:

Paul Price

Message:

Flag an image as bad if its transparency comes out as NAN (e.g., if no
stars matched). This allows the process to finish, and bad images to
be identified.

File:

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

trunk/psModules/src/objects/pmSourceMatch.c

@@ -63 +63 @@
         pmSource *source = sources->data[i]; // Source of interest
         if (!source) continue;
-        if (source->mode & SOURCE_MASK) continue;
-        if (!isfinite(source->psfMag)) continue;
-        if (!isfinite(source->errMag)) continue;
-        if (source->psfMag > SOURCE_FAINTEST) continue;
+        if (source->mode & SOURCE_MASK) continue;
+        if (!isfinite(source->psfMag)) continue;
+        if (!isfinite(source->errMag)) continue;
+        if (source->psfMag > SOURCE_FAINTEST) continue;
 
         float xSrc, ySrc;               // Coordinates of source
@@ -290 +290 @@
 static float sourceMatchRelphotIterate(psVector *trans, // Transparencies
                                        psVector *stars, // Star magnitudes
+                                       psVector *badImage, // Bad image mask
                                        const psArray *matches, // Array of matches
                                        const psVector *zp, // Zero points for each image (incl. airmass term)
@@ -363 +364 @@
     for (int i = 0; i < numImages; i++) {
         trans->data.F32[i] = accum->data.F64[i] / accumErr->data.F64[i];
-
+        if (!isfinite(trans->data.F32[i])) {
+            badImage->data.U8[i] = 0xFF;
+        }
         psTrace("psModules.objects", 3, "Transparency for image %d: %f\n", i, trans->data.F32[i]);
     }
@@ -379 +382 @@
             }
             int index = match->image->data.U32[j]; // Image index
+            if (badImage->data.U8[index]) {
+                continue;
+            }
             float mag = match->mag->data.F32[j]; // Measured magnitude
             float magErr2 = PS_SQR(match->magErr->data.F32[j]) + sysErr2; // Error in measured magnitude
@@ -402 +408 @@
 static int sourceMatchRelphotPhotometric(psVector *photo, // Photometric determination
                                          const psVector *trans, // Estimated transparencies
+                                         const psVector *badImage, // Bad image?
                                          int transIter, // Iterations for transparency
                                          float transClip, // Clipping level for transparency
@@ -409 +416 @@
     psAssert(photo && photo->type.type == PS_TYPE_U8, "Need photometric determination");
     psAssert(trans && trans->type.type == PS_TYPE_F32, "Need transparencies");
+    psAssert(badImage && badImage->type.type == PS_TYPE_U8, "Need bad image determination");
 
     int numImages = photo->n;              // Number of images
 
     psAssert(trans->n == numImages, "Not enough transparencies: %ld", trans->n);
+    psAssert(badImage->n == numImages, "Not enough bad image determinations: %ld", badImage->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);
@@ -421 +430 @@
     stats->clipSigma = transClip;
 
-    if (!psVectorStats(stats, trans, NULL, NULL, 0)) {
+    if (!psVectorStats(stats, trans, NULL, badImage, 0xFF)) {
         psError(PS_ERR_UNKNOWN, false, "Unable to perform statistics on transparencies.");
         psFree(stats);
@@ -432 +441 @@
     int numPhoto = 0;                   // Number of photometric images
     for (int i = 0; i < numImages; i++) {
+        if (badImage->data.U8[i]) {
+            continue;
+        }
         if (trans->data.F32[i] < thresh) {
             photo->data.U8[i] = 0xFF;
@@ -451 +463 @@
                                       const psVector *zp, // Zero points for each image
                                       const psVector *photo, // Photometric image?
+                                      const psVector *badImage, // Bad image?
                                       float starClip, // Clipping for stars
                                       float sysErr2 // Systematic error squared
@@ -468 +481 @@
     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);
+    psAssert(!badImage || badImage->type.type == PS_TYPE_U8, "Photometric determination is wrong type");
+    psAssert(!badImage || badImage->n == numImages, "Not enough bad determinations: %ld", badImage->n);
 
     starClip = PS_SQR(starClip);
@@ -481 +496 @@
             numMeasurements++;
             int index = match->image->data.U32[j]; // Image index
+            if (badImage->data.U8[index]) {
+                continue;
+            }
             float mag = match->mag->data.F32[j]; // Measured magnitude
             float magErr = match->magErr->data.F32[j]; // Error in measured magnitude
@@ -489 +507 @@
             float dev = mag + cal - stars->data.F32[i]; // Deviation
 
-            // only reject detections from photometric images (non-photometric images can
-            // have large errors.  XXX Or: allow a much higher rejection threshold
-            if (photo->data.U8[index]) {
-              if (PS_SQR(dev) > starClip * (PS_SQR(magErr) + sysErr2)) {
-                numRejected++;
-                match->mask->data.PS_TYPE_VECTOR_MASK_DATA[j] = 0xFF;
-              }
-            }
+            // only reject detections from photometric images (non-photometric images can
+            // have large errors.  XXX Or: allow a much higher rejection threshold
+            if (photo->data.U8[index]) {
+                if (PS_SQR(dev) > starClip * (PS_SQR(magErr) + sysErr2)) {
+                    numRejected++;
+                    match->mask->data.PS_TYPE_VECTOR_MASK_DATA[j] = 0xFF;
+                }
+            }
         }
     }
@@ -529 +547 @@
     psVector *photo = psVectorAlloc(numImages, PS_TYPE_U8); // Photometric determination for each image
     psVectorInit(photo, 0);
+    psVector *badImage = psVectorAlloc(numImages, PS_TYPE_U8); // Bad image?
+    psVectorInit(badImage, 0);
     psVector *stars = psVectorAlloc(numStars, PS_TYPE_F32); // Magnitudes for each star
 
-    float chi2 = sourceMatchRelphotIterate(trans, stars, matches, zp, photo, sysErr); // chi^2 for solution
+    float chi2 = sourceMatchRelphotIterate(trans, stars, badImage, matches, zp,
+                                           photo, sysErr); // chi^2 for solution
     psTrace("psModules.objects", 1, "Initial: chi^2 = %f\n", chi2);
     float lastChi2 = INFINITY;          // chi^2 on last iteration
@@ -540 +561 @@
 
         // 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, "Pass 1: Determined %d/%d are photometric", numPhoto, numImages);
-
-        // XXX use 20 sigma rejection and 0.1 mag systematic error (move these to the recipe)
-        fracRej = sourceMatchRelphotReject(trans, stars, matches, zp, photo, 20.0, PS_SQR(0.1));
-        psTrace("psModules.objects", 3, "Pass 1: %f%% of measurements rejected", fracRej * 100);
-
-        // XXX use 0.05 mag systematic error (move these to the recipe)
-        chi2 = sourceMatchRelphotIterate(trans, stars, matches, zp, photo, PS_SQR(0.1)); // chi^2 for solution
-        psTrace("psModules.objects", 1, "Pass 1: iter = %d: chi^2 = %f rejected = %f\n", i, chi2, fracRej);
-    }
-
-    for (int i = 0; i < maxIter && (fabsf(lastChi2 - chi2) > tol * chi2 || fracRej > rejLimit); i++) {
-        lastChi2 = chi2;
-
-        // Identify photometric nights
-        int numPhoto = sourceMatchRelphotPhotometric(photo, trans, transIter, transClip,
+        int numPhoto = sourceMatchRelphotPhotometric(photo, trans, badImage, transIter, transClip,
                                                      photoLevel); // Number of photometric images
         if (numPhoto < 0) {
@@ -572 +570 @@
             return NULL;
         }
+        psTrace("psModules.objects", 3, "Pass 1: Determined %d/%d are photometric", numPhoto, numImages);
+
+        // XXX use 20 sigma rejection and 0.1 mag systematic error (move these to the recipe)
+        fracRej = sourceMatchRelphotReject(trans, stars, matches, zp, photo, badImage, 20.0, PS_SQR(0.1));
+        psTrace("psModules.objects", 3, "Pass 1: %f%% of measurements rejected", fracRej * 100);
+
+        // XXX use 0.05 mag systematic error (move these to the recipe)
+        chi2 = sourceMatchRelphotIterate(trans, stars, badImage, matches, zp, photo, PS_SQR(0.1));
+        psTrace("psModules.objects", 1, "Pass 1: iter = %d: chi^2 = %f rejected = %f\n", i, chi2, fracRej);
+    }
+
+    for (int i = 0; i < maxIter && (fabsf(lastChi2 - chi2) > tol * chi2 || fracRej > rejLimit); i++) {
+        lastChi2 = chi2;
+
+        // Identify photometric nights
+        int numPhoto = sourceMatchRelphotPhotometric(photo, trans, badImage, 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);
+        fracRej = sourceMatchRelphotReject(trans, stars, matches, zp, photo, badImage, starClip, sysErr);
         psTrace("psModules.objects", 3, "%f%% of measurements rejected", fracRej * 100);
 
-        chi2 = sourceMatchRelphotIterate(trans, stars, matches, zp, photo, sysErr); // chi^2 for solution
+        chi2 = sourceMatchRelphotIterate(trans, stars, badImage, matches, zp, photo, sysErr);
         psTrace("psModules.objects", 1, "iter = %d: chi^2 = %f rejected = %f\n", i, chi2, fracRej);
     }
 
     psFree(photo);
+    psFree(badImage);
     psFree(stars);