Changeset 26893 for trunk/psModules/src/astrom/pmAstrometryObjects.c

Timestamp:

Feb 10, 2010, 7:34:39 PM (16 years ago)

Author:

eugene

Message:

updates from eam_branches/20091201 (substantially changes to the kernel matching code; updates to pmDetection container, added pmPattern, etc)

File:

• trunk/psModules/src/astrom/pmAstrometryObjects.c (modified) (16 diffs)

trunk/psModules/src/astrom/pmAstrometryObjects.c

@@ -108 +108 @@
         }
 
-        if (found1->data.S8[i]) {
-            i++;
-            continue;
-        }
-        if (found2->data.S8[j]) {
-            j++;
-            continue;
-        }
+        if (found1->data.S8[i]) {
+            i++;
+            continue;
+        }
+        if (found2->data.S8[j]) {
+            j++;
+            continue;
+        }
 
         jStart = j;
@@ -128 +128 @@
                 continue;
             }
-            if (found2->data.S8[j]) {
-                j++;
-                continue;
-            }
+            if (found2->data.S8[j]) {
+                j++;
+                continue;
+            }
 
             // got a match; add to output list
@@ -138 +138 @@
             psFree (match);
 
-            found1->data.S8[i] = 1;
-            found2->data.S8[j] = 1;
+            found1->data.S8[i] = 1;
+            found2->data.S8[j] = 1;
 
             j++;
@@ -196 +196 @@
     psArray *matches = match_lists(x1, y1, x2, y2, sorted1, sorted2, RADIUS); \
     \
-    psFree((void *)sorted1); \
-    psFree((void *)sorted2); \
+    psFree(sorted1); \
+    psFree(sorted2); \
     psFree(x1); \
     psFree(y1); \
@@ -314 +314 @@
     if (!xRes) abort();
     psFree (xFit);
-    
+
     psVector *yFit = psPolynomial2DEvalVector (map->y, X, Y);
     if (!yFit) abort();
@@ -322 +322 @@
 
     // extract a high-quality subset (unmasked, S/N > XXX) and position errors
-    // XXX for now, generate a position error based on the magnitude error 
+    // XXX for now, generate a position error based on the magnitude error
     psVector *xErr     = psVectorAllocEmpty (match->n, PS_TYPE_F32);
     psVector *yErr     = psVectorAllocEmpty (match->n, PS_TYPE_F32);
@@ -329 +329 @@
 
     for (int i = 0; i < match->n; i++) {
-        if (mask->data.PS_TYPE_VECTOR_MASK_DATA[i]) continue;
+        if (mask->data.PS_TYPE_VECTOR_MASK_DATA[i]) continue;
         pmAstromMatch *pair = match->data[i];
         pmAstromObj *rawStar = raw->data[pair->raw];
@@ -335 +335 @@
         if (rawStar->dMag > 0.02) continue;
 
-        // two likely failure values: NAN or 0.0 --> use dMag in this case
-        float xErrValue, yErrValue;
-        if (isfinite(rawStar->chip->xErr) && (rawStar->chip->xErr > 0.0)) {
-            xErrValue = rawStar->chip->xErr;
-        } else {
-            xErrValue = PS_MAX(0.005, rawStar->dMag);
-        }
-        if (isfinite(rawStar->chip->yErr) && (rawStar->chip->yErr > 0.0)) {
-            yErrValue = rawStar->chip->yErr;
-        } else {
-            yErrValue = PS_MAX(0.005, rawStar->dMag);
-        }
-
-        psVectorAppend (xErr,     xErrValue);
-        psVectorAppend (yErr,     yErrValue);
-        psVectorAppend (xResGood, xRes->data.F32[i]);
-        psVectorAppend (yResGood, yRes->data.F32[i]);
+        // two likely failure values: NAN or 0.0 --> use dMag in this case
+        float xErrValue, yErrValue;
+        if (isfinite(rawStar->chip->xErr) && (rawStar->chip->xErr > 0.0)) {
+            xErrValue = rawStar->chip->xErr;
+        } else {
+            xErrValue = PS_MAX(0.005, rawStar->dMag);
+        }
+        if (isfinite(rawStar->chip->yErr) && (rawStar->chip->yErr > 0.0)) {
+            yErrValue = rawStar->chip->yErr;
+        } else {
+            yErrValue = PS_MAX(0.005, rawStar->dMag);
+        }
+
+        psVectorAppend (xErr,     xErrValue);
+        psVectorAppend (yErr,     yErrValue);
+        psVectorAppend (xResGood, xRes->data.F32[i]);
+        psVectorAppend (yResGood, yRes->data.F32[i]);
     }
 
