Changeset 41285 for trunk

Timestamp:

Feb 24, 2020, 4:01:13 PM (6 years ago)

Author:

eugene

Message:

adding recipe values PSASTRO.REFSTAR.CLUMP.NITER, PSASTRO.REFSTAR.CLUMP.SCALE (defaulting to 150, 3 as previously defined); adding log messages to track timing; adding chip-only iteration-depended fitting order; chip-only fitting needs to use bilevel astrometry if the toTPA transformation has order > 1

Location:

trunk/psastro/src

Files:

• psastro.c (modified) (2 diffs)
• psastroAnalysis.c (modified) (2 diffs)
• psastroChipAstrom.c (modified) (4 diffs)
• psastroChooseRefstars.c (modified) (2 diffs)
• psastroOneChipFit.c (modified) (4 diffs)
• psastroRemoveClumps.c (modified) (4 diffs)

trunk/psastro/src/psastro.c

@@ -51 +51 @@
     }
 
+    psLogMsg("psastro", 3, "TIMEMARK: psastroDataLoad: %f sec\n", psTimerMark ("complete"));
+
     psMetadata *stats = psMetadataAlloc(); // Statistics, for output
     psMetadataAddS32(stats, PS_LIST_TAIL, "QUALITY", 0, "No problems", 0);
@@ -61 +63 @@
         exit(PS_EXIT_SYS_ERROR);
     }
+
+    psLogMsg("psastro", 3, "TIMEMARK: psastroAnalysis: %f sec\n", psTimerMark ("complete"));
 
     // write out the results

trunk/psastro/src/psastroAnalysis.c

@@ -135 +135 @@
         }
     }
+
+    psLogMsg("psastro", 3, "TIMEMARK: psastroChipAstrom: %f sec\n", psTimerMark ("complete"));
+
     if (mosastro) {
         if (!psastroMosaicAstrom (config)) {
@@ -150 +153 @@
         }
     }
+
+    psLogMsg("psastro", 3, "TIMEMARK: psastroMosaicAstrom: %f sec\n", psTimerMark ("complete"));
 
     if (!skipastro) {

trunk/psastro/src/psastroChipAstrom.c

@@ -28 +28 @@
     }
 
+    // In this function, we are fitting the chip->FPA model, keeping the toTPA static.  
+    // If the toTPA model has order > 1, then the WCS terms need to use the BiLevel
+    // model.  It is possible for this conversion to fail.  However, if we are going to try
+    // again with the mosaic model (fitting both), then we need to react to failure differently
+    bool mosastro  = psMetadataLookupBool (&status, config->arguments, "PSASTRO.MOSAIC.MODE");
+    if (!status) {
+        mosastro  = psMetadataLookupBool (&status, recipe, "PSASTRO.MOSAIC.MODE");
+    }
+
     // select the input data sources
     pmFPAfile *input = psMetadataLookupPtr (&status, config->files, "PSASTRO.INPUT");
@@ -39 +48 @@
 
     int numGoodChips = 0;               // Number of chips for which astrometry succeeds
+
+    bool useBilevelWCS = (fpa->toTPA->x->nX > 1) || (fpa->toTPA->x->nY > 1);
 
     while ((chip = pmFPAviewNextChip (view, fpa, 1)) != NULL) {
@@ -145 +156 @@
                 psU64 astrom_chip_val = 1;
                 astrom_chip_val <<= view->chip;
-                psMetadataAddU64 (updates, PS_LIST_TAIL, "ASTROM_CHIPS", PS_META_REPLACE, 
-                    "chips that passed astrometry", astrom_chip_val);
+                psMetadataAddU64 (updates, PS_LIST_TAIL, "ASTROM_CHIPS", PS_META_REPLACE, "chips that passed astrometry", astrom_chip_val);
 
                 // write the elapsed time here; this will be updated in psastroMosaicAstrometry, if called
                 psMetadataAddF32 (updates, PS_LIST_TAIL, "DT_ASTR", PS_META_REPLACE, "elapsed psastro time", psTimerMark ("psastroAnalysis"));
                 
-                fpa->wcsCDkeys = psMetadataLookupBool(&status, recipe , "PSASTRO.WCS.USECDKEYS");
-                pmAstromWriteWCS (updates, fpa, chip, NONLIN_TOL);
+                // if the toTPA distortion model uses higher-order terms, then we need to use BiLevel astrometry 
+                if (useBilevelWCS) {
+                    // create the header keywords to descripe the results
+                    if (!pmAstromWriteBilevelChip (updates, chip, NONLIN_TOL)) {
+                        if (!mosastro) { readout->data_exists = false; } // unless we are going to try again, give up on this chip
+                        psError(PS_ERR_UNKNOWN, false, "invalid solution for %d,%d,%d\n", view->chip, view->cell, view->readout);
+                        psErrorStackPrint(stderr, "failure to generate WCS keywords for one chip\n");
+                        psErrorClear();
+                        continue;
+                    }
+                } else {
+                    fpa->wcsCDkeys = psMetadataLookupBool(&status, recipe , "PSASTRO.WCS.USECDKEYS");
+                    pmAstromWriteWCS (updates, fpa, chip, NONLIN_TOL);
+                }
 
                 if (psTraceGetLevel("psastro.dump") > 0) {
@@ -179 +201 @@
     }
 
+    if (useBilevelWCS) {
+        // save WCS and analysis metadata in update header.
+        // (pull or create local view to entry on readout->analysis)
+        psMetadata *updates = psMetadataLookupMetadata (&status, fpa->analysis, "PSASTRO.HEADER");
+        if (!updates) {
+            updates = psMetadataAlloc ();
+            psMetadataAddMetadata (fpa->analysis, PS_LIST_TAIL, "PSASTRO.HEADER",  PS_META_REPLACE, "psastro header stats", updates);
+            psFree (updates);
+        }
+        if (!pmAstromWriteBilevelMosaic (updates, fpa, NONLIN_TOL)) {
+            psError(psErrorCodeLast(), false, "Failed to save header terms");
+            return false;
+        }
+    }
+
     if (fpa->chips->n == 1 && numGoodChips == 0) {
         psError(PSASTRO_ERR_UNKNOWN, false, "Failed to fit single chip.");

trunk/psastro/src/psastroChooseRefstars.c

@@ -46 +46 @@
 
     bool matchLumFunc = psMetadataLookupBool (&status, recipe, "PSASTRO.MATCH.LUMFUNC");
+
+    int nIter = psMetadataLookupS32 (&status, recipe, "PSASTRO.REFSTAR.CLUMP.NITER");
+    if (!status) nIter = 3;
+
+    psF32 clumpScale = psMetadataLookupS32 (&status, recipe, "PSASTRO.REFSTAR.CLUMP.SCALE");
+    if (!status) clumpScale = 150;
 
     pmFPAview *view = pmFPAviewAlloc (0);
@@ -178 +184 @@
                 // generate a reduced subset excluding the clumps
                 // XXX do we need both REFSTARS and SUBSET? 
-                psArray *subset = psastroRemoveClumpsIterate(refstars, 150, 3);
+                psArray *subset = psastroRemoveClumpsIterate(refstars, clumpScale, nIter);
                 psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSASTRO.REFSTARS.SUBSET", PS_DATA_ARRAY, "astrometry objects", subset);
 

trunk/psastro/src/psastroOneChipFit.c

@@ -36 +36 @@
 
     // select the desired chip order
-    REQUIRED_RECIPE_VALUE (int order, "PSASTRO.CHIP.ORDER", S32);
+    REQUIRED_RECIPE_VALUE (int defaultOrder, "PSASTRO.CHIP.ORDER", S32);
 
     // allowed limits for valid solutions
@@ -57 +57 @@
         }
 
+        sprintf (name, "PSASTRO.ONE.CHIP.ORDER.N%d", iter);
+        int order = psMetadataLookupS32 (&status, recipe, name);
+        if (!status) {
+            order = defaultOrder;
+        }
 
         // use small radius to match stars
@@ -76 +81 @@
             match = unique;
         }
