IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 28354


Ignore:
Timestamp:
Jun 16, 2010, 12:54:06 PM (16 years ago)
Author:
eugene
Message:

fix a bug in merge_catalogs_old, some optimization to speed things up

Location:
trunk/Ohana/src/dvomerge/src
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • trunk/Ohana/src/dvomerge/src/dvo_image_merge_dbs.c

    r28241 r28354  
    11# include "dvomerge.h"
     2
     3void sort_IDmap (IDmapType *IDmap);
    24
    35// merge db2 into db1
    … …  
    2931  }
    3032
     33  // sort IDmap->old,new on the basis of IDmap->old:
     34  sort_IDmap (IDmap);
     35
    3136  if (!out[0].swapped) {
    3237    gfits_convert_Image ((Image *) out[0].ftable.buffer, sizeof(Image), Nout);
    … …  
    4348}
    4449
    45 // optimize with sort and bisection
     50// XXX isn't the map just ID_new = ID_old + offset ??
    4651off_t dvo_map_image_ID (IDmapType *IDmap, off_t oldID) {
    4752
    48   off_t i;
     53  // off_t i;
     54  //
     55  // for (i = 0; i < IDmap->Nmap; i++) {
     56  //   if (IDmap->old[i] != oldID) continue;
     57  //   return (IDmap->new[i]);
     58  // }
    4959
    50   for (i = 0; i < IDmap->Nmap; i++) {
    51     if (IDmap->old[i] != oldID) continue;
    52     return (IDmap->new[i]);
     60  off_t Nlo, Nhi, N;
     61
     62  // find the a close entry below desired ID
     63  Nlo = 0; Nhi = IDmap->Nmap;
     64  while (Nhi - Nlo > 10) {
     65    N = 0.5*(Nlo + Nhi);
     66    if (IDmap->old[N] < oldID) {
     67      Nlo = MAX(N, 0);
     68    } else {
     69      Nhi = MIN(N + 1, IDmap->Nmap);
     70    }
     71  }
     72
     73  // search for the desired ID starting from Nlo, give up at Nhi
     74  for (N = Nlo; N < Nhi; N++) {
     75    if (IDmap->old[N] < oldID) continue;
     76    if (IDmap->old[N] > oldID) return 0;
     77    return (IDmap->new[N]);
    5378  }
    5479  return 0;
    … …  
    7297  return TRUE;
    7398}
     99
     100// sort two times vectors and an index by first time vector
     101void sort_IDmap (IDmapType *IDmap) {
     102
     103# define SWAPFUNC(A,B){ off_t tmp_old, tmp_new;    \
     104  tmp_old = IDmap->old[A]; IDmap->old[A] = IDmap->old[B]; IDmap->old[B] = tmp_old; \
     105  tmp_new = IDmap->new[A]; IDmap->new[A] = IDmap->new[B]; IDmap->new[B] = tmp_new; \
     106}
     107# define COMPARE(A,B)(IDmap->old[A] < IDmap->old[B])
     108
     109  OHANA_SORT (IDmap->Nmap, COMPARE, SWAPFUNC);
     110
     111# undef SWAPFUNC
     112# undef COMPARE
     113
     114}
     115

  • trunk/Ohana/src/dvomerge/src/dvomergeUpdate.c

    r28342 r28354  
    6969     
    7070    // load in all of the tables from input for this region
    71     inlist = SkyListByBounds (insky, depth, outsky[0].regions[i].Rmin, outsky[0].regions[i].Rmax, outsky[0].regions[i].Dmin, outsky[0].regions[i].Dmax);
     71    // SkyListByBounds will return neighbor catalogs if the boundaries exactly match (due to rounding).  Since the regions are not infinitely small,
     72    // compare to a slightly reduced footprint
     73    float dPos = 2.0/3600.0;
     74    inlist = SkyListByBounds (insky, depth, outsky[0].regions[i].Rmin + dPos, outsky[0].regions[i].Rmax - dPos, outsky[0].regions[i].Dmin + dPos, outsky[0].regions[i].Dmax - dPos);
    7275    for (j = 0; j < inlist[0].Nregions; j++) {
    7376      if (VERBOSE) fprintf (stderr, "input : %s\n", inlist[0].regions[j][0].name);
    … …  
    8689      dvo_catalog_unlock (&incatalog);
    8790      dvo_catalog_free (&incatalog);
     91
     92      fprintf (stderr, "merged %s into %s\n", outsky[0].regions[i].name, inlist[0].regions[j][0].name);
    8893    }
    8994    SkyListFree (inlist);

  • trunk/Ohana/src/dvomerge/src/merge_catalogs_old.c

    r28328 r28354  
    22# define PSPS_ID TRUE
    33
    4 # define IN_REGION(R,D) ( \
    5 ((D) >= region[0].Dmin) && ((D) < region[0].Dmax) && \
    6 ((R) >= region[0].Rmin)  && ((R) < region[0].Rmax))
     4# define MARKTIME(MSG,...) {                    \
     5    float dtime;                                \
     6    gettimeofday (&stop, (void *) NULL);        \
     7    dtime = DTIME (stop, start);                \
     8    fprintf (stderr, MSG, __VA_ARGS__);         \
     9    gettimeofday (&start, (void *) NULL);       \
     10  }
     11
     12# define IN_REGION(R,D) (                                       \
     13    ((D) >= region[0].Dmin) && ((D) < region[0].Dmax) &&        \
     14    ((R) >= region[0].Rmin)  && ((R) < region[0].Rmax))
    715
    816// merge the input data into the output catalog
    … …  
    1018int merge_catalogs_old (SkyRegion *region, Catalog *output, Catalog *input, double RADIUS) {
    1119
    12   off_t i, j, Nin, offset, J, Jmin, status, Nstars;
     20  off_t i, j, k, Nin, offset, J, Jmin, status, Nstars;
    1321  double RADIUS2, Rmin, Rin, Din;
    1422  double *X1, *Y1, *X2, *Y2;
    … …  
    1927  unsigned int objID, catID;
    2028  Coords tcoords;
     29
     30  // struct timeval start, stop;
     31  // gettimeofday (&start, (void *) NULL);
    2132
    2233  Nsecfilt = output[0].Nsecfilt;
    … …  
    108119  /* choose a radius for matches */
    109120  RADIUS2 = RADIUS*RADIUS;
     121
     122  // MARKTIME("set up structures: %f sec\n", dtime);
    110123
    111124  /** find matched stars **/
    … …  
    215228    // update the average properties to reflect the incoming entries:
    216229    // if the original value is NAN but the input value is not, accept the input:
    217     for (j = 0; j < Nsecfilt; j++) {
    218       if ( isfinite(output[0].secfilt[n*Nsecfilt+j].M)) continue;
    219       if (!isfinite( input[0].secfilt[N*Nsecfilt+j].M)) continue;
    220       output[0].secfilt[n*Nsecfilt+j].M     = input[0].secfilt[N*Nsecfilt+j].M;
    221       output[0].secfilt[n*Nsecfilt+j].dM    = input[0].secfilt[N*Nsecfilt+j].dM;
    222       output[0].secfilt[n*Nsecfilt+j].Xm    = input[0].secfilt[N*Nsecfilt+j].Xm;
    223       output[0].secfilt[n*Nsecfilt+j].M_20  = input[0].secfilt[N*Nsecfilt+j].M_20;
    224       output[0].secfilt[n*Nsecfilt+j].M_80  = input[0].secfilt[N*Nsecfilt+j].M_80;
    225       output[0].secfilt[n*Nsecfilt+j].Ncode = input[0].secfilt[N*Nsecfilt+j].Ncode;
    226       output[0].secfilt[n*Nsecfilt+j].Nused = input[0].secfilt[N*Nsecfilt+j].Nused;
     230    for (k = 0; k < Nsecfilt; k++) {
     231      if ( isfinite(output[0].secfilt[n*Nsecfilt+k].M)) continue;
     232      if (!isfinite( input[0].secfilt[N*Nsecfilt+k].M)) continue;
     233      output[0].secfilt[n*Nsecfilt+k] = input[0].secfilt[N*Nsecfilt+k];
    227234    }
    228235
    … …  
    232239    i++;
    233240  }
    234 
    235 # if (0)
    236   fprintf (stderr, "--- 1 ---\n");
    237   for (i = 0; i < Nmeas; i++) {
    238     fprintf (stderr, "Nave : %d, Nmeas : "OFF_T_FMT", dR: %f, dD: %f, catID: %d\n", output[0].measure[i].averef,  i, output[0].measure[i].dR, output[0].measure[i].dD, output[0].measure[i].catID);
    239   }
    240 # endif
     241  // MARKTIME("find matched stars: %f sec for "OFF_T_FMT","OFF_T_FMT" stars\n", dtime, Nstars, Nave);
    241242
    242243  /** incorporate unmatched image stars, if this star is in field of this catalog **/
    … …  
    260261    if (!IN_REGION (input[0].average[N].R, input[0].average[N].D)) continue;
    261262
    262     // XXX should we accept the input measurements for the fields?
     263    // XXX should we accept the input measurements for these fields?
    263264
    264265    output[0].average[Nave].R              = input[0].average[N].R;
    … …  
    344345    int Ngroup = input[0].average[N].Nmeasure;
    345346    for (j = 0; j < Ngroup - 1; j++) {
    346         next_meas[Nmeas - Ngroup + j] = Nmeas - Ngroup + j + 1;
     347      next_meas[Nmeas - Ngroup + j] = Nmeas - Ngroup + j + 1;
    347348    }
    348349    Nave ++;
    349350  }
    350351     
     352  // MARKTIME("save unmatched stars: %f sec\n", dtime);
     353
    351354  REALLOCATE (output[0].average, Average, Nave);
    352355  REALLOCATE (output[0].measure, Measure, Nmeas);
    353356 
    354 # if (0)
    355   fprintf (stderr, "--- 2 ---\n");
    356   for (i = 0; i < Nmeas; i++) {
    357     fprintf (stderr, "Nave : %d, Nmeas : "OFF_T_FMT", dR: %f, dD: %f, catID: %d\n", output[0].measure[i].averef,  i, output[0].measure[i].dR, output[0].measure[i].dD, output[0].measure[i].catID);
    358   }
    359 # endif
    360 
    361357# define NOSORT 0
    362358  if (NOSORT) {
    … …  
    366362    output[0].measure = sort_measure (output[0].average, Nave, output[0].measure, Nmeas, next_meas);
    367363  }
    368 
    369 # if (0)
    370   fprintf (stderr, "--- 3 ---\n");
    371   for (i = 0; i < Nmeas; i++) {
    372     fprintf (stderr, "Nave : %d, Nmeas : "OFF_T_FMT", dR: %f, dD: %f, catID: %d\n", output[0].measure[i].averef,  i, output[0].measure[i].dR, output[0].measure[i].dD, output[0].measure[i].catID);
    373   }
    374 # endif
    375364
    376365  /* note stars which have been found in this catalog */
    … …  
    400389  free (Y1);
    401390  free (N1);
     391
     392  // MARKTIME("cleanup: %f sec\n", dtime);
    402393  return (Nmatch);
    403394}
Note: See TracChangeset for help on using the changeset viewer.