@@ -434 +434 @@
     // Get the minimum and maximum values
     if (!psVectorStats(stats, myVector, NULL, NULL, 0)) {
-        psFree(stats);
-        return NAN;
+        psFree(stats);
+        return NAN;
     }
     min = stats->min;
     max = stats->max;
     if (isnan(min) || isnan(max)) {
-        psFree(stats);
-        return NAN;
+        psFree(stats);
+        return NAN;
     }
 
@@ -448 +448 @@
     // If all data points have the same value, then we set the appropriate members of stats and return.
     if (fabs(max - min) <= FLT_EPSILON) {
-        psFree (stats);
-        return 0.0;
-    }
-    
-    // Define the histogram bin size.  
+        psFree (stats);
+        return 0.0;
+    }
+
+    // Define the histogram bin size.
     float binSize = 0.001;
     long numBins = PS_MIN(100000, (max - min) / binSize); // Number of bins
@@ -463 +463 @@
 
     if (!psVectorHistogram(histogram, myVector, NULL, NULL, 0)) {
-        // if psVectorHistogram returns false, we have a programming error
-        psAbort ("Unable to generate histogram");
+        // if psVectorHistogram returns false, we have a programming error
+        psAbort ("Unable to generate histogram");
     }
     if (psTraceGetLevel("psModules.astrom") >= 8) {
-        PS_VECTOR_PRINT_F32(histogram->bounds);
-        PS_VECTOR_PRINT_F32(histogram->nums);
+        PS_VECTOR_PRINT_F32(histogram->bounds);
+        PS_VECTOR_PRINT_F32(histogram->nums);
     }
 
@@ -475 +475 @@
     cumulative->nums->data.F32[0] = histogram->nums->data.F32[0];
     for (long i = 1; i < histogram->nums->n; i++) {
-        cumulative->nums->data.F32[i] = cumulative->nums->data.F32[i-1] + histogram->nums->data.F32[i];
-        cumulative->bounds->data.F32[i-1] = histogram->bounds->data.F32[i];
+        cumulative->nums->data.F32[i] = cumulative->nums->data.F32[i-1] + histogram->nums->data.F32[i];
+        cumulative->bounds->data.F32[i-1] = histogram->bounds->data.F32[i];
     }
     if (psTraceGetLevel("psModules.astrom") >= 8) {
-        PS_VECTOR_PRINT_F32(cumulative->bounds);
-        PS_VECTOR_PRINT_F32(cumulative->nums);
+        PS_VECTOR_PRINT_F32(cumulative->bounds);
+        PS_VECTOR_PRINT_F32(cumulative->nums);
     }
 
@@ -837 +837 @@
         }
 
