Changeset 5560 for trunk/psastro/src/pmAstrom.c

Timestamp:

Nov 21, 2005, 11:03:53 AM (21 years ago)

Author:

eugene

Message:

various cleanups and fixes

File:

• trunk/psastro/src/pmAstrom.c (modified) (7 diffs)

trunk/psastro/src/pmAstrom.c

@@ -1 +1 @@
 # include "pmAstrom.h"
 
+// sort by mag (descending)
+int pmAstromObjSortByFPX (const void **a, const void **b)
+{
+    pmAstromObj *A = *(pmAstromObj **)a;
+    pmAstromObj *B = *(pmAstromObj **)b;
+
+    psF32 diff = A->FP.x - B->FP.x;
+    if (diff > FLT_EPSILON) return (+1);
+    if (diff < FLT_EPSILON) return (-1);
+    return (0);
+}
+
+psPlane *psCoordReadoutToCell (psPlane *cellpix, psPlane *readpix, pmReadout *readout) {
+
+    if (cellpix == NULL) {
+        cellpix = psPlaneAlloc ();
+    }
+
+    cellpix->x = readpix->x*readout->colBins + readout->col0;
+    cellpix->y = readpix->y*readout->rowBins + readout->row0;
+
+    return (cellpix);
+}
+
+psPlane *psCoordCellToReadout (psPlane *readpix, psPlane *cellpix, pmReadout *readout) {
+
+    if (readpix == NULL) {
+        readpix = psPlaneAlloc ();
+    }
+
+    readpix->x = (cellpix->x - readout->col0) / readout->colBins;
+    readpix->y = (cellpix->y - readout->row0) / readout->rowBins;
+
+    return (readpix);
+}
+
+psPlane* psCoordChipToCell_EAM(psPlane* cellCoord,
+                           const psPlane* chipCoord,
+                           const pmCell* cell)
+{
+    PS_ASSERT_PTR_NON_NULL(chipCoord, NULL);
+    PS_ASSERT_PTR_NON_NULL(cell, NULL);
+    PS_ASSERT_PTR_NON_NULL(cell->parent, NULL);
+
+    // XXX EAM : why was this being done?
+    // pmCell *tmpCell = pmCellInChip(chipCoord, cell->parent);
+    PS_ASSERT_PTR_NON_NULL(cell->toChip, NULL);
+    psPlaneTransform *tmpChipToCell = p_psPlaneTransformLinearInvert(cell->toChip);
+    PS_ASSERT_PTR_NON_NULL(tmpChipToCell, NULL);
+    cellCoord = psPlaneTransformApply(cellCoord, tmpChipToCell, chipCoord);
+    psFree(tmpChipToCell);
+    return(cellCoord);
+}
+
+bool psastroProjectFPA (pmFPA *fpa, char *starlist, bool toSky) {
+
+    bool status;
+
+    for (int i = 0; i < fpa->chips->n; i++) {
+        pmChip *chip = fpa->chips->data[i];
+        for (int j = 0; j < chip->cells->n; j++) {
+            pmCell *cell = chip->cells->data[j];
+            for (int k = 0; k < cell->readouts->n; k++) {
+                pmReadout *readout = cell->readouts->data[k];
+                psArray *stars = psMetadataLookupPtr (&status, readout->analysis, starlist);
+                if (toSky) {
+                    psastroProjectRawstars (stars, readout);
+                } else {
+                    psastroProjectRefstars (stars, readout);
+                }
+            }
+        }
+    }
+    return true;
+}
+ 
+bool psastroProjectRawstars (psArray *stars, pmReadout *readout) {
+
+    pmCell *cell = readout->parent;
+    pmChip *chip = cell->parent;
+    pmFPA  *fpa  = chip->parent;
+
+    for (int i = 0; i < stars->n; i++) {
+        pmAstromObj *star = stars->data[i];
+        psCoordReadoutToCell (&star->cell, &star->pix, readout);
+        psCoordCellToChip (&star->chip, &star->cell, cell);
+        psCoordChipToFPA (&star->FP, &star->chip, chip);
+        psCoordFPAToTP (&star->TP, &star->FP, 0.0, 0.0, fpa);
+        psCoordTPToSky (&star->sky, &star->TP, fpa->projection);
+    }
+    return true;
+}
+ 
+bool psastroProjectRefstars (psArray *stars, pmReadout *readout) {
+
+    pmCell *cell = readout->parent;
+    pmChip *chip = cell->parent;
+    pmFPA  *fpa  = chip->parent;
+
+    for (int i = 0; i < stars->n; i++) {
+        pmAstromObj *star = stars->data[i];
+        psCoordSkyToTP (&star->TP, &star->sky, fpa->projection);
+        psCoordTPToFPA (&star->FP, &star->TP, 0.0, 0.0, fpa);
+        psCoordFPAToChip (&star->chip, &star->FP, chip);
+        psCoordChipToCell_EAM (&star->cell, &star->chip, cell);
+        psCoordCellToReadout (&star->pix, &star->cell, readout);
+    }
+    return true;
+}
+ 
 psArray *pmAstromRadiusMatch (psArray *st1, psArray *st2, psMetadata *config) {
 
-    // assume the lists are sorted, or sort in place?
-    // sort st1 by P
-    // sort st2 by P
-
-    double dX, dY;
+    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, pmAstromObjSortByFPX);
+    st2 = psArraySort (st2, pmAstromObjSortByFPX);
+
+    psArray *matches = psArrayAlloc (100);
+    matches->n = 0;
 
     int i = 0;
     int j = 0;
