IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Mar 23, 2010, 5:26:17 PM (16 years ago)
Author:
eugene
Message:

add timing marks; use sorted indexes and bracketting to speed the detection->image and image->mosaic matches

Location:
branches/eam_branches/largefiles.20100314/Ohana/src/relastro
Files:
6 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/largefiles.20100314/Ohana/src/relastro/doc/notes.txt

    r17043 r27425  
     1
     22010.03.23
     3
     4  I've been working on optimization.  I've updated ImageOps to make
     5  the measure->image relationship generation use the image IDs along
     6  with a bracket search -- this goes much faster as a result.  More
     7  work is needed:
     8
     9  * select_images needs to apply astrometry to all image corners:
     10    -> cache the image center and radius and use this to narrow down
     11       the searches.
     12
     13  * MosaicOps uses the time to match mosaics.  I've added sorting to
     14    speed this up, but this should be done based on the image and
     15    parent index.
     16
     17    -> make a tool to create the parent ID for existing DBs.
     18    -> make sure addstar is populating the parent IDs
     19    -> use the parent IDs to make the link.
     20
     21load image data: 0.006347 sec
     22  setup sky: 0.091475 sec
     23  convert image table: 0.091498 sec
     24  select images: 2.925978 sec
     25  init images: 2.926957 sec
     26  init mosaics: 2.927319 sec
     27load images: 2.933762 sec
     28load catalog data: 9.772259 sec
     29make image bins: 9.780618 sec
     30set up image indexes: 10.196519 sec
     31
     32load image data: 0.006746 sec
     33  setup sky: 0.091440 sec
     34  convert image table: 0.091467 sec
     35  select images: 0.374193 sec
     36  init images: 0.375355 sec
     37  init mosaics: 0.375753 sec
     38load images: 0.382606 sec
     39load catalog data: 7.260094 sec
     40make image bins: 7.268564 sec
     41set up image indexes: 7.687485 sec
     42
    143
    2442008.03.01

  • branches/eam_branches/largefiles.20100314/Ohana/src/relastro/src/ImageOps.c

    r27312 r27425  
    11# include "relastro.h"
    2 #define TESTING
    3 
    4 static unsigned int *start;
    5 static unsigned int *stop;
     2# define TESTING
     3# define USE_IMAGE_ID 1
    64
    75static off_t       **bin;     // link from catalog,measure to image
    … …  
    1412static off_t        Nimage;   // number of available images
    1513
     14# if USE_IMAGE_ID
     15// we sort (imageIDs, imageIdx) by imageIDs to get a sorted index
     16static off_t        *imageIDs; // list of all image IDs
     17static off_t        *imageIdx; // list of index for image IDs
     18# else
     19static unsigned int *start;
     20static unsigned int *stop;
     21# endif
     22
    1623Image *getimages (off_t *N) {
    1724  *N = Nimage;
    … …  
    3037  Nimage = N;
    3138
    32   ALLOCATE (start,   unsigned, Nimage);
    33   ALLOCATE (stop,    unsigned, Nimage);
     39# if USE_IMAGE_ID
     40  ALLOCATE (imageIDs, off_t, Nimage);
     41  ALLOCATE (imageIdx, off_t, Nimage);
     42
     43  for (i = 0; i < Nimage; i++) {
     44    imageIdx[i] = i;
     45    imageIDs[i] = image[i].imageID;
     46  }
     47  llsortpair (imageIDs, imageIdx, Nimage);
     48# else
     49  ALLOCATE (start, unsigned, Nimage);
     50  ALLOCATE (stop, unsigned, Nimage);
    3451
    3552  for (i = 0; i < Nimage; i++) {
    … …  
    3754    stop[i]  = image[i].tzero + MAX(1.01*image[i].trate*image[i].NY, 1);
    3855  }
     56# endif
     57
     58}
     59
     60off_t getImageByID (off_t ID) {
     61
     62  // we have a pair of vectors (imageIDs, imageIdx) sorted by imageIDs
     63  // use bisection to find the specified image ID
     64
     65# if USE_IMAGE_ID
     66  off_t Nlo, Nhi, N;
     67
     68  Nlo = 0; Nhi = Nimage;
     69  while (Nhi - Nlo > 10) {
     70    N = 0.5*(Nlo + Nhi);
     71    if (imageIDs[N] < ID) {
     72      Nlo = MAX(N, 0);
     73    } else {
     74      Nhi = MIN(N + 1, Nimage);
     75    }
     76  }
     77
     78  for (N = Nlo; N < Nhi; N++) {
     79    if (imageIDs[N] == ID)
     80      return (imageIdx[N]);
     81  }
     82# endif
     83
     84  return (-1);
    3985}
    4086
    … …  
    97143}
    98144
     145# if USE_IMAGE_ID
    99146/* modify this function to use the measure->imageID field */
    100147void matchImage (Catalog *catalog, off_t meas, int cat) {
    101148
     149  off_t idx, ID;
     150  Measure *measure;
     151
     152  measure = &catalog[cat].measure[meas];
     153
     154  ID = measure[0].imageID;
     155  idx = getImageByID (ID);
     156  if (idx == -1) {
     157    fprintf (stderr, "can't match detection to image?\n");
     158    abort();
     159  }
     160
     161  // index for (catalog, measure) -> image
     162  bin[cat][meas] = idx;
     163
     164  // index for image, Nentry -> catalog
     165  clist[idx][Nlist[idx]] = cat;
     166
     167  // index for image, Nentry -> measure
     168  mlist[idx][Nlist[idx]] = meas;
     169  Nlist[idx] ++;
     170
     171  if (Nlist[idx] == NLIST[idx]) {
     172    NLIST[idx] += 100;
     173    REALLOCATE (clist[idx], int,   NLIST[idx]);
     174    REALLOCATE (mlist[idx], off_t, NLIST[idx]);
     175  }
     176  return;
     177}
     178
     179# else
     180/* modify this function to use the measure->imageID field */
     181void matchImage (Catalog *catalog, off_t meas, int cat) {
     182
    102183  off_t i;
    103184  Measure *measure;
    … …  
    107188  /* find the image that supplied this measurement */
    108189  for (i = 0; i < Nimage; i++) {
     190    // let's try the very slow method first before adding a bisection search
     191    // if (image[i].imageID != measure[0].imageID) continue;
    109192    if (image[0].photcode == -1) continue;
    110193    if (measure[0].photcode != image[i].photcode) continue;
    … …  
    129212    return;
    130213  }
    131 }
     214  fprintf (stderr, "can't match detection to image?\n");
     215  abort();
     216}
     217# endif
    132218
    133219Coords *getCoords (off_t meas, int cat) {
    … …  
    199285  mosaic = getMosaicForImage (im);
    200286  if (mosaic != NULL) {
    201       moscoords = &mosaic[0].coords;
     287    moscoords = &mosaic[0].coords;
    202288  }
    203289  imcoords = &image[im].coords;
    … …  
    270356  moscoords = NULL;
    271357  if (mode == MODE_MOSAIC) {
    272       mosaic = getMosaicForImage (im);
    273       if (mosaic == NULL) {
    274         fprintf (stderr, "mosaic not found for image %s\n", image[im].name);
    275         exit (1);
    276       }
    277       moscoords = &mosaic[0].coords;
     358    mosaic = getMosaicForImage (im);
     359    if (mosaic == NULL) {
     360      fprintf (stderr, "mosaic not found for image %s\n", image[im].name);
     361      exit (1);
     362    }
     363    moscoords = &mosaic[0].coords;
    278364  }
    279365
    … …  
    319405        LM_to_RD (&raw[i].R, &raw[i].D, raw[i].P, raw[i].Q, moscoords);
    320406        break;
    321     default:
    322       fprintf (stderr, "error: invalid mode in getImageRaw");
    323       abort ();
     407      default:
     408        fprintf (stderr, "error: invalid mode in getImageRaw");
     409        abort ();
    324410    }
    325411  }
    … …  
    371457    switch (mode) {
    372458      case MODE_SIMPLE:
    373       RD_to_LM (&ref[i].P, &ref[i].Q, ref[i].R, ref[i].D, &image[im].coords);
    374       ref[i].L = ref[i].P;
    375       ref[i].M = ref[i].Q;
    376       LM_to_XY (&ref[i].X, &ref[i].Y, ref[i].L, ref[i].M, &image[im].coords);
    377       break;
     459        RD_to_LM (&ref[i].P, &ref[i].Q, ref[i].R, ref[i].D, &image[im].coords);
     460        ref[i].L = ref[i].P;
     461        ref[i].M = ref[i].Q;
     462        LM_to_XY (&ref[i].X, &ref[i].Y, ref[i].L, ref[i].M, &image[im].coords);
     463        break;
    378464      case MODE_MOSAIC:
    379       RD_to_LM (&ref[i].P, &ref[i].Q, ref[i].R, ref[i].D, moscoords);
    380       LM_to_XY (&ref[i].L, &ref[i].M, ref[i].P, ref[i].Q, moscoords);
    381       LM_to_XY (&ref[i].X, &ref[i].Y, ref[i].L, ref[i].M, &image[im].coords);
    382       break;
     465        RD_to_LM (&ref[i].P, &ref[i].Q, ref[i].R, ref[i].D, moscoords);
     466        LM_to_XY (&ref[i].L, &ref[i].M, ref[i].P, ref[i].Q, moscoords);
     467        LM_to_XY (&ref[i].X, &ref[i].Y, ref[i].L, ref[i].M, &image[im].coords);
     468        break;
    383469      default:
    384470        fprintf (stderr, "invalid case");

  • branches/eam_branches/largefiles.20100314/Ohana/src/relastro/src/MosaicOps.c

    r27312 r27425  
    1919}
    2020
     21off_t getMosaicByTimes (unsigned int start, unsigned int stop, unsigned int *startMos, unsigned int *stopMos, off_t *indexMos) {
     22
     23  // use bisection to find the overlapping mosaic
     24
     25  off_t Nlo, Nhi, N;
     26
     27  // find the last mosaic before start
     28  Nlo = 0; Nhi = Nmosaic;
     29  while (Nhi - Nlo > 10) {
     30    N = 0.5*(Nlo + Nhi);
     31    if (startMos[N] < start) {
     32      Nlo = MAX(N, 0);
     33    } else {
     34      Nhi = MIN(N + 1, Nmosaic);
     35    }
     36  }
     37
     38  // check for the matched mosaic starting from Nlo
     39  // we may have to go much beyond Nlo since stop is not sorted
     40  // can we use a sorted version of stop to check when we are beyond the valid range??
     41  for (N = Nlo; N < Nmosaic; N++) {
     42    if (stop  < stopMos[N]) continue;
     43    if (start > startMos[N])  continue;
     44    if (stop  < startMos[N]) return (-1);
     45    return (indexMos[N]);
     46  }
     47
     48  return (-1);
     49}
     50
     51// sort two times vectors and an index by first time vector
     52void sort_mosaic_times (unsigned int *S, unsigned int *E, off_t *I, off_t N) {
     53
     54# define SWAPFUNC(A,B){ unsigned int tmp_t; off_t tmp_i; \
     55  tmp_t = S[A]; S[A] = S[B]; S[B] = tmp_t; \
     56  tmp_t = E[A]; E[A] = E[B]; E[B] = tmp_t; \
     57  tmp_i = I[A]; I[A] = I[B]; I[B] = tmp_i; \
     58}
     59# define COMPARE(A,B)(S[A] < S[B])
     60
     61  OHANA_SORT (N, COMPARE, SWAPFUNC);
     62
     63# undef SWAPFUNC
     64# undef COMPARE
     65
     66}
     67
     68// first, let's continue to use the time to make the match, but use bracketing to make it faster:
     69
    2170// find mosaic frames (unique time periods & photcode name matches mosaic)
    2271void initMosaics (Image *image, off_t Nimage) {
    2372
    24   off_t i, j, found, NMOSAIC;
     73  off_t i, Nmos, NMOSAIC;
    2574  unsigned int start, stop;
     75
     76  unsigned int *startMos, *stopMos;
     77  off_t *indexMos;
    2678
    2779  Nmosaic = 0;
    … …  
    3284  ALLOCATE (Amosaic_own_images, off_t, NMOSAIC);
    3385  ALLOCATE (mosaic_own_images, off_t *, NMOSAIC);
     86  ALLOCATE (startMos, unsigned int, NMOSAIC);
     87  ALLOCATE (stopMos, unsigned int, NMOSAIC);
     88  ALLOCATE (indexMos, off_t, NMOSAIC);
    3489
    3590  /* find the mosaic images (coords.ctype = DIS); generate list of unique mosaics */
    … …  
    56111    ALLOCATE (mosaic_own_images[Nmosaic], off_t, Amosaic_own_images[Nmosaic]);
    57112
     113    startMos[Nmosaic] = start;
     114    stopMos[Nmosaic] = stop;
     115    indexMos[Nmosaic] = Nmosaic;
     116
    58117    Nmosaic ++;
    59118    if (Nmosaic == NMOSAIC) {
    … …  
    63122      REALLOCATE (Nmosaic_own_images, off_t, NMOSAIC);
    64123      REALLOCATE (Amosaic_own_images, off_t, NMOSAIC);
    65     }
    66   }
    67 
     124      REALLOCATE (startMos, unsigned int, NMOSAIC);
     125      REALLOCATE (stopMos, unsigned int, NMOSAIC);
     126      REALLOCATE (indexMos, off_t, NMOSAIC);
     127    }
     128  }
     129
     130  // sort the index, start, and stop by the start times:
     131  sort_mosaic_times (startMos, stopMos, indexMos, Nmosaic);
     132 
    68133  // array to store image->mosaic index
    69134  Nmosaic_for_images = Nimage;
    … …  
    80145    stop  = image[i].tzero + MAX(1.01*image[i].trate*image[i].NY, 1);
    81146
    82     /* find existing mosaic with this time range */
    83     found = FALSE;
    84     for (j = 0; !found && (j < Nmosaic); j++) {
    85       if (stop  < mosaic[j].start) continue;
    86       if (start > mosaic[j].stop)  continue;
    87       found = TRUE;
    88       break;
    89     }
    90     /* if no matching mosaic exists, skip this image */
    91     if (!found) continue;
     147    Nmos = getMosaicByTimes (start, stop, startMos, stopMos, indexMos);
     148    if (Nmos == -1) continue;
    92149
    93150    // mosaic corresponding to this image
    94     mosaic_for_images[i] = j;
     151    mosaic_for_images[i] = Nmos;
    95152
    96153    // add image to mosaic_own_image list
    97     mosaic_own_images[j][Nmosaic_own_images[j]] = i;
    98     Nmosaic_own_images[j] ++;
    99     if (Nmosaic_own_images[j] == Amosaic_own_images[j]) {
    100       Amosaic_own_images[j] += 10;
    101       REALLOCATE (mosaic_own_images[j], off_t, Amosaic_own_images[j]);
    102     }
    103     assert (Nmosaic_own_images[j] < Amosaic_own_images[j]);
     154    mosaic_own_images[Nmos][Nmosaic_own_images[Nmos]] = i;
     155    Nmosaic_own_images[Nmos] ++;
     156    if (Nmosaic_own_images[Nmos] == Amosaic_own_images[Nmos]) {
     157      Amosaic_own_images[Nmos] += 10;
     158      REALLOCATE (mosaic_own_images[Nmos], off_t, Amosaic_own_images[Nmos]);
     159    }
     160    assert (Nmosaic_own_images[Nmos] < Amosaic_own_images[Nmos]);
    104161  }
    105162

  • branches/eam_branches/largefiles.20100314/Ohana/src/relastro/src/load_images.c

    r27312 r27425  
    11# include "relastro.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, SkyRegion *region) {
    … …  
    612  off_t      Nimage, Nsubset;
    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)
    … …  
    1625  // determine the populated SkyRegions overlapping the requested area
    1726  skylist = SkyListByPatch (sky, -1, region);
     27  MARKTIME("  setup sky: %f sec\n", dtime);
    1828
    1929  // convert database table to internal structure
    2030  image = gfits_table_get_Image (&db[0].ftable, &Nimage, &db[0].swapped);
     31  MARKTIME("  convert image table: %f sec\n", dtime);
    2132
    2233  // select the images which overlap the selected sky regions
    2334  subset = select_images (skylist, image, Nimage, &LineNumber, &Nsubset);
     35  MARKTIME("  select images: %f sec\n", dtime);
    2436
    2537  gfits_vtable_from_ftable (&db[0].ftable, &db[0].vtable, LineNumber, Nsubset);
    2638
    2739  initImages (subset, Nsubset);
     40  MARKTIME("  init images: %f sec\n", dtime);
     41
    2842  initMosaics (subset, Nsubset);
     43  MARKTIME("  init mosaics: %f sec\n", dtime);
    2944 
    3045  /* unlock, if we can (else, unlocked below) */

  • branches/eam_branches/largefiles.20100314/Ohana/src/relastro/src/relastro.c

    r17243 r27425  
    11# include "relastro.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
    39int main (int argc, char **argv) {
    … …  
    612  Catalog *catalog;
    713  FITS_DB db;
     14  struct timeval start, stop;
    815
    916  SkyList *skylist = NULL;
     17
     18  gettimeofday (&start, (void *) NULL);
    1019
    1120  /* get configuration info, args */
    … …  
    2635  if (db.dbstate == LCK_EMPTY) Shutdown ("ERROR: No images in catalog %s (1)", db.filename);
    2736  if (!dvo_image_load (&db, VERBOSE, FALSE)) Shutdown ("can't read image catalog %s", db.filename);
     37  MARKTIME("load image data: %f sec\n", dtime);
    2838
    2939  /* load regions and images based on specified sky patch (default depth) */
    3040  skylist = load_images (&db, &UserPatch);
     41  MARKTIME("load images: %f sec\n", dtime);
    3142
    3243  /* load catalog data from region files : subselect high-quality measurements */
    3344  catalog = load_catalogs (skylist, &Ncatalog, TRUE);
     45  MARKTIME("load catalog data: %f sec\n", dtime);
    3446
    3547  /* match measurements with images */
    3648  initImageBins (catalog, Ncatalog);
     49  MARKTIME("make image bins: %f sec\n", dtime);
     50
    3751  findImages (catalog, Ncatalog);
     52  MARKTIME("set up image indexes: %f sec\n", dtime);
    3853
    3954  if (PLOTSTUFF) {

  • branches/eam_branches/largefiles.20100314/Ohana/src/relastro/src/select_images.c

    r27312 r27425  
    1313} SkyRegionCoords;
    1414
     15void dsortindex (double *X, off_t *Y, int N);
     16off_t getRegionStartByRA (double R, double *Rref, off_t Nregions);
     17
    1518Image *select_images (SkyList *skylist, Image *timage, off_t Ntimage, off_t **LineNumber, off_t *Nimage) {
    1619 
    1720  Image *image;
    18   off_t i, j, k, m, nimage, NIMAGE;
     21  off_t i, j, k, m, nStart, iSky, nimage, NIMAGE;
    1922  int InRange, found;
    2023  double Ri[5], Di[5], Xi[5], Yi[5], dx, dy;
    … …  
    2326  SkyRegionCoords *skycoords;
    2427 
     28  double *RmaxSky;
     29  off_t *index;
     30
    2531  if (skylist[0].Nregions < 1) {
    2632    *Nimage = 0;
    … …  
    3844
    3945  ALLOCATE (skycoords, SkyRegionCoords, skylist[0].Nregions);
     46
     47  ALLOCATE (RmaxSky, double, skylist[0].Nregions);
     48  ALLOCATE (index, off_t, skylist[0].Nregions);
    4049
    4150  /* compare with each region file */
    … …  
    5665    skycoords[i].Yc[4] = skycoords[i].Yc[0];   
    5766
     67    RmaxSky[i] = skylist[0].regions[i][0].Rmax;
     68    index[i] = i;
     69
    5870    dx = 0.02*(skycoords[i].Xc[2] - skycoords[i].Xc[0]);
    5971    dy = 0.02*(skycoords[i].Yc[2] - skycoords[i].Yc[0]);
    … …  
    6476    skycoords[i].Xc[4] -= dx; skycoords[i].Yc[4] -= dy;
    6577  }
     78
     79  dsortindex (RmaxSky, index, skylist[0].Nregions);
    6680
    6781  if (VERBOSE) fprintf (stderr, "finding images\n");
    … …  
    100114    }
    101115   
    102     if (!FindMosaicForImage (timage, Ntimage, i)) continue;
     116    if (!FindMosaicForImage (timage, Ntimage, i)) {
     117      fprintf (stderr, "cannot find mosaic for %lld\n", (long long) i);
     118      continue;
     119    }
    103120
    104121    /* define image corners */
    … …  
    111128
    112129    /* transform to ra,dec */
     130    double RminImage = 360.0;
     131    // double RmaxImage =   0.0;
    113132    for (j = 0; j < 5; j++) {
    114133      XY_to_RD (&Ri[j], &Di[j], Xi[j], Yi[j], &timage[i].coords);
    115     }
     134      RminImage = MIN(RminImage, Ri[j]);
     135      // RmaxImage = MIN(RmaxImage, Ri[j]);
     136    }
     137
     138    // RA(nStart) is guaranteed to be < RminImage:
     139    nStart = getRegionStartByRA (RminImage, RmaxSky, skylist[0].Nregions);
    116140
    117141    /* compare with each region file */
    118     for (m = 0; (m < skylist[0].Nregions) && !found; m++) {
     142    for (iSky = nStart; (iSky < skylist[0].Nregions) && !found; iSky++) {
     143
     144      m = index[iSky];
     145      // if (RmaxImage < skylist[0].regions[m][0].Rmin) {
     146      // break;
     147      // }
    119148
    120149      /* we make positional comparisons in the projection of catalog */
    … …  
    122151      tcoords.crval2 = skycoords[m].Dc;
    123152
    124       /* transform to ra,dec */
     153      /* transform corner coords to X,Y in this catalog system */
    125154      InRange = TRUE;
    126155      for (j = 0; (j < 5) && InRange; j++) {
    … …  
    253282}
    254283
     284void dsortindex (double *X, off_t *Y, int N) {
     285
     286# define SWAPFUNC(A,B){ double tmpf; off_t tmpi; \
     287  tmpf = X[A]; X[A] = X[B]; X[B] = tmpf; \
     288  tmpi = Y[A]; Y[A] = Y[B]; Y[B] = tmpi; \
     289}
     290# define COMPARE(A,B)(X[A] < X[B])
     291
     292  OHANA_SORT (N, COMPARE, SWAPFUNC);
     293
     294# undef SWAPFUNC
     295# undef COMPARE
     296
     297}
     298
     299off_t getRegionStartByRA (double R, double *Rref, off_t Nregions) {
     300
     301  // use bisection to find the overlapping mosaic
     302
     303  off_t Nlo, Nhi, N;
     304
     305  // find the last mosaic before start
     306  Nlo = 0; Nhi = Nregions;
     307  while (Nhi - Nlo > 10) {
     308    N = 0.5*(Nlo + Nhi);
     309    if (Rref[N] < R) {
     310      Nlo = MAX(N, 0);
     311    } else {
     312      Nhi = MIN(N, Nregions);
     313    }
     314  }
     315  return (Nlo);
     316}
Note: See TracChangeset for help on using the changeset viewer.