Changeset 23239 for branches/cnb_branches/cnb_branch_20090215/psModules/src/objects/pmSourceMatch.c

Timestamp:

Mar 9, 2009, 3:41:26 PM (17 years ago)

Author:

beaumont

Message:

sync with trunk

Location:

branches/cnb_branches/cnb_branch_20090215

Files:

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

branches/cnb_branches/cnb_branch_20090215/psModules/src/objects/pmSourceMatch.c

@@ -152 +152 @@
 
 
-psArray *pmSourceMatchSources(const psArray *sourceArrays, float radius)
+psArray *pmSourceMatchSources(const psArray *sourceArrays, float radius, bool cullSingles)
 {
     PS_ASSERT_ARRAY_NON_NULL(sourceArrays, NULL);
@@ -220 +220 @@
 
         // Match with the master list
-
-        psTree *tree = psTreePlant(2, SOURCES_MAX_LEAF, xMaster, yMaster); // kd Tree with sources
-        long numMatch = 0;              // Number of matches
-
-        long size = numMaster + numSources; // New size
-        xMaster = psVectorRealloc(xMaster, size);
-        yMaster = psVectorRealloc(yMaster, size);
-        matches = psArrayRealloc(matches, size);
-
-        psVector *coords = psVectorAlloc(2, PS_TYPE_F32); // Coordinates for tree lookup
-        for (int j = 0; j < numSources; j++) {
-            coords->data.F32[0] = xImage->data.F32[j];
-            coords->data.F32[1] = yImage->data.F32[j];
-            long index = psTreeNearestWithin(tree, coords, radius); // Match index
-            if (index >= 0) {
-                // Record the match
-                pmSourceMatch *match = matches->data[index]; // Match data
-                pmSourceMatchAdd(match, magImage->data.F32[j], magErrImage->data.F32[j], i, j);
-                numMatch++;
-            } else {
-                // Add to the master list
-                pmSourceMatch *match = pmSourceMatchAlloc(numImages); // Match data
-                pmSourceMatchAdd(match, magImage->data.F32[j], magErrImage->data.F32[j], i, j);
-                xMaster->data.F32[numMaster] = xImage->data.F32[j];
-                yMaster->data.F32[numMaster] = yImage->data.F32[j];
-                matches->data[numMaster] = match;
-                numMaster++;
-                xMaster->n = yMaster->n = matches->n = numMaster;
-            }
-
-        }
-        psFree(coords);
-        psFree(tree);
-
+        {
+            psTree *tree = psTreePlant(2, SOURCES_MAX_LEAF, xMaster, yMaster); // kd Tree with sources
+            long numMatch = 0;              // Number of matches
+
+            long size = numMaster + numSources; // New size
+            xMaster = psVectorRealloc(xMaster, size);
+            yMaster = psVectorRealloc(yMaster, size);
+            matches = psArrayRealloc(matches, size);
+
+            psVector *coords = psVectorAlloc(2, PS_TYPE_F32); // Coordinates for tree lookup
+            for (int j = 0; j < numSources; j++) {
+                coords->data.F32[0] = xImage->data.F32[j];
+                coords->data.F32[1] = yImage->data.F32[j];
+                long index = psTreeNearestWithin(tree, coords, radius); // Match index
+                if (index >= 0) {
+                    // Record the match
+                    pmSourceMatch *match = matches->data[index]; // Match data
+                    pmSourceMatchAdd(match, magImage->data.F32[j], magErrImage->data.F32[j], i, j);
+                    numMatch++;
+                } else {
+                    // Add to the master list
+                    pmSourceMatch *match = pmSourceMatchAlloc(numImages); // Match data
+                    pmSourceMatchAdd(match, magImage->data.F32[j], magErrImage->data.F32[j], i, j);
+                    xMaster->data.F32[numMaster] = xImage->data.F32[j];
+                    yMaster->data.F32[numMaster] = yImage->data.F32[j];
+                    matches->data[numMaster] = match;
+                    numMaster++;
+                    xMaster->n = yMaster->n = matches->n = numMaster;
+                }
+
+            }
+            psFree(coords);
+            psFree(tree);
+        }
 
     match_image_done:
@@ -263 +263 @@
     }
 
