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

Timestamp:

Oct 3, 2007, 4:32:39 PM (19 years ago)

Author:

eugene

Message:

iterate on the match as well as the fit, added user-set radii for iterations

File:

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

trunk/psastro/src/psastroOneChip.c

@@ -1 +1 @@
 # include "psastroInternal.h"
 
-# define REQUIRED_RECIPE_VALUE(VALUE, NAME, TYPE, MESSAGE)\
+# define REQUIRED_RECIPE_VALUE(VALUE, NAME, TYPE)\
   VALUE = psMetadataLookup##TYPE (&status, recipe, NAME); \
   if (!status) { \
-   psError(PSASTRO_ERR_CONFIG, false, MESSAGE); \
-   return false; } 
+   psAbort ("Failed to find %s in recipe", NAME); }
 
 bool psastroOneChip (pmFPA *fpa, pmChip *chip, psArray *refstars, psArray *rawstars, psMetadata *recipe, psMetadata *updates) {
 
     bool status;
-    pmAstromStats *gridStats = NULL;
     pmAstromStats *stats = NULL;
 
-    // supplied radius is in pixels
-    REQUIRED_RECIPE_VALUE (double RADIUS, "PSASTRO.MATCH.RADIUS", F32, "Failed to lookup matching radius"); 
+    // default value for match/fit : radius is in pixels
+    REQUIRED_RECIPE_VALUE (double RADIUS, "PSASTRO.MATCH.RADIUS", F32); 
+
+    // run the match/fit sequence NITER times
+    REQUIRED_RECIPE_VALUE (int nIter, "PSASTRO.MATCH.FIT.NITER", S32); 
 
     // 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"); 
+    REQUIRED_RECIPE_VALUE (double pixelScale, "PSASTRO.PIXEL.SCALE", F32); 
     RADIUS *= pixelScale;
 
     // select the desired chip order
-    REQUIRED_RECIPE_VALUE (int order, "PSASTRO.CHIP.ORDER", S32, "failed to find single-chip fit order\n");
+    REQUIRED_RECIPE_VALUE (int order, "PSASTRO.CHIP.ORDER", S32);
 
     // 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");
+    REQUIRED_RECIPE_VALUE (float maxError, "PSASTRO.MAX.ERROR", F32);
+    REQUIRED_RECIPE_VALUE (int minNstar, "PSASTRO.MIN.NSTAR", S32);
 
     // do we need to get a rough initial match?
-    REQUIRED_RECIPE_VALUE (bool gridSearch, "PSASTRO.GRID.SEARCH", Bool, "failed to find chip grid-search option\n");
+    REQUIRED_RECIPE_VALUE (bool gridSearch, "PSASTRO.GRID.SEARCH", Bool);
+
+    // do we need to get a rough initial match?
+    REQUIRED_RECIPE_VALUE (int maxNstar, "PSASTRO.GRID.NSTAR.MAX", S32);
 
     if (gridSearch) {
+        // generate the bright subset of maxNstar entries (note: rawstars is sorted by S/N)
+        psArray *gridStars = psArrayAlloc (PS_MIN (maxNstar, rawstars->n));
+        for (int i = 0; (i < maxNstar) && (i < rawstars->n); i++) {
+            gridStars->data[i] = psMemIncrRefCounter (rawstars->data[i]);
+        }
+
         // find initial offset / rotation
-        gridStats = pmAstromGridMatch (rawstars, refstars, recipe);
+        pmAstromStats *gridStats = pmAstromGridMatch (gridStars, refstars, recipe);
         if (gridStats == NULL) {
             psLogMsg ("psastro", 3, "failed to find a grid match solution\n");
+            psFree (gridStars);
             return false;
         }
@@ -40 +51 @@
 
         // tweak the position by finding peak of matches stars
-        stats = pmAstromGridTweak (rawstars, refstars, recipe, gridStats);
+        stats = pmAstromGridTweak (gridStars, refstars, recipe, gridStats);
         if (stats == NULL) {
             psLogMsg ("psastro", 3, "failed to measure tweaked grid solution\n");
             psFree (gridStats);
+            psFree (gridStars);
             return false;
         }
@@ -51 +63 @@
         pmAstromGridApply (chip->toFPA, stats);
         psastroUpdateChipToFPA (fpa, chip, rawstars, refstars);
-    }
-
-    // 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");
-        psFree (stats);
         psFree (gridStats);
-        return false;
-    }
-
-    // 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; 
-    } 
-
-    // create output toFPA; set masks appropriate to the Elixir DVO astrometry format
-    psFree (chip->toFPA);
-    chip->toFPA = psPlaneTransformAlloc (order, order);
-    for (int i = 0; i <= chip->toFPA->x->nX; i++) {
-        for (int j = 0; j <= chip->toFPA->x->nY; j++) {
-            if (i + j > order) {
-                chip->toFPA->x->mask[i][j] = 1;
-                chip->toFPA->y->mask[i][j] = 1;
-            }
-        }
-    }
-
-    // XXX allow statitic to be set by the user
-    psStats *fitStats = psStatsAlloc (PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV);
-    // psStats *fitStats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
-    fitStats->clipSigma = psMetadataLookupF32 (&status, recipe, "PSASTRO.CHIP.NSIGMA");
-    fitStats->clipIter = psMetadataLookupS32 (&status, recipe, "PSASTRO.CHIP.NITER");
-
-    // improved fit for astrometric terms
-    pmAstromFitResults *results = pmAstromMatchFit (chip->toFPA, rawstars, refstars, match, fitStats);
-    if (!results) {
-        psLogMsg ("psastro", 3, "failed to perform the matched fit\n");
-        psFree (match);
-        psFree (stats);
-        psFree (fitStats);
-        psFree (gridStats);
-        return false;
-    }
+        psFree (gridStars);
+    }
+
+    psArray *match = NULL;
+    psStats *fitStats = NULL;
+    pmAstromFitResults *results = NULL;
+
+    for (int iter = 0; iter < nIter; iter++) {
+        
+        char name[128];
+
+        sprintf (name, "PSASTRO.MATCH.RADIUS.N%d", iter);
+        float radius = psMetadataLookupF32 (&status, recipe, name);
+        radius *= pixelScale;
+        if (!status || (radius == 0.0)) {
+            radius = RADIUS;
+        }
+
+
+        // use small radius to match stars
+        match = pmAstromRadiusMatchFP (rawstars, refstars, radius);
+        if (match == NULL) {
+            psLogMsg ("psastro", 3, "failed to find radius-matched sources\n");
+            psFree (stats);
+            return false;
+        }
+
+        // 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);
+            return false; 
+        } 
+
+        // create output toFPA; set masks appropriate to the Elixir DVO astrometry format
+        psFree (chip->toFPA);
+        chip->toFPA = psPlaneTransformAlloc (order, order);
+        for (int i = 0; i <= chip->toFPA->x->nX; i++) {
+            for (int j = 0; j <= chip->toFPA->x->nY; j++) {
+                if (i + j > order) {
+                    chip->toFPA->x->mask[i][j] = 1;
+                    chip->toFPA->y->mask[i][j] = 1;
+                }
+            }
+        }
+
+        // XXX allow statitic 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);
+        fitStats->clipSigma = psMetadataLookupF32 (&status, recipe, "PSASTRO.CHIP.NSIGMA");
+        fitStats->clipIter = psMetadataLookupS32 (&status, recipe, "PSASTRO.CHIP.NITER");
+
+        // improved fit for astrometric terms
+        results = pmAstromMatchFit (chip->toFPA, rawstars, refstars, match, fitStats);
+        if (!results) {
+            psLogMsg ("psastro", 3, "failed to perform the matched fit\n");
+            psFree (match);
+            psFree (stats);
+            psFree (fitStats);
+            return false;
+        }
     
