Changeset 12492

Timestamp:

Mar 18, 2007, 12:28:16 PM (19 years ago)

Author:

eugene

Message:

a variety of fixes to get psastro working in chip and mosaic mode

• various memory leaks
• when loading mosaic data, the modification of the WCS for the common focal-plane scale was getting the offset wrong (fixed in psModules)
• added some debugging dump; more plots still needed
• error handling in psastroMosaicAstrom
• loading chip and fpa dimensions from concepts
• now reading chip concepts when reading objects (in psModules)
• now supplying region and subdivisions for gradient measurements
• grid search is now optional

-

Location:

trunk/psastro/src

Files:

• psastro.c (modified) (3 diffs)
• psastro.h (modified) (1 diff)
• psastroArguments.c (modified) (1 diff)
• psastroAstromGuess.c (modified) (1 diff)
• psastroChipAstrom.c (modified) (1 diff)
• psastroCleanup.c (modified) (1 diff)
• psastroDemoDump.c (modified) (1 diff)
• psastroLoadRefstars.c (modified) (1 diff)
• psastroMosaicAstrom.c (modified) (4 diffs)
• psastroMosaicGetGrads.c (modified) (2 diffs)
• psastroMosaicGradients.c (modified) (3 diffs)
• psastroMosaicRescaleChips.c (modified) (1 diff)
• psastroOneChip.c (modified) (4 diffs)
• psastroRefstarSubset.c (modified) (1 diff)
• psastroUtils.c (modified) (2 diffs)

trunk/psastro/src/psastro.c

