IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 30443


Ignore:
Timestamp:
Jan 31, 2011, 2:47:28 PM (15 years ago)
Author:
eugene
Message:

speed up mosaic match to detections by limiting the search region and sorting

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

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ipp-20101205/Ohana/src/relphot/Makefile

    r17242 r30443  
    2929$(SRC)/StarOps.$(ARCH).o         \
    3030$(SRC)/args.$(ARCH).o            \
     31$(SRC)/help.$(ARCH).o            \
    3132$(SRC)/bcatalog.$(ARCH).o        \
    3233$(SRC)/global_stats.$(ARCH).o    \

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

    r27586 r30443  
    215215void          write_coords        PROTO((Header *header, Coords *coords));
    216216int relphot_objects (void);
     217
     218void relphot_usage (void);
     219void relphot_help (int argc, char **argv);

  • branches/eam_branches/ipp-20101205/Ohana/src/relphot/src/MosaicOps.c

    r29001 r30443  
    400400      mark = (N < IMAGE_TOOFEW) || (N < IMAGE_GOOD_FRACTION*Nlist[i]);
    401401      if (mark) {
    402         fprintf (stderr, "marked mosaic %s ("OFF_T_FMT"), ("OFF_T_FMT" < %d) || ("OFF_T_FMT" < %f*"OFF_T_FMT")\n", image[imlist[i][0]].name,  i,  N, IMAGE_TOOFEW,  N, IMAGE_GOOD_FRACTION,  Nlist[i]);
     402        if (VERBOSE2) { fprintf (stderr, "marked mosaic %s ("OFF_T_FMT"), ("OFF_T_FMT" < %d) || ("OFF_T_FMT" < %f*"OFF_T_FMT")\n", image[imlist[i][0]].name,  i,  N, IMAGE_TOOFEW,  N, IMAGE_GOOD_FRACTION,  Nlist[i]); }
    403403        mosaic[i].flags |= ID_IMAGE_PHOTOM_FEW;
    404404        Nfew ++;

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

    r27586 r30443  
    11# include "relphot.h"
    2 void usage (void);
    32
    43int args (int argc, char **argv) {
    … …  
    228227  if (UpdateAverages && (argc == 1)) return TRUE;
    229228  if (UserPatchSelect && (argc == 2)) return TRUE;
    230   if (argc != 3) usage ();
     229  if (argc != 3) relphot_usage ();
    231230
    232231  return TRUE;
    233232}
    234 
    235 void usage () {
    236   fprintf (stderr, "ERROR: USAGE: relphot (region) (photcode)\n");
    237   fprintf (stderr, "       or:    relphot (photcode) -region RA RA DEC DEC\n");
    238   fprintf (stderr, "       or:    relphot -averages -region RA RA DEC DEC\n");
    239   fprintf (stderr, "  options: \n");
    240   fprintf (stderr, "  -time (start) (stop)\n");
    241   fprintf (stderr, "  -v : verbose output\n");
    242   fprintf (stderr, "  -vv : more verbose output\n");
    243   fprintf (stderr, "  -outroot (outroot)\n");
    244   fprintf (stderr, "  -plot\n");
    245   fprintf (stderr, "  -plotdelay (seconds)\n");
    246   fprintf (stderr, "  -statmode (mode)\n");
    247   fprintf (stderr, "  -refcode (name) : give extra weight to this photcode\n");
    248   fprintf (stderr, "  -n (nloop)\n");
    249   fprintf (stderr, "  -reset\n");
    250   fprintf (stderr, "  -update\n");
    251   fprintf (stderr, "  -params\n");
    252   fprintf (stderr, "  -mosaic (mosaic)\n");
    253   fprintf (stderr, "  -imfreeze\n");
    254   fprintf (stderr, "  -grid\n");
    255   fprintf (stderr, "  -area Xmin Xmax Ymin Ymax\n");
    256   fprintf (stderr, "  -instmag min max\n");
    257   fprintf (stderr, "  \n");
    258   exit (2);
    259 }

  • branches/eam_branches/ipp-20101205/Ohana/src/relphot/src/bcatalog.c

    r28660 r30443  
    119119
    120120  if (VERBOSE) {
    121     fprintf (stderr, "using "OFF_T_FMT" stars ("OFF_T_FMT" measures) of "OFF_T_FMT" for catalog\n",
    122               subcatalog[0].Naverage,  subcatalog[0].Nmeasure,  i);
     121    fprintf (stderr, "using "OFF_T_FMT" stars ("OFF_T_FMT" measures) of "OFF_T_FMT" for catalog %s\n",
     122             subcatalog[0].Naverage,  subcatalog[0].Nmeasure,  i, catalog[0].filename);
    123123    fprintf (stderr, "rejections: %d code, %d time, %d dophot, %d mag, %d sigma, %d imag, %d few\n",
    124124             Ncode, Ntime, Ndophot, Nmag, Nsigma, Nimag, Nfew);

  • branches/eam_branches/ipp-20101205/Ohana/src/relphot/src/initialize.c

    r29001 r30443  
    55  int N;
    66
     7  relphot_help (argc, argv);
    78  ConfigInit (&argc, argv);
    89  args (argc, argv);

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

    r29754 r30443  
    66  Catalog *catalog, tcatalog;
    77
    8   if (VERBOSE) fprintf (stderr, "loading catalog data\n");
     8  if (VERBOSE2) fprintf (stderr, "loading catalog data\n");
    99
    1010  ALLOCATE (catalog, Catalog, skylist[0].Nregions);
    … …  
    2222    tcatalog.Nsecfilt  = GetPhotcodeNsecfilt ();               // set the desired number in case we need to create the catalog
    2323
    24     if (!dvo_catalog_open (&tcatalog, skylist[0].regions[i], VERBOSE, "r")) {
     24    if (!dvo_catalog_open (&tcatalog, skylist[0].regions[i], VERBOSE2, "r")) {
    2525      fprintf (stderr, "ERROR: failure reading catalog %s\n", tcatalog.filename);
    2626      exit (1);
    2727    }
    28     if (VERBOSE && !tcatalog.Naves_disk) fprintf (stderr, "no data in %s, skipping\n", tcatalog.filename);
     28    if (!tcatalog.Naves_disk) {
     29        if (VERBOSE2) { fprintf (stderr, "no data in %s, skipping\n", tcatalog.filename); }
     30        dvo_catalog_unlock (&tcatalog);
     31        dvo_catalog_free (&tcatalog);
     32        continue;
     33    }
     34
    2935    Nstar_total += tcatalog.Naverage;
    3036    Nmeas_total += tcatalog.Nmeasure;

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

    r29001 r30443  
    11# include "relphot.h"
     2
     3# define MARKTIME(MSG,...) { \
     4  float dtime; \
     5  gettimeofday (&stop, (void *) NULL); \
     6  dtime = DTIME (stop, start); \
     7  fprintf (stderr, MSG, __VA_ARGS__); }
    28
    39SkyList *load_images (FITS_DB *db, char *regionName, SkyRegion *region, int RegionSelect) {
    … …  
    612  off_t      Nimage, Nsubset, Nchar;
    713  off_t     *LineNumber;
     14  struct timeval start, stop;
    815
    916  SkyTable *sky = NULL;
    1017  SkyList *skylist = NULL;
     18
     19  gettimeofday (&start, (void *) NULL);
    1120
    1221  // load the current sky table (layout of all SkyRegions)
    … …  
    2938      exit (2);
    3039  }
     40  MARKTIME("read image table: %f sec\n", dtime);
    3141
    3242  // select the images which overlap the selected sky regions
    3343  subset = select_images (skylist, image, Nimage, &LineNumber, &Nsubset);
     44  MARKTIME("selected images: %f sec\n", dtime);
    3445
    3546  gfits_vtable_from_ftable (&db[0].ftable, &db[0].vtable, LineNumber, Nsubset);
     47  MARKTIME("converted ftable to vtable: %f sec\n", dtime);
    3648
    3749  initImages (subset, Nsubset);
     50  MARKTIME("init images: %f sec\n", dtime);
     51
    3852  initMosaics (subset, Nsubset);
     53  MARKTIME("init mosaics: %f sec\n", dtime);
    3954 
    4055  return (skylist);

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

    r29001 r30443  
    1515void dsortindex (double *X, off_t *Y, int N);
    1616off_t getRegionStartByRA (double R, double *Rref, off_t Nregions);
     17
     18# define MARKTIME(MSG,...) { \
     19  float dtime; \
     20  gettimeofday (&stop, (void *) NULL); \
     21  dtime = DTIME (stop, start); \
     22  fprintf (stderr, MSG, __VA_ARGS__); }
    1723
    1824Image *select_images (SkyList *skylist, Image *timage, off_t Ntimage, off_t **LineNumber, off_t *Nimage) {
    … …  
    2531  Coords tcoords;
    2632  SkyRegionCoords *skycoords;
    27  
     33  struct timeval start, stop;
     34 
     35  double RmaxSkyRegion, RminSkyRegion, DminSkyRegion, DmaxSkyRegion;
     36
    2837  double *RmaxSky;
    2938  off_t *index;
    … …  
    3544    return NULL;
    3645  }
     46
     47  gettimeofday (&start, (void *) NULL);
    3748
    3849  // the comparison is made in the catalog local projection. below we set crval1,2
    … …  
    4758  ALLOCATE (RmaxSky, double, skylist[0].Nregions);
    4859  ALLOCATE (index, off_t, skylist[0].Nregions);
     60
     61  RminSkyRegion = +360.0;
     62  RmaxSkyRegion = -360.0;
     63  DminSkyRegion = +90.0;
     64  DmaxSkyRegion = -90.0;
    4965
    5066  /* compare with each region file */
    … …  
    7591    skycoords[i].Xc[3] -= dx; skycoords[i].Yc[3] += dy;
    7692    skycoords[i].Xc[4] -= dx; skycoords[i].Yc[4] -= dy;
    77   }
     93
     94    RminSkyRegion = MIN(RminSkyRegion, skylist[0].regions[i][0].Rmin);
     95    RmaxSkyRegion = MAX(RmaxSkyRegion, skylist[0].regions[i][0].Rmax);
     96    DminSkyRegion = MIN(DminSkyRegion, skylist[0].regions[i][0].Dmin);
     97    DmaxSkyRegion = MAX(DmaxSkyRegion, skylist[0].regions[i][0].Dmax);
     98  }
     99  MARKTIME("create sky region coords: %f sec\n", dtime);
    78100
    79101  dsortindex (RmaxSky, index, skylist[0].Nregions);
     102  MARKTIME("sort sky coords: %f sec\n", dtime);
    80103
    81104  if (VERBOSE) fprintf (stderr, "finding images\n");
    82105  BuildChipMatch (timage, Ntimage);
     106  MARKTIME("build chip match: %f sec\n", dtime);
    83107
    84108  nimage = 0;
    … …  
    100124    }
    101125   
     126    // this adds 1.3 sec for 3M images
    102127    if (!FindMosaicForImage (timage, Ntimage, i)) {
    103128      fprintf (stderr, "cannot find mosaic for "OFF_T_FMT"\n", i);
    … …  
    113138    found = FALSE;
    114139
    115     /* transform corners to ra,dec */
     140    /* transform corners to ra,dec -- costs ~3sec for 3M images */
    116141    double RminImage = 360.0;
     142    double RmaxImage =   0.0;
     143    double DminImage = +90.0;
     144    double DmaxImage = -90.0;
    117145    for (j = 0; j < 5; j++) {
    118146      XY_to_RD (&Ri[j], &Di[j], Xi[j], Yi[j], &timage[i].coords);
    119147      RminImage = MIN(RminImage, Ri[j]);
    120     }
    121 
    122     // RA(nStart) is guaranteed to be < RminImage:
     148      RmaxImage = MAX(RmaxImage, Ri[j]);
     149      DminImage = MIN(DminImage, Di[j]);
     150      DmaxImage = MAX(DmaxImage, Di[j]);
     151    }
     152    if (RmaxImage - RminImage > 180.0) {
     153        double tmp = RminImage;
     154        RmaxImage = RminImage;
     155        RminImage = tmp - 360.0;
     156    }
     157   
     158    // check that this image is even in range of the searched region
     159    if (DminImage > DmaxSkyRegion) continue;
     160    if (DmaxImage < DminSkyRegion) continue;
     161   
     162    // the sky region RA is defined to be 0 - 360.0
     163    if (RminImage > RmaxSkyRegion) continue;
     164    if (RmaxImage < RminSkyRegion) continue;
     165
     166    // RA(nStart) is guaranteed to be < RminImage: -- costs 0.5sec for 3M images
    123167    nStart = getRegionStartByRA (RminImage, RmaxSky, skylist[0].Nregions);
    124168
    125169    /* compare with each region file */
    126     for (iSky = 0; (iSky < skylist[0].Nregions) && !found; iSky++) {
     170    for (iSky = nStart; (iSky < skylist[0].Nregions) && !found; iSky++) {
    127171
    128172      m = index[iSky];
    … …  
    174218    }
    175219  }
    176      
     220  MARKTIME("finish image selection: %f sec\n", dtime);
     221
    177222  if (VERBOSE) fprintf (stderr, "found "OFF_T_FMT" images\n", nimage);
    178223
    … …  
    295340  return (Nlo);
    296341}
     342
     343off_t getRegionStopByRA (double R, double *Rref, off_t Nregions) {
     344
     345  // use bisection to find the overlapping mosaic
     346
     347  off_t Nlo, Nhi, N;
     348
     349  // find the last mosaic before start
     350  Nlo = 0; Nhi = Nregions;
     351  while (Nhi - Nlo > 10) {
     352    N = 0.5*(Nlo + Nhi);
     353    if (Rref[N] < R) {
     354      Nlo = MAX(N, 0);
     355    } else {
     356      Nhi = MIN(N, Nregions);
     357    }
     358  }
     359  return (Nlo);
     360}
Note: See TracChangeset for help on using the changeset viewer.