Changeset 10855

Timestamp:

Dec 29, 2006, 8:36:33 AM (20 years ago)

Author:

eugene

Message:

extensive work to catch errors and clear leaks; implement mosaics

Location:

trunk/psastro/src

Files:

• psastro.c (modified) (1 diff)
• psastroArguments.c (modified) (1 diff)
• psastroAstromGuess.c (modified) (5 diffs)
• psastroChipAstrom.c (modified) (1 diff)
• psastroChooseRefstars.c (modified) (3 diffs)
• psastroConvert.c (modified) (1 diff)
• psastroLuminosityFunction.c (modified) (1 diff)
• psastroOneChip.c (modified) (6 diffs)
• psastroParseCamera.c (modified) (1 diff)
• psastroRefstarSubset.c (modified) (1 diff)
• psastroUtils.c (modified) (2 diffs)

trunk/psastro/src/psastro.c

@@ -53 +53 @@
     }
 
-    if (!psastroChipAstrom (config, refs)) {
-      psErrorStackPrint(stderr, "failed to perform single chip astrometry\n");
-      exit (1);
+    bool chipastro = psMetadataLookupBool (NULL, config->arguments, "PSASTRO.CHIP.MODE");
+    if (chipastro) {
+        if (!psastroChipAstrom (config, refs)) {
+            psErrorStackPrint(stderr, "failed to perform single chip astrometry\n");
+            exit (1);
+        }
     }
 
-    bool mosastro = psMetadataLookupBool (NULL, config->arguments, "MOSAIC.MODE");
-    if (!mosastro) {
+    bool mosastro = psMetadataLookupBool (NULL, config->arguments, "PSASTRO.MOSAIC.MODE");
+    if (mosastro) {
         if (!psastroMosaicAstrom (config, refs)) {
             psErrorStackPrint(stderr, "failed to perform mosaic camera astrometry\n");

trunk/psastro/src/psastroArguments.c

@@ -33 +33 @@
     if ((N = psArgumentGet (argc, argv, "-chip"))) {
         psArgumentRemove (N, &argc, argv);
-        psMetadataAddStr (config->arguments, PS_LIST_TAIL, "CHIP_SELECTIONS", PS_DATA_STRING, "", psStringCopy(argv[N]));
+        psMetadataAddStr (config->arguments, PS_LIST_TAIL, "CHIP_SELECTIONS", PS_META_REPLACE, "", psStringCopy(argv[N]));
         psArgumentRemove (N, &argc, argv);
     }
 
-    // mosastro mode also specifies output header file
+    // apply mosastro mode?
     if ((N = psArgumentGet (argc, argv, "-mosastro"))) {
         psArgumentRemove (N, &argc, argv);
-        psMetadataAddStr (config->arguments, PS_LIST_TAIL, "MOSAIC.MODE", PS_DATA_BOOL, "", psStringCopy(argv[N]));
+        psMetadataAddBool (config->arguments, PS_LIST_TAIL, "PSASTRO.MOSAIC.MODE", PS_META_REPLACE, "", false);
+    }
+    if ((N = psArgumentGet (argc, argv, "+mosastro"))) {
         psArgumentRemove (N, &argc, argv);
+        psMetadataAddBool (config->arguments, PS_LIST_TAIL, "PSASTRO.MOSAIC.MODE", PS_META_REPLACE, "", true);
+    }
+
+    // apply chipastro mode?
+    if ((N = psArgumentGet (argc, argv, "-chipastro"))) {
+        psArgumentRemove (N, &argc, argv);
+        psMetadataAddBool (config->arguments, PS_LIST_TAIL, "PSASTRO.CHIP.MODE", PS_META_REPLACE, "", false);
+    }
+    if ((N = psArgumentGet (argc, argv, "+chipastro"))) {
+        psArgumentRemove (N, &argc, argv);
+        psMetadataAddBool (config->arguments, PS_LIST_TAIL, "PSASTRO.CHIP.MODE", PS_META_REPLACE, "", true);
     }
 

trunk/psastro/src/psastroAstromGuess.c

@@ -11 +11 @@
     bool newFPA = true;
     bool status = false;
-    double RAmin = NAN;
-    double RAmax = NAN;
+    double RAmin = FLT_MAX;
+    double RAmax = FLT_MIN;
+    double DECmin = FLT_MAX;
+    double DECmax = FLT_MIN;
+
     double RAminSky = NAN;
     double RAmaxSky = NAN;
-    double DECmin = NAN;
-    double DECmax = NAN;
 
     pmChip *chip = NULL;
     pmCell *cell = NULL;
     pmReadout *readout = NULL;
-
-    DECmin = DECmax = -90;
-    RAmin = RAmax = RAminSky = RAmaxSky = 0;
 
     // select the current recipe
@@ -63 +61 @@
             RAminSky = fpa->toSky->R - M_PI;
             RAmaxSky = fpa->toSky->R + M_PI;
-            RAmin = RAmax = fpa->toSky->R;
-            DECmin = DECmax = fpa->toSky->D;
         }
 
@@ -79 +75 @@
                 if (rawstars == NULL) { continue; }
 
+                FILE *f = fopen ("rawstars.dat", "w");
                 for (int i = 0; i < rawstars->n; i++) {
                     pmAstromObj *raw = rawstars->data[i];
@@ -85 +82 @@
                     psPlaneTransformApply (raw->TP, fpa->toTPA, raw->FP);
                     psDeproject (raw->sky, raw->TP, fpa->toSky);
+
+                    fprintf (f, "%d  %f %f  %f %f  %f %f  %f %f\n", i,
+                             raw->sky->r, raw->sky->d, 
+                             raw->TP->x, raw->TP->y, 
+                             raw->FP->x, raw->FP->y, 
+                             raw->chip->x, raw->chip->y);
 
                     if (i < 10) {
@@ -122 +125 @@
                     DECmax = PS_MAX (raw->sky->d, DECmax);
                 }
+                fclose (f);
             }
         }

trunk/psastro/src/psastroChipAstrom.c

@@ -49 +49 @@
                 if (refstars == NULL) { continue; }
 
+                // the absolute minimum number of stars is 4 (for order = 1)
+                if ((rawstars->n < 4) || (refstars->n < 4)) {
+                    psLogMsg ("psastro", 3, "insufficient rawstars (%ld) or refstars (%ld)", 
+                            rawstars->n, refstars->n);
+                    continue;
+                } 
+
                 // save WCS and analysis metadata in update header
                 psMetadata *updates = psMetadataAlloc();
 
+                // XXX update the header with info to reflect the failure
                 if (!psastroOneChip (fpa, chip, refstars, rawstars, recipe, updates)) {
-                    psError (PSASTRO_ERR_DATA, false, "failed to find a solution\n");
-                    return false;
+                    psLogMsg ("psastro", 3, "failed to find a solution\n");
+                    psFree (updates);
+                    continue;
                 }
                 pmAstromWriteWCS (updates, fpa, chip, NONLIN_TOL);

trunk/psastro/src/psastroChooseRefstars.c

@@ -56 +56 @@
                 psArray *refstars = psArrayAllocEmpty (100);
 
+                FILE *f = fopen ("refstars.dat", "w");
+
                 // select the reference objects within range of this readout
                 // project the reference objects to this chip
@@ -64 +66 @@
                     psPlaneTransformApply (ref->FP, fpa->fromTPA, ref->TP);
                     psPlaneTransformApply (ref->chip, chip->fromFPA, ref->FP);
+
+                    fprintf (f, "%d  %f %f  %f %f  %f %f  %f %f\n", i,
+                             ref->sky->r, ref->sky->d, 
+                             ref->TP->x, ref->TP->y, 
+                             ref->FP->x, ref->FP->y, 
+                             ref->chip->x, ref->chip->y);
 
                     // limit the X,Y range of the refs to the selected chip
@@ -75 +83 @@
                     psFree (ref);
                 }
+                fclose (f);
                 psTrace ("psastro", 4, "Added %ld refstars\n", refstars->n);
 

trunk/psastro/src/psastroConvert.c

@@ -39 +39 @@
 
     psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSASTRO.RAWSTARS", PS_DATA_ARRAY, "astrometry objects", rawstars);
+    psLogMsg ("psastro", 4, "loaded %ld sources\n", rawstars->n);
+
     psFree (rawstars);
 

trunk/psastro/src/psastroLuminosityFunction.c

@@ -109 +109 @@
   }
 
