Changeset 7014 for trunk/psastro/src/psastroUtils.c

Timestamp:

Apr 30, 2006, 12:15:03 PM (20 years ago)

Author:

eugene

Message:

first complete working version (chip only)

File:

• trunk/psastro/src/psastroUtils.c (modified) (2 diffs)

trunk/psastro/src/psastroUtils.c

@@ -1 +1 @@
 # include "psastro.h"
 # define RENORM 0
+
+bool psastroUpdateChipToFPA (pmFPA *fpa, pmChip *chip, psArray *rawstars, psArray *refstars) {
+
+    // XXX this region needs to be defined more intelligently
+    psRegion region = {0, 0, 2000, 4500};
+    psFree (chip->fromFPA);
+    chip->fromFPA = psPlaneTransformInvert (NULL, chip->toFPA, region, 20);
+
+    for (int i = 0; i < rawstars->n; i++) {
+        pmAstromObj *raw = rawstars->data[i];
+        
+        psPlaneTransformApply (raw->FP, chip->toFPA, raw->chip);
+        psPlaneDistortApply (raw->TP, fpa->toTangentPlane, raw->FP, 0.0, 0.0);
+        p_psDeproject (raw->sky, raw->TP, fpa->projection);
+    }
+
+    for (int i = 0; i < refstars->n; i++) {
+        pmAstromObj *ref = refstars->data[i];
+        psPlaneTransformApply (ref->chip, chip->fromFPA, ref->FP);
+    }
+
+    return true;
+}
+
+# if 0
 
 bool psastroSelectBrightStars (pmFPA *fpa, psMetadata *config) {
@@ -143 +168 @@
     return true; 
 }
