Changeset 39926 for trunk/psastro/src/psastroMosaicAstrom.c

Timestamp:

Jan 6, 2017, 11:30:10 AM (10 years ago)

Author:

eugene

Message:

merging changes from czw dev branch (compare with r39924)

Location:

trunk

Files:

• . (modified) (1 prop)
• psastro/src/psastroMosaicAstrom.c (modified) (3 diffs)

trunk/psastro/src/psastroMosaicAstrom.c

@@ -15 +15 @@
 
 bool psastroMosaicFit (pmFPA *fpa, psMetadata *recipe, const char *rootname, int pass);
+bool psastroProjectionRefit (pmFPA *fpa, psMetadata *recipe);
 
 // XXX require this fpa to have multiple chip extensions and a PHU?
@@ -43 +44 @@
     if (!status) psAbort ("missing config value");
 
+    // if projection is not TAN, fit the measured projection here and combine toTPA/fromTPA functions
+    bool fitMosaicDistortion = true;
+    if ((fpa->toSky->type != PS_PROJ_TAN) && (fpa->toSky->type != PS_PROJ_DIS)) {
+      if (!psastroProjectionRefit (fpa, recipe)) psAbort ("failed to refit distortion");
+      fitMosaicDistortion = false;
+    }
+
     // this should be in a loop with nIter =
-    for (int iter = 0; iter < nIter; iter++) {
+    for (int iter = 0; fitMosaicDistortion && (iter < nIter); iter++) {
         if (!psastroMosaicFit (fpa, recipe, outroot, iter)) return false;
     }
 
-    // now fit the chips under the common distortion with higher-order terms
-    // first, re-perform the match with a slightly tighter circle
-    if (!psastroMosaicSetMatch (fpa, recipe, nIter)) {
-        psError(PSASTRO_ERR_UNKNOWN, false, "failed to match raw and ref stars for mosaic (4th pass)");
-        return false;
-    }
-    if (!psastroMosaicChipAstrom (fpa, recipe, nIter)) {
-        psError(PSASTRO_ERR_UNKNOWN, false, "failed to measure chip astrometry in mosaic mode (4th pass)");
-        return false;
-    }
-
-    if (psTraceGetLevel("psastro.dump") > 0) {
+    for (int iter = 0; (iter < nIter); iter++) {
+      // now fit the chips under the common distortion with higher-order terms
+      // first, re-perform the match with a slightly tighter circle
+      if (!psastroMosaicSetMatch (fpa, recipe, iter)) {
+        psError(PSASTRO_ERR_UNKNOWN, false, "failed to match raw and ref stars for mosaic (4th pass %d)", iter);
+        return false;
+      }
+      if (!psastroMosaicChipAstrom (fpa, recipe, iter)) {
+        psError(PSASTRO_ERR_UNKNOWN, false, "failed to measure chip astrometry in mosaic mode (pass %d)", iter);
+        return false;
+      }
+      
+      if (psTraceGetLevel("psastro.dump") > 0) {
         // the last filename (see filenames in psastroMosaicFit)
         char filename[256];
         snprintf (filename, 256, "%s.%d.dat", outroot, 2*nIter + 2);
         psastroDumpMatches (fpa, filename);
-    }
-
+      }
+    }
+    
     // save WCS and analysis metadata in update header.
     // (pull or create local view to entry on readout->analysis)
@@ -154 +164 @@
     return true;
 }
