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

Timestamp:

Dec 24, 2006, 3:56:53 PM (20 years ago)

Author:

eugene

Message:

finished up the mosaic astrometry code

File:

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

trunk/psastro/src/psastroOneChip.c

@@ -25 +25 @@
     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"); 
+        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 = pmAstromRadiusMatch (rawstars, refstars, recipe);
+    psArray *match = pmAstromRadiusMatchFP (rawstars, refstars, RADIUS);
     if (stats == NULL) {
         psError(PSASTRO_ERR_UNKNOWN, false, "failed to find radius-matched sources\n");
@@ -32 +47 @@
     }
 
-    int nX = psMetadataLookupS32 (&status, recipe, "PSASTRO.CHIP.NX");
-    if (!status) nX = 1;
-
-    int nY = psMetadataLookupS32 (&status, recipe, "PSASTRO.CHIP.NY");
-    if (!status) nY = 1;
-
-    psPlaneTransform *toFPAout = psPlaneTransformAlloc (nX, nY);
-
-    // improved fit for astrometric terms
-    if (!pmAstromMatchFit (toFPAout, rawstars, refstars, match, recipe, updates)) {
-        psError(PSASTRO_ERR_DATA, false, "failed to find a valid fitted match solution\n");
+    // 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;
     }
     psFree (chip->toFPA);
-    chip->toFPA = toFPAout;
+    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;
+            }
+        }
+    }
 
-    // write results
+    // XXX allow statitic to be set by the user
+    psStats *fitStats = psStatsAlloc (PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_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) {
+        psError(PSASTRO_ERR_DATA, false, "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");
+
+    // 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;
+    int astNstar = results->yStats->clippedNvalues;
+
+    bool validSolution = true;
+
+    // XXX should these result in errors or be handled another way?
+    psLogMsg ("psastro", PS_LOG_INFO, "astrometry solution: error: %f, Nstars: %d", astError, astNstar);
+    if (astError > maxError) {
+        psLogMsg("psastro", PS_LOG_INFO, "residual error is too large, failed to find a solution: %f > %f", astError, maxError);
+        validSolution = false;
+    }
+    if (astNstar < minNstar) {
+        psLogMsg("psastro", PS_LOG_INFO, "solution uses too few stars: %d < %d", astNstar, minNstar);
+        validSolution = false;
+    }
+
+    // DVO expects NASTRO = 0 if we fail to find a solution
+    if (validSolution) {
+        psMetadataAddS32 (updates, PS_LIST_TAIL, "NASTRO",   PS_META_REPLACE, "", astNstar);
+        psMetadataAddF32 (updates, PS_LIST_TAIL, "CPRECISE", PS_META_REPLACE, "", astError/sqrt(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, "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
+
+// write results
     psastroUpdateChipToFPA (fpa, chip, rawstars, refstars);
 
     psFree (match);
     psFree (stats);
+    psFree (results);
+    psFree (fitStats);
     psFree (gridStats);
     return true;