@@ -7 +7 @@
 
 int main (int argc, char **argv) {
+
+    pmConfig *config = NULL;
+    psArray *refs = NULL;
 
     psTimerStart ("complete");
@@ -17 +20 @@
 
     // load configuration information
-    pmConfig *config = psastroArguments (argc, argv);
+    config = psastroArguments (argc, argv);
     if (!config) usage ();
 
@@ -41 +44 @@
 
     // load the reference stars overlapping the data stars
-    psArray *refs = psastroLoadRefstars (config);
+    refs = psastroLoadRefstars (config);
     if (!refs) {
         psErrorStackPrint(stderr, "failed to load reference data\n");

trunk/psastro/src/psastro.h

@@ -78 +78 @@
 bool psastroDumpRefstars (psArray *refstars, char *filename);
 bool psastroDumpMatches (pmFPA *fpa, char *filename);
+bool psastroDumpStars (psArray *stars, char *filename);
 
 bool psastroMosaicSetAstrom_tmp (pmFPA *fpa);

trunk/psastro/src/psastroArguments.c

@@ -32 +32 @@
     if ((N = psArgumentGet (argc, argv, "-chip"))) {
         psArgumentRemove (N, &argc, argv);
-        psMetadataAddStr (config->arguments, PS_LIST_TAIL, "CHIP_SELECTIONS", PS_META_REPLACE, "", psStringCopy(argv[N]));
+        psMetadataAddStr (config->arguments, PS_LIST_TAIL, "CHIP_SELECTIONS", PS_META_REPLACE, "", argv[N]);
         psArgumentRemove (N, &argc, argv);
     }

trunk/psastro/src/psastroAstromGuess.c

@@ -40 +40 @@
     double pixelScale = psMetadataLookupF32 (&status, recipe, "PSASTRO.PIXEL.SCALE");
     if (!status) {
-        psError(PS_ERR_IO, false, "Failed to lookup pixel scale"); 
+        psError(PS_ERR_IO, true, "Failed to lookup pixel scale"); 
         return false; 
     } 

trunk/psastro/src/psastroChipAstrom.c

@@ -52 +52 @@
                 } 
 
+                char *filename = NULL;
+                char *chipname = psMetadataLookupStr (NULL, chip->concepts, "CHIP.NAME");
+
+                psStringAppend (&filename, "rawstars.ch.%s.dat", chipname);
+                psastroDumpStars (rawstars, filename);
+                psFree (filename);
+                filename = NULL;
+
+                psStringAppend (&filename, "refstars.ch.%s.dat", chipname);
+                psastroDumpStars (refstars, filename);
+                psFree (filename);
+                filename = NULL;
+
                 // save WCS and analysis metadata in update header
                 psMetadata *updates = psMetadataAlloc();

trunk/psastro/src/psastroCleanup.c

@@ -12 +12 @@
     pmConceptsDone ();
     pmConfigDone ();
-    // fprintf (stderr, "found %d leaks at %s\n", psMemCheckLeaks (0, NULL, stdout, false), "psastro");
-    fprintf (stderr, "found %d leaks at %s\n", psMemCheckLeaks (0, NULL, NULL, false), "psastro");
+    fprintf (stderr, "found %d leaks at %s\n", psMemCheckLeaks (0, NULL, stdout, false), "psastro");
+    // fprintf (stderr, "found %d leaks at %s\n", psMemCheckLeaks (0, NULL, NULL, false), "psastro");
 
     return;

trunk/psastro/src/psastroDemoDump.c

@@ -1 +1 @@
 # include "psastro.h"
+
+// this function is used for test purposes (-trace psastro.dump.psastroAstromGuess 1)
+bool psastroDumpStars (psArray *stars, char *filename) {
+
+    FILE *f = fopen (filename, "w");
+
+    for (int i = 0; i < stars->n; i++) {
+        pmAstromObj *obj = stars->data[i];
+
+        // write out the upward projections
+        fprintf (f, "%d  %f %f  %f  %f %f  %f %f  %f %f\n", i,
+                 obj->sky->r, obj->sky->d, obj->Mag, 
+                 obj->TP->x, obj->TP->y, 
+                 obj->FP->x, obj->FP->y, 
+                 obj->chip->x, obj->chip->y);
+    }
+    fclose (f);
+    return true;
+}
 
 // this function is used for test purposes (-trace psastro.dump.psastroAstromGuess 1)
 bool psastroDumpRawstars (psArray *rawstars, pmFPA *fpa, pmChip *chip) {
 
-    FILE *f1 = fopen ("rawstars.dat", "w");
-    FILE *f2 = fopen ("rawstars.down.dat", "w");
+    char *filename = NULL;
+    char *chipname = psMetadataLookupStr (NULL, chip->concepts, "CHIP.NAME");
+
+    psStringAppend (&filename, "rawstars.up.%s.dat", chipname);
+    FILE *f1 = fopen (filename, "w");
+    psFree (filename);
+    filename = NULL;
+
+    psStringAppend (&filename, "rawstars.dn.%s.dat", chipname);
+    FILE *f2 = fopen (filename, "w");
+    psFree (filename);
+    filename = NULL;
 
     for (int i = 0; i < rawstars->n; i++) {

trunk/psastro/src/psastroLoadRefstars.c

@@ -54 +54 @@
     // XXX set getstar in config?
     // psStringAppend (&getstarLine, "/home/kiawe/eugene/psconfig/dev.lin64/bin/getstar %s -D CATMODE mef -maglim %f -region %f %f %f %f -o %s", catformat, MAGmax, RAmin, DECmin, RAmax, DECmax, tempFile);
-    psStringAppend (&getstarLine, "getstar %s -D CATMODE mef -maglim %f -region %f %f %f %f -o %s", catformat, MAGmax, RAmin, DECmin, RAmax, DECmax, tempFile);
+    psStringAppend (&getstarLine, "/home/kiawe/eugene/psconfig/dev.lin64/bin/getstar %s -D CATMODE mef -maglim %f -region %f %f %f %f -o %s", catformat, MAGmax, RAmin, DECmin, RAmax, DECmax, tempFile);
     psTrace ("psastro", 3, "%s\n", getstarLine);
 

trunk/psastro/src/psastroMosaicAstrom.c

@@ -5 +5 @@
 bool psastroMosaicAstrom (pmConfig *config, psArray *refs) {
 
+    bool status;
+
     // select the current recipe
     psMetadata *recipe  = psMetadataLookupPtr (NULL, config->recipes, PSASTRO_RECIPE);
     if (!recipe) {
-        psError(PSASTRO_ERR_CONFIG, true, "Can't find PSASTRO recipe!\n");
+        psError(PSASTRO_ERR_CONFIG, false, "Can't find PSASTRO recipe!\n");
         return false;
     }
+
+    // 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; 
+    } 
 
     // select the input data sources
     pmFPAfile *input = psMetadataLookupPtr (NULL, config->files, "PSASTRO.INPUT");
     if (!input) {
-        psError(PSASTRO_ERR_CONFIG, true, "Can't find input data!\n");
+        psError(PSASTRO_ERR_CONFIG, false, "Can't find input data!\n");
         return false;
     }
@@ -27 +36 @@
     // given the existing per-chip astrometry, determine matches between raw and ref stars
     // is this needed? yes, if we didn't do SingleChip astrometry first
-    psastroMosaicSetMatch (fpa, recipe, 0);
+    if (!psastroMosaicSetMatch (fpa, recipe, 0)) {
+        psError(PSASTRO_ERR_UNKNOWN, false, "failed to match raw and ref stars for mosaic");
+        return false;
+    }
     if (psTraceGetLevel("psastro.dump") > 0) { psastroDumpMatches (fpa, "match.0.dat"); }
 
@@ -33 +45 @@
     // modify the chip->toFPA scaling to match knowledge about pixel scale,
     // then recalculate raw and ref positions
-    psastroMosaicCommonScale (fpa, recipe);
+    if (!psastroMosaicCommonScale (fpa, recipe)) {
+        psError(PSASTRO_ERR_UNKNOWN, false, "failed to set a common scale for the chips");
+        return false;
+    }
     if (psTraceGetLevel("psastro.dump") > 0) { psastroDumpMatches (fpa, "match.1.dat"); }
 
@@ -39 +54 @@
     // apply the new distortion terms up and down
     // refit the per-chip terms with linear fits only
-    psastroMosaicGradients (fpa, recipe);
+    if (!psastroMosaicGradients (fpa, recipe)) {
+        psError(PSASTRO_ERR_UNKNOWN, false, "failed to measure mosaic gradients");
+        return false;
+    }
     if (psTraceGetLevel("psastro.dump") > 0) { psastroDumpMatches (fpa, "match.2.dat"); }
 
-    psastroMosaicChipAstrom (fpa, recipe, false);
+    // measure the astrometry for the chips under the distortion term
+    if (!psastroMosaicChipAstrom (fpa, recipe, false)) {
+        psError(PSASTRO_ERR_UNKNOWN, false, "failed to measure chip astrometry in mosaic mode");
+        return false;
+    }
     if (psTraceGetLevel("psastro.dump") > 0) { psastroDumpMatches (fpa, "match.3.dat"); }
 
-    // re-perform the match with a slightly tighter circle
-    psastroMosaicSetMatch (fpa, recipe, 1);
+    // do a second pass on the distortion with improved chip positions and rotations
+    // first, re-perform the match with a slightly tighter circle
+    if (!psastroMosaicSetMatch (fpa, recipe, 1)) {
+        psError(PSASTRO_ERR_UNKNOWN, false, "failed to match raw and ref stars for mosaic (2nd pass)");
+        return false;
+    }
+    if (!psastroMosaicCommonScale (fpa, recipe)) {
+        psError(PSASTRO_ERR_UNKNOWN, false, "failed to set a common scale for the chips (2nd pass)");
+        return false;
+    }
+    if (!psastroMosaicGradients (fpa, recipe)) {
+        psError(PSASTRO_ERR_UNKNOWN, false, "failed to measure mosaic gradients (2nd pass)");
+        return false;
+    }
+    if (psTraceGetLevel("psastro.dump") > 0) { psastroDumpMatches (fpa, "match.4.dat"); }
 
-    // do a second pass on the distortion with improved chip positions and rotations
-    psastroMosaicCommonScale (fpa, recipe);
-    psastroMosaicGradients (fpa, recipe);
-    psastroMosaicChipAstrom (fpa, recipe, false);
-    psastroMosaicSetMatch (fpa, recipe, 1);
-    
+    if (!psastroMosaicChipAstrom (fpa, recipe, false)) {
+        psError(PSASTRO_ERR_UNKNOWN, false, "failed to measure chip astrometry in mosaic mode (2nd pass)");
+        return false;
+    }
+    if (psTraceGetLevel("psastro.dump") > 0) { psastroDumpMatches (fpa, "match.5.dat"); }
+
     // now fit the chips under the common distortion with higher-order terms
     // first, re-perform the match with a slightly tighter circle
-    psastroMosaicChipAstrom (fpa, recipe, true);
+    if (!psastroMosaicSetMatch (fpa, recipe, 1)) {
+        psError(PSASTRO_ERR_UNKNOWN, false, "failed to match raw and ref stars for mosaic (3rd pass)");
+        return false;
+    }
+    if (!psastroMosaicChipAstrom (fpa, recipe, true)) {
+        psError(PSASTRO_ERR_UNKNOWN, false, "failed to measure chip astrometry in mosaic mode (3rd pass)");
+        return false;
+    }
+    if (psTraceGetLevel("psastro.dump") > 0) { psastroDumpMatches (fpa, "match.6.dat"); }
 
-    // save WCS and analysis metadata in update header
+    // save WCS and analysis metadata in update header XXX this is still wrong: the pmFPAExtent
+    // function returns the FP dimensions in pixel units relative to the 0,0 corner of chip
+    // 0,0.  I need to have a function which loops over all cells, determines the pixel
+    // dimensions for each cell, converts them to chip pixels, then uses the fpa astrometry
+    // structures to get the total dimensions of the fpa in fpa units.  equivalent to
+    // pmFPAExtent
+    psRegion *region = pmFPAExtent (fpa);
+    region->x0 *= pixelScale;
+    region->x1 *= pixelScale;
+    region->y0 *= pixelScale;
+    region->y1 *= pixelScale;
+
     // XXX also need to add DATE/TIME info and NAXIS1, NAXIS2
     psMetadata *updates = psMetadataAlloc();
-    psMetadataAddS32 (updates, PS_LIST_TAIL, "NAXIS1",   PS_META_REPLACE, "fpa dimensions (um)", 30000);
-    psMetadataAddS32 (updates, PS_LIST_TAIL, "NAXIS2",   PS_META_REPLACE, "fpa dimensions (um)", 30000);
-    psMetadataAddStr (updates, PS_LIST_TAIL, "DATE-OBS", PS_META_REPLACE, "date", "2003-08-28");
-    psMetadataAddStr (updates, PS_LIST_TAIL, "UTC-OBS",  PS_META_REPLACE, "date", "7:15:46.47");
+    psMetadataAddS32 (updates, PS_LIST_TAIL, "NAXIS1",   PS_META_REPLACE, "fpa dimensions (mm)", region->x1 - region->x0);
+    psMetadataAddS32 (updates, PS_LIST_TAIL, "NAXIS2",   PS_META_REPLACE, "fpa dimensions (mm)", region->y1 - region->y0);
+    psFree (region);
     
-    pmAstromWriteBilevelMosaic (updates, fpa, NONLIN_TOL);
+    if (!pmAstromWriteBilevelMosaic (updates, fpa, NONLIN_TOL)) {
+        psAbort ("failed to save header terms");
+    }
+
     psMetadataAddMetadata (fpa->analysis, PS_LIST_TAIL, "PSASTRO.HEADER",  PS_META_REPLACE, "psastro header stats", updates);
     psFree (updates);

trunk/psastro/src/psastroMosaicGetGrads.c

@@ -15 +15 @@
         if (!chip->process || !chip->file_exists) { continue; }
         
+        psRegion *region = pmChipExtent (chip);
+
         while ((cell = pmFPAviewNextCell (view, fpa, 1)) != NULL) {
             psTrace ("psastro", 4, "Cell %d: %x %x\n", view->cell, cell->file_exists, cell->process);
@@ -36 +38 @@
 
                 // measure the local gradients for this set of stars
-                // XXX update the function prototype to accept an incoming gradient structure to which the new elements are added
-                grads = pmAstromMeasureGradients (grads, rawstars, refstars, match, recipe);
+                // the new elements are added to the incoming gradient structure 
+                grads = pmAstromMeasureGradients (grads, rawstars, refstars, match, region, 2, 2);
             }
         }
+        psFree (region);
     }
     return (grads);

trunk/psastro/src/psastroMosaicGradients.c

@@ -16 +16 @@
         if (!chip->process || !chip->file_exists) { continue; }
         
+        psRegion *region = pmChipExtent (chip);
+
         while ((cell = pmFPAviewNextCell (view, fpa, 1)) != NULL) {
             psTrace ("psastro", 4, "Cell %d: %x %x\n", view->cell, cell->file_exists, cell->process);
@@ -37 +39 @@
 
                 // measure the local gradients for this set of stars
-                gradients = pmAstromMeasureGradients (gradients, rawstars, refstars, match, recipe);
+                // XXX 2,2 are the number of test boxes on the chip.  put this in the recipe
+                gradients = pmAstromMeasureGradients (gradients, rawstars, refstars, match, region, 2, 2);
             }
         }
