IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 33605


Ignore:
Timestamp:
Mar 24, 2012, 5:15:54 PM (14 years ago)
Author:
eugene
Message:

adding mean ap mag, reset all secfilt.flags (no non-phot flags), skip catalog load and related for NLOOP = 0

Location:
branches/eam_branches/ipp-20111122/Ohana/src/relphot/src
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ipp-20111122/Ohana/src/relphot/src/StarOps.c

    r33564 r33605  
    167167  int i, N;
    168168  off_t j, k, m, Nmax;
    169   float Msys, Mcal, Mmos, Mgrid;
    170   double *list, *dlist;
    171   StatType stats;
     169  float Msys, Mcal, Mmos, Mgrid, Map;
     170  double *list, *dlist, *aplist, *daplist;
     171  StatType stats, apstats;
    172172  int Nsec;
    173173
    … …  
    183183  ALLOCATE (list, double, MAX (1, Nmax));
    184184  ALLOCATE (dlist, double, MAX (1, Nmax));
     185
     186  ALLOCATE (aplist, double, MAX (1, Nmax));
     187  ALLOCATE (daplist, double, MAX (1, Nmax));
    185188
    186189  for (i = 0; i < Ncatalog; i++) {
    … …  
    231234          list[N] = Msys - Mcal - Mmos - Mgrid;
    232235
     236          Map = PhotAper (&catalog[i].measure[m]);
     237          aplist[N] = Map - Mcal - Mmos - Mgrid;
     238
    233239          int myUbercalDist = getUbercalDist(m,i);
    234240          minUbercalDist = MIN(minUbercalDist, myUbercalDist);
    … …  
    286292        if (N < 1) continue;
    287293
     294        for (m = 0; m < N; m++) {
     295          daplist[m] = dlist[m];
     296        }
     297
    288298        // XXX force WT_MEAN or MEAN here?
    289299        liststats (list, dlist, N, &stats);
     300        liststats (aplist, daplist, N, &apstats);
    290301        catalog[i].found[Nsecfilt*j+Nsec] = TRUE;
    291302        switch (pass) {
    … …  
    316327        /* use sigma or error in dM for output? */
    317328        catalog[i].secfilt[Nsecfilt*j+Nsec].M  = stats.mean;
     329        catalog[i].secfilt[Nsecfilt*j+Nsec].Map = apstats.mean;
    318330        catalog[i].secfilt[Nsecfilt*j+Nsec].dM = MAX (stats.error, stats.sigma);
    319331        catalog[i].secfilt[Nsecfilt*j+Nsec].Xm = (stats.Nmeas > 1) ? 100.0*log10(stats.chisq) : NAN_S_SHORT;
    … …  
    334346  free (list);
    335347  free (dlist);
     348
     349  free (aplist);
     350  free (daplist);
    336351  return (TRUE);
    337352}

  • branches/eam_branches/ipp-20111122/Ohana/src/relphot/src/relphot.c

    r33405 r33605  
    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__); }
     3# define MARKTIME(MSG,...) {                    \
     4    float dtime;                                \
     5    gettimeofday (&stop, (void *) NULL);        \
     6    dtime = DTIME (stop, start);                \
     7    fprintf (stderr, MSG, __VA_ARGS__); }
    88
    99int main (int argc, char **argv) {
    … …  
    6666  FlatCorrectionTable *flatcorr = FlatCorrectionLoad (flatcorrfile, VERBOSE);
    6767
    68   /* load catalog data from region files (hostID is 0 since we are not a client */
    69   catalog = load_catalogs (skylist, &Ncatalog, 0, NULL);
    70   MARKTIME("-- load catalog data: %f sec\n", dtime);
     68  if (NLOOP > 0) {
     69    /* load catalog data from region files (hostID is 0 since we are not a client */
     70    catalog = load_catalogs (skylist, &Ncatalog, 0, NULL);
     71    MARKTIME("-- load catalog data: %f sec\n", dtime);
    7172 
    72   /* add in a loop over the catalogs calling dvo_catalog_chipcoords */
    73 
    74   /* match measurements with images, mosaics */
    75   initImageBins  (catalog, Ncatalog, TRUE);
    76   MARKTIME("-- make image bins: %f sec\n", dtime);
    77 
    78   initMosaicBins (catalog, Ncatalog, TRUE);
    79   initGridBins   (catalog, Ncatalog);
    80   initMrel (catalog, Ncatalog);
    81 
    82   findImages (catalog, Ncatalog, TRUE);
    83   MARKTIME("-- set up image indexes: %f sec\n", dtime);
    84 
    85   findMosaics (catalog, Ncatalog, TRUE);  /* also sets Grid values */
    86   MARKTIME("-- set up mosaic indexes: %f sec\n", dtime);
    87 
    88   SAVEPLOT = FALSE;
    89 
    90   setExclusions (catalog, Ncatalog, TRUE);
    91 
    92   global_stats (catalog, Ncatalog, flatcorr);
    93 
    94   if (PLOTSTUFF) {
    95     plot_star_coords (catalog, Ncatalog);
    96     // plot_mosaic_fields (catalog);
    97   }
    98 
    99   // if we are measuring the flat-field correction grid, we need to perform a number of iterations first:
    100   if (USE_GRID) {
     73    /* add in a loop over the catalogs calling dvo_catalog_chipcoords */
     74
     75    /* match measurements with images, mosaics */
     76    initImageBins  (catalog, Ncatalog, TRUE);
     77    MARKTIME("-- make image bins: %f sec\n", dtime);
     78
     79    initMosaicBins (catalog, Ncatalog, TRUE);
     80    initGridBins   (catalog, Ncatalog);
     81    initMrel (catalog, Ncatalog);
     82
     83    findImages (catalog, Ncatalog, TRUE);
     84    MARKTIME("-- set up image indexes: %f sec\n", dtime);
     85
     86    findMosaics (catalog, Ncatalog, TRUE);  /* also sets Grid values */
     87    MARKTIME("-- set up mosaic indexes: %f sec\n", dtime);
     88
     89    SAVEPLOT = FALSE;
     90
     91    setExclusions (catalog, Ncatalog, TRUE);
     92
     93    global_stats (catalog, Ncatalog, flatcorr);
     94
     95    if (PLOTSTUFF) {
     96      plot_star_coords (catalog, Ncatalog);
     97      // plot_mosaic_fields (catalog);
     98    }
     99
     100    // if we are measuring the flat-field correction grid, we need to perform a number of iterations first:
     101    if (USE_GRID) {
    101102      int star_toofew;
    102103
    … …  
    106107      STAR_TOOFEW = 0;
    107108      for (i = 0; i < NGRID; i++) {
    108           STAR_BAD = ID_STAR_POOR;
    109           setMrel  (catalog, Ncatalog, flatcorr);
    110           STAR_BAD  = ID_STAR_POOR | ID_STAR_FEW;
    111           setMgrid (catalog, flatcorr);
     109        STAR_BAD = ID_STAR_POOR;
     110        setMrel  (catalog, Ncatalog, flatcorr);
     111        STAR_BAD  = ID_STAR_POOR | ID_STAR_FEW;
     112        setMgrid (catalog, flatcorr);
    112113      }
    113114      STAR_BAD  = ID_STAR_POOR | ID_STAR_FEW;
    114115      STAR_TOOFEW = star_toofew;
    115   }
    116 
    117   /* determine fit values */
    118   for (i = 0; i < NLOOP; i++) {
    119     if (PLOTSTUFF) {
    120       plot_scatter (catalog, Ncatalog, flatcorr);
    121     }
    122     setMrel  (catalog, Ncatalog, flatcorr);
    123     if (PLOTSTUFF) {
    124       plot_scatter (catalog, Ncatalog, flatcorr);
    125     }
    126     setMcal  (catalog, FALSE, flatcorr);
    127     setMmos  (catalog, FALSE, flatcorr);
    128     setMgrid (catalog, flatcorr);
     116    }
     117
     118    /* determine fit values */
     119    for (i = 0; i < NLOOP; i++) {
     120      if (PLOTSTUFF) {
     121        plot_scatter (catalog, Ncatalog, flatcorr);
     122      }
     123      setMrel  (catalog, Ncatalog, flatcorr);
     124      if (PLOTSTUFF) {
     125        plot_scatter (catalog, Ncatalog, flatcorr);
     126      }
     127      setMcal  (catalog, FALSE, flatcorr);
     128      setMmos  (catalog, FALSE, flatcorr);
     129      setMgrid (catalog, flatcorr);
    129130   
     131      if (PLOTSTUFF) {
     132        plot_scatter (catalog, Ncatalog, flatcorr);
     133        plot_grid (catalog, flatcorr);
     134        plot_mosaics ();
     135        plot_images ();
     136        plot_stars (catalog, Ncatalog);
     137        plot_chisq (catalog, Ncatalog);
     138      }
     139      // if (i < NLOOP - 1) rationalize_mosaics (catalog, Ncatalog);
     140      // if (i % 6 == 1) rationalize_images ();
     141      if (i % 6 == 2) clean_measures (catalog, Ncatalog, FALSE, flatcorr);
     142      if (i % 6 == 3) clean_stars (catalog, Ncatalog);
     143      if (i % 6 == 5) clean_mosaics ();
     144      if (i % 6 == 5) clean_images ();
     145
     146      // if ((i == 1) || (i == 5) || (i ==  9) || (i == 13)) clean_measures (catalog, Ncatalog, FALSE);
     147      // if ((i == 2) || (i == 6) || (i == 10) || (i == 14)) clean_stars (catalog, Ncatalog);
     148      // if ((i == 4) || (i == 8) || (i == 12) || (i == 16)) clean_mosaics ();
     149      // if ((i == 4) || (i == 8) || (i == 12) || (i == 16)) clean_images ();
     150      global_stats (catalog, Ncatalog, flatcorr);
     151      MARKTIME("-- finished loop %d: %f sec\n", i, dtime);
     152    }
     153
    130154    if (PLOTSTUFF) {
    131155      plot_scatter (catalog, Ncatalog, flatcorr);
    … …  
    136160      plot_chisq (catalog, Ncatalog);
    137161    }
    138     // if (i < NLOOP - 1) rationalize_mosaics (catalog, Ncatalog);
    139     // if (i % 6 == 1) rationalize_images ();
    140     if (i % 6 == 2) clean_measures (catalog, Ncatalog, FALSE, flatcorr);
    141     if (i % 6 == 3) clean_stars (catalog, Ncatalog);
    142     if (i % 6 == 5) clean_mosaics ();
    143     if (i % 6 == 5) clean_images ();
    144 
    145     // if ((i == 1) || (i == 5) || (i ==  9) || (i == 13)) clean_measures (catalog, Ncatalog, FALSE);
    146     // if ((i == 2) || (i == 6) || (i == 10) || (i == 14)) clean_stars (catalog, Ncatalog);
    147     // if ((i == 4) || (i == 8) || (i == 12) || (i == 16)) clean_mosaics ();
    148     // if ((i == 4) || (i == 8) || (i == 12) || (i == 16)) clean_images ();
    149     global_stats (catalog, Ncatalog, flatcorr);
    150     MARKTIME("-- finished loop %d: %f sec\n", i, dtime);
    151   }
    152 
    153   if (PLOTSTUFF) {
    154     plot_scatter (catalog, Ncatalog, flatcorr);
    155     plot_grid (catalog, flatcorr);
    156     plot_mosaics ();
    157     plot_images ();
    158     plot_stars (catalog, Ncatalog);
    159     plot_chisq (catalog, Ncatalog);
    160   }
    161162 
    162   if (USE_GRID) dump_grid ();
    163 
    164   /* set Mcal & Mmos for bad images */
    165   setMcal  (catalog, TRUE, flatcorr);
    166   setMmos  (catalog, TRUE, flatcorr);
    167   MARKTIME("-- finalize Mcal values: %f sec\n", dtime);
    168 
    169   setMcalFinal (); // copy per-mosaic calibrations to the images
    170 
    171   if (SAVE_IMAGE_UPDATES) {
    172 
    173     FITS_DB dbX;
    174     off_t *LineNumber;
    175     Image *image, *subset;
    176     off_t Nimage, Nsubset, i, Nx;
    177     char filename[1024];
    178 
    179     gfits_db_init (&dbX);
    180     dbX.lockstate = LCK_XCLD;
    181     dbX.timeout   = 60.0;
    182     dbX.mode      = db.mode;
    183     dbX.format    = db.format;
    184 
    185     snprintf (filename, 1024, "%s.bck", ImageCat);
    186     if (!gfits_db_lock (&dbX, filename)) {
    187       fprintf (stderr, "can't lock backup image image catalog\n");
    188       return (FALSE);
    189     }
     163    if (USE_GRID) dump_grid ();
     164
     165    /* set Mcal & Mmos for bad images */
     166    setMcal  (catalog, TRUE, flatcorr);
     167    setMmos  (catalog, TRUE, flatcorr);
     168    MARKTIME("-- finalize Mcal values: %f sec\n", dtime);
     169
     170    setMcalFinal (); // copy per-mosaic calibrations to the images
     171
     172    if (SAVE_IMAGE_UPDATES) {
     173
     174      FITS_DB dbX;
     175      off_t *LineNumber;
     176      Image *image, *subset;
     177      off_t Nimage, Nsubset, i, Nx;
     178      char filename[1024];
     179
     180      gfits_db_init (&dbX);
     181      dbX.lockstate = LCK_XCLD;
     182      dbX.timeout   = 60.0;
     183      dbX.mode      = db.mode;
     184      dbX.format    = db.format;
     185
     186      snprintf (filename, 1024, "%s.bck", ImageCat);
     187      if (!gfits_db_lock (&dbX, filename)) {
     188        fprintf (stderr, "can't lock backup image image catalog\n");
     189        return (FALSE);
     190      }
    190191   
    191     // apply the changes from the image subset to the full table:
    192 
    193     // convert database table to internal structure (binary to Image)
    194     // 'image' points to the same memory as db->ftable->buffer
    195     image = gfits_table_get_Image (&db.ftable, &Nimage, &db.swapped);
    196     if (!image) {
     192      // apply the changes from the image subset to the full table:
     193
     194      // convert database table to internal structure (binary to Image)
     195      // 'image' points to the same memory as db->ftable->buffer
     196      image = gfits_table_get_Image (&db.ftable, &Nimage, &db.swapped);
     197      if (!image) {
    197198        fprintf (stderr, "ERROR: failed to read images\n");
    198199        exit (2);
    199     }
    200     gfits_scan (db.ftable.header, "NAXIS1", OFF_T_FMT, 1,  &Nx);
    201 
    202     subset = getimages (&Nsubset, &LineNumber);
    203     for (i = 0; i < Nsubset; i++) {
    204       if (LineNumber[i] == -1) continue;
    205       memcpy (&image[LineNumber[i]], &subset[i], Nx);
    206     }
    207 
    208     // copy Images.dat data (all or only the subset / vtable elements?)  I think I need to dump
    209     // the entire Image table, but with the updates in place I think this says: re-work the
    210     // ftable/vtable usage in this program to be more sensible...
    211     gfits_copy_header (&db.header,  &dbX.header);
    212     gfits_copy_matrix (&db.matrix,  &dbX.matrix);
    213     gfits_copy_header (&db.theader, &dbX.theader);
    214     gfits_copy_ftable (&db.ftable,  &dbX.ftable);
    215 
    216     dbX.ftable.header = &dbX.theader;
    217     dbX.virtual = FALSE;
    218 
    219     // save Images.dat using the copied structure
    220     if (UPDATE_CATFORMAT) {
    221       // ensure the db format is updated
    222       dbX.format = dvo_catalog_catformat (UPDATE_CATFORMAT);
    223       char photcodeFile[1024];
    224       sprintf (photcodeFile, "%s/Photcodes.dat", CATDIR);
    225       SavePhotcodesFITS (photcodeFile);
    226     }
    227     dvo_image_save (&dbX, VERBOSE);
    228     dvo_image_unlock (&dbX);
    229     MARKTIME("-- save Image.dat.bck: %f sec\n", dtime);
    230   }
     200      }
     201      gfits_scan (db.ftable.header, "NAXIS1", OFF_T_FMT, 1,  &Nx);
     202
     203      subset = getimages (&Nsubset, &LineNumber);
     204      for (i = 0; i < Nsubset; i++) {
     205        if (LineNumber[i] == -1) continue;
     206        memcpy (&image[LineNumber[i]], &subset[i], Nx);
     207      }
     208
     209      // copy Images.dat data (all or only the subset / vtable elements?)  I think I need to dump
     210      // the entire Image table, but with the updates in place I think this says: re-work the
     211      // ftable/vtable usage in this program to be more sensible...
     212      gfits_copy_header (&db.header,  &dbX.header);
     213      gfits_copy_matrix (&db.matrix,  &dbX.matrix);
     214      gfits_copy_header (&db.theader, &dbX.theader);
     215      gfits_copy_ftable (&db.ftable,  &dbX.ftable);
     216
     217      dbX.ftable.header = &dbX.theader;
     218      dbX.virtual = FALSE;
     219
     220      // save Images.dat using the copied structure
     221      if (UPDATE_CATFORMAT) {
     222        // ensure the db format is updated
     223        dbX.format = dvo_catalog_catformat (UPDATE_CATFORMAT);
     224        char photcodeFile[1024];
     225        sprintf (photcodeFile, "%s/Photcodes.dat", CATDIR);
     226        SavePhotcodesFITS (photcodeFile);
     227      }
     228      dvo_image_save (&dbX, VERBOSE);
     229      dvo_image_unlock (&dbX);
     230      MARKTIME("-- save Image.dat.bck: %f sec\n", dtime);
     231    }
     232  } // (NLOOP > 0) : this loop determines the offsets per chip
    231233
    232234  // only change the real database if -update is requested

  • branches/eam_branches/ipp-20111122/Ohana/src/relphot/src/setMrelFinal.c

    r33564 r33605  
    4848        catalog[0].secfilt[N].Ncode  = 0;
    4949        catalog[0].secfilt[N].Nused  = 0;
    50         catalog[0].secfilt[N].flags &= ~secfiltFlags;
     50        catalog[0].secfilt[N].flags  = 0; // XXX there are no *astrometry* bits in secfilt.flags
    5151        catalog[0].secfilt[N].ubercalDist = 1000;
    5252
Note: See TracChangeset for help on using the changeset viewer.