Changeset 12536

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

Location:

trunk/psastro/src

Files:

• psastroAstromGuess.c (modified) (2 diffs)
• psastroLoadRefstars.c (modified) (1 diff)
• psastroMosaicAstrom.c (modified) (3 diffs)
• psastroMosaicRescaleChips.c (deleted)
• psastroUtils.c (modified) (3 diffs)
• psastroWCS.c (deleted)

trunk/psastro/src/psastroAstromGuess.c

@@ -47 +47 @@
     pmFPA *fpa = input->fpa;
 
+    // load mosaic-level astrometry?
+    bool bilevelAstrometry = false;
+    pmHDU *phu = pmFPAviewThisPHU (view, fpa);
+    if (phu) {
+      char *ctype = psMetadataLookupStr (NULL, phu->header, "CTYPE1");
+      if (ctype) {
+        bilevelAstrometry = !strcmp (&ctype[4], "-DIS");
+      }
+    }
+    if (bilevelAstrometry) {
+      pmAstromReadBilevelMosaic (fpa, phu->header);
+    } 
+
     while ((chip = pmFPAviewNextChip (view, fpa, 1)) != NULL) {
         psTrace ("psastro", 4, "Chip %d: %x %x\n", view->chip, chip->file_exists, chip->process);
@@ -53 +66 @@
         // read WCS data from the corresponding header
         pmHDU *hdu = pmFPAviewThisHDU (view, fpa);
+        if (bilevelAstrometry) {
+          pmAstromReadBilevelChip (chip, hdu->header); 
+        } else {
+          pmAstromReadWCS (fpa, chip, hdu->header, pixelScale);
+        }
 
-        pmAstromReadWCS (fpa, chip, hdu->header, pixelScale);
         if (newFPA) {
             newFPA = false;

trunk/psastro/src/psastroLoadRefstars.c

@@ -83 +83 @@
 
     unlink (tempFile);
+
+    if (table == NULL) {
+        psError(PSASTRO_ERR_REFSTARS, true, "failure to load astrometric reference\n");
+        return NULL;
+    }
+
     psLogMsg ("psastro", 3, "read getstar output table : %f sec\n", psTimerMark ("psastro"));
 

trunk/psastro/src/psastroMosaicAstrom.c

@@ -4 +4 @@
 // XXX require this fpa to have multiple chip extensions and a PHU?
 bool psastroMosaicAstrom (pmConfig *config, psArray *refs) {
-
-    bool status;
 
     // select the current recipe
@@ -13 +11 @@
         return false;
     }
-
-    // physical pixel scale in microns per pixel
-    double pixelScale = psMetadataLookupF32 (&status, recipe, "PSASTRO.PIXEL.SCALE");
-    if (!status) {
-        psError(PS_ERR_IO, false, "Failed to lookup pixel scale"); 
-        return false; 
-    } 
 
     // select the input data sources
@@ -101 +92 @@
     if (psTraceGetLevel("psastro.dump") > 0) { psastroDumpMatches (fpa, "match.6.dat"); }
 
-    // save WCS and analysis metadata in update header XXX this is still wrong: the pmFPAExtent
-    // function returns the FP dimensions in pixel units relative to the 0,0 corner of chip
-    // 0,0.  I need to have a function which loops over all cells, determines the pixel
-    // dimensions for each cell, converts them to chip pixels, then uses the fpa astrometry
-    // structures to get the total dimensions of the fpa in fpa units.  equivalent to
-    // pmFPAExtent
-    psRegion *region = pmFPAExtent (fpa);
-    region->x0 *= pixelScale;
-    region->x1 *= pixelScale;
-    region->y0 *= pixelScale;
-    region->y1 *= pixelScale;
-
-    // XXX also need to add DATE/TIME info and NAXIS1, NAXIS2
+    // save WCS and analysis metadata in update header.
     psMetadata *updates = psMetadataAlloc();
-    psMetadataAddS32 (updates, PS_LIST_TAIL, "NAXIS1",   PS_META_REPLACE, "fpa dimensions (mm)", region->x1 - region->x0);
-    psMetadataAddS32 (updates, PS_LIST_TAIL, "NAXIS2",   PS_META_REPLACE, "fpa dimensions (mm)", region->y1 - region->y0);
-    psFree (region);
-    
     if (!pmAstromWriteBilevelMosaic (updates, fpa, NONLIN_TOL)) {
         psAbort ("failed to save header terms");

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) {