IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 33104


Ignore:
Timestamp:
Jan 15, 2012, 4:28:36 PM (15 years ago)
Author:
eugene
Message:

finish adding flat-field correction to libdvo & uniphot (setphot) -- it builds, but is not tested yet

Location:
branches/eam_branches/ipp-20111122/Ohana/src
Files:
1 added
9 edited

Legend:

Unmodified
Added
Removed

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

    r33103 r33104  
    22EXTNAME FLAT_CORRECTION
    33TYPE    BINTABLE
    4 SIZE    20
     4SIZE    16
    55
    66FIELD   ID,        ID,             int

  • branches/eam_branches/ipp-20111122/Ohana/src/libdvo/Makefile

    r31636 r33104  
    7070$(SRC)/dvo_convert_PS1_V2.$(ARCH).o \
    7171$(SRC)/dvo_convert_PS1_REF.$(ARCH).o \
     72$(SRC)/flatcorr_io.$(ARCH).o    \
    7273$(SRC)/skyregion_io.$(ARCH).o    \
    7374$(SRC)/skyregion_gsc.$(ARCH).o    \

  • branches/eam_branches/ipp-20111122/Ohana/src/libdvo/include/dvo.h

    r31663 r33104  
    205205} SkyList;
    206206
     207typedef struct {
     208    off_t Nimage;
     209    FlatCorrectionImage *image;
     210    off_t Ncorr;
     211    int *IDtoSeq;
     212    FlatCorrection *corr;
     213    float ***offset; // the correction images represented as a set of arrays (same sequence as *image)
     214} FlatCorrectionTable;
     215
    207216// special-case function:
    208217CMF_PS1_V2 *gfits_table_get_CMF_PS1_V1_Alt (FTable *ftable, off_t *Ndata, char *swapped);
    … …  
    484493void dvo_image_create (FITS_DB *db, double ZeroPoint);
    485494
     495/* flatcorr APIs */
     496FlatCorrectionTable *FlatCorrectionLoad (char *filename, int VERBOSE);
     497int FlatCorrectionInternal(FlatCorrectionTable *flatcorrTable);
     498int FlatCorrectionSave (FlatCorrectionTable *flatcorrTable, char *filename);
     499float FlatCorrectionOffset (FlatCorrectionTable *flatcorr, int ID, int X, int Y);
     500
    486501/* skyregion APIs */
    487502int        SkyTableSave            PROTO((SkyTable *table, char *filename));

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

    r33103 r33104  
    1010 */
    1111
    12 typedef struct {
    13     int Nimage;
    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)
    18 } FlatCorrectionTable;
    19 
    2012FlatCorrectionTable *FlatCorrectionLoad (char *filename, int VERBOSE) {
    2113
    … …  
    2517  Header theader;
    2618  FTable ftable;
    27   int i;
    2819  FlatCorrectionTable *flatcorrTable;
    2920 
    … …  
    8576
    8677  // convert the table format to map format:
     78  FlatCorrectionInternal (flatcorrTable);
    8779
    88   // generate the arrays to hold the corrections
     80  return (flatcorrTable);
     81}
     82
     83int FlatCorrectionInternal(FlatCorrectionTable *flatcorrTable) {
     84
     85  int i, j;
     86
     87  // assert that the internal arrays are not yet allocated?
     88
     89  // generate the arrays to hold the corrections & initialize
    8990  ALLOCATE (flatcorrTable->offset, float **, flatcorrTable->Nimage);
    9091  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       }
     92    ALLOCATE (flatcorrTable->offset[i], float *, flatcorrTable->image[i].Nx);
     93    for (j = 0; j < flatcorrTable->image[i].Nx; j++) {
     94      ALLOCATE (flatcorrTable->offset[i][j], float, flatcorrTable->image[i].Ny);
     95      memset (flatcorrTable->offset[i][j], 0.0, flatcorrTable->image[i].Ny*sizeof(float));
     96    }
    9697  }
    9798
    9899  // create the lookup table (ID -> Seq)
     100
     101  // find the max value of ID
     102  int MaxID = 0;
     103  for (i = 0; i < flatcorrTable->Nimage; i++) {
     104    MaxID = MAX(flatcorrTable->image[i].ID, MaxID);
     105  }
     106
     107  // generate the index and init values to -1
    99108  ALLOCATE (flatcorrTable->IDtoSeq, int, MaxID);
     109  for (i = 0; i < MaxID; i++) {
     110    flatcorrTable->IDtoSeq[i] = -1;
     111  }
     112
     113  // assign the ID values for each image
    100114  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;
     115    int ID = flatcorrTable->image[i].ID;
     116    assert (flatcorrTable->IDtoSeq[ID] == -1);
     117    flatcorrTable->IDtoSeq[ID] = i;
    107118  }
    108119
    109120  // assign the known specific correction values
    110121  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);
     122    int ID = flatcorrTable->corr[i].ID;
     123    int x = flatcorrTable->corr[i].x;
     124    int y = flatcorrTable->corr[i].y;
     125    int seq = flatcorrTable->IDtoSeq[ID];
     126    assert (seq != -1);
    116127
    117       flatcorrTable->offset[i][x][y] = flatcorrTable->corr[i].offset;     
     128    flatcorrTable->offset[i][x][y] = flatcorrTable->corr[i].offset;     
    118129  }
    119130
    120   return (flatcorrTable);
     131  return TRUE;
    121132}
    122133
    … …  
    142153
    143154  ftableCorr.header = &theaderCorr;
    144   gfits_table_set_FlatCorrectionCorr (&ftableCorr, flatcorrTable->corr, flatcorrTable->Ncorr);
     155  gfits_table_set_FlatCorrection (&ftableCorr, flatcorrTable->corr, flatcorrTable->Ncorr);
    145156
    146157  f = fopen (filename, "w");
    … …  
    163174float FlatCorrectionOffset (FlatCorrectionTable *flatcorr, int ID, int X, int Y) {
    164175
    165     // validate the flat_id (not out of range?)
    166     seq = flatcorr->IDtoSeq[ID];
     176  // validate the flat_id (not out of range?)
     177  int seq = flatcorr->IDtoSeq[ID];
    167178
    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);
     179  // convert X,Y to Xbin, Ybin:
     180  int Xbin = MAX(MIN(X * flatcorr->image[seq].Nx / flatcorr->image[seq].DX, flatcorr->image[seq].Nx - 1), 0);
     181  int Ybin = MAX(MIN(Y * flatcorr->image[seq].Ny / flatcorr->image[seq].DY, flatcorr->image[seq].Ny - 1), 0);
    171182
    172     // XXX warn if X,Y are out of range? not super important unless *way* out of range
     183  // XXX warn if X,Y are out of range? not super important unless *way* out of range
    173184
    174     offset = flatcorr->offset[seq][Xbin][Ybin];
    175     return offset;
     185  float offset = flatcorr->offset[seq][Xbin][Ybin];
     186  return offset;
    176187}
    177188

  • branches/eam_branches/ipp-20111122/Ohana/src/uniphot/Makefile

    r31160 r33104  
    5252$(SRC)/load_images.$(ARCH).o        \
    5353$(SRC)/match_zpts_to_images.$(ARCH).o       \
     54$(SRC)/match_flatcorr_to_images.$(ARCH).o           \
    5455$(SRC)/update_catalog_setphot.$(ARCH).o \
    5556$(SRC)/SetSignals.$(ARCH).o         \

  • branches/eam_branches/ipp-20111122/Ohana/src/uniphot/include/uniphot.h

    r31635 r33104  
    6666char         STATMODE[64];
    6767int          VERBOSE;
     68int          UBERCAL; // load the supplied ubercal zero point fits table (with flat-field corrections)
    6869int          NLOOP;
    6970int          TimeSelect;
    … …  
    132133Image        *load_images_setphot    PROTO((FITS_DB *db, off_t *Nimage));
    133134int           match_zpts_to_images   PROTO((Image *image, off_t Nimage, ZptTable *zpts, int Nzpts));
    134 int           update_setphot         PROTO((Image *image, off_t Nimage));
    135 void          update_catalog_setphot PROTO((Catalog *catalog, Image *image, off_t *index, off_t Nimage));
     135int           update_setphot         PROTO((Image *image, off_t Nimage, FlatCorrectionTable *flatcorr));
     136void          update_catalog_setphot PROTO((Catalog *catalog, Image *image, off_t *index, off_t Nimage, FlatCorrectionTable *flatcorr));
     137
     138ZptTable *load_zpt_ubercal (char *filename, int *nzpts, FlatCorrectionTable *flatcorrTable);
     139int match_flatcorr_to_images (Image *image, off_t Nimage, FlatCorrectionTable *flatcorrTable);
    136140
    137141/*** time/coord conversion functions not supplied by libohana ***/

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

    r33103 r33104  
    4242}
    4343
    44 ZptTable *load_zpt_ubercal (char *filename, int *nzpts) {
     44# define NFILTER 5
     45# define NSEASON 4
    4546
     47char   filter[NFILTER] = {'g', 'r', 'i', 'z', 'y'};
     48double tstart[NSEASON] = {55000.0, 55296.0, 55327.0, 55662.0};
     49double tstop[NSEASON]  = {55296.0, 55327.0, 55662.0, 60000.0};
     50
     51# define NCHIP_X 8
     52# define NCHIP_Y 8
     53# define NCELL_X 2
     54# define NCELL_Y 2
     55
     56ZptTable *load_zpt_ubercal (char *filename, int *nzpts, FlatCorrectionTable *flatcorrTable) {
     57
     58  int i, nfilter, nseason, ix, iy, ixc, iyc, Ncol, seq;
     59  off_t Nrow;
     60  char type[16];
    4661  int Nzpts, NZPTS;
    4762  ZptTable *zpts;
    48   double zpt, mjd, zpt_err;
     63  Header header;
     64  Header theader;
     65  Matrix matrix;
     66  FTable ftable;
    4967
    5068  FILE *f;
    … …  
    5775
    5876  /* load in PHU segment (ignore) */
    59   ftable.header = &theader;
    6077  if (!gfits_fread_header (f, &header)) {
    6178    if (VERBOSE) fprintf (stderr, "can't read Flat Correction header\n");
    … …  
    7491  ALLOCATE (zpts, ZptTable, NZPTS);
    7592
    76   ALLOCATE (flatcorr, FlatCorrection, 4*5*X);
     93  // we have 5 filters, and 4 flat-field correction sets for each
     94  flatcorrTable->Ncorr = NFILTER*NSEASON*NCHIP_X*NCHIP_Y*NCELL_X*NCELL_Y;
     95  flatcorrTable->Nimage = NFILTER*NSEASON*NCHIP_X*NCHIP_Y;
    7796
    78   for (i = 0; i < 5; i++) {
    79       // load data for this header
    80       if (!gfits_load_header (f, &theader)) return (FALSE);
     97  ALLOCATE (flatcorrTable->corr, FlatCorrection, flatcorrTable->Ncorr);
     98  ALLOCATE (flatcorrTable->image, FlatCorrectionImage, flatcorrTable->Nimage);
    8199
    82       // read the fits table bytes
    83       if (gfits_fread_ftable_data (f, &table, FALSE)) return (TRUE);
     100  int corrID = 1;
     101  int Nimage = 0;
     102  ftable.header = &theader;
     103  for (nfilter = 0; nfilter < NFILTER; nfilter++) {
     104    // load data for this header
     105    if (!gfits_load_header (f, &theader)) return (NULL);
     106
     107    // read the fits table bytes
     108    if (gfits_fread_ftable_data (f, &ftable, FALSE)) return (NULL);
    84109 
    85       // skip over remaining bytes in data segment
    86       fseeko (f, table[0].datasize - table[0].validsize, SEEK_CUR);
     110    // skip over remaining bytes in data segment
     111    fseeko (f, ftable.datasize - ftable.validsize, SEEK_CUR);
    87112
    88       // load data for this header
    89       if (!gfits_load_header (f, &header)) return (FALSE);
     113    // need to create and assign to flat-field correction
     114    double *mjd = gfits_get_bintable_column_data (&theader, &ftable, "mjd_obs", type, &Nrow, &Ncol);
     115    assert (!strcmp(type, "double"));
    90116
    91       // read the fits table bytes
    92       if (gfits_fread_matrix (f, &matrix, FALSE)) return (TRUE);
     117    double *zp = gfits_get_bintable_column_data (&theader, &ftable, "zp", type, &Nrow, &Ncol);
     118    assert (!strcmp(type, "double"));
     119     
     120    double *zperr = gfits_get_bintable_column_data (&theader, &ftable, "resid", type, &Nrow, &Ncol);
     121    assert (!strcmp(type, "double"));
     122     
     123    NZPTS += Nrow;
     124    REALLOCATE (zpts, ZptTable, NZPTS);
     125    for (i = 0; i < Nrow; i++) {
     126      zpts[i+Nzpts].time = ohana_mjd_to_sec(mjd[i]);
     127      zpts[i+Nzpts].zpt = zp[i];
     128      zpts[i+Nzpts].zpt_err = zperr[i];
     129    }
     130    Nzpts += Nrow;
    93131
    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"));
     132    // the image contains the flat-field corrections
    97133
    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];
     134    // load data for this header
     135    if (!gfits_load_header (f, &header)) return (NULL);
     136
     137    // read the fits table bytes
     138    if (gfits_fread_matrix (f, &matrix, FALSE)) return (NULL);
     139
     140    for (nseason = 0; nseason < NSEASON; nseason++) { // seasons
     141      for (iy = 0; iy < NCHIP_Y; iy++) { // y-chip
     142        for (ix = 0; ix < NCHIP_X; ix++) { // x-chip
     143          flatcorrTable->image[Nimage].Nx = NCELL_X;
     144          flatcorrTable->image[Nimage].Ny = NCELL_Y;
     145          flatcorrTable->image[Nimage].ID = corrID;
     146          flatcorrTable->image[Nimage].DX = 4880;
     147          flatcorrTable->image[Nimage].DY = 4864;
     148          flatcorrTable->image[Nimage].tstart = ohana_mjd_to_sec(tstart[nseason]);
     149          flatcorrTable->image[Nimage].tstop  = ohana_mjd_to_sec(tstop[nseason]);
     150         
     151          // photcode name
     152          char photname[64];
     153          snprintf (photname, 64, "GPC1.%c.XY%d%d", filter[nfilter], ix, iy);
     154          // note that the XY00, XY07, etc, chips will have photcode values of 0
     155          flatcorrTable->image[Nimage].photcode = GetPhotcodeCodebyName(photname);
     156          flatcorrTable->image[Nimage].tstart = tstart[nseason];
     157          flatcorrTable->image[Nimage].tstop  = tstop[nseason];
     158         
     159          for (iyc = 0; iyc < NCELL_Y; iyc++) {
     160            for (ixc = 0; ixc < NCELL_X; ixc++) {
     161              if (ix > 3) {
     162                seq = nfilter*NSEASON*NCHIP_X*NCHIP_Y*NCELL_X*NCELL_Y + nseason*NCHIP_X*NCHIP_Y*NCELL_X*NCELL_Y + (iy*NCELL_Y + NCELL_Y - 1 - iyc)*NCHIP_X + (ix*NCELL_X + NCELL_X - 1 - ixc);
     163              } else {
     164                seq = nfilter*NSEASON*NCHIP_X*NCHIP_Y*NCELL_X*NCELL_Y + nseason*NCHIP_X*NCHIP_Y*NCELL_X*NCELL_Y + (iy*NCELL_Y + iyc)*NCHIP_X + (ix*NCELL_X + ixc);
     165              }
     166              flatcorrTable->corr[seq].x = ixc;
     167              flatcorrTable->corr[seq].y = iyc;
     168              flatcorrTable->corr[seq].offset = matrix.buffer[seq];
     169              flatcorrTable->corr[seq].ID = corrID;
     170            }
     171          }
     172          corrID ++;
     173          Nimage ++;
     174        }
    110175      }
    111       Nzpts += Nrow;
    112 
    113       // the image is the flat-field correction. 
    114       // need to byteswap or whatever
     176    }         
    115177  }
    116178
    117   if (status != EOF) {
    118     fprintf (stderr, "unexpected formatting on line %d\n", Nzpts);
    119     exit (2);
    120   }
     179  /*** convert from corr,image format to offsets ***/
     180  FlatCorrectionInternal (flatcorrTable);
    121181
    122182  fprintf (stderr, "loaded %d zero points\n", Nzpts);

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

    r33103 r33104  
    88  ZptTable *zpts;
    99  Image *image;
     10  FlatCorrectionTable flatcorrTable;
    1011
    1112  /* get configuration info, args, lockfile */
    … …  
    1819
    1920  if (UBERCAL) {
    20     zpts = load_zpt_ubercal (argv[1], &Nzpts);
     21    zpts = load_zpt_ubercal (argv[1], &Nzpts, &flatcorrTable);
     22    char flatcorrfile[64];
     23    snprintf (flatcorrfile, 64, "%s/flatcorr.fits", CATDIR);
     24    FlatCorrectionSave(&flatcorrTable, flatcorrfile);
    2125  } else {
    2226    zpts = load_zpt_table (argv[1], &Nzpts);
    … …  
    2933  match_zpts_to_images (image, Nimage, zpts, Nzpts);
    3034
    31   update_setphot(image, Nimage);
     35  if (UBERCAL) {
     36    match_flatcorr_to_images (image, Nimage, &flatcorrTable);
     37  }
     38
     39  update_setphot(image, Nimage, &flatcorrTable);
    3240
    3341  // write image table

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

    r33103 r33104  
    22# include <glob.h>
    33
    4 int update_setphot (Image *image, off_t Nimage) {
     4int update_setphot (Image *image, off_t Nimage, FlatCorrectionTable *flatcorr) {
    55
    66  SkyRegion UserPatch;
    … …  
    5454    }
    5555
    56     update_catalog_setphot (&catalog, image, index, Nimage);
     56    update_catalog_setphot (&catalog, image, index, Nimage, flatcorr);
    5757
    5858    if (!UPDATE) {
    … …  
    7373void update_catalog_setphot (Catalog *catalog, Image *image, off_t *index, off_t Nimage, FlatCorrectionTable *flatcorr) {
    7474
    75   off_t i, j, m, id, idx, found;
    76   float Mcal, dMcal;
     75  off_t i, j, found;
    7776
    7877  found = 0;   
    7978  for (i = 0; i < catalog[0].Naverage; i++) {
    8079
    81     m = catalog[0].average[i].measureOffset;
     80    off_t m = catalog[0].average[i].measureOffset;
    8281    for (j = 0; j < catalog[0].average[i].Nmeasure; j++, m++) {
    83       idx = catalog[0].measure[m].imageID;
     82      off_t idx = catalog[0].measure[m].imageID;
    8483      if (idx <= 0) continue; // detections with imageID == 0 do not have a valid image (eg, ref photcode)
    8584
    86       id = index[idx];
    87       Mcal = image[id].Mcal;
    88       dMcal = image[id].dMcal;
    89       Mcal_offset = 0.0;
     85      off_t id = index[idx];
     86      float Mcal = image[id].Mcal;
     87      float dMcal = image[id].dMcal;
     88      float Mcal_offset = 0.0;
    9089
    9190      // 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;
     91      int flat_id = image[id].photom_map_id;
    9392      if (flat_id > 0) {
    9493          Mcal_offset = FlatCorrectionOffset (flatcorr, flat_id, catalog[0].measure[m].Xccd, catalog[0].measure[m].Yccd);
Note: See TracChangeset for help on using the changeset viewer.