IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 27480


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

add detection->image, image->mosaic optimizations from relastro

Location:
trunk/Ohana/src/relphot/src
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • trunk/Ohana/src/relphot/src/ImageOps.c

    r27435 r27480  
    11# include "relphot.h"
    22
     3# define USE_IMAGE_ID 1
     4
     5static off_t       **bin;     // link from catalog,measure to image         
     6static off_t       **clist;   // catalog which supplied measurement on image
     7static off_t       **mlist;   // measure reference for measurement on image 
     8static off_t        *Nlist;   // number of measurements on image             
     9static off_t        *NLIST;   // allocated number of measurements on image   
     10
     11static Image        *image;   // list of available images 
     12static off_t        Nimage;   // number of available images
     13
     14// if we search by image ID, we sort (imageIDs, imageIdx) by imageIDs to get a sorted
     15// index
     16
     17# if USE_IMAGE_ID
     18static off_t        *imageIDs; // list of all image IDs
     19static off_t        *imageIdx; // list of index for image IDs
     20# else
    321static unsigned int *start;
    422static unsigned int *stop;
    5 static off_t       **bin;
    6 
    7 static off_t       **clist;
    8 static off_t       **mlist;
    9 static off_t        *Nlist;
    10 static off_t        *NLIST;
    11 
    12 static Image        *image;
    13 static off_t        Nimage;
     23# endif
     24
     25Image *getimages (off_t *N) {
     26
     27  *N = Nimage;
     28  return (image);
     29}
     30
     31Image *getimage (off_t N) {
     32  return (&image[N]);
     33}
    1434
    1535void initImages (Image *input, off_t N) {
    … …  
    2040  Nimage = N;
    2141
     42# if USE_IMAGE_ID
     43  ALLOCATE (imageIDs, off_t, Nimage);
     44  ALLOCATE (imageIdx, off_t, Nimage);
     45
     46  for (i = 0; i < Nimage; i++) {
     47    imageIdx[i] = i;
     48    imageIDs[i] = image[i].imageID;
     49  }
     50  llsortpair (imageIDs, imageIdx, Nimage);
     51# else
    2252  ALLOCATE (start,   unsigned, Nimage);
    2353  ALLOCATE (stop,    unsigned, Nimage);
    … …  
    2757    stop[i]  = image[i].tzero + MAX(1.01*image[i].trate*image[i].NY, 1);
    2858  }
     59# endif
     60}
     61
     62off_t getImageByID (off_t ID) {
     63
     64  // we have a pair of vectors (imageIDs, imageIdx) sorted by imageIDs
     65  // use bisection to find the specified image ID
     66
     67# if USE_IMAGE_ID
     68  off_t Nlo, Nhi, N;
     69
     70  Nlo = 0; Nhi = Nimage;
     71  while (Nhi - Nlo > 10) {
     72    N = 0.5*(Nlo + Nhi);
     73    if (imageIDs[N] < ID) {
     74      Nlo = MAX(N, 0);
     75    } else {
     76      Nhi = MIN(N + 1, Nimage);
     77    }
     78  }
     79
     80  for (N = Nlo; N < Nhi; N++) {
     81    if (imageIDs[N] == ID)
     82      return (imageIdx[N]);
     83  }
     84# endif
     85
     86  return (-1);
    2987}
    3088
    … …  
    67125  free (mlist);
    68126}
    69 
    70 /* **** use this to mimic the Mosaic version ***/
    71 # if (0)
    72 void initImageGrid (Image *image, off_t Nimage) {
    73 
    74   /* find max dR, dD range for all mosaics */
    75   /* define mosaic.coords to cover dR, dD */
    76   /* send results to initGridBins */
    77 
    78   off_t i, j, m, NX, NY;
    79   int dXmax, dYmax;
    80   double dS, dX, dY;
    81   double R, D, Rmin, Rmax, Dmin, Dmax;
    82   double Mcal, dMcal, Xm;
    83 
    84   dXmax = dYmax = 0.0;
    85   for (i = 0; i < Nmosaic; i++) {
    86     Dmin = Rmin =  360.0;
    87     Dmax = Rmax = -360.0;
    88     dS = 0.0;
    89     Mcal = dMcal = Xm = 0;
    90     for (j = 0; j < Nimlist[i]; j++) {
    91       m = imlist[i][j];
    92       NX = image[m].NX;
    93       NY = image[m].NY;
    94       dS += hypot (image[m].coords.cdelt1*image[m].coords.pc1_1, image[m].coords.cdelt1*image[m].coords.pc2_1);
    95       XY_to_RD (&R, &D, 0.0, 0.0, &image[m].coords);
    96       Rmin = MIN (Rmin, R);
    97       Rmax = MAX (Rmax, R);
    98       Dmin = MIN (Dmin, D);
    99       Dmax = MAX (Dmax, D);
    100       XY_to_RD (&R, &D, (double) NX, 0.0, &image[m].coords);
    101       Rmin = MIN (Rmin, R);
    102       Rmax = MAX (Rmax, R);
    103       Dmin = MIN (Dmin, D);
    104       Dmax = MAX (Dmax, D);
    105       XY_to_RD (&R, &D, (double) NX, (double) NY, &image[m].coords);
    106       Rmin = MIN (Rmin, R);
    107       Rmax = MAX (Rmax, R);
    108       Dmin = MIN (Dmin, D);
    109       Dmax = MAX (Dmax, D);
    110       XY_to_RD (&R, &D, 0.0, (double) NY, &image[m].coords);
    111       Rmin = MIN (Rmin, R);
    112       Rmax = MAX (Rmax, R);
    113       Dmin = MIN (Dmin, D);
    114       Dmax = MAX (Dmax, D);
    115       Mcal  += image[m].Mcal;
    116       dMcal += image[m].dMcal;
    117       Xm    += image[m].Xm;
    118       /* we are using mosaic.Mcal, not image.Mcal. reset image.Mcal */
    119       image[m].Mcal  = 0.0;
    120       image[m].dMcal = NAN;
    121       image[m].Xm    = NAN_S_SHORT;
    122     }
    123     dS /= Nimlist[i];
    124     strcpy (mosaic[i].coords.ctype, "RA---TAN");
    125     mosaic[i].coords.crval1 = Rmin;
    126     mosaic[i].coords.crval2 = Dmin;
    127     mosaic[i].coords.crpix1 = 0.0;
    128     mosaic[i].coords.crpix2 = 0.0;
    129     mosaic[i].coords.cdelt1 = dS;
    130     mosaic[i].coords.cdelt2 = dS;
    131     mosaic[i].coords.pc1_1  = 1.0;
    132     mosaic[i].coords.pc2_2  = 1.0;
    133     mosaic[i].coords.pc1_2  = 0.0;
    134     mosaic[i].coords.pc2_1  = 0.0;
    135     RD_to_XY (&dX, &dY, Rmax, Dmax, &mosaic[i].coords);
    136 
    137     mosaic[i].Mcal  = Mcal / Nimlist[i];
    138     mosaic[i].dMcal = dMcal / Nimlist[i];
    139     mosaic[i].Xm    = Xm / Nimlist[i];
    140   }
    141   if (!USE_GRID) return;
    142 
    143   dXmax = MAX (dXmax, dX);
    144   dYmax = MAX (dYmax, dY);
    145   initGrid (dXmax, dYmax);
    146   return;
    147 }
    148 # endif
    149127
    150128/* select all image equivalent to the current photcode */
    … …  
    163141}
    164142
    165 /* XXX convert this to use the values of measure->Xccd,Yccd */
    166 /* XXX need to apply the conversion for ELIXIR and LONEOS formats */
    167 void matchImage (Catalog *catalog, off_t meas, int cat) {
    168 
    169   off_t i;
     143int findCCD (off_t idx, off_t meas, int cat, Measure *measure) {
     144
    170145  int ccdnum;
    171146  double X, Y;
    172147  char *pname, *filter, *p, base[256];
     148
     149  /* identify the ccd on the basis of the photcode name */
     150  pname = GetPhotcodeNamebyCode (image[idx].photcode);
     151  filter = photcode[0].name;
     152  sprintf (base, "%s.%s.", MOSAICNAME, filter);
     153  if (strncmp (pname, base, strlen (base))) return (FALSE);
     154  p = pname + strlen(base);
     155
     156  /* p is pointing at the DETECTOR ID part of the photcode: CAMERA.FILTER.DETECTOR
     157     for now, we will let this be of the form SSSSnn where SSS is an arbitrary string
     158     and nn is an integer -> ccdnum */
     159     
     160  /// XXX we have an inconsistency here wrt dvo.layout and dvo.photcode
     161
     162  while (!isdigit(*p) && *p) p++;
     163  if (*p == 0) return (FALSE);
     164  ccdnum = atoi (p);
     165
     166  /* ccdnum is an integer, but not necessarily a sequence number.  when the camera layout is
     167     constructed, there will be null values for undefined ccdnums */
     168
     169  // old code to add this measurement to the grid cell for this chip
     170  // ave = measure[0].averef;
     171  // ra  = catalog[cat].average[ave].R - measure[0].dR / 3600.0;
     172  // dec = catalog[cat].average[ave].D - measure[0].dD / 3600.0;
     173  // RD_to_XY (&X, &Y, ra, dec, &image[i].coords);
     174
     175  // XXX we can now use these values (but need to be careful about old formats)
     176  X = measure[0].Xccd;
     177  Y = measure[0].Yccd;
     178  setGridMeasure (meas, cat, X, Y, ccdnum);
     179
     180  return (TRUE);
     181}
     182
     183# if USE_IMAGE_ID
     184void matchImage (Catalog *catalog, off_t meas, int cat) {
     185
     186  off_t idx, ID;
     187  int status;
     188  Measure *measure;
     189 
     190  measure = &catalog[cat].measure[meas];
     191
     192  ID = measure[0].imageID;
     193  idx = getImageByID (ID);
     194  if (idx == -1) {
     195    fprintf (stderr, "can't match detection to image?\n");
     196    abort();
     197  }
     198
     199  if (USE_GRID) {
     200    status = findCCD (idx, meas, cat, measure);
     201    if (!status) {
     202      fprintf (stderr, "failed to determine CCD for %lld, %d\n", (long long) meas, cat);
     203      abort();
     204    }
     205  }
     206
     207  // index for (catalog, measure) -> image
     208  bin[cat][meas] = idx;
     209
     210  // index for image, Nentry -> catalog
     211  clist[idx][Nlist[idx]] = cat;
     212
     213  // index for image, Nentry -> measure
     214  mlist[idx][Nlist[idx]] = meas;
     215  Nlist[idx] ++;
     216
     217  if (Nlist[idx] == NLIST[idx]) {
     218    NLIST[idx] += 100;
     219    REALLOCATE (clist[idx], off_t, NLIST[idx]);
     220    REALLOCATE (mlist[idx], off_t, NLIST[idx]);
     221  }     
     222
     223  return;
     224}
     225
     226# else
     227// this is the time-based match
     228void matchImage (Catalog *catalog, off_t meas, int cat) {
     229
     230  off_t i;
     231  int status;
    173232  Measure *measure;
    174233 
    … …  
    181240   
    182241    if (USE_GRID) {
    183 
    184       /* identify the ccd on the basis of the photcode name */
    185       pname = GetPhotcodeNamebyCode (image[i].photcode);
    186       filter = photcode[0].name;
    187       sprintf (base, "%s.%s.", MOSAICNAME, filter);
    188       if (strncmp (pname, base, strlen (base))) continue;
    189       p = pname + strlen(base);
    190 
    191       /* p is pointing at the DETECTOR ID part of the photcode: CAMERA.FILTER.DETECTOR
    192          for now, we will let this be of the form SSSSnn where SSS is an arbitrary string
    193          and nn is an integer -> ccdnum */
    194      
    195       /// XXX we have an inconsistency here wrt dvo.layout and dvo.photcode
    196 
    197       while (!isdigit(*p) && *p) p++;
    198       if (*p == 0) continue;
    199       ccdnum = atoi (p);
    200 
    201       /* ccdnum is an integer, but not necessarily a sequence number.  when the camera layout is
    202          constructed, there will be null values for undefined ccdnums */
    203 
    204       // old code to add this measurement to the grid cell for this chip
    205       // ave = measure[0].averef;
    206       // ra  = catalog[cat].average[ave].R - measure[0].dR / 3600.0;
    207       // dec = catalog[cat].average[ave].D - measure[0].dD / 3600.0;
    208       // RD_to_XY (&X, &Y, ra, dec, &image[i].coords);
    209 
    210       // XXX we can now use these values (but need to be careful about old formats)
    211       X = measure[0].Xccd;
    212       Y = measure[0].Yccd;
    213       setGridMeasure (meas, cat, X, Y, ccdnum);
    214     }
    215 
     242      status = findCCD(i, meas, cat);
     243      if (!status) continue;
     244    }
     245
     246    // index for (catalog, measure) -> image
    216247    bin[cat][meas] = i;
    217248
     249    // index for image, Nentry -> catalog
    218250    clist[i][Nlist[i]] = cat;
     251
     252    // index for image, Nentry -> measure
    219253    mlist[i][Nlist[i]] = meas;
    220254    Nlist[i] ++;
    … …  
    229263  /*  fprintf (stderr, "can't find source image for this measurement: %d (%d)\n", measure[0].t, measure[0].photcode); */
    230264}
     265# endif
    231266
    232267off_t getImageEntry (off_t meas, int cat) {
    … …  
    585620}
    586621
    587 Image *getimages (off_t *N) {
    588 
    589   *N = Nimage;
    590   return (image);
    591 }
    592 
    593 Image *getimage (off_t N) {
    594   return (&image[N]);
    595 }
    596 

  • trunk/Ohana/src/relphot/src/MosaicOps.c

    r27435 r27480  
    11# include "relphot.h"
    22
     3// array of mosaic definition structures
    34static off_t  Nmosaic;
    45static Mosaic *mosaic;
    56
    6 static off_t  *Nimlist;
     7// list of all images associated with a mosaic
     8static off_t   *Nimlist;
    79static off_t   **imlist; /* mosaic -> image[] */
    810static off_t   **bin;    /* catalog, measure -> mosaic */
    911
     12// list of mosaic associated with each image 
    1013static int     **clist;  /* mosaic -> catalog[] */
    1114static off_t   **mlist;  /* mosiac -> measure[] */
    … …  
    2932  ALLOCATE (Nimlist, off_t,   NMOSAIC);
    3033  ALLOCATE (NIMLIST, off_t,   NMOSAIC);
     34
     35  /* a 'mosaic' in relphot is (unlike relastro) a virtual concept: there is no
     36   * entry in the image table that represents this mosaic.  Instead, it is an
     37   * internal construct that defines a group of related images
     38   */
    3139
    3240  /* generate list of unique mosaics */

  • trunk/Ohana/src/relphot/src/relphot.c

    r21153 r27480  
    11# include "relphot.h"
    22# define USE_DIRECT 0
     3
     4# define MARKTIME(MSG,...) { \
     5  float dtime; \
     6  gettimeofday (&stop, (void *) NULL); \
     7  dtime = DTIME (stop, start); \
     8  fprintf (stderr, MSG, __VA_ARGS__); }
    39
    410int main (int argc, char **argv) {
    … …  
    713  Catalog *catalog;
    814  FITS_DB db;
     15  struct timeval start, stop;
     16  SkyList *skylist = NULL;
    917
    10   SkyList *skylist = NULL;
     18  gettimeofday (&start, (void *) NULL);
    1119
    1220  /* get configuration info, args */
    … …  
    3947    if (!dvo_image_load (&db, VERBOSE, FALSE)) Shutdown ("can't read image catalog %s", db.filename);
    4048  }
     49  MARKTIME("load image data: %f sec\n", dtime);
    4150
    4251  /* load regions and images based on specified sky patch */
    … …  
    4554  // XXX this is fairly lame: argv[1] is photcode if UserPatchSelect is true
    4655  skylist = load_images (&db, argv[1], &UserPatch, UserPatchSelect);
     56  MARKTIME("load images: %f sec\n", dtime);
    4757
    4858  /* unlock, if we can (else, unlocked below) */
    … …  
    5161  /* load catalog data from region files */
    5262  catalog = load_catalogs (skylist, &Ncatalog);
     63  MARKTIME("load catalog data: %f sec\n", dtime);
    5364 
    5465  /* add in a loop over the catalogs calling dvo_catalog_chipcoords */
    … …  
    5667  /* match measurements with images, mosaics */
    5768  initImageBins  (catalog, Ncatalog);
     69  MARKTIME("make image bins: %f sec\n", dtime);
     70
    5871  initMosaicBins (catalog, Ncatalog);
    5972  initGridBins   (catalog, Ncatalog);
    … …  
    6174
    6275  findImages (catalog, Ncatalog);
     76  MARKTIME("set up image indexes: %f sec\n", dtime);
     77
    6378  findMosaics (catalog, Ncatalog);  /* also sets Grid values */
     79  MARKTIME("set up mosaic indexes: %f sec\n", dtime);
     80
    6481  SAVEPLOT = FALSE;
    6582

  • trunk/Ohana/src/relphot/src/select_images.c

    r27435 r27480  
    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;
    19   off_t *line_number, nimage, NIMAGE;
     21  off_t i, j, k, m, nStart, iSky, nimage, NIMAGE;
     22  off_t *line_number;
    2023  int InRange, ecode, found;
    2124  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]);
    … …  
    6577  }
    6678
     79  dsortindex (RmaxSky, index, skylist[0].Nregions);
     80
    6781  if (VERBOSE) fprintf (stderr, "finding images\n");
    6882  BuildChipMatch (timage, Ntimage);
    … …  
    86100    }
    87101   
    88     if (!FindMosaicForImage (timage, Ntimage, i)) continue;
     102    if (!FindMosaicForImage (timage, Ntimage, i)) {
     103      fprintf (stderr, "cannot find mosaic for %lld\n", (long long) i);
     104      continue;
     105    }
    89106
    90107    /* define image corners */
    … …  
    96113    found = FALSE;
    97114
    98     /* transform to ra,dec */
     115    /* transform corners to ra,dec */
     116    double RminImage = 360.0;
    99117    for (j = 0; j < 5; j++) {
    100118      XY_to_RD (&Ri[j], &Di[j], Xi[j], Yi[j], &timage[i].coords);
    101     }
     119      RminImage = MIN(RminImage, Ri[j]);
     120    }
     121
     122    // RA(nStart) is guaranteed to be < RminImage:
     123    nStart = getRegionStartByRA (RminImage, RmaxSky, skylist[0].Nregions);
    102124
    103125    /* compare with each region file */
    104     for (m = 0; (m < skylist[0].Nregions) && !found; m++) {
     126    for (iSky = 0; (iSky < skylist[0].Nregions) && !found; iSky++) {
     127
     128      m = index[iSky];
    105129
    106130      /* we make positional comparisons in the projection of catalog */
    … …  
    238262}
    239263
     264void dsortindex (double *X, off_t *Y, int N) {
     265
     266# define SWAPFUNC(A,B){ double tmpf; off_t tmpi; \
     267  tmpf = X[A]; X[A] = X[B]; X[B] = tmpf; \
     268  tmpi = Y[A]; Y[A] = Y[B]; Y[B] = tmpi; \
     269}
     270# define COMPARE(A,B)(X[A] < X[B])
     271
     272  OHANA_SORT (N, COMPARE, SWAPFUNC);
     273
     274# undef SWAPFUNC
     275# undef COMPARE
     276
     277}
     278
     279off_t getRegionStartByRA (double R, double *Rref, off_t Nregions) {
     280
     281  // use bisection to find the overlapping mosaic
     282
     283  off_t Nlo, Nhi, N;
     284
     285  // find the last mosaic before start
     286  Nlo = 0; Nhi = Nregions;
     287  while (Nhi - Nlo > 10) {
     288    N = 0.5*(Nlo + Nhi);
     289    if (Rref[N] < R) {
     290      Nlo = MAX(N, 0);
     291    } else {
     292      Nhi = MIN(N, Nregions);
     293    }
     294  }
     295  return (Nlo);
     296}
Note: See TracChangeset for help on using the changeset viewer.