-
-    double RADIUS = psMetadataLookupR32 (&status, myHeader, "ASTROM_MATCH_RADIUS");
-    double RADIUS_SQR = PS_SQR (RADIUS);
-
-    psArray *matches = psArrayAlloc (100);
-    matches->n = 0;
-
     while ((i < st1->n) && (j < st2->n)) {
 
@@ -34 +145 @@
             
             dX = ((pmAstromObj *)st1->data[i])->FP.x - ((pmAstromObj *)st2->data[j])->FP.x;
-            dQ = ((pmAstromObj *)st1->data[i])->FP.y - ((pmAstromObj *)st2->data[j])->FP.y;
-            dR = dX*dX + dQ*dQ;
+            dY = ((pmAstromObj *)st1->data[i])->FP.y - ((pmAstromObj *)st2->data[j])->FP.y;
+            dR = dX*dX + dY*dY;
 
             if (dR > RADIUS_SQR) {
@@ -44 +155 @@
             // got a match; add to output list
             pmAstromMatch *match = pmAstromMatchAlloc (i, j);
-            psArrayAdd (matches, match, 100);
+            psArrayAdd (matches, 100, match);
 
             j++;
@@ -51 +162 @@
         i++;
     }
-    return (match);
+    return (matches);
 }
 
 // take two matched star lists and fit a psPlaneTransform between them
 // 
