IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Mar 24, 2010, 11:22:25 AM (16 years ago)
Author:
eugene
Message:

large update merging in changes for Ohana to support large files

Location:
trunk/Ohana
Files:
15 edited

Legend:

Unmodified
Added
Removed

  • trunk/Ohana

  • trunk/Ohana/src/relastro/doc/notes.txt

    r17043 r27435  
     1
     22010.03.24
     3
     4  I have finished only part of the optimizations listed below: I have
     5  made the detection->image and image->mosaic lookups use bracket
     6  searches.  I have NOT converted the image->mosaic lookups to use the
     7  parentID, so these matches are still based on time ranges.  I have
     8  also added a bit of positional constraint to the image match
     9  function in 'select_images' (comparing Rmax of the sky regions to
     10  Rmin of the images).  These have had a huge improvement.  Here is
     11  the execution time for the SAS data (~450 images) with the old vs
     12  the new code on pikake:
     13
     14  old : 39min
     15  new :  7min
     16
     17  The bulk of the time for the new code is now disk I/O.
     18
     192010.03.23
     20
     21  I've been working on optimization.  I've updated ImageOps to make
     22  the measure->image relationship generation use the image IDs along
     23  with a bracket search -- this goes much faster as a result.  More
     24  work is needed:
     25
     26  * select_images needs to apply astrometry to all image corners:
     27    -> cache the image center and radius and use this to narrow down
     28       the searches.
     29
     30  * MosaicOps uses the time to match mosaics.  I've added sorting to
     31    speed this up, but this should be done based on the image and
     32    parent index.
     33
     34    -> make a tool to create the parent ID for existing DBs.
     35    -> make sure addstar is populating the parent IDs
     36    -> use the parent IDs to make the link.
     37
     38load image data: 0.006347 sec
     39  setup sky: 0.091475 sec
     40  convert image table: 0.091498 sec
     41  select images: 2.925978 sec
     42  init images: 2.926957 sec
     43  init mosaics: 2.927319 sec
     44load images: 2.933762 sec
     45load catalog data: 9.772259 sec
     46make image bins: 9.780618 sec
     47set up image indexes: 10.196519 sec
     48
     49load image data: 0.006746 sec
     50  setup sky: 0.091440 sec
     51  convert image table: 0.091467 sec
     52  select images: 0.374193 sec
     53  init images: 0.375355 sec
     54  init mosaics: 0.375753 sec
     55load images: 0.382606 sec
     56load catalog data: 7.260094 sec
     57make image bins: 7.268564 sec
     58set up image indexes: 7.687485 sec
     59
    160
    2612008.03.01

  • trunk/Ohana/src/relastro/include/relastro.h

    r24308 r27435  
    140140char         *GetPhotnamebyCode   PROTO((PhotCodeData *photcodes, int code));
    141141void          InterpolateGrid     PROTO((float *buffer, int Nx, int Ny, Coords *ccd, Coords *gcoords));
    142 int          *SelectRefMosaic     PROTO((Mosaic **refmosaic, int *Nimage));
     142off_t          *SelectRefMosaic     PROTO((Mosaic **refmosaic, off_t *Nimage));
    143143int           args                PROTO((int argc, char **argv));
    144144int           bcatalog            PROTO((Catalog *subcatalog, Catalog *catalog));
    … …  
    153153void          findImages          PROTO((Catalog *catalog, int Ncatalog));
    154154int           findMosaics         PROTO((Catalog *catalog, int Ncatalog));
    155 Image        *find_images         PROTO((FITS_DB *db, GSCRegion *region, int Nregion, int *Nimage, int **LineNum));
     155Image        *find_images         PROTO((FITS_DB *db, GSCRegion *region, off_t Nregion, off_t *Nimage, off_t **LineNum));
    156156void set_db (FITS_DB *in);
    157 int Shutdown (char *format, ...);
     157int Shutdown (char *format, ...) OHANA_FORMAT(printf, 1, 2);
    158158void TrapSignal (int sig);
    159159void SetProtect (int mode);
    160160int SetSignals (void);
    161161
    162 GSCRegion    *find_regions        PROTO((Image *image, int Nimage, int *Nregions, GSCRegion *fullregion));
     162GSCRegion    *find_regions        PROTO((Image *image, off_t Nimage, int *Nregions, GSCRegion *fullregion));
    163163void          freeGridBins        PROTO((int Ncatalog));
    164164void          freeImageBins       PROTO((int Ncatalog));
    … …  
    166166void          free_catalogs       PROTO((Catalog *catalog, int Ncatalog));
    167167int           gcatalog            PROTO((Catalog *catalog, int FINAL));
    168 Coords       *getCoords           PROTO((int meas, int cat));
    169 float         getMcal             PROTO((int meas, int cat));
    170 float         getMgrid            PROTO((int meas, int cat));
    171 float         getMmos             PROTO((int meas, int cat));
    172 float         getMrel             PROTO((Catalog *catalog, int meas, int cat));
    173 GSCRegion    *get_regions         PROTO((double minRa, double maxRa, double minDec, double maxDec, int *Nregions));
    174 void          getfullregion       PROTO((Image *image, int Nimage, GSCRegion *fullregion));
    175 Image        *getimage            PROTO((int N));
    176 Image        *getimages           PROTO((int *N));
     168Coords       *getCoords           PROTO((off_t meas, int cat));
     169float         getMcal             PROTO((off_t meas, int cat));
     170float         getMgrid            PROTO((off_t meas, int cat));
     171float         getMmos             PROTO((off_t meas, int cat));
     172float         getMrel             PROTO((Catalog *catalog, off_t meas, int cat));
     173GSCRegion    *get_regions         PROTO((double minRa, double maxRa, double minDec, double maxDec, off_t *Nregions));
     174void          getfullregion       PROTO((Image *image, off_t Nimage, GSCRegion *fullregion));
     175Image        *getimage            PROTO((off_t N));
     176Image        *getimages           PROTO((off_t *N));
    177177void          global_stats        PROTO((Catalog *catalog, int Ncatalog));
    178178void          initGrid            PROTO((int dX, int dY));
    179179void          initGridBins        PROTO((Catalog *catalog, int Ncatalog));
    180180void          initImageBins       PROTO((Catalog *catalog, int Ncatalog));
    181 void          initImages          PROTO((Image *input, int N));
     181void          initImages          PROTO((Image *input, off_t N));
    182182void          initMosaicBins      PROTO((Catalog *catalog, int Ncatalog));
    183 void          initMosaicGrid      PROTO((Image *image, int Nimage));
    184 void          initMosaics         PROTO((Image *image, int Nimage));
     183void          initMosaicGrid      PROTO((Image *image, off_t Nimage));
     184void          initMosaics         PROTO((Image *image, off_t Nimage));
    185185void          initMrel            PROTO((Catalog *catalog, int Ncatalog));
    186186void          initialize          PROTO((int argc, char **argv));
    … …  
    189189Catalog      *load_catalogs       PROTO((SkyList *skylist, int *Ncatalog, int subselect));
    190190SkyList      *load_images         PROTO((FITS_DB *db, SkyRegion *region));
    191 Image        *select_images       PROTO((SkyList *skylist, Image *timage, int Ntimage, int **LineNumber, int *Nimage));
     191Image        *select_images       PROTO((SkyList *skylist, Image *timage, off_t Ntimage, off_t **LineNumber, off_t *Nimage));
    192192
    193193void check_permissions (char *basefile);
    … …  
    198198
    199199int           main                PROTO((int argc, char **argv));
    200 void          mark_images         PROTO((Image *image, int Nimage, Image *timage, int Ntimage));
    201 void          matchImage          PROTO((Catalog *catalog, int meas, int cat));
    202 void          matchMosaics        PROTO((Catalog *catalog, int meas, int cat));
    203 GSCRegion    *name_region         PROTO((char *name, int *Nregions));
     200void          mark_images         PROTO((Image *image, off_t Nimage, Image *timage, off_t Ntimage));
     201void          matchImage          PROTO((Catalog *catalog, off_t meas, int cat));
     202void          matchMosaics        PROTO((Catalog *catalog, off_t meas, int cat));
     203GSCRegion    *name_region         PROTO((char *name, off_t *Nregions));
    204204double        opening_angle       PROTO((double x1, double y1, double x2, double y2, double x3, double y3));
    205205void          plot_chisq          PROTO((Catalog *catalog, int Ncatalog));
    … …  
    224224int           setMrelOutput       PROTO((Catalog *catalog, int Ncatalog, int mark));
    225225void          set_ZP              PROTO((double ZERO));
    226 int           setrefcode          PROTO((Image *image, int Nimage));
     226int           setrefcode          PROTO((Image *image, off_t Nimage));
    227227void          skip_measurements   PROTO((Catalog *catalog, int pass));
    228228void          sortA               PROTO((double *X, int N));
    … …  
    268268int UpdateMosaic (Catalog *catalog, int Ncatalog);
    269269int UpdateMeasures (Catalog *catalog, int Ncatalog);
    270 void fixImageRaw (Catalog *catalog, int Ncatalog, int im);
     270void fixImageRaw (Catalog *catalog, int Ncatalog, off_t im);
    271271void FlagOutliers(Catalog *catalog);
    272272int MeasFilterTest(Measure *measure);
    … …  
    278278int FitPMandPar (PMFit *fit, double *X, double *dX, double *Y, double *dY, double *T, double *pR, double *pD, int Npts);
    279279
    280 Mosaic *getMosaicForImage (int N);
    281 
    282 StarData *getImageRef (Catalog *catalog, int Ncatalog, int im, int *Nstars, CoordMode mode);
    283 StarData *getImageRaw (Catalog *catalog, int Ncatalog, int im, int *Nstars, CoordMode mode);
    284 
    285 Mosaic *getmosaics (int *N);
    286 void initMosaics (Image *image, int Nimage);
    287 StarData *getMosaicRaw (Catalog *catalog, int Ncatalog, int mos, int *Nstars);
    288 StarData *getMosaicRef (Catalog *catalog, int Ncatalog, int mos, int *Nstars);
    289 Mosaic *getMosaicForImage (int im);
     280Mosaic *getMosaicForImage (off_t N);
     281
     282StarData *getImageRef (Catalog *catalog, int Ncatalog, off_t im, off_t *Nstars, CoordMode mode);
     283StarData *getImageRaw (Catalog *catalog, int Ncatalog, off_t im, off_t *Nstars, CoordMode mode);
     284
     285Mosaic *getmosaics (off_t *N);
     286void initMosaics (Image *image, off_t Nimage);
     287StarData *getMosaicRaw (Catalog *catalog, int Ncatalog, off_t mos, off_t *Nstars);
     288StarData *getMosaicRef (Catalog *catalog, int Ncatalog, off_t mos, off_t *Nstars);
     289Mosaic *getMosaicForImage (off_t im);
    290290
    291291double getMeanR (Measure *measure, Average *average, SecFilt *secfilt);
    … …  
    297297int relastro_objects (void);
    298298int UpdateObjectOffsets (SkyList *skylist);
     299
     300int relastroVisualPlotRawRef(StarData *raw, StarData *ref, double dRmax, int numObj);
     301int relastroVisualPlotScatter(double values[], double thresh, int npts);
     302int relastroVisualPlotOutliers(Catalog *catalog, int offset, int Nmeasure,
     303                               StatType statsR, StatType statsD, double thresh);
     304

  • trunk/Ohana/src/relastro/src/FitSimple.c

    r16810 r27435  
    1717
    1818    // XXX the value needs to be set in a more intelligent way
    19     if (dR > 0.15) continue;
     19    if (dR > 0.50) continue;
    2020   
    2121    fit_add (fit, raw[i].X, raw[i].Y, ref[i].P, ref[i].Q, 1.0);

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

    r24308 r27435  
    11# include "relastro.h"
    2 #define TESTING
    3 
     2# define TESTING
     3# define USE_IMAGE_ID 1
     4
     5static off_t       **bin;     // link from catalog,measure to image
     6static int         **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 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
    419static unsigned int *start;
    520static unsigned int *stop;
    6 static int         **bin;
    7 
    8 static int         **clist;
    9 static int         **mlist;
    10 static int          *Nlist;
    11 static int          *NLIST;
    12 
    13 static Image        *image;
    14 static int          Nimage;
    15 
    16 Image *getimages (int *N) {
     21# endif
     22
     23Image *getimages (off_t *N) {
    1724  *N = Nimage;
    1825  return (image);
    1926}
    2027
    21 Image *getimage (int N) {
     28Image *getimage (off_t N) {
    2229  return (&image[N]);
    2330}
    2431
    25 void initImages (Image *input, int N) {
    26 
    27   int i;
     32void initImages (Image *input, off_t N) {
     33
     34  off_t i;
    2835
    2936  image = input;
    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
    4187void initImageBins (Catalog *catalog, int Ncatalog) {
    4288
    43   int i, j;
    44 
    45   ALLOCATE (bin, int *, Ncatalog);
     89  off_t i, j;
     90
     91  ALLOCATE (bin, off_t *, Ncatalog);
    4692  for (i = 0; i < Ncatalog; i++) {
    47     ALLOCATE (bin[i], int, MAX (catalog[i].Nmeasure, 1));
     93    ALLOCATE (bin[i], off_t, MAX (catalog[i].Nmeasure, 1));
    4894    for (j = 0; j < catalog[i].Nmeasure; j++) bin[i][j] = -1;
    4995  }
    5096
    51   ALLOCATE (Nlist, int, Nimage);
    52   ALLOCATE (NLIST, int, Nimage);
    53   ALLOCATE (clist, int *, Nimage);
    54   ALLOCATE (mlist, int *, Nimage);
     97  ALLOCATE (Nlist, off_t,  Nimage);
     98  ALLOCATE (NLIST, off_t,  Nimage);
     99  ALLOCATE (clist, int *,   Nimage);
     100  ALLOCATE (mlist, off_t *, Nimage);
    55101
    56102  for (i = 0; i < Nimage; i++) {
    … …  
    58104    NLIST[i] = 100;
    59105    ALLOCATE (clist[i], int, NLIST[i]);
    60     ALLOCATE (mlist[i], int, NLIST[i]);
     106    ALLOCATE (mlist[i], off_t, NLIST[i]);
    61107  }
    62108}
    … …  
    64110void freeImageBins (int Ncatalog) {
    65111
    66   int i;
     112  off_t i;
    67113
    68114  for (i = 0; i < Ncatalog; i++) {
    … …  
    81127void findImages (Catalog *catalog, int Ncatalog) {
    82128
    83   int i, j;
     129  off_t i, j;
    84130  char *name;
    85131
    … …  
    92138  for (i = 0; VERBOSE && (i < Nimage); i++) {
    93139    name = GetPhotcodeNamebyCode (image[i].photcode);
    94     fprintf (stderr, "image %d has %d measures (%s, %s)\n", i, Nlist[i],
     140    fprintf (stderr, "image %lld has %lld measures (%s, %s)\n", (long long) i, (long long) Nlist[i],
    95141             ohana_sec_to_date(image[i].tzero), name);
    96142  }
    97143}
    98144
     145# if USE_IMAGE_ID
    99146/* modify this function to use the measure->imageID field */
    100 void matchImage (Catalog *catalog, int meas, int cat) {
    101 
    102   int i;
     147void matchImage (Catalog *catalog, off_t meas, int cat) {
     148
     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
     183  off_t i;
    103184  Measure *measure;
    104185
    … …  
    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;
    … …  
    124207    if (Nlist[i] == NLIST[i]) {
    125208      NLIST[i] += 100;
    126       REALLOCATE (clist[i], int, NLIST[i]);
    127       REALLOCATE (mlist[i], int, NLIST[i]);
     209      REALLOCATE (clist[i], int,   NLIST[i]);
     210      REALLOCATE (mlist[i], off_t, NLIST[i]);
    128211    }
    129212    return;
    130213  }
    131 }
    132 
    133 Coords *getCoords (int meas, int cat) {
    134 
    135   int i;
     214  fprintf (stderr, "can't match detection to image?\n");
     215  abort();
     216}
     217# endif
     218
     219Coords *getCoords (off_t meas, int cat) {
     220
     221  off_t i;
    136222
    137223  i = bin[cat][meas];
    … …  
    142228void plot_images () {
    143229
    144   int i, bin;
     230  off_t i, bin;
    145231  double *xlist, *Mlist, *dlist;
    146232  Graphdata graphdata;
    … …  
    188274// return StarData values for detections in the specified image, converting coordinates from the
    189275// chip positions: X,Y -> L,M -> P,Q -> R,D
    190 void fixImageRaw (Catalog *catalog, int Ncatalog, int im) {
    191 
    192   int i, m, c, n;
     276void fixImageRaw (Catalog *catalog, int Ncatalog, off_t im) {
     277
     278  off_t i, m, c, n;
    193279  double X, Y, L, M, P, Q, R, D, dR, dD;
    194280
    … …  
    199285  mosaic = getMosaicForImage (im);
    200286  if (mosaic != NULL) {
    201       moscoords = &mosaic[0].coords;
     287    moscoords = &mosaic[0].coords;
    202288  }
    203289  imcoords = &image[im].coords;
    … …  
    257343// chip positions: X,Y -> L,M -> P,Q -> R,D.  This function is used by the image fitting steps, for
    258344// which the detections have already been filtered when they were loaded (bcatalog)
    259 StarData *getImageRaw (Catalog *catalog, int Ncatalog, int im, int *Nstars, CoordMode mode) {
    260 
    261   int i, m, c, n;
     345StarData *getImageRaw (Catalog *catalog, int Ncatalog, off_t im, off_t *Nstars, CoordMode mode) {
     346
     347  off_t i, m, c, n;
    262348
    263349  Mosaic *mosaic;
    … …  
    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  }
    … …  
    332418// the sky positions: R,D -> P,Q -> L,M -> X,Y
    333419
    334 StarData *getImageRef (Catalog *catalog, int Ncatalog, int im, int *Nstars, CoordMode mode) {
    335 
    336   int i, m, c, n;
     420StarData *getImageRef (Catalog *catalog, int Ncatalog, off_t im, off_t *Nstars, CoordMode mode) {
     421
     422  off_t i, m, c, n;
    337423
    338424  Mosaic *mosaic;
    … …  
    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");
    … …  
    399485  return;
    400486
    401   int i, j, k, m, N, Ndel, Nave, Nmax, TOOFEW, Nsecfilt;
     487  int Ndel, Nave;
     488  off_t i, j, k, m, N, Nmax, TOOFEW, Nsecfilt;
    402489  double Ns, theta, x, y;
    403490  double *R, *D, *dR, *dD;
    … …  
    490577void FlagOutliers2D (Catalog *catalog) {
    491578
    492   int i, j, k, m, N, Ndel, Nave, Nmax, TOOFEW, Nsecfilt;
     579  int Ndel, Nave;
     580  off_t i, j, k, m, N, Nmax, TOOFEW, Nsecfilt;
    493581  double *index;
    494582  double Ns, theta, x, y;
    … …  
    569657    // recalculate image center, sigma based on closest 50% of points
    570658    for(k = 0;  k < N; k++) {
    571       int ind = (int) index[k];
     659      off_t ind = (off_t) index[k];
    572660      R[k] = catalog[0].measure[ind].dR;
    573661      D[k] = catalog[0].measure[ind].dD;

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

    r21508 r27435  
    22
    33// array of mosaic definition structures
    4 static int    Nmosaic;
     4static off_t    Nmosaic;
    55static Mosaic *mosaic;
    66
    77// list of all images associated with a mosaic
    8 static int   *Nmosaic_own_images; // number of images for this mosaic
    9 static int   *Amosaic_own_images; // size of allocated array
    10 static int   **mosaic_own_images; // array of arrays: mosaic -> images
     8static off_t   *Nmosaic_own_images; // number of images for this mosaic
     9static off_t   *Amosaic_own_images; // size of allocated array
     10static off_t   **mosaic_own_images; // array of arrays: mosaic -> images
    1111
    1212// list of mosaic associated with each image 
    13 static int    Nmosaic_for_images; // number of images (for internal checks)
    14 static int    *mosaic_for_images; // array of: image -> mosaic
    15 
    16 Mosaic *getmosaics (int *N) {
     13static off_t    Nmosaic_for_images; // number of images (for off_ternal checks)
     14static off_t    *mosaic_for_images; // array of: image -> mosaic
     15
     16Mosaic *getmosaics (off_t *N) {
    1717  *N = Nmosaic;
    1818  return (mosaic);
    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)
    22 void initMosaics (Image *image, int Nimage) {
    23 
    24   int i, j, found, NMOSAIC;
     71void initMosaics (Image *image, off_t Nimage) {
     72
     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;
    … …  
    2981  ALLOCATE (mosaic, Mosaic, NMOSAIC);
    3082
    31   ALLOCATE (Nmosaic_own_images, int, NMOSAIC);
    32   ALLOCATE (Amosaic_own_images, int, NMOSAIC);
    33   ALLOCATE (mosaic_own_images, int *, NMOSAIC);
     83  ALLOCATE (Nmosaic_own_images, off_t, NMOSAIC);
     84  ALLOCATE (Amosaic_own_images, off_t, NMOSAIC);
     85  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 */
    … …  
    54109    Nmosaic_own_images[Nmosaic] = 0;
    55110    Amosaic_own_images[Nmosaic] = 10;
    56     ALLOCATE (mosaic_own_images[Nmosaic], int, Amosaic_own_images[Nmosaic]);
     111    ALLOCATE (mosaic_own_images[Nmosaic], off_t, Amosaic_own_images[Nmosaic]);
     112
     113    startMos[Nmosaic] = start;
     114    stopMos[Nmosaic] = stop;
     115    indexMos[Nmosaic] = Nmosaic;
    57116
    58117    Nmosaic ++;
    … …  
    60119      NMOSAIC += 10;
    61120      REALLOCATE (mosaic, Mosaic, NMOSAIC);
    62       REALLOCATE (mosaic_own_images, int *, NMOSAIC);
    63       REALLOCATE (Nmosaic_own_images, int, NMOSAIC);
    64       REALLOCATE (Amosaic_own_images, int, NMOSAIC);
    65     }
    66   }
    67 
     121      REALLOCATE (mosaic_own_images, off_t *, NMOSAIC);
     122      REALLOCATE (Nmosaic_own_images, off_t, NMOSAIC);
     123      REALLOCATE (Amosaic_own_images, off_t, NMOSAIC);
     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;
    70   ALLOCATE (mosaic_for_images, int, Nmosaic_for_images);
     135  ALLOCATE (mosaic_for_images, off_t, Nmosaic_for_images);
    71136
    72137  /* now assign the WRP images to these mosaics */
    … …  
    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], int, 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
    … …  
    109166// return StarData values for detections in the specified image, converting coordinates from the
    110167// chip positions: X,Y -> L,M -> P,Q -> R,D
    111 StarData *getMosaicRaw (Catalog *catalog, int Ncatalog, int mos, int *Nstars) {
    112 
    113   int i, j, im, Nraw, Nnew;
     168StarData *getMosaicRaw (Catalog *catalog, int Ncatalog, off_t mos, off_t *Nstars) {
     169
     170  off_t i, j, im, Nraw, Nnew;
    114171  StarData *raw, *new;
    115172
    … …  
    142199// return StarData values for averages positions in the specified image, converting coordinates from
    143200// the sky positions: R,D -> P,Q -> L,M -> X,Y
    144 StarData *getMosaicRef (Catalog *catalog, int Ncatalog, int mos, int *Nstars) {
    145 
    146   int i, j, im, Nref, Nnew;
     201StarData *getMosaicRef (Catalog *catalog, int Ncatalog, off_t mos, off_t *Nstars) {
     202
     203  off_t i, j, im, Nref, Nnew;
    147204  StarData *ref, *new;
    148205 
    … …  
    172229}
    173230
    174 Mosaic *getMosaicForImage (int im) {
    175 
    176   int mos;
     231Mosaic *getMosaicForImage (off_t im) {
     232
     233  off_t mos;
    177234
    178235  if (im < 0) abort();

  • trunk/Ohana/src/relastro/src/UpdateChips.c

    r17151 r27435  
    44
    55  /* we can measure new image parameters for each non-mosaic chip independently */
    6   int i, Nimage, Nraw, Nref;
     6  off_t i, Nimage, Nraw, Nref;
    77  Image *image;
    88  StarData *raw, *ref;
    … …  
    2525
    2626    // FitChip does iterative, clipped fitting
    27     fprintf (stderr, "image %d : Nstars: %d\n", i, Nraw);
     27    fprintf (stderr, "image %lld : Nstars: %lld\n", (long long) i, (long long) Nraw);
    2828    FitChip (raw, ref, Nraw, &image[i].coords);
    2929

  • trunk/Ohana/src/relastro/src/UpdateMeasures.c

    r16810 r27435  
    33int UpdateMeasures (Catalog *catalog, int Ncatalog) {
    44
    5   int i, Nimage;
     5  off_t i, Nimage;
    66  Image *image;
    77

  • trunk/Ohana/src/relastro/src/UpdateMosaic.c

    r15600 r27435  
    44
    55  /* we can measure new image parameters for each mosaic independently */
    6   int i, Nmosaic, Nstars;
     6  off_t i, Nmosaic, Nstars;
    77  Mosaic *mosaic;
    88  StarData *raw, *ref;

  • trunk/Ohana/src/relastro/src/UpdateObjects.c

    r24308 r27435  
    11# include "relastro.h"
    22
    3 static int Nmax;
     3static off_t  Nmax;
    44static double *X, *dX;
    55static double *Y, *dY;
    … …  
    1313void initObjectData (Catalog *catalog, int Ncatalog) {
    1414
    15   int i, j;
     15  off_t i, j;
    1616 
    1717  Nmax = 0;
    … …  
    4040int UpdateObjects (Catalog *catalog, int Ncatalog) {
    4141
    42   int i, j, k, m, N, Nsecfilt, found, kp;
     42  int i, j, k, m, N, Nsecfilt;
    4343  StatType statsR, statsD;
    4444  Coords coords;
    … …  
    4747  int mode, Nave, Npm, Npar, Nskip;
    4848  double Tmin, Tmax;
    49   float mag;
    50   int mask;
    51   PhotCode *code;
    5249
    5350  initObjectData (catalog, Ncatalog);
    … …  
    7471  for (i = 0; i < Ncatalog; i++) {
    7572
    76     if (VERBOSE) fprintf (stderr, "astrometrize catalog %d : %d ave, %d meas\n", i, catalog[i].Naverage, catalog[i].Nmeasure);
     73    if (VERBOSE) fprintf (stderr, "astrometrize catalog %d : %lld ave, %lld meas\n", i, (long long) catalog[i].Naverage, (long long) catalog[i].Nmeasure);
    7774
    7875    Nskip = 0;

  • trunk/Ohana/src/relastro/src/UpdateSimple.c

    r15590 r27435  
    44
    55  /* we can measure new image parameters for each non-mosaic chip independently */
    6   int i, Nimage, Nstars;
     6  off_t i, Nimage, Nstars;
    77  Image *image;
    88  StarData *raw, *ref;

  • trunk/Ohana/src/relastro/src/bcatalog.c

    r24308 r27435  
    33int bcatalog (Catalog *subcatalog, Catalog *catalog) {
    44
    5   int i, j, k, offset, found;
    6   int NAVERAGE, NMEASURE, Naverage, Nmeasure, Nm, Nsecfilt;
    7   float mag;
    8   int mask;
    9   PhotCode *code;
     5  off_t i, j, offset;
     6  off_t NAVERAGE, NMEASURE, Naverage, Nmeasure, Nm;
     7  int Nsecfilt;
    108
    119  // XXX in the future, use catalog[0].Nsecfilt only?  allow catalogs to have variable Nsecfilt?
    … …  
    106104
    107105  if (VERBOSE) {
    108     fprintf (stderr, "%d: using %d stars (%d measures) for catalog\n", i,
    109              subcatalog[0].Naverage, subcatalog[0].Nmeasure);
     106    fprintf (stderr, "%lld: using %lld stars (%lld measures) for catalog\n", (long long) i, (long long) subcatalog[0].Naverage, (long long) subcatalog[0].Nmeasure);
    110107   }
    111108  return (TRUE);

  • trunk/Ohana/src/relastro/src/load_images.c

    r25757 r27435  
    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) {
    410
    511  Image     *image, *subset;
    6   int        Nimage, Nsubset;
    7   int       *LineNumber;
     12  off_t      Nimage, Nsubset;
     13  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) */

  • trunk/Ohana/src/relastro/src/relastro.c

    r17243 r27435  
    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) {

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

    r21508 r27435  
    1313} SkyRegionCoords;
    1414
    15 Image *select_images (SkyList *skylist, Image *timage, int Ntimage, int **LineNumber, int *Nimage) {
     15void dsortindex (double *X, off_t *Y, int N);
     16off_t getRegionStartByRA (double R, double *Rref, off_t Nregions);
     17
     18Image *select_images (SkyList *skylist, Image *timage, off_t Ntimage, off_t **LineNumber, off_t *Nimage) {
    1619 
    1720  Image *image;
    18   int i, j, k, m, found, nimage, NIMAGE;
    19   int InRange;
     21  off_t i, j, k, m, nStart, iSky, nimage, NIMAGE;
     22  int InRange, found;
    2023  double Ri[5], Di[5], Xi[5], Yi[5], dx, dy;
    21   int *line_number;
     24  off_t *line_number;
    2225  Coords tcoords;
    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);
    … …  
    7185  NIMAGE = 100;
    7286  ALLOCATE (image, Image, NIMAGE);
    73   ALLOCATE (line_number, int, NIMAGE);
     87  ALLOCATE (line_number, off_t, NIMAGE);
    7488 
    7589  // go through the complete list of images, selecting ones which overlap any region
    … …  
    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++) {
    … …  
    161190        NIMAGE += 100;
    162191        REALLOCATE (image, Image, NIMAGE);
    163         REALLOCATE (line_number, int, NIMAGE);
     192        REALLOCATE (line_number, off_t, NIMAGE);
    164193      }
    165194    }
    166195  }
    167196     
    168   if (VERBOSE) fprintf (stderr, "found %d images\n", nimage);
     197  if (VERBOSE) fprintf (stderr, "found %lld images\n", (long long) nimage);
    169198
    170199  REALLOCATE (image, Image, MAX (nimage, 1));
    171   REALLOCATE (line_number, int, MAX (nimage, 1));
     200  REALLOCATE (line_number, off_t, MAX (nimage, 1));
    172201  free (skycoords);
    173202
    … …  
    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.