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

Timestamp:

Jun 1, 2006, 4:15:35 PM (20 years ago)

Author:

magnier

Message:

finish up the basic astrometry functions

File:

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

trunk/psModules/src/astrom/pmAstrometryObjects.c

@@ -8 +8 @@
 *  @author EAM, IfA
 *
-*  @version $Revision: 1.6 $ $Name: not supported by cvs2svn $
-*  @date $Date: 2006-05-19 15:49:37 $
+*  @version $Revision: 1.7 $ $Name: not supported by cvs2svn $
+*  @date $Date: 2006-06-02 02:15:35 $
 *
 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -34 +34 @@
  ******************************************************************************/
 
+// XXX an example (might be faster)
+void pmAstromObjSortByChipX2 (psArray *objects)
+{
+
+    int l,j,ir,i;
+    float v1, v2;
+    pmAstromObj *tmp;
+    pmAstromObj *obj;
+
+    int N = objects->n;
+
+    if (N < 2)
+        return;
+    l = N >> 1;
+    ir = N - 1;
+    for (;;) {
+        if (l > 0) {
+            tmp = objects->data[--l];
+        } else {
+            tmp = objects->data[ir];
+            objects->data[ir] = objects->data[0];
+            if (--ir == 0) {
+                objects->data[0] = tmp;
+                return;
+            }
+        }
+        i = l;
+        j = (l << 1) + 1;
+        while (j <= ir) {
+            obj = objects->data[j];
+            v1 = obj->chip->x;
+            obj = objects->data[j+1];
+            v2 = obj->chip->x;
+            if (j < ir && v1 < v2)
+                ++j;
+            obj = objects->data[j];
+            v1 = obj->chip->x;
+            v2 = tmp->chip->x;
+            if (v2 < v1) {
+                objects->data[i] = objects->data[j];
+                j += (i=j) + 1;
+            } else
+                j = ir + 1;
+        }
+        objects->data[i] = tmp;
+    }
+}
+
 int pmAstromObjSortByFPX(
     const void **a,
@@ -59 +107 @@
 }
 
-
+int pmAstromObjSortByChipX(
+    const void **a,
+    const void **b)
+{
+    if (PM_ASTROMETRYOBJECTS_DEBUG) {
+        PS_ASSERT_PTR_NON_NULL(a, 0);
+        PS_ASSERT_PTR_NON_NULL(*a, 0);
+        PS_ASSERT_PTR_NON_NULL(b, 0);
+        PS_ASSERT_PTR_NON_NULL(*b, 0);
+    }
+
+    pmAstromObj *A = *(pmAstromObj **)a;
+    pmAstromObj *B = *(pmAstromObj **)b;
+
+    psF32 diff = A->chip->x - B->chip->x;
+    if (diff > FLT_EPSILON) {
+        return (+1);
+    }
+
+    if (diff < FLT_EPSILON) {
+        return (-1);
+    }
+    return (0);
+}
 
 /******************************************************************************
@@ -116 +187 @@
     st2 = psArraySort (st2, pmAstromObjSortByFPX);
 
+    // pmAstromObjSortByChipX2 (st1);
+    // pmAstromObjSortByChipX2 (st2);
+
     psArray *matches = psArrayAlloc (100);
     matches->n = 0;
@@ -124 +198 @@
 
         dX = ((pmAstromObj *)st1->data[i])->FP->x - ((pmAstromObj *)st2->data[j])->FP->x;
-        if (dX < -RADIUS) {
+        if (dX <= -RADIUS) {
             i++;
             continue;
         }
-        if (dX > +RADIUS) {
+        if (dX >= +RADIUS) {
             j++;
             continue;
@@ -134 +208 @@
 
         jStart = j;
-        while ((dX > -RADIUS) && (j < st2->n)) {
+        while ((fabs(dX) < RADIUS) && (j < st2->n)) {
 
             dX = ((pmAstromObj *)st1->data[i])->FP->x - ((pmAstromObj *)st2->data[j])->FP->x;
@@ -155 +229 @@
         i++;
     }
-    psLogMsg (__func__, 3, "radius match: %d pairs\n", matches->n);
+    psLogMsg (__func__, 3, "radius match: %d pairs (radius: %f)\n", matches->n, RADIUS);
+    return (matches);
+}
+
+
+
+/******************************************************************************
+pmAstromRadiusMatchChip(st1, st2, config) : match objects on the chip
+ ******************************************************************************/
+psArray *pmAstromRadiusMatchChip(
+    psArray *st1,
+    psArray *st2,
+    psMetadata *config)
+{
+    PS_ASSERT_PTR_NON_NULL(st1, NULL);
+    PS_ASSERT_PTR_NON_NULL(st2, NULL);
+    PS_ASSERT_PTR_NON_NULL(config, NULL);
+    bool status;
+    int jStart;
+    double dX, dY, dR;
+
+    double RADIUS = psMetadataLookupF32 (&status, config, "PSASTRO.MATCH.RADIUS");
+    double RADIUS_SQR = PS_SQR (RADIUS);
+
+    // sort both lists by Focal Plane X coord
+    st1 = psArraySort (st1, pmAstromObjSortByChipX);
+    st2 = psArraySort (st2, pmAstromObjSortByChipX);
+
+    psArray *matches = psArrayAlloc (100);
+    matches->n = 0;
+
+    int i = 0;
+    int j = 0;
+    while ((i < st1->n) && (j < st2->n)) {
+
+        dX = ((pmAstromObj *)st1->data[i])->chip->x - ((pmAstromObj *)st2->data[j])->chip->x;
+        if (dX <= -RADIUS) {
+            i++;
+            continue;
+        }
+        if (dX >= +RADIUS) {
+            j++;
+            continue;
+        }
+
+        jStart = j;
+        while ((dX > -RADIUS) && (j < st2->n)) {
+
+            dX = ((pmAstromObj *)st1->data[i])->chip->x - ((pmAstromObj *)st2->data[j])->chip->x;
+            dY = ((pmAstromObj *)st1->data[i])->chip->y - ((pmAstromObj *)st2->data[j])->chip->y;
+            dR = dX*dX + dY*dY;
+
+            if (dR > RADIUS_SQR) {
+                j++;
+                continue;
+            }
+
+            // got a match; add to output list
+            pmAstromMatch *match = pmAstromMatchAlloc (i, j);
+            psArrayAdd (matches, 100, match);
+            psFree (match);
+
+            j++;
+        }
+        j = jStart;
+        i++;
+    }
+    psLogMsg (__func__, 3, "radius match: %d pairs (radius: %f)\n", matches->n, RADIUS);
     return (matches);
 }
@@ -235 +376 @@
     // fit chip-to-FPA transformation
     psVectorClipFitPolynomial2D (map->x, stats, xMask, 0, x, wt, X, Y);
+    psLogMsg ("psModules.astrom", 3, "x resid: %f +/- %f (%d of %d)\n", stats->clippedMean, stats->clippedStdev, stats->clippedNvalues, x->n);
+
     psVectorClipFitPolynomial2D (map->y, stats, yMask, 0, y, wt, X, Y);
+    psLogMsg ("psModules.astrom", 3, "y resid: %f +/- %f (%d of %d)\n", stats->clippedMean, stats->clippedStdev, stats->clippedNvalues, y->n);
 
     psFree (x);