Changeset 39686

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

Location:

branches/czw_branch/20160809/psastro/src

Files:

• psastroAstromGuess.c (modified) (3 diffs)
• psastroDemoDump.c (modified) (1 diff)
• psastroModelAdjust.c (modified) (3 diffs)
• psastroMosaicAstrom.c (modified) (3 diffs)
• psastroMosaicGradients.c (modified) (1 diff)
• psastroMosaicOneChip.c (modified) (1 diff)
• psastroOneChipFit.c (modified) (1 diff)

branches/czw_branch/20160809/psastro/src/psastroAstromGuess.c

@@ -75 +75 @@
     pmFPA *fpa = input->fpa;
 
-    if (DEBUG) psastroDumpCorners ("corners.up.guess1.dat", "corners.dn.guess1.dat", fpa);
+    // this call only works if we have loaded a model : seg fault if not
+    if (DEBUG && useModel) psastroDumpCorners ("corners.up.guess1.dat", "corners.dn.guess1.dat", fpa);
 
     // load mosaic-level astrometry?
@@ -91 +92 @@
             if (!psastroAstromGuessSetChip (fpa, chip, view, pixelScale, bilevelAstrometry)) continue;
         }
+
+        if (!chip->toFPA || !chip->fromFPA) {
+          char *name = psMetadataLookupStr (NULL, chip->concepts, "CHIP.NAME");
+          fprintf (stderr, "no astrom model for %s, skipping\n", name);
+          continue;
+        }
 
         if (newFPA) {
@@ -312 +319 @@
     while ((chip = pmFPAviewNextChip (view, fpa, 1)) != NULL) {
         if (!chip->process || !chip->file_exists || !chip->data_exists) { continue; }
+
+        if (!chip->toFPA || !chip->fromFPA) {
+          char *name = psMetadataLookupStr (NULL, chip->concepts, "CHIP.NAME");
+          fprintf (stderr, "no astrom model for %s, skipping\n", name);
+          continue;
+        }
 
         // XXX we are currently inconsistent with marking the good vs the bad data

branches/czw_branch/20160809/psastro/src/psastroDemoDump.c

@@ -234 +234 @@
         if (!chip->process || !chip->file_exists || !chip->data_exists) { continue; }
 
+        if (!chip->toFPA || !chip->fromFPA) {
+          char *name = psMetadataLookupStr (NULL, chip->concepts, "CHIP.NAME");
+          fprintf (stderr, "no astrom model for %s, skipping\n", name);
+          continue;
+        }
+
         // XXX write out the four corners for a test
         psRegion *region = pmChipPixels (chip);

branches/czw_branch/20160809/psastro/src/psastroModelAdjust.c

@@ -52 +52 @@
         return false; 
     } 
+    float refChipAngleNominal = PS_RAD_DEG*psMetadataLookupF32 (&status, recipe, "PSASTRO.MODEL.REF.CHIP.ANGLE");
+    if (!refChipName) {
+        psError(PS_ERR_IO, true, "reference chip is missing from recipe"); 
+        return false; 
+    } 
 
     // get reference chip from name
@@ -87 +92 @@
     // get the current posangle of the ref chip
     float chipAngle = atan2 (refChip->toFPA->y->coeff[1][0], refChip->toFPA->x->coeff[1][0]);
-    fprintf (stderr, "chipAngle: %f\n", chipAngle*PS_DEG_RAD);
+
+    // chipAngle should be refChipAngleNominal @ POSANGLE = 0.0
+    float posAngleOffset = chipAngle - refChipAngleNominal;
+
+    fprintf (stderr, "chipAngle is: %f, at PA = 0.0, it should be %f, rotating model by %f\n", 
+             chipAngle*PS_DEG_RAD, refChipAngleNominal*PS_DEG_RAD, posAngleOffset*PS_DEG_RAD);
+
     // psMetadataAddF32 (output->fpa->concepts, PS_LIST_TAIL, "FPA.POSANGLE", PS_META_REPLACE, "boresite parameter", posangle);
 
@@ -99 +110 @@
         psRegion *region = pmChipPixels (chip);
 
-        psPlaneTransform *toFPA = psPlaneTransformRotate (NULL, chip->toFPA, chipAngle);
+        psPlaneTransform *toFPA = psPlaneTransformRotate (NULL, chip->toFPA, posAngleOffset);
         psFree (chip->toFPA);
         chip->toFPA = toFPA;

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;
+}
+

branches/czw_branch/20160809/psastro/src/psastroMosaicGradients.c

@@ -79 +79 @@
     int order = psMetadataLookupF32 (&status, recipe, "PSASTRO.MOSAIC.ORDER");
     if (!status) {
-        psError(PSASTRO_ERR_UNKNOWN, false, "failed to find single-chip fit order\n");
+        psError(PSASTRO_ERR_UNKNOWN, false, "failed to find mosaic distortion fit order\n");
         psFree (gradients);
         psFree (view);

branches/czw_branch/20160809/psastro/src/psastroMosaicOneChip.c

@@ -100 +100 @@
     }
 
-    // XXX allow statitic to be set by the user
+    // XXX allow statistic to be set by the user
     // only clip if we are fitting the chip parameters.
     psStats *fitStats = NULL;
-//    if (order == 0) {
-//      fitStats = psStatsAlloc (PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV);
-//      fitStats->clipSigma = psMetadataLookupF32 (&status, recipe, "PSASTRO.MOSAIC.CHIP.NSIGMA");
-//      fitStats->clipIter = psMetadataLookupS32 (&status, recipe, "PSASTRO.MOSAIC.CHIP.NITER");
-//    } else {
-        fitStats = psStatsAlloc (PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV);
-        fitStats->clipSigma = psMetadataLookupF32 (&status, recipe, "PSASTRO.MOSAIC.CHIP.NSIGMA");
-        fitStats->clipIter = psMetadataLookupS32 (&status, recipe, "PSASTRO.MOSAIC.CHIP.NITER");
-//    }
+    if (FALSE && (order == 0)) {
+      fitStats = psStatsAlloc (PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV);
+      fitStats->clipSigma = psMetadataLookupF32 (&status, recipe, "PSASTRO.MOSAIC.CHIP.NSIGMA");
+      fitStats->clipIter = psMetadataLookupS32 (&status, recipe, "PSASTRO.MOSAIC.CHIP.NITER");
+    } else {
+      fitStats = psStatsAlloc (PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV);
+      fitStats->clipSigma = psMetadataLookupF32 (&status, recipe, "PSASTRO.MOSAIC.CHIP.NSIGMA");
+      fitStats->clipIter = psMetadataLookupS32 (&status, recipe, "PSASTRO.MOSAIC.CHIP.NITER");
+    }
 
     // need to pass in an update header, sent in from above
-    pmAstromFitResults *results = pmAstromMatchFit (chip->toFPA, rawstars, refstars, match, fitStats);
+    pmAstromFitResults *results = pmAstromMatchFit (chip->toFPA, rawstars, refstars, match, fitStats, recipe);
     if (!results) {
         psError(PSASTRO_ERR_DATA, false, "failed to perform the matched fit\n");

branches/czw_branch/20160809/psastro/src/psastroOneChipFit.c

@@ -124 +124 @@
 
         // improved fit for astrometric terms
-        results = pmAstromMatchFit (chip->toFPA, rawstars, refstars, match, fitStats);
+        results = pmAstromMatchFit (chip->toFPA, rawstars, refstars, match, fitStats, recipe);
         if (!results) {
             psLogMsg ("psastro", 3, "failed to perform the matched fit\n");