IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Apr 19, 2013, 4:31:05 PM (13 years ago)
Author:
eugene
Message:

merge changes from eam_branches/ipp-20130307 : parallelize dvosecfilt, use dvo_average,measure,secfilt_init functions more broadly; move image calibration to relphot_images; add synthetic photometry referece (w-band); put relphot MARKTIME & INITTIME in macros; in dvo_clients, mextract fields which need image data use a subset image metadata table (instead of loading full Images.dat table); create FIELD and MOSAIC fields for mextract; somewhat better rules for photcode:ave & similar selections in mextract; add hpm_* concept in relastro (high-speed proper motions); fix precision errors in dvopsps; check for dvopsps exit conditions; error-bar clipping to limits of plotting window; parallelize delstar -dup-images; do NOT delete parents (because parent IDs are broken after dvomerge); minor handshake when setting up KAPA connection; avextract and related do not need to launch remote jobs on all clients if region does not include tables from all hosts; add fitplx and xsection functions to mana; enable addstar of diff cmfs (at least PS1_DV3)

Location:
trunk/Ohana
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/Ohana

  • trunk/Ohana/src/relphot/src/relphot.c

    r33963 r35416  
    11# include "relphot.h"
    22
    3 # define MARKTIME(MSG,...) {                    \
    4     float dtime;                                \
    5     gettimeofday (&stop, (void *) NULL);        \
    6     dtime = DTIME (stop, start);                \
    7     fprintf (stderr, MSG, __VA_ARGS__); }
    8 
    93int main (int argc, char **argv) {
    10 
    11   int i, status, Ncatalog;
    12   Catalog *catalog = NULL;
    13   FITS_DB db;
    14   struct timeval start, stop;
    15   SkyList *skylist = NULL;
    16 
    17   gettimeofday (&start, (void *) NULL);
    184
    195  /* get configuration info, args */
    … …  
    3622  }
    3723
    38   /* register database handle with shutdown procedure */
    39   set_db (&db);
    40   db.mode   = dvo_catalog_catmode (CATMODE);
    41   db.format = dvo_catalog_catformat (CATFORMAT);
    42 
    43   /* lock and load the image db table */
    44   status = dvo_image_lock (&db, ImageCat, 60.0, (UPDATE ? LCK_XCLD : LCK_SOFT));
    45   if (!status && UPDATE) {
    46     fprintf (stderr, "error\n");
    47     Shutdown ("ERROR: failure to lock image catalog %s", db.filename);
    48   }
    49   // if the file is missing, db.dbstate will have a value of either:
    50   // LCK_EMPTY (if UPDATE) or LCK_MISSING (if !UPDATE)
    51   if ((db.dbstate == LCK_EMPTY) || (db.dbstate == LCK_MISSING)) {
    52 
    53     // 25.0 is the nominal zero point for measurements in the database
    54     // for all measurements except PHOT_REF:
    55     // measure.M contains: -2.5*log(DN) + 2.5*log(exptime) + 25.0
    56     // for measurements using PHOT_REF:
    57     // measure.M contains: -2.5*log(DN) + 2.5*log(exptime) + True Zero Point (but exptime and ZP are known)
    58     dvo_image_create (&db, 25.0);
    59 
    60     // XXX why are we running relphot on an empty database?
    61   } else {
    62     if (!dvo_image_load (&db, VERBOSE, FALSE)) Shutdown ("can't read image catalog %s", db.filename);
    63   }
    64   MARKTIME("-- load image data: %f sec\n", dtime);
    65 
    66   /* load regions and images based on specified sky patch and/or catalog */
    67   skylist = load_images (&db, UserCatalog, &UserPatch);
    68   MARKTIME("-- load images: %f sec\n", dtime);
    69 
    70   /* unlock, if we can (else, unlocked below) */
    71   if (!UPDATE) dvo_image_unlock (&db);
    72 
    73   // load the flat correction table (if defined)
    74   char flatcorrfile[256];
    75   sprintf (flatcorrfile, "%s/flatcorr.fits", CATDIR);
    76   FlatCorrectionTable *flatcorr = FlatCorrectionLoad (flatcorrfile, VERBOSE);
    77 
    78   if (NLOOP > 0) {
    79     /* load catalog data from region files (hostID is 0 since we are not a client */
    80     catalog = load_catalogs (skylist, &Ncatalog, 0, NULL);
    81     MARKTIME("-- load catalog data: %f sec\n", dtime);
    82  
    83     /* add in a loop over the catalogs calling dvo_catalog_chipcoords */
    84 
    85     /* match measurements with images, mosaics */
    86     initImageBins  (catalog, Ncatalog, TRUE);
    87     MARKTIME("-- make image bins: %f sec\n", dtime);
    88 
    89     initMosaicBins (catalog, Ncatalog, TRUE);
    90     initGridBins   (catalog, Ncatalog);
    91     initMrel (catalog, Ncatalog);
    92 
    93     findImages (catalog, Ncatalog, TRUE);
    94     MARKTIME("-- set up image indexes: %f sec\n", dtime);
    95 
    96     findMosaics (catalog, Ncatalog, TRUE);  /* also sets Grid values */
    97     MARKTIME("-- set up mosaic indexes: %f sec\n", dtime);
    98 
    99     SAVEPLOT = FALSE;
    100 
    101     setExclusions (catalog, Ncatalog, TRUE);
    102 
    103     global_stats (catalog, Ncatalog, flatcorr);
    104 
    105     if (PLOTSTUFF) {
    106       plot_star_coords (catalog, Ncatalog);
    107       // plot_mosaic_fields (catalog);
    108     }
    109 
    110     // if we are measuring the flat-field correction grid, we need to perform a number of iterations first:
    111     if (USE_GRID) {
    112       int star_toofew;
    113 
    114       // until we finish the grid analysis, do not reject stars out-of-hand based on ID_STAR_FEW
    115       // XXX this is kind of poor: need to have a better distinctions about STAR_BAD in setMrel vs getMrel
    116       star_toofew = STAR_TOOFEW;
    117       STAR_TOOFEW = 0;
    118       for (i = 0; i < NGRID; i++) {
    119         STAR_BAD = ID_STAR_POOR;
    120         setMrel  (catalog, Ncatalog, flatcorr);
    121         STAR_BAD  = ID_STAR_POOR | ID_STAR_FEW;
    122         setMgrid (catalog, flatcorr);
    123       }
    124       STAR_BAD  = ID_STAR_POOR | ID_STAR_FEW;
    125       STAR_TOOFEW = star_toofew;
    126     }
    127 
    128     /* determine fit values */
    129     for (i = 0; i < NLOOP; i++) {
    130       if (PLOTSTUFF) {
    131         plot_scatter (catalog, Ncatalog, flatcorr);
    132       }
    133       setMrel  (catalog, Ncatalog, flatcorr); // threaded
    134       if (PLOTSTUFF) {
    135         plot_scatter (catalog, Ncatalog, flatcorr);
    136       }
    137       setMcal  (catalog, FALSE, flatcorr);
    138       setMmos  (catalog, FALSE, flatcorr);
    139       setMgrid (catalog, flatcorr);
    140       MARKTIME("-- set Mrel, Mcal, Mmos, Mgrid : %f sec\n", dtime);
    141    
    142       if (PLOTSTUFF) {
    143         plot_scatter (catalog, Ncatalog, flatcorr);
    144         plot_grid (catalog, flatcorr);
    145         plot_mosaics ();
    146         plot_images ();
    147         plot_stars (catalog, Ncatalog);
    148         plot_chisq (catalog, Ncatalog);
    149       }
    150       // if (i < NLOOP - 1) rationalize_mosaics (catalog, Ncatalog);
    151       // if (i % 6 == 1) rationalize_images ();
    152 
    153       // NOTE : in the past, I was not iterating enough before cleaning.  make sure we do
    154       // at least 8 loops first -- that should get the systematic errors down to the ~1%
    155       // level, even in cases where we have an even split between photometric data and
    156       // data with 1 mag of extinction.
    157 
    158       if ((i > 8) && (i % 8 == 2)) clean_measures (catalog, Ncatalog, FALSE, flatcorr);
    159       if ((i > 8) && (i % 8 == 3)) clean_stars (catalog, Ncatalog);
    160       if ((i > 8) && (i % 8 == 5)) clean_mosaics ();
    161       if ((i > 8) && (i % 8 == 5)) clean_images ();
    162 
    163       // if ((i == 1) || (i == 5) || (i ==  9) || (i == 13)) clean_measures (catalog, Ncatalog, FALSE);
    164       // if ((i == 2) || (i == 6) || (i == 10) || (i == 14)) clean_stars (catalog, Ncatalog);
    165       // if ((i == 4) || (i == 8) || (i == 12) || (i == 16)) clean_mosaics ();
    166       // if ((i == 4) || (i == 8) || (i == 12) || (i == 16)) clean_images ();
    167       if (i % 3 == 2) global_stats (catalog, Ncatalog, flatcorr);
    168       MARKTIME("-- finished loop %d: %f sec\n", i, dtime);
    169     }
    170 
    171     if (PLOTSTUFF) {
    172       plot_scatter (catalog, Ncatalog, flatcorr);
    173       plot_grid (catalog, flatcorr);
    174       plot_mosaics ();
    175       plot_images ();
    176       plot_stars (catalog, Ncatalog);
    177       plot_chisq (catalog, Ncatalog);
    178     }
    179  
    180     if (USE_GRID) dump_grid ();
    181 
    182     /* set Mcal & Mmos for bad images */
    183     setMcal  (catalog, TRUE, flatcorr);
    184     setMmos  (catalog, TRUE, flatcorr);
    185     MARKTIME("-- finalize Mcal values: %f sec\n", dtime);
    186 
    187     setMcalFinal (); // copy per-mosaic calibrations to the images
    188 
    189     if (SAVE_IMAGE_UPDATES) {
    190 
    191       FITS_DB dbX;
    192       off_t *LineNumber;
    193       Image *image, *subset;
    194       off_t Nimage, Nsubset, i, Nx;
    195       char filename[1024];
    196 
    197       gfits_db_init (&dbX);
    198       dbX.lockstate = LCK_XCLD;
    199       dbX.timeout   = 60.0;
    200       dbX.mode      = db.mode;
    201       dbX.format    = db.format;
    202 
    203       snprintf (filename, 1024, "%s.bck", ImageCat);
    204       if (!gfits_db_lock (&dbX, filename)) {
    205         fprintf (stderr, "can't lock backup image image catalog\n");
    206         return (FALSE);
    207       }
    208    
    209       // apply the changes from the image subset to the full table:
    210 
    211       // convert database table to internal structure (binary to Image)
    212       // 'image' points to the same memory as db->ftable->buffer
    213       image = gfits_table_get_Image (&db.ftable, &Nimage, &db.swapped);
    214       if (!image) {
    215         fprintf (stderr, "ERROR: failed to read images\n");
    216         exit (2);
    217       }
    218       gfits_scan (db.ftable.header, "NAXIS1", OFF_T_FMT, 1,  &Nx);
    219 
    220       subset = getimages (&Nsubset, &LineNumber);
    221       for (i = 0; i < Nsubset; i++) {
    222         if (LineNumber[i] == -1) continue;
    223         memcpy (&image[LineNumber[i]], &subset[i], Nx);
    224       }
    225 
    226       // copy Images.dat data (all or only the subset / vtable elements?)  I think I need to dump
    227       // the entire Image table, but with the updates in place I think this says: re-work the
    228       // ftable/vtable usage in this program to be more sensible...
    229       gfits_copy_header (&db.header,  &dbX.header);
    230       gfits_copy_matrix (&db.matrix,  &dbX.matrix);
    231       gfits_copy_header (&db.theader, &dbX.theader);
    232       gfits_copy_ftable (&db.ftable,  &dbX.ftable);
    233 
    234       dbX.ftable.header = &dbX.theader;
    235       dbX.virtual = FALSE;
    236 
    237       // save Images.dat using the copied structure
    238       if (UPDATE_CATFORMAT) {
    239         // ensure the db format is updated
    240         dbX.format = dvo_catalog_catformat (UPDATE_CATFORMAT);
    241         char photcodeFile[1024];
    242         sprintf (photcodeFile, "%s/Photcodes.dat", CATDIR);
    243         SavePhotcodesFITS (photcodeFile);
    244       }
    245       dvo_image_save (&dbX, VERBOSE);
    246       dvo_image_unlock (&dbX);
    247       MARKTIME("-- save Image.dat.bck: %f sec\n", dtime);
    248     }
    249 
    250     /* at this point, we have correct cal coeffs in the image/mosaic structures */
    251     for (i = 0; i < Ncatalog; i++) {
    252       free_tiny_values (&catalog[i]);
    253       dvo_catalog_free (&catalog[i]);
    254     }
    255     freeImageBins (Ncatalog, TRUE);
    256     freeMosaicBins (Ncatalog, TRUE);
    257     freeGridBins (Ncatalog);
    258 
    259   } // (NLOOP > 0) : this loop determines the offsets per chip
    260 
    261   reload_images (&db);
    262 
    263   // only change the real database files if -update is requested
    264   if (!UPDATE) exit (0);
    265  
    266   /* load catalog data from region files, update Mrel include all data */
    267   reload_catalogs (skylist, flatcorr, 0, NULL);
    268   MARKTIME("-- updated all catalogs: %f sec\n", dtime);
    269 
    270   if (UPDATE_CATFORMAT) {
    271     // ensure the db format is updated
    272     db.format = dvo_catalog_catformat (UPDATE_CATFORMAT);
    273   }
    274   dvo_image_update (&db, VERBOSE);
    275   dvo_image_unlock (&db);
     24  relphot_images ();
    27625
    27726  exit (0);
Note: See TracChangeset for help on using the changeset viewer.