-        if (psTraceGetLevel("psModules.astrom") >= 5) {
-            char line[16];
-            psFits *fits = psFitsOpen ("grid.image.fits", "w");
-            psFitsWriteImage (fits, NULL, gridNP, 0, NULL);
-            psFitsClose (fits);
-            fprintf (stderr, "wrote grid image, press return to continue\n");
-            if (!fgets(line, 15, stdin)) {
-                fprintf(stderr, "Something went wrong; continuing.");
+        if (psTraceGetLevel("psModules.astrom") >= 5) {
+            char line[16];
+            psFits *fits = psFitsOpen ("grid.image.fits", "w");
+            psFitsWriteImage (fits, NULL, gridNP, 0, NULL);
+            psFitsClose (fits);
+            fprintf (stderr, "wrote grid image, press return to continue\n");
+            if (!fgets (line, 15, stdin)) {
+                fprintf(stderr, "Error waiting for RETURN.");
             }
-        }
+        }
 
         // only check bins with at least 1/2 of max bin
@@ -1191 +1191 @@
 psArray *pmAstromRadiusMatchUniq (psArray *rawstars, psArray *refstars, psArray *matches) {
 
-    // I have the matches between the refstars and the rawstars.  
+    // I have the matches between the refstars and the rawstars.
     // For each refstar, find the single match which has the smallest radius and reject
-    // all others.  
+    // all others.
 
     // create an array of refstars->n arrays, each containing all of the matches for the
-    // given refstar.  
+    // given refstar.
 
     psArray *refstarMatches = psArrayAlloc (refstars->n);
@@ -1202 +1202 @@
     for (int i = 0; i < matches->n; i++) {
 
-        pmAstromMatch *match = matches->data[i];
-        
-        // accumulate this refstar match on the array for this refstar (create if needed)
-        psArray *refSet = refstarMatches->data[match->ref];
-        if (!refSet) {
-            refstarMatches->data[match->ref] = psArrayAllocEmpty (8);
-            refSet = refstarMatches->data[match->ref];
-        }
-
-        pmAstromMatchInfo *matchInfo = pmAstromMatchInfoAlloc();
-
-        pmAstromObj *refStar = refstars->data[match->ref];
-        pmAstromObj *rawStar = rawstars->data[match->raw];
-        
-        matchInfo->match = match; // reference to the match of interest
-        matchInfo->radius = hypot (refStar->FP->x - rawStar->FP->x, refStar->FP->y - rawStar->FP->y);
-
-        psArrayAdd (refSet, 8, matchInfo); // matchInfo->match is just a reference
-        psFree (matchInfo);
+        pmAstromMatch *match = matches->data[i];
+
+        // accumulate this refstar match on the array for this refstar (create if needed)
+        psArray *refSet = refstarMatches->data[match->ref];
+        if (!refSet) {
+            refstarMatches->data[match->ref] = psArrayAllocEmpty (8);
+            refSet = refstarMatches->data[match->ref];
+        }
+
+        pmAstromMatchInfo *matchInfo = pmAstromMatchInfoAlloc();
+
+        pmAstromObj *refStar = refstars->data[match->ref];
+        pmAstromObj *rawStar = rawstars->data[match->raw];
+
+        matchInfo->match = match; // reference to the match of interest
+        matchInfo->radius = hypot (refStar->FP->x - rawStar->FP->x, refStar->FP->y - rawStar->FP->y);
+
+        psArrayAdd (refSet, 8, matchInfo); // matchInfo->match is just a reference
+        psFree (matchInfo);
     }
 
@@ -1229 +1229 @@
     for (int i = 0; i < refstars->n; i++) {
 
-        psArray *refSet = refstarMatches->data[i];
-        if (!refSet) continue;
-        if (refSet->n == 0) continue; // not certain how this can happen...
-
-        if (refSet->n == 1) {
-            pmAstromMatchInfo *matchInfo = refSet->data[0];
-            psArrayAdd (unique, 32, matchInfo->match);
-            continue;
-        }
-
-        pmAstromMatchInfo *matchInfo = refSet->data[0];
-        float minRadius = matchInfo->radius;
-        pmAstromMatch *minMatch = matchInfo->match;
-        for (int j = 1; j < refSet->n; j++) {
-            pmAstromMatchInfo *matchInfo = refSet->data[j];
-            if (minRadius < matchInfo->radius) continue;
-            minMatch = matchInfo->match;
-            minRadius = matchInfo->radius;
-        }
-        
-        psArrayAdd (unique, 32, minMatch); // minMatch is just a reference to a match on matches,
+        psArray *refSet = refstarMatches->data[i];
+        if (!refSet) continue;
+        if (refSet->n == 0) continue; // not certain how this can happen...
+
+        if (refSet->n == 1) {
+            pmAstromMatchInfo *matchInfo = refSet->data[0];
+            psArrayAdd (unique, 32, matchInfo->match);
+            continue;
+        }
+
+        pmAstromMatchInfo *matchInfo = refSet->data[0];
+        float minRadius = matchInfo->radius;
+        pmAstromMatch *minMatch = matchInfo->match;
+        for (int j = 1; j < refSet->n; j++) {
+            pmAstromMatchInfo *matchInfo = refSet->data[j];
+            if (minRadius < matchInfo->radius) continue;
+            minMatch = matchInfo->match;
+            minRadius = matchInfo->radius;
+        }
+
+        psArrayAdd (unique, 32, minMatch); // minMatch is just a reference to a match on matches,
     }