IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Feb 13, 2011, 11:30:02 AM (15 years ago)
Author:
eugene
Message:

various fixes to photdbc, relastro and relphot based on analysis of the first reference database

Location:
trunk/Ohana/src/relphot
Files:
12 edited
1 copied

Legend:

Unmodified
Added
Removed

  • trunk/Ohana/src/relphot/Makefile

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

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

    r27586 r30616  
    1313  float Mcal;
    1414  float dMcal;
     15  float dMsys;
     16  short nFitPhotom;
    1517  short Xm;
    1618  float secz;
    … …  
    5052double STAR_CHISQ;
    5153double MIN_ERROR;
     54double IMFIT_SYS_SIGMA_LIM;
    5255
    5356int    VERBOSE;
    … …  
    103106SkyRegion UserPatch;
    104107int UserPatchSelect;
     108
     109int USE_BASIC_CHECK;
    105110
    106111# ifdef GRID_V1
    … …  
    215220void          write_coords        PROTO((Header *header, Coords *coords));
    216221int relphot_objects (void);
     222
     223void relphot_usage (void);
     224void relphot_help (int argc, char **argv);

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

    r25757 r30616  
    2020  GetConfig (config, "MAG_LIM",                "%lf", 0, &MAG_LIM);
    2121  GetConfig (config, "SIGMA_LIM",              "%lf", 0, &SIGMA_LIM);
     22
     23  if (!ScanConfig (config, "RELPHOT_IMFIT_SYS_SIGMA_LIM", "%lf", 0, &IMFIT_SYS_SIGMA_LIM)) IMFIT_SYS_SIGMA_LIM = 0.01;
     24
    2225  GetConfig (config, "STAR_SCATTER",           "%lf", 0, &STAR_SCATTER);
    2326  GetConfig (config, "IMAGE_SCATTER",          "%lf", 0, &IMAGE_SCATTER);

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

    r29001 r30616  
    124124  free (clist);
    125125  free (mlist);
     126  free (Nlist);
     127  free (NLIST);
    126128}
    127129
    … …  
    299301
    300302  off_t i, j, m, c, n, N, Nmax, mark, bad;
    301   int Nfew, Nbad, Nmos, Ngrid, Nrel, Nsys;
    302   float Msys, Mrel, Mmos, Mgrid;
     303  int Nfew, Nbad, Nmos, Ngrid, Nrel, Nsys, Nbright;
     304  float Msys, Mrel, Mmos, Mgrid, McalBright, McalBright2;
    303305  double *list, *dlist;
    304306  StatType stats;
    … …  
    318320
    319321  Nfew = Nbad = Nmos = Ngrid = Nrel = Nsys = 0;
     322
     323  // counters to measure the bright-end scatter
     324  McalBright = McalBright2 = 0.0;
     325  Nbright = 0;
    320326
    321327  for (i = 0; i < Nimage; i++) {
    … …  
    361367      list[N] = Msys - Mrel - Mmos - Mgrid;
    362368      dlist[N] = MAX (catalog[c].measure[m].dM, MIN_ERROR);
     369      if (catalog[c].measure[m].dM < IMFIT_SYS_SIGMA_LIM) {
     370        McalBright += list[N];
     371        McalBright2 += SQ(list[N]);
     372        Nbright ++;
     373      }
    363374      N++;
    364375    }
    … …  
    378389    liststats (list, dlist, N, &stats);
    379390    image[i].Mcal  = stats.mean;
    380     image[i].dMcal = stats.sigma;
     391    image[i].dMcal = stats.error;
     392    image[i].dMagSys = stats.sigma;
     393    image[i].nFitPhotom = N;
    381394    image[i].Xm    = 100.0*log10(stats.chisq);
    382395  }

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

    r29001 r30616  
    7373    mosaic[Nmosaic].Mcal  = 0.0;
    7474    mosaic[Nmosaic].dMcal = 0.0;
     75    mosaic[Nmosaic].dMsys = 0.0;
    7576    mosaic[Nmosaic].Xm    = 0.0;
    7677    mosaic[Nmosaic].flags  = image[i].flags;
    … …  
    189190      image[im].dMcal = mosaic[i].dMcal;
    190191      image[im].Xm = mosaic[i].Xm;
     192      image[im].dMagSys = mosaic[i].dMsys;
     193      image[im].nFitPhotom = mosaic[i].nFitPhotom;
    191194      image[im].flags |= (mosaic[i].flags & ID_IMAGE_PHOTOM_FEW);
    192195      image[im].flags |= (mosaic[i].flags & ID_IMAGE_PHOTOM_POOR);
    … …  
    400403      mark = (N < IMAGE_TOOFEW) || (N < IMAGE_GOOD_FRACTION*Nlist[i]);
    401404      if (mark) {
    402         fprintf (stderr, "marked mosaic %s ("OFF_T_FMT"), ("OFF_T_FMT" < %d) || ("OFF_T_FMT" < %f*"OFF_T_FMT")\n", image[imlist[i][0]].name,  i,  N, IMAGE_TOOFEW,  N, IMAGE_GOOD_FRACTION,  Nlist[i]);
     405        if (VERBOSE2) { fprintf (stderr, "marked mosaic %s ("OFF_T_FMT"), ("OFF_T_FMT" < %d) || ("OFF_T_FMT" < %f*"OFF_T_FMT")\n", image[imlist[i][0]].name,  i,  N, IMAGE_TOOFEW,  N, IMAGE_GOOD_FRACTION,  Nlist[i]); }
    403406        mosaic[i].flags |= ID_IMAGE_PHOTOM_FEW;
    404407        Nfew ++;
    … …  
    410413    if (PoorImages) fprintf (stderr, "Mmos: %f %f %d "OFF_T_FMT"\n", stats.mean, stats.sigma, stats.Nmeas,  N);
    411414    mosaic[i].Mcal  = stats.mean;
    412     mosaic[i].dMcal = stats.sigma;
     415    mosaic[i].dMcal = stats.error;
     416    mosaic[i].dMsys = stats.sigma;
     417    mosaic[i].nFitPhotom = N;
    413418    mosaic[i].Xm    = 100.0*log10(stats.chisq);
    414419  }

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

    r27586 r30616  
    11# include "relphot.h"
    2 void usage (void);
    32
    43int args (int argc, char **argv) {
    … …  
    5352  }
    5453
     54  USE_BASIC_CHECK = FALSE;
     55  if ((N = get_argument (argc, argv, "-basic-image-search"))) {
     56    remove_argument (N, &argc, argv);
     57    USE_BASIC_CHECK = TRUE;
     58  }
     59
    5560  VERBOSE = VERBOSE2 = FALSE;
    5661  if ((N = get_argument (argc, argv, "-v"))) {
    … …  
    228233  if (UpdateAverages && (argc == 1)) return TRUE;
    229234  if (UserPatchSelect && (argc == 2)) return TRUE;
    230   if (argc != 3) usage ();
     235  if (argc != 3) relphot_usage ();
    231236
    232237  return TRUE;
    233238}
    234 
    235 void usage () {
    236   fprintf (stderr, "ERROR: USAGE: relphot (region) (photcode)\n");
    237   fprintf (stderr, "       or:    relphot (photcode) -region RA RA DEC DEC\n");
    238   fprintf (stderr, "       or:    relphot -averages -region RA RA DEC DEC\n");
    239   fprintf (stderr, "  options: \n");
    240   fprintf (stderr, "  -time (start) (stop)\n");
    241   fprintf (stderr, "  -v : verbose output\n");
    242   fprintf (stderr, "  -vv : more verbose output\n");
    243   fprintf (stderr, "  -outroot (outroot)\n");
    244   fprintf (stderr, "  -plot\n");
    245   fprintf (stderr, "  -plotdelay (seconds)\n");
    246   fprintf (stderr, "  -statmode (mode)\n");
    247   fprintf (stderr, "  -refcode (name) : give extra weight to this photcode\n");
    248   fprintf (stderr, "  -n (nloop)\n");
    249   fprintf (stderr, "  -reset\n");
    250   fprintf (stderr, "  -update\n");
    251   fprintf (stderr, "  -params\n");
    252   fprintf (stderr, "  -mosaic (mosaic)\n");
    253   fprintf (stderr, "  -imfreeze\n");
    254   fprintf (stderr, "  -grid\n");
    255   fprintf (stderr, "  -area Xmin Xmax Ymin Ymax\n");
    256   fprintf (stderr, "  -instmag min max\n");
    257   fprintf (stderr, "  \n");
    258   exit (2);
    259 }

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

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

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

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

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

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

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

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

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

    r27480 r30616  
    8585  if (PLOTSTUFF) {
    8686    plot_star_coords (catalog, Ncatalog);
    87     plot_mosaic_fields (catalog);
     87    // plot_mosaic_fields (catalog);
    8888  }
    8989

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

    r29001 r30616  
    1515void dsortindex (double *X, off_t *Y, int N);
    1616off_t getRegionStartByRA (double R, double *Rref, off_t Nregions);
     17
     18# define MARKTIME(MSG,...) { \
     19  float dtime; \
     20  gettimeofday (&stop, (void *) NULL); \
     21  dtime = DTIME (stop, start); \
     22  fprintf (stderr, MSG, __VA_ARGS__); }
    1723
    1824Image *select_images (SkyList *skylist, Image *timage, off_t Ntimage, off_t **LineNumber, off_t *Nimage) {
    … …  
    2531  Coords tcoords;
    2632  SkyRegionCoords *skycoords;
    27  
     33  struct timeval start, stop;
     34 
     35  double RmaxSkyRegion, RminSkyRegion, DminSkyRegion, DmaxSkyRegion, RmidSkyRegion;
     36
    2837  double *RmaxSky;
    2938  off_t *index;
    … …  
    3544    return NULL;
    3645  }
     46
     47  gettimeofday (&start, (void *) NULL);
    3748
    3849  // the comparison is made in the catalog local projection. below we set crval1,2
    … …  
    4758  ALLOCATE (RmaxSky, double, skylist[0].Nregions);
    4859  ALLOCATE (index, off_t, skylist[0].Nregions);
     60
     61  RminSkyRegion = +360.0;
     62  RmaxSkyRegion = -360.0;
     63  DminSkyRegion = +90.0;
     64  DmaxSkyRegion = -90.0;
    4965
    5066  /* compare with each region file */
    … …  
    7591    skycoords[i].Xc[3] -= dx; skycoords[i].Yc[3] += dy;
    7692    skycoords[i].Xc[4] -= dx; skycoords[i].Yc[4] -= dy;
    77   }
     93
     94    RminSkyRegion = MIN(RminSkyRegion, skylist[0].regions[i][0].Rmin);
     95    RmaxSkyRegion = MAX(RmaxSkyRegion, skylist[0].regions[i][0].Rmax);
     96    DminSkyRegion = MIN(DminSkyRegion, skylist[0].regions[i][0].Dmin);
     97    DmaxSkyRegion = MAX(DmaxSkyRegion, skylist[0].regions[i][0].Dmax);
     98  }
     99  RmidSkyRegion = 0.5*(RminSkyRegion + RmaxSkyRegion);
     100  MARKTIME("create sky region coords: %f sec\n", dtime);
    78101
    79102  dsortindex (RmaxSky, index, skylist[0].Nregions);
     103  MARKTIME("sort sky coords: %f sec\n", dtime);
    80104
    81105  if (VERBOSE) fprintf (stderr, "finding images\n");
    82106  BuildChipMatch (timage, Ntimage);
     107  MARKTIME("build chip match: %f sec\n", dtime);
    83108
    84109  nimage = 0;
    … …  
    100125    }
    101126   
     127    // this adds 1.3 sec for 3M images
    102128    if (!FindMosaicForImage (timage, Ntimage, i)) {
    103129      fprintf (stderr, "cannot find mosaic for "OFF_T_FMT"\n", i);
    … …  
    113139    found = FALSE;
    114140
    115     /* transform corners to ra,dec */
     141    /* transform corners to ra,dec -- costs ~3sec for 3M images */
    116142    double RminImage = 360.0;
     143    double RmaxImage =   0.0;
     144    double DminImage = +90.0;
     145    double DmaxImage = -90.0;
    117146    for (j = 0; j < 5; j++) {
    118147      XY_to_RD (&Ri[j], &Di[j], Xi[j], Yi[j], &timage[i].coords);
     148      Ri[j] = ohana_normalize_angle_to_midpoint (Ri[j], RmidSkyRegion);
     149     
    119150      RminImage = MIN(RminImage, Ri[j]);
    120     }
    121 
    122     // RA(nStart) is guaranteed to be < RminImage:
     151      RmaxImage = MAX(RmaxImage, Ri[j]);
     152      DminImage = MIN(DminImage, Di[j]);
     153      DmaxImage = MAX(DmaxImage, Di[j]);
     154    }
     155    if (RmaxImage - RminImage > 180.0) {
     156        double tmp = RminImage;
     157        RmaxImage = RminImage;
     158        RminImage = tmp - 360.0;
     159    }
     160   
     161    // check that this image is even in range of the searched region
     162    if (DminImage > DmaxSkyRegion) continue;
     163    if (DmaxImage < DminSkyRegion) continue;
     164   
     165    // the sky region RA is defined to be 0 - 360.0
     166    if (RminImage > RmaxSkyRegion) continue;
     167    if (RmaxImage < RminSkyRegion) continue;
     168
     169    // image overlaps region, keep it
     170    if (USE_BASIC_CHECK) goto found_it;
     171
     172    // RA(nStart) is guaranteed to be < RminImage: -- costs 0.5sec for 3M images
    123173    nStart = getRegionStartByRA (RminImage, RmaxSky, skylist[0].Nregions);
    124174
    125175    /* compare with each region file */
    126     for (iSky = 0; (iSky < skylist[0].Nregions) && !found; iSky++) {
     176    for (iSky = nStart; (iSky < skylist[0].Nregions) && !found; iSky++) {
    127177
    128178      m = index[iSky];
    … …  
    142192      for (j = 0; (j < 4) && !found; j++) {
    143193        found = corner_check (&Xi[j], &Yi[j], &skycoords[m].Xc[0], &skycoords[m].Yc[0]);
     194        if (found) goto found_it;
    144195      }
    145196      /* check if catalog corner inside image */
    146197      for (j = 0; (j < 4) && !found; j++) {
    147198        found = corner_check (&skycoords[m].Xc[j], &skycoords[m].Yc[j], &Xi[0], &Yi[0]);
     199        if (found) goto found_it;
    148200      }
    149201      /* check if edges cross */
    … …  
    151203        for (k = 0; (k < 4) && !found; k++) {
    152204          found = edge_check (&Xi[j], &Yi[j], &skycoords[m].Xc[k], &skycoords[m].Yc[k]);
     205          if (found) goto found_it;
    153206        }
    154207      }
    155       if (!found) continue;
    156 
    157       image[nimage] = timage[i];
    158       /* always allow 'few' images to succeed, if possible */
    159       if (image[nimage].flags & ID_IMAGE_PHOTOM_FEW) {
    160         image[nimage].flags &= ~(ID_IMAGE_PHOTOM_FEW | ID_IMAGE_PHOTOM_POOR);
    161       }
    162       if (RESET) {
    163         assignMcal (&image[nimage], (double *) NULL, -1);
    164         image[nimage].dMcal = NAN;
    165         image[nimage].flags &= ~ID_IMAGE_PHOTOM_POOR;
    166       }
    167       line_number[nimage] = i;
    168       nimage ++;
    169       if (nimage == NIMAGE) {
    170         NIMAGE += 100;
    171         REALLOCATE (image, Image, NIMAGE);
    172         REALLOCATE (line_number, off_t, NIMAGE);
    173       }
    174     }
    175   }
    176      
     208    }
     209    if (!found) continue;
     210
     211  found_it:
     212    image[nimage] = timage[i];
     213    /* always allow 'few' images to succeed, if possible */
     214    if (image[nimage].flags & ID_IMAGE_PHOTOM_FEW) {
     215      image[nimage].flags &= ~(ID_IMAGE_PHOTOM_FEW | ID_IMAGE_PHOTOM_POOR);
     216    }
     217    if (RESET) {
     218      // XXX assignMcal (&image[nimage], (double *) NULL, -1);
     219      // XXX this needs to be thought through a bit more
     220      image[nimage].Mcal = 0.0;
     221      image[nimage].dMcal = NAN;
     222      image[nimage].flags &= ~ID_IMAGE_PHOTOM_POOR;
     223    }
     224    line_number[nimage] = i;
     225    nimage ++;
     226    if (nimage == NIMAGE) {
     227      NIMAGE += 100;
     228      REALLOCATE (image, Image, NIMAGE);
     229      REALLOCATE (line_number, off_t, NIMAGE);
     230    }
     231  }
     232  MARKTIME("finish image selection: %f sec\n", dtime);
     233
    177234  if (VERBOSE) fprintf (stderr, "found "OFF_T_FMT" images\n", nimage);
    178235
    … …  
    295352  return (Nlo);
    296353}
     354
     355off_t getRegionStopByRA (double R, double *Rref, off_t Nregions) {
     356
     357  // use bisection to find the overlapping mosaic
     358
     359  off_t Nlo, Nhi, N;
     360
     361  // find the last mosaic before start
     362  Nlo = 0; Nhi = Nregions;
     363  while (Nhi - Nlo > 10) {
     364    N = 0.5*(Nlo + Nhi);
     365    if (Rref[N] < R) {
     366      Nlo = MAX(N, 0);
     367    } else {
     368      Nhi = MIN(N, Nregions);
     369    }
     370  }
     371  return (Nlo);
     372}
Note: See TracChangeset for help on using the changeset viewer.