-    // Now cull the matches that contain only a single star
-    int numGood = 0;                    // Number of good matches
+    if (cullSingles) {
+        // Now cull the matches that contain only a single star
+        int numGood = 0;                    // Number of good matches
+        for (int i = 0; i < matches->n; i++) {
+            pmSourceMatch *match = matches->data[i]; // Match of interest
+            if (match->num > 1) {
+                if (i != numGood) {
+                    psFree(matches->data[numGood]);
+                    matches->data[numGood] = psMemIncrRefCounter(match);
+                }
+                numGood++;
+            }
+        }
+        matches->n = numGood;
+        for (int i = numGood; i < numMaster; i++) {
+            psFree(matches->data[i]);
+            matches->data[i] = NULL;
+        }
+    }
+
+    return matches;
+}
+
+
+psArray *pmSourceMatchMerge(psArray *sourceArrays, float radius)
+{
+    PS_ASSERT_ARRAY_NON_NULL(sourceArrays, NULL);
+    PS_ASSERT_FLOAT_LARGER_THAN(radius, 0.0, NULL);
+
+    psArray *matches = pmSourceMatchSources(sourceArrays, radius, false); // Source matches
+    if (!matches) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to match source lists.");
+        return NULL;
+    }
+
+    int index = 0;                      // Index to current position on list
     for (int i = 0; i < matches->n; i++) {
+        pmSource *source = NULL;        // Source to put in merged list
         pmSourceMatch *match = matches->data[i]; // Match of interest
-        if (match->num > 1) {
-            if (i != numGood) {
-                psFree(matches->data[numGood]);
-                matches->data[numGood] = psMemIncrRefCounter(match);
-                //                psFree(match);
-            }
-            numGood++;
-        }
-    }
-    matches->n = numGood;
-    for (int i = numGood; i < numMaster; i++) {
-        psFree(matches->data[i]);
-        matches->data[i] = NULL;
-    }
+        for (int j = 0; j < match->num && !source; j++) {
+            if (!isfinite(match->mag->data.F32[j]) || !isfinite(match->magErr->data.F32[j])) {
+                continue;
+            }
+            int image = match->image->data.S32[j]; // Index of image
+            int index = match->index->data.S32[j]; // Index of source for image
+            psArray *list = sourceArrays->data[image]; // List of interest
+            source = list->data[index];
+            break;
+        }
+
+        if (source) {
+            psFree(matches->data[index]);
+            matches->data[index] = psMemIncrRefCounter(source);
+            index++;
+        }
+    }
+
+    // Clear out the rest of the list
+    int num = index;                    // Number of good sources
+    for (; index < matches->n; index++) {
+        psFree(matches->data[index]);
+        matches->data[index] = NULL;
+    }
+    matches->n = num;
 
     return matches;
 }
-
 
 // Iterate on the star magnitudes and image transparencies
@@ -524 +569 @@
 psVector *pmSourceMatchRelphot(const psArray *matches, // Array of matches
                                const psVector *zp, // Zero points for each image (including airmass term)
-                               int maxIter, // Maximum number of iterations
                                float tol, // Relative tolerance for convergence
+                               int iter1, // Number of iterations for pass 1
+                               float rej1, // Limit on rejection between iterations for pass 1
+                               float sys1, // Systematic error in measurements for pass 1
+                               int iter2, // Number of iterations for pass 2
+                               float rej2, // Limit on rejection between iterations for pass 2
+                               float sys2, // Systematic error in measurements for pass 2
                                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
+                               float photoLevel // Level at which we declare image is photometric
                                )
 {
@@ -539 +587 @@
     PS_ASSERT_FLOAT_LARGER_THAN(transClip, 0.0, NULL);
 
-    sysErr *= sysErr;
+    sys1 *= sys1;
+    sys2 *= sys2;
 
     int numImages = zp->n;              // Number of images
@@ -552 +601 @@
 
     float chi2 = sourceMatchRelphotIterate(trans, stars, badImage, matches, zp,
-                                           photo, sysErr); // chi^2 for solution
+                                           photo, sys1); // chi^2 for solution
     psTrace("psModules.objects", 1, "Initial: chi^2 = %f\n", chi2);
     float lastChi2 = INFINITY;          // chi^2 on last iteration
     float fracRej = INFINITY;        // Fraction of measurements rejected
 
-    // in the first passes, the transparencies are not well deteremined: use high systematic error and the rejection thresholds
-    for (int i = 0; i < 5; i++) {
+    // In the first passes, the transparencies are not well deteremined: use high systematic error and
+    // rejection thresholds
+    for (int i = 0; i < iter1; i++) {
 
         // Identify photometric nights
@@ -572 +622 @@
         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));
+        fracRej = sourceMatchRelphotReject(trans, stars, matches, zp, photo, badImage, rej1, sys1);
         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));
+        chi2 = sourceMatchRelphotIterate(trans, stars, badImage, matches, zp, photo, sys1);
         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++) {
+    for (int i = 0; i < iter2 && (fabsf(lastChi2 - chi2) > tol * chi2 || fracRej > rejLimit); i++) {
         lastChi2 = chi2;
 
@@ -594 +642 @@
             return NULL;
         }
-        psTrace("psModules.objects", 3, "Determined %d/%d are photometric", numPhoto, numImages);
-
-        fracRej = sourceMatchRelphotReject(trans, stars, matches, zp, photo, badImage, starClip, sysErr);
-        psTrace("psModules.objects", 3, "%f%% of measurements rejected", fracRej * 100);
-
-        chi2 = sourceMatchRelphotIterate(trans, stars, badImage, matches, zp, photo, sysErr);
-        psTrace("psModules.objects", 1, "iter = %d: chi^2 = %f rejected = %f\n", i, chi2, fracRej);
+        psTrace("psModules.objects", 3, "Pass 2: Determined %d/%d are photometric", numPhoto, numImages);
+
+        fracRej = sourceMatchRelphotReject(trans, stars, matches, zp, photo, badImage, rej2, sys2);
+        psTrace("psModules.objects", 3, "Pass 2: %f%% of measurements rejected", fracRej * 100);
+
+        chi2 = sourceMatchRelphotIterate(trans, stars, badImage, matches, zp, photo, sys2);
+        psTrace("psModules.objects", 1, "Pass 2: iter = %d: chi^2 = %f rejected = %f\n", i, chi2, fracRej);
     }