Changeset 28386

Timestamp:

Jun 17, 2010, 8:27:02 AM (16 years ago)

Author:

Paul Price

Message:

Remove 'CDi_jX' WCS keywords

File:

• trunk/psModules/src/astrom/pmAstrometryWCS.c (modified) (10 diffs)

trunk/psModules/src/astrom/pmAstrometryWCS.c

@@ -378 +378 @@
     // XXX make it optional to write out CDi_j terms, or other versions
     // apply CDELT1,2 (degrees / pixel) to yield PCi,j terms of order unity
-    if (!wcs->wcsCDkeys) { 
+    if (!wcs->wcsCDkeys) {
 
       double cdelt1 = wcs->cdelt1;
@@ -384 +384 @@
       psMetadataAddF64 (header, PS_LIST_TAIL, "CDELT1", PS_META_REPLACE, "", cdelt1);
       psMetadataAddF64 (header, PS_LIST_TAIL, "CDELT2", PS_META_REPLACE, "", cdelt2);
-      
+
       // test the PC00i00j varient:
       psMetadataAddF64 (header, PS_LIST_TAIL, "PC001001", PS_META_REPLACE, "", wcs->trans->x->coeff[1][0] / cdelt1); // == PC1_1
@@ -390 +390 @@
       psMetadataAddF64 (header, PS_LIST_TAIL, "PC002001", PS_META_REPLACE, "", wcs->trans->y->coeff[1][0] / cdelt1); // == PC2_1
       psMetadataAddF64 (header, PS_LIST_TAIL, "PC002002", PS_META_REPLACE, "", wcs->trans->y->coeff[0][1] / cdelt2); // == PC2_2
-      
+
       // Elixir-style polynomial terms
       // XXX currently, Elixir/DVO cannot accept mixed orders
@@ -398 +398 @@
       if (fitOrder > 1) {
         for (int i = 0; i <= fitOrder; i++) {
-          for (int j = 0; j <= fitOrder; j++) {
-            if (i + j < 2)
-              continue;
-            if (i + j > fitOrder)
-              continue;
-            sprintf (name, "PCA1X%1dY%1d", i, j);
-            psMetadataAddF64 (header, PS_LIST_TAIL, name, PS_META_REPLACE, "", wcs->trans->x->coeff[i][j] / pow(cdelt1, i) / pow(cdelt2, j));
-            sprintf (name, "PCA2X%1dY%1d", i, j);
-            psMetadataAddF64 (header, PS_LIST_TAIL, name, PS_META_REPLACE, "", wcs->trans->y->coeff[i][j] / pow(cdelt1, i) / pow(cdelt2, j));
-          }
+          for (int j = 0; j <= fitOrder; j++) {
+            if (i + j < 2)
+              continue;
+            if (i + j > fitOrder)
+              continue;
+            sprintf (name, "PCA1X%1dY%1d", i, j);
+            psMetadataAddF64 (header, PS_LIST_TAIL, name, PS_META_REPLACE, "", wcs->trans->x->coeff[i][j] / pow(cdelt1, i) / pow(cdelt2, j));
+            sprintf (name, "PCA2X%1dY%1d", i, j);
+            psMetadataAddF64 (header, PS_LIST_TAIL, name, PS_META_REPLACE, "", wcs->trans->y->coeff[i][j] / pow(cdelt1, i) / pow(cdelt2, j));
+          }
         }
         psMetadataAddS32 (header, PS_LIST_TAIL, "NPLYTERM", PS_META_REPLACE, "", fitOrder);
       }
-      
+
       // remove any existing 'CDi_j style' wcs keywords
       if (psMetadataLookup(header, "CD1_1")) {
@@ -419 +419 @@
         psMetadataRemoveKey(header, "CD2_2");
       }
+
+      // Remove 'CDi_jX' WCS keywords
+      psString cd11 = psStringCopy("CD1_1 ");
+      psString cd12 = psStringCopy("CD1_2 ");
+      psString cd21 = psStringCopy("CD2_1 ");
+      psString cd22 = psStringCopy("CD2_2 ");
+      for (char extra = 'A'; extra <= 'Z'; extra++) {
+          cd11[strlen(cd11)-1] = extra;
+          if (psMetadataLookup(header, cd11)) {
+              cd12[strlen(cd12)-1] = extra;
+              cd21[strlen(cd21)-1] = extra;
+              cd22[strlen(cd22)-1] = extra;
+              psMetadataRemoveKey(header, cd11);
+              psMetadataRemoveKey(header, cd12);
+              psMetadataRemoveKey(header, cd21);
+              psMetadataRemoveKey(header, cd22);
+          }
+      }
+      psFree(cd11);
+      psFree(cd12);
+      psFree(cd21);
+      psFree(cd22);
+
+
     } else {
 
@@ -427 +451 @@
 
       if (psMetadataLookup(header, "PC001001")) {
-        psMetadataRemoveKey(header, "PC001001");
-        psMetadataRemoveKey(header, "PC001002");
-        psMetadataRemoveKey(header, "PC002001");
-        psMetadataRemoveKey(header, "PC002002");
+        psMetadataRemoveKey(header, "PC001001");
+        psMetadataRemoveKey(header, "PC001002");
+        psMetadataRemoveKey(header, "PC002001");
+        psMetadataRemoveKey(header, "PC002002");
       }
     }
@@ -685 +709 @@
     }
 