+
+psPolynomial2D *psPolynomial2DCopy (psPolynomial2D *input) {
+
+    psPolynomial2D *output = psPolynomial2DAlloc (input->nX, input->nY, input->type);
+
+    for (int i = 0; i < input->nX; i++) {
+        for (int j = 0; j < input->nY; j++) {
+            output->mask[i][j]     = input->mask[i][j];
+            output->coeff[i][j]    = input->coeff[i][j];
+            output->coeffErr[i][j] = input->coeffErr[i][j];
+        }
+    }
+    return (output);
+}
+
+psPolynomial4D *psPolynomial4DCopy (psPolynomial4D *input) {
+
+    psPolynomial4D *output = psPolynomial4DAlloc (input->nX, input->nY, input->nZ, input->nT, input->type);
+
+    for (int i = 0; i < input->nX; i++) {
+        for (int j = 0; j < input->nY; j++) {
+            for (int k = 0; k < input->nZ; k++) {
+                for (int m = 0; m < input->nT; m++) {
+                    output->mask[i][j][k][m]     = input->mask[i][j][k][m];
+                    output->coeff[i][j][k][m]    = input->coeff[i][j][k][m];
+                    output->coeffErr[i][j][k][m] = input->coeffErr[i][j][k][m];
+                }
+            }
+        }
+    }
+    return (output);
+}
+
+psPlaneDistort *psPlaneDistortCopy (psPlaneDistort *input) {
+
+    psPlaneDistort *output = psPlaneDistortAlloc (input->x->nX, input->x->nY, input->x->nZ, input->x->nT);
+
+    for (int i = 0; i < input->x->nX; i++) {
+        for (int j = 0; j < input->x->nY; j++) {
+            for (int k = 0; k < input->x->nZ; k++) {
+                for (int m = 0; m < input->x->nT; m++) {
+                    // x-terms
+                    output->x->mask[i][j][k][m]     = input->x->mask[i][j][k][m];
+                    output->x->coeff[i][j][k][m]    = input->x->coeff[i][j][k][m];
+                    output->x->coeffErr[i][j][k][m] = input->x->coeffErr[i][j][k][m];
+                    // y-terms
+                    output->y->mask[i][j][k][m]     = input->y->mask[i][j][k][m];
+                    output->y->coeff[i][j][k][m]    = input->y->coeff[i][j][k][m];
+                    output->y->coeffErr[i][j][k][m] = input->y->coeffErr[i][j][k][m];
+                }
+            }
+        }
+    }
+    return (output);
+}
+
+psPlaneTransform *psPlaneTransformCopy (psPlaneTransform *input) {
+
+    psPlaneTransform *output = psPlaneTransformAlloc (input->x->nX, input->x->nY);
+
+    for (int i = 0; i < input->x->nX; i++) {
+        for (int j = 0; j < input->x->nY; j++) {
+            // x-terms
+            output->x->mask[i][j]     = input->x->mask[i][j];
+            output->x->coeff[i][j]    = input->x->coeff[i][j];
+            output->x->coeffErr[i][j] = input->x->coeffErr[i][j];
+            // y-terms
+            output->y->mask[i][j]     = input->y->mask[i][j];
+            output->y->coeff[i][j]    = input->y->coeff[i][j];
+            output->y->coeffErr[i][j] = input->y->coeffErr[i][j];
+        }
+    }
+    return (output);
+}
+
+psProjection *psProjectionCopy (psProjection *input) {
+
+    psProjection *output = psProjectionAlloc (input->R, input->D, input->Xs, input->Ys, input->type);
+    return (output);
+}
+
+// very crude distortion inversion: assumes 0 order in z and t, linear in x and y:
+psPlaneDistort *psPlaneDistortInvert(psPlaneDistort *distort) {
+    PS_ASSERT_PTR_NON_NULL(distort, 0);
+    PS_ASSERT_PTR_NON_NULL(distort->x, 0);
+    PS_ASSERT_PTR_NON_NULL(distort->y, 0);
+
+    psPlaneDistort *out = psPlaneDistortAlloc(1, 1, 0, 0);
+
+    /* simple matrix inversion code */
+
+    psF64 r11 = distort->x->coeff[1][0][0][0];
+    psF64 r12 = distort->x->coeff[0][1][0][0];
+    psF64 r21 = distort->y->coeff[1][0][0][0];
+    psF64 r22 = distort->y->coeff[0][1][0][0];
+    psF64 xo  = distort->x->coeff[0][0][0][0];
+    psF64 yo  = distort->y->coeff[0][0][0][0];
+
+    psF64 invDet = 1.0 / (r11 * r22 - r12 * r21); // Inverse of the determinant
+
+    out->x->coeff[1][0][0][0] = +invDet * r22;
+    out->x->coeff[0][1][0][0] = -invDet * r12;
+    out->y->coeff[1][0][0][0] = -invDet * r21;
+    out->y->coeff[0][1][0][0] = +invDet * r11;
+
+    out->x->coeff[0][0][0][0] = - invDet * (r22 * xo - r12 * yo);
+    out->y->coeff[0][0][0][0] = - invDet * (r11 * yo - r21 * xo);
+
+    return(out);
+}
+
+// returns the rotation term, forcing positive parity
+double psPlaneTransformGetRotation (psPlaneTransform *map) {
+
+    if (map->x->nX < 1) return 0;
+    if (map->x->nY < 1) return 0;
+
+    if (map->y->nX < 1) return 0;
+    if (map->y->nY < 1) return 0;
+    
+    double pc1_1 = map->x->coeff[1][0];
+    double pc1_2 = map->x->coeff[0][1];
+    double pc2_1 = map->y->coeff[1][0];
+    double pc2_2 = map->y->coeff[0][1];
+
+    double px = SIGN (pc1_1);
+    double py = SIGN (pc2_2);
+
+    // both x and y terms imply an angle. take the average
+    double t1 = -atan2 (px*pc1_2, px*pc1_1);
+    double t2 = +atan2 (py*pc2_1, py*pc2_2);
+    
+    // careful near -pi,+pi boundary...
+    if (t1 - t2 > M_PI/2) t2 += 2*M_PI;
+    if (t2 - t1 > M_PI/2) t1 += 2*M_PI;
+
+    double theta = 0.5*(t1 + t2);
+    while (theta < M_PI) theta += 2*M_PI;
+    while (theta > M_PI) theta -= 2*M_PI;
+    
+    return (theta);
+}
+
+# endif