-  psLogMsg ("psastro", 4, "logN %f + %f mag\n", line->coeff[0], line->coeff[1]);
+  psLogMsg ("psastro", 4, "logN = %f + %f mag\n", line->coeff[0], line->coeff[1]);
   psLogMsg ("psastro", 4, "logN vs mag residuals: %f +/- %f\n", stats->sampleMedian, stats->sampleStdev);
 

trunk/psastro/src/psastroOneChip.c

@@ -1 +1 @@
 # include "psastro.h"
+
+# define REQUIRED_RECIPE_VALUE(VALUE, NAME, TYPE, MESSAGE)\
+  VALUE = psMetadataLookup##TYPE (&status, recipe, NAME); \
+  if (!status) { \
+   psError(PSASTRO_ERR_CONFIG, false, MESSAGE); \
+   return false; } 
 
 bool psastroOneChip (pmFPA *fpa, pmChip *chip, psArray *refstars, psArray *rawstars, psMetadata *recipe, psMetadata *updates) {
@@ -5 +11 @@
     bool status;
 
+    // supplied radius is in pixels
+    REQUIRED_RECIPE_VALUE (double RADIUS, "PSASTRO.MATCH.RADIUS", F32, "Failed to lookup matching radius"); 
+
+    // correct radius to FP units (physical pixel scale in microns per pixel)
+    REQUIRED_RECIPE_VALUE (double pixelScale, "PSASTRO.PIXEL.SCALE", F32, "Failed to lookup pixel scale"); 
+    RADIUS *= pixelScale;
+
+    // select the desired chip order
+    REQUIRED_RECIPE_VALUE (int order, "PSASTRO.CHIP.ORDER", S32, "failed to find single-chip fit order\n");
+
+    // allowed limits for valid solutions
+    REQUIRED_RECIPE_VALUE (float maxError, "PSASTRO.MAX.ERROR", F32, "failed to find single-chip max allowed error\n");
+
+    REQUIRED_RECIPE_VALUE (int minNstar, "PSASTRO.MIN.NSTAR", S32, "failed to find single-chip min allowed stars\n");
+
     // find initial offset / rotation
     pmAstromStats *gridStats = pmAstromGridMatch (rawstars, refstars, recipe);
     if (gridStats == NULL) {
-        psError(PSASTRO_ERR_DATA, false, "failed to find a grid match solution\n");
+        psLogMsg ("psastro", 3, "failed to find a grid match solution\n");
         return false;
     }
@@ -16 +37 @@
     pmAstromStats *stats = pmAstromGridTweak (rawstars, refstars, recipe, gridStats);
     if (stats == NULL) {
-        psError(PSASTRO_ERR_DATA, false, "failed to measure tweaked grid solution\n");
+        psLogMsg ("psastro", 3, "failed to measure tweaked grid solution\n");
         return false;
     }
@@ -25 +46 @@
     psastroUpdateChipToFPA (fpa, chip, rawstars, refstars);
 
-    // supplied radius is in pixels
-    double RADIUS = psMetadataLookupF32 (&status, recipe, "PSASTRO.MATCH.RADIUS"); 
-    if (!status) { 
-        psError(PS_ERR_IO, false, "Failed to lookup matching radius"); 
+    // use small radius to match stars
+    psArray *match = pmAstromRadiusMatchFP (rawstars, refstars, RADIUS);
+    if (match == NULL) {
+        psLogMsg ("psastro", 3, "failed to find radius-matched sources\n");
+        return false;
+    }
+
+    // create the output fit model, modify the order to correspond to the actual number of
+    // matched stars:
+    if ((match->n < 11) && (order >= 3)) order = 2;
+    if ((match->n <  7) && (order >= 2)) order = 1;
+    if ((match->n <  4) && (order >= 1)) order = 0;
+    if (order < 1) {
+        psLogMsg ("psastro", 3, "insufficient stars or invalid order: %ld stars", match->n); 
+        psFree (match);
+        psFree (stats);
+        psFree (gridStats);
         return false; 
     } 
 
-    // correct to FP units (physical pixel scale in microns per pixel)
-    double pixelScale = psMetadataLookupF32 (&status, recipe, "PSASTRO.PIXEL.SCALE");
-    if (!status) { 
-        psError(PS_ERR_IO, false, "Failed to lookup pixel scale"); 
-        return false; 
-    } 
-    RADIUS *= pixelScale;
-
-    // use small radius to match stars
-    psArray *match = pmAstromRadiusMatchFP (rawstars, refstars, RADIUS);
-    if (stats == NULL) {
-        psError(PSASTRO_ERR_UNKNOWN, false, "failed to find radius-matched sources\n");
-        return false;
-    }
-
-    // create the output fit model
-    int order = psMetadataLookupS32 (&status, recipe, "PSASTRO.CHIP.ORDER");
-    if (!status) {
-        psError(PSASTRO_ERR_UNKNOWN, false, "failed to find single-chip fit order\n");
-        return false;
-    }
+    // set masks appropriate to the Elixir DVO astrometry format
     psFree (chip->toFPA);
     chip->toFPA = psPlaneTransformAlloc (order, order);
@@ -72 +86 @@
     pmAstromFitResults *results = pmAstromMatchFit (chip->toFPA, rawstars, refstars, match, fitStats);
     if (!results) {
-        psError(PSASTRO_ERR_DATA, false, "failed to perform the matched fit\n");
+        psLogMsg ("psastro", 3, "failed to perform the matched fit\n");
         return false;
     }
     
-    // allowed limits for valid solutions
-    float maxError = psMetadataLookupF32 (&status, recipe, "PSASTRO.MAX.ERROR");
-    int minNstar = psMetadataLookupS32 (&status, recipe, "PSASTRO.MIN.NSTAR");
+    // toTPA converts from FPA units (microns) to TPA units (linear degrees)
+    float plateScale1 = hypot (fpa->toTPA->x->coeff[1][0], fpa->toTPA->x->coeff[0][1]);
+    float plateScale2 = hypot (fpa->toTPA->y->coeff[1][0], fpa->toTPA->y->coeff[0][1]);
+    float plateScale = 0.5*(plateScale1 + plateScale2)*3600.0;
 
-    // astError is the average 1D scatter in pixels ('results' are in FPA units = microns)
-    float astError = 0.5*(results->xStats->clippedStdev + results->yStats->clippedStdev) / pixelScale;
+    // pixError is the average 1D scatter in pixels ('results' are in FPA units = microns)
+    float pixError = 0.5*(results->xStats->clippedStdev + results->yStats->clippedStdev) / pixelScale;
+
+    // astError is the average 1D scatter in arcsec ('results' are in FPA units = microns)
+    float astError = 0.5*(results->xStats->clippedStdev + results->yStats->clippedStdev) * plateScale;
     int astNstar = results->yStats->clippedNvalues;
 
@@ -98 +116 @@
 
     // DVO expects NASTRO = 0 if we fail to find a solution
+    psMetadataAddF32 (updates, PS_LIST_TAIL, "PERROR",   PS_META_REPLACE, "astrometry error (pixels)", pixError);
+    psMetadataAddF32 (updates, PS_LIST_TAIL, "CERROR",   PS_META_REPLACE, "astrometry error (arcsec)", astError);
     if (validSolution) {
-        psMetadataAddS32 (updates, PS_LIST_TAIL, "NASTRO",   PS_META_REPLACE, "", astNstar);
-        psMetadataAddF32 (updates, PS_LIST_TAIL, "CPRECISE", PS_META_REPLACE, "", astError/sqrt(astNstar));
+        psMetadataAddF32 (updates, PS_LIST_TAIL, "CPRECISE", PS_META_REPLACE, "astrometry precision (arcsec)", astError/sqrt(astNstar));
+        psMetadataAddS32 (updates, PS_LIST_TAIL, "NASTRO",   PS_META_REPLACE, "number of astrometry stars", astNstar);
     } else {
-        psMetadataAddS32 (updates, PS_LIST_TAIL, "NASTRO",   PS_META_REPLACE, "", 0);
-        psMetadataAddF32 (updates, PS_LIST_TAIL, "CPRECISE", PS_META_REPLACE, "", 0.0);
+        psMetadataAddF32 (updates, PS_LIST_TAIL, "CPRECISE", PS_META_REPLACE, "astrometry precision (arcsec)", 0.0);
+        psMetadataAddS32 (updates, PS_LIST_TAIL, "NASTRO",   PS_META_REPLACE, "number of astrometry stars", 0);
     }
-    psMetadataAddF32 (updates, PS_LIST_TAIL, "CERROR",   PS_META_REPLACE, "", astError);
-    psMetadataAddF32 (updates, PS_LIST_TAIL, "EQUINOX",  PS_META_REPLACE, "", 2000.0); // XXX this is bogus: should be defined based on equinox of refstars
+    psMetadataAddF32 (updates, PS_LIST_TAIL, "EQUINOX",  PS_META_REPLACE, "equinox of ref catalog", 2000.0); // XXX this is bogus: should be defined based on equinox of refstars
 
 // write results

trunk/psastro/src/psastroParseCamera.c

@@ -22 +22 @@
     pmFPAfileDefineOutput (config, input->fpa, "PSASTRO.OUTPUT");
 
-    // Chip selection: turn on only the chips specified
+    // Chip selection: turn on only the chips specified (option is not required)
     char *chipLine = psMetadataLookupStr(NULL, config->arguments, "CHIP_SELECTIONS"); 
+    if (chipLine == NULL) psErrorClear();
+
     psArray *chips = psStringSplitArray (chipLine, ",", false);
     if (chips->n > 0) {

trunk/psastro/src/psastroRefstarSubset.c

@@ -61 +61 @@
   return true;
 }
+
+/* this test is a bit sensitive to the total number of refstars or rawstars available
+   watch out if:
+   - the fitted slopes are extremely different 
+   - the average number of stars per bin is ~1
+   
+   skip the cut in these cases?
+*/

trunk/psastro/src/psastroUtils.c

@@ -37 +37 @@
         for (int i = 0; i < fpa->chips->n; i++) {
             pmChip *chip = fpa->chips->data[i];
+            if (!chip->process || !chip->file_exists) { continue; }
             
             pixelScale1 = hypot (chip->toFPA->x->coeff[1][0], chip->toFPA->x->coeff[0][1]);
@@ -50 +51 @@
     for (int i = 0; i < fpa->chips->n; i++) {
         pmChip *chip = fpa->chips->data[i];
+        if (!chip->process || !chip->file_exists) { continue; }
 
         psPlaneTransform *toFPA = chip->toFPA;