IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 33103


Ignore:
Timestamp:
Jan 15, 2012, 1:04:13 PM (15 years ago)
Author:
eugene
Message:

adding flat-correction features to dvo for ubercal (and other)

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

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ipp-20111122/Ohana/src/libautocode/def/FlatCorrection.d

    r33100 r33103  
    22EXTNAME FLAT_CORRECTION
    33TYPE    BINTABLE
    4 SIZE    16
     4SIZE    20
    55
     6FIELD   ID,        ID,             int
    67FIELD   x,         X,              int
    78FIELD   y,         Y,              int
    8 FIELD   offst,     OFFSET,         float
    9 FIELD   photcode,  PHOTCODE,       short
    10 FIELD   padding,   PADDING,        short
     9FIELD   offset,    OFFSET,         float

  • branches/eam_branches/ipp-20111122/Ohana/src/libdvo/src/flatcorr_io.c

    r33100 r33103  
    1212typedef struct {
    1313    int Nimage;
    14     FlatCorrectionImage *flatcorrImage;
    15 
    16     int Nflatcorr;
    17     FlatCorrection *flatcorr;
     14    FlatCorrectionImage *image;
     15    int Ncorr;
     16    FlatCorrection *corr;
     17    float ***offset; // the correction images represented as a set of arrays (same sequence as *image)
    1818} FlatCorrectionTable;
    1919
    … …  
    2525  Header theader;
    2626  FTable ftable;
    27   int i, Nflatcorr;
    28   FlatCorrection *flatcorr;
     27  int i;
     28  FlatCorrectionTable *flatcorrTable;
    2929 
    3030  f = fopen (filename, "r");
    … …  
    4747    return (NULL);
    4848  }
    49   if (!gfits_fread_ftable (f, &ftable, "FLAT_CORRECTION")) {
     49
     50  ALLOCATE (flatcorrTable, FlatCorrectionTable, 1);
     51
     52  if (!gfits_fread_ftable (f, &ftable, "FLAT_CORRECTION_IMAGE")) {
    5053    if (VERBOSE) fprintf (stderr, "can't read Flat Correction table\n");
    5154    gfits_free_header (&header);
    5255    gfits_free_matrix (&matrix);
     56    free (flatcorrTable);
    5357    fclose (f);
    5458    return (NULL);
    5559  }
    56 
    57   ALLOCATE (flatcorrTable, FlatCorrectionTable, 1);
    58 
    59   flatcorrTable.images = gfits_table_get_FlatCorrectionImage (&ftable, &flatcorrTable.Nimage, NULL);
    60   if (!skytable[0].regions) {
    61     fprintf (stderr, "ERROR: failed to read sky regions\n");
     60  flatcorrTable->image = gfits_table_get_FlatCorrectionImage (&ftable, &flatcorrTable->Nimage, NULL);
     61  if (!flatcorrTable->image) {
     62    fprintf (stderr, "ERROR: failed to read Flat Correction Images\n");
    6263    exit (2);
    6364  }
     65  gfits_free_header (&theader);
     66  // NOTE: do not free ftable or the data will be freed
    6467
    6568  if (!gfits_fread_ftable (f, &ftable, "FLAT_CORRECTION")) {
    … …  
    7073    return (NULL);
    7174  }
    72 
    73   flatcorrTable.images = gfits_table_get_FlatCorrectionImage (&ftable, &flatcorrTable.Nimage, NULL);
    74   if (!skytable[0].regions) {
    75     fprintf (stderr, "ERROR: failed to read sky regions\n");
     75  flatcorrTable->corr = gfits_table_get_FlatCorrection (&ftable, &flatcorrTable->Ncorr, NULL);
     76  if (!flatcorrTable->corr) {
     77    fprintf (stderr, "ERROR: failed to read Flat Corrections\n");
    7678    exit (2);
    7779  }
    … …  
    8082  gfits_free_matrix (&matrix);
    8183  gfits_free_header (&theader);
     84  // NOTE: do not free ftable or the data will be freed
    8285
    83   return (flatcorr);
     86  // convert the table format to map format:
     87
     88  // generate the arrays to hold the corrections
     89  ALLOCATE (flatcorrTable->offset, float **, flatcorrTable->Nimage);
     90  for (i = 0; i < flatcorrTable->Nimage; i++) {
     91      ALLOCATE (flatcorrTable->offset[i], float *, flatcorrTable->image[i]->Nx);
     92      for (j = 0; j < flatcorrTable->image[i]->Nx; j++) {
     93          ALLOCATE (flatcorrTable->offset[i][j], float, flatcorrTable->image[i]->Ny);
     94          memset (flatcorrTable->offset[i][j], 0.0, flatcorrTable->image[i]->Ny*sizeof(float));
     95      }
     96  }
     97
     98  // create the lookup table (ID -> Seq)
     99  ALLOCATE (flatcorrTable->IDtoSeq, int, MaxID);
     100  for (i = 0; i < flatcorrTable->Nimage; i++) {
     101      flatcorrTable->IDtoSeq[i] = -1;
     102  }
     103  for (i = 0; i < flatcorrTable->Nimage; i++) {
     104      ID = flatcorrTable->image[i].ID;
     105      assert (flatcorrTable->IDtoSeq[ID] == -1);
     106      flatcorrTable->IDtoSeq[ID] = i;
     107  }
     108
     109  // assign the known specific correction values
     110  for (i = 0; i < flatcorrTable->Ncorr; i++) {
     111      ID = flatcorrTable->corr[i].ID;
     112      x = flatcorrTable->corr[i].x;
     113      y = flatcorrTable->corr[i].y;
     114      seq = flatcorrTable->IDtoSeq[ID];
     115      assert (seq != -1);
     116
     117      flatcorrTable->offset[i][x][y] = flatcorrTable->corr[i].offset;     
     118  }
     119
     120  return (flatcorrTable);
    84121}
    85122
    86 int SkyTableSave (SkyTable *skytable, char *filename) {
     123int FlatCorrectionSave (FlatCorrectionTable *flatcorrTable, char *filename) {
    87124
    88125  Header header;
    89126  Matrix matrix;
    90   Header theader;
    91   FTable ftable;
     127  Header theaderImage;
     128  FTable ftableImage;
     129
     130  Header theaderCorr;
     131  FTable ftableCorr;
    92132  FILE *f;
    93133
    94134  /* make phu header (no matrix needed) */
    95   ftable.header = &theader;
    96135  gfits_init_header (&header);
    97136  header.extend = TRUE;
    … …  
    99138  gfits_create_matrix (&header, &matrix);
    100139
    101   gfits_table_set_SkyRegion (&ftable, skytable[0].regions, skytable[0].Nregions);
     140  ftableImage.header = &theaderImage;
     141  gfits_table_set_FlatCorrectionImage (&ftableImage, flatcorrTable->image, flatcorrTable->Nimage);
     142
     143  ftableCorr.header = &theaderCorr;
     144  gfits_table_set_FlatCorrectionCorr (&ftableCorr, flatcorrTable->corr, flatcorrTable->Ncorr);
    102145
    103146  f = fopen (filename, "w");
    … …  
    109152  gfits_fwrite_header  (f, &header);
    110153  gfits_fwrite_matrix  (f, &matrix);
    111   gfits_fwrite_Theader (f, &theader);
    112   gfits_fwrite_table  (f, &ftable);
     154  gfits_fwrite_Theader (f, &theaderImage);
     155  gfits_fwrite_table  (f, &ftableImage);
     156  gfits_fwrite_Theader (f, &theaderCorr);
     157  gfits_fwrite_table  (f, &ftableCorr);
    113158  fclose (f);
    114159
    … …  
    116161}
    117162
    118 SkyTable *SkyTableLoadOptimal (char *catdir, char *skyfile, char *gscfile, int readwrite, int depth, int verbose) {
     163float FlatCorrectionOffset (FlatCorrectionTable *flatcorr, int ID, int X, int Y) {
    119164
    120   char filename[256];
    121   struct stat filestat;
    122   SkyTable *sky;
     165    // validate the flat_id (not out of range?)
     166    seq = flatcorr->IDtoSeq[ID];
    123167
    124   /* first option: CATDIR/SkyTable.fits */
    125   sprintf (filename, "%s/SkyTable.fits", catdir);
    126   if (stat (filename, &filestat)) goto SKYFILE;
    127   if (!check_file_access (filename, FALSE, readwrite, verbose)) goto SKYFILE;
    128   sky = SkyTableLoad (filename, verbose);
    129   if (sky == NULL) {
    130     fprintf (stderr, "error loading sky table\n");
    131     exit (1);
    132   }
    133   return (sky);
     168    // convert X,Y to Xbin, Ybin:
     169    Xbin = MAX(MIN(X / flatcorr->image[seq].DX, flatcorr->image[seq].NX - 1), 0);
     170    Ybin = MAX(MIN(Y / flatcorr->image[seq].DY, flatcorr->image[seq].NY - 1), 0);
    134171
    135 SKYFILE:
    136   /* second option: SKYFILE */
    137   if (skyfile == NULL) goto GSCFILE;
    138   if (skyfile[0] != 0) goto GSCFILE;
    139   if (stat (skyfile, &filestat)) goto GSCFILE;
    140   if (!check_file_access (skyfile, FALSE, readwrite, verbose)) goto GSCFILE;
    141   sky = SkyTableLoad (skyfile, verbose);
    142   if (sky == NULL) {
    143     fprintf (stderr, "error loading sky table\n");
    144     return (NULL);
    145   }
    146   /* set the depths to the default depth */
    147   SkyTableSetDepth (sky, depth);
     172    // XXX warn if X,Y are out of range? not super important unless *way* out of range
    148173
    149   /* create CATDIR copy */
    150   sprintf (filename, "%s/SkyTable.fits", catdir);
    151   check_file_access (filename, FALSE, readwrite, verbose);
    152   if (!SkyTableSave (sky, filename)) return NULL;
    153 
    154   gfits_convert_SkyRegion (sky[0].regions, sizeof (SkyTable), sky[0].Nregions);
    155   return (sky);
    156 
    157 GSCFILE:
    158 
    159   /* third option: GSCFILE */
    160 
    161   if (gscfile == NULL) {
    162     fprintf (stderr, "error loading sky table from existing db (%s)\n", catdir);
    163     return (NULL);
    164   }
    165 
    166   sky = SkyTableFromGSC (gscfile, depth, verbose);
    167   if (sky == NULL) {
    168     fprintf (stderr, "error loading sky table\n");
    169     return (NULL);
    170   }
    171 
    172   /* create CATDIR copy */
    173   sprintf (filename, "%s/SkyTable.fits", catdir);
    174   check_file_access (filename, FALSE, readwrite, verbose);
    175   if (!SkyTableSave (sky, filename)) return NULL;
    176 
    177   gfits_convert_SkyRegion (sky[0].regions, sizeof (SkyRegion), sky[0].Nregions);
    178   return (sky);
     174    offset = flatcorr->offset[seq][Xbin][Ybin];
     175    return offset;
    179176}
    180177
    181 int SkyListSetFilenames (SkyList *list, char *path, char *ext) {
    182 
    183   int i;
    184   char line[256];
    185 
    186   // this generates the names, be sure to free when not needed
    187   for (i = 0; i < list[0].Nregions; i++) {
    188     sprintf (line, "%s/%s.%s", path, list[0].regions[i][0].name, ext);
    189     list[0].filename[i] = strcreate (line);
    190   }
    191 
    192   return (TRUE);
    193 }
    194 
    195 int SkyTableSetFilenames (SkyTable *sky, char *path, char *ext) {
    196 
    197   int i;
    198   char line[256];
    199 
    200   // this generates the names, be sure to free when not needed
    201   for (i = 0; i < sky[0].Nregions; i++) {
    202     sprintf (line, "%s/%s.%s", path, sky[0].regions[i].name, ext);
    203     sky[0].filename[i] = strcreate (line);
    204   }
    205 
    206   return (TRUE);
    207 }
    208 
     178// other needed operations:
     179// * assign an image to a flatcorr (based on photcode and obs_date)

  • branches/eam_branches/ipp-20111122/Ohana/src/uniphot/src/args.c

    r31160 r33103  
    8282  }
    8383
     84  UBERCAL = FALSE;
     85  if ((N = get_argument (argc, argv, "-ubercal"))) {
     86    remove_argument (N, &argc, argv);
     87    UBERCAL = TRUE;
     88  }
     89
    8490  UPDATE = FALSE;
    8591  if ((N = get_argument (argc, argv, "-update"))) {

  • branches/eam_branches/ipp-20111122/Ohana/src/uniphot/src/load_zpt_table.c

    r28717 r33103  
    4141  return zpts;
    4242}
     43
     44ZptTable *load_zpt_ubercal (char *filename, int *nzpts) {
     45
     46  int Nzpts, NZPTS;
     47  ZptTable *zpts;
     48  double zpt, mjd, zpt_err;
     49
     50  FILE *f;
     51
     52  f = fopen (filename, "r");
     53  if (!f) {
     54    fprintf (stderr, "ERROR: cannot open zpt table file %s\n", filename);
     55    exit (1);
     56  }
     57
     58  /* load in PHU segment (ignore) */
     59  ftable.header = &theader;
     60  if (!gfits_fread_header (f, &header)) {
     61    if (VERBOSE) fprintf (stderr, "can't read Flat Correction header\n");
     62    fclose (f);
     63    return (NULL);
     64  }
     65  if (!gfits_fread_matrix (f, &matrix, &header)) {
     66    if (VERBOSE) fprintf (stderr, "can't read Flat Correction matrix\n");
     67    gfits_free_header (&header);
     68    fclose (f);
     69    return (NULL);
     70  }
     71
     72  Nzpts = 0;
     73  NZPTS = 0;
     74  ALLOCATE (zpts, ZptTable, NZPTS);
     75
     76  ALLOCATE (flatcorr, FlatCorrection, 4*5*X);
     77
     78  for (i = 0; i < 5; i++) {
     79      // load data for this header
     80      if (!gfits_load_header (f, &theader)) return (FALSE);
     81
     82      // read the fits table bytes
     83      if (gfits_fread_ftable_data (f, &table, FALSE)) return (TRUE);
     84 
     85      // skip over remaining bytes in data segment
     86      fseeko (f, table[0].datasize - table[0].validsize, SEEK_CUR);
     87
     88      // load data for this header
     89      if (!gfits_load_header (f, &header)) return (FALSE);
     90
     91      // read the fits table bytes
     92      if (gfits_fread_matrix (f, &matrix, FALSE)) return (TRUE);
     93
     94      // need to create and assign to flat-field correction
     95      double *mjd = gfits_get_bintable_column_data (header, table, "mjd_obs", &type, &Nrow, &Ncol);
     96      assert (!strcmp(type, "double"));
     97
     98      double *zpts = gfits_get_bintable_column_data (header, table, "zp", &type, &Nrow, &Ncol);
     99      assert (!strcmp(type, "double"));
     100     
     101      double *zperr = gfits_get_bintable_column_data (header, table, "resid", &type, &Nrow, &Ncol);
     102      assert (!strcmp(type, "double"));
     103     
     104      NZPTS += Nrow;
     105      REALLOCATE (zpts, ZptTable, NZPTS);
     106      for (i = 0; i < Nrow; i++) {
     107          zpts[i+Nzpts].time = ohana_mjd_to_sec(mjd[i]);
     108          zpts[i+Nzpts].zpt = zp[i];
     109          zpts[i+Nzpts].zpt_err = zperr[i];
     110      }
     111      Nzpts += Nrow;
     112
     113      // the image is the flat-field correction. 
     114      // need to byteswap or whatever
     115  }
     116
     117  if (status != EOF) {
     118    fprintf (stderr, "unexpected formatting on line %d\n", Nzpts);
     119    exit (2);
     120  }
     121
     122  fprintf (stderr, "loaded %d zero points\n", Nzpts);
     123
     124  *nzpts = Nzpts;
     125  return zpts;
     126}

  • branches/eam_branches/ipp-20111122/Ohana/src/uniphot/src/setphot.c

    r27790 r33103  
    1717  if (db.dbstate == LCK_EMPTY) Shutdown ("ERROR: No images in catalog %s (1)", db.filename);
    1818
    19   zpts = load_zpt_table (argv[1], &Nzpts);
     19  if (UBERCAL) {
     20    zpts = load_zpt_ubercal (argv[1], &Nzpts);
     21  } else {
     22    zpts = load_zpt_table (argv[1], &Nzpts);
     23  }
    2024
    2125  // load images
    … …  
    3943   setphot (zpt_table)
    4044
    41  * load text table of zpts, time (photcode?)
    42  * load images
    43  * match images to zpts
    44  * set zpts
    45  * load catalogs
    46  * update detection (& averages?)
     45   setphot has 2 modes : basic & ubercal
     46
     47   * outline of basic mode:
     48
     49   ** load text table of zpts, time (photcode?)
     50   ** load images
     51   ** match images to zpts
     52   ** set zpts
     53   ** load catalogs
     54   ** update detection (& averages?)
     55
     56   * outline of ubercal mode:
     57
     58   ** load fits table of zpts & time + table of flat-field corrections
     59   ** load images
     60   ** match images to zpts
     61   ** set zpts
     62   ** load catalogs
     63   ** update detection (& averages?)
    4764
    4865 */

  • branches/eam_branches/ipp-20111122/Ohana/src/uniphot/src/update_catalog_setphot.c

    r28717 r33103  
    7171}     
    7272
    73 void update_catalog_setphot (Catalog *catalog, Image *image, off_t *index, off_t Nimage) {
     73void update_catalog_setphot (Catalog *catalog, Image *image, off_t *index, off_t Nimage, FlatCorrectionTable *flatcorr) {
    7474
    7575  off_t i, j, m, id, idx, found;
    … …  
    8787      Mcal = image[id].Mcal;
    8888      dMcal = image[id].dMcal;
     89      Mcal_offset = 0.0;
    8990
    90       catalog[0].measure[m].Mcal = Mcal;
     91      // if we know about a flat-field correction, then we need to apply the sub-chip correction
     92      flat_id = image[id].photom_map_id;
     93      if (flat_id > 0) {
     94          Mcal_offset = FlatCorrectionOffset (flatcorr, flat_id, catalog[0].measure[m].Xccd, catalog[0].measure[m].Yccd);
     95      }
     96
     97      catalog[0].measure[m].Mcal = Mcal + Mcal_offset;
    9198      catalog[0].measure[m].dMcal = dMcal;
    9299      found ++;
Note: See TracChangeset for help on using the changeset viewer.