IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 29779


Ignore:
Timestamp:
Nov 14, 2010, 4:02:44 PM (15 years ago)
Author:
eugene
Message:

updates to dvorepair; add ability to plot mosaics by photcode (requires FindMosaicForImage to return the mosaic sequence number)

Location:
branches/eam_branches/ipp-20101103/Ohana/src
Files:
5 added
8 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ipp-20101103/Ohana/src/dvomerge/Makefile

    r29694 r29779  
    7272DVOREPAIR = \
    7373$(SRC)/dvorepair.$(ARCH).o \
     74$(SRC)/dvorepairFixCPT.$(ARCH).o \
     75$(SRC)/dvorepairImagesVsMeasures.$(ARCH).o \
     76$(SRC)/dvorepairDeleteImageList.$(ARCH).o \
     77$(SRC)/dvorepairFixImages.$(ARCH).o \
    7478$(SRC)/psps_ids.$(ARCH).o \
    7579$(SRC)/LoadImages.$(ARCH).o \
     80$(SRC)/ReadDeleteList.$(ARCH).o \
    7681$(SRC)/match_image.$(ARCH).o \
     82$(SRC)/help.$(ARCH).o \
    7783$(SRC)/ImageOps.$(ARCH).o
    7884

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

    r29694 r29779  
    1616# include <arpa/inet.h>
    1717# include <glob.h>
     18
     19# define myAbort(MSG) { fprintf (stderr, "%s\n", MSG); abort(); }
     20# define myAssert(LOGIC,MSG) { if (!(LOGIC)) { fprintf (stderr, "%s\n", MSG); abort(); } }
    1821
    1922int    VERBOSE;
    … …  
    8790int        dvo_update_image_IDs   PROTO((IDmapType *IDmap, Catalog *catalog));
    8891
    89 Image *LoadImages (char *filename, off_t *Nimage);
    90 Image *MatchImage (Image *image, off_t Nimage, unsigned int time, short int source, unsigned int imageID);
    91 off_t match_image (Image *image, off_t Nimage, unsigned int T, short int S);
     92// dvorepair prototypes
     93Image     *LoadImages             PROTO((FITS_DB *db, char *filename, off_t *Nimage));
     94Image     *MatchImage             PROTO((Image *image, off_t Nimage, unsigned int time, short int source, unsigned int imageID));
     95off_t      match_image            PROTO((Image *image, off_t Nimage, unsigned int T, short int S));
     96
     97int        SaveImages             PROTO((FITS_DB *oldDB, char *filename, Image *imagesOut, off_t Nout));
     98
     99int        dvorepairFixCPT        PROTO((int argc, char **argv));
     100int        dvorepairImagesVsMeasures PROTO((int argc, char **argv));
     101int        dvorepairDeleteImageList PROTO((int argc, char **argv));
     102int        dvorepairFixImages     PROTO((int argc, char **argv));
     103
     104int       *ReadDeleteList         PROTO((char *filename, int *nindex));
     105int        RepairTableCPT         PROTO((char *cptFilenameSrc, char *cptFilenameTgt, char *cpsFilenameSrc, char *cpsFilenameTgt, Measure *measure, off_t Nmeasure, Image *image, off_t Nimage, char catformat));
     106void       dvorepair_help         PROTO((int argc, char **argv));

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

    r29694 r29779  
    11# include "dvomerge.h"
    22
    3 int dvoUseImageCache = 1;
    4 
    5 static char *lastFilename = NULL;
    6 static time_t lastModified = 0;
    7 static Image *imageCache = NULL;
    8 static off_t cacheNimage = 0;
    9 
    10 static time_t getLastModified(char *filename);
    11 
    12 Image *LoadImages (char *filename, off_t *Nimage) {
     3Image *LoadImages (FITS_DB *db, char *filename, off_t *Nimage) {
    134
    145  int status;
    156  Image *image;
    16   FITS_DB db;
     7  double ZeroPoint;
    178 
    18   if (lastFilename) {
    19     if (dvoUseImageCache && !strcmp(lastFilename, filename)) {
    20       // Make sure the file hasn't changed since we loaded it
    21       if (getLastModified(filename) == lastModified) {
    22         *Nimage = cacheNimage;
    23         return  imageCache;
    24       }
    25     }
    26     free(lastFilename);
    27     lastFilename = NULL;
    28     free(imageCache);
    29     imageCache = NULL;
    30     cacheNimage = 0;
    31     lastModified = 0;
    32   }
     9  myAssert(filename, "programming errror");
     10  myAssert(Nimage, "programming errror");
     11   
     12  gfits_db_init (db);
     13  db->lockstate = LCK_SOFT;
     14  db->timeout   = 120.0;
    3315
    34 
    35   gfits_db_init (&db);
    36   db.lockstate = LCK_SOFT;
    37   db.timeout   = 120.0;
    38 
    39   if (!gfits_db_lock (&db, filename)) {
     16  if (!gfits_db_lock (db, filename)) {
    4017    fprintf (stderr, "error opening image catalog %s (1)\n", filename);
    4118    return (NULL);
    4219  }
    4320
    44   if (db.dbstate == LCK_EMPTY) {
     21  if (db->dbstate == LCK_EMPTY) {
    4522    fprintf (stderr, "note: image catalog is empty\n");
    4623    ALLOCATE (image, Image, 1);
    … …  
    4926  }
    5027
    51   status = dvo_image_load (&db, TRUE, FALSE);
    52   gfits_db_close (&db);
     28  status = dvo_image_load (db, TRUE, FALSE);
     29  gfits_db_close (db);
    5330
    5431  if (!status) {
    … …  
    5734  }
    5835
    59   image = gfits_table_get_Image (&db.ftable, Nimage, &db.swapped);
     36  image = gfits_table_get_Image (&db->ftable, Nimage, &db->swapped);
    6037  if (!image) {
    6138    fprintf (stderr, "ERROR: failed to read images\n");
    6239    return (NULL);
    6340  }
    64   if (dvoUseImageCache && image) {
    65     cacheNimage = *Nimage;
    66     imageCache = image;
    67     lastFilename = strdup(filename);
    68     lastModified = getLastModified(filename);
    69   }
     41
     42  gfits_scan (&db->header, "ZERO_PT", "%lf", 1, &ZeroPoint);
     43  SetZeroPoint (ZeroPoint);
    7044
    7145  return (image);
    7246}
    73 
    74 static time_t getLastModified(char *filename) {
    75   struct stat statbuf;
    76   if (!stat(filename, &statbuf)) {
    77     return statbuf.st_mtime;
    78   } else {
    79     return 0;
    80   }
    81 }
    82 
    83 void FreeImages(Image *images) {
    84   if (!dvoUseImageCache && (images != NULL)) {
    85     free(images);
    86   } else {
    87     // defer free until next LoadImages with a different or modified Images.dat
    88   }
    89 }
    90 

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

    r29735 r29779  
    11# include "dvomerge.h"
    2 
    3 # define myAbort(MSG) { fprintf (stderr, "%s\n", MSG); abort(); }
    4 # define myAssert(LOGIC,MSG) { if (!(LOGIC)) { fprintf (stderr, "%s\n", MSG); abort(); } }
    5 
    6 // delete images based on a text table of the target IDs
    7 // * load the delete table file (validate format)
    8 // * load the Images.dat file
    9 // * create an index for imageID (seq = imageIDindex[i])
    10 // * create an index of the imageIDs to be deleted (delete (T/F) = imageDELindex[i])
    11 // * determine the full RA & DEC range of the images to be deleted
    12 // * scan over all catalog files in the specified RA & DEC range
    13 // * load the cpm file (pad if short, identify the padded section)
    14 // * create a new cpm file
    15 // * loop over detections : only keep those not from images to be deleted
    16 // * load the cpt file
    17 // * rebuild the cpt file or delete the matching entries?
    18 
    19 int *ReadDeleteList(char *filename, int *nindex);
    20 int RepairTableCPT(char *cptFilenameSrc, char *cptFilenameTgt, char *cpsFilenameSrc, char *cpsFilenameTgt, Measure *measure, off_t Nmeasure, Image *image, off_t Nimage, char catformat);
    212
    223int main (int argc, char **argv) {
    234
    24   off_t i, j, Nmeasure, NmeasureNew, Ndelete, Nvalid, Nimage, Nindex, *imageIdx, index;
    25   int seq, Nbytes, NbytesPerRow, Ncheck, nPass, raPass;
    26   double Rthis, Dthis, Rmin, Rmax, Dmin, Dmax, Qthis, Qmin, Qmax, dR, dQ;
    27 
    28   Image *image;
    29   Measure *measure;
    30   Measure *measureNew;
    31 
    32   Matrix matrix;
    33 
    34   SkyTable *insky;
    35   SkyList *inlist;
    36 
    37   char *cpmFilenameSrc = NULL;
    38   char *cptFilenameSrc = NULL;
    39   char *cpsFilenameSrc = NULL;
    40   char *cpmFilenameTgt = NULL;
    41   char *cptFilenameTgt = NULL;
    42   char *cpsFilenameTgt = NULL;
    43 
    44   char *imageFilename = NULL;
    45 
    46   Header cpmHeaderPHU;
    47   Header cpmHeaderTBL;
    48   FTable cpmFtable;
    49 
    50   FILE *cpmFile = NULL;
    51 
    52   char catformat;
    53 
    54   if (argc != 3) {
    55     fprintf (stderr, "USAGE: dvorepair (catdir) (deleteList)\n");
    56     fprintf (stderr, "  catdir : database of interest\n");
    57     fprintf (stderr, "  deleteList : table from dvorepair giving images to be deleted\n");
    58     exit (2);
    59   }
    60 
    61   char *catdir = argv[1];
    62   char *delList = argv[2];
    63 
    64   // load the image data
    65   ALLOCATE(imageFilename, char, strlen(catdir) + 12);
    66   sprintf (imageFilename, "%s/Images.dat", catdir);
    67   if ((image = LoadImages (imageFilename, &Nimage)) == NULL) return (FALSE);
    68   BuildChipMatch (image, Nimage);
    69 
    70   // generate an index for imageIDs
    71 
    72   // find the full range of the imageID values
    73   Nindex = 0;
    74   for (i = 0; i < Nimage; i++) {
    75     Nindex = MAX(image[i].imageID, Nindex);
    76   }
    77 
    78   // create the index vector
    79   ALLOCATE (imageIdx, off_t, (Nindex + 1));
    80   memset (imageIdx, 0, (Nindex + 1)*sizeof(off_t));
    81 
    82   // assign the imageIDs to the index vector
    83   for (i = 0; i < Nimage; i++) {
    84     index = image[i].imageID;
    85     myAssert(index, "image index cannot be 0");
    86     myAssert(!imageIdx[index], "stomping on an earlier index");
    87     imageIdx[index] = i;
    88   }
    89 
    90   // generate an index of the deletion image
    91   int *deleteIndex;
    92   ALLOCATE (deleteIndex, int, Nindex + 1);
    93   memset (deleteIndex, 0, (Nindex + 1)*sizeof(int));
    94 
    95   int NdeleteIDs;
    96   int *deleteIDs = ReadDeleteList(delList, &NdeleteIDs);
    97 
    98   for (i = 0; i < NdeleteIDs; i++) {
    99     index = deleteIDs[i];
    100     myAssert(index > 0, "image index cannot be <= 0");
    101     myAssert(index <= Nindex, "delete index out of range of real data");
    102     myAssert(!deleteIndex[index], "stomping on an earlier index");
    103     deleteIndex[index] = TRUE;
    104   }
    105 
    106   // find the RA & DEC range of the images we want to delete
    107   Rmin = 360.0;
    108   Rmax = 0.0;
    109   Dmin = +90.0;
    110   Dmax = -90.0;
    111   for (i = 0; i < NdeleteIDs; i++) {
    112     index = deleteIDs[i];
    113     seq = imageIdx[index];
    114     if (!FindMosaicForImage (image, Nimage, seq)) {
    115       fprintf (stderr, "cannot find mosaic for %s, skipping\n", image[seq].name);
    116       continue;
    117     }
    118 
    119     XY_to_RD(&Rthis, &Dthis, 0, 0, &image[seq].coords);
    120     Rmin = MIN(Rthis, Rmin);
    121     Rmax = MAX(Rthis, Rmax);
    122     Dmin = MIN(Dthis, Dmin);
    123     Dmax = MAX(Dthis, Dmax);
    124     Qthis = (Rthis < 180.0) ? Rthis : Rthis - 360.0;
    125     Qmin = MIN(Qthis, Qmin);
    126     Qmax = MAX(Qthis, Qmax);
    127 
    128     // fprintf (stderr, "image @ %f,%f\n", Rthis, Dthis);
    129 
    130     XY_to_RD(&Rthis, &Dthis, image[seq].NX, 0, &image[seq].coords);
    131     Rmin = MIN(Rthis, Rmin);
    132     Rmax = MAX(Rthis, Rmax);
    133     Dmin = MIN(Dthis, Dmin);
    134     Dmax = MAX(Dthis, Dmax);
    135     Qthis = (Rthis < 180.0) ? Rthis : Rthis - 360.0;
    136     Qmin = MIN(Qthis, Qmin);
    137     Qmax = MAX(Qthis, Qmax);
    138 
    139     XY_to_RD(&Rthis, &Dthis, 0, image[seq].NY, &image[seq].coords);
    140     Rmin = MIN(Rthis, Rmin);
    141     Rmax = MAX(Rthis, Rmax);
    142     Dmin = MIN(Dthis, Dmin);
    143     Dmax = MAX(Dthis, Dmax);
    144     Qthis = (Rthis < 180.0) ? Rthis : Rthis - 360.0;
    145     Qmin = MIN(Qthis, Qmin);
    146     Qmax = MAX(Qthis, Qmax);
    147 
    148     XY_to_RD(&Rthis, &Dthis, image[seq].NX, image[seq].NY, &image[seq].coords);
    149     Rmin = MIN(Rthis, Rmin);
    150     Rmax = MAX(Rthis, Rmax);
    151     Dmin = MIN(Dthis, Dmin);
    152     Dmax = MAX(Dthis, Dmax);
    153     Qthis = (Rthis < 180.0) ? Rthis : Rthis - 360.0;
    154     Qmin = MIN(Qthis, Qmin);
    155     Qmax = MAX(Qthis, Qmax);
    156   }
    157 
    158   dQ = Qmax - Qmin;
    159   dR = Rmax - Rmin;
    160 
    161   // load the sky table for the existing database
    162   insky = SkyTableLoadOptimal (catdir, NULL, NULL, FALSE, SKY_DEPTH_HST, VERBOSE);
    163   myAssert(insky, "can't read SkyTable");
    164   SkyTableSetFilenames (insky, catdir, "cpt");
     5  dvorepair_help(argc, argv);
    1656 
    166   // XXX apply this...generate the subset matching the user-selected region
    167   SkyRegion UserPatch;
    168 
    169   if (dR < dQ) {
    170     fprintf (stderr, "R,D range: %f - %f, %f - %f\n", Rmin, Rmax, Dmin, Dmax);
    171     nPass = 1;
    172   } else {
    173     fprintf (stderr, "R,D range: %f - %f, %f - %f\n", Qmin, Qmax, Dmin, Dmax);
    174     nPass = 2;
    175   }
    176  
    177   ALLOCATE(cpmFilenameSrc, char, strlen(catdir) + 64);
    178   ALLOCATE(cptFilenameSrc, char, strlen(catdir) + 64);
    179   ALLOCATE(cpsFilenameSrc, char, strlen(catdir) + 64);
    180   ALLOCATE(cpmFilenameTgt, char, strlen(catdir) + 64);
    181   ALLOCATE(cptFilenameTgt, char, strlen(catdir) + 64);
    182   ALLOCATE(cpsFilenameTgt, char, strlen(catdir) + 64);
    183 
    184   for (raPass = 0; raPass < nPass; raPass++) {
    185    
    186     UserPatch.Dmin = Dmin + 0.1;
    187     UserPatch.Dmax = Dmax + 0.1;
    188 
    189     if (dR < dQ) {
    190       UserPatch.Rmin = Rmin;
    191       UserPatch.Rmax = Rmax;
    192     } else {
    193       if (raPass == 0) {
    194         UserPatch.Rmin = 0.0;
    195         UserPatch.Rmax = Qmax;
    196       } else {
    197         UserPatch.Rmin = Qmin + 360.0;
    198         UserPatch.Rmax = 360.0;
    199       }
    200     }
    201 
    202     inlist = SkyListByPatch (insky, -1, &UserPatch);
    203  
    204     // SkyListPopulatedRange (&Ns, &Ne, inlist, 0);
    205     // depth = inlist[0].regions[Ns][0].depth;
    206  
    207     // loop over the populated input regions
    208     Ncheck = 0;
    209     for (i = 0; i < inlist[0].Nregions; i++) {
    210       if (!inlist[0].regions[i][0].table) continue;
    211 
    212       sprintf (cpmFilenameSrc, "%s/%s.cpm", catdir, inlist[0].regions[i][0].name);
    213       sprintf (cptFilenameSrc, "%s/%s.cpt", catdir, inlist[0].regions[i][0].name);
    214       sprintf (cpsFilenameSrc, "%s/%s.cps", catdir, inlist[0].regions[i][0].name);
    215       sprintf (cpmFilenameTgt, "%s/%s.cpm.fixed", catdir, inlist[0].regions[i][0].name);
    216       sprintf (cptFilenameTgt, "%s/%s.cpt.fixed", catdir, inlist[0].regions[i][0].name);
    217       sprintf (cpsFilenameTgt, "%s/%s.cps.fixed", catdir, inlist[0].regions[i][0].name);
    218       cpmFtable.header = &cpmHeaderTBL;
    219 
    220       /*** read and examine the CPM file ***/
    221       {
    222         // fprintf (stderr, "check %s\n", cpmFilenameSrc);
    223 
    224         // open cpm file
    225         cpmFile = fopen(cpmFilenameSrc, "r");
    226         if (!cpmFile) continue;
    227         // myAssert(cpmFile, "failed to open cpm file");
    228    
    229         // load the cpm header
    230         if (!gfits_fread_header (cpmFile, &cpmHeaderPHU)) {
    231           myAbort("failure to cpm header");
    232         }
    233 
    234         // move to TBL header
    235         Nbytes = cpmHeaderPHU.datasize + gfits_data_size (&cpmHeaderPHU);
    236         fseeko (cpmFile, Nbytes, SEEK_SET);
    237 
    238         // read cpm TBL header
    239         if (!gfits_fread_header (cpmFile, &cpmHeaderTBL)) {
    240           myAbort("can't read header for cpm table");
    241         }
    242 
    243         // read Measure table data : format is irrelevant here */
    244         if (!gfits_fread_ftable_data (cpmFile, &cpmFtable, TRUE)) {
    245           myAbort("can't read data for cpm table");
    246         }
    247 
    248         gfits_scan(&cpmHeaderTBL, "NAXIS1", "%d", 1, &NbytesPerRow);
    249 
    250         measure = FtableToMeasure (&cpmFtable, &Nmeasure, &catformat);
    251         myAssert(measure, "failed to convert ftable to measure data");
    252    
    253         Nvalid = (int)(cpmFtable.validsize / NbytesPerRow);
    254         Nvalid = MIN(Nmeasure, Nvalid);
    255 
    256         // close the input cpm file
    257         fclose(cpmFile);
    258 
    259         // allocate an output array of measures (to replace, if needed)
    260         ALLOCATE (measureNew, Measure, Nvalid);
    261 
    262         NmeasureNew = 0;
    263         Ndelete = 0;
    264 
    265         // examine all measurements: find ones that need to be deleted
    266         for (j = 0; j < Nvalid; j++) {
    267           index = measure[j].imageID;
    268           myAssert(index, "measure is missing an image ID");
    269 
    270           if (deleteIndex[index]) {
    271             Ndelete ++;
    272             continue;
    273           }
    274           measureNew[NmeasureNew] = measure[j];
    275           NmeasureNew ++;
    276         }
    277       }
    278 
    279       fprintf (stderr, "deleting %d of %d (keep %d) detections from %s -> %s\n", (int) Ndelete, (int) Nvalid, (int) NmeasureNew, cpmFilenameSrc, cpmFilenameTgt);
    280 
    281       // if we actually want to delete any measurements, write out a new cpm file
    282       if (Ndelete > 0) {
    283 
    284         // the CPT and CPS tables need to be regenerated.  This must happen first because, in the process, we also update measure->averef
    285         RepairTableCPT(cptFilenameSrc, cptFilenameTgt, cpsFilenameSrc, cpsFilenameTgt, measureNew, NmeasureNew, image, Nimage, catformat);
    286 
    287         // convert internal to external format
    288         if (!MeasureToFtable (&cpmFtable, measureNew, NmeasureNew, catformat)) {
    289           myAbort("trouble converting format");
    290         }
    291 
    292         // create and write the output file
    293         cpmFile = fopen(cpmFilenameTgt, "w");
    294         myAssert(cpmFile, "failed to open cpt file");
    295        
    296         // write PHU header
    297         if (!gfits_fwrite_header (cpmFile, &cpmHeaderPHU)) {
    298           myAbort("can't write primary header");
    299         }
    300 
    301         // write the PHU matrix; this is probably a NOP, do I have to keep it in?
    302         gfits_create_matrix (&cpmHeaderPHU, &matrix);
    303         if (!gfits_fwrite_matrix  (cpmFile, &matrix)) {
    304           myAbort("can't write primary matrix");
    305         }
    306         gfits_free_matrix (&matrix);
    307        
    308         // write the table data
    309         if (!gfits_fwrite_ftable_range (cpmFile, &cpmFtable, 0, NmeasureNew, 0, NmeasureNew)) {
    310           myAbort("can't write table data");
    311         }
    312         fclose (cpmFile);
    313       }
    314       gfits_free_header (&cpmHeaderPHU);
    315       gfits_free_header (&cpmHeaderTBL);
    316       free (measure);
    317       free (measureNew);
    318 
    319       Ncheck ++;
    320       if (Ncheck % 1000 == 0) {
    321         fprintf (stderr, "%s...", inlist[0].regions[i][0].name);
    322       }
    323     }
    324     fprintf (stderr, "\n");
    325     SkyListFree(inlist);
    326   }
    327 
    328   free (imageFilename);
    329   free (imageIdx);
    330   free (deleteIndex);
    331   free (deleteIDs);
    332 
    333   free(cpmFilenameSrc);
    334   free(cptFilenameSrc);
    335   free(cpsFilenameSrc);
    336   free(cpmFilenameTgt);
    337   free(cptFilenameTgt);
    338   free(cpsFilenameTgt);
    339 
    340   SkyTableFree(insky);
    341 
    342   exit (0);
     7  if (!strcmp(argv[1], "-fix-cpt")) dvorepairFixCPT(argc, argv);
     8  if (!strcmp(argv[1], "-images-vs-measures")) dvorepairImagesVsMeasures(argc, argv);
     9  if (!strcmp(argv[1], "-delete-image-list")) dvorepairDeleteImageList(argc, argv);
     10  if (!strcmp(argv[1], "-fix-images")) dvorepairFixImages(argc, argv);
     11  dvorepair_help(0, NULL);
     12  exit (2);
    34313}
    344 
    345 int *ReadDeleteList(char *filename, int *nindex) {
    346 
    347   int j, index, Nindex, NINDEX, *indexList, Nline;
    348   int Nstart, Nbytes, Nread, status;
    349   char *c0, *c1, *space, *indexPoint, *buffer;
    350 
    351   FILE *f = fopen (filename, "r");
    352   myAssert(f, "failed to open delete list");
    353  
    354   Nindex = 0;
    355   NINDEX = 1000;
    356   ALLOCATE(indexList, int, NINDEX);
    357 
    358   ALLOCATE (buffer, char, 0x10001);
    359   bzero (buffer, 0x10001);
    360 
    361   Nstart = 0;
    362   Nline = 0;
    363   while (1) {
    364     Nbytes = 0x10000 - Nstart;
    365     bzero (&buffer[Nstart], Nbytes);
    366     Nread = fread (&buffer[Nstart], 1, Nbytes, f);
    367     if (ferror (f)) {
    368       perror ("error reading data file");
    369       exit (1);
    370     }
    371     if (Nread == 0) break; // end of data
    372     // nbytes = Nread + Nstart;
    373    
    374     status = TRUE;
    375     c0 = buffer;
    376     while (status) {
    377       // identify the range of the line
    378       c1 = strchr (c0, '\n');
    379       if (c1 == (char *) NULL) {
    380         Nstart = strlen (c0);
    381         memmove (buffer, c0, Nstart);
    382         break;
    383       } else {
    384         *c1 = 0;
    385       }     
    386       if (*c0 == '#') goto next_line;
    387 
    388       // confirm 'image' as first token:
    389       if (strncmp(c0, "image", 5)) {
    390         fprintf (stderr, "error on line %d: missing 'image' token\n", Nline);
    391         goto next_line;
    392       }
    393 
    394       // confirm we have 7 fields broken by 6 spaces:
    395       space = c0;
    396       for (j = 0; j < 6; j++) {
    397         space = strchr(space, ' ');
    398         if (!space) {
    399           fprintf (stderr, "line only has %d word(s)\n", j + 1);
    400           goto next_line;
    401         }
    402         space ++;
    403         if (j == 1) {
    404           indexPoint = space + 1;
    405         }
    406       }
    407 
    408       // save the value we have found:
    409       index = atoi(indexPoint);
    410       indexList[Nindex] = index;
    411 
    412       // fprintf (stderr, "index: %d, line: %s\n", index, c0);
    413 
    414       Nindex ++;
    415       if (Nindex == NINDEX) {
    416         NINDEX += 1000;
    417         REALLOCATE (indexList, int, NINDEX);
    418       }
    419      
    420     next_line:
    421       Nline ++;
    422       c0 = c1 + 1;
    423     }
    424   }
    425 
    426   *nindex = Nindex;
    427   return (indexList);
    428 }
    429 
    430 /* Delete List Format:
    431    image o5252g0035o.140669.cm.51017.smf[XY11.hdr] (219974) : 368 vs 375
    432    image o5252g0051o.140685.cm.51123.smf[XY11.hdr] (226417) : 350 vs 356
    433    image o5399g0207o.194950.cm.98351.smf[XY46.hdr] (1683277) : 10341 vs 10347
    434    image o5464g0284o.230438.cm.123774.smf[XY21.hdr] (2634913) : 214 vs 242
    435    image o5464g0305o.230459.cm.123795.smf[XY21.hdr] (2635758) : 232 vs 257
    436    image o5465g0306o.231171.cm.124201.smf[XY01.hdr] (2654941) : 0 vs 494
    437    image o5465g0306o.231171.cm.124201.smf[XY02.hdr] (2654942) : 10 vs 813
    438    image o5465g0306o.231171.cm.124201.smf[XY10.hdr] (2654947) : 0 vs 534
    439    image o5465g0306o.231171.cm.124201.smf[XY11.hdr] (2654948) : 0 vs 447
    440 */
    441 
    442 int RepairTableCPT(char *cptFilenameSrc, char *cptFilenameTgt, char *cpsFilenameSrc, char *cpsFilenameTgt, Measure *measure, off_t Nmeasure, Image *image, off_t Nimage, char catformat) {
    443 
    444   off_t *averefMatch;
    445   off_t i, NaveMax, Naverage, NAVERAGE, NaverageOut, Nave, Nout, Nold;
    446   int *found, Nsecfilt;
    447 
    448   Image *thisImage;
    449   Average *average, *averageOut;
    450 
    451   Matrix matrix;
    452 
    453   Header cptHeaderPHU;
    454   Header cptHeaderTBL;
    455   FTable cptFtable;
    456 
    457   FILE *cptFileSrc = NULL;
    458   FILE *cptFileTgt = NULL;
    459 
    460   cptFtable.header = &cptHeaderTBL;
    461 
    462   NaveMax = 0;
    463   NAVERAGE = 1000;
    464   ALLOCATE (average, Average, NAVERAGE);
    465   memset (average, 0, NAVERAGE*sizeof(Average));
    466 
    467   ALLOCATE (found, int, NAVERAGE);
    468   memset (found, 0, NAVERAGE*sizeof(int));
    469 
    470   // examine all measurements and new objects as needed
    471   for (i = 0; i < Nmeasure; i++) {
    472     Nave = measure[i].averef;
    473 
    474     if (Nave >= NAVERAGE) {
    475       Nold = NAVERAGE;
    476       NAVERAGE = MAX(Nave + 1000, NAVERAGE + 1000);
    477       REALLOCATE (average, Average, NAVERAGE);
    478       memset (&average[Nold], 0, (NAVERAGE - Nold)*sizeof(Average));
    479 
    480       REALLOCATE (found, int, NAVERAGE);
    481       memset (&found[Nold], 0, (NAVERAGE - Nold)*sizeof(int));
    482     }
    483 
    484     if (found[Nave]) {
    485       average[Nave].Nmeasure ++;
    486       myAssert(average[Nave].objID == measure[i].objID, "objIDs do not match!");
    487       myAssert(average[Nave].catID == measure[i].catID, "catIDs do not match!");
    488       continue;
    489     }
    490 
    491     NaveMax = MAX(Nave, NaveMax);
    492 
    493     found[Nave] = TRUE;
    494 
    495     // we are going to leave most of the elements of average unset: they are the result of
    496     // the relastro analysis for this object and can be recreated with a call to relastro
    497 
    498     // need to find image so we can use ccd coordinates to determine RA & DEC
    499     thisImage = MatchImage (image, Nimage, measure[i].t, measure[i].photcode, measure[i].imageID);
    500     XY_to_RD (&average[Nave].R, &average[Nave].D, measure[i].Xccd, measure[i].Yccd, &thisImage[0].coords);
    501 
    502     average[Nave].Nmeasure = 1;
    503     average[Nave].Nmissing = 0;
    504 
    505     // assume the resulting table set is unsorted
    506     average[Nave].measureOffset = -1;
    507     average[Nave].missingOffset = -1;
    508     average[Nave].extendOffset = -1;
    509 
    510     average[Nave].objID = measure[i].objID;
    511     average[Nave].catID = measure[i].catID;
    512     average[Nave].extID = CreatePSPSObjectID(average[Nave].R, average[Nave].D);
    513   }
    514   Naverage = NaveMax + 1;
    515 
    516   // we now have an average table, but there will be holes due to deleted measurements
    517   // create a new average table with only existing entries
    518 
    519   ALLOCATE (averageOut, Average, Naverage);
    520   memset (averageOut, 0, Naverage*sizeof(Average));
    521 
    522   ALLOCATE (averefMatch, off_t, Naverage);
    523   memset (averefMatch, 0, Naverage*sizeof(int));
    524 
    525   Nave = 0;
    526   for (i = 0; i < Naverage; i++) {
    527     if (!found[i]) continue;
    528     averageOut[Nave] = average[i]; // use a memcpy?
    529     averefMatch[i] = Nave;
    530     Nave ++;
    531   }
    532   NaverageOut = Nave;
    533 
    534   for (i = 0; i < Nmeasure; i++) {
    535     Nave = measure[i].averef;
    536     Nout = averefMatch[Nave];
    537     myAssert(Nout < NaverageOut, "output averef is wrong");
    538    
    539     myAssert(average[Nave].objID == measure[i].objID, "objIDs do not match");
    540     myAssert(average[Nave].catID == measure[i].catID, "objIDs do not match");
    541     myAssert(averageOut[Nout].objID == measure[i].objID, "objIDs do not match");
    542     myAssert(averageOut[Nout].catID == measure[i].catID, "objIDs do not match");
    543 
    544     measure[i].averef = Nout;
    545   }
    546 
    547   fprintf (stderr, "cpt file : %d obj -> %d obj (%s -> %s)\n", (int) Naverage, (int) NaverageOut, cptFilenameSrc, cptFilenameTgt);
    548 
    549   // open source cpt file
    550   cptFileSrc = fopen(cptFilenameSrc, "r");
    551   myAssert(cptFileSrc, "failed to open cpt file");
    552 
    553   // load the cpt header (use for CATID, RA,DEC range, filenames)
    554   if (!gfits_fread_header (cptFileSrc, &cptHeaderPHU)) {
    555     myAbort("failure to cpt header");
    556   }
    557 
    558   // update the output header
    559   gfits_modify (&cptHeaderPHU, "NSTARS",     OFF_T_FMT, 1,  NaverageOut);
    560   gfits_modify (&cptHeaderPHU, "NMEAS",      OFF_T_FMT, 1,  Nmeasure);
    561   gfits_modify (&cptHeaderPHU, "NMISS",      "%d",      1,  0);
    562   gfits_modify_alt (&cptHeaderPHU, "SORTED", "%t",      1,  FALSE);
    563 
    564   gfits_scan (&cptHeaderPHU, "NSECFILT",     "%d",      1,  &Nsecfilt);
    565 
    566   /* convert internal to external format */
    567   if (!AverageToFtable (&cptFtable, averageOut, NaverageOut, catformat, NULL)) {
    568     myAbort("trouble converting format");
    569   }
    570 
    571   // create and write the output file
    572   cptFileTgt = fopen(cptFilenameTgt, "w");
    573   myAssert(cptFileTgt, "failed to open cpt file");
    574    
    575   // write PHU header
    576   if (!gfits_fwrite_header (cptFileTgt, &cptHeaderPHU)) {
    577     myAbort("can't write primary header");
    578   }
    579 
    580   // write the PHU matrix; this is probably a NOP, do I have to keep it in?
    581   gfits_create_matrix (&cptHeaderPHU, &matrix);
    582   if (!gfits_fwrite_matrix  (cptFileTgt, &matrix)) {
    583     myAbort("can't write primary matrix");
    584   }
    585   gfits_free_matrix (&matrix);
    586 
    587   // write the table data
    588   if (!gfits_fwrite_ftable_range (cptFileTgt, &cptFtable, 0, NaverageOut, 0, NaverageOut)) {
    589     myAbort("can't write table data");
    590   }
    591  
    592   fclose(cptFileTgt);
    593   fclose(cptFileSrc);
    594 
    595   gfits_free_header (&cptHeaderPHU);
    596   gfits_free_header (&cptHeaderTBL);
    597   free (average);
    598   free (averageOut);
    599 
    600   free (found);
    601   free (averefMatch);
    602 
    603   {
    604     Header cpsHeaderPHU;
    605     Header cpsHeaderTBL;
    606     FTable cpsFtable;
    607 
    608     FILE *cpsFileSrc = NULL;
    609     FILE *cpsFileTgt = NULL;
    610 
    611     SecFilt *secfilt = NULL;
    612 
    613     cpsFtable.header = &cpsHeaderTBL;
    614 
    615     // open source cpt file
    616     cpsFileSrc = fopen(cpsFilenameSrc, "r");
    617     myAssert(cpsFileSrc, "failed to open cps file");
    618    
    619     // load the cps header (use for CATID, RA,DEC range, filenames)
    620     if (!gfits_fread_header (cpsFileSrc, &cpsHeaderPHU)) {
    621       myAbort("failure to cps header");
    622     }
    623 
    624     int Nrows = Nsecfilt*NaverageOut;
    625     ALLOCATE (secfilt, SecFilt, Nrows);
    626 
    627     /* convert internal to external format */
    628     if (!SecFiltToFtable (&cpsFtable, secfilt, Nrows, catformat)) {
    629       myAbort("trouble converting format");
    630     }
    631 
    632     // create and write the output file
    633     cpsFileTgt = fopen(cpsFilenameTgt, "w");
    634     myAssert(cpsFileTgt, "failed to open cps file");
    635    
    636     // write PHU header
    637     if (!gfits_fwrite_header (cpsFileTgt, &cpsHeaderPHU)) {
    638       myAbort("can't write primary header");
    639     }
    640 
    641     // write the PHU matrix; this is probably a NOP, do I have to keep it in?
    642     gfits_create_matrix (&cpsHeaderPHU, &matrix);
    643     if (!gfits_fwrite_matrix  (cpsFileTgt, &matrix)) {
    644       myAbort("can't write primary matrix");
    645     }
    646     gfits_free_matrix (&matrix);
    647 
    648     // write the table data
    649     if (!gfits_fwrite_ftable_range (cpsFileTgt, &cpsFtable, 0, Nrows, 0, Nrows)) {
    650       myAbort("can't write table data");
    651     }
    652  
    653     fclose(cpsFileTgt);
    654     fclose(cpsFileSrc);
    655 
    656     gfits_free_header (&cpsHeaderPHU);
    657     gfits_free_header (&cpsHeaderTBL);
    658     free (secfilt);
    659   }
    660 
    661   return (TRUE);
    662 }
    663 

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

    r29001 r29779  
    8686}
    8787
     88void dvorepair_help (int argc, char **argv) {
     89
     90  /* check for help request */
     91  if (!argv) goto show_help;
     92  if (get_argument (argc, argv, "-help")) goto show_help;
     93  if (get_argument (argc, argv, "-h"))    goto show_help;
     94  if (argc < 2) goto show_help;
     95  return;
     96
     97show_help:
     98
     99  fprintf (stderr, "USAGE\n");
     100  fprintf (stderr, "  dvorepair -fix-cpt (images) (cpm) (cptInput) (cptOutput) - regenerate a specific cpt file from the cpm file\n");
     101  fprintf (stderr, "  dvorepair -images-vs-measures (catdir) (Ntol) - find images with too many missing detections\n");
     102  fprintf (stderr, "  dvorepair -delete-image-list (catdir) (deleteList) - delete a set of images based on image IDs (output from -images-vs-measures)\n");
     103  fprintf (stderr, "  dvorepair -fix-images (catdir) (deleteList) - delete a set of images based on image IDs (output from -images-vs-measures)\n");
     104  fprintf (stderr, "\n");
     105
     106  fprintf (stderr, "  optional flags:\n");
     107  fprintf (stderr, "  -help                       : this list\n");
     108  fprintf (stderr, "  -h                          : this list\n\n");
     109  exit (2);
     110}
     111

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

    r29001 r29779  
    293293int isRegisteredMosaic (void);
    294294off_t GetRegisteredMosaic (void);
     295off_t *GetChipMatch ();
    295296int GetMosaicCoords (Coords *coords);
    296297int FindMosaicForImage (Image *images, off_t Nimages, off_t entry);

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

    r27435 r29779  
    1717}
    1818
     19/* what is the currently registered mosaic image? */
     20off_t *GetChipMatch () {
     21  return (ChipMatch);
     22}
     23
    1924/* given an image array and a current entry of type WRP, find matching DIS & Register it */
    2025int FindMosaicForImage (Image *images, off_t Nimages, off_t entry) {
    … …  
    4045  if (strcmp(&images[entry].coords.ctype[4], "-WRP")) {
    4146    /* not a wrp image, do nothing */
    42     return (TRUE); /* error or not */
     47    return (entry + 1); /* error or not */
    4348  }
    4449
    … …  
    5257    RegisterMosaic (&images[i].coords);
    5358    iDIS = i;
    54     return (TRUE);
     59    return (i + 1);
    5560  }
    5661
    … …  
    6267    RegisterMosaic (&images[i].coords);
    6368    iDIS = i;
    64     return (TRUE);
     69    return (i + 1);
    6570  }
    6671  return (FALSE);
    … …  
    7479  if (strcmp(&images[entry].coords.ctype[4], "-WRP")) {
    7580    /* not a wrp image, do nothing */
    76     return (TRUE); /* error or not? */
     81    return (entry + 1); /* error or not? */
    7782  }
    7883
    … …  
    8590  RegisterMosaic (&images[N].coords);
    8691  iDIS = N;
    87   return (TRUE);
     92  return (N + 1);
    8893}
    8994
    … …  
    165170  SortDISindex (DIStzero, DISentry, Ndis);
    166171
     172  // ChipMatch has a few possible values:
     173  // -3 : image is a mosaic (DIS)
     174  // -2 : image is an unassigned WRP image
     175  // -1 : image is not a WRP or DIS images
     176  // >= 0 : value is the mosaic seq number for this WRP image
     177
    167178  /* find all matched WRP images */
    168179  ALLOCATE (ChipMatch, off_t, Nimages);
    169180  for (i = 0; i < Nimages; i++) {
    170181    ChipMatch[i] = -1;
     182    if (!strcmp(&images[i].coords.ctype[4], "-DIS")) {
     183      ChipMatch[i] = -3;
     184      continue;
     185    }
    171186    if (strcmp(&images[i].coords.ctype[4], "-WRP")) continue;
    172187

  • branches/eam_branches/ipp-20101103/Ohana/src/opihi/dvo/images.c

    r28958 r29779  
    1111int images (int argc, char **argv) {
    1212
    13   off_t i, Nimage;
     13  off_t i, Nimage, Nmosaic;
    1414  int j, status, InPic, leftside, *plist, TimeSelect, ByName;
    1515  int WITH_MOSAIC, SOLO_MOSAIC, HIDDEN;
    1616  time_t tzero, tend;
    17   int N, NPTS, n, npts, Npts, kapa;
     17  int N, NPTS, n, npts, Npts, kapa, *foundMosaic;
    1818  Vector Xvec, Yvec;
    1919  double r[8], d[8], x[8], y[8], Rmin, Rmax, Rmid, trange, Radius;
    … …  
    132132  BuildChipMatch (image, Nimage);
    133133
     134  if (SOLO_MOSAIC && photcode) {
     135    ALLOCATE(foundMosaic, int, Nimage);
     136    memset(foundMosaic, 0, Nimage*sizeof(int));
     137  }
     138
    134139  Rmin = graphmode.coords.crval1 - 180.0;
    135140  Rmax = graphmode.coords.crval1 + 180.0;
    … …  
    137142 
    138143  int DistortImage = wordhash ("-DIS");
     144  int ChipImage    = wordhash ("-WRP");
    139145  int TriangleUp   = wordhash ("TRP-");
    140146  int TriangleDn   = wordhash ("TRM-");
    … …  
    150156    if (ByName && strncmp (image[i].name, name, strlen(name))) continue;
    151157    if (TimeSelect && ((image[i].tzero < tzero) || (image[i].tzero+image[i].trate*image[i].NY > tzero + trange))) continue;
    152     if (!FindMosaicForImage (image, Nimage, i)) continue;
     158    if (!(Nmosaic = FindMosaicForImage (image, Nimage, i))) continue;
     159    Nmosaic --; // XXX kind of a hack: FindMosaicForImage returns 0 or the mosaic seq number + 1
    153160    if (photcode) {
    154161      if ( photcodeEquiv && (photcode[0].code != GetPhotcodeEquivCodebyCode(image[i].photcode))) continue;
    … …  
    161168
    162169    typehash = wordhash (&image[i].coords.ctype[4]);
     170
     171    if (photcode && SOLO_MOSAIC) {
     172      // mosaic (DIS) images are not currently given a photcode : plot these via the WRP entries
     173      /* DIS images represent a field, not a chip */
     174      if (typehash != ChipImage) continue;
     175      if (foundMosaic[Nmosaic]) continue;
     176      x[0] = -0.5*image[Nmosaic].NX; y[0] = -0.5*image[Nmosaic].NY;
     177      x[1] = +0.5*image[Nmosaic].NX; y[1] = -0.5*image[Nmosaic].NY;
     178      x[2] = +0.5*image[Nmosaic].NX; y[2] = +0.5*image[Nmosaic].NY;
     179      x[3] = -0.5*image[Nmosaic].NX; y[3] = +0.5*image[Nmosaic].NY;
     180      for (j = 0; j < Npts; j++) {
     181        status = XY_to_RD (&r[j], &d[j], x[j], y[j], &image[Nmosaic].coords);
     182        if (!status) break;
     183        r[j] = ohana_normalize_angle (r[j]);
     184        while (r[j] < Rmin) { r[j] += 360.0; }
     185        while (r[j] > Rmax) { r[j] -= 360.0; }
     186        if (j == 0) {
     187          leftside = (r[j] < Rmid);
     188        }
     189        if (j > 0) {
     190          if ( leftside && (r[j] > Rmid + 90)) { r[j] -= 360.0; }
     191          if (!leftside && (r[j] < Rmid - 90)) { r[j] += 360.0; }
     192        }
     193      }
     194      foundMosaic[Nmosaic] = TRUE;
     195      goto plot_points;
     196    }
    163197
    164198    /* DIS images represent a field, not a chip */
    … …  
    270304    if (Npts == 0) continue;
    271305
     306  plot_points:
     307
    272308    status = FALSE;
    273309    for (j = 0; j < Npts; j++) {
    … …  
    320356  free (Yvec.elements.Flt);
    321357  FreeImages (image);
     358  free (foundMosaic);
    322359  return (TRUE);
    323360
Note: See TracChangeset for help on using the changeset viewer.