-    // generate transform with the original orientation (does this rotate about 'center'?) 
+    // generate transform with the original orientation (does this rotate about 'center'?)
     psPlaneTransform *tpa2 = psPlaneTransformRotate (NULL, tpa1, -1.0*angle);
 
@@ -967 +991 @@
 
     // outFPA projection must be defined as the goal
-    
+
     // the output transformations are:
     // chip -> FPA : standard linear trans with needed rotation, etc
@@ -990 +1014 @@
         for (int i =  0; i < nSamples; i++) {
 
-            psSphere srcSky;
-            psPlane *srcChip = psPlaneAlloc();
-            psPlane *dstTP = psPlaneAlloc();
+            psSphere srcSky;
+            psPlane *srcChip = psPlaneAlloc();
+            psPlane *dstTP = psPlaneAlloc();
 
             srcChip->x = bounds->x0 + (i * deltaX / nSamples);
             srcChip->y = y;
 
-            psPlaneTransformApply (&srcFP, inChip->toFPA, srcChip);
-            psPlaneTransformApply (&srcTP, inFPA->toTPA, &srcFP);
-            psDeproject (&srcSky, &srcTP, inFPA->toSky);
-            
-            // fprintf (stderr, "%f %f | %f %f | %f %f | %f %f\n", srcChip->x, srcChip->y, srcFP.x, srcFP.y, srcTP.x, srcTP.y, srcSky.r*PS_DEG_RAD, srcSky.d*PS_DEG_RAD);
-
-            psProject (dstTP, &srcSky, outFPA->toSky);
+            psPlaneTransformApply (&srcFP, inChip->toFPA, srcChip);
+            psPlaneTransformApply (&srcTP, inFPA->toTPA, &srcFP);
+            psDeproject (&srcSky, &srcTP, inFPA->toSky);
+
+            // fprintf (stderr, "%f %f | %f %f | %f %f | %f %f\n", srcChip->x, srcChip->y, srcFP.x, srcFP.y, srcTP.x, srcTP.y, srcSky.r*PS_DEG_RAD, srcSky.d*PS_DEG_RAD);
+
+            psProject (dstTP, &srcSky, outFPA->toSky);
 
             srcChip->x -= bounds->x0;
             srcChip->y -= bounds->y0;
-            psArrayAdd (src, 100, srcChip);
-            psArrayAdd (dst, 100, dstTP);
+            psArrayAdd (src, 100, srcChip);
+            psArrayAdd (dst, 100, dstTP);
 
             psFree(srcChip);  // drop our refs to s and d
@@ -1021 +1045 @@
     if (!psPlaneTransformFit(newToFPA, src, dst, 0, 0)) {
         psError(PS_ERR_UNKNOWN, false, "linear fit to transform failed");
-        psFree(src);
-        psFree(dst);
+        psFree(src);
+        psFree(dst);
         return NULL;
     }
-    
+
 # if (0)
     for (int i = 0; i < src->n; i++) {
-        
-        psSphere srcSky, dstSky;
-        psPlane *srcChip = src->data[i];
-        psPlane *dstTP   = dst->data[i];
-
-        psPlaneTransformApply (&srcFP, newToFPA, srcChip);
-        psDeproject (&srcSky, &srcFP, outFPA->toSky);
-        psDeproject (&dstSky, dstTP, outFPA->toSky);
-
-        double dX = (srcSky.r*PS_DEG_RAD - dstSky.r*PS_DEG_RAD)*3600.0;
-        double dY = (srcSky.d*PS_DEG_RAD - dstSky.d*PS_DEG_RAD)*3600.0;
-        fprintf (stderr, "%f %f | %f %f | %f %f | %f %f | %f %f | %f %f\n", dX, dY, srcChip->x, srcChip->y, srcFP.x, srcFP.y, dstTP->x, dstTP->y, srcSky.r*PS_DEG_RAD, srcSky.d*PS_DEG_RAD, dstSky.r*PS_DEG_RAD, dstSky.d*PS_DEG_RAD);
+
+        psSphere srcSky, dstSky;
+        psPlane *srcChip = src->data[i];
+        psPlane *dstTP   = dst->data[i];
+
+        psPlaneTransformApply (&srcFP, newToFPA, srcChip);
+        psDeproject (&srcSky, &srcFP, outFPA->toSky);
+        psDeproject (&dstSky, dstTP, outFPA->toSky);
+
+        double dX = (srcSky.r*PS_DEG_RAD - dstSky.r*PS_DEG_RAD)*3600.0;
+        double dY = (srcSky.d*PS_DEG_RAD - dstSky.d*PS_DEG_RAD)*3600.0;
+        fprintf (stderr, "%f %f | %f %f | %f %f | %f %f | %f %f | %f %f\n", dX, dY, srcChip->x, srcChip->y, srcFP.x, srcFP.y, dstTP->x, dstTP->y, srcSky.r*PS_DEG_RAD, srcSky.d*PS_DEG_RAD, dstSky.r*PS_DEG_RAD, dstSky.d*PS_DEG_RAD);
 
     }