IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
May 16, 2012, 4:30:08 PM (14 years ago)
Author:
eugene
Message:

fix selection of images, especially at 0,360 boundary; catch failures to write image table; add ApplyOffsets mode [-apply-offsets]

Location:
branches/eam_branches/ipp-20120405/Ohana/src/relphot
Files:
8 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ipp-20120405/Ohana/src/relphot/include/relphot.h

    r33874 r33885  
    145145int    PLOTDELAY;
    146146int    UpdateAverages;
     147int    ApplyOffsets;
    147148
    148149char  *PhotcodeList;
    … …  
    250251
    251252SkyList      *load_images         PROTO((FITS_DB *db, char *regionName, SkyRegion *region));
    252 Image        *select_images       PROTO((SkyList *skylist, Image *timage, off_t Ntimage, char *inSubset, off_t **LineNumber, off_t *Nimage));
     253Image        *select_images       PROTO((SkyList *skylist, Image *timage, off_t Ntimage, char *inSubset, off_t **LineNumber, off_t *Nimage, SkyRegion *region));
    253254
    254255int           main                PROTO((int argc, char **argv));

  • branches/eam_branches/ipp-20120405/Ohana/src/relphot/src/ImageSubset.c

    r33651 r33885  
    9292}
    9393
     94// STATUS is value expected for success
     95# define CHECK_STATUS(STATUS,MSG,...)                                   \
     96  if (!(STATUS)) {                                                      \
     97    fprintf (stderr, MSG, __VA_ARGS__);                                 \
     98    return FALSE;                                                       \
     99  }
     100
    94101int ImageSubsetSave(char *filename, ImageSubset *image, off_t Nimage) {
    95102
    … …  
    167174  }
    168175
    169   gfits_fwrite_header  (f, &header);
    170   gfits_fwrite_matrix  (f, &matrix);
    171   gfits_fwrite_Theader (f, &theader);
    172   gfits_fwrite_table  (f, &ftable);
     176  int status;
     177  status = gfits_fwrite_header  (f, &header);
     178  CHECK_STATUS (status, "ERROR: cannot write header for image subset %s\n", filename);
     179
     180  status = gfits_fwrite_matrix  (f, &matrix);
     181  CHECK_STATUS (status, "ERROR: cannot write matrix for image subset %s\n", filename);
     182
     183  status = gfits_fwrite_Theader (f, &theader);
     184  CHECK_STATUS (status, "ERROR: cannot write table header for image subset %s\n", filename);
     185
     186  status = gfits_fwrite_table  (f, &ftable);
     187  CHECK_STATUS (status, "ERROR: cannot write table data for image subset %s\n", filename);
    173188
    174189  gfits_free_header (&header);
    … …  
    177192  gfits_free_table (&ftable);
    178193
    179   fclose (f);
    180   fflush (f);
     194  int fd = fileno (f);
     195
     196  status = fflush (f);
     197  CHECK_STATUS (!status, "ERROR: cannot flush file image subset %s\n", filename);
     198
     199  status = fsync (fd);
     200  CHECK_STATUS (!status, "ERROR: cannot flush file image subset %s\n", filename);
     201
     202  status = fclose (f);
     203  CHECK_STATUS (!status, "ERROR: problem closing image subset file %s\n", filename);
    181204
    182205  return TRUE;

  • branches/eam_branches/ipp-20120405/Ohana/src/relphot/src/args.c

    r33797 r33885  
    307307  }
    308308
     309  ApplyOffsets = FALSE;
     310  if ((N = get_argument (argc, argv, "-apply-offsets"))) {
     311    remove_argument (N, &argc, argv);
     312    ApplyOffsets = TRUE;
     313  }
     314
    309315  if (UpdateAverages && (argc == 1)) return TRUE;
    310316  if (argc != 2) relphot_usage ();

  • branches/eam_branches/ipp-20120405/Ohana/src/relphot/src/load_catalogs.c

    r33874 r33885  
    190190  }
    191191  if (!PARALLEL_MANUAL && !PARALLEL_SERIAL) {
    192     HostTableWaitJobsGetIO (table, __FILE__, __LINE__, VERBOSE);
     192    // if any one of the remote jobs fails, we should fail and exit
     193    int status = HostTableWaitJobsGetIO (table, __FILE__, __LINE__, VERBOSE);
     194    if (!status) {
     195      fprintf (stderr, "at least one remote client job failed to load data, exiting\n");
     196      exit (3);
     197    }
    193198  }
    194199

  • branches/eam_branches/ipp-20120405/Ohana/src/relphot/src/load_images.c

    r33797 r33885  
    5050    // select the images which overlap the selected sky regions
    5151    // 'subset' points to a new copy of the data (different from 'image')
    52     subset = select_images (skylist, image, Nimage, inSubset, &LineNumber, &Nsubset);
     52    subset = select_images (skylist, image, Nimage, inSubset, &LineNumber, &Nsubset, NULL);
    5353    MARKTIME("selected %d overlapping images: %f sec\n", (int) Nsubset, dtime);
    5454  } else {
    … …  
    5959      skylist = SkyListByPatch (sky, -1, &subregion);
    6060
    61       subregion.Rmin = region[0].Rmin; subregion.Rmax = 360.0;
    62       SkyList *extraList = SkyListByPatch (sky, -1, &subregion);
    63 
    6461      // select the images which overlap the selected sky regions
    6562      // 'subset' points to a new copy of the data (different from 'image')
    66       subset = select_images (skylist, image, Nimage, inSubset, &LineNumber, &Nsubset);
     63      subset = select_images (skylist, image, Nimage, inSubset, &LineNumber, &Nsubset, &subregion);
    6764      MARKTIME("selected %d overlapping images: %f sec\n", (int) Nsubset, dtime);
    6865
    … …  
    7168      off_t *LineNumberExtra, i;
    7269
    73       subsetExtra = select_images (extraList, image, Nimage, inSubset, &LineNumberExtra, &NsubsetExtra);
     70      subregion.Rmin = region[0].Rmin; subregion.Rmax = 360.0;
     71      SkyList *extraList = SkyListByPatch (sky, -1, &subregion);
     72
     73      subsetExtra = select_images (extraList, image, Nimage, inSubset, &LineNumberExtra, &NsubsetExtra, &subregion);
    7474      MARKTIME("selected %d overlapping images: %f sec\n", (int) Nsubset, dtime);
    7575
    … …  
    8888      // select the images which overlap the selected sky regions
    8989      // 'subset' points to a new copy of the data (different from 'image')
    90       subset = select_images (skylist, image, Nimage, inSubset, &LineNumber, &Nsubset);
     90      subset = select_images (skylist, image, Nimage, inSubset, &LineNumber, &Nsubset, region);
    9191      MARKTIME("selected %d overlapping images: %f sec\n", (int) Nsubset, dtime);
    9292    }
    9393  }
     94
     95# if (1)
     96  FILE *ftest = fopen ("skydump.dat", "w");
     97  int i;
     98  for (i = 0; i < skylist[0].Nregions; i++) {
     99    fprintf (ftest, "%s : %f %f : %f %f\n", skylist[0].regions[i][0].name, skylist[0].regions[i][0].Rmin, skylist[0].regions[i][0].Rmax, skylist[0].regions[i][0].Dmin, skylist[0].regions[i][0].Dmax);
     100  }
     101  fclose (ftest);
     102# endif
    94103
    95104  // match chips to mosaics (if applicable)

  • branches/eam_branches/ipp-20120405/Ohana/src/relphot/src/reload_catalogs.c

    r33874 r33885  
    120120int reload_catalogs_parallel (SkyList *sky) {
    121121
    122   off_t Nimage;
    123   ImageSubset *image = getimages_subset (&Nimage);
    124 
    125   // write out the subset table of image information
     122  // name of image subset file:
    126123  char imageFile[512];
    127124  snprintf (imageFile, 512, "%s/Images.subset.dat", CATDIR);
    128   if (!ImageSubsetSave (imageFile, image, Nimage)) {
    129     fprintf (stderr, "failed to write image subset\n");
    130     exit (1);
    131   }
    132   free (image);
     125
     126  // the ApplyOffsets option re-uses an existing Images.subset.dat file
     127  if (!ApplyOffsets) {
     128    off_t Nimage;
     129    ImageSubset *image = getimages_subset (&Nimage);
     130    if (!ImageSubsetSave (imageFile, image, Nimage)) {
     131      fprintf (stderr, "failed to write image subset\n");
     132      exit (1);
     133    }
     134    free (image);
     135  } else {
     136    // load the current sky table (layout of all SkyRegions)
     137    SkyTable *skytable = SkyTableLoadOptimal (CATDIR, SKY_TABLE, GSCFILE, TRUE, SKY_DEPTH, VERBOSE);
     138    // XXX kind of hackish...
     139    ALLOCATE (sky, SkyList, 1);
     140    strcpy (sky->hosts, skytable->hosts);
     141  }
    133142
    134143  // now launch the relphot_client jobs to the parallel hosts
    … …  
    202211  }
    203212  if (!PARALLEL_MANUAL && !PARALLEL_SERIAL) {
    204     HostTableWaitJobsGetIO (table, __FILE__, __LINE__, VERBOSE);
     213    int status = HostTableWaitJobsGetIO (table, __FILE__, __LINE__, VERBOSE);
     214    if (!status) {
     215      fprintf (stderr, "at least one remote client job failed to load data, exiting\n");
     216      exit (3);
     217    }
    205218  }
    206219  return (TRUE);

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

    r33827 r33885  
    2323  if (UpdateAverages) {
    2424    relphot_objects (0, NULL);
     25    exit (0);
     26  }
     27
     28  if (ApplyOffsets) {
     29    // re-run this step from a previous attempt (assumes an existing Images.subset.dat file)
     30    if (!PARALLEL) {
     31      fprintf (stderr, "-apply-offsets only makes sense in an parallel context\n");
     32      exit (2);
     33    }
     34    reload_catalogs (NULL, NULL, 0, NULL);
    2535    exit (0);
    2636  }

  • branches/eam_branches/ipp-20120405/Ohana/src/relphot/src/select_images.c

    r33827 r33885  
    1313} SkyRegionCoords;
    1414
     15void MakeSkyCoordIndex (SkyRegionCoords **skycoords_out, double **RmaxSky_out, off_t **index_out, SkyList *skylist);
    1516void dsortindex (double *X, off_t *Y, int N);
    1617off_t getRegionStartByRA (double R, double *Rref, off_t Nregions);
    … …  
    2223  fprintf (stderr, MSG, __VA_ARGS__); }
    2324
    24 Image *select_images (SkyList *skylist, Image *timage, off_t Ntimage, char *inSubset, off_t **LineNumber, off_t *Nimage) {
     25Image *select_images (SkyList *skylist, Image *timage, off_t Ntimage, char *inSubset, off_t **LineNumber, off_t *Nimage, SkyRegion *region) {
    2526 
    2627  Image *image;
    … …  
    2829  off_t *line_number;
    2930  int InRange, ecode, found;
    30   double Ri[5], Di[5], Xi[5], Yi[5], dx, dy;
     31  double Ri[5], Di[5], Xi[5], Yi[5];
    3132  Coords tcoords;
    32   SkyRegionCoords *skycoords;
    3333  struct timeval start, stop;
    3434 
    35   double *RmaxSky;
    36   off_t *index;
    37 
    3835  if (skylist[0].Nregions < 1) {
    3936    *Nimage = 0;
    … …  
    5249  strcpy (tcoords.ctype, "RA---TAN");
    5350
    54   ALLOCATE (skycoords, SkyRegionCoords, skylist[0].Nregions);
    55 
    56   ALLOCATE (RmaxSky, double, skylist[0].Nregions);
    57   ALLOCATE (index, off_t, skylist[0].Nregions);
    58 
    59   double RminSkyRegion = +360.0;
    60   double RmaxSkyRegion = -360.0;
    61   double DminSkyRegion = +90.0;
    62   double DmaxSkyRegion = -90.0;
    63 
    64   D_NIMAGE = 6000;
    65 
    66   // FILE *ftest = fopen ("relphot.dump.dat", "w");
    67 
    68   // The regions have RA values in the range 0-360, which implies an RA midpoint of 180.
    69   // A patch which covers the 0,360 boundary (but not much more) needs to treat the values
    70   // as being centered on 0.0.  Run the following test with both and see which yields a
    71   // smaller contiguious patch
    72 
    73   double RminSkyRegionAlt = +360.0;
    74   double RmaxSkyRegionAlt = -360.0;
    75 
    76   /* compare with each region file */
    77   for (i = 0; i < skylist[0].Nregions; i++) {
    78 
    79     /* we make positional comparisons in the projection of catalog */
    80     skycoords[i].Rc = 0.5*(skylist[0].regions[i][0].Rmax + skylist[0].regions[i][0].Rmin);
    81     skycoords[i].Dc = 0.5*(skylist[0].regions[i][0].Dmax + skylist[0].regions[i][0].Dmin);
    82     tcoords.crval1 = skycoords[i].Rc;
    83     tcoords.crval2 = skycoords[i].Dc;
    84 
    85     /* define catalog corners */
    86     RD_to_XY (&skycoords[i].Xc[0], &skycoords[i].Yc[0], skylist[0].regions[i][0].Rmin, skylist[0].regions[i][0].Dmin, &tcoords);
    87     RD_to_XY (&skycoords[i].Xc[1], &skycoords[i].Yc[1], skylist[0].regions[i][0].Rmax, skylist[0].regions[i][0].Dmin, &tcoords);
    88     RD_to_XY (&skycoords[i].Xc[2], &skycoords[i].Yc[2], skylist[0].regions[i][0].Rmax, skylist[0].regions[i][0].Dmax, &tcoords);
    89     RD_to_XY (&skycoords[i].Xc[3], &skycoords[i].Yc[3], skylist[0].regions[i][0].Rmin, skylist[0].regions[i][0].Dmax, &tcoords);
    90     skycoords[i].Xc[4] = skycoords[i].Xc[0];   
    91     skycoords[i].Yc[4] = skycoords[i].Yc[0];   
    92 
    93     RmaxSky[i] = skylist[0].regions[i][0].Rmax;
    94     index[i] = i;
    95 
    96     dx = 0.02*(skycoords[i].Xc[2] - skycoords[i].Xc[0]);
    97     dy = 0.02*(skycoords[i].Yc[2] - skycoords[i].Yc[0]);
    98     skycoords[i].Xc[0] -= dx; skycoords[i].Yc[0] -= dy;
    99     skycoords[i].Xc[1] += dx; skycoords[i].Yc[1] -= dy;
    100     skycoords[i].Xc[2] += dx; skycoords[i].Yc[2] += dy;
    101     skycoords[i].Xc[3] -= dx; skycoords[i].Yc[3] += dy;
    102     skycoords[i].Xc[4] -= dx; skycoords[i].Yc[4] -= dy;
    103 
    104     RminSkyRegion = MIN(RminSkyRegion, skylist[0].regions[i][0].Rmin);
    105     RmaxSkyRegion = MAX(RmaxSkyRegion, skylist[0].regions[i][0].Rmax);
    106     DminSkyRegion = MIN(DminSkyRegion, skylist[0].regions[i][0].Dmin);
    107     DmaxSkyRegion = MAX(DmaxSkyRegion, skylist[0].regions[i][0].Dmax);
    108 
    109     double RminAlt   = ohana_normalize_angle_to_midpoint (skylist[0].regions[i][0].Rmin, 0.0);
    110     double RmaxAlt   = ohana_normalize_angle_to_midpoint (skylist[0].regions[i][0].Rmax, 0.0);
    111     RminSkyRegionAlt = MIN(RminSkyRegionAlt, RminAlt);
    112     RmaxSkyRegionAlt = MAX(RmaxSkyRegionAlt, RmaxAlt);
    113 
    114     RminSkyRegion = MIN(RminSkyRegion, skylist[0].regions[i][0].Rmin);
    115     RmaxSkyRegion = MAX(RmaxSkyRegion, skylist[0].regions[i][0].Rmax);
    116     DminSkyRegion = MIN(DminSkyRegion, skylist[0].regions[i][0].Dmin);
    117     DmaxSkyRegion = MAX(DmaxSkyRegion, skylist[0].regions[i][0].Dmax);
    118   }
    119   // use the option which gives the smaller RA range
    120   if ((RmaxSkyRegionAlt - RminSkyRegionAlt) < (RmaxSkyRegion - RminSkyRegion)) {
    121     RminSkyRegion = RminSkyRegionAlt;
    122     RmaxSkyRegion = RmaxSkyRegionAlt;
    123   }
     51  double RminSkyRegion = region[0].Rmin;
     52  double RmaxSkyRegion = region[0].Rmax;
     53  double DminSkyRegion = region[0].Dmin;
     54  double DmaxSkyRegion = region[0].Dmax;
     55
    12456  double RmidSkyRegion = 0.5*(RminSkyRegion + RmaxSkyRegion);
    125   MARKTIME("create sky region coords: %f sec\n", dtime);
    126 
    127   dsortindex (RmaxSky, index, skylist[0].Nregions);
    128   MARKTIME("sort sky coords: %f sec\n", dtime);
    12957
    13058  // generate a set of Rmin/Rmax values using this RmidSkyRegion for DEC bands of width
    … …  
    15482  }
    15583
    156 # if (0)
    157   // XXX quick test
     84# if (1)
     85  // XXX quick test of the dec bands:
    15886  FILE *fout = fopen ("dec.bands.dat", "w");
    15987  for (i = 0; i < NDecBands; i++) {
    … …  
    16391# endif
    16492
     93  SkyRegionCoords *skycoords = NULL;
     94  double *RmaxSky = NULL;
     95  off_t *index = NULL;
     96
     97  if (!USE_BASIC_CHECK) {
     98    MakeSkyCoordIndex (&skycoords, &RmaxSky, &index, skylist);
     99  }
     100
    165101  if (VERBOSE) fprintf (stderr, "finding images\n");
    166102  BuildChipMatch (timage, Ntimage);
    167103  MARKTIME("build chip match for %d images: %f sec\n", (int) Ntimage, dtime);
     104
     105  D_NIMAGE = 6000;
    168106
    169107  nimage = 0;
    … …  
    266204      continue;
    267205    }
    268 
    269     // fprintf (ftest, "%s %lf %lf - %lf %lf : %lf %lf\n", timage[i].name, RminImage, DminImage, RmaxImage, DmaxImage, 0.5*(RminImage + RmaxImage), 0.5*(DminImage + DmaxImage));
    270206
    271207    // image overlaps region, keep it
    … …  
    339275  MARKTIME("finish image selection: %f sec\n", dtime);
    340276
    341   // fclose (ftest);
    342 
    343277  if (VERBOSE) fprintf (stderr, "found "OFF_T_FMT" images\n", nimage);
    344278
    345279  REALLOCATE (image, Image, MAX (nimage, 1));
    346280  REALLOCATE (line_number, off_t, MAX (nimage, 1));
    347   free (skycoords);
    348   free (RmaxSky);
    349   free (index);
     281
     282  if (!USE_BASIC_CHECK) {
     283    free (skycoords);
     284    free (RmaxSky);
     285    free (index);
     286  }
    350287
    351288  *Nimage  = nimage;
    … …  
    463400  return (Nlo);
    464401}
     402
     403void MakeSkyCoordIndex (SkyRegionCoords **skycoords_out, double **RmaxSky_out, off_t **index_out, SkyList *skylist) {
     404
     405  off_t i;
     406  SkyRegionCoords *skycoords;
     407  Coords tcoords;
     408  double *RmaxSky;
     409  off_t *index;
     410  double dx, dy;
     411
     412  // the comparison is made in the catalog local projection. below we set crval1,2
     413  tcoords.crpix1 = tcoords.crpix2 = 0.0;
     414  tcoords.cdelt1 = tcoords.cdelt2 = 1.0 / 3600.0;
     415  tcoords.pc1_1  = tcoords.pc2_2 = 1.0;
     416  tcoords.pc1_2  = tcoords.pc2_1 = 0.0;
     417  strcpy (tcoords.ctype, "RA---TAN");
     418
     419  /* compare with each region file */
     420  ALLOCATE (skycoords, SkyRegionCoords, skylist[0].Nregions);
     421
     422  ALLOCATE (RmaxSky, double, skylist[0].Nregions);
     423  ALLOCATE (index, off_t, skylist[0].Nregions);
     424
     425  for (i = 0; i < skylist[0].Nregions; i++) {
     426
     427    /* we make positional comparisons in the projection of catalog */
     428    skycoords[i].Rc = 0.5*(skylist[0].regions[i][0].Rmax + skylist[0].regions[i][0].Rmin);
     429    skycoords[i].Dc = 0.5*(skylist[0].regions[i][0].Dmax + skylist[0].regions[i][0].Dmin);
     430    tcoords.crval1 = skycoords[i].Rc;
     431    tcoords.crval2 = skycoords[i].Dc;
     432
     433    /* define catalog corners */
     434    RD_to_XY (&skycoords[i].Xc[0], &skycoords[i].Yc[0], skylist[0].regions[i][0].Rmin, skylist[0].regions[i][0].Dmin, &tcoords);
     435    RD_to_XY (&skycoords[i].Xc[1], &skycoords[i].Yc[1], skylist[0].regions[i][0].Rmax, skylist[0].regions[i][0].Dmin, &tcoords);
     436    RD_to_XY (&skycoords[i].Xc[2], &skycoords[i].Yc[2], skylist[0].regions[i][0].Rmax, skylist[0].regions[i][0].Dmax, &tcoords);
     437    RD_to_XY (&skycoords[i].Xc[3], &skycoords[i].Yc[3], skylist[0].regions[i][0].Rmin, skylist[0].regions[i][0].Dmax, &tcoords);
     438    skycoords[i].Xc[4] = skycoords[i].Xc[0];   
     439    skycoords[i].Yc[4] = skycoords[i].Yc[0];   
     440
     441    RmaxSky[i] = skylist[0].regions[i][0].Rmax;
     442    index[i] = i;
     443
     444    dx = 0.02*(skycoords[i].Xc[2] - skycoords[i].Xc[0]);
     445    dy = 0.02*(skycoords[i].Yc[2] - skycoords[i].Yc[0]);
     446    skycoords[i].Xc[0] -= dx; skycoords[i].Yc[0] -= dy;
     447    skycoords[i].Xc[1] += dx; skycoords[i].Yc[1] -= dy;
     448    skycoords[i].Xc[2] += dx; skycoords[i].Yc[2] += dy;
     449    skycoords[i].Xc[3] -= dx; skycoords[i].Yc[3] += dy;
     450    skycoords[i].Xc[4] -= dx; skycoords[i].Yc[4] -= dy;
     451  }
     452
     453  dsortindex (RmaxSky, index, skylist[0].Nregions);
     454
     455  *RmaxSky_out = RmaxSky;
     456  *index_out = index;
     457  *skycoords_out = skycoords;
     458
     459  return;
     460}
Note: See TracChangeset for help on using the changeset viewer.