+
+        // XXX check if we correctly applied the new transformation:
+        if (psTraceGetLevel("psastro.dump") > 0) {
+          char *filename = NULL;
+          char *chipname = psMetadataLookupStr (&status, chip->concepts, "CHIP.NAME");
+          psStringAppend (&filename, "match.pref.%s.%d.dat", chipname, iter);
+          psastroDumpMatchedStars (filename, rawstars, refstars, match);
+          psFree (filename);
+          filename = NULL;
+        }
+
 
         // modify the order to correspond to the actual number of matched stars:
@@ -137 +153 @@
         psastroUpdateChipToFPA (fpa, chip); // updates PSASTRO.RAWSTARS and PSASTRO.REFSTARS
 
+        // XXX check if we correctly applied the new transformation:
+        if (psTraceGetLevel("psastro.dump") > 0) {
+          char *filename = NULL;
+          char *chipname = psMetadataLookupStr (&status, chip->concepts, "CHIP.NAME");
+          psStringAppend (&filename, "match.post.%s.%d.dat", chipname, iter);
+          psastroDumpMatchedStars (filename, rawstars, refstars, match);
+          psFree (filename);
+          filename = NULL;
+        }
+
         // toSky converts from FPA & TPA units (microns) to sky units (radians)
         float plateScale = 0.5*(fpa->toSky->Xs + fpa->toSky->Ys)*3600.0*PM_DEG_RAD;

trunk/psastro/src/psastroRemoveClumps.c

@@ -35 +35 @@
     }
 
+    int nIter = psMetadataLookupS32 (&status, recipe, "PSASTRO.REFSTAR.CLUMP.NITER");
+    if (!status) nIter = 3;
+
+    psF32 clumpScale = psMetadataLookupS32 (&status, recipe, "PSASTRO.REFSTAR.CLUMP.SCALE");
+    if (!status) clumpScale = 150;
+
     pmFPAview *view = pmFPAviewAlloc (0);
     pmFPA *fpa = input->fpa;
@@ -58 +64 @@
                 // XXX do we need both RAWSTARS and SUBSET? 
                 // XXX put these parameters in the recipe, please
-                psArray *subset = psastroRemoveClumpsIterate(rawstars, 150, 3);
+                psArray *subset = psastroRemoveClumpsIterate(rawstars, clumpScale, nIter);
                 psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSASTRO.RAWSTARS.SUBSET", PS_DATA_ARRAY, "astrometry objects", subset);
                 psFree (subset);
@@ -65 +71 @@
                 if ((gridstars == rawstars) || (gridstars == NULL)) { continue; }
 
-                psArray *gridstars_subset = psastroRemoveClumpsIterate(gridstars, 150, 3);
+                psArray *gridstars_subset = psastroRemoveClumpsIterate(gridstars, clumpScale, nIter);
                 psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSASTRO.GRID.RAWSTARS.SUBSET", PS_DATA_ARRAY, "astrometry objects", gridstars_subset);
                 psFree (gridstars_subset);
@@ -79 +85 @@
     psArray *newset = psMemIncrRefCounter (input);
     for (int i = 0; i < nIter; i++) {
-        psArray *subset = psastroRemoveClumps (newset, 150);
+        psArray *subset = psastroRemoveClumps (newset, scale);
         psFree (newset);
         newset = subset;