+        psFree (region);
     }
 
@@ -59 +63 @@
     }
 
+    // physical pixel scale in microns per pixel (FP is in physical units, chip is in pixels)
+    double pixelScale = psMetadataLookupF32 (&status, recipe, "PSASTRO.PIXEL.SCALE");
+    if (!status) {
+        psError(PS_ERR_IO, false, "Failed to lookup pixel scale"); 
+        return false; 
+    } 
+
     // fit the measured gradients with the telescope distortion model (polynomial order based on toTPA)
-    pmAstromFitDistortion (fpa, gradients, recipe);
-    psastroMosaicSetAstrom (fpa);
+    if (!pmAstromFitDistortion (fpa, gradients, pixelScale)) {
+        psError(PSASTRO_ERR_UNKNOWN, false, "failed to fit the distortion terms\n");
+        psFree (gradients);
+        psFree (view);
+        return false;
+    }
+        
+    if (!psastroMosaicSetAstrom (fpa)) {
+        psError(PSASTRO_ERR_UNKNOWN, false, "failed to apply mosaic distortion terms\n");
+        psFree (gradients);
+        psFree (view);
+        return false;
+    }
 
     psFree (gradients);

trunk/psastro/src/psastroMosaicRescaleChips.c

@@ -1 +1 @@
 # include "psastro.h"
