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

Timestamp:

Mar 21, 2007, 9:30:37 PM (19 years ago)

Author:

eugene

Message:

update to new psAstromWCS APIs; drop old files, old functions from psastroUtils.c

File:

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

trunk/psastro/src/psastroUtils.c

@@ -79 +79 @@
 bool psastroUpdateChipToFPA (pmFPA *fpa, pmChip *chip, psArray *rawstars, psArray *refstars) {
 
-    psRegion *region = pmChipExtent (chip);
-    region->x1 -= region->x0;
-    region->y1 -= region->y0;
-    region->x0 = 0;
-    region->y0 = 0;
+    psRegion *region = pmChipPixels (chip);
+
     psFree (chip->fromFPA);
-    chip->fromFPA = psPlaneTransformInvert (NULL, chip->toFPA, *region, 20);
+    chip->fromFPA = psPlaneTransformInvert (NULL, chip->toFPA, *region, 50);
     psFree (region);
 
@@ -134 +131 @@
     }
     return true;
-}
-
-bool psastroNormFPA (pmFPA *fpa, psMetadata *config) {
-
-    // save the raw astrometry for later reference
-    pmFPA *raw = pmFPACopyAstrom (fpa);
-
-    // first pass: measure the per-chip solutions, modify the chip.toFPA terms
-    for (int i = 0; i < fpa->chips->n; i++) {
-        pmChip *chip = fpa->chips->data[i];
-        psastroChipAstrom (chip, config);
-    }
-
-    // second stage: re-normalize the chip terms, passing the
-    // average rotation and offset values to the fpa.toSky
-    if (RENORM) {
-
-        // this code is needed for the mosastro stage, with multiple chip solutions
-
-        double dX, dY, dT, dN;
-        dX = dY = dT = dN = 0;
-
-        psPlane origin, P1, P2;
-        origin.x = 0;
-        origin.y = 0;
-
-        // calculate the average rotation and boresite offset relative to raw
-        for (int i = 0; i < fpa->chips->n; i++) {
-            pmChip *iChip = raw->chips->data[i];
-            pmChip *oChip = fpa->chips->data[i];
-
-            // offset of chip
-            psCoordChipToFPA (&P1, &origin, iChip);
-            psCoordChipToFPA (&P2, &origin, oChip);
-            dX += (P2.x - P1.x);
-            dY += (P2.y - P1.y);
-
-            // get parity-independent rotations for old and new solutions
-            double T1 = psPlaneTransformGetRotation (iChip->toFPA);
-            double T2 = psPlaneTransformGetRotation (oChip->toFPA);
-            dT += T2 - T1;
-            dN ++;
-        }
-
-        dT /= dN;
-        dX /= dN;
-        dY /= dN;
-
-        // R(T)
-        double PC1_1 = fpa->toTPA->x->coeff[1][0][0][0];
-        double PC1_2 = fpa->toTPA->x->coeff[0][1][0][0];
-        double PC2_1 = fpa->toTPA->y->coeff[1][0][0][0];
-        double PC2_2 = fpa->toTPA->y->coeff[0][1][0][0];
-
-        // R(dT)
-        double dPC1_1 = +cos (dT);
-        double dPC1_2 = +sin (dT);
-        double dPC2_1 = -sin (dT);
-        double dPC2_2 = +cos (dT);
-
-        // R'(T) = R(T) * R(dT)
-        double pc1_1 = PC1_1*dPC1_1 + PC1_2*dPC2_1;
-        double pc1_2 = PC1_1*dPC1_2 + PC1_2*dPC2_2;
-        double pc2_1 = PC2_1*dPC1_1 + PC2_2*dPC2_1;
-        double pc2_2 = PC2_1*dPC1_2 + PC2_2*dPC2_2;
-
-        double det = 1.0 / (pc1_1*pc2_2 - pc1_2*pc2_1);
-
-        // R'(-T)  (matrix inverse, not just rotation inverse -- keeps parity)
-        double pi1_1 = +pc2_2 * det;
-        double pi1_2 = -pc1_2 * det;
-        double pi2_1 = -pc2_1 * det;
-        double pi2_2 = +pc1_1 * det;
-
-        // apply the new modifcations in rotation and boresite
-        for (int i = 0; i < fpa->chips->n; i++) {
-            pmChip *oChip = fpa->chips->data[i];
-
-            // r(T)
-            double pr1_1 = oChip->toFPA->x->coeff[1][0];
-            double pr1_2 = oChip->toFPA->x->coeff[0][1];
-            double pr2_1 = oChip->toFPA->y->coeff[1][0];
-            double pr2_2 = oChip->toFPA->y->coeff[0][1];
-
-            // ri'(T) = R(T) r(t)
-            double ri1_1 = PC1_1*pr1_1 + PC1_2*pr2_1;
-            double ri1_2 = PC1_1*pr1_2 + PC1_2*pr2_2;
-            double ri2_1 = PC2_1*pr1_1 + PC2_2*pr2_1;
-            double ri2_2 = PC2_1*pr1_2 + PC2_2*pr2_2;
-
-            // r'(T) = R'(-T) ri'(T)
-            oChip->toFPA->x->coeff[1][0] = pi1_1*ri1_1 + pi1_2*ri2_1;
-            oChip->toFPA->x->coeff[0][1] = pi1_1*ri1_2 + pi1_2*ri2_2;
-            oChip->toFPA->y->coeff[1][0] = pi2_1*ri1_1 + pi2_2*ri2_1;
-            oChip->toFPA->y->coeff[0][1] = pi2_1*ri1_2 + pi2_2*ri2_2;
-
-            double dx = PC1_1*oChip->toFPA->x->coeff[0][0] + PC1_2*oChip->toFPA->y->coeff[0][0] + dX;
-            double dy = PC2_1*oChip->toFPA->x->coeff[0][0] + PC2_2*oChip->toFPA->y->coeff[0][0] + dY;
-
-            oChip->toFPA->x->coeff[0][0] = pi1_1*dx + pi1_2*dy;
-            oChip->toFPA->y->coeff[0][0] = pi2_1*dx + pi2_2*dy;
-        }
-
-        fpa->toTPA->x->coeff[0][0][0][0] -= dX;
-        fpa->toTPA->y->coeff[0][0][0][0] -= dY;
-
-        fpa->toTPA->x->coeff[1][0][0][0] = pc1_1;
-        fpa->toTPA->x->coeff[0][1][0][0] = pc1_2;
-        fpa->toTPA->y->coeff[1][0][0][0] = pc2_1;
-        fpa->toTPA->y->coeff[0][1][0][0] = pc2_2;
-    }
-    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);
 }
 
@@ -328 +181 @@
 }
 
-// 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) {