IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 29714


Ignore:
Timestamp:
Nov 8, 2010, 6:20:37 PM (16 years ago)
Author:
eugene
Message:

add the ability to read a broken fits table (pad the missing areas); adjust dvorepair to count missing detections per image

Location:
branches/eam_branches/ipp-20101103/Ohana/src
Files:
3 added
11 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ipp-20101103/Ohana/src/addstar/src/ReadStarsFITS.c

    r29541 r29714  
    2222
    2323  /* load the table data */
    24   if (!gfits_fread_ftable_data (f, &table)) {
     24  if (!gfits_fread_ftable_data (f, &table, FALSE)) {
    2525    fprintf (stderr, "ERROR: can't read table header\n");
    2626    exit (1);

  • branches/eam_branches/ipp-20101103/Ohana/src/addstar/src/ReadStarsSDSS.c

    r28939 r29714  
    6060
    6161  /* load the table data */
    62   if (!gfits_fread_ftable_data (f, &table)) {
     62  if (!gfits_fread_ftable_data (f, &table, FALSE)) {
    6363    fprintf (stderr, "ERROR: can't read table header\n");
    6464    exit (1);

  • branches/eam_branches/ipp-20101103/Ohana/src/dvomerge/src/dvorepair.c

    r29694 r29714  
    44# define myAssert(LOGIC,MSG) { if (!LOGIC) { fprintf (stderr, "%s\n", MSG); abort(); } }
    55
    6 // broken cpt file, valid cpm file: we can recover everything in the cpt file from the cpm file:
    7 // * load the full cpm file
    8 // * loop over detections
    9 // * if ave_ref is new: create a new object
    10 //   * determine RA & DEC from ext_id
    11 //   * obj_id, cat_id are defined in detection
    12 //   *
     6// find images which are missing detections:
     7// * load the Images.dat file
     8// * create an index for imageID (seq = imageIDindex[i])
     9// * scan over all catalog files
     10// * load the cpm file (pad if short, identify the padded section)
     11// * loop over detections: increment detection count for each imageID
    1312
    1413int main (int argc, char **argv) {
    1514
    16   off_t Nmeasure, Nimage;
    17   int i, Nbytes, NaveMax, Naverage, NAVERAGE, Nave, Nold;
    18   int *found;
     15  off_t i, j, Nmeasure, Nvalid, Nimage, Nindex, *imageIdx, index, imageSeq;
     16  int Nbytes, NbytesPerRow, Nbad, Ncheck, Ntol;
     17  int *detCounts;
    1918
    20   Image *image, *thisImage;
    21   Average *average;
     19  Image *image;
    2220  Measure *measure;
    2321
    24   Matrix matrix;
     22  SkyTable *insky;
     23  SkyList *inlist;
    2524
    26   char *cpmFilename, *cptFilenameSrc, *cptFilenameTgt, *imageFilename;
     25  char *catdir = NULL;
     26  char *cpmFilename = NULL;
     27  char *imageFilename = NULL;
    2728
    28   Header cptHeaderPHU, cptHeaderTBL;
    29   Header cpmHeaderPHU, cpmHeaderTBL;
    30   FTable cpmFtable, cptFtable;
     29  Header cpmHeaderPHU;
     30  Header cpmHeaderTBL;
     31  FTable cpmFtable;
    3132
    32   FILE *cptFileSrc = NULL;
    33   FILE *cptFileTgt = NULL;
    3433  FILE *cpmFile = NULL;
    3534
    3635  char catformat;
    3736
    38   if (argc != 5) {
    39     fprintf (stderr, "USAGE: dvorepair (images) (cpm) (cptInput) (cptOutput)\n");
     37  if (argc != 3) {
     38    fprintf (stderr, "USAGE: dvorepair (catdir) (Ntol)\n");
     39    fprintf (stderr, "  catdir : database of interest\n");
     40    fprintf (stderr, "  Ntol : allow Ntol missing detections\n");
    4041    exit (2);
    4142  }
    4243
    43   imageFilename  = argv[1];
    44   cpmFilename    = argv[2];
    45   cptFilenameSrc = argv[3];
    46   cptFilenameTgt = argv[4];
    47 
    48   cpmFtable.header = &cpmHeaderTBL;
    49   cptFtable.header = &cptHeaderTBL;
    50 
    51   // XXX don't bother locking for now: this is totally manual..
     44  catdir = argv[1];
     45  Ntol = atoi(argv[2]);
    5246
    5347  // load the image data
     48  ALLOCATE(imageFilename, char, strlen(catdir) + 12);
     49  sprintf (imageFilename, "%s/Images.dat", catdir);
    5450  if ((image = LoadImages (imageFilename, &Nimage)) == NULL) return (FALSE);
    5551  BuildChipMatch (image, Nimage);
    5652
    57   // open cpm file
    58   cpmFile = fopen(cpmFilename, "r");
    59   myAssert(cpmFile, "failed to open cpm file");
     53  // generate an index for imageIDs:
     54  Nindex = 0;
     55  for (i = 0; i < Nimage; i++) {
     56    Nindex = MAX(image[i].imageID, Nindex);
     57  }
     58  ALLOCATE (imageIdx, off_t, (Nindex + 1));
     59  memset (imageIdx, 0, (Nindex + 1)*sizeof(off_t));
    6060
    61   // load the cpm header
    62   if (!gfits_fread_header (cpmFile, &cpmHeaderPHU)) {
    63     myAbort("failure to cpm header");
     61  for (i = 0; i < Nimage; i++) {
     62    index = image[i].imageID;
     63    if (index == 0) {
     64      fprintf (stderr, "?");
     65      continue;
     66    }
     67    if (imageIdx[index]) {
     68      fprintf (stderr, "!");
     69      continue;
     70    }
     71    imageIdx[index] = i;
    6472  }
    6573
    66   // move to TBL header
    67   Nbytes = cpmHeaderPHU.datasize + gfits_data_size (&cpmHeaderPHU);
    68   fseeko (cpmFile, Nbytes, SEEK_SET);
     74  // generate a list of the detection counts:
     75  ALLOCATE (detCounts, int, Nimage);
     76  memset (detCounts, 0, Nimage*sizeof(int));
    6977
    70   // read cpm TBL header
    71   if (!gfits_fread_header (cpmFile, &cpmHeaderTBL)) {
    72     myAbort("can't read header for cpm table");
     78  // load the sky table for the existing database
     79  insky = SkyTableLoadOptimal (catdir, NULL, NULL, FALSE, SKY_DEPTH_HST, VERBOSE);
     80  myAssert(insky, "can't read SkyTable");
     81  SkyTableSetFilenames (insky, catdir, "cpt");
     82
     83  // XXX apply this...generate the subset matching the user-selected region
     84  SkyRegion UserPatch;
     85  UserPatch.Rmin = 0.0;
     86  UserPatch.Rmax = 360.0;
     87  UserPatch.Dmin = -90.0;
     88  UserPatch.Dmax = +90.0;
     89  inlist = SkyListByPatch (insky, -1, &UserPatch);
     90
     91  // SkyListPopulatedRange (&Ns, &Ne, inlist, 0);
     92  // depth = inlist[0].regions[Ns][0].depth;
     93 
     94  ALLOCATE(cpmFilename, char, strlen(catdir) + 64);
     95
     96  // loop over the populated input regions
     97  Ncheck = 0;
     98  for (i = 0; i < inlist[0].Nregions; i++) {
     99    if (!inlist[0].regions[i][0].table) continue;
     100
     101    if (1) {
     102      sprintf (cpmFilename, "%s/%s.cpm", catdir, inlist[0].regions[i][0].name);
     103      cpmFtable.header = &cpmHeaderTBL;
     104
     105      // open cpm file
     106      cpmFile = fopen(cpmFilename, "r");
     107      if (!cpmFile) continue;
     108      // myAssert(cpmFile, "failed to open cpm file");
     109   
     110      // fprintf (stderr, "input: %s\n", inlist[0].regions[i][0].name);
     111
     112      // load the cpm header
     113      if (!gfits_fread_header (cpmFile, &cpmHeaderPHU)) {
     114        myAbort("failure to cpm header");
     115      }
     116
     117      // move to TBL header
     118      Nbytes = cpmHeaderPHU.datasize + gfits_data_size (&cpmHeaderPHU);
     119      fseeko (cpmFile, Nbytes, SEEK_SET);
     120
     121      // read cpm TBL header
     122      if (!gfits_fread_header (cpmFile, &cpmHeaderTBL)) {
     123        myAbort("can't read header for cpm table");
     124      }
     125
     126      // read Measure table data : format is irrelevant here */
     127      if (!gfits_fread_ftable_data (cpmFile, &cpmFtable, TRUE)) {
     128        myAbort("can't read data for cpm table");
     129      }
     130
     131      gfits_scan(&cpmHeaderTBL, "NAXIS1", "%d", 1, &NbytesPerRow);
     132
     133      measure = FtableToMeasure (&cpmFtable, &Nmeasure, &catformat);
     134      myAssert(measure, "failed to convert ftable to measure data");
     135   
     136      Nvalid = (int)(cpmFtable.validsize / NbytesPerRow);
     137      Nvalid = MIN(Nmeasure, Nvalid);
     138
     139      // examine all measurements and new objects as needed
     140      for (j = 0; j < Nvalid; j++) {
     141        index = measure[j].imageID;
     142        if (!index) {
     143          fprintf (stderr, "?");
     144          continue;
     145        }
     146
     147        imageSeq = imageIdx[index];
     148        // XXX check the range?
     149
     150        detCounts[imageSeq] ++;
     151      }
     152      fclose(cpmFile);
     153      gfits_free_header (&cpmHeaderPHU);
     154      gfits_free_header (&cpmHeaderTBL);
     155      free (measure);
     156
     157      Ncheck ++;
     158      if (Ncheck % 1000 == 0) {
     159        fprintf (stderr, "%s...", inlist[0].regions[i][0].name);
     160      }
     161    }
    73162  }
    74   // read Measure table data : format is irrelevant here */
    75   if (!gfits_fread_ftable_data (cpmFile, &cpmFtable)) {
    76     myAbort("can't read data for cpm table");
     163  fprintf (stderr, "\n");
     164
     165  Nbad = 0;
     166  for (i = 0; i < Nimage; i++) {
     167    if (detCounts[i] < image[i].nstar - Ntol) {
     168      fprintf (stdout, "image %s : %d vs %d\n", image[i].name, detCounts[i], image[i].nstar);
     169      Nbad ++;
     170    }
    77171  }
    78 
    79   measure = FtableToMeasure (&cpmFtable, &Nmeasure, &catformat);
    80   myAssert(measure, "failed to convert ftable to measure data");
    81 
    82   NaveMax = 0;
    83   NAVERAGE = 1000;
    84   ALLOCATE (average, Average, NAVERAGE);
    85   memset (average, 0, NAVERAGE*sizeof(Average));
    86 
    87   ALLOCATE (found, int, NAVERAGE);
    88   memset (found, 0, NAVERAGE*sizeof(int));
    89 
    90   // examine all measurements and new objects as needed
    91   for (i = 0; i < Nmeasure; i++) {
    92     Nave = measure[i].averef;
    93     if (found[Nave]) {
    94       average[Nave].Nmeasure ++;
    95       myAssert(average[Nave].objID == measure[i].objID, "objIDs do not match!");
    96       myAssert(average[Nave].catID == measure[i].catID, "catIDs do not match!");
    97       continue;
    98     }
    99 
    100     if (Nave >= NAVERAGE) {
    101       Nold = NAVERAGE;
    102       NAVERAGE = MAX(Nave + 1000, NAVERAGE + 1000);
    103       REALLOCATE (average, Average, NAVERAGE);
    104       memset (&average[Nold], 0, (NAVERAGE - Nold)*sizeof(Average));
    105 
    106       REALLOCATE (found, int, NAVERAGE);
    107       memset (&found[Nold], 0, (NAVERAGE - Nold)*sizeof(int));
    108     }
    109 
    110     NaveMax = MAX(Nave, NaveMax);
    111 
    112     found[Nave] = TRUE;
    113 
    114     // we are going to leave most of the elements of average unset: they are the result of
    115     // the relastro analysis for this object and can be recreated with a call to relastro
    116 
    117     // fields we have to set:
    118 
    119     // need to find image so we can use ccd coordinates to determine RA & DEC
    120     thisImage = MatchImage (image, Nimage, measure[i].t, measure[i].photcode, measure[i].imageID);
    121     XY_to_RD (&average[Nave].R, &average[Nave].D, measure[i].Xccd, measure[i].Yccd, &thisImage[0].coords);
    122 
    123     average[Nave].Nmeasure = 1;
    124     average[Nave].Nmissing = 0;
    125 
    126     // assume the resulting table set is unsorted
    127     average[Nave].measureOffset = -1;
    128     average[Nave].missingOffset = -1;
    129     average[Nave].extendOffset = -1;
    130 
    131     average[Nave].objID = measure[i].objID;
    132     average[Nave].catID = measure[i].catID;
    133     average[Nave].extID = CreatePSPSObjectID(average[Nave].R, average[Nave].D);
     172  if (!Nbad) {
     173    fprintf (stderr, "no bad images found\n");
    134174  }
    135   Naverage = NaveMax + 1;
    136 
    137   // have we created all objects in the range 0 - Naverage?
    138   for (i = 0; i < Naverage; i++) {
    139     myAssert(found[i], "failed to find one");
    140   }
    141 
    142   // open source cpt file
    143   cptFileSrc = fopen(cptFilenameSrc, "r");
    144   myAssert(cptFileSrc, "failed to open cpt file");
    145 
    146   // load the cpt header (use for CATID, RA,DEC range, filenames)
    147   if (!gfits_fread_header (cptFileSrc, &cptHeaderPHU)) {
    148     myAbort("failure to cpt header");
    149   }
    150 
    151   // update the output header
    152   gfits_modify (&cptHeaderPHU, "NSTARS",     "%d",      1,  Naverage);
    153   gfits_modify (&cptHeaderPHU, "NMEAS",      OFF_T_FMT, 1,  Nmeasure);
    154   gfits_modify (&cptHeaderPHU, "NMISS",      "%d",      1,  0);
    155   gfits_modify_alt (&cptHeaderPHU, "SORTED", "%t",      1,  FALSE);
    156 
    157   /* convert internal to external format */
    158   if (!AverageToFtable (&cptFtable, average, Naverage, catformat, NULL)) {
    159     myAbort("trouble converting format");
    160   }
    161 
    162   // create and write the output file
    163   cptFileTgt = fopen(cptFilenameTgt, "w");
    164   myAssert(cptFileTgt, "failed to open cpt file");
    165    
    166   // write PHU header
    167   if (!gfits_fwrite_header (cptFileTgt, &cptHeaderPHU)) {
    168     myAbort("can't write primary header");
    169   }
    170 
    171   // write the PHU matrix; this is probably a NOP, do I have to keep it in?
    172   gfits_create_matrix (&cptHeaderPHU, &matrix);
    173   if (!gfits_fwrite_matrix  (cptFileTgt, &matrix)) {
    174     myAbort("can't write primary matrix");
    175   }
    176   gfits_free_matrix (&matrix);
    177 
    178   // write the table data
    179   if (!gfits_fwrite_ftable_range (cptFileTgt, &cptFtable, 0, Naverage, 0, Naverage)) {
    180     myAbort("can't write table data");
    181   }
    182  
    183   fclose(cptFileTgt);
    184   fclose(cptFileSrc);
    185   fclose(cpmFile);
    186175
    187176  exit (0);

  • branches/eam_branches/ipp-20101103/Ohana/src/imregister/imphot/rfits.c

    r7080 r29714  
    1717    return (FALSE);
    1818  }
    19   if (!gfits_fread_ftable_data (db[0].f, &db[0].ftable)) {
     19  if (!gfits_fread_ftable_data (db[0].f, &db[0].ftable, FALSE)) {
    2020    fprintf (stderr, "can't read table data");
    2121    return (FALSE);

  • branches/eam_branches/ipp-20101103/Ohana/src/libdvo/src/dvo_catalog_mef.c

    r29001 r29714  
    6666  /* read Average table data (or skip) */
    6767  if (catalog[0].catflags & LOAD_AVES) {
    68     if (!gfits_fread_ftable_data (f, &ftable)) {
     68    if (!gfits_fread_ftable_data (f, &ftable, FALSE)) {
    6969      if (VERBOSE) fprintf (stderr, "can't read table average data");
    7070      return (FALSE);
    … …  
    9494  /* read Measure table data */
    9595  if (catalog[0].catflags & LOAD_MEAS) {
    96     if (!gfits_fread_ftable_data (f, &ftable)) {
     96    if (!gfits_fread_ftable_data (f, &ftable, FALSE)) {
    9797      if (VERBOSE) fprintf (stderr, "can't read table measure data");
    9898      return (FALSE);
    … …  
    119119  /* read Missing table data */
    120120  if (catalog[0].catflags & LOAD_MISS) {
    121     if (!gfits_fread_ftable_data (f, &ftable)) {
     121    if (!gfits_fread_ftable_data (f, &ftable, FALSE)) {
    122122      if (VERBOSE) fprintf (stderr, "can't read table missing data");
    123123      return (FALSE);
    … …  
    149149  /* read secfilt table data */
    150150  if (catalog[0].catflags & LOAD_SECF) {
    151     if (!gfits_fread_ftable_data (f, &ftable)) {
     151    if (!gfits_fread_ftable_data (f, &ftable, FALSE)) {
    152152      if (VERBOSE) fprintf (stderr, "can't read table secfilt data");
    153153      return (FALSE);

  • branches/eam_branches/ipp-20101103/Ohana/src/libdvo/src/dvo_catalog_split.c

    r29001 r29714  
    218218    }
    219219    /* read Average table data : format is irrelevant here */
    220     if (!gfits_fread_ftable_data (catalog[0].f, &ftable)) {
     220    if (!gfits_fread_ftable_data (catalog[0].f, &ftable, FALSE)) {
    221221      if (VERBOSE) fprintf (stderr, "can't read table average data");
    222222      return (FALSE);
    … …  
    249249    // XXX this allows an empty Measure catalog with non-empty Average catalog : is that OK?
    250250    /* read Measure table data */
    251     if (!gfits_fread_ftable_data (catalog[0].measure_catalog[0].f, &ftable)) {
     251    if (!gfits_fread_ftable_data (catalog[0].measure_catalog[0].f, &ftable, FALSE)) {
    252252      if (VERBOSE) fprintf (stderr, "can't read table measure data\n");
    253253      return (FALSE);
    … …  
    279279  if ((status != DVO_CAT_OPEN_EMPTY) && (catalog[0].catflags & LOAD_MISS)) {
    280280    /* read Missing table data */
    281     if (!gfits_fread_ftable_data (catalog[0].missing_catalog[0].f, &ftable)) {
     281    if (!gfits_fread_ftable_data (catalog[0].missing_catalog[0].f, &ftable, FALSE)) {
    282282      if (VERBOSE) fprintf (stderr, "can't read table missing data\n");
    283283      return (FALSE);
    … …  
    313313  if ((status != DVO_CAT_OPEN_EMPTY) && (catalog[0].catflags & LOAD_SECF)) {
    314314    /* read secfilt table data */
    315     if (!gfits_fread_ftable_data (catalog[0].secfilt_catalog[0].f, &ftable)) {
     315    if (!gfits_fread_ftable_data (catalog[0].secfilt_catalog[0].f, &ftable, FALSE)) {
    316316      if (VERBOSE) fprintf (stderr, "can't read table secfilt data\n");
    317317      return (FALSE);

  • branches/eam_branches/ipp-20101103/Ohana/src/libdvo/src/fits_db.c

    r28246 r29714  
    8383    return (FALSE);
    8484  }
    85   if (!gfits_fread_ftable_data (db[0].f, &db[0].ftable)) {
     85  if (!gfits_fread_ftable_data (db[0].f, &db[0].ftable, FALSE)) {
    8686    fprintf (stderr, "can't read table data");
    8787    gfits_db_free (db);

  • branches/eam_branches/ipp-20101103/Ohana/src/libfits/include/gfitsio.h

    r29537 r29714  
    6464  char                   *buffer;
    6565  off_t                   datasize;
     66  off_t                   validsize;
    6667} FTable;
    6768
    … …  
    172173int     gfits_define_table_column      PROTO((Header *header, char *format, char *label, char *comment, char *unit));
    173174int     gfits_fread_ftable             PROTO((FILE *f, FTable *ftable, char *extname));
    174 int     gfits_fread_ftable_data        PROTO((FILE *f, FTable *ftable));
     175int     gfits_fread_ftable_data        PROTO((FILE *f, FTable *ftable, int padIfShort));
    175176int     gfits_fread_ftable_range       PROTO((FILE *f, FTable *ftable, off_t start, off_t Nrows));
    176177int     gfits_fread_vtable             PROTO((FILE *f, VTable *vtable, char *extname, off_t Nrow, off_t *row));

  • branches/eam_branches/ipp-20101103/Ohana/src/libfits/table/F_read_T.c

    r28241 r29714  
    3636    gfits_scan (header, "EXTNAME", "%s", 1, tname);
    3737    if (!strcmp (tname, extname)) {
    38       if (gfits_fread_ftable_data (f, table)) return (TRUE);
     38      if (gfits_fread_ftable_data (f, table, FALSE)) return (TRUE);
    3939      gfits_free_header (header);
    4040      return (FALSE);
    … …  
    5050
    5151/*********************** fits read ftable data ***********************************/
    52 int gfits_fread_ftable_data (FILE *f, FTable *table) {
     52int gfits_fread_ftable_data (FILE *f, FTable *table, int padIfShort) {
    5353
    5454  off_t Nbytes, Nread;
    … …  
    6565    if (Nread < gfits_data_min_size (table[0].header)) {
    6666      fprintf (stderr, "error: fits read error in %s, read "OFF_T_FMT", need "OFF_T_FMT"\n", __func__,  Nread,  gfits_data_min_size (table[0].header));
    67       gfits_free_table  (table);
    68       return (FALSE);
     67      if (padIfShort) {
     68        memset (&table[0].buffer[Nread], 0, Nbytes - Nread);
     69        fprintf (stderr, "warning: file missing data, padding with zeros: USE AT YOUR OWN RISK!\n");
     70        table[0].validsize = Nread;
     71        table[0].datasize = Nbytes;
     72        return (TRUE);
     73      } else {
     74        gfits_free_table  (table);
     75        return (FALSE);
     76      }
    6977    }
    7078    fprintf (stderr, "warning: file missing pad\n");
    7179  }
     80  table[0].validsize = Nread;
    7281  table[0].datasize = Nbytes;
    7382  return (TRUE);

  • branches/eam_branches/ipp-20101103/Ohana/src/opihi/cmd.data/rd.c

    r28241 r29714  
    184184    ftable.header[0].buffer = NULL;
    185185    gfits_copy_header (&buf[0].header, ftable.header);
    186     status = gfits_fread_ftable_data (f, &ftable);  // this just reads the bytes (not even a SWAP)
     186    status = gfits_fread_ftable_data (f, &ftable, FALSE);  // this just reads the bytes (not even a SWAP)
    187187    status = gfits_uncompress_image (&buf[0].header, &buf[0].matrix, &ftable);
    188188    // uncompressing the image leaves the format as an extension

  • branches/eam_branches/ipp-20101103/Ohana/src/opihi/cmd.data/read_vectors.c

    r29540 r29714  
    160160
    161161  off_t Nbytes;
    162   int i, j, k, N, Nextend, Ny, Binary, vecType;
     162  int i, j, k, N, Nextend, Ny, Binary, vecType, padIfShort;
    163163  char type[16], ID[80], *CCDKeyword;
    164164  FTable table;
    … …  
    174174    CCDKeyword = strcreate (argv[N]);
    175175    remove_argument (N, &argc, argv);
     176  }
     177
     178  padIfShort = FALSE;
     179  if ((N = get_argument (argc, argv, "-pad-if-short"))) {
     180    remove_argument (N, &argc, argv);
     181    padIfShort = TRUE;
    176182  }
    177183
    … …  
    205211    }
    206212    if (!gfits_load_header (f, &header)) ESCAPE ("error reading header for extension");
    207     if (!gfits_fread_ftable_data (f, &table)) ESCAPE ("error reading table for extension");
     213    if (!gfits_fread_ftable_data (f, &table, padIfShort)) ESCAPE ("error reading table for extension");
    208214
    209215  } else {
    … …  
    236242        continue;
    237243      }
    238       if (!gfits_fread_ftable_data (f, &table)) ESCAPE ("error reading table for extension");
     244      if (!gfits_fread_ftable_data (f, &table, padIfShort)) ESCAPE ("error reading table for extension");
    239245      break;
    240246    }
Note: See TracChangeset for help on using the changeset viewer.