Changeset 21208 for branches/cnb_branch_20090113/psastro/src/psastroOneChipFit.c

Timestamp:

Jan 28, 2009, 1:36:37 PM (18 years ago)

Author:

beaumont

Message:

Added QA plots to ppSub. Migrated psastro QA plots to psModules/extras.

File:

• branches/cnb_branch_20090113/psastro/src/psastroOneChipFit.c (modified) (5 diffs)

branches/cnb_branch_20090113/psastro/src/psastroOneChipFit.c

@@ -11 +11 @@
 
     // default value for match/fit : radius is in pixels
-    REQUIRED_RECIPE_VALUE (double RADIUS, "PSASTRO.MATCH.RADIUS", F32); 
+    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); 
+    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); 
+    REQUIRED_RECIPE_VALUE (double pixelScale, "PSASTRO.PIXEL.SCALE", F32);
     RADIUS *= pixelScale;
 
@@ -32 +32 @@
 
     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;
-        }
+        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");
-            return false;
-        }
+        // use small radius to match stars
+        match = pmAstromRadiusMatchFP (rawstars, refstars, radius);
+        if (match == NULL) {
+            psLogMsg ("psastro", 3, "failed to find radius-matched sources\n");
+            return false;
+        }
 
-        // modify the order to correspond to the actual number of matched stars:
-        int Ndof_min = 3;
-        int order_max = 0.5*(3 + sqrt(4*match->n - 4*Ndof_min + 1));
-        order = PS_MIN (order, order_max);
+        // modify the order to correspond to the actual number of matched stars:
+        int Ndof_min = 3;
+        int order_max = 0.5*(3 + sqrt(4*match->n - 4*Ndof_min + 1));
+        order = PS_MIN (order, order_max);
 
-        // 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 ((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);
-            return false; 
-        } 
+        if (order < 1) {
+            psLogMsg ("psastro", 3, "insufficient stars or invalid order: %ld stars", match->n);
+            psFree (match);
+            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->coeffMask[i][j] = PS_POLY_MASK_SET;
-                    chip->toFPA->y->coeffMask[i][j] = PS_POLY_MASK_SET;
-                }
-            }
-        }
+        // 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->coeffMask[i][j] = PS_POLY_MASK_SET;
+                    chip->toFPA->y->coeffMask[i][j] = PS_POLY_MASK_SET;
+                }
+            }
+        }
 
-        // 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");
+        // 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 (fitStats);
-            return false;
-        }
-    
-        // determine fromFPA transformation and apply new transformation to raw & ref stars
-        psastroUpdateChipToFPA (fpa, chip, rawstars, 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;
-        int astNstar = results->yStats->clippedNvalues;
-        psLogMsg ("psastro", PS_LOG_INFO, "pass %d, error: %f arcsec, Nstars: %d", iter, astError, astNstar);
+        // 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 (fitStats);
+            return false;
+        }
 
-        if (iter < nIter - 1) {
-            psFree (fitStats);
-            psFree (results);
-            psFree (match);
-        }
+        // determine fromFPA transformation and apply new transformation to raw & ref stars
+        psastroUpdateChipToFPA (fpa, chip, rawstars, 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;
+        int astNstar = results->yStats->clippedNvalues;
+        psLogMsg ("psastro", PS_LOG_INFO, "pass %d, error: %f arcsec, Nstars: %d", iter, astError, astNstar);
+
+        if (iter < nIter - 1) {
+            psFree (fitStats);
+            psFree (results);
+            psFree (match);
+        }
     }
 
@@ -127 +127 @@
     if (astError > maxError) {
         psLogMsg("psastro", PS_LOG_INFO, "residual error is too large, failed to find a solution: %f > %f", astError, maxError);
-        validSolution = false;
+        validSolution = false;
     }
     if (astNstar < minNstar) {
         psLogMsg("psastro", PS_LOG_INFO, "solution uses too few stars: %d < %d", astNstar, minNstar);
-        validSolution = false;
+        validSolution = false;
     }
 
@@ -138 +138 @@
     psMetadataAddF32 (updates, PS_LIST_TAIL, "CERROR",   PS_META_REPLACE, "astrometry error (arcsec)", astError);
     if (validSolution) {
-        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);
+        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 {
-        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, "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, "EQUINOX",  PS_META_REPLACE, "equinox of ref catalog", 2000.0); // XXX this is bogus: should be defined based on equinox of refstars
@@ -148 +148 @@
     // 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);
-        psastroDumpMatchedStars ("match.dat", rawstars, refstars, match);
-        psastroDumpStars (refstars, "refstars.cal.dat");
+        psastroDumpRawstars (rawstars, fpa, chip);
+        psastroDumpMatchedStars ("match.dat", rawstars, refstars, match);
+        psastroDumpStars (refstars, "refstars.cal.dat");
     }
 
-    psastroVisualPlotOneChipFit (rawstars, refstars, match, recipe);
+    pmAstromVisualPlotOneChipFit (rawstars, refstars, match, recipe);
 
     if (psTraceGetLevel("psastro.plot") > 0) {
-        psastroPlotOneChipFit (rawstars, refstars, match, recipe);
+        psastroPlotOneChipFit (rawstars, refstars, match, recipe);
     }