+        // determine fromFPA transformation and apply new transformation to raw & ref stars
+        psastroUpdateChipToFPA (fpa, chip, rawstars, refstars);
+    
+        if (iter < nIter - 1) {
+            psFree (fitStats);
+            psFree (results);
+            psFree (match);
+        }
+
+        // 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;
+        int astNstar = results->yStats->clippedNvalues;
+        psLogMsg ("psastro", PS_LOG_INFO, "pass %d, error: %f arcsec, Nstars: %d", iter, astError, astNstar);
+    }
+
+
+
     // 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;
 
     // 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*(results->xStats->clippedStdev + results->yStats->clippedStdev) / pixelScale;
+    float pixError = 0.5*(results->xStats->robustStdev + results->yStats->robustStdev) / 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*(results->xStats->clippedStdev + results->yStats->clippedStdev) * plateScale;
+    float astError = 0.5*(results->xStats->robustStdev + results->yStats->robustStdev) * plateScale;
     int astNstar = results->yStats->clippedNvalues;
 
@@ -140 +186 @@
     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
 
-    // determine fromFPA transformation and apply new transformation to raw & ref stars
-    psastroUpdateChipToFPA (fpa, chip, rawstars, refstars);
+    // XXX drop from here : determine fromFPA transformation and apply new transformation to raw & ref stars
+    // psastroUpdateChipToFPA (fpa, chip, rawstars, refstars);
     
     // XXX check if we correctly applied the new transformation:
     if (psTraceGetLevel("psastro.dump") > 0) {
         psastroDumpRawstars (rawstars, fpa, chip);
-    }
-    if (psTraceGetLevel("psastro.dump") > 0) {
+        psastroDumpMatchedStars ("match.dat", rawstars, refstars, match);
         psastroDumpStars (refstars, "refstars.cal.dat");
     }
 
     if (psTraceGetLevel("psastro.plot") > 0) {
-        psastroPlotOneChipFit (rawstars, refstars, match, results);
+        psastroPlotOneChipFit (rawstars, refstars, match, recipe);
     }
 
@@ -159 +204 @@
     psFree (results);
     psFree (fitStats);
-    psFree (gridStats);
     return validSolution;
 }