Changeset 10855 for trunk/psastro/src/psastroOneChip.c

Timestamp:

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

Author:

eugene

Message:

extensive work to catch errors and clear leaks; implement mosaics

File:

• trunk/psastro/src/psastroOneChip.c (modified) (6 diffs)

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