Changeset 12505

Timestamp:

Mar 19, 2007, 3:58:02 PM (19 years ago)

Author:

eugene

Message:

updates to fix mosaic astrometry

Location:

trunk/pswarp/src

Files:

• Makefile.am (modified) (1 diff)
• pswarp.c (modified) (1 diff)
• pswarp.h (modified) (2 diffs)
• pswarpArguments.c (modified) (1 diff)
• pswarpCleanup.c (modified) (1 diff)
• pswarpDataLoad.c (modified) (5 diffs)
• pswarpDefine.c (modified) (2 diffs)
• pswarpMapGrid.c (modified) (2 diffs)
• pswarpParseCamera.c (modified) (3 diffs)
• pswarpTransformReadout_Opt.c (modified) (5 diffs)

trunk/pswarp/src/Makefile.am

@@ -17 +17 @@
 pswarpDataSave.c           \
 pswarpDefine.c             \
+pswarpDefineSkycell.c      \
 pswarpErrorCodes.c         \
 pswarpMapGrid.c            \

trunk/pswarp/src/pswarp.c

@@ -14 +14 @@
 
     psTimerStart ("complete");
+
+    // model inits are needed in pmSourceIO
+    // models defined in psphot/src/models are not available in psastro
+    pmModelGroupInit ();
 
     pmConfig *config = pswarpArguments (argc, argv);

trunk/pswarp/src/pswarp.h

@@ -11 +11 @@
 
 # include "pswarpErrorCodes.h"
-# define PSWARP_RECIPE "PSWARP" // Name of the recipe to use
+# define PSWARP_RECIPE  "PSWARP" // Name of the recipe to use
+# define PSASTRO_RECIPE "PSASTRO" // Name of the recipe to use
 
 // a single pswarpMap converts coordinates from one image to a second image 
@@ -47 +48 @@
 pswarpMapGrid *pswarpMapGridFromImage (pmReadout *dest, pmReadout *src, int nXpix, int nYpix);
 bool pswarpMapGridSetGrid (pswarpMapGrid *grid, int ix, int iy, int *gridX, int *gridY);
-bool pswarpMapGridNextGrid (pswarpMapGrid *grid, int gridX, int gridY, int *nextX, int *nextY);
+bool pswarpMapGridNextGrid_X (pswarpMapGrid *grid, int gridX, int *nextX);
+bool pswarpMapGridNextGrid_Y (pswarpMapGrid *grid, int gridY, int *nextY);
 double pswarpMapGridMaxError (pswarpMapGrid *grid);
 bool pswarpMapApply (double *outX, double *outY, pswarpMap *map, double inX, double inY);
 bool pswarpMapSetLocalModel (pswarpMap *map, pmReadout *dest, pmReadout *src, int ix, int iy);
+
+bool pswarpDefineSkycell (pmFPAfile **outFile, pmConfig **outConfig, pmConfig *config, const char *filename, const char *argname);

trunk/pswarp/src/pswarpArguments.c

@@ -23 +23 @@
     psMetadata *options = pmConfigRecipeOptions (config, PSWARP_RECIPE);
 
-    // astrom : used to supply astrometric parameters
-    if ((N = psArgumentGet (argc, argv, "-astrom"))) {
-        psArgumentRemove (N, &argc, argv);
-        psMetadataAddStr (options, PS_LIST_TAIL, "ASTROM_SOURCE", PS_META_REPLACE, "", argv[N]);
-        psArgumentRemove (N, &argc, argv);
-    }
+    pmConfigFileSetsMD (config->arguments, config, "ASTROM",   "-astrom", "-astromlist");
 
     // chip selection is used to limit chips to be processed
     if ((N = psArgumentGet (argc, argv, "-chip"))) {
         psArgumentRemove (N, &argc, argv);
-        psMetadataAddStr (config->arguments, PS_LIST_TAIL, "CHIP_SELECTIONS", PS_DATA_STRING, "", psStringCopy(argv[N]));
+        psMetadataAddStr (config->arguments, PS_LIST_TAIL, "CHIP_SELECTIONS", PS_DATA_STRING, "", argv[N]);
         psArgumentRemove (N, &argc, argv);
     }
+
 
     // drop the local view on the options (saved in config->arguments)

trunk/pswarp/src/pswarpCleanup.c

@@ -7 +7 @@
     psTimerStop ();
     psMemCheckCorruption(stderr, true);
-    // XXX not used in pswarp pmModelGroupCleanup ();
+    pmModelGroupCleanup ();
     psTimeFinalize ();
     pmConceptsDone ();
     pmConfigDone ();
-    // fprintf (stderr, "found %d leaks at %s\n", psMemCheckLeaks (0, NULL, stdout, false), "pswarp");
-    fprintf (stderr, "found %d leaks at %s\n", psMemCheckLeaks (0, NULL, NULL, false), "pswarp");
+    fprintf (stderr, "found %d leaks at %s\n", psMemCheckLeaks (0, NULL, stdout, false), "pswarp");
+    // fprintf (stderr, "found %d leaks at %s\n", psMemCheckLeaks (0, NULL, NULL, false), "pswarp");
 
     return;