+
+// XXX this file is no longer used
 
 // XXX what is this doing? 

trunk/psastro/src/psastroOneChip.c

@@ -10 +10 @@
 
     bool status;
+    pmAstromStats *gridStats = NULL;
+    pmAstromStats *stats = NULL;
 
     // supplied radius is in pixels
@@ -25 +27 @@
     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) {
-        psLogMsg ("psastro", 3, "failed to find a grid match solution\n");
-        return false;
+    // 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");
+
+    if (gridSearch) {
+        // find initial offset / rotation
+        gridStats = pmAstromGridMatch (rawstars, refstars, recipe);
+        if (gridStats == NULL) {
+            psLogMsg ("psastro", 3, "failed to find a grid match solution\n");
+            return false;
+        }
+        psLogMsg ("psastro", 3, "basic grid search result - offset: %f,%f pixels, rotation: %f deg\n", gridStats->offset.x, gridStats->offset.y, DEG_RAD*gridStats->angle);
+
+        // tweak the position by finding peak of matches stars
+        stats = pmAstromGridTweak (rawstars, refstars, recipe, gridStats);
+        if (stats == NULL) {
+            psLogMsg ("psastro", 3, "failed to measure tweaked grid solution\n");
+            psFree (gridStats);
+            return false;
+        }
+        psLogMsg ("psastro", 3, "tweak grid search result - offset: %f,%f pixels, rotation: %f deg\n", stats->offset.x, stats->offset.y, DEG_RAD*stats->angle);
+
+        // adjust the chip.toFPA terms only
+        pmAstromGridApply (chip->toFPA, stats);
+        psastroUpdateChipToFPA (fpa, chip, rawstars, refstars);
     }
-    psLogMsg ("psastro", 3, "basic grid search result - offset: %f,%f pixels, rotation: %f deg\n", gridStats->offset.x, gridStats->offset.y, DEG_RAD*gridStats->angle);
-
-    // tweak the position by finding peak of matches stars
-    pmAstromStats *stats = pmAstromGridTweak (rawstars, refstars, recipe, gridStats);
-    if (stats == NULL) {
-        psLogMsg ("psastro", 3, "failed to measure tweaked grid solution\n");
-        return false;
-    }
-    psLogMsg ("psastro", 3, "tweak grid search result - offset: %f,%f pixels, rotation: %f deg\n", stats->offset.x, stats->offset.y, DEG_RAD*stats->angle);
-
-    // adjust the chip.toFPA terms only
-    pmAstromGridApply (chip->toFPA, stats);
-    psastroUpdateChipToFPA (fpa, chip, rawstars, refstars);
 
     // use small radius to match stars
@@ -49 +57 @@
     if (match == NULL) {
         psLogMsg ("psastro", 3, "failed to find radius-matched sources\n");
+        psFree (stats);
+        psFree (gridStats);
         return false;
     }
@@ -85 +95 @@
     if (!results) {
         psLogMsg ("psastro", 3, "failed to perform the matched fit\n");
+        psFree (match);
+        psFree (stats);
+        psFree (fitStats);
+        psFree (gridStats);
         return false;
     }

trunk/psastro/src/psastroRefstarSubset.c

@@ -56 +56 @@
 
   if (psTraceGetLevel("psastro.dump") > 0) {
-      psastroDumpRefstars (subset, "refstars.subset.dat");
+      pmChip *chip = readout->parent->parent;
+
+      char *filename = NULL;
+      char *chipname = psMetadataLookupStr (NULL, chip->concepts, "CHIP.NAME");
+      psStringAppend (&filename, "refstars.%s.dat", chipname);
+      psastroDumpRefstars (subset, filename);
+      psFree (filename);
   }
 

trunk/psastro/src/psastroUtils.c

@@ -79 +79 @@
 bool psastroUpdateChipToFPA (pmFPA *fpa, pmChip *chip, psArray *rawstars, psArray *refstars) {
 
-    // XXX this region needs to be defined more intelligently
-    // psRegion *region = pmChipExtent (chip);
-    psRegion *region = psRegionAlloc (0, 2000, 0, 2000);
-
+    psRegion *region = pmChipExtent (chip);
+    region->x1 -= region->x0;
+    region->y1 -= region->y0;
+    region->x0 = 0;
+    region->y0 = 0;
     psFree (chip->fromFPA);
     chip->fromFPA = psPlaneTransformInvert (NULL, chip->toFPA, *region, 20);
-
-//    PS_POLY_PRINT_2D(chip->toFPA->x);
-//    PS_POLY_PRINT_2D(chip->toFPA->y);
-
-//    PS_POLY_PRINT_2D(chip->fromFPA->x);
-//    PS_POLY_PRINT_2D(chip->fromFPA->y);
+    psFree (region);
 
     for (int i = 0; i < rawstars->n; i++) {
@@ -105 +101 @@
     }
 
-    psFree (region);
     return true;
 }