-psPlaneTransform *pmAstromMatchedListFit (psPlaneTransform *map, psArray *st1, psArray *st2, psArray *match, psMetadata *config) {
+psPlaneTransform *pmAstromMatchFit (psPlaneTransform *map, psArray *raw, psArray *ref, psArray *match, psMetadata *config) {
+
+    bool status;
+    pmAstromObj *rawStar, *refStar;
+    pmAstromMatch *pair;
 
     if (map == NULL) {
-        // int nX = psMetadataLookupS32 (&status, myHeader, "CHIP.NX");
-        // int nY = psMetadataLookupS32 (&status, myHeader, "CHIP.NY");
-        map = psPlaneTransform (2, 2);
-    }
-
-    psStats *stats = psStatsAlloc (CLIPPED_MEAN);
-    stats->nSigma = psMetadataLookupS32 (&status, myHeader, "CHIP.NITER");
-    stats->nSigma = psMetadataLookupS32 (&status, myHeader, "CHIP.NSIGMA");
-
-    psVector *X = psVectorAlloc (match->n);
-    psVector *Y = psVectorAlloc (match->n);
-    psVector *x = psVectorAlloc (match->n);
-    psVector *y = psVectorAlloc (match->n);
+        int nX = psMetadataLookupS32 (&status, config, "PSASTRO.CHIP.NX");
+        if (!status) nX = 1;
+        int nY = psMetadataLookupS32 (&status, config, "PSASTRO.CHIP.NY");
+        if (!status) nY = 1;
+        map = psPlaneTransformAlloc (nX, nY);
+        // XXX EAM : not clear that we are allowed to use orders > 1
+    }
+
+    psStats *stats = psStatsAlloc (PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV);
+    stats->clipSigma = psMetadataLookupF32 (&status, config, "PSASTRO.CHIP.NSIGMA");
+    stats->clipIter  = psMetadataLookupS32 (&status, config, "PSASTRO.CHIP.NITER");
+
+    // XXX EAM : clip fit seems to only work for F32!
+    // XXX EAM : clip fit fails to handle NULL error
+    psVector *X = psVectorAlloc (match->n, PS_TYPE_F32);
+    psVector *Y = psVectorAlloc (match->n, PS_TYPE_F32);
+    psVector *x = psVectorAlloc (match->n, PS_TYPE_F32);
+    psVector *y = psVectorAlloc (match->n, PS_TYPE_F32);
+    psVector *wt = psVectorAlloc (match->n, PS_TYPE_F32);
   
     // take the matched stars, first fit 
     for (int i = 0; i < match->n; i++) {
 
-        pair = match->data[i];
-        ob1 = st1->data[pair->i1];
-        ob2 = st2->data[pair->i2];
-
-        X->data.F64[i] = ob1->chip.x;
-        Y->data.F64[i] = ob1->chip.y;
-
-        x->data.F64[i] = ob2->FP.x;
-        y->data.F64[i] = ob2->FP.y;
+        pair    = match->data[i];
+        rawStar = raw->data[pair->i1];
+        refStar = ref->data[pair->i2];
+
+        X->data.F32[i] = rawStar->chip.x;
+        Y->data.F32[i] = rawStar->chip.y;
+
+        x->data.F32[i] = refStar->FP.x;
+        y->data.F32[i] = refStar->FP.y;
+
+        wt->data.F32[i] = 1.0;
     }
 
     // no masking, no errors
-    psVectorClipFitPolynomial2D (map->x, stats, NULL, 0, X, NULL, x, y);
-    psVectorClipFitPolynomial2D (map->y, stats, NULL, 0, Y, NULL, x, y);
+    psVectorClipFitPolynomial2D (map->x, stats, NULL, 0, X, wt, x, y);
+    psVectorClipFitPolynomial2D (map->y, stats, NULL, 0, Y, wt, x, y);
     psFree (x);
     psFree (y);
     psFree (X);
     psFree (Y);
+    psFree (wt);
     psFree (stats);
 
@@ -103 +227 @@
 psArray *pmAstromRotateObj (psArray *old, psPlane center, double angle) {
 
+    double X, Y;
     pmAstromObj *newObj;
+    pmAstromObj *oldObj;
 
     psArray *new = psArrayAlloc (old->n);
 
-    cs = cos(angle);
-    sn = sin(angle);
-    xCenter = center.x;
-    yCenter = center.y;
+    double cs = cos(angle);
+    double sn = sin(angle);
+    double xCenter = center.x;
+    double yCenter = center.y;
     
     for (int i = 0; i < old->n; i++) {
@@ -120 +246 @@
         Y = oldObj->FP.y - yCenter;
         
-        newObj->FP.x = X*cs + Y*sn;
-        newObj->FP.y = Y*cs - X*sn;
+        newObj->FP.x = X*cs + Y*sn + xCenter;
+        newObj->FP.y = Y*cs - X*sn + yCenter;
         
         new->data[i] = newObj;
@@ -148 +274 @@
   obj->TP.y = 0;
 
-  obj->sky.x = 0;
-  obj->sky.y = 0;
-
-  obj->mag = 0;
+  obj->sky.r = 0;
+  obj->sky.d = 0;
+
+  obj->Mag = 0;
+  obj->dMag = 0;
 
   return (obj);