Changeset 39686 for branches/czw_branch/20160809/psastro/src/psastroMosaicAstrom.c

Timestamp:

Sep 11, 2016, 2:33:18 PM (10 years ago)

Author:

eugene

Message:

add skips for chips without a valid solution (do not just crash); add recipe value for PSASTRO.MODEL.REF.CHIP.ANGLE, the position angle expected for the reference chip for an exposure with PA = 0.0; posAngle needs to be adjusted to match reference chip target value; if the projection is not TAN or DIS when attempting bilevel mosaic fitting, re-fit the projection + distortion to be TAN + distortion

File:

• branches/czw_branch/20160809/psastro/src/psastroMosaicAstrom.c (modified) (3 diffs)

branches/czw_branch/20160809/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?
@@ -42 +43 @@
     int nIter = psMetadataLookupS32 (&status, recipe, "PSASTRO.MOSAIC.CHIP.NITER");
     if (!status) psAbort ("missing config value");
+
+    // XXX if projection is not TAN, fit the measured projection here and combine toTPA/fromTPA functions
+    if ((fpa->toSky->type != PS_PROJ_TAN) && (fpa->toSky->type != PS_PROJ_DIS)) {
+      if (!psastroProjectionRefit (fpa, recipe)) psAbort ("failed to refit distortion");
+    }
 
     // this should be in a loop with nIter =
@@ -154 +160 @@
     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
+    fpa->fromTPA = psPlaneTransformInvert (fpa->fromTPA, fpa->toTPA, fitRegion, 100);
+
+    psFree (fpa->toSky);
+    fpa->toSky = toSkyTan;
+
+    psFree (L);
+    psFree (M);
+    psFree (P);
+    psFree (Q);
+
+    psFree (view);
+    return true;
+}
+