trunk/pswarp/src/pswarpDataLoad.c

@@ -9 +9 @@
 bool pswarpDataLoad (pmConfig *config) {
 
+    bool status;
     pmChip *chip;
     pmCell *cell;
     pmReadout *readout;
     pmFPAview *view;
+
+    // select the current recipe
+    psMetadata *recipe  = psMetadataLookupPtr (NULL, config->recipes, PSASTRO_RECIPE);
+    if (!recipe) {
+        psError(PSWARP_ERR_CONFIG, false, "Can't find PSASTRO recipe needed for pixel scale!\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
@@ -21 +36 @@
     }
 
+    // use the external astrometry source if supplied
+    pmFPAfile *astrom = psMetadataLookupPtr (NULL, config->files, "PSWARP.ASTROM");
+    if (!astrom) {
+        fprintf (stderr, "*** using header astrometry\n");
+        astrom = input;
+    }
+ 
     // select the output readout
     view = pmFPAviewAlloc (0);
@@ -44 +66 @@
     // XXX need to optionally load the astrometry datafile
     bool bilevelAstrometry = false;
-    pmHDU *phu = pmFPAviewThisPHU (view, input->fpa);
+    pmHDU *phu = pmFPAviewThisPHU (view, astrom->fpa);
     if (phu) {
       char *ctype = psMetadataLookupStr (NULL, phu->header, "CTYPE1");
@@ -52 +74 @@
     }
     if (bilevelAstrometry) {
-      pmAstromReadBilevelMosaic (input->fpa, phu->header);
+      pmAstromReadBilevelMosaic (input->fpa, phu->header, pixelScale);
     } 
 
@@ -64 +86 @@
 
         // read WCS data from the corresponding header
-        pmHDU *hdu = pmFPAviewThisHDU (view, input->fpa);
+        pmHDU *hdu = pmFPAviewThisHDU (view, astrom->fpa);
         if (bilevelAstrometry) {
           pmAstromReadBilevelChip (chip, hdu->header); 
         } else {
-          // XXX get pixelScale from recipes.  does it matter?
-          float pixelScale = 13.5;
           pmAstromReadWCS (input->fpa, chip, hdu->header, pixelScale);
         }

trunk/pswarp/src/pswarpDefine.c

@@ -4 +4 @@
 bool pswarpDefine (pmConfig *config) {
 
-    // select the current recipe
-    psMetadata *recipe  = psMetadataLookupPtr (NULL, config->recipes, PSWARP_RECIPE);
+    bool status;
+    psMetadata *recipe = NULL;
+
+    // load the pixel scale from the PSASTRO recipe
+    recipe  = psMetadataLookupPtr (NULL, config->recipes, PSASTRO_RECIPE);
+    if (!recipe) {
+        psError(PSWARP_ERR_CONFIG, false, "Can't find PSASTRO recipe needed for pixel scale!\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; 
+    } 
+
+    // load the PSWARP recipe
+    recipe  = psMetadataLookupPtr (NULL, config->recipes, PSWARP_RECIPE);
     if (!recipe) {
         psError(PSWARP_ERR_CONFIG, true, "Can't find PSWARP recipe!\n");
@@ -36 +52 @@
     }
     if (bilevelAstrometry) {
-      pmAstromReadBilevelMosaic (skycell->fpa, phu->header);
+      pmAstromReadBilevelMosaic (skycell->fpa, phu->header, pixelScale);
       pmAstromReadBilevelChip (chip, hdu->header); 
     } else {
-      // XXX get pixelScale from recipes.  does it matter?
-      float pixelScale = 13.5;
       pmAstromReadWCS (skycell->fpa, chip, hdu->header, pixelScale);
     }

trunk/pswarp/src/pswarpMapGrid.c

@@ -8 +8 @@
     int j, nj;
 
+    // XXX I was trying to match the grids too closely
     // split the difference of the remainder
-    int xMin = 0.5*(src->image->numCols % nXpix);
-    int yMin = 0.5*(src->image->numRows % nYpix);
-
-    int nXpts = src->image->numCols / nXpix;
-    if (src->image->numCols % nXpix == 0) nXpts ++;
-    int nYpts = src->image->numRows / nYpix;
-    if (src->image->numRows % nYpix == 0) nYpts ++;
+    // int xMin = 0.5*(src->image->numCols % nXpix);
+    // int yMin = 0.5*(src->image->numRows % nYpix);
+    int xMin = 0;
+    int yMin = 0;
+
+    int nXpts = src->image->numCols / nXpix + 1;
+    if (src->image->numCols % nXpix) nXpts ++;
+    int nYpts = src->image->numRows / nYpix + 1;
+    if (src->image->numRows % nYpix) nYpts ++;
 
     pswarpMapGrid *grid = pswarpMapGridAlloc (nXpts, nYpts);
@@ -39 +42 @@
 }
 
-bool pswarpMapGridNextGrid (pswarpMapGrid *grid, int gridX, int gridY, int *nextX, int *nextY) {
+bool pswarpMapGridNextGrid_X (pswarpMapGrid *grid, int gridX, int *nextX) {
 
     *nextX = gridX*grid->nXpix + grid->xMin + 0.5*grid->nXpix;
+    return true;
+}
+
+bool pswarpMapGridNextGrid_Y (pswarpMapGrid *grid, int gridY, int *nextY) {
+
     *nextY = gridY*grid->nYpix + grid->yMin + 0.5*grid->nYpix;
     return true;

trunk/pswarp/src/pswarpParseCamera.c

@@ -4 +4 @@
 
     bool status;
+    pmFPAfile *skycell = NULL;
+    pmConfig *skyConfig = NULL;
 
     // the input image(s) are required arguments; they define the camera
@@ -13 +15 @@
     }
 
+    // the input image(s) are required arguments; they define the camera
+    status = false;
+    pmFPAfileDefineFromArgs (&status, config, "PSWARP.ASTROM", "ASTROM");
+    if (status) {
+        fprintf (stderr, "using supplied astrometry\n");
+    } else {
+        fprintf (stderr, "using header astrometry\n");
+    }
+
+    // the mask is not required - but must conform to input camera
+    if (!pmFPAfileBindFromArgs (NULL, input, config, "PSWARP.MASK", "MASK")) {
+        fprintf (stderr, "no mask supplied\n");
+    } 
+
+    if (!pmFPAfileBindFromArgs (NULL, input, config, "PSWARP.WEIGHT", "WEIGHT")) {
+        fprintf (stderr, "no weight supplied\n");
+    } 
+
     // the input skycell is a required argument: it defines the output image
-    status = false;
-    pmFPAfile *skycell = pmFPAfileDefineFromArgs (&status, config, "PSWARP.SKYCELL", "SKYCELL");
-    if (!skycell) {
+    // XXX we may need a different skycell structure here
+    status = pswarpDefineSkycell (&skycell, &skyConfig, config, "PSWARP.SKYCELL", "SKYCELL");
+    if (!status) {
         psError(PSWARP_ERR_CONFIG, false, "Failed to build FPA from PSWARP.SKYCELL");
         return false;
     }
 
-    // the mask is not required; the  skycell is a required argument: it defines the output image
-    status = false;
-    pmFPAfileDefineFromArgs (&status, config, "PSWARP.MASK", "MASK");
-    if (!status) {
-      // build the basic mask concept here
-    }
-
-    // the weight is not required; the  skycell is a required argument: it defines the output image
-    status = false;
-    pmFPAfileDefineFromArgs (&status, config, "PSWARP.WEIGHT", "WEIGHT");
-    if (!status) {
-      // build the basic weight concept here
-    }
-
     // these calls bind the I/O handle to the specified fpa
-    if (!pmFPAfileDefineOutput (config, skycell->fpa, "PSWARP.OUTPUT")) {
+    if (!pmFPAfileDefineOutput (skyConfig, skycell->fpa, "PSWARP.OUTPUT")) {
         psError(PSWARP_ERR_CONFIG, false, "Failed to build FPA from PSWARP.OUTPUT");
         return false;
@@ -55 +61 @@
     }
     psFree (chips);
+    psFree (skyConfig);
 
     psTrace("pswarp", 1, "Done with pswarpParseCamera...\n");

trunk/pswarp/src/pswarpTransformReadout_Opt.c

@@ -5 +5 @@
 
     int minX, minY, maxX, maxY;
-    int gridX, gridY, nextGridX, nextGridY;
+    int gridXo, gridX, gridY, nextGridX, nextGridY;
     pswarpMap *map = NULL;
 
@@ -19 +19 @@
 
     // we might want to do the rectangular regions outside of the selection independently
-    psImageInit (output->image, NAN);
+    // psImageInit (output->image, NAN);
 
     // find the output pixel range
@@ -31 +31 @@
 
     pswarpMapGridSetGrid (grid, minX, minY, &gridX, &gridY); 
-    pswarpMapGridNextGrid (grid, gridX, gridY, &nextGridX, &nextGridY);
+    pswarpMapGridNextGrid_Y (grid, gridY, &nextGridY);
     map = grid->maps[gridX][gridY];
 
@@ -38 +38 @@
         if (y >= nextGridY) {
             gridY ++;
-            pswarpMapGridNextGrid (grid, gridX, gridY, &nextGridX, &nextGridY);
+            pswarpMapGridNextGrid_Y (grid, gridY, &nextGridY);
             map = grid->maps[gridX][gridY];
         }
+
+        gridXo = gridX;
         for (int x = minX; x < maxX; x++) {
             if (x >= nextGridX) {
                 gridX ++;
-                pswarpMapGridNextGrid (grid, gridX, gridY, &nextGridX, &nextGridY);
+                pswarpMapGridNextGrid_X (grid, gridX, &nextGridX);
                 map = grid->maps[gridX][gridY];
             }
@@ -60 +62 @@
             outData[y][x] = (psF32)psImagePixelInterpolate(inImage, inPix->x, inPix->y, NULL, 1, NAN, PS_INTERPOLATE_BILINEAR);
         }
+        gridX = gridXo;
+        pswarpMapGridNextGrid_X (grid, gridX, &nextGridX);
+        map = grid->maps[gridX][gridY];
     }