+
+// we have a fpa->toSky projection which is NOT of type TAN along with a toTPA which is
+// the Identity transform.  we want to convert this to TAN projection plus a non-identity
+// transformation to take the non-TAN components.
+
+// we are going to generate a set of FP->(x,y) using the current projection + fromTPA
+// transformations along with a set of TP->(x,y) values using the desired TAN projection,
+// then we will fit TP vs FP
+
+bool psastroProjectionRefit (pmFPA *fpa, psMetadata *recipe) {
+
+  bool status;
+
+    // allocate mosaic-level polynomial transformation and set masks needed by DVO
+    int order = psMetadataLookupF32 (&status, recipe, "PSASTRO.MOSAIC.ORDER");
+    if (!status) {
+        psError(PSASTRO_ERR_UNKNOWN, false, "failed to find mosaic distortion fit order\n");
+        return false;
+    }
+
+    pmChip *chip = NULL;
+    pmCell *cell = NULL;
+    pmReadout *readout = NULL;
+    pmFPAview *view = pmFPAviewAlloc (0);
+
+    psVector *L = psVectorAllocEmpty (1000, PS_TYPE_F32);
+    psVector *M = psVectorAllocEmpty (1000, PS_TYPE_F32);
+
+    psVector *P = psVectorAllocEmpty (1000, PS_TYPE_F32);
+    psVector *Q = psVectorAllocEmpty (1000, PS_TYPE_F32);
+
+    psProjection *toSkyTan = psProjectionAlloc (fpa->toSky->R, fpa->toSky->D, fpa->toSky->Xs, fpa->toSky->Ys, PS_PROJ_TAN);
+
+    float xMin = NAN;
+    float xMax = NAN;
+    float yMin = NAN;
+    float yMax = NAN;
+
+    bool firstObject = true;
+
+    // this loop selects the matched stars for all chips
+    while ((chip = pmFPAviewNextChip (view, fpa, 1)) != NULL) {
+        psTrace ("psastro", 4, "Chip %d: %x %x\n", view->chip, chip->file_exists, chip->process);
+        if (!chip->process || !chip->file_exists) continue;
+        
+        while ((cell = pmFPAviewNextCell (view, fpa, 1)) != NULL) {
+            psTrace ("psastro", 4, "Cell %d: %x %x\n", view->cell, cell->file_exists, cell->process);
+            if (!cell->process || !cell->file_exists) continue;
+
+            // process each of the readouts
+            // XXX there can only be one readout per chip, right?
+            while ((readout = pmFPAviewNextReadout (view, fpa, 1)) != NULL) {
+                if (! readout->data_exists) continue;
+
+                // select the raw objects for this readout
+                psArray *refstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.REFSTARS.SUBSET");
+                if (refstars == NULL) continue;
+                psTrace ("psastro", 4, "Trying %ld refstars\n", refstars->n);
+
+                psArray *matches = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.MATCH");
+                if (matches == NULL) continue;
+
+                // we are looking over the matched refstars only to be sure we are
+                // covering the valid space of the projection + transformation, and not
+                // beyond
+
+                for (int i = 0; i < matches->n; i++) {
+                    pmAstromMatch *match = matches->data[i];
+                    pmAstromObj *ref = refstars->data[match->ref];
+
+                    psPlane fpOld, tpNew, tpOld;
+
+                    psProject (&tpOld, ref->sky, fpa->toSky); // find the focal-plane coord of this RA,DEC coord using the ref chip projection
+                    psPlaneTransformApply (&fpOld, fpa->fromTPA, &tpOld);
+
+                    psProject (&tpNew, ref->sky, toSkyTan); // find the focal-plane coord of this RA,DEC coord using the ref chip projection
+
+                    psVectorAppend (L, fpOld.x);
+                    psVectorAppend (M, fpOld.y);
+
+                    psVectorAppend (P, tpNew.x);
+                    psVectorAppend (Q, tpNew.y);
+
+                    if (firstObject) {
+                      xMin = xMax = fpOld.x;
+                      yMin = yMax = fpOld.y;
+                      firstObject = false;
+                    }
+
+                    xMin = PS_MIN (xMin, fpOld.x);
+                    xMax = PS_MAX (xMax, fpOld.x);
+                    yMin = PS_MIN (yMin, fpOld.y);
+                    yMax = PS_MAX (yMax, fpOld.y);
+                }
+            }
+        }
+    }
+
+    // the original transforms are (1, 1), but we will need higher order to fit the distortions
+    psFree (fpa->toTPA);
+
+    // the original transforms are (1, 1), but we will need higher order to fit the distortions
+    fpa->toTPA = psPlaneTransformAlloc (order, order);
+    for (int i = 0; i <= fpa->toTPA->x->nX; i++) {
+        for (int j = 0; j <= fpa->toTPA->x->nY; j++) {
+            if (i + j > order) {
+                fpa->toTPA->x->coeffMask[i][j] = PS_POLY_MASK_SET;
+                fpa->toTPA->y->coeffMask[i][j] = PS_POLY_MASK_SET;
+            }
+        }
+    }
+
+    psVectorFitPolynomial2D (fpa->toTPA->x, NULL, 0, P, NULL, L, M);
+    psVectorFitPolynomial2D (fpa->toTPA->y, NULL, 0, Q, NULL, L, M);
+    
+    psRegion fitRegion = psRegionSet (xMin, xMax, yMin, yMax);
+
+    // psPlaneTransformInvert will generate a new fromTPA with order to match toTPA
+    psFree (fpa->fromTPA);
+    psPlaneTransform *myPT = psPlaneTransformAlloc(fpa->toTPA->x->nX+4, fpa->toTPA->x->nY+4);
+    fpa->fromTPA = psPlaneTransformInvert (myPT, fpa->toTPA, fitRegion, 100);
+    psFree (myPT);
+
+    psFree (fpa->toSky);
+    fpa->toSky = toSkyTan;
+
+    psFree (L);
+    psFree (M);
+    psFree (P);
+    psFree (Q);
+
+    psFree (view);
+
+    if (!psastroMosaicSetAstrom (fpa)) {
+        psError(PSASTRO_ERR_UNKNOWN, false, "failed to apply mosaic distortion terms\n");
+        return false;
+    }
+
+    return true;
+}
+