IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 36833


Ignore:
Timestamp:
Jun 7, 2014, 6:33:11 AM (12 years ago)
Author:
eugene
Message:

merge changes from eam_branches/ipp-20140423: add dcr-color option, ref_color, bad detection clipping to relastro; fix RA,DEC for mosaics in imextract; allow photdbc to restrict images to output region; change PHOTCODE_SKIP_LIST to PHOTCODE_KEEP_LIST in photdbc; add -reset-to-photcode option to relastro; ensure relastro child processes get the complete set of options; allow relastro region hosts to continue after load-object clients are done, but before the data is loaded; invalidate extreme proper motion; limit relastro catalog density on a grid; require more detections for higher order fits; option to delete duplicate images based on exp time/date and photcode; have kapa use XSync instead of XFlush to avoid flashing; handle some cases to avoid crashing kapa; extend kapa connection timeout; allow kapa to handle > 2GB images; ensure dvomerge works even if incoming db has no images; add dvo fuctions coordmosaic & psastro_model

Location:
trunk/Ohana/src
Files:
68 edited
24 copied

Legend:

Unmodified
Added
Removed
  • trunk/Ohana/src/addstar/src/ReadImageHeader.c

    r31395 r36833  
    4444  image[0].NX = Nx;
    4545  image[0].NY = Ny;
     46
     47  image[0].refColor = NAN;
    4648
    4749  if (!gfits_scan (header, "TZERO",   "%d",  1, &image[0].tzero) && !ACCEPT_TIME) {
  • trunk/Ohana/src/delstar/include/delstar.h

    r35837 r36833  
    22# include <dvo.h>
    33# include <signal.h>
     4
     5# define MARKTIME(MSG,...) {                    \
     6    gettimeofday (&stopTimer, (void *) NULL);   \
     7    float dtime = DTIME (stopTimer, startTimer);        \
     8    fprintf (stderr, MSG, __VA_ARGS__); }
     9
     10# define INITTIME \
     11  struct timeval startTimer, stopTimer; \
     12  gettimeofday (&startTimer, (void *) NULL);
    413
    514// options for generating the IndexArray used to select images for deletion (delete_duplicate_images.c)
     
    6776int    UPDATE;
    6877int    IMAGE_DETAILS;
     78int    IMAGE_DUPLICATES_BY_OBSTIME;
    6979int    IMAGE_ONLY;
    7080int    ORPHAN;
  • trunk/Ohana/src/delstar/src/args.c

    r35758 r36833  
    1818  fprintf (stderr, "  -image-details : list info about the deleted images (-dup-images only)\n");
    1919  fprintf (stderr, "  -image-only : only examine the image table (changes are NOT saved; -dup-images only)\n");
     20  fprintf (stderr, "  -image-only-force : modify only the image table (-dup-images only)\n");
     21  fprintf (stderr, "  -image-by-obstime : use date/time and photcode (not externID) to find duplicates\n");
    2022  fprintf (stderr, "  -region Rmin Rmax Dmin Dmax : apply changes to this part of the sky\n");
    2123 
     
    6870  if ((N = get_argument (argc, argv, "-image-details"))) {
    6971    IMAGE_DETAILS = TRUE;
     72    remove_argument (N, &argc, argv);
     73  }
     74
     75  IMAGE_DUPLICATES_BY_OBSTIME = FALSE;
     76  if ((N = get_argument (argc, argv, "-image-by-obstime"))) {
     77    IMAGE_DUPLICATES_BY_OBSTIME = TRUE;
    7078    remove_argument (N, &argc, argv);
    7179  }
  • trunk/Ohana/src/delstar/src/delete_duplicate_images.c

    r35758 r36833  
    11# include "delstar.h"
     2
     3IndexArray *find_duplicates_obstime (Image *image, off_t Nimage, off_t *Nduplicates);
     4off_t find_obstime_range (Image *image, off_t Nimage, off_t firstEntry);
     5void sort_by_obstime (e_time *T, short *P, off_t *I, off_t N);
     6void sort_by_photcode (short *P, off_t *I, off_t N);
    27
    38// this function identifies the images to be deleted based on duplication of the
     
    1823
    1924  off_t Nduplicates = 0;
    20   IndexArray *imageID = find_duplicates (image, Nimage, &Nduplicates);
     25  IndexArray *imageID = NULL;
     26
     27  if (IMAGE_DUPLICATES_BY_OBSTIME) {
     28    imageID = find_duplicates_obstime (image, Nimage, &Nduplicates);
     29  } else {
     30    imageID = find_duplicates (image, Nimage, &Nduplicates);
     31  }
    2132
    2233  if (Nduplicates == 0) {
     
    494505}
    495506
     507// sort by increasing obstime
     508void sort_by_obstime (e_time *T, short *P, off_t *I, off_t N) {
     509
     510# define SWAPFUNC(A,B){ e_time tmpT; short tmpP; off_t tmpI; \
     511  tmpT = T[A]; T[A] = T[B]; T[B] = tmpT; \
     512  tmpP = P[A]; P[A] = P[B]; P[B] = tmpP; \
     513  tmpI = I[A]; I[A] = I[B]; I[B] = tmpI; \
     514}
     515# define COMPARE(A,B)(T[A] < T[B])
     516
     517  OHANA_SORT (N, COMPARE, SWAPFUNC);
     518
     519# undef SWAPFUNC
     520# undef COMPARE
     521
     522}
     523
     524// sort by increasing photcode
     525void sort_by_photcode (short *P, off_t *I, off_t N) {
     526
     527# define SWAPFUNC(A,B){ short tmpP; off_t tmpI; \
     528  tmpP = P[A]; P[A] = P[B]; P[B] = tmpP; \
     529  tmpI = I[A]; I[A] = I[B]; I[B] = tmpI; \
     530}
     531# define COMPARE(A,B)(P[A] < P[B])
     532
     533  OHANA_SORT (N, COMPARE, SWAPFUNC);
     534
     535# undef SWAPFUNC
     536# undef COMPARE
     537
     538}
     539
     540static e_time *obstime  = NULL;
     541static short  *photcode = NULL;
     542static off_t  *primary  = NULL;
     543static off_t  *idx      = NULL;
     544static char   *keep     = NULL;
     545
     546static short  *photcode_subset = NULL;
     547static off_t  *idx_subset      = NULL;
     548
     549static off_t Nsubset = 0;
     550static off_t NSUBSET = 300;
     551
     552off_t find_obstime_range (Image *image, off_t Nimage, off_t firstEntry) {
     553
     554  Nsubset = 0;
     555
     556  // find all entries with the same obstime:
     557  e_time firstTime = obstime[firstEntry];
     558
     559  idx_subset[Nsubset] = idx[firstEntry];
     560  photcode_subset[Nsubset] = photcode[firstEntry];
     561  Nsubset++;
     562
     563  off_t i = firstEntry + Nsubset;
     564  while ((i < Nimage) && (obstime[i] == firstTime)) {
     565    idx_subset[Nsubset] = idx[i];
     566    photcode_subset[Nsubset] = photcode[i];
     567    Nsubset++;
     568    i++;
     569    if (Nsubset >= NSUBSET) {
     570      NSUBSET += 1000;
     571      REALLOCATE (photcode_subset, short, NSUBSET);
     572      REALLOCATE (idx_subset, off_t, NSUBSET);
     573    }
     574  }     
     575
     576  sort_by_photcode (photcode_subset, idx_subset, Nsubset);
     577
     578  return i;
     579}
     580
     581// alternative version to find duplicates based on obstime and photcode
     582IndexArray *find_duplicates_obstime (Image *image, off_t Nimage, off_t *Nduplicates) {
     583
     584  if (Nimage < 1) return NULL;
     585
     586  // how to find duplicates:
     587  // generate a set of arrays (idx, obstime, photcode, keep)
     588  // sort the arrays by obstime
     589  // loop over obstime. 
     590  // for a given new obstime, scan through the entries with the same value
     591  // create a sub-array of photcodes, idx
     592  // sort by photcode
     593  // loop over entries
     594  // find matching photcodes
     595  // mark duplicate enties
     596
     597  ALLOCATE (obstime, e_time, Nimage);
     598  ALLOCATE (photcode, short, Nimage);
     599  ALLOCATE (primary, off_t, Nimage);
     600  ALLOCATE (idx, off_t, Nimage);
     601  ALLOCATE (keep, char, Nimage);
     602
     603  off_t i;
     604
     605  INITTIME;
     606
     607  // skip entries with photcode == 0?
     608  for (i = 0; i < Nimage; i++) {
     609    idx[i] = i;
     610    primary[i] = -1; // only duplicates get a value for primary
     611    keep[i] = TRUE;
     612    obstime[i] = image[i].tzero;
     613    photcode[i] = image[i].photcode;
     614  }
     615  MARKTIME("  generate index arrays: %f sec\n", dtime);
     616
     617  // sort the 4 arrays
     618  sort_by_obstime (obstime, photcode, idx, Nimage);
     619  MARKTIME("  sort index arrays: %f sec\n", dtime);
     620
     621  // image[idx[i]].tzero = obstime[i]
     622  // keep[i] -> keep image[i] ('keep' is NOT resorted)
     623
     624  // allocate arrays to store the subsets (these are global static)
     625  // These get reallocated if necessary in find_obstime_range()
     626  ALLOCATE (photcode_subset, short, NSUBSET);
     627  ALLOCATE (idx_subset, off_t, NSUBSET);
     628
     629  // entries of idx_subset correspond to the original image sequence:
     630  // image[idx_subset[i]].tzero = obstime[i]
     631
     632  off_t firstEntry = 0;
     633  off_t nextEntry = 0;
     634 
     635  while (nextEntry < Nimage) {
     636    if (firstEntry >= Nimage) {
     637      fprintf (stderr, "error, too far?\n");
     638    }
     639   
     640    // generate photcode_subset, idx_subset in order of photcode for this unique obstime[firstEntry]
     641    // returned value is the first value of the next entry
     642    nextEntry = find_obstime_range (image, Nimage, firstEntry);
     643   
     644    // step through the photcodes and find duplicates
     645    int j;
     646    int firstCodeEntry = 0;
     647    short firstCode = photcode_subset[firstCodeEntry];
     648    for (j = 1; j < Nsubset; j++) {
     649      if (photcode_subset[j] == firstCode) {
     650        // mark as duplicate
     651        off_t dupIndex = idx_subset[j];
     652        keep[dupIndex] = FALSE;
     653        primary[dupIndex] = idx_subset[firstCodeEntry];
     654      } else {
     655        // new value of photcode, call if the first one
     656        firstCodeEntry = j;
     657        firstCode = photcode_subset[firstCodeEntry];
     658      }
     659    }
     660    firstEntry = nextEntry;
     661  }
     662  MARKTIME("  find duplicates: %f sec\n", dtime);
     663
     664  IndexArray *imageID = make_index_array (image, Nimage, IMAGE_ID);
     665  MARKTIME("  make index array: %f sec\n", dtime);
     666
     667  // set imageID->value to TRUE for images we want to delete
     668  off_t Ndup = 0;
     669  for (i = 0; i < Nimage; i++) {
     670    if (keep[i]) continue;
     671    off_t Ni = image[i].imageID - imageID->minID;
     672    myAssert (Ni >= 0, "oops");
     673    myAssert (Ni < imageID->range, "oops");
     674    imageID->value[Ni] = TRUE;
     675    Ndup ++;
     676  }
     677  MARKTIME("  mark duplicates: %f sec\n", dtime);
     678
     679  for (i = 0; IMAGE_DETAILS && (i < Nimage); i++) {
     680    off_t Ni = image[i].imageID - imageID->minID;
     681    if (!imageID->value[Ni]) continue;
     682   
     683    char *date = NULL;
     684    date = ohana_sec_to_date (image[i].tzero);
     685    fprintf (stderr, "delete image : (" OFF_T_FMT "), extID = %d : %30s : %20s %5d  ==  ", i, image[i].externID, image[i].name, date, image[i].photcode);
     686    free (date);
     687
     688    off_t myPrimary = primary[i];
     689    if (myPrimary < 0) {
     690      fprintf (stderr, "ERROR: this should never happen\n");
     691      abort();
     692    }
     693    date = ohana_sec_to_date (image[myPrimary].tzero);
     694    fprintf (stderr, "parent image : (" OFF_T_FMT "), extID = %d : %30s : %20s %5d\n", myPrimary, image[myPrimary].externID, image[myPrimary].name, date, image[myPrimary].photcode);
     695    free (date);
     696  }
     697
     698  *Nduplicates = Ndup;
     699
     700  free (photcode_subset);
     701  free (idx_subset);
     702 
     703  free (obstime);
     704  free (photcode);
     705  free (primary);
     706  free (idx);
     707  free (keep);
     708
     709  return imageID;;
     710}
     711
    496712// find the min & max values of the given ID (externID or imageID)
    497713// construct an empty array with length needed to fit IDs
  • trunk/Ohana/src/dvomerge/src/dvo_image_merge_dbs.c

    r36528 r36833  
    168168
    169169  if (!IDmap->old) {
     170    fprintf (stderr, "input database has image IDs, but no Image table\n");
     171    return FALSE;
     172  }
     173
     174  if (!IDmap->Nmap) {
    170175    fprintf (stderr, "input database has image IDs, but no Image table\n");
    171176    return FALSE;
  • trunk/Ohana/src/dvomerge/src/dvomergeImageIDs.c

    r35765 r36833  
    1818  if (inDB.dbstate == LCK_EMPTY) {
    1919    dvo_image_unlock (&inDB); // unlock input
     20    IDmap->old = NULL;
     21    IDmap->new = NULL;
    2022    IDmap->Nmap = 0;
    2123    return TRUE;
     
    99101  if (inDB.dbstate == LCK_EMPTY) {
    100102    dvo_image_unlock (&inDB); // unlock input
     103    IDmap->old = NULL;
     104    IDmap->new = NULL;
    101105    IDmap->Nmap = 0;
    102106    return TRUE;
  • trunk/Ohana/src/dvomerge/src/dvorepairCPT.c

    r33657 r36833  
    126126    average[Nave].measureOffset = -1;
    127127    average[Nave].missingOffset = -1;
    128     average[Nave].extendOffset = -1;
     128    average[Nave].refColor = NAN;
    129129
    130130    average[Nave].objID = measure[i].objID;
  • trunk/Ohana/src/dvomerge/src/dvorepairDeleteImageList.c

    r29938 r36833  
    411411    average[Nave].measureOffset = -1;
    412412    average[Nave].missingOffset = -1;
    413     average[Nave].extendOffset = -1;
     413    average[Nave].refColor = NAN;
    414414
    415415    average[Nave].objID = measure[i].objID;
  • trunk/Ohana/src/dvomerge/src/dvorepairFixTables.c

    r29938 r36833  
    266266    average[Nave].measureOffset = -1;
    267267    average[Nave].missingOffset = -1;
    268     average[Nave].extendOffset = -1;
     268    average[Nave].refColor = NAN;
    269269
    270270    average[Nave].objID = measure[i].objID;
  • trunk/Ohana/src/dvopsps/src/insert_objects_dvopsps_catalog.c

    r35450 r36833  
    292292  PrintIOBuffer (ave_buffer, "%d, ",  average->measureOffset); 
    293293  PrintIOBuffer (ave_buffer, "%d, ",  average->missingOffset); 
    294   PrintIOBuffer (ave_buffer, "%d, ",  average->extendOffset);   
     294  PrintIOBuffer (ave_buffer, "%d, ",  0); // was : average->extendOffset);   
    295295  PrintIOBuffer (ave_buffer, "%u, ",  average->flags);         
    296296  PrintIOBuffer (ave_buffer, "%u, ",  average->photFlagsUpper);
  • trunk/Ohana/src/kapa2/src/CheckPipe.c

    r29938 r36833  
    11# include "Ximage.h"
     2# include <errno.h>
    23# define STRCONST(A) ((int)(0x1000000*A[0] + 0x10000*A[1] + 0x100*A[2] + 0x1*A[3]))
    34
     
    4142  word[4] = 0;
    4243  switch (status) {
    43   case -1:                          /* no input from pipe: continue */
    44     return (TRUE);
     44  case -1:
     45    if (errno == EAGAIN) {
     46      /* no input from pipe: continue */
     47      return (TRUE);
     48    }
     49    perror ("exiting due to problem with socket connection in CheckPipe");
     50    return (FALSE);
    4551    break;
    4652
  • trunk/Ohana/src/kapa2/src/EventLoop.c

    r21153 r36833  
    11# include "Ximage.h"
     2# define DEBUG 0
    23
    34/* list events being selected below, all other masks are ignored */
     
    67int LastEvent (Display *display, int type, XEvent *event) {
    78
    8   int found;
     9  int found, Nfound;
    910
    1011  found = FALSE;
     12  Nfound = 0;
    1113  while (XCheckTypedEvent (display, type, event)) {
     14    // If the link is slow, then I should wait a little while for some config events to
     15    // build up (typically, the window is being dragged around, so we get a whole series of
     16    // config events.  the flush / sync time is slow on a slow link, so we should try to wait
     17    // for the config events to stop, but on a timescale deteremined by the flush/sync time
     18    // if (type == ConfigureNotify) fprintf (stderr, "config (%d)\n", Nfound);
     19    Nfound ++;
    1220    found = TRUE;
     21  }
     22  if (DEBUG && found) {
     23    if (type == ConfigureNotify) fprintf (stderr, "config (%d)\n", Nfound);
     24    if (type == CirculateNotify) fprintf (stderr, "circul (%d)\n", Nfound);
     25    if (type == Expose) fprintf (stderr, "expose (%d)\n", Nfound);
     26  }
     27
     28  // if I have a ConfigureNotify event, I should purge all Expose events as well:
     29  if (found && (type == ConfigureNotify)) {
     30    XEvent discard;
     31    while (XCheckTypedEvent (display, Expose, &discard));
    1332  }
    1433  return (found);
     
    3857
    3958    if (XEventsQueued (display, QueuedAfterFlush) < 1) {
    40       /* fprintf (stderr, "."); */
    4159      usleep (50000);
    4260      continue;
    4361    }
    4462
     63    // If I have a config event, I want to also purge all expose events
     64
    4565    /* grab the last entry for these events */
    46     if (LastEvent (display, ConfigureNotify, &event)) Reconfig (&event);
    47     if (LastEvent (display, CirculateNotify, &event)) Reconfig (&event);
    48     if (LastEvent (display, Expose,          &event)) Refresh ();
     66    if (LastEvent (display, ConfigureNotify, &event)) { Reconfig (&event); continue; }
     67    if (LastEvent (display, CirculateNotify, &event)) { Reconfig (&event); continue; }
     68    if (LastEvent (display, Expose,          &event)) { Refresh (); continue; }
    4969    if (LastEvent (display, MappingNotify,   &event)) XRefreshKeyboardMapping ((XMappingEvent *) &event);
    5070    if (LastEvent (display, MotionNotify,    &event)) UpdatePointer (graphic, (XMotionEvent *) &event);
  • trunk/Ohana/src/kapa2/src/FlushDisplay.c

    r21153 r36833  
    3030
    3131  if (flush) {
    32     XFlush (graphic->display);
     32    // I changed XFlush to XSync to avoid lag
     33    // of the display (build up of config events)
     34    XSync (graphic->display, FALSE);
    3335    reftime = now;
    3436  }
  • trunk/Ohana/src/kapa2/src/Refresh.c

    r27790 r36833  
    4747  Y0 = graphic[0].dy - graphic[0].dy * section[0].y - dY;
    4848  X0 = graphic[0].dx * section[0].x;
    49   // fprintf (stderr, "section: (%s) %d %d - %d %d\n", section[0].name, X0, Y0, dX, dY);
    5049
    5150  XFillRectangle (graphic[0].display,  graphic[0].window, graphic[0].gc, X0, Y0, dX, dY);
  • trunk/Ohana/src/libautocode/def/average-ps1-dev-2.d

    r16810 r36833  
    2525FIELD measureOffset,  OFF_MEASURE, int,             offset to first psf measurement
    2626FIELD missingOffset,  OFF_MISSING, int,             offset to first missing obs
    27 FIELD extendOffset,   OFF_EXTEND,  int,             offset to first extended measurement
     27FIELD refColor,       REF_COLOR,   float,           color of astrometry ref stars
    2828
    2929FIELD code,           code,       unsigned short,   ID code (star; ghost; etc)
  • trunk/Ohana/src/libautocode/def/average-ps1-v1.d

    r21508 r36833  
    3030FIELD measureOffset,  OFF_MEASURE, uint32_t,        offset to first psf measurement
    3131FIELD missingOffset,  OFF_MISSING, uint32_t,        offset to first missing obs
    32 FIELD extendOffset,   OFF_EXTEND,  uint32_t,        offset to first extended measurement
     32FIELD refColor,       REF_COLOR,   float,           color of astrometry ref stars
    3333
    3434# 'flags' was called 'code' prior to 2009.02.07
  • trunk/Ohana/src/libautocode/def/average-ps1-v2.d

    r33647 r36833  
    3535FIELD measureOffset,  OFF_MEASURE, uint32_t,        offset to first psf measurement
    3636FIELD missingOffset,  OFF_MISSING, uint32_t,        offset to first missing obs
    37 FIELD extendOffset,   OFF_EXTEND,  uint32_t,        offset to first extended measurement
     37FIELD refColor,       REF_COLOR,   float,           color of astrometry ref stars
    3838
    3939# 'flags' was called 'code' prior to 2009.02.07
  • trunk/Ohana/src/libautocode/def/average-ps1-v3.d

    r33647 r36833  
    3535FIELD measureOffset,  OFF_MEASURE, uint32_t,        offset to first psf measurement
    3636FIELD missingOffset,  OFF_MISSING, uint32_t,        offset to first missing obs
    37 FIELD extendOffset,   OFF_EXTEND,  uint32_t,        offset to first extended measurement
     37FIELD refColor,       REF_COLOR,   float,           color of astrometry ref stars
    3838
    3939# 'flags' was called 'code' prior to 2009.02.07
  • trunk/Ohana/src/libautocode/def/average-ps1-v4.d

    r34405 r36833  
    3838FIELD measureOffset,  OFF_MEASURE, uint32_t,        offset to first psf measurement
    3939FIELD missingOffset,  OFF_MISSING, uint32_t,        offset to first missing obs
    40 FIELD extendOffset,   OFF_EXTEND,  uint32_t,        offset to first extended measurement
     40FIELD refColor,       REF_COLOR,   float,           color of reference stars
    4141
    4242# 'flags' was called 'code' prior to 2009.02.07
  • trunk/Ohana/src/libautocode/def/average.d

    r34405 r36833  
    3838FIELD measureOffset,  OFF_MEASURE, int,             offset to first psf measurement
    3939FIELD missingOffset,  OFF_MISSING, int,             offset to first missing obs
    40 FIELD extendOffset,   OFF_EXTEND,  int,             offset to first extended measurement
     40FIELD refColor,       REF_COLOR,   float,           color of astrometry ref stars
    4141
    4242# 'flags' was called 'code' prior to 2009.02.07
  • trunk/Ohana/src/libautocode/def/image-ps1-v2.d

    r30604 r36833  
    4242FIELD     DECo,             DEC_CENTER,           float,          image center,              degrees
    4343FIELD     Radius,           RADIUS,               float,          image radius,              degrees
    44 FIELD     DUMMY,            DUMMY,                float,          dummy
     44FIELD     refColor,         REF_COLOR,            float,          dummy
    4545
    4646# should we define the max length of name as a macro?
  • trunk/Ohana/src/libautocode/def/image-ps1-v3.d

    r33647 r36833  
    4242FIELD     DECo,             DEC_CENTER,           float,          image center,              degrees
    4343FIELD     Radius,           RADIUS,               float,          image radius,              degrees
    44 FIELD     DUMMY,            DUMMY,                float,          dummy
     44FIELD     refColor,         REF_COLOR,            float,          dummy
    4545
    4646# should we define the max length of name as a macro?
  • trunk/Ohana/src/libautocode/def/image-ps1-v4.d

    r34260 r36833  
    4242FIELD     DECo,             DEC_CENTER,           float,          image center,              degrees
    4343FIELD     Radius,           RADIUS,               float,          image radius,              degrees
    44 FIELD     DUMMY,            DUMMY,                float,          dummy
     44FIELD     refColor,         REF_COLOR,            float,          dummy
    4545
    4646# should we define the max length of name as a macro?
  • trunk/Ohana/src/libautocode/def/image.d

    r33963 r36833  
    4444FIELD     DECo,             DEC_CENTER,           float,          image center,              degrees
    4545FIELD     Radius,           RADIUS,               float,          image radius,              degrees
    46 FIELD     DUMMY,            DUMMY,                float,          dummy
     46FIELD     refColor,         REF_COLOR,            float,          dummy
    4747
    4848FIELD     name,             NAME,                 char[121],      name of original image
  • trunk/Ohana/src/libdvo/include/dvo.h

    r36630 r36833  
    171171  ID_STAR_USE_PAR   = 0x00200000, // parallax used (not AVE or PM)
    172172  ID_STAR_NO_ASTROM = 0x00400000, // mean astrometry could not be measured
     173  ID_STAR_BAD_PM    = 0x00800000, // mean astrometry could not be measured
    173174  ID_OBJ_EXT        = 0x01000000, // extended in our data (eg, PS)
    174175  ID_OBJ_EXT_ALT    = 0x02000000, // extended in external data (eg, 2MASS)
     
    481482  uint32_t         measureOffset;        // offset to first psf measurement
    482483  uint32_t         missingOffset;        // offset to first missing obs
    483   uint32_t         extendOffset;         // offset to first extended measurement
     484  float            refColor;            // offset to first extended measurement
    484485  uint32_t         flags;                // average object flags (star; ghost; etc)
    485486  uint32_t         photFlagsUpper;       // upper bit of 2 bit summary of per-measure photflags
  • trunk/Ohana/src/libdvo/include/dvodb.h

    r36630 r36833  
    222222      AVE_EXTID_HI,
    223223      AVE_EXTID_LO,
     224      AVE_REF_COLOR,
    224225};
    225226
     
    289290      IMAGE_NFIT_ASTROM,
    290291      IMAGE_NLINK_PHOTOM,
    291       IMAGE_NLINK_ASTROM
     292      IMAGE_NLINK_ASTROM,
     293      IMAGE_REF_COLOR
    292294};
    293295
  • trunk/Ohana/src/libdvo/src/dbExtractAverages.c

    r36630 r36833  
    161161    case AVE_STARGAL:
    162162      value.Flt = average[0].stargal;
     163      break;
     164    case AVE_REF_COLOR:
     165      value.Flt = average[0].refColor;
    163166      break;
    164167
  • trunk/Ohana/src/libdvo/src/dbExtractImages.c

    r35102 r36833  
    6969      if (!haveCelestial) {
    7070        if (!FindMosaicForImage (image, Nimage, N)) return value;
    71         x = 0.5*image[N].NX;
    72         y = 0.5*image[N].NY;
     71        if (!strcmp(&image[N].coords.ctype[4], "-DIS")) {
     72          x = 0.0;
     73          y = 0.0;
     74        } else {
     75          x = 0.5*image[N].NX;
     76          y = 0.5*image[N].NY;
     77        }
    7378        XY_to_RD (&RAo, &DECo, x, y, &image[N].coords);
    7479        haveCelestial = TRUE;
     
    7984      if (!haveCelestial) {
    8085        if (!FindMosaicForImage (image, Nimage, N)) return value;
    81         x = 0.5*image[N].NX;
    82         y = 0.5*image[N].NY;
     86        if (!strcmp(&image[N].coords.ctype[4], "-DIS")) {
     87          x = 0.0;
     88          y = 0.0;
     89        } else {
     90          x = 0.5*image[N].NX;
     91          y = 0.5*image[N].NY;
     92        }
    8393        XY_to_RD (&RAo, &DECo, x, y, &image[N].coords);
    8494        haveCelestial = TRUE;
     
    90100        if (!haveCelestial) {
    91101          if (!FindMosaicForImage (image, Nimage, N)) return value;
    92           x = 0.5*image[N].NX;
    93           y = 0.5*image[N].NY;
     102          if (!strcmp(&image[N].coords.ctype[4], "-DIS")) {
     103            x = 0.0;
     104            y = 0.0;
     105          } else {
     106            x = 0.5*image[N].NX;
     107            y = 0.5*image[N].NY;
     108          }
    94109          XY_to_RD (&RAo, &DECo, x, y, &image[N].coords);
    95110          haveCelestial = TRUE;
     
    104119        if (!haveCelestial) {
    105120          if (!FindMosaicForImage (image, Nimage, N)) return value;
    106           x = 0.5*image[N].NX;
    107           y = 0.5*image[N].NY;
     121          if (!strcmp(&image[N].coords.ctype[4], "-DIS")) {
     122            x = 0.0;
     123            y = 0.0;
     124          } else {
     125            x = 0.5*image[N].NX;
     126            y = 0.5*image[N].NY;
     127          }
    108128          XY_to_RD (&RAo, &DECo, x, y, &image[N].coords);
    109129          haveCelestial = TRUE;
     
    118138        if (!haveCelestial) {
    119139          if (!FindMosaicForImage (image, Nimage, N)) return value;
    120           x = 0.5*image[N].NX;
    121           y = 0.5*image[N].NY;
     140          if (!strcmp(&image[N].coords.ctype[4], "-DIS")) {
     141            x = 0.0;
     142            y = 0.0;
     143          } else {
     144            x = 0.5*image[N].NX;
     145            y = 0.5*image[N].NY;
     146          }
    122147          XY_to_RD (&RAo, &DECo, x, y, &image[N].coords);
    123148          haveCelestial = TRUE;
     
    132157        if (!haveCelestial) {
    133158          if (!FindMosaicForImage (image, Nimage, N)) return value;
    134           x = 0.5*image[N].NX;
    135           y = 0.5*image[N].NY;
     159          if (!strcmp(&image[N].coords.ctype[4], "-DIS")) {
     160            x = 0.0;
     161            y = 0.0;
     162          } else {
     163            x = 0.5*image[N].NX;
     164            y = 0.5*image[N].NY;
     165          }
    136166          XY_to_RD (&RAo, &DECo, x, y, &image[N].coords);
    137167          haveCelestial = TRUE;
     
    344374      value.Int = image[N].nLinkAstrom;
    345375      break;
     376
     377    case IMAGE_REF_COLOR:
     378      value.Flt = image[N].refColor;
     379      break;
    346380  }
    347381  return (value);
  • trunk/Ohana/src/libdvo/src/dbFields.c

    r36630 r36833  
    381381  if (!strcasecmp (fieldName, "EXTID_HI"))    ESCAPE (AVE_EXTID_HI,    MAG_NONE, OPIHI_INT);
    382382  if (!strcasecmp (fieldName, "EXTID_LO"))    ESCAPE (AVE_EXTID_LO,    MAG_NONE, OPIHI_INT);
     383  if (!strcasecmp (fieldName, "REF_COLOR"))   ESCAPE (AVE_REF_COLOR,   MAG_NONE, OPIHI_FLT);
    383384
    384385  // check for code:mode in photcode name
     
    522523  if (!strcasecmp (fieldName, "UCDIST")) ESCAPE (IMAGE_UBERCAL_DIST,MAG_NONE, OPIHI_INT);
    523524
    524   if (!strcasecmp (fieldName, "NFIT_PHOTOM")) ESCAPE (IMAGE_NFIT_PHOTOM,MAG_NONE, OPIHI_INT);
    525   if (!strcasecmp (fieldName, "NFIT_ASTROM")) ESCAPE (IMAGE_NFIT_ASTROM,MAG_NONE, OPIHI_INT);
    526   if (!strcasecmp (fieldName, "NLINK_PHOTOM")) ESCAPE (IMAGE_NLINK_PHOTOM,MAG_NONE, OPIHI_INT);
    527   if (!strcasecmp (fieldName, "NLINK_ASTROM")) ESCAPE (IMAGE_NLINK_ASTROM,MAG_NONE, OPIHI_INT);
     525  if (!strcasecmp (fieldName, "NFIT_PHOTOM"))  ESCAPE (IMAGE_NFIT_PHOTOM,  MAG_NONE, OPIHI_INT);
     526  if (!strcasecmp (fieldName, "NFIT_ASTROM"))  ESCAPE (IMAGE_NFIT_ASTROM,  MAG_NONE, OPIHI_INT);
     527  if (!strcasecmp (fieldName, "NLINK_PHOTOM")) ESCAPE (IMAGE_NLINK_PHOTOM, MAG_NONE, OPIHI_INT);
     528  if (!strcasecmp (fieldName, "NLINK_ASTROM")) ESCAPE (IMAGE_NLINK_ASTROM, MAG_NONE, OPIHI_INT);
     529  if (!strcasecmp (fieldName, "REF_COLOR"))    ESCAPE (IMAGE_REF_COLOR,    MAG_NONE, OPIHI_FLT);
    528530
    529531  // for words that don't parse, try a photcode
  • trunk/Ohana/src/libdvo/src/dvo_catalog.c

    r36630 r36833  
    100100  average->measureOffset   = -1;
    101101  average->missingOffset   = -1;
    102   average->extendOffset    = -1;
     102  average->refColor        = NAN;
    103103
    104104  average->flags           = 0;
  • trunk/Ohana/src/libdvo/src/dvo_convert_PS1_DEV_2.c

    r35162 r36833  
    119119    out[i].measureOffset = in[i].measureOffset;
    120120    out[i].missingOffset = in[i].missingOffset;
    121     out[i].extendOffset  = in[i].extendOffset;
     121    out[i].refColor  = in[i].refColor;
    122122    out[i].objID         = in[i].objID;
    123123    out[i].catID         = in[i].catID;
     
    153153    out[i].measureOffset = in[i].measureOffset;
    154154    out[i].missingOffset = in[i].missingOffset;
    155     out[i].extendOffset  = in[i].extendOffset;
     155    out[i].refColor  = in[i].refColor;
    156156    out[i].objID         = in[i].objID;
    157157    out[i].catID         = in[i].catID;
  • trunk/Ohana/src/libdvo/src/dvo_convert_PS1_V1.c

    r35162 r36833  
    140140    out[i].measureOffset = in[i].measureOffset;
    141141    out[i].missingOffset = in[i].missingOffset;
    142     out[i].extendOffset  = in[i].extendOffset;
     142    out[i].refColor  = in[i].refColor;
    143143    out[i].flags         = in[i].flags;   
    144144    out[i].objID         = in[i].objID;
     
    175175    out[i].measureOffset = in[i].measureOffset;
    176176    out[i].missingOffset = in[i].missingOffset;
    177     out[i].extendOffset  = in[i].extendOffset;
     177    out[i].refColor  = in[i].refColor;
    178178    out[i].flags         = in[i].flags;   
    179179    out[i].objID         = in[i].objID;
  • trunk/Ohana/src/libdvo/src/dvo_convert_PS1_V2.c

    r35162 r36833  
    145145    out[i].measureOffset = in[i].measureOffset;
    146146    out[i].missingOffset = in[i].missingOffset;
    147     out[i].extendOffset  = in[i].extendOffset;
     147    out[i].refColor  = in[i].refColor;
    148148    out[i].flags         = in[i].flags;   
    149149    out[i].objID         = in[i].objID;
     
    184184    out[i].measureOffset = in[i].measureOffset;
    185185    out[i].missingOffset = in[i].missingOffset;
    186     out[i].extendOffset  = in[i].extendOffset;
     186    out[i].refColor  = in[i].refColor;
    187187    out[i].flags         = in[i].flags;   
    188188    out[i].objID         = in[i].objID;
     
    269269    out[i].DECo             = in[i].DECo;
    270270    out[i].Radius           = in[i].Radius;
     271    out[i].refColor         = in[i].refColor;
    271272
    272273    out[i].detection_limit  = in[i].detection_limit;
     
    331332    out[i].DECo             = in[i].DECo;
    332333    out[i].Radius           = in[i].Radius;
     334    out[i].refColor         = in[i].refColor;
    333335
    334336    out[i].detection_limit  = in[i].detection_limit;
  • trunk/Ohana/src/libdvo/src/dvo_convert_PS1_V3.c

    r35162 r36833  
    145145    out[i].measureOffset = in[i].measureOffset;
    146146    out[i].missingOffset = in[i].missingOffset;
    147     out[i].extendOffset  = in[i].extendOffset;
     147    out[i].refColor  = in[i].refColor;
    148148    out[i].flags         = in[i].flags;   
    149149    out[i].photFlagsUpper = in[i].photFlagsUpper;   
     
    186186    out[i].measureOffset = in[i].measureOffset;
    187187    out[i].missingOffset = in[i].missingOffset;
    188     out[i].extendOffset  = in[i].extendOffset;
     188    out[i].refColor  = in[i].refColor;
    189189    out[i].flags         = in[i].flags;   
    190190    out[i].photFlagsUpper = in[i].photFlagsUpper;   
     
    279279    out[i].DECo             = in[i].DECo;
    280280    out[i].Radius           = in[i].Radius;
     281    out[i].refColor         = in[i].refColor;
    281282
    282283    out[i].detection_limit  = in[i].detection_limit;
     
    341342    out[i].DECo             = in[i].DECo;
    342343    out[i].Radius           = in[i].Radius;
     344    out[i].refColor         = in[i].refColor;
    343345
    344346    out[i].detection_limit  = in[i].detection_limit;
  • trunk/Ohana/src/libdvo/src/dvo_convert_PS1_V4.c

    r36630 r36833  
    166166    out[i].measureOffset = in[i].measureOffset;
    167167    out[i].missingOffset = in[i].missingOffset;
    168     out[i].extendOffset  = in[i].extendOffset;
     168    out[i].refColor      = in[i].refColor;
    169169    out[i].flags         = in[i].flags;   
    170170    out[i].photFlagsUpper = in[i].photFlagsUpper;   
     
    210210    out[i].measureOffset = in[i].measureOffset;
    211211    out[i].missingOffset = in[i].missingOffset;
    212     out[i].extendOffset  = in[i].extendOffset;
     212    out[i].refColor      = in[i].refColor;
    213213    out[i].flags         = in[i].flags;   
    214214    out[i].photFlagsUpper = in[i].photFlagsUpper;   
     
    317317    out[i].DECo             = in[i].DECo;
    318318    out[i].Radius           = in[i].Radius;
     319    out[i].refColor         = in[i].refColor;
    319320
    320321    out[i].detection_limit  = in[i].detection_limit;
     
    379380    out[i].DECo             = in[i].DECo;
    380381    out[i].Radius           = in[i].Radius;
     382    out[i].refColor         = in[i].refColor;
    381383
    382384    out[i].detection_limit  = in[i].detection_limit;
     
    680682    out[i].measureOffset = in[i].measureOffset;
    681683    out[i].missingOffset = in[i].missingOffset;
    682     out[i].extendOffset  = in[i].extendOffset;
     684    out[i].refColor      = in[i].refColor;
    683685    out[i].flags         = in[i].flags;   
    684686    out[i].photFlagsUpper = in[i].photFlagsUpper;   
  • trunk/Ohana/src/libdvo/src/skyregion_io.c

    r33649 r36833  
    7878  gfits_free_matrix (&matrix);
    7979  gfits_free_header (&theader);
     80  fclose (f);
    8081
    8182  return (skytable);
  • trunk/Ohana/src/libkapa/src/KapaOpen.c

    r35761 r36833  
    218218  // apparently, I can connect on someone else's port (eg GoogleTalkPlugin)
    219219  // do a simple handshake before we set !NONBLOCK:
    220   char line[5];
    221   int Nout = read (InitSocket, line, 4);
    222   if (Nout != 4) {
    223     if (DEBUG) fprintf (stderr, "connection failed\n");
    224     close (InitSocket);
    225     return (-1);
    226   }
    227   if (strncmp (line, "KAPA", 4)) {
    228     if (DEBUG) fprintf (stderr, "connection to the wrong server\n");
    229     close (InitSocket);
    230     return (-1);
     220
     221  // ensure the socket is NONBLOCK first...
     222  fcntl (InitSocket, F_SETFL, O_NONBLOCK);
     223
     224  int Ntry = 0;
     225retry_message:
     226  {
     227    char line[5];
     228    int Nout = read (InitSocket, line, 4);
     229    if ((Nout == -1) && (errno == EAGAIN)) {
     230      Ntry ++;
     231      if (Ntry > 500) {
     232        if (DEBUG) fprintf (stderr, "handshake failure\n");
     233        close (InitSocket);
     234        return (-1);
     235      }
     236      usleep (10000);
     237      goto retry_message;
     238    }
     239    if (Nout != 4) {
     240      if (DEBUG) fprintf (stderr, "connection failed\n");
     241      close (InitSocket);
     242      return (-1);
     243    }
     244    if (strncmp (line, "KAPA", 4)) {
     245      if (DEBUG) fprintf (stderr, "connection to the wrong server\n");
     246      close (InitSocket);
     247      return (-1);
     248    }
    231249  }
    232250
     
    265283  }
    266284
    267 # define NTRY 500
     285# define NTRY 1000
    268286  Ntry = 0;
    269287  while (Ntry < NTRY) {
    270288    sock = KapaClientSocket ("localhost");
    271289    if (sock != -1) break;
    272     // if (errno != EAGAIN) break;
    273     if (errno != ECONNREFUSED) break;
     290    if (errno != ECONNREFUSED) {
     291      perror ("KapaOpen");
     292      break;
     293    }
    274294    // no connection yet. try again, but first check
    275295    // if the kapa job has exited
     
    330350  fcntl (InitSocket, F_SETFL, O_NONBLOCK);
    331351
    332 # define NTRY 500
     352# define NTRY 1000
    333353  Ntry = 0;
    334354  while (Ntry < NTRY) {
  • trunk/Ohana/src/libkapa/src/KiiPicture.c

    r34088 r36833  
    2121int KiiNewPicture1D (int fd, KiiImage *image, KapaImageData *data, Coords *coords) {
    2222
    23   int Nwrite, Npix, Ncolors, size;
     23  int Ncolors;
     24  off_t Nwrite, Npix, size;
    2425  float min, max;
    2526
     
    3637  KiiSendMessage (fd, "%8d %8d", image[0].Nx, image[0].Ny);
    3738  KiiSendMessage (fd, "%f %f %s %s", data[0].zero, data[0].range, data[0].name, data[0].file);
    38   KiiSendMessage (fd, "%f %f %d ", min, max, size);
     39  KiiSendMessage (fd, "%f %f "OFF_T_FMT" ", min, max, size);
    3940  KiiSendMessage (fd, "%f %f %g %g %g ", coords[0].crval1, coords[0].crpix1, coords[0].cdelt1, coords[0].pc1_1, coords[0].pc1_2);
    4041  KiiSendMessage (fd, "%f %f %g %g %g ", coords[0].crval2, coords[0].crpix2, coords[0].cdelt2, coords[0].pc2_1, coords[0].pc2_2);
     
    4243
    4344  /* send the image data */
    44   Nwrite = write (fd, image[0].data1d, size);
    45   if (Nwrite != size) {
    46     fprintf (stderr, "error sending picture to kapa\n");
    47     return (FALSE);
     45  off_t bytes_left = size;
     46  while (bytes_left > 0) {
     47    Nwrite = write (fd, image[0].data1d, bytes_left);
     48    if (Nwrite == 0) {
     49      fprintf (stderr, "unable to send more data to kapa?\n");
     50      return FALSE;
     51    }
     52    if (Nwrite < 0) {
     53      if (errno == EAGAIN) continue;
     54      if (errno == EWOULDBLOCK) continue;
     55      perror ("KiiNewPicture1D:");
     56      return FALSE;
     57    }
     58    bytes_left -= Nwrite;
    4859  }
    4960
     
    114125int KiiNewPicture2D (int fd, KiiImage *image, KapaImageData *data, Coords *coords) {
    115126
    116   int j;
    117   int Nwrite, Nbytes, Npix, Ncolors, size;
     127  int Ncolors;
     128  off_t j, Nwrite, Nbytes, Npix, size;
    118129  float min, max;
    119130
     
    130141  KiiSendMessage (fd, "%8d %8d", image[0].Nx, image[0].Ny);
    131142  KiiSendMessage (fd, "%f %f %s %s", data[0].zero, data[0].range, data[0].name, data[0].file);
    132   KiiSendMessage (fd, "%f %f %d ", min, max, size);
     143  KiiSendMessage (fd, "%f %f "OFF_T_FMT" ", min, max, size);
    133144  KiiSendMessage (fd, "%f %f %g %g %g ", coords[0].crval1, coords[0].crpix1, coords[0].cdelt1, coords[0].pc1_1, coords[0].pc1_2);
    134145  KiiSendMessage (fd, "%f %f %g %g %g ", coords[0].crval2, coords[0].crpix2, coords[0].cdelt2, coords[0].pc2_1, coords[0].pc2_2);
  • trunk/Ohana/src/libohana/src

    • Property svn:mergeinfo deleted
  • trunk/Ohana/src/opihi/dvo/Makefile

    r36375 r36833  
    4141$(SRC)/cmpread.$(ARCH).o                \
    4242$(SRC)/coordimage.$(ARCH).o             \
     43$(SRC)/coordmosaic.$(ARCH).o            \
    4344$(SRC)/ddmags.$(ARCH).o         \
    4445$(SRC)/detrend.$(ARCH).o                \
  • trunk/Ohana/src/opihi/dvo/init.c

    r36375 r36833  
    1414int cmpread         PROTO((int, char **));
    1515int coordimage      PROTO((int, char **));
     16int coordmosaic     PROTO((int, char **));
     17int psastro_model   PROTO((int, char **));
    1618int ddmags          PROTO((int, char **));
    1719int detrend         PROTO((int, char **));
     
    7375  {1, "cmpread",     cmpread,      "read data from cmp format files"},
    7476  {1, "coordimage",  coordimage,   "generate a map of the transformation residuals"},
     77  {1, "coordmosaic", coordmosaic,  "generate a map of the distortion"},
     78  {1, "psastro_model", psastro_model, "save psastro-format astrometry model"},
    7579  {1, "ddmags",      ddmags,       "plot magnitude differences"},
    7680  {1, "detrend",     detrend,      "extract from detrend database?"},
  • trunk/Ohana/src/photdbc/Makefile

    r33655 r36833  
    3131$(SRC)/args.$(ARCH).o              \
    3232$(SRC)/copy_images.$(ARCH).o       \
     33$(SRC)/select_images.$(ARCH).o     \
    3334$(SRC)/Shutdown.$(ARCH).o          \
    3435$(SRC)/join_stars.$(ARCH).o        \
     
    4142$(SRC)/ConfigInit.$(ARCH).o        \
    4243$(SRC)/args.$(ARCH).o              \
    43 $(SRC)/copy_images.$(ARCH).o       \
    4444$(SRC)/Shutdown.$(ARCH).o          \
    4545$(SRC)/join_stars.$(ARCH).o        \
  • trunk/Ohana/src/photdbc/include/photdbc.h

    r33963 r36833  
    2020  int    Nmeas;
    2121} StatType;
     22
     23typedef struct {
     24  double Xc[5];
     25  double Yc[5];
     26  double Rc;
     27  double Dc;
     28} SkyRegionCoords;
     29
     30# define MARKTIME(MSG,...) {                    \
     31    gettimeofday (&stopTimer, (void *) NULL);   \
     32    float dtime = DTIME (stopTimer, startTimer);        \
     33    fprintf (stderr, MSG, __VA_ARGS__); }
     34
     35# define INITTIME \
     36  struct timeval startTimer, stopTimer; \
     37  gettimeofday (&startTimer, (void *) NULL);
    2238
    2339int    PARALLEL;
     
    8096PhotCodeData photcodes;
    8197
    82 char          *PHOTCODE_DROP_LIST, *PHOTCODE_SKIP_LIST;
    83 int           NphotcodesDrop,      NphotcodesSkip;
    84 PhotCode     **photcodesDrop,     **photcodesSkip;
     98char          *PHOTCODE_DROP_LIST, *PHOTCODE_KEEP_LIST;
     99int           NphotcodesDrop,      NphotcodesKeep;
     100PhotCode     **photcodesDrop,     **photcodesKeep;
    85101
    86102# define FLAG_AREA            0X0001
     
    131147void SetProtect (int mode);
    132148int SetSignals (void);
    133 int copy_images (char *outdir);
     149int copy_images (char *outdir, SkyList *skylist);
    134150void usage();
     151
     152void dsortindex (double *X, off_t *Y, int N);
     153off_t getRegionStartByRA (double R, double *Rref, off_t Nregions);
     154
     155Image *select_images (SkyList *skylist, Image *timage, off_t Ntimage, off_t **LineNumber, off_t *Nimage, int UseFullOverlap);
    135156
    136157int photdbc_catalogs (char *outroot, SkyList *skylist, int hostID);
  • trunk/Ohana/src/photdbc/src/args.c

    r33655 r36833  
    8585  }
    8686
     87  // measurements with these photcodes are not copied to the output
    8788  PHOTCODE_DROP_LIST = NULL;
    8889  if ((N = get_argument (argc, argv, "-photcode-drop"))) {
     
    9293  }
    9394
    94   PHOTCODE_SKIP_LIST = NULL;
    95   if ((N = get_argument (argc, argv, "-photcode-skip"))) {
    96     remove_argument (N, &argc, argv);
    97     PHOTCODE_SKIP_LIST = strcreate(argv[N]);
     95  // measurements with these photcodes are kept ***regardless of quality***
     96  // -photcode-keep J will keep all J-band measurements of all kinds
     97  // -photcode-keep GPC1.02.g will keep all g-band measurements from chip XY02
     98  PHOTCODE_KEEP_LIST = NULL;
     99  if ((N = get_argument (argc, argv, "-photcode-keep"))) {
     100    remove_argument (N, &argc, argv);
     101    PHOTCODE_KEEP_LIST = strcreate(argv[N]);
    98102    remove_argument (N, &argc, argv);
    99103  }
     
    254258  }
    255259
    256   PHOTCODE_SKIP_LIST = NULL;
    257   if ((N = get_argument (argc, argv, "-photcode-skip"))) {
    258     remove_argument (N, &argc, argv);
    259     PHOTCODE_SKIP_LIST = strcreate(argv[N]);
     260  PHOTCODE_KEEP_LIST = NULL;
     261  if ((N = get_argument (argc, argv, "-photcode-keep"))) {
     262    remove_argument (N, &argc, argv);
     263    PHOTCODE_KEEP_LIST = strcreate(argv[N]);
    260264    remove_argument (N, &argc, argv);
    261265  }
     
    290294 
    291295  fprintf (stderr, " -photcode-drop   : remove these photcodes from the output (REF or DEP only)\n");
    292   fprintf (stderr, " -photcode-skip  : ignore these photcodes when assessing the validity (keep unless object is dropped)\n");
     296  fprintf (stderr, " -photcode-keep  : ignore these photcodes when assessing the validity (keep unless object is dropped)\n");
    293297
    294298  fprintf (stderr, " -instmag (min) (max) : range of valid instrumental magnitudes (or measurements are dropped)\n");
  • trunk/Ohana/src/photdbc/src/copy_images.c

    r33655 r36833  
    11# include "photdbc.h"
    22
    3 int copy_images (char *outdir) {
     3int copy_images (char *outdir, SkyList *skylist) {
    44
    55  int status;
    6   off_t Nimage;
     6  off_t Nimage, Nsubset;
     7  off_t *LineNumber;
    78  char *ImageOut;
    89  unsigned int imageID;
    910  FITS_DB in;
    1011  FITS_DB out;
    11   Image *image;
     12  Image *image, *subset;
    1213  struct stat filestat;
    1314  char *path;
     
    5051      exit (2);
    5152  }
    52   dvo_image_addrows (&out, image, Nimage);
     53
     54  subset = select_images (skylist, image, Nimage, &LineNumber, &Nsubset, FALSE);
     55
     56  dvo_image_addrows (&out, subset, Nsubset);
    5357
    5458  // note that imageID is unsigned int
  • trunk/Ohana/src/photdbc/src/initialize.c

    r33655 r36833  
    1313
    1414  photcodesDrop = ParsePhotcodeList (PHOTCODE_DROP_LIST, &NphotcodesDrop, FALSE);
    15   photcodesSkip = ParsePhotcodeList (PHOTCODE_SKIP_LIST, &NphotcodesSkip, FALSE);
     15  photcodesKeep = ParsePhotcodeList (PHOTCODE_KEEP_LIST, &NphotcodesKeep, FALSE);
    1616
    1717  if (SHOW_PARAMS) {
     
    5454
    5555  photcodesDrop = ParsePhotcodeList (PHOTCODE_DROP_LIST, &NphotcodesDrop, FALSE);
    56   photcodesSkip = ParsePhotcodeList (PHOTCODE_SKIP_LIST, &NphotcodesSkip, FALSE);
     56  photcodesKeep = ParsePhotcodeList (PHOTCODE_KEEP_LIST, &NphotcodesKeep, FALSE);
    5757
    5858  if (SHOW_PARAMS) {
  • trunk/Ohana/src/photdbc/src/make_subcatalog.c

    r31635 r36833  
    99  off_t NAVERAGE, NMEASURE, Naverage, Nmeasure, Nm, Nsecfilt;
    1010  double mag, minMag, minSigma;
    11   int keep, *secSkip;
     11  int keep, *secKeep;
    1212  PhotCode *photcode;
    1313 
     
    1616
    1717  // set up a list of SEC entries to ignore when evaluating a source
    18   ALLOCATE (secSkip, int, Nsecfilt);
     18  ALLOCATE (secKeep, int, Nsecfilt);
    1919  for (i = 0; i < Nsecfilt; i++) {
    20       secSkip[i] = FALSE;
     20      secKeep[i] = FALSE;
    2121      photcode = GetPhotcodebyNsec(i);
    22       for (k = 0; k < NphotcodesSkip; k++) {
    23           if (photcodesSkip[k][0].code != photcode[0].code) continue;
    24           secSkip[i] = TRUE;
     22      for (k = 0; k < NphotcodesKeep; k++) {
     23          if (photcodesKeep[k][0].code != photcode[0].code) continue;
     24          secKeep[i] = TRUE;
    2525      }
    2626  }
     
    6767      keep = FALSE;
    6868      for (j = 0; !keep && (j < Nsecfilt); j++) {
    69           if (secSkip[j]) continue;
     69          if (secKeep[j]) continue;
    7070          if (catalog[0].secfilt[Nsecfilt*i+j].Ncode >= NCODE_MIN) {
    7171              keep = TRUE;
     
    101101      if (NphotcodesDrop > 0) {
    102102          found = FALSE;
    103           for (k = 0; (k < NphotcodesSkip) && !found; k++) {
    104               if (photcodesSkip[k][0].code == catalog[0].measure[offset].photcode) found = TRUE;
    105               if (photcodesSkip[k][0].code == GetPhotcodeEquivCodebyCode(catalog[0].measure[offset].photcode)) found = TRUE;
     103          for (k = 0; (k < NphotcodesDrop) && !found; k++) {
     104              if (photcodesDrop[k][0].code == catalog[0].measure[offset].photcode) found = TRUE;
     105              if (photcodesDrop[k][0].code == GetPhotcodeEquivCodebyCode(catalog[0].measure[offset].photcode)) found = TRUE;
    106106          }
    107107          if (found) continue;
     
    109109 
    110110      // ignore certain photcodes to assess the measurements
    111       if (NphotcodesSkip > 0) {
     111      if (NphotcodesKeep > 0) {
    112112          found = FALSE;
    113           for (k = 0; (k < NphotcodesSkip) && !found; k++) {
    114               if (photcodesSkip[k][0].code == catalog[0].measure[offset].photcode) found = TRUE;
    115               if (photcodesSkip[k][0].code == GetPhotcodeEquivCodebyCode(catalog[0].measure[offset].photcode)) found = TRUE;
     113          for (k = 0; (k < NphotcodesKeep) && !found; k++) {
     114              if (photcodesKeep[k][0].code == catalog[0].measure[offset].photcode) found = TRUE;
     115              if (photcodesKeep[k][0].code == GetPhotcodeEquivCodebyCode(catalog[0].measure[offset].photcode)) found = TRUE;
    116116          }
    117117          if (found) goto keep;
  • trunk/Ohana/src/photdbc/src/photdbc.c

    r33655 r36833  
    99  initialize (argc, argv);
    1010
    11   // load and copy the image table
    12   copy_images (argv[1]);
    13 
    1411  // the output catalog needs to inherit the SKY_DEPTH of the input catalog
    1512  sky = SkyTableLoadOptimal (CATDIR, NULL, NULL, TRUE, 0, VERBOSE);
    1613  SkyTableSetFilenames (sky, CATDIR, "cpt");
    1714  skylist = SkyListByPatch (sky, -1, &REGION);
     15
     16  // load and copy the image table
     17  copy_images (argv[1], skylist);
    1818
    1919  // hostID is 0 for master program
  • trunk/Ohana/src/photdbc/src/photdbc_catalogs.c

    r33963 r36833  
    120120    if (ExcludeByMaxMinMag) { snprintf (tmpline, DVO_MAX_PATH, "%s -maxminmag %f",     command, MAX_MIN_MAG);                   strcpy (command, tmpline); }
    121121    if (PHOTCODE_DROP_LIST) { snprintf (tmpline, DVO_MAX_PATH, "%s -photcode-drop %s", command, PHOTCODE_DROP_LIST);    strcpy (command, tmpline); }
    122     if (PHOTCODE_SKIP_LIST) { snprintf (tmpline, DVO_MAX_PATH, "%s -photcode-skip %s", command, PHOTCODE_SKIP_LIST);    strcpy (command, tmpline); }
     122    if (PHOTCODE_KEEP_LIST) { snprintf (tmpline, DVO_MAX_PATH, "%s -photcode-keep %s", command, PHOTCODE_KEEP_LIST);    strcpy (command, tmpline); }
    123123    if (CATFORMAT)          { snprintf (tmpline, DVO_MAX_PATH, "%s -set-format %s",    command, CATFORMAT);                     strcpy (command, tmpline); }
    124124    if (CATMODE)            { snprintf (tmpline, DVO_MAX_PATH, "%s -set-mode %s",      command, CATMODE);                       strcpy (command, tmpline); }
  • trunk/Ohana/src/relastro/Makefile

    r36630 r36833  
    127127$(SRC)/plotstuff.$(ARCH).o           \
    128128$(SRC)/relastroVisual.$(ARCH).o \
     129$(SRC)/syncfile.$(ARCH).o \
    129130$(SRC)/BrightCatalog.$(ARCH).o
    130131
  • trunk/Ohana/src/relastro/include/relastro.h

    r36630 r36833  
    8383  double L, M;  /* Focal Plane   - pixels  */
    8484  double X, Y;  /* Chip Coords   - pixels  */
    85   double Mag;
    86   double dMag;
    87   double dPos;
     85  float Mag;
     86  float Color;
     87  float dMag;
     88  float dPos;
    8889  int mask;
    8990  int Nmeas;
     
    207208double MaxDensityValue;
    208209
    209 char          *PHOTCODE_KEEP_LIST, *PHOTCODE_SKIP_LIST;
    210 int           NphotcodesKeep,      NphotcodesSkip;
    211 PhotCode     **photcodesKeep,     **photcodesSkip;
     210char          *PHOTCODE_KEEP_LIST, *PHOTCODE_SKIP_LIST, *PHOTCODE_RESET_LIST;
     211int           NphotcodesKeep,      NphotcodesSkip,      NphotcodesReset;
     212PhotCode     **photcodesKeep,     **photcodesSkip,     **photcodesReset;
    212213
    213214char          *PHOTCODE_A_LIST,  *PHOTCODE_B_LIST;
     
    217218SkyRegionSelection SELECTION;
    218219
     220char         *DCR_COLOR_POS,    *DCR_COLOR_NEG;
     221PhotCode     *DCR_PHOTCODE_POS, *DCR_PHOTCODE_NEG;
     222int           DCR_NSEC_POS,      DCR_NSEC_NEG;
     223
    219224int ImagSelect;
    220225double ImagMin, ImagMax;
     
    233238int    CLIP_THRESH;
    234239int USE_BASIC_CHECK;
     240
     241int ExcludeBogus;
     242double ExcludeBogusRadius;
    235243
    236244FitMode FIT_MODE;
     
    307315int           liststats           PROTO((double *value, double *dvalue, int N, StatType *stats));
    308316int           liststats_pos       PROTO((double *value, double *dvalue, int N, StatType *stats, int XVERB));
    309 Catalog      *load_catalogs       PROTO((SkyList *skylist, int *Ncatalog, int subselect, int hostID, char *hostpath));
     317Catalog      *load_catalogs       PROTO((SkyList *skylist, int *Ncatalog, int subselect, int hostID, char *hostpath, char *syncfile));
    310318int           load_images         PROTO((FITS_DB *db, SkyList *skylist, int UseFullOverlap));
    311319Image        *select_images       PROTO((SkyList *skylist, Image *timage, off_t Ntimage, off_t **LineNumber, off_t *Nimage, int UseFullOverlap));
     
    417425int setMeanR (double ra_fit, MeasureTiny *measure, Average *average, SecFilt *secfilt);
    418426int setMeanD (double dec_fit, MeasureTiny *measure, Average *average, SecFilt *secfilt);
     427double getMeanR_Big (Measure *measure, Average *average, SecFilt *secfilt);
     428double getMeanD_Big (Measure *measure, Average *average, SecFilt *secfilt);
    419429int setMeanR_Big (double ra_fit, Measure *measure, Average *average, SecFilt *secfilt);
    420430int setMeanD_Big (double dec_fit, Measure *measure, Average *average, SecFilt *secfilt);
     
    521531int ImageTableSave (char *filename, Image *images, off_t Nimages);
    522532int select_mosaics_hostregion (RegionHostTable *regionHosts, Image *image, off_t Nimage);
     533
     534float getColor (off_t meas, int cat);
     535
     536int strextend (char *input, char *format,...);
     537
     538int areImagesLoaded ();
     539int areImagesMatched ();
  • trunk/Ohana/src/relastro/src/BrightCatalog.c

    r35105 r36833  
    167167    GET_COLUMN(measureOffset,  "OFF_MEASURE", int);
    168168    GET_COLUMN(missingOffset,  "OFF_MISSING", int);
    169     GET_COLUMN(extendOffset,   "OFF_EXTEND",  int);
     169    GET_COLUMN(refColor,   "REF_COLOR",  float);
    170170    GET_COLUMN(flags,          "FLAGS",       int);
    171171    GET_COLUMN(photFlagsUpper, "PHOTFLAGS_U", int);
     
    202202      average[i].measureOffset   = measureOffset[i]   ;
    203203      average[i].missingOffset   = missingOffset[i]   ;
    204       average[i].extendOffset    = extendOffset[i]    ; 
     204      average[i].refColor    = refColor[i]    ; 
    205205      average[i].flags           = flags[i]           ;
    206206      average[i].photFlagsUpper  = photFlagsUpper[i]  ;
     
    234234    free (measureOffset);
    235235    free (missingOffset);
    236     free (extendOffset);
     236    free (refColor);
    237237    free (flags);
    238238    free (photFlagsUpper);
     
    476476    gfits_define_bintable_column (&theader, "J", "OFF_MEASURE", "offset to first psf measurement",                        "", 1.0, 1.0*0x8000);
    477477    gfits_define_bintable_column (&theader, "J", "OFF_MISSING", "offset to first missing obs",                            "", 1.0, 1.0*0x8000);
    478     gfits_define_bintable_column (&theader, "J", "OFF_EXTEND",  "offset to first extended measurement",                   "", 1.0, 1.0*0x8000);
     478    gfits_define_bintable_column (&theader, "E", "REF_COLOR",   "reference color",                                        "", 1.0, 0.0);
    479479    gfits_define_bintable_column (&theader, "J", "FLAGS",       "average object flags (star; ghost; etc)",                "", 1.0, 1.0*0x8000);
    480480    gfits_define_bintable_column (&theader, "J", "PHOTFLAGS_U", "upper bit of 2 bit summary of per-measure photflags",    "", 1.0, 1.0*0x8000);
     
    510510    int      *measureOffset ; ALLOCATE (measureOffset , int     , catalog->Naverage);
    511511    int      *missingOffset ; ALLOCATE (missingOffset , int     , catalog->Naverage);
    512     int      *extendOffset  ; ALLOCATE (extendOffset  , int     , catalog->Naverage);
     512    float    *refColor      ; ALLOCATE (refColor  , float     , catalog->Naverage);
    513513    int      *flags         ; ALLOCATE (flags         , int     , catalog->Naverage);
    514514    int      *photFlagsUpper; ALLOCATE (photFlagsUpper, int     , catalog->Naverage);
     
    543543      measureOffset[i]   = average[i].measureOffset   ;
    544544      missingOffset[i]   = average[i].missingOffset   ;
    545       extendOffset[i]    = average[i].extendOffset    ; 
     545      refColor[i]    = average[i].refColor    ; 
    546546      flags[i]           = average[i].flags           ;
    547547      photFlagsUpper[i]  = average[i].photFlagsUpper  ;
     
    575575    gfits_set_bintable_column (&theader, &ftable, "OFF_MEASURE", measureOffset,   catalog->Naverage);
    576576    gfits_set_bintable_column (&theader, &ftable, "OFF_MISSING", missingOffset,   catalog->Naverage);
    577     gfits_set_bintable_column (&theader, &ftable, "OFF_EXTEND",  extendOffset,    catalog->Naverage);
     577    gfits_set_bintable_column (&theader, &ftable, "REF_COLOR",   refColor,        catalog->Naverage);
    578578    gfits_set_bintable_column (&theader, &ftable, "FLAGS",       flags,           catalog->Naverage);
    579579    gfits_set_bintable_column (&theader, &ftable, "PHOTFLAGS_U", photFlagsUpper,  catalog->Naverage);
     
    605605    free (measureOffset);
    606606    free (missingOffset);
    607     free (extendOffset);
     607    free (refColor);
    608608    free (flags);
    609609    free (photFlagsUpper);
  • trunk/Ohana/src/relastro/src/FitChip.c

    r33652 r36833  
    44int FitChip (StarData *raw, StarData *ref, int Nmatch, Image *image) {
    55
    6   int i, NstatFull, Nstat, Niter;
     6  int i, NstatFull, Nstat, Niter, Nkeep;
    77  float dLsig, dMsig, dRsig;
    88  float dLsigFull, dMsigFull, dRsigFull;
     
    2424
    2525    // measure the scatter distribution (use only the bright end detections)
     26    Nkeep = 0;
    2627    for (i = 0; i < Nmatch; i++) {
    2728      if (raw[i].mask) continue;
     
    4243        continue;
    4344      }
    44     }
     45      Nkeep ++;
     46    }
     47
     48    // I'm rejecting some points from the fit above; I count the remainders and actually
     49    // use that count to set order_use below
    4550
    4651    // figures to assess the fitting process:
     
    5055    relastroVisualPlotChipFit(raw, ref, dRmax, Nmatch);
    5156
    52     // fit the requested order polynomial
    53     if (CHIPORDER > 0) {
    54       int Ndof_min = 3;
    55       int order_max = 0.5*(sqrt(4*Nmatch - 4*Ndof_min + 1) - 3);
    56       int order_use = MIN (CHIPORDER, order_max);
    57       if (order_use < 1) {
    58         if (VERBOSE2) fprintf (stderr, "insufficient measurements (%d) for linear fit\n", Nmatch);
    59         image[0].flags |= ID_IMAGE_ASTROM_FEW;
    60         return FALSE;
    61       }
    62       image[0].coords.Npolyterms = order_use;
    63     }
     57    // set the maximum order for the polynomial (based on number of stars kept above)
     58    int order_use = 0;
     59    if (Nkeep >  5) order_use = 1; // 4 stars per polynomial term (per dimension)
     60    if (Nkeep > 30) order_use = 2; // 5 stars per polynomial term (per dimension)
     61    if (Nkeep > 60) order_use = 3; // 6 stars per polynomial term (per dimension)
     62    if (order_use < 1) {
     63      if (VERBOSE2) fprintf (stderr, "insufficient measurements (%d) for linear fit\n", Nkeep);
     64      image[0].flags |= ID_IMAGE_ASTROM_FEW;
     65      return FALSE;
     66    }
     67    image[0].coords.Npolyterms = order_use;
    6468
    6569    if (fit) fit_free (fit);
     
    7175      fit_add (fit, raw[i].X, raw[i].Y, ref[i].L, ref[i].M, raw[i].dPos);
    7276    }
    73 
    74 # if (0)
    75     // check if the fit has enough data points for the polynomial order
    76     skip = FALSE;
    77     switch (image[0].coords.Npolyterms) {
    78       case 0:
    79       case 1:
    80         skip = (fit[0].Npts < 8);
    81         break;
    82       case 2:
    83         skip = (fit[0].Npts < 11);
    84         break;
    85       case 3:
    86         skip = (fit[0].Npts < 15);
    87         break;
    88       default:
    89         fprintf (stderr, "invalid chip order %d\n", image[0].coords.Npolyterms);
    90         skip = TRUE;
    91     }
    92     if (skip) {
    93       if (VERBOSE2) fprintf (stderr, "insufficient measurements (%d) for requested order (%d)\n", fit[0].Npts, image[0].coords.Npolyterms);
    94       fit_free (fit);
    95       image[0].flags |= ID_IMAGE_ASTROM_FEW;
    96       return FALSE;
    97     }
    98 # endif
    9977
    10078    // measure the fit, update the coords & object coordinates
     
    134112  }
    135113
     114  int Ncolor = 0;
     115  float *colorList = NULL;
     116  ALLOCATE (colorList, float, Nmatch);
     117
     118  for (i = 0; i < Nmatch; i++) {
     119    if (!raw[i].mask) continue;
     120    if (!isnan(ref[i].Color)) {
     121      colorList[Ncolor] = ref[i].Color;
     122      Ncolor ++;
     123    }
     124  }
     125  fsort (colorList, Ncolor);
     126  float colorMedian = (Ncolor > 0) ? colorList[(int)(0.5*Ncolor)] : NAN;
     127  image[0].refColor = colorMedian;
     128  free (colorList);
     129
    136130  GetScatterRawRef(&dLsigFull, &dMsigFull, &dRsigFull, &NstatFull, raw, ref, Nmatch, SIGMA_LIM);
    137131  GetScatterRawRef(&dLsig,     &dMsig,     &dRsig,     &Nstat,     raw, ref, Nmatch, IMFIT_SYS_SIGMA_LIM);
     
    151145  image[0].dYpixSys = dMsig;
    152146  image[0].nFitAstrom = fit[0].Npts;
     147
     148  // fprintf (stderr, "%s %6.3f %4d %4d\n", image[0].name, image[0].refColor, Ncolor, image[0].nFitAstrom);
    153149
    154150  if (fit) fit_free (fit);
  • trunk/Ohana/src/relastro/src/GetAstromError.c

    r33652 r36833  
    3030
    3131  code  = GetPhotcodebyCode (measure[0].photcode);
     32  if (!code) return NAN;
    3233
    3334  // do not raise an exception, just send back the result
     
    7273
    7374  code  = GetPhotcodebyCode (measure[0].photcode);
     75  if (!code) return NAN;
    7476
    7577  // do not raise an exception, just send back the result
  • trunk/Ohana/src/relastro/src/ImageOps.c

    r36630 r36833  
    3939// N_onImage was 'Nlist'
    4040// N_ONIMAGE was 'NLIST'
     41
     42int areImagesLoaded () {
     43
     44  if (image) return TRUE;
     45  return FALSE;
     46}
     47
     48int areImagesMatched () {
     49
     50  if (MeasureToImage) return TRUE;
     51  return FALSE;
     52}
    4153
    4254Image *getimages (off_t *N, off_t **line_number) {
     
    333345*/
    334346
     347float getColor (off_t meas, int cat) {
     348
     349  off_t i;
     350
     351  i = MeasureToImage[cat][meas];
     352  if (i == -1) return (NAN);
     353  return (image[i].refColor);
     354}
     355
    335356void plot_images () {
    336357
     
    468489    }
    469490
     491    X = measure[0].Xccd;
     492    Y = measure[0].Yccd;
    470493    if (USE_FIXED_PIXCOORDS) {
    471       X = isfinite(measure[0].Xfix) ? measure[0].Xfix : measure[0].Xccd;
    472       Y = isfinite(measure[0].Yfix) ? measure[0].Yfix : measure[0].Yccd;
    473     } else {
    474       X = measure[0].Xccd;
    475       Y = measure[0].Yccd;
     494      if (isfinite(measure[0].Xfix) && isfinite(measure[0].Yfix)) {
     495        float dX = measure[0].Xfix - measure[0].Xccd;
     496        float dY = measure[0].Yfix - measure[0].Yccd;
     497        if (hypot(dX,dY) < 2.0) {
     498          X = measure[0].Xfix;
     499          Y = measure[0].Yfix;
     500        }
     501      }
    476502    }
    477503    n = measure[0].averef;
     
    560586
    561587  int Noff = NoffRAave + NoffDECave + NoffRAori + NoffDECori;
    562   if (VERBOSE && (Noff > 0)) fprintf (stderr, "Noff ave RA %d, Noff ave DEC %d, Noff ori RA %d, Noff ori DEC %d\n", NoffRAave, NoffDECave, NoffRAori, NoffDECori);
     588  if (VERBOSE2 && (Noff > 0)) fprintf (stderr, "Noff ave RA %d, Noff ave DEC %d, Noff ori RA %d, Noff ori DEC %d\n", NoffRAave, NoffDECave, NoffRAori, NoffDECori);
    563589  saveOffsets (dPos, nPos, im);
    564590
     
    609635    Measure *measure = &catalog[c].measure[m];
    610636
     637    X = measure[0].Xccd;
     638    Y = measure[0].Yccd;
    611639    if (USE_FIXED_PIXCOORDS) {
    612       X = isfinite(measure[0].Xfix) ? measure[0].Xfix : measure[0].Xccd;
    613       Y = isfinite(measure[0].Yfix) ? measure[0].Yfix : measure[0].Yccd;
    614     } else {
    615       X = measure[0].Xccd;
    616       Y = measure[0].Yccd;
     640      if (isfinite(measure[0].Xfix) && isfinite(measure[0].Yfix)) {
     641        float dX = measure[0].Xfix - measure[0].Xccd;
     642        float dY = measure[0].Yfix - measure[0].Yccd;
     643        if (hypot(dX,dY) < 2.0) {
     644          X = measure[0].Xfix;
     645          Y = measure[0].Yfix;
     646        }
     647      }
    617648    }
    618649    n = measure[0].averef;
     
    679710
    680711    /* apply the current image transformation or use the current value of R+dR, D+dD? */
     712    raw[i].X = measure[0].Xccd;
     713    raw[i].Y = measure[0].Yccd;
    681714    if (USE_FIXED_PIXCOORDS) {
    682       raw[i].X = isfinite(measure[0].Xfix) ? measure[0].Xfix : measure[0].Xccd;
    683       raw[i].Y = isfinite(measure[0].Yfix) ? measure[0].Yfix : measure[0].Yccd;
    684     } else {
    685       raw[i].X = measure[0].Xccd;
    686       raw[i].Y = measure[0].Yccd;
    687     }
    688 
     715      if (isfinite(measure[0].Xfix) && isfinite(measure[0].Yfix)) {
     716        float dX = measure[0].Xfix - measure[0].Xccd;
     717        float dY = measure[0].Yfix - measure[0].Yccd;
     718        if (hypot(dX,dY) < 2.0) {
     719          raw[i].X = measure[0].Xfix;
     720          raw[i].Y = measure[0].Yfix;
     721        }
     722      }
     723    }
    689724    raw[i].Mag  = measure[0].M;
    690725    raw[i].dMag = measure[0].dM;
     
    738773  StarData *ref;
    739774
     775  int Nsecfilt = GetPhotcodeNsecfilt();
     776
    740777  mosaic = NULL;
    741778  moscoords = NULL;
     
    765802    ref[i].dMag = measure[0].dM;
    766803    ref[i].dPos = GetAstromErrorTiny (&measure[0], ERROR_MODE_POS);
     804
     805    if ((DCR_NSEC_POS >= 0) && (DCR_NSEC_NEG >= -1)) {
     806      ref[i].Color = catalog[c].secfilt[n*Nsecfilt + DCR_NSEC_POS].M - catalog[c].secfilt[n*Nsecfilt + DCR_NSEC_NEG].M;
     807    }
    767808
    768809    ref[i].mask = FALSE;
     
    859900    N = 0;
    860901    for (k = 0; k < catalog[0].average[j].Nmeasure; k++) {
    861       //reset flag on each invocation
     902      // reset flag on each invocation
    862903      measure[k].dbFlags &= ~ID_MEAS_POOR_ASTROM;
    863904
     
    898939  double *R, *D, *dR, *dD, *d2;
    899940  StatType statsR, statsD;
     941
     942  // XXX we are not going to use this for now
     943  return;
    900944
    901945  Nsecfilt = GetPhotcodeNsecfilt();
     
    10631107    } else {
    10641108      code = GetPhotcodebyCode (measure[0].photcode);
     1109      if (!code) return FALSE;
    10651110      mask = code[0].astromBadMask;
    10661111    }
     
    11291174    } else {
    11301175      code = GetPhotcodebyCode (measure[0].photcode);
     1176      if (!code) return FALSE;
    11311177      mask = code[0].astromBadMask;
    11321178    }
  • trunk/Ohana/src/relastro/src/UpdateObjectOffsets.c

    r35763 r36833  
    126126              table->hosts[i].hostID, CATDIR, table->hosts[i].pathname, UserPatch.Rmin, UserPatch.Rmax, UserPatch.Dmin, UserPatch.Dmax, STATMODE, MIN_ERROR);
    127127
    128     char tmpline[1024];
    129     if (FIT_MODE == FIT_PM_ONLY)         { snprintf (tmpline, 1024, "%s -pm",    command);           strcpy (command, tmpline); }
    130     if (FIT_MODE == FIT_PAR_ONLY)        { snprintf (tmpline, 1024, "%s -par",   command);           strcpy (command, tmpline); }
    131     if (FIT_MODE == FIT_PM_AND_PAR)      { snprintf (tmpline, 1024, "%s -pmpar", command);           strcpy (command, tmpline); }
    132 
    133     if (VERBOSE)       { snprintf (tmpline, 1024, "%s -v",              command);                    strcpy (command, tmpline); }
    134     if (VERBOSE2)      { snprintf (tmpline, 1024, "%s -vv",             command);                    strcpy (command, tmpline); }
    135     if (RESET)         { snprintf (tmpline, 1024, "%s -reset",          command);                    strcpy (command, tmpline); }
    136     if (UPDATE)        { snprintf (tmpline, 1024, "%s -update",         command);                    strcpy (command, tmpline); }
    137 
    138     if (RESET_BAD_IMAGES) { snprintf (tmpline, 1024, "%s -reset-bad-images", command);               strcpy (command, tmpline); }
    139 
    140     if (ImagSelect)    { snprintf (tmpline, 1024, "%s -instmag %f %f",  command, ImagMin, ImagMax);  strcpy (command, tmpline); }
    141     if (MaxDensityUse) { snprintf (tmpline, 1024, "%s -max-density %f", command, MaxDensityValue);   strcpy (command, tmpline); }
     128    if (FIT_MODE == FIT_PM_ONLY)         strextend (command, "-pm");
     129    if (FIT_MODE == FIT_PAR_ONLY)        strextend (command, "-par");
     130    if (FIT_MODE == FIT_PM_AND_PAR)      strextend (command, "-pmpar");
     131
     132    if (VERBOSE)       strextend (command, "-v");
     133    if (VERBOSE2)      strextend (command, "-vv");
     134    if (RESET)         strextend (command, "-reset");
     135    if (UPDATE)        strextend (command, "-update");
     136
     137    if (RESET_BAD_IMAGES) strextend (command, "-reset-bad-images");
     138
     139    if (ImagSelect)    strextend (command, "-instmag %f %f", ImagMin, ImagMax);
     140    if (MaxDensityUse) strextend (command, "-max-density %f", MaxDensityValue);
    142141   
    143     if (USE_BASIC_CHECK) { snprintf (tmpline, 1024, "%s -basic-image-search", command);              strcpy (command, tmpline); }
    144     if (FlagOutlier)   { snprintf (tmpline, 1024, "%s -clip %d",        command, CLIP_THRESH);       strcpy (command, tmpline); }
     142    if (USE_BASIC_CHECK) strextend (command, "-basic-image-search");
     143    if (FlagOutlier)     strextend (command, "-clip %d", CLIP_THRESH);
     144    if (ExcludeBogus)    strextend (command, "-exclude-bogus %f", ExcludeBogusRadius);
    145145   
    146     if (USE_FIXED_PIXCOORDS) { snprintf (tmpline, 1024, "%s -D USE_FIXED_PIXCOORDS 1", command);     strcpy (command, tmpline); }
    147 
    148     if (PHOTCODE_KEEP_LIST) { snprintf (tmpline, 1024, "%s +photcode %s", command, PHOTCODE_KEEP_LIST); strcpy (command, tmpline); }
    149     if (PHOTCODE_SKIP_LIST) { snprintf (tmpline, 1024, "%s -photcode %s", command, PHOTCODE_SKIP_LIST); strcpy (command, tmpline); }
    150     if (PhotFlagSelect)    { snprintf (tmpline, 1024, "%s +photflags",   command);                     strcpy (command, tmpline); }
    151     if (PhotFlagBad)       { snprintf (tmpline, 1024, "%s +photflagbad %d", command, PhotFlagBad);     strcpy (command, tmpline); }
    152     if (PhotFlagPoor)      { snprintf (tmpline, 1024, "%s +photflagpoor %d", command, PhotFlagPoor);   strcpy (command, tmpline); }
     146    if (USE_FIXED_PIXCOORDS) strextend (command, "-D USE_FIXED_PIXCOORDS 1");
     147
     148    if (PHOTCODE_KEEP_LIST) strextend (command, "+photcode %s", PHOTCODE_KEEP_LIST);
     149    if (PHOTCODE_SKIP_LIST) strextend (command, "-photcode %s", PHOTCODE_SKIP_LIST);
     150    if (PhotFlagSelect)     strextend (command, "+photflags");
     151    if (PhotFlagBad)        strextend (command, "+photflagbad %d", PhotFlagBad);
     152    if (PhotFlagPoor)       strextend (command, "+photflagpoor %d", PhotFlagPoor);
    153153    // XXX note that the above pass in the flag as decimal -- also note that args.c cannot handle 0xHEX values
    154154
    155     if (MinBadQF > 0.0)    { snprintf (tmpline, 1024, "%s -min-bad-psfqf %f",        command, MinBadQF);      strcpy (command, tmpline); }
    156     if (MaxMeanOffset != 10.0) { snprintf (tmpline, 1024, "%s -max-mean-offset  %f", command, MaxMeanOffset); strcpy (command, tmpline); }
     155    if (DCR_COLOR_POS && DCR_COLOR_NEG) {
     156      strextend (command, "-dcr-color %s %s", DCR_COLOR_POS, DCR_COLOR_NEG);
     157    }
     158
     159    if (MinBadQF > 0.0)        strextend (command, "-min-bad-psfqf %f", MinBadQF);
     160    if (MaxMeanOffset != 10.0) strextend (command, "-max-mean-offset  %f", MaxMeanOffset);
    157161
    158162    if (TimeSelect) {
    159163      char *tstart = ohana_sec_to_date (TSTART);
    160164      char *tstop  = ohana_sec_to_date (TSTOP);
    161       snprintf (tmpline, 1024, "%s -time %s %s", command, tstart, tstop);
     165      strextend (command, "-time %s %s", tstart, tstop);
    162166      free (tstart);
    163167      free (tstop);
    164       strcpy (command, tmpline);
    165168    }
    166169    fprintf (stderr, "command: %s\n", command);
  • trunk/Ohana/src/relastro/src/UpdateObjects.c

    r36630 r36833  
    1010static double *T;
    1111static double *dT;
     12static double *C;
    1213
    1314void initObjectData (Catalog *catalog, int Ncatalog) {
     
    3637  ALLOCATE (pX, double, MAX (1, Nmax));
    3738  ALLOCATE (pY, double, MAX (1, Nmax));
     39
     40  ALLOCATE (C, double, MAX (1, Nmax));
    3841
    3942
     
    5457  free (pX);
    5558  free (pY);
     59
     60  free (C);
    5661
    5762
     
    6166
    6267  off_t j, k, m;
    63   int i, N, Nsecfilt, mode, result, status, XVERB;
     68  int i, N, Nc, Nsecfilt, mode, result, status, XVERB;
    6469  StatType statsR, statsD;
    6570  Coords coords;
     
    7580  memset (&fitPAR, 0, sizeof(fitPAR));
    7681  initObjectData (catalog, Ncatalog);
     82
     83  int setRefColor = areImagesMatched();
    7784
    7885  /* project coordinates to a plane centered on the object with units of arcsec */
     
    112119      fitPAR.chisq = NAN;
    113120
     121      // if we fail to fit the astrometry for some reason, we need to set/reset these
     122      catalog[i].average[j].flags |= ID_STAR_NO_ASTROM;
     123      catalog[i].average[j].ChiSqAve  = NAN;
     124      catalog[i].average[j].ChiSqPM   = NAN;
     125      catalog[i].average[j].ChiSqPar  = NAN;
     126      catalog[i].average[j].Npos = 0;
     127
    114128      if (catalog[i].average[j].Nmeasure == 0) {
    115           continue;
    116       }
    117 
     129        continue;
     130      }
     131
     132      Nc = 0;
    118133      N = 0;
    119134      m = catalog[i].average[j].measureOffset;
     
    138153        }
    139154
     155        // reset the bit to note that a detection was used (or not)
     156        measure[k].dbFlags &= ~ID_MEAS_USED_OBJ;
     157        if (measureBig) { measureBig[k].dbFlags &= ~ID_MEAS_USED_OBJ; }
     158
    140159        // does the measurement pass the supplied filtering constraints?
    141160        // MeasFilterTestTiny does not test psfQF
     161        // exclude bad detections based on: photcodes, psfQF, time range, photflags & astromBadMask, mag_inst
    142162        int keepMeasure = measureBig ? MeasFilterTest(&measureBig[k], FALSE) : MeasFilterTestTiny(&measure[k], FALSE);
    143163        if (!keepMeasure) {
    144           measure[k].dbFlags &= ~ID_MEAS_USED_OBJ;
    145           if (measureBig) { measureBig[k].dbFlags &= ~ID_MEAS_USED_OBJ; }
    146164          continue;
    147165        }
    148166
     167        // mark (as POOR) any measurements which are deviant from the mean by > ExcludeBogusRadius
     168        if (ExcludeBogus) {
     169          double Ri = getMeanR (&measure[k], &catalog[i].average[j], &catalog[i].secfilt[j*Nsecfilt]);
     170          double Di = getMeanD (&measure[k], &catalog[i].average[j], &catalog[i].secfilt[j*Nsecfilt]);
     171          coords.crval1 = catalog[i].average[j].R;
     172          coords.crval2 = catalog[i].average[j].D;
     173          double Xi, Yi;
     174          RD_to_XY (&Xi, &Yi, Ri, Di, &coords);
     175          double radius = hypot(Xi, Yi);
     176          if (radius > ExcludeBogusRadius) {
     177              measure[k].dbFlags |= ID_MEAS_POOR_ASTROM;
     178              if (measureBig) { measureBig[k].dbFlags |= ID_MEAS_POOR_ASTROM; }
     179              continue;
     180          }
     181          measure[k].dbFlags &= ~ID_MEAS_POOR_ASTROM;
     182          if (measureBig) { measureBig[k].dbFlags &= ~ID_MEAS_POOR_ASTROM; }
     183        }
     184
    149185        // outlier rejection
    150         if (FlagOutlier && (measure[k].dbFlags & ID_MEAS_POOR_ASTROM)) {
    151           measure[k].dbFlags &= ~ID_MEAS_USED_OBJ;
    152           if (measureBig) { measureBig[k].dbFlags &= ~ID_MEAS_USED_OBJ; }
     186        if (FALSE && FlagOutlier && (measure[k].dbFlags & ID_MEAS_POOR_ASTROM)) {
    153187          continue;
    154188        }
    155 
    156         measure[k].dbFlags |= ID_MEAS_USED_OBJ;
    157         if (measureBig) { measureBig[k].dbFlags |= ID_MEAS_USED_OBJ; }
    158189
    159190        R[N] = getMeanR (&measure[k], &catalog[i].average[j], &catalog[i].secfilt[j*Nsecfilt]);
     
    170201        // allow a given photcode or measurement to be
    171202        // ignored if the error is NAN (for photcode, set astromErrSys to NaN)
    172         if (isnan(dX[N])) {
    173           measure[k].dbFlags &= ~ID_MEAS_USED_OBJ;
    174           if (measureBig) { measureBig[k].dbFlags &= ~ID_MEAS_USED_OBJ; }
    175           continue;
    176         }
    177         if (isnan(dY[N])) {
    178           measure[k].dbFlags &= ~ID_MEAS_USED_OBJ;
    179           if (measureBig) { measureBig[k].dbFlags &= ~ID_MEAS_USED_OBJ; }
    180           continue;
    181         }
     203        if (isnan(dX[N])) continue;
     204        if (isnan(dY[N])) continue;
    182205
    183206        // add systematic error in quadrature, if desired
     
    201224        dD[N] = dY[N] / 3600.0;
    202225
     226        if (setRefColor) {
     227          float color = getColor (m+k, i);
     228          if (!isnan(color)) {
     229            C[Nc] = color;
     230            Nc++;
     231          }
     232        }
     233
     234        measure[k].dbFlags |= ID_MEAS_USED_OBJ;
     235        if (measureBig) { measureBig[k].dbFlags |= ID_MEAS_USED_OBJ; }
     236
    203237        N++;
    204       }
     238      } // loop over measurements : catalog[i].average[j].Nmeasure
    205239
    206240      // if we have too few good detections for the desired fit, or too limited a
    207241      // baseline, use a fit with fewer parameters.  XXX if we have too few measurements
    208242      // for even the average position, consider including the lower-quality detections?
    209 
    210       catalog[i].average[j].flags &= ~ID_STAR_NO_ASTROM;
    211243
    212244      // find Tmin & Tmax from the list of accepted measurements
     
    231263      // too few measurements for average position (require 2 values)
    232264      if (N < SRC_MEAS_TOOFEW) {
    233         // XXX need to define PHOTOM and ASTROM object flags
    234         // XXX reset the average value fields?
    235         catalog[i].average[j].flags |= ID_STAR_NO_ASTROM;
    236         catalog[i].average[j].ChiSqAve  = NAN;
    237         catalog[i].average[j].ChiSqPM   = NAN;
    238         catalog[i].average[j].ChiSqPar  = NAN;
    239265        if (N < 2) continue;
    240266      }
     
    268294        catalog[i].average[j].flags |= ID_STAR_FIT_PM;
    269295        Npm ++;
     296
     297        // XXX a hard-wired hack...
     298        if ((fabs(fitPM.uR) > 2.0) || (fabs(fitPM.uD) > 2.0)) {
     299          mode = FIT_AVERAGE;
     300          catalog[i].average[j].flags |= ID_STAR_BAD_PM;
     301        }
    270302      }
    271303
    272304      if (mode == FIT_PM_AND_PAR) {
    273         // fprintf (stderr, "parallax fitting is still untested (%s, %d)\n", __FILE__, __LINE__);
    274 
    275305        float pXmin = +2.0;
    276306        float pXmax = -2.0;
     
    296326          catalog[i].average[j].flags |= ID_STAR_FIT_PAR;
    297327          Npar ++;
     328
     329          // XXX a hard-wired hack...
     330          if ((fabs(fitPAR.uR) > 2.0) || (fabs(fitPAR.uD) > 2.0)) {
     331            mode = FIT_AVERAGE;
     332            catalog[i].average[j].flags |= ID_STAR_BAD_PM;
     333          }
    298334        } else {
    299335          // need to set mode = FIT_PM_ONLY if we do not fit for parallax
     
    321357        catalog[i].average[j].flags |= ID_STAR_FIT_AVE;
    322358        Nave ++;
     359      }
     360
     361      if (setRefColor) {
     362        dsort (C, Nc);
     363        float colorMedian = (Nc > 0) ? C[(int)(0.5*Nc)] : NAN;
     364        catalog[i].average[j].refColor = colorMedian;
    323365      }
    324366
     
    411453      catalog[i].average[j].Trange = (Trange * 86400 * 365.25);
    412454      catalog[i].average[j].Npos = fit.Nfit;
     455
     456      // unset the NO_ASTROM bit (not(NO_ASTROM) == HAVE_ASTROM)
     457      catalog[i].average[j].flags &= ~ID_STAR_NO_ASTROM;
     458
    413459      if (XVERB) fprintf (stderr, "%f %f -> %f %f (%f,%f) pm=(%f %f) chisq=(%f, %f, %f)\n",
    414460                          catalog[i].average[j].R,
  • trunk/Ohana/src/relastro/src/args.c

    r36630 r36833  
    179179  FlagOutlier = FALSE;
    180180  if ((N = get_argument (argc, argv, "-clip"))) {
     181    fprintf (stderr, "-clip is currently disabled\n");
    181182    remove_argument (N, &argc, argv);
    182183    CLIP_THRESH = atof (argv[N]);
    183184    remove_argument (N, &argc, argv);
    184185    FlagOutlier = TRUE;
     186  }
     187
     188  ExcludeBogus = FALSE;
     189  ExcludeBogusRadius = 0.0;
     190  if ((N = get_argument (argc, argv, "-exclude-bogus"))) {
     191    remove_argument (N, &argc, argv);
     192    ExcludeBogusRadius = atof (argv[N]);
     193    remove_argument (N, &argc, argv);
     194    ExcludeBogus = TRUE;
    185195  }
    186196
     
    284294  }
    285295
     296  DCR_COLOR_POS = NULL;
     297  DCR_COLOR_NEG = NULL;
     298  if ((N = get_argument (argc, argv, "-dcr-color"))) {
     299    remove_argument (N, &argc, argv);
     300    DCR_COLOR_POS = strcreate(argv[N]);
     301    remove_argument (N, &argc, argv);
     302    DCR_COLOR_NEG = strcreate(argv[N]);
     303    remove_argument (N, &argc, argv);
     304  }
     305
     306  PHOTCODE_RESET_LIST = NULL;
     307  if ((N = get_argument (argc, argv, "-reset-to-photcode"))) {
     308    remove_argument (N, &argc, argv);
     309    PHOTCODE_RESET_LIST = strcreate(argv[N]);
     310    remove_argument (N, &argc, argv);
     311  }
     312
    286313  PHOTCODE_KEEP_LIST = NULL;
    287314  if ((N = get_argument (argc, argv, "+photcode"))) {
     
    291318  }
    292319
    293   PHOTCODE_SKIP_LIST = strcreate("SCOS.103a.E,SCOS.4414.OG590,SCOS.4415.OG590,SCOS.IIIaF.OG590,SCOS.IIIaF.RG610,SCOS.IIIaF.RG630,SCOS.IIIaJ.GG385,SCOS.IIIaJ.GG395,SCOS.IVN.RG715,SCOS.IVN.RG9");
     320  char *SuperCOSMOS_SKIP = strcreate("SCOS.103a.E,SCOS.4414.OG590,SCOS.4415.OG590,SCOS.IIIaF.OG590,SCOS.IIIaF.RG610,SCOS.IIIaF.RG630,SCOS.IIIaJ.GG385,SCOS.IIIaJ.GG395,SCOS.IVN.RG715,SCOS.IVN.RG9");
     321
     322  PHOTCODE_SKIP_LIST = NULL;
    294323  if ((N = get_argument (argc, argv, "-photcode"))) {
    295324    remove_argument (N, &argc, argv);
    296     char *tmp1 = strcreate(argv[N]);
    297 
    298     int Ntotal = strlen(tmp1) + strlen(PHOTCODE_SKIP_LIST) + 5;
    299 
    300     char *tmp2 = NULL;
    301     ALLOCATE (tmp2, char, Ntotal);
    302     snprintf (tmp2, Ntotal, "%s,%s", PHOTCODE_SKIP_LIST, tmp1);
    303 
    304     free (tmp1);
    305     free (PHOTCODE_SKIP_LIST);
    306 
    307     PHOTCODE_SKIP_LIST = tmp2;
     325    char *RawSkip = strcreate(argv[N]);
     326
     327    char *GotSkip = strstr (RawSkip, SuperCOSMOS_SKIP);
     328    if (!GotSkip) {
     329      int Ntotal = strlen(RawSkip) + strlen(SuperCOSMOS_SKIP) + 5;
     330      ALLOCATE (PHOTCODE_SKIP_LIST, char, Ntotal);
     331      snprintf (PHOTCODE_SKIP_LIST, Ntotal, "%s,%s", SuperCOSMOS_SKIP, RawSkip);
     332      free (RawSkip);
     333      free (SuperCOSMOS_SKIP);
     334    } else {
     335      PHOTCODE_SKIP_LIST = RawSkip;
     336    }
    308337    remove_argument (N, &argc, argv);
    309338  }
     
    582611  FlagOutlier = FALSE;
    583612  if ((N = get_argument (argc, argv, "-clip"))) {
     613    fprintf (stderr, "-clip is currently disabled\n");
    584614    remove_argument (N, &argc, argv);
    585615    CLIP_THRESH = atof (argv[N]);
    586616    remove_argument (N, &argc, argv);
    587617    FlagOutlier = TRUE;
     618  }
     619
     620  ExcludeBogus = FALSE;
     621  ExcludeBogusRadius = 0.0;
     622  if ((N = get_argument (argc, argv, "-exclude-bogus"))) {
     623    remove_argument (N, &argc, argv);
     624    ExcludeBogusRadius = atof (argv[N]);
     625    remove_argument (N, &argc, argv);
     626    ExcludeBogus = TRUE;
    588627  }
    589628
     
    655694    remove_argument (N, &argc, argv);
    656695    TimeSelect = TRUE;
     696  }
     697
     698  DCR_COLOR_POS = NULL;
     699  DCR_COLOR_NEG = NULL;
     700  if ((N = get_argument (argc, argv, "-dcr-color"))) {
     701    remove_argument (N, &argc, argv);
     702    DCR_COLOR_POS = strcreate(argv[N]);
     703    remove_argument (N, &argc, argv);
     704    DCR_COLOR_NEG = strcreate(argv[N]);
     705    remove_argument (N, &argc, argv);
     706  }
     707
     708  PHOTCODE_RESET_LIST = NULL;
     709  if ((N = get_argument (argc, argv, "-reset-to-photcode"))) {
     710    remove_argument (N, &argc, argv);
     711    PHOTCODE_RESET_LIST = strcreate(argv[N]);
     712    remove_argument (N, &argc, argv);
    657713  }
    658714
     
    861917  exit (2);
    862918}
     919
     920int strextend (char *input, char *format,...) {
     921
     922  char tmpextra[1024], tmpline[1024];
     923  va_list argp;
     924
     925  va_start (argp, format);
     926  vsnprintf (tmpextra, 1024, format, argp);
     927  snprintf (tmpline, 1024, "%s %s", input, tmpextra);
     928  strcpy (input, tmpline);
     929
     930  return TRUE;
     931}
  • trunk/Ohana/src/relastro/src/bcatalog.c

    r34429 r36833  
    11# include "relastro.h"
    22
     3static int Nkeep1 = 0;
     4static int Nkeep2 = 0;
    35static int Nskip1 = 0;
    46static int Nskip2 = 0;
    5 static unsigned int Tref = 1323353985;
    6 static short Cref = 10355;
     7
     8FILE *fbogus = NULL;
     9static int NskipBogus = 0;
     10
     11// test image: 2013/06/15,13:25:51, GPC1.r.XY50
     12static int CHECK_TEST_IMAGE = FALSE;
     13// static unsigned int Tref = 1378812312;
     14// static short Cref = 10001;
     15
     16static unsigned int Tref = 1379570672;
     17static short Cref = 10341;
     18
     19int LimitDensityCatalog_ByNmeasureGrid (Catalog *subcatalog, Catalog *oldcatalog);
    720
    821int bcatalog (Catalog *subcatalog, Catalog *catalog) {
     
    1124  off_t NAVERAGE, NMEASURE, Naverage, Nmeasure, Nm;
    1225  int Nsecfilt;
     26  Coords coords;
     27
     28  /* for outlier rejection, project coordinates to a plane centered on the object with units of arcsec */
     29  coords.crval1 = 0;
     30  coords.crval2 = 0;
     31  coords.crpix1 = 0;
     32  coords.crpix2 = 0;
     33  coords.cdelt1 = coords.cdelt2 = 1.0 / 3600.0;
     34  coords.pc1_1  = coords.pc2_2 = 1.0;
     35  coords.pc1_2  = coords.pc2_1 = 0.0;
     36  coords.Npolyterms = 1;
     37  strcpy (coords.ctype, "DEC--SIN");
    1338
    1439  // XXX in the future, use catalog[0].Nsecfilt only?  allow catalogs to have variable Nsecfilt?
     
    3459    ID_STAR_NO_ASTROM ;
    3560
     61  if (VERBOSE2 && ExcludeBogus && (fbogus == NULL)) {
     62    char name[1024];
     63    snprintf (name, 1024, "%s/bogus.%02d.dat", CATDIR, HOST_ID);
     64    fbogus = fopen (name, "w");
     65    if (!fbogus) {
     66      fprintf (stderr, "trouble opening bogus detection dump : %s\n", name);
     67    }
     68  }
     69
    3670  /* exclude stars not in range or with too few measurements */
    3771  for (i = 0; i < catalog[0].Naverage; i++) {
     
    6094      if (!MeasFilterTest(&catalog[0].measure[offset], TRUE)) {
    6195        catalog[0].measure[offset].dbFlags &= ~ID_MEAS_USED_CHIP;
    62         if (FALSE && (abs(catalog[0].measure[offset].t - Tref) < 10) && (catalog[0].measure[offset].photcode == Cref)) {
     96        if (CHECK_TEST_IMAGE && (abs(catalog[0].measure[offset].t - Tref) < 10) && (catalog[0].measure[offset].photcode == Cref)) {
    6397          Nskip1 ++;
    6498        }
    6599        continue;
    66100      }
    67 
    68       // filter out outliers
    69       if (FlagOutlier && (catalog[0].measure[offset].dbFlags & ID_MEAS_POOR_ASTROM)) {
     101      if (CHECK_TEST_IMAGE && (abs(catalog[0].measure[offset].t - Tref) < 10) && (catalog[0].measure[offset].photcode == Cref)) {
     102        Nkeep1 ++;
     103      }
     104
     105      // filter out outliers - these are detections inconsistent with the offset distribution
     106      // XXX disable this for now
     107      if (FALSE && FlagOutlier && (catalog[0].measure[offset].dbFlags & ID_MEAS_POOR_ASTROM)) {
    70108        catalog[0].measure[offset].dbFlags &= ~ID_MEAS_USED_CHIP;
    71109        if (FALSE && (abs(catalog[0].measure[offset].t - Tref) < 10) && (catalog[0].measure[offset].photcode == Cref)) {
     
    76114      catalog[0].measure[offset].dbFlags |= ID_MEAS_USED_CHIP;
    77115
     116      // exclude bogus
     117      if (ExcludeBogus) {
     118          double Ri = getMeanR_Big (&catalog[0].measure[offset], &catalog[0].average[i], &catalog[0].secfilt[i*Nsecfilt]);
     119          double Di = getMeanD_Big (&catalog[0].measure[offset], &catalog[0].average[i], &catalog[0].secfilt[i*Nsecfilt]);
     120          coords.crval1 = catalog[0].average[i].R;
     121          coords.crval2 = catalog[0].average[i].D;
     122          double Xi, Yi;
     123          RD_to_XY (&Xi, &Yi, Ri, Di, &coords);
     124          double radius = hypot(Xi, Yi);
     125          if (radius > ExcludeBogusRadius) {
     126              NskipBogus ++;
     127              if (VERBOSE2) {
     128                FILE *foutput = fbogus ? fbogus : stderr;
     129                char *date = ohana_sec_to_date(catalog[0].measure[offset].t);
     130                fprintf (foutput, "exclude bogus: %10.6f %10.6f : %10.6f %10.6f : %6.2f %6.2f : %6.2f : %5d %s\n", catalog[0].average[i].R, catalog[0].average[i].D, Ri, Di, Xi, Yi, radius, catalog[0].measure[offset].photcode, date);
     131                free (date);
     132              }
     133              continue;
     134          }
     135      }
     136
    78137      // re-assess on each run of relastro if a measurement should be used
    79138
     
    87146      // but before the final average properties are calculated, these measurements may be
    88147      // allowed.
     148
     149      if (CHECK_TEST_IMAGE && (abs(catalog[0].measure[offset].t - Tref) < 10) && (catalog[0].measure[offset].photcode == Cref)) {
     150        Nkeep2 ++;
     151      }
    89152
    90153      CopyMeasureToTiny (&subcatalog[0].measureT[Nmeasure], &catalog[0].measure[offset]);
     
    125188  // limit the total number of stars in the catalog
    126189  if (MaxDensityUse) {
    127     LimitDensityCatalog_ByNmeasure (subcatalog, catalog);
     190    LimitDensityCatalog_ByNmeasureGrid (subcatalog, catalog);
    128191  } else {
    129192    if (VERBOSE2) {
     
    134197    }
    135198  }
     199  if (CHECK_TEST_IMAGE && (Nkeep1 + Nkeep2 + Nskip1 + Nskip2 > 0)) {
     200    fprintf (stderr, "kept %d %d, skipped %d %d\n", Nkeep1, Nkeep2, Nskip1, Nskip2);
     201  }
     202
    136203  return (TRUE);
    137204}
    138205
    139206void bcatalog_show_skips () {
    140   fprintf (stderr, "Nskip: %d, %d\n", Nskip1, Nskip2);
     207  if (ExcludeBogus) {
     208    if (fbogus) fclose (fbogus);
     209    fprintf (stderr, "NskipBogus: %d\n", NskipBogus);
     210    // fprintf (stderr, "Nskip: %d, %d\n", Nskip1, Nskip2);
     211    // fprintf (stderr, "Nkeep: %d, %d\n", Nkeep1, Nkeep2);
     212  }
    141213}
    142214
     
    155227# undef COMPARE
    156228
     229}
     230
     231/* generate a grid in a locally-projected coordinate system, try to select average entries
     232   from each grid cell in decending Nmeasure order.
     233 */
     234int LimitDensityCatalog_ByNmeasureGrid (Catalog *subcatalog, Catalog *oldcatalog) {
     235
     236  off_t i, j;
     237  int ix, iy;
     238
     239  Catalog tmpcatalog;
     240
     241  double Rmin, Rmax, Dmin, Dmax;
     242
     243  int Nsecfilt = GetPhotcodeNsecfilt ();
     244
     245  gfits_scan (&oldcatalog[0].header, "RA0",  "%lf", 1, &Rmin);
     246  gfits_scan (&oldcatalog[0].header, "DEC0", "%lf", 1, &Dmin);
     247  gfits_scan (&oldcatalog[0].header, "RA1",  "%lf", 1, &Rmax);
     248  gfits_scan (&oldcatalog[0].header, "DEC1", "%lf", 1, &Dmax);
     249
     250  if (VERBOSE2) fprintf (stderr, "extracting from catalog covering region %f,%f to %f,%f\n", Rmin, Dmin, Rmax, Dmax);
     251
     252  float AREA = fabs(Dmax - Dmin) * fabs(Rmax - Rmin) * cos (0.5*RAD_DEG*(Dmax + Dmin));
     253  assert (AREA > 0);
     254
     255  off_t Nmax = MaxDensityValue * AREA;
     256  if (subcatalog[0].Naverage <= Nmax) {
     257    if (VERBOSE) {
     258      fprintf (stderr, "subcatalog has less than the max density\n");
     259    }
     260    return (TRUE);
     261  }
     262
     263  off_t Naverage = subcatalog[0].Naverage;
     264
     265  // generate a grid in locally projected space
     266  double Rc = 0.5*(Rmin + Rmax);
     267  double Dc = 0.5*(Dmin + Dmax);
     268
     269  /* project coordinates to a plane centered on the object with units of arcsec */
     270  Coords coords;
     271  coords.crval1 = Rc;
     272  coords.crval2 = Dc;
     273  coords.crpix1 = 0;
     274  coords.crpix2 = 0;
     275  coords.cdelt1 = coords.cdelt2 = 1.0 / 3600.0;
     276  coords.pc1_1  = coords.pc2_2 = 1.0;
     277  coords.pc1_2  = coords.pc2_1 = 0.0;
     278  coords.Npolyterms = 1;
     279  strcpy (coords.ctype, "DEC--SIN");
     280
     281  // convert all average R,D values to X,Y:
     282  double *X, *Y;
     283  ALLOCATE (X, double, Naverage);
     284  ALLOCATE (Y, double, Naverage);
     285  float Xmin = +10000.0, Ymin = +10000.0;
     286  float Xmax = -10000.0, Ymax = -10000.0;
     287  for (i = 0; i < Naverage; i++) {
     288    X[i] = NAN;
     289    Y[i] = NAN;
     290    // skip any stars which are outside of nominal catalog range
     291    if (subcatalog[0].average[i].R < Rmin) continue;
     292    if (subcatalog[0].average[i].R > Rmax) continue;
     293    if (subcatalog[0].average[i].D < Dmin) continue;
     294    if (subcatalog[0].average[i].D > Dmax) continue;
     295    RD_to_XY (&X[i], &Y[i], subcatalog[0].average[i].R, subcatalog[0].average[i].D, &coords);
     296    Xmin = MIN (Xmin, X[i]);
     297    Xmax = MAX (Xmax, X[i]);
     298    Ymin = MIN (Ymin, Y[i]);
     299    Ymax = MAX (Ymax, Y[i]);
     300  }
     301
     302  // how many grid cells? what is the grid spacing?
     303  float dX = Xmax - Xmin;
     304  float dY = Ymax - Ymin;
     305
     306  // *** XXX for the moment, I'm using a hard-wired cell size (200 arcsec ~ 3.3 arcmin)
     307  int NX = (int)(dX / 200) + 1;
     308  int NY = (int)(dY / 200) + 1;
     309  // fprintf (stderr, "Density Grid: %d x %d\n", NX, NY);
     310  // XXX check that NX,NY are sensible (5 degrees / 200 arcsec seems like the absolute max)
     311  if (NX > 1000) {
     312    fprintf (stderr, "serious problem with %s: NX = %d\n", subcatalog[0].filename, NX);
     313    exit (3);
     314  }
     315  if (NY > 1000) {
     316    fprintf (stderr, "serious problem with %s: NY = %d\n", subcatalog[0].filename, NY);
     317    exit (3);
     318  }
     319 
     320  // kind of ugly : generate a grid of index, Nmeasure arrays
     321  // to be filled below (I also need NN and Nn to track the number of
     322  // entries in each).
     323  int    **NN_grid;
     324  int    **Nn_grid;
     325  int   ***Nm_grid;
     326  off_t ***idxgrid;
     327  ALLOCATE (NN_grid, int *, NX);
     328  ALLOCATE (Nn_grid, int *, NX);
     329  ALLOCATE (Nm_grid, int **, NX);
     330  ALLOCATE (idxgrid, off_t **, NX);
     331
     332  for (ix = 0; ix < NX; ix++) {
     333    ALLOCATE (NN_grid[ix], int, NY);
     334    ALLOCATE (Nn_grid[ix], int, NY);
     335    ALLOCATE (Nm_grid[ix], int *, NY);
     336    ALLOCATE (idxgrid[ix], off_t *, NY);
     337    for (iy = 0; iy < NY; iy++) {
     338      Nn_grid[ix][iy] = 0;
     339      NN_grid[ix][iy] = 100;
     340      ALLOCATE (Nm_grid[ix][iy], int,   NN_grid[ix][iy]);
     341      ALLOCATE (idxgrid[ix][iy], off_t, NN_grid[ix][iy]);
     342    }
     343  }
     344
     345  // assign all of the average entries to a grid cell
     346  for (i = 0; i < Naverage; i++) {
     347    if (isnan(X[i])) continue;
     348    if (isnan(Y[i])) continue;
     349    ix = (X[i] - Xmin) / 200.0;
     350    iy = (Y[i] - Ymin) / 200.0;
     351    int Nn = Nn_grid[ix][iy];
     352    Nm_grid[ix][iy][Nn] = subcatalog[0].average[i].Nmeasure;
     353   
     354    // if we are resetting to a given photcode, we need to have that photcode...
     355    if (NphotcodesReset) {
     356      int k;
     357      int foundReset = FALSE;
     358      int m = subcatalog[0].average[i].measureOffset;
     359      MeasureTiny *measure = &subcatalog[0].measureT[m];
     360      for (j = 0; (j < subcatalog[0].average[i].Nmeasure) && !foundReset; j++) {
     361        if (CHECK_TEST_IMAGE && (abs(measure[j].t - Tref) < 10) && (measure[j].photcode == Cref)) {
     362          fprintf (stderr, ".");
     363        }
     364        for (k = 0; (k < NphotcodesReset) && !foundReset; k++) {
     365          if (photcodesReset[k][0].code == measure[j].photcode) foundReset = TRUE;
     366        }
     367      }
     368      if (!foundReset) {
     369        Nm_grid[ix][iy][Nn] = 0;
     370      }
     371    }
     372
     373    idxgrid[ix][iy][Nn] = i;
     374    Nn_grid[ix][iy] ++;
     375    if (Nn_grid[ix][iy] >= NN_grid[ix][iy]) {
     376      NN_grid[ix][iy] += 100;
     377      REALLOCATE (Nm_grid[ix][iy], int,   NN_grid[ix][iy]);
     378      REALLOCATE (idxgrid[ix][iy], off_t, NN_grid[ix][iy]);
     379    }
     380  }
     381   
     382  // sort all of the grid cells
     383  for (ix = 0; ix < NX; ix++) {
     384    for (iy = 0; iy < NY; iy++) {
     385      sort_by_Nmeasure (Nm_grid[ix][iy], idxgrid[ix][iy], Nn_grid[ix][iy]);
     386      NN_grid[ix][iy] = 0; // I'm going to use this array to track which element I've already selected
     387    }
     388  }
     389
     390  // cycle over the grid until we ready Nmax
     391  off_t *keepidx = NULL;
     392  ALLOCATE (keepidx, off_t, Naverage);
     393  memset (keepidx, 0, Naverage*sizeof(off_t));
     394  int Nkeep = 0;
     395
     396  for (i = 0; (i < 20) && (Nkeep < Nmax); i++) {
     397    for (ix = 0; (ix < NX) && (Nkeep < Nmax); ix++) {
     398      for (iy = 0; (iy < NY) && (Nkeep < Nmax); iy++) {
     399        if (NN_grid[ix][iy] >= Nn_grid[ix][iy]) continue; // all used up!
     400        int Nn = NN_grid[ix][iy];
     401        keepidx[Nkeep] = idxgrid[ix][iy][Nn];
     402        Nkeep ++;
     403        NN_grid[ix][iy] ++;
     404      }
     405    }
     406  }
     407
     408  // count the number of measurements this selection will yield
     409  off_t ave, NMEASURE = 0;
     410  for (i = 0; i < Nkeep; i++) {
     411    ave = keepidx[i];
     412    NMEASURE += subcatalog[0].average[ave].Nmeasure;
     413  }
     414
     415  // test catID : 37262 37261 37257 37258
     416  int dumpit = FALSE;
     417  dumpit |= (oldcatalog[0].catID == 37007);
     418  // dumpit |= (oldcatalog[0].catID == 37261);
     419  // dumpit |= (oldcatalog[0].catID == 37257);
     420  // dumpit |= (oldcatalog[0].catID == 37258);
     421  if (dumpit) {
     422    char name[64];
     423    snprintf (name, 64, "cat.%05d.dump.dat", oldcatalog[0].catID);
     424    FILE *fdump = fopen (name, "w");
     425    for (i = 0; i < Nkeep; i++) {
     426      ave = keepidx[i];
     427      fprintf (fdump, "%10.6f %10.6f %d\n", subcatalog[0].average[ave].R, subcatalog[0].average[ave].D, subcatalog[0].average[ave].Nmeasure);
     428    }
     429    fclose (fdump);
     430  }
     431
     432  // allocate the output data
     433  ALLOCATE (tmpcatalog.average,  Average,     Nkeep);
     434  ALLOCATE (tmpcatalog.measureT, MeasureTiny, NMEASURE);
     435  ALLOCATE (tmpcatalog.secfilt,  SecFilt,     Nkeep * Nsecfilt);
     436
     437  off_t Nmeasure = 0;
     438
     439  // copy the Nkeep selected entries from subcatalog to tmpcatalog (adjusting links)
     440  for (i = 0; i < Nkeep; i++) {
     441    ave = keepidx[i];
     442    tmpcatalog.average[i] = subcatalog[0].average[ave];
     443    tmpcatalog.average[i].measureOffset = Nmeasure;
     444    for (j = 0; j < tmpcatalog.average[i].Nmeasure; j++) {
     445      off_t offset = subcatalog[0].average[ave].measureOffset + j;
     446      tmpcatalog.measureT[Nmeasure] = subcatalog[0].measureT[offset];
     447      tmpcatalog.measureT[Nmeasure].averef = i;
     448      Nmeasure ++;
     449    }
     450    for (j = 0; j < Nsecfilt; j++) {
     451      tmpcatalog.secfilt[i*Nsecfilt + j] = subcatalog[0].secfilt[ave*Nsecfilt + j];
     452    }
     453  }
     454
     455  if (VERBOSE2) {
     456    char *basename = filebasename (oldcatalog[0].filename);
     457    fprintf (stderr, "limited to %d ("OFF_T_FMT" subset, "OFF_T_FMT" total) stars, "OFF_T_FMT" ("OFF_T_FMT" subset, "OFF_T_FMT" total) measures for catalog %s\n",
     458             Nkeep, subcatalog[0].Naverage, oldcatalog[0].Naverage, Nmeasure, subcatalog[0].Nmeasure,  oldcatalog[0].Nmeasure, basename);
     459    free (basename);
     460  }
     461
     462  free (X);
     463  free (Y);
     464
     465  for (ix = 0; ix < NX; ix++) {
     466    for (iy = 0; iy < NY; iy++) {
     467      free (Nm_grid[ix][iy]);
     468      free (idxgrid[ix][iy]);
     469    }
     470    free (NN_grid[ix]);
     471    free (Nn_grid[ix]);
     472    free (Nm_grid[ix]);
     473    free (idxgrid[ix]);
     474  }
     475  free (NN_grid);
     476  free (Nn_grid);
     477  free (Nm_grid);
     478  free (idxgrid);
     479
     480  free (keepidx);
     481
     482  free (subcatalog[0].average);
     483  free (subcatalog[0].measureT);
     484  free (subcatalog[0].secfilt);
     485
     486  subcatalog[0].average = tmpcatalog.average;
     487  subcatalog[0].measureT = tmpcatalog.measureT;
     488  subcatalog[0].secfilt = tmpcatalog.secfilt;
     489  subcatalog[0].Naverage = Nkeep;
     490  subcatalog[0].Nmeasure = Nmeasure;
     491  subcatalog[0].Nsecfilt = oldcatalog[0].Nsecfilt;
     492  subcatalog[0].Nsecf_mem = Naverage * oldcatalog[0].Nsecfilt;
     493
     494  return (TRUE);
    157495}
    158496
     
    205543    NMEASURE += subcatalog[0].average[ave].Nmeasure;
    206544  }
     545
     546# if (0)
     547  if (oldcatalog[0].catID == 59962) {
     548    FILE *fdump = fopen ("cat.dump.dat", "w");
     549    for (i = 0; i < Nmax; i++) {
     550      ave = index[i];
     551      fprintf (fdump, "%10.6f %10.6f %d\n", subcatalog[0].average[ave].R, subcatalog[0].average[ave].D, subcatalog[0].average[ave].Nmeasure);
     552    }
     553    fclose (fdump);
     554  }
     555# endif
    207556
    208557  // allocate the output data
  • trunk/Ohana/src/relastro/src/dvo_astrom_ops.c

    r33652 r36833  
    9191}
    9292
     93double getMeanR_Big (Measure *measure, Average *average, SecFilt *secfilt) {
     94
     95  double ra;
     96
     97  /* the measure carries the instantaneous mean position at the epoch t */
     98  ra = average[0].R - measure[0].dR / 3600.0;
     99
     100  /* possible corrections to mean ra:
     101
     102  - proper-motion and parallax
     103  - abberation
     104  - precession and nutation, etc
     105  - refraction
     106  - DCR
     107
     108  */
     109
     110  return (ra);
     111}
     112
     113double getMeanD_Big (Measure *measure, Average *average, SecFilt *secfilt) {
     114
     115  double dec;
     116
     117  /* the measure carries the instantaneous mean position at the epoch t */
     118  dec = average[0].D - measure[0].dD / 3600.0;
     119
     120  /* possible corrections to mean ra:
     121
     122  - proper-motion and parallax
     123  - abberation
     124  - precession and nutation, etc
     125  - refraction
     126  - DCR
     127
     128  */
     129
     130  return (dec);
     131}
     132
    93133int setMeanR_Big (double ra_fit, Measure *measure, Average *average, SecFilt *secfilt) {
    94134
     
    144184  return (TRUE);
    145185}
     186
  • trunk/Ohana/src/relastro/src/initialize.c

    r36630 r36833  
    88  if (RELASTRO_OP == OP_MERGE_SOURCE) return;
    99
    10   fprintf (stderr, "PHOTCODE_KEEP_LIST: %s\n", PHOTCODE_KEEP_LIST);
    11   fprintf (stderr, "PHOTCODE_SKIP_LIST: %s\n", PHOTCODE_SKIP_LIST);
    12   fprintf (stderr, "PHOTCODE_A_LIST: %s\n", PHOTCODE_A_LIST);
    13   fprintf (stderr, "PHOTCODE_B_LIST: %s\n", PHOTCODE_B_LIST);
     10  if (DCR_COLOR_POS)  fprintf (stderr, "DCR_COLOR_POS:  %s\n", DCR_COLOR_POS);
     11  if (DCR_COLOR_NEG)  fprintf (stderr, "DCR_COLOR_NEG:  %s\n", DCR_COLOR_NEG);
    1412
    15   photcodesKeep = ParsePhotcodeList (PHOTCODE_KEEP_LIST, &NphotcodesKeep, FALSE);
    16   photcodesSkip = ParsePhotcodeList (PHOTCODE_SKIP_LIST, &NphotcodesSkip, FALSE);
    17   photcodesGroupA = ParsePhotcodeList (PHOTCODE_A_LIST, &NphotcodesGroupA, TRUE);
    18   photcodesGroupB = ParsePhotcodeList (PHOTCODE_B_LIST, &NphotcodesGroupB, TRUE);
     13  if (PHOTCODE_KEEP_LIST)  fprintf (stderr, "PHOTCODE_KEEP_LIST:  %s\n", PHOTCODE_KEEP_LIST);
     14  if (PHOTCODE_SKIP_LIST)  fprintf (stderr, "PHOTCODE_SKIP_LIST:  %s\n", PHOTCODE_SKIP_LIST);
     15  if (PHOTCODE_RESET_LIST) fprintf (stderr, "PHOTCODE_RESET_LIST: %s\n", PHOTCODE_RESET_LIST);
     16  if (PHOTCODE_A_LIST)     fprintf (stderr, "PHOTCODE_A_LIST:     %s\n", PHOTCODE_A_LIST);
     17  if (PHOTCODE_B_LIST)     fprintf (stderr, "PHOTCODE_B_LIST:     %s\n", PHOTCODE_B_LIST);
     18
     19  photcodesKeep   = ParsePhotcodeList (PHOTCODE_KEEP_LIST,  &NphotcodesKeep,   FALSE);
     20  photcodesSkip   = ParsePhotcodeList (PHOTCODE_SKIP_LIST,  &NphotcodesSkip,   FALSE);
     21  photcodesReset  = ParsePhotcodeList (PHOTCODE_RESET_LIST, &NphotcodesReset,  FALSE);
     22  photcodesGroupA = ParsePhotcodeList (PHOTCODE_A_LIST,     &NphotcodesGroupA, TRUE);
     23  photcodesGroupB = ParsePhotcodeList (PHOTCODE_B_LIST,     &NphotcodesGroupB, TRUE);
     24
     25  DCR_NSEC_POS = DCR_NSEC_NEG = -1;
     26  if (DCR_COLOR_POS) {
     27    DCR_PHOTCODE_POS = GetPhotcodebyName (DCR_COLOR_POS);
     28    if (!DCR_PHOTCODE_POS) {
     29      fprintf (stderr, "ERROR: photcode %s not found in photcode table\n", DCR_COLOR_POS);
     30      exit (1);
     31    }
     32    DCR_NSEC_POS = GetPhotcodeNsec (DCR_PHOTCODE_POS[0].code);
     33  }
     34  if (DCR_COLOR_NEG) {
     35    DCR_PHOTCODE_NEG = GetPhotcodebyName (DCR_COLOR_NEG);
     36    if (!DCR_PHOTCODE_NEG) {
     37      fprintf (stderr, "ERROR: photcode %s not found in photcode table\n", DCR_COLOR_NEG);
     38      exit (1);
     39    }
     40    DCR_NSEC_NEG = GetPhotcodeNsec (DCR_PHOTCODE_NEG[0].code);
     41  }
    1942
    2043  initstats (STATMODE);
     
    4265  args_client (argc, argv);
    4366
    44   fprintf (stderr, "PHOTCODE_KEEP_LIST: %s\n", PHOTCODE_KEEP_LIST);
    45   fprintf (stderr, "PHOTCODE_SKIP_LIST: %s\n", PHOTCODE_SKIP_LIST);
    46   fprintf (stderr, "PHOTCODE_A_LIST: %s\n", PHOTCODE_A_LIST);
    47   fprintf (stderr, "PHOTCODE_B_LIST: %s\n", PHOTCODE_B_LIST);
     67  if (PHOTCODE_KEEP_LIST)  fprintf (stderr, "PHOTCODE_KEEP_LIST:  %s\n", PHOTCODE_KEEP_LIST);
     68  if (PHOTCODE_SKIP_LIST)  fprintf (stderr, "PHOTCODE_SKIP_LIST:  %s\n", PHOTCODE_SKIP_LIST);
     69  if (PHOTCODE_RESET_LIST) fprintf (stderr, "PHOTCODE_RESET_LIST: %s\n", PHOTCODE_RESET_LIST);
     70  if (PHOTCODE_A_LIST)     fprintf (stderr, "PHOTCODE_A_LIST:     %s\n", PHOTCODE_A_LIST);
     71  if (PHOTCODE_B_LIST)     fprintf (stderr, "PHOTCODE_B_LIST:     %s\n", PHOTCODE_B_LIST);
    4872
    49   photcodesKeep = ParsePhotcodeList (PHOTCODE_KEEP_LIST, &NphotcodesKeep, FALSE);
    50   photcodesSkip = ParsePhotcodeList (PHOTCODE_SKIP_LIST, &NphotcodesSkip, FALSE);
    51   photcodesGroupA = ParsePhotcodeList (PHOTCODE_A_LIST, &NphotcodesGroupA, TRUE);
    52   photcodesGroupB = ParsePhotcodeList (PHOTCODE_B_LIST, &NphotcodesGroupB, TRUE);
     73  photcodesKeep   = ParsePhotcodeList (PHOTCODE_KEEP_LIST,  &NphotcodesKeep,   FALSE);
     74  photcodesSkip   = ParsePhotcodeList (PHOTCODE_SKIP_LIST,  &NphotcodesSkip,   FALSE);
     75  photcodesReset  = ParsePhotcodeList (PHOTCODE_RESET_LIST, &NphotcodesReset,  FALSE);
     76  photcodesGroupA = ParsePhotcodeList (PHOTCODE_A_LIST,     &NphotcodesGroupA, TRUE);
     77  photcodesGroupB = ParsePhotcodeList (PHOTCODE_B_LIST,     &NphotcodesGroupB, TRUE);
    5378
    5479  initstats (STATMODE);
  • trunk/Ohana/src/relastro/src/launch_region_hosts.c

    r36630 r36833  
    11# include "relastro.h"
    22# define DEBUG 0
    3 
    4 int strextend (char *input, char *format,...) {
    5 
    6   char tmpextra[1024], tmpline[1024];
    7   va_list argp;
    8 
    9   va_start (argp, format);
    10   vsnprintf (tmpextra, 1024, format, argp);
    11   snprintf (tmpline, 1024, "%s %s", input, tmpextra);
    12   strcpy (input, tmpline);
    13 
    14   return TRUE;
    15 }
    163
    174int launch_region_hosts (RegionHostTable *regionHosts) {
     
    8875    if (PHOTCODE_KEEP_LIST) strextend (command, "+photcode %s", PHOTCODE_KEEP_LIST);
    8976    if (PHOTCODE_SKIP_LIST) strextend (command, "-photcode %s", PHOTCODE_SKIP_LIST);
     77    if (PHOTCODE_RESET_LIST) strextend (command, "-reset-to-photcode %s", PHOTCODE_RESET_LIST);
    9078
    9179    if (MaxDensityUse)      strextend (command, "-max-density %f", MaxDensityValue);
    9280    if (ImagSelect)         strextend (command, "-instmag %f %f", ImagMin, ImagMax);
     81    if (ExcludeBogus)       strextend (command, "-exclude-bogus %f", ExcludeBogusRadius);
     82
     83    if (DCR_COLOR_POS && DCR_COLOR_NEG) {
     84      strextend (command, "-dcr-color %s %s", DCR_COLOR_POS, DCR_COLOR_NEG);
     85    }
     86
    9387    if (PhotFlagSelect)     strextend (command, "+photflags");
    9488    if (PhotFlagBad)        strextend (command, "+photflagbad %d", PhotFlagBad);
     
    9791    if (MinBadQF > 0.0)        strextend (command, "-min-bad-psfqf %f", MinBadQF);
    9892    if (MaxMeanOffset != 10.0) strextend (command, "-max-mean-offset  %f", MaxMeanOffset);
     93
     94    strextend (command, "-D RELASTRO_SRC_MEAS_TOOFEW %d", SRC_MEAS_TOOFEW);
     95    strextend (command, "-D RELASTRO_SIGMA_LIM %f", SIGMA_LIM);
     96    strextend (command, "-D RELASTRO_DPOS_MAX %f", DPOS_MAX);
     97    strextend (command, "-D ADDSTAR_RADIUS %f", ADDSTAR_RADIUS);
     98
     99    if (USE_FIXED_PIXCOORDS) strextend (command, "-D USE_FIXED_PIXCOORDS 1");
    99100
    100101    if (TimeSelect) {
  • trunk/Ohana/src/relastro/src/load_catalogs.c

    r36630 r36833  
    11# include "relastro.h"
    22
    3 Catalog *load_catalogs_parallel (SkyList *sky, int *Ncatalog);
     3Catalog *load_catalogs_parallel (SkyList *sky, int *Ncatalog, char *syncfile);
    44void bcatalog_show_skips ();
    55
    6 Catalog *load_catalogs (SkyList *skylist, int *Ncatalog, int subselect, int hostID, char *hostpath) {
    7 
    8   int i, j, k, m, Nstar;
     6Catalog *load_catalogs (SkyList *skylist, int *Ncatalog, int subselect, int hostID, char *hostpath, char *syncfile) {
     7
     8  int i, j, Nstar;
     9  // int k, m;
    910  Catalog *catalog, *pcatalog, tcatalog;
    1011
    1112  // XXX need to decide how to determine PARALLEL mode...
    1213  if (PARALLEL && !hostID) {
    13     catalog = load_catalogs_parallel (skylist, Ncatalog);
     14    catalog = load_catalogs_parallel (skylist, Ncatalog, syncfile);
    1415    return catalog;
    1516  }
     
    7374      if (RESET) {
    7475        for (j = 0; j < catalog[Ncat].Naverage; j++) {
     76# if (0)         
    7577          catalog[Ncat].average[j].flags = 0;
    7678          m = catalog[Ncat].average[j].measureOffset;
     
    7880            catalog[Ncat].measure[m+k].dbFlags = 0;
    7981          }
     82# endif
    8083        }
    8184      }
     
    8487  }
    8588
    86   // XXX TEST : bcatalog_show_skips();
     89  bcatalog_show_skips();
    8790
    8891  Nstar = 0;
     
    113116// CATDIR is supplied globally
    114117# define DEBUG 1
    115 Catalog *load_catalogs_parallel (SkyList *sky, int *Ncatalog) {
     118Catalog *load_catalogs_parallel (SkyList *sky, int *Ncatalog, char *syncfile) {
    116119
    117120  char uniquer[12];
     
    150153              table->hosts[i].results, table->hosts[i].hostID, CATDIR, table->hosts[i].pathname, UserPatch.Rmin, UserPatch.Rmax, UserPatch.Dmin, UserPatch.Dmax, STATMODE, MIN_ERROR, SIGMA_LIM);
    151154
    152     char tmpline[1024];
    153     if (FIT_MODE == FIT_PM_ONLY)         { snprintf (tmpline, 1024, "%s -pm",    command);           strcpy (command, tmpline); }
    154     if (FIT_MODE == FIT_PAR_ONLY)        { snprintf (tmpline, 1024, "%s -par",   command);           strcpy (command, tmpline); }
    155     if (FIT_MODE == FIT_PM_AND_PAR)      { snprintf (tmpline, 1024, "%s -pmpar", command);           strcpy (command, tmpline); }
    156 
    157     if (VERBOSE)       { snprintf (tmpline, 1024, "%s -v",              command);                    strcpy (command, tmpline); }
    158     if (VERBOSE2)      { snprintf (tmpline, 1024, "%s -vv",             command);                    strcpy (command, tmpline); }
    159     if (RESET)         { snprintf (tmpline, 1024, "%s -reset",          command);                    strcpy (command, tmpline); }
    160     if (ImagSelect)    { snprintf (tmpline, 1024, "%s -instmag %f %f",  command, ImagMin, ImagMax);  strcpy (command, tmpline); }
    161     if (MaxDensityUse) { snprintf (tmpline, 1024, "%s -max-density %f", command, MaxDensityValue);   strcpy (command, tmpline); }
     155    if (FIT_MODE == FIT_PM_ONLY)         strextend (command, "-pm");
     156    if (FIT_MODE == FIT_PAR_ONLY)        strextend (command, "-par");
     157    if (FIT_MODE == FIT_PM_AND_PAR)      strextend (command, "-pmpar");
     158
     159    if (VERBOSE)       strextend (command, "-v");
     160    if (VERBOSE2)      strextend (command, "-vv");
     161    if (RESET)         strextend (command, "-reset");
     162    if (ImagSelect)    strextend (command, "-instmag %f %f", ImagMin, ImagMax);
     163    if (MaxDensityUse) strextend (command, "-max-density %f", MaxDensityValue);
     164    if (FlagOutlier)   strextend (command, "-clip %d", CLIP_THRESH);
     165    if (ExcludeBogus)  strextend (command, "-exclude-bogus %f", ExcludeBogusRadius);
    162166   
    163     if (FlagOutlier)   { snprintf (tmpline, 1024, "%s -clip %d",        command, CLIP_THRESH);       strcpy (command, tmpline); }
    164    
    165     if (USE_FIXED_PIXCOORDS) { snprintf (tmpline, 1024, "%s -D USE_FIXED_PIXCOORDS 1", command);     strcpy (command, tmpline); }
    166 
    167     if (PHOTCODE_KEEP_LIST) { snprintf (tmpline, 1024, "%s +photcode %s", command, PHOTCODE_KEEP_LIST); strcpy (command, tmpline); }
    168     if (PHOTCODE_SKIP_LIST) { snprintf (tmpline, 1024, "%s -photcode %s", command, PHOTCODE_SKIP_LIST); strcpy (command, tmpline); }
    169     if (PhotFlagSelect)    { snprintf (tmpline, 1024, "%s +photflags",   command);                     strcpy (command, tmpline); }
    170     if (PhotFlagBad)       { snprintf (tmpline, 1024, "%s +photflagbad %d", command, PhotFlagBad);     strcpy (command, tmpline); }
    171     if (PhotFlagPoor)      { snprintf (tmpline, 1024, "%s +photflagpoor %d", command, PhotFlagPoor);   strcpy (command, tmpline); }
     167    if (DCR_COLOR_POS && DCR_COLOR_NEG) {
     168      strextend (command, "-dcr-color %s %s", DCR_COLOR_POS, DCR_COLOR_NEG);
     169    }
     170
     171    if (USE_FIXED_PIXCOORDS) strextend (command, "-D USE_FIXED_PIXCOORDS 1");
     172    if (PHOTCODE_KEEP_LIST)  strextend (command, "+photcode %s", PHOTCODE_KEEP_LIST);
     173    if (PHOTCODE_SKIP_LIST)  strextend (command, "-photcode %s", PHOTCODE_SKIP_LIST);
     174    if (PhotFlagSelect)      strextend (command, "+photflags");
     175    if (PhotFlagBad)         strextend (command, "+photflagbad %d", PhotFlagBad);
     176    if (PhotFlagPoor)        strextend (command, "+photflagpoor %d", PhotFlagPoor);
    172177    // XXX note that the above pass in the flag as decimal -- also note that args.c cannot handle 0xHEX values
    173178
     
    175180      char *tstart = ohana_sec_to_date (TSTART);
    176181      char *tstop  = ohana_sec_to_date (TSTOP);
    177       snprintf (tmpline, 1024, "%s -time %s %s", command, tstart, tstop);
     182      strextend (command, "-time %s %s", tstart, tstop);
    178183      free (tstart);
    179184      free (tstop);
    180       strcpy (command, tmpline);
    181185    }
    182186
     
    211215  }
    212216
     217  // update syncfile here (save lots of I/O time)
     218
     219  // NOTE: if I let all hosts load blindly, I saturate the data clients with too many
     220  // relastro_client requests.  I need to have the master mediate this.  the master
     221  // will not launch the next remote job until this one says it is done
     222  if (syncfile) {
     223    update_sync_file (syncfile, 1);
     224  }
     225
    213226  // each host generates a BrightCatalog structure, with the measure, average, etc value
    214227  // loaded into a single set of arrays (of MeasureTiny, AverageTiny, Secfilt).  I need to
  • trunk/Ohana/src/relastro/src/relastro_client.c

    r36630 r36833  
    3333      // USAGE: relastro_client -load-objects
    3434      int Ncatalog;
    35       Catalog *catalog = load_catalogs (skylist, &Ncatalog, TRUE, HOST_ID, HOSTDIR);
     35      Catalog *catalog = load_catalogs (skylist, &Ncatalog, TRUE, HOST_ID, HOSTDIR, NULL);
    3636      if (!catalog) {
    3737        fprintf (stderr, "ERROR loading catalogs from %s\n", CATDIR);
  • trunk/Ohana/src/relastro/src/relastro_images.c

    r36630 r36833  
    2828  // XXX pass in the image table
    2929  // XXX who carries the image grid?
    30   catalog = load_catalogs (skylist, &Ncatalog, TRUE, 0, NULL);
     30
     31  // photcodesKeep is used here to allow measurements from the images being calibrated
     32  // note if -reset-to-photcode is selected, photocodesKeep is replaced with below with photcodesReset
     33  catalog = load_catalogs (skylist, &Ncatalog, TRUE, 0, NULL, NULL);
    3134  MARKTIME("load catalog data: %f sec\n", dtime);
     35
     36  if (photcodesReset) {
     37    photcodesKeep  = photcodesReset;
     38    NphotcodesKeep = NphotcodesReset;
     39  }
    3240
    3341  if (Ncatalog == 0) {
     
    5058  // set test points based on the starmap
    5159  createStarMap (catalog, Ncatalog);
     60
     61  // XXX NOTE : for 2mass reset, photcodesKeep should now limit to 2MASS measurements
    5262
    5363  /* major modes */
  • trunk/Ohana/src/relastro/src/relastro_objects.c

    r36630 r36833  
    4040    snprintf (hostfile, 1024, "%s/%s.cpt", hostpath, skylist[0].regions[i]->name);
    4141    catalog.filename  = hostID ? hostfile : skylist[0].filename[i];
     42
     43    // set up the basic catalog info
    4244    catalog.catformat = dvo_catalog_catformat (CATFORMAT);    // set the default catformat from config data
    4345    catalog.catmode   = dvo_catalog_catmode (CATMODE);        // set the default catmode from config data
  • trunk/Ohana/src/relastro/src/relastro_parallel_images.c

    r36630 r36833  
    4545
    4646  /* load catalog data from region files (hostID is 0 since we are not a client */
    47   catalog = load_catalogs (skylist, &Ncatalog, TRUE, 0, NULL);
     47  char *syncfile = make_filename (CATDIR, regionHosts->hosts[myHost].hostname, REGION_HOST_ID, "loadcat.sync");
     48  catalog = load_catalogs (skylist, &Ncatalog, TRUE, 0, NULL, syncfile);
    4849  MARKTIME("-- load catalog data: %f sec\n", dtime);
     50  free (syncfile);
    4951 
    50   // NOTE: if I let all hosts load blindly, I saturate the data clients with too many
    51   // relastro_client requests.  I need to have the master mediate this.  the master
    52   // will not launch the next remote job until this one says it is done
    53   char *syncfile = make_filename (CATDIR, regionHosts->hosts[myHost].hostname, REGION_HOST_ID, "loadcat.sync");
    54   update_sync_file (syncfile, 1);
    55  
     52  if (photcodesReset) {
     53    photcodesKeep  = photcodesReset;
     54    NphotcodesKeep = NphotcodesReset;
     55  }
     56
    5657  // generate tables go from catID,objID -> catSeq,objSeq
    5758  indexCatalogs (catalog, Ncatalog);
  • trunk/Ohana/src/uniphot/src/update_catalog_fiximids.c

    r35764 r36833  
    11# include "fiximids.h"
     2
     3// test image: 2013/06/15,13:25:51, GPC1.r.XY50
     4static int CHECK_TEST_IMAGE = FALSE;
     5static unsigned int Tref = 1378812312;
     6static short Cref = 10001;
    27
    38void update_catalog_fiximids (Catalog *catalog) {
     
    1823    short photcode = catalog[0].measure[i].photcode;
    1924    e_time time = catalog[0].measure[i].t;
     25
     26    if (CHECK_TEST_IMAGE && (abs(time - Tref) < 10) && (photcode == Cref)) {
     27      fprintf (stderr, ".");
     28    }
    2029
    2130    // skip detections with no valid imageID (eg, ref photcode)
Note: See TracChangeset for help on using the changeset viewer.