IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 35287


Ignore:
Timestamp:
Mar 9, 2013, 3:00:20 PM (13 years ago)
Author:
eugene
Message:

adding hpm search

Location:
branches/eam_branches/ipp-20130307/Ohana/src/relastro
Files:
1 added
6 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ipp-20130307/Ohana/src/relastro/Makefile

    r33652 r35287  
    6464$(SRC)/high_speed_objects.$(ARCH).o  \
    6565$(SRC)/high_speed_utils.$(ARCH).o  \
     66$(SRC)/hpm_catalogs.$(ARCH).o  \
     67$(SRC)/hpm_objects.$(ARCH).o  \
    6668$(SRC)/relastro_merge_source.$(ARCH).o  \
    6769$(SRC)/resort_catalog.$(ARCH).o  \
    … …  
    106108$(SRC)/high_speed_objects.$(ARCH).o  \
    107109$(SRC)/high_speed_utils.$(ARCH).o  \
     110$(SRC)/hpm_catalogs.$(ARCH).o  \
     111$(SRC)/hpm_objects.$(ARCH).o  \
    108112$(SRC)/relastro_merge_source.$(ARCH).o  \
    109113$(SRC)/resort_catalog.$(ARCH).o  \

  • branches/eam_branches/ipp-20130307/Ohana/src/relastro/include/relastro.h

    r35105 r35287  
    2020typedef enum {FIT_NONE, FIT_AVERAGE, FIT_PM_ONLY, FIT_PAR_ONLY, FIT_PM_AND_PAR} FitMode;
    2121
    22 typedef enum {TARGET_NONE, TARGET_SIMPLE, TARGET_CHIPS, TARGET_MOSAICS, TARGET_HIGH_SPEED, TARGET_MERGE_SOURCE, TARGET_UPDATE_OBJECTS, TARGET_UPDATE_OFFSETS, TARGET_LOAD_OBJECTS} FitTarget;
     22typedef enum {TARGET_NONE, TARGET_SIMPLE, TARGET_CHIPS, TARGET_MOSAICS, TARGET_HIGH_SPEED, TARGET_MERGE_SOURCE, TARGET_UPDATE_OBJECTS, TARGET_UPDATE_OFFSETS, TARGET_LOAD_OBJECTS, TARGET_HPM} FitTarget;
    2323
    2424typedef enum {
    … …  
    442442
    443443PhotCode **ParsePhotcodeList (char *rawlist, int *nphotcodes, int needAve);
     444
     445int hpm_catalogs (SkyTable *sky, SkyList *skylist, int hostID, char *hostpath);
     446int hpm_catalogs_parallel (SkyList *skylist);
     447int hpm_objects (SkyRegion *region, Catalog *catalog);

  • branches/eam_branches/ipp-20130307/Ohana/src/relastro/src/args.c

    r34749 r35287  
    8989  if ((N = get_argument (argc, argv, "-high-speed"))) {
    9090    // XXX include a parallax / no-parallax option
     91    if (N >= argc - 5) usage();
     92    FIT_TARGET = TARGET_HIGH_SPEED;
     93    remove_argument (N, &argc, argv);
     94    PHOTCODE_A_LIST = strcreate(argv[N]);
     95    remove_argument (N, &argc, argv);
     96    PHOTCODE_B_LIST = strcreate(argv[N]);
     97    remove_argument (N, &argc, argv);
     98    RADIUS = atof(argv[N]);
     99    remove_argument (N, &argc, argv);
     100    HIGH_SPEED_DIR = abspath(argv[N], DVO_MAX_PATH);
     101    remove_argument (N, &argc, argv);
     102  }
     103
     104  if ((N = get_argument (argc, argv, "-hpm"))) {
    91105    if (N >= argc - 3) usage();
    92     FIT_TARGET = TARGET_HIGH_SPEED;
    93     remove_argument (N, &argc, argv);
    94     PHOTCODE_A_LIST = strcreate(argv[N]);
    95     remove_argument (N, &argc, argv);
    96     PHOTCODE_B_LIST = strcreate(argv[N]);
     106    FIT_TARGET = TARGET_HPM;
    97107    remove_argument (N, &argc, argv);
    98108    RADIUS = atof(argv[N]);
    … …  
    106116    remove_argument (N, &argc, argv);
    107117    PARALLEL_OUTPUT = TRUE;
    108     if (FIT_TARGET != TARGET_HIGH_SPEED) {
    109       fprintf (stderr, "-parallel-output only valid for -high-speed mode\n");
     118    if ((FIT_TARGET != TARGET_HIGH_SPEED) && (FIT_TARGET != TARGET_HPM)) {
     119      fprintf (stderr, "-parallel-output only valid for -high-speed or -hpm modes\n");
    110120      exit (1);
    111121    }
    … …  
    456466  if ((N = get_argument (argc, argv, "-high-speed"))) {
    457467    // XXX include a parallax / no-parallax option
    458     if (N >= argc - 3) usage_client();
     468    if (N >= argc - 5) usage_client();
    459469    FIT_TARGET = TARGET_HIGH_SPEED;
    460470    remove_argument (N, &argc, argv);
    … …  
    466476    remove_argument (N, &argc, argv);
    467477    HIGH_SPEED_DIR = strcreate(argv[N]);
     478    remove_argument (N, &argc, argv);
     479  }
     480
     481  if ((N = get_argument (argc, argv, "-hpm"))) {
     482    if (N >= argc - 3) usage();
     483    FIT_TARGET = TARGET_HPM;
     484    remove_argument (N, &argc, argv);
     485    RADIUS = atof(argv[N]);
     486    remove_argument (N, &argc, argv);
     487    HIGH_SPEED_DIR = abspath(argv[N], DVO_MAX_PATH);
    468488    remove_argument (N, &argc, argv);
    469489  }
    … …  
    668688  fprintf (stderr, "       OR:    relastro -update-objects [options]\n");
    669689  fprintf (stderr, "       OR:    relastro -high-speed [options]\n");
     690  fprintf (stderr, "       OR:    relastro -hpm [options]\n");
    670691  fprintf (stderr, "       OR:    relastro -merge-source [options]\n\n");
    671692
    … …  
    679700  fprintf (stderr, "  -update-mosaics\n");
    680701  fprintf (stderr, "  -high-speed (code[,code,code]) (code[,code,code]) (radius) (output catdir)\n");
     702  fprintf (stderr, "  -hpm (radius) (output catdir)\n");
    681703  fprintf (stderr, "  -merge-source (objID) (catID) into (objID) (catID)\n\n");
    682704

  • branches/eam_branches/ipp-20130307/Ohana/src/relastro/src/hpm_objects.c

    r35286 r35287  
    11# include "relastro.h"
     2# define MAX_TRANGE 180.0
     3# define MIN_PS1_DET 3
     4# define MAX_CHISQ_PM 10.0
     5# define MIN_NPOS_FRAC 0.8
    26
    37# define NEXT_I { if (Ngroup < 2) slowMoving[ni] = TRUE; newI = TRUE; i++; continue; }
    … …  
    610int hpm_objects (SkyRegion *region, Catalog *catalog) {
    711
    8   off_t i, j, m, J, ni, nj, *N1;
    9   off_t Nslow, Ninvalid, NgroupA, NgroupB, NgroupAbad, NgroupBbad;
    10   off_t l, i1, NAVERAGE, NMEASURE;
    11   int *slowMoving, *groupA, *groupB, status, Nmatch, Nmatchmeas, Nepoch, nv[2], Nmatchmeasobj;
    12   int foundA, foundB, XVERB;
     12  off_t i, j, J, ni, nj, *N1, k, m;
     13  int *slowMoving;
    1314  double *X1, *Y1;
    1415  double dX, dY, dR, RADIUS2;
    1516  Coords tcoords;
    16   Catalog catalogOut;
    17   Catalog testcat;
     17  Catalog catalogOut, testcat;
    1818
    1919  int Nsecfilt;
    2020  char filename[1024];
    2121
    22   // XXX are we saving these in an hpm dvodb?
     22  // we need at least 2 objects if we are going to match anything...
     23  if (catalog[0].Naverage < 2) return (TRUE);
     24  if (VERBOSE) fprintf (stderr, "checking "OFF_T_FMT" objects\n",  catalog[0].Naverage);
     25
     26  // we save the best objects in an hpm dvodb
    2327  snprintf (filename, 1024, "%s/%s.cpt", HIGH_SPEED_DIR, region[0].name);
    24   fprintf (stderr, "%s\n",filename);
    25 
    26   dvo_catalog_init (&catalogOut, TRUE); /*initialise new catalogue*/
     28
     29  dvo_catalog_init (&catalogOut, TRUE); /* init new catalog */
    2730  catalogOut.filename = strcreate(filename);
    2831
    29   SIGMA_LIM = 0.0;
    30 
    3132  Nsecfilt = GetPhotcodeNsecfilt();
    32   // off_t Naverage = catalog[0].Naverage;
    33   // off_t Nmeasure = catalog[0].Nmeasure;
    34   catalogOut.filename = filename; // based on the input name, need to keep everything below the catdir portion
     33  catalogOut.filename = filename;
    3534  catalogOut.Nsecfilt = Nsecfilt;
    3635  catalogOut.catflags = LOAD_AVES | LOAD_MEAS | LOAD_MISS | LOAD_SECF; // load all data
    … …  
    4544  }
    4645
    47   NAVERAGE = 1000;
    48   NMEASURE = 10000;
     46  off_t NAVERAGE = 1000;
     47  off_t NMEASURE = 10000;
    4948  REALLOCATE (catalogOut.average, Average, NAVERAGE);
    5049  REALLOCATE (catalogOut.measure, Measure, NMEASURE);
    5150  REALLOCATE (catalogOut.secfilt, SecFilt, NAVERAGE*Nsecfilt);
     51  catalogOut.Naverage = 0;
     52  catalogOut.Nmeasure = 0;
    5253
    5354  // testcat is used to determine the fit for a single object group
    … …  
    5556  dvo_catalog_init (&testcat, TRUE);
    5657  testcat.Naverage = 1; // this is fixed -- only one obj in testcat
    57   NMEASURE_TEST = 1000;
     58  off_t NMEASURE_TESTCAT = 1000;
    5859  REALLOCATE (catalogOut.average, Average, 1);
    59   REALLOCATE (catalogOut.measure, Measure, NMEASURE_TEST);
     60  REALLOCATE (catalogOut.measure, Measure, NMEASURE_TESTCAT);
    6061  REALLOCATE (catalogOut.secfilt, SecFilt, Nsecfilt);
    61 
    62   // we need at least 2 objects if we are going to match anything...
    63   if (catalog[0].Naverage < 2) return (TRUE);
    6462
    6563  // mask with which to mark objects to be ignored
    6664  ALLOCATE (slowMoving, int, catalog[0].Naverage);
    6765  memset (slowMoving, 0, catalog[0].Naverage*sizeof(int));
    68 
    69   if (VERBOSE) fprintf (stderr, "checking "OFF_T_FMT" objects\n",  catalog[0].Naverage);
    70   Nslow = 0;
    71   Ninvalid = 0;
    72   NgroupAbad = 0;
    73   NgroupBbad = 0;
    74 
    75   // mark (exclude) objects with both sets of target photcodes
     66  off_t Nslow = 0;
     67
     68  // mark (exclude) objects which do not meet the selection criterion
    7669  for (i = 0; i < catalog[0].Naverage; i++) {
    7770
    78     XVERB  = (catalog[0].average[i].objID == OBJ_ID_SRC) && (catalog[0].average[i].catID == CAT_ID_SRC);
     71    int XVERB = FALSE;
     72    XVERB |= (catalog[0].average[i].objID == OBJ_ID_SRC) && (catalog[0].average[i].catID == CAT_ID_SRC);
    7973    XVERB |= (catalog[0].average[i].objID == OBJ_ID_DST) && (catalog[0].average[i].catID == CAT_ID_DST);
    8074    if (XVERB) {
    … …  
    8276    }
    8377
    84     // selection criteria:
    85     // (Nps1 > 2) && (Trange < 180)
     78    // selection criteria : (Nps1 > 2) && (Trange < 180)
    8679    if (catalog[0].average[i].Trange > MAX_TRANGE) {
    8780      slowMoving[i] = TRUE;
    … …  
    129122  /* build spatial index (RA sort) referencing input array sequence */
    130123  for (i = 0; i < catalog[0].Naverage; i++) {
    131     status = RD_to_XY (&X1[i], &Y1[i], catalog[0].average[i].R, catalog[0].average[i].D, &tcoords);
     124    int status = RD_to_XY (&X1[i], &Y1[i], catalog[0].average[i].R, catalog[0].average[i].D, &tcoords);
    132125    N1[i] = i;
    133126    assert (status);
    … …  
    135128  sort_coords_index (X1, Y1, N1, catalog[0].Naverage);
    136129
    137   Nmatch = 0;
    138   Nmatchmeas = 0;
    139130  RADIUS2 = SQ(RADIUS);
    140131
    141132  // group is a list of objects that are within a clump.  this set will be tested
    142133  // via a clipped fit to the measurements after the group is identified
    143   Ngroup = 0;
    144   NGROUP = 100;
     134  int Ngroup = 0;
     135  int NGROUP = 100;
     136  int *group = NULL;
    145137  ALLOCATE (group, int, NGROUP);
    146138
    147139  // in the loop below, we need to do a bunch of things when we go to the next main object
    148   newI = TRUE;
     140  int newI = TRUE;
    149141
    150142  // mark (skip) objects with both sets of target photcodes
    … …  
    159151
    160152    if (newI) {
    161       Nmatch = 0;
    162153      Ngroup = 1;
    163154      group[0] = ni;
    … …  
    165156    }
    166157
    167     XVERB  = (catalog[0].average[ni].objID == OBJ_ID_SRC) && (catalog[0].average[ni].catID == CAT_ID_SRC);
     158    int XVERB = FALSE;
     159    XVERB |= (catalog[0].average[ni].objID == OBJ_ID_SRC) && (catalog[0].average[ni].catID == CAT_ID_SRC);
    168160    XVERB |= (catalog[0].average[nj].objID == OBJ_ID_DST) && (catalog[0].average[nj].catID == CAT_ID_DST);
    169161    if (XVERB) {
    … …  
    186178    if (dX >= +1.02*RADIUS) NEXT_J; // positive dX: j is too small
    187179
    188     // within match range; look for valid matches
     180    // within match range; look for valid matches & accumulate the group
    189181    for (J = j; (dX > -1.02*RADIUS) && (J < catalog[0].Naverage); J++) {     
    190182      if (J == i) continue;  // avoid auto-matches
    … …  
    219211    // the mean object will start with info from the primary object
    220212    // remember: testcat.Naverage = 1 -- does not change
    221     testcat.average[0] = catalog[0].average[group[0]];
     213    ni = group[0];
     214    int Nmatchmeas = 0;
     215    testcat.average[0] = catalog[0].average[ni];
    222216    testcat.average[0].measureOffset = 0;
    223     Nmatchmeas = 0;
    224     I = group[0];
    225217    for (J = 0; J < Ngroup; J++) {
    226       J = group[J];
    227       m = catalog[0].average[Nj].measureOffset;
    228       for (k = 0; k < catalog[0].average[Nj].Nmeasure; k++) {
     218      nj = group[J];
     219      m = catalog[0].average[nj].measureOffset;
     220      for (k = 0; k < catalog[0].average[nj].Nmeasure; k++) {
    229221        testcat.measure[Nmatchmeas] = catalog[0].measure[m+k];
    230222        /* Set offset RA and Dec wrt correct average value*/
    231         testcat.measure[Nmatchmeas].dR = catalog[0].measure[m+k].dR + 3600.0*(catalog[0].average[I].R - catalog[0].average[J].R);
    232         testcat.measure[Nmatchmeas].dD = catalog[0].measure[m+k].dD + 3600.0*(catalog[0].average[I].D - catalog[0].average[J].D);
     223        testcat.measure[Nmatchmeas].dR = catalog[0].measure[m+k].dR + 3600.0*(catalog[0].average[ni].R - catalog[0].average[nj].R);
     224        testcat.measure[Nmatchmeas].dD = catalog[0].measure[m+k].dD + 3600.0*(catalog[0].average[ni].D - catalog[0].average[nj].D);
    233225        testcat.measure[Nmatchmeas].averef = 0;
    234226        Nmatchmeas++;
    235         CHECK_REALLOCATE (testcat.measure, Measure, NMEASURE_TEST, Nmatchmeas, 1000);
     227        CHECK_REALLOCATE (testcat.measure, Measure, NMEASURE_TESTCAT, Nmatchmeas, 1000);
    236228      }
    237229    }
    … …  
    244236    UpdateObjects (&testcat, 1);
    245237
    246     // logic for keeping the fit:
    247     if (testcat.average[0].ChiSqPM < XXX) good = TRUE;
    248     if (testcat.average[0].Npos > 0.X * testcat.average[0].Nmeasure) good = TRUE;
    249     // what else?
     238    // logic for keeping the fit (anything else?)
     239    int good = FALSE;
     240    good |= (testcat.average[0].ChiSqPM < MAX_CHISQ_PM);
     241    good |= (testcat.average[0].Npos > MIN_NPOS_FRAC * testcat.average[0].Nmeasure);
    250242
    251243    if (good) {
    252244      // save the new object on catalogOut
    253       {
    254         catalogOut.average[Nout] = testcat.average[0];
    255         catalogOut.average[Nout].measureOffset = Nmatchmeas;
    256         m = testcat.average[0].measureOffset;
    257         for (k = 0; k < testcat.average[0].Nmeasure; k++) {
    258           catalogOut.measure[Noutmeas] = testcat.measure[m+k];
    259           catalogOut.measure[Noutmeas].averef = Nout;
    260         }
    261         Nout ++;
    262         CHECK_REALLOCATE (catalogOut.average, Average, NAVERAGE, Nout, 100);
     245      // (note: there is only one object in testcat and thus measureOffset = 0
     246      catalogOut.average[catalogOut.Naverage] = testcat.average[0];
     247      catalogOut.average[catalogOut.Naverage].measureOffset = catalogOut.Nmeasure;
     248      for (k = 0; k < testcat.average[0].Nmeasure; k++) {
     249        catalogOut.measure[catalogOut.Nmeasure] = testcat.measure[k];
     250        catalogOut.measure[catalogOut.Nmeasure].averef = catalogOut.Naverage;
     251        catalogOut.Nmeasure ++;
     252        CHECK_REALLOCATE (catalogOut.measure, Measure, NMEASURE, catalogOut.Nmeasure, 100);
    263253      }
     254      catalogOut.Naverage ++;
     255      CHECK_REALLOCATE (catalogOut.average, Average, NAVERAGE, catalogOut.Naverage, 100);
    264256    }
    265257    NEXT_I;
    266258  }
    267 
    268   catalogOut.Naverage=Nmatch;
    269   catalogOut.Nmeasure=Nmatchmeas;
    270   catalogOut.Nsecfilt=Nsecfilt;
    271   catalogOut.Nsecf_mem=Nmatch*Nsecfilt;
     259  catalogOut.Nsecfilt = Nsecfilt;
     260  catalogOut.Nsecf_mem = Nsecfilt * catalogOut.Naverage;
    272261
    273262  populate_tiny_values (&catalogOut, DVO_TV_MEASURE);
    274263  UpdateObjects (&catalogOut, 1);
    275264
    276   fprintf (stderr, "found %d matches\n", Nmatch);
     265  fprintf (stderr, "found "OFF_T_FMT" matches with "OFF_T_FMT" detections\n", catalogOut.Naverage, catalogOut.Nmeasure);
    277266  dvo_catalog_save (&catalogOut, VERBOSE2);
    278267  dvo_catalog_unlock (&catalogOut);
    279268  dvo_catalog_free (&catalogOut);
    280269  free (slowMoving);
    281   free (groupA);
    282   free (groupB);
    283270  free (X1);
    284271  free (Y1);
    … …  
    287274  return TRUE;
    288275}
    289 
    290 // return TRUE if measure->photcode match any in the photcodeSet
    291 // (but, return FALSE if the measurement is bad)
    292 int MeasMatchesPhotcode(Measure *measure, PhotCode **photcodeSet, int Nset) {
    293 
    294   int found, k;
    295 
    296   if (!Nset) return TRUE;
    297   // XXX require a set or not?  assert (Nset > 0);
    298 
    299   if (!finite(measure[0].dR)) return FALSE;
    300   if (!finite(measure[0].dD)) return FALSE;
    301   if (!finite(measure[0].M))  return FALSE;
    302  
    303   float dX = GetAstromError (measure, ERROR_MODE_RA);
    304   if (isnan(dX)) return FALSE;
    305 
    306   float dY = GetAstromError (measure, ERROR_MODE_DEC);
    307   if (isnan(dY)) return FALSE;
    308 
    309   /* select measurements by photcode, or equiv photcode, if specified */
    310   found = FALSE;
    311   for (k = 0; (k < Nset) && !found; k++) {
    312     if (photcodeSet[k][0].code == measure[0].photcode) found = TRUE;
    313     if (photcodeSet[k][0].code == GetPhotcodeEquivCodebyCode(measure[0].photcode)) found = TRUE;
    314   }
    315 
    316   return found;
    317 }

  • branches/eam_branches/ipp-20130307/Ohana/src/relastro/src/relastro.c

    r33652 r35287  
    1919      /* high-speed is a 2pt cross-correlation process for linking moving objects (high PM) */
    2020      high_speed_catalogs (sky, skylist, 0, NULL);
     21      exit (0);
     22
     23    case TARGET_HPM:
     24      hpm_catalogs (sky, skylist, 0, NULL);
    2125      exit (0);
    2226

  • branches/eam_branches/ipp-20130307/Ohana/src/relastro/src/relastro_client.c

    r33652 r35287  
    5858    }
    5959
     60    case TARGET_HPM: {
     61      // USAGE: relastro_client -high-speed
     62      hpm_catalogs (sky, skylist, HOST_ID, HOSTDIR);
     63      break;
     64    }
     65
    6066      // XXX loading the images is fairly costly -- see if we can do an image subset
    6167    case TARGET_UPDATE_OFFSETS: {
Note: See TracChangeset for help on using the changeset viewer.