Changeset 26259 for trunk/psastro/src/psastroOneChipFit.c

Timestamp:

Nov 22, 2009, 2:57:41 PM (16 years ago)

Author:

eugene

Message:

various fixes to psastro:

1) added bootstrap resampling to zero point error analysis
2) added iterative clump removal from refstars and rawstars
3) added unique reference match option
4) some improved visualizations
5) improved mosaic iterations

File:

• trunk/psastro/src/psastroOneChipFit.c (modified) (7 diffs)

trunk/psastro/src/psastroOneChipFit.c

@@ -28 +28 @@
     REQUIRED_RECIPE_VALUE (int nIter, "PSASTRO.MATCH.FIT.NITER", S32);
 
+    // for iterations >= uniqIter, require a single match per reference
+    REQUIRED_RECIPE_VALUE (int uniqIter, "PSASTRO.MATCH.UNIQ.ITER", S32);
+
     // correct radius to FP units (physical pixel scale in microns per pixel)
     REQUIRED_RECIPE_VALUE (double pixelScale, "PSASTRO.PIXEL.SCALE", F32);
@@ -61 +64 @@
             return false;
         }
+
+        if (iter >= uniqIter) {
+            psArray *unique = pmAstromRadiusMatchUniq (rawstars, refstars, match);
+            if (!unique) {
+                psLogMsg ("psastro", 3, "failed to generate a uniq set of matched sources\n");
+                return false;
+            }
+            psFree (match);
+            match = unique;
+        }
 
         // modify the order to correspond to the actual number of matched stars:
@@ -101 +114 @@
         }
 
-        // XXX allow statitic to be set by the user
+        // XXX allow statistic to be set by the user
         // fitStats = psStatsAlloc (PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV);
         fitStats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
@@ -117 +130 @@
 
         // determine fromFPA transformation and apply new transformation to raw & ref stars
-        psastroUpdateChipToFPA (fpa, chip, rawstars, refstars);
+        // this applies the transformation to all stars (including subset)
+        psastroUpdateChipToFPA (fpa, chip); // updates PSASTRO.RAWSTARS and PSASTRO.REFSTARS
 
         // 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;
-        // float astError = 0.5*(results->xStats->clippedStdev + results->yStats->clippedStdev) * plateScale;
-        float astError = 0.5*(results->xStats->robustStdev + results->yStats->robustStdev) * plateScale;
+
+        float rawXstdev = psStatsGetValue (results->xStats, psStatsStdevOption(results->xStats->options));
+        float rawYstdev = psStatsGetValue (results->yStats, psStatsStdevOption(results->yStats->options));
+
+        float astError = 0.5*(rawXstdev + rawYstdev) * plateScale;
         int astNstar = results->yStats->clippedNvalues;
         psLogMsg ("psastro", PS_LOG_INFO, "pass %d, error: %f arcsec, Nstars: %d", iter, astError, astNstar);
@@ -136 +153 @@
     float plateScale = 0.5*(fpa->toSky->Xs + fpa->toSky->Ys)*3600.0*PM_DEG_RAD;
 
+    float rawXstdev = psStatsGetValue (results->xStats, psStatsStdevOption(results->xStats->options));
+    float rawYstdev = psStatsGetValue (results->yStats, psStatsStdevOption(results->yStats->options));
+
     // 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;
-    float pixError = 0.5*(results->xStats->robustStdev + results->yStats->robustStdev) / pixelScale;
+    float pixError = 0.5*(rawXstdev + rawYstdev) / 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;
-    float astError = 0.5*(results->xStats->robustStdev + results->yStats->robustStdev) * plateScale;
+    float astError = 0.5*(rawXstdev + rawYstdev) * plateScale;
+
+    // x and y are forced to use the same subset of values:
     int astNstar = results->yStats->clippedNvalues;
 
@@ -170 +190 @@
     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
 
-    // XXX drop from here : determine fromFPA transformation and apply new transformation to raw & ref stars
-    // psastroUpdateChipToFPA (fpa, chip, rawstars, refstars);
+    // additional error measurements:
+    psMetadataAddF32 (updates, PS_LIST_TAIL, "AST_CDX",   PS_META_REPLACE, "chip astrometry X stdev (arcsec)", rawXstdev * plateScale);
+    psMetadataAddF32 (updates, PS_LIST_TAIL, "AST_CSX",   PS_META_REPLACE, "chip astrometry X systematic err (arcsec)", results->dXsys * plateScale);
+    psMetadataAddF32 (updates, PS_LIST_TAIL, "AST_CRX",   PS_META_REPLACE, "chip astrometry X 10-90 percentile (arcsec)", results->dXrange * plateScale);
+
+    psMetadataAddF32 (updates, PS_LIST_TAIL, "AST_CDY",   PS_META_REPLACE, "chip astrometry Y stdev (arcsec)", rawYstdev * plateScale);
+    psMetadataAddF32 (updates, PS_LIST_TAIL, "AST_CSY",   PS_META_REPLACE, "chip astrometry Y systematic err (arcsec)", results->dYsys * plateScale);
+    psMetadataAddF32 (updates, PS_LIST_TAIL, "AST_CRY",   PS_META_REPLACE, "chip astrometry Y 10-90 percentile (arcsec)", results->dYrange * plateScale);
 
     // XXX check if we correctly applied the new transformation:
@@ -189 +215 @@
     psFree (results);
     psFree (fitStats);
+
     return validSolution;
 }