IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 28855


Ignore:
Timestamp:
Aug 5, 2010, 2:29:00 PM (16 years ago)
Author:
bills
Message:

Enhance dvomerge to handle merging catdirs with different (but compatible) secondary filter lists

Location:
trunk/Ohana/src
Files:
9 edited

Legend:

Unmodified
Added
Removed

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

    r27583 r28855  
    2525float  RADIUS;
    2626int    SKY_DEPTH;
     27char   *ALTERNATE_PHOTCODE_FILE;
    2728
    2829SkyRegion UserPatch;  // used by MODE CAT
    … …  
    6162int        LoadCatalog            PROTO((Catalog *catalog, SkyRegion *region, char *filename, char *mode));
    6263
    63 int        merge_catalogs_new     PROTO((SkyRegion *region, Catalog *output, Catalog *input));
    64 int        merge_catalogs_old     PROTO((SkyRegion *region, Catalog *output, Catalog *input, double RADIUS));
     64int        merge_catalogs_new     PROTO((SkyRegion *region, Catalog *output, Catalog *input, int *secflitMap));
     65int        merge_catalogs_old     PROTO((SkyRegion *region, Catalog *output, Catalog *input, double RADIUS, int *secflitMap));
    6566
    6667off_t     *init_measure_links     PROTO((Average *average, off_t Naverage, Measure *measure, off_t Nmeasure));

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

    r27583 r28855  
    1010  if ((N = get_argument (*argc, argv, "-v"))) {
    1111    VERBOSE = TRUE;
     12    remove_argument (N, argc, argv);
     13  }
     14  if ((N = get_argument (*argc, argv, "-photcode-file"))) {
     15    remove_argument (N, argc, argv);
     16    ALTERNATE_PHOTCODE_FILE = strdup(argv[N]);
    1217    remove_argument (N, argc, argv);
    1318  }

  • trunk/Ohana/src/dvomerge/src/dvoconvert.c

    r28241 r28855  
    8585        continue;
    8686      }
    87       merge_catalogs_new (&outsky[0].regions[i], &outcatalog, &incatalog);
     87      merge_catalogs_new (&outsky[0].regions[i], &outcatalog, &incatalog, NULL);
    8888      dvo_catalog_unlock (&incatalog);
    8989      dvo_catalog_free (&incatalog);

  • trunk/Ohana/src/dvomerge/src/dvomergeCreate.c

    r27583 r28855  
    1010  SkyList *inlist;
    1111  Catalog incatalog, outcatalog;
    12   char filename[256], *input1, *input2, *output;
     12  char filename[256];
     13  char filename1[256];
     14  char *codesFilename;
     15  char *input1, *input2, *output;
    1316  IDmapType IDmap1, IDmap2;
    1417
    … …  
    2225  dvomergeImagesCreate (&IDmap1, input1, &IDmap2, input2, output);
    2326
    24   // the first input defines the photcode table & db layout
    25   sprintf (filename, "%s/Photcodes.dat", input1);
     27  PhotCodeData *input1Photcodes = NULL;
     28  PhotCodeData *input2Photcodes = NULL;
     29  PhotCodeData *outputPhotcodes = NULL;
     30  int *secfiltMap1 = NULL;
     31  int *secfiltMap2 = NULL;
     32
     33  // Read the input1 photcodes
     34  SetPhotcodeTable(NULL);
     35  sprintf (filename1, "%s/Photcodes.dat", input1);
     36  if (!LoadPhotcodes (filename1, NULL, FALSE)) {
     37    fprintf (stderr, "error loading photcode table %s\n", filename1);
     38    exit (1);
     39  }
     40  input1Photcodes = GetPhotcodeTable();
     41
     42  // Read the input2 photcodes
     43  SetPhotcodeTable(NULL);
     44  sprintf (filename, "%s/Photcodes.dat", input2);
    2645  if (!LoadPhotcodes (filename, NULL, FALSE)) {
    2746    fprintf (stderr, "error loading photcode table %s\n", filename);
    2847    exit (1);
    2948  }
    30   // save the photcodes in the output catdir
     49  input2Photcodes = GetPhotcodeTable();
     50
     51  if (ALTERNATE_PHOTCODE_FILE) {
     52    SetPhotcodeTable(NULL);
     53    // First attempt to read file as fits table ...
     54    if (!LoadPhotcodes(ALTERNATE_PHOTCODE_FILE, NULL, FALSE)) {
     55      // and if that fails try it as a text file
     56      if (!LoadPhotcodesText(ALTERNATE_PHOTCODE_FILE)) {
     57        fprintf (stderr, "failed to read photcode file %s\n", ALTERNATE_PHOTCODE_FILE);
     58        exit (1);
     59      }
     60    }
     61    outputPhotcodes = GetPhotcodeTable();
     62
     63    secfiltMap1 = GetSecFiltMap(outputPhotcodes, input1Photcodes);
     64    if (!secfiltMap1) {
     65      fprintf (stderr, "failed to map %s secfilt photcodes to alternate photcodes table %s\n", input1,
     66        ALTERNATE_PHOTCODE_FILE);
     67      exit (1);
     68    }
     69
     70    codesFilename = ALTERNATE_PHOTCODE_FILE;
     71
     72  } else {
     73    // Use the first input file's photcodes
     74    outputPhotcodes = input1Photcodes;
     75    codesFilename = filename1;
     76    // secfitMap1 can remain NULL since input1 defines the set of codes
     77  }
     78
     79  secfiltMap2 = GetSecFiltMap(outputPhotcodes, input2Photcodes);
     80  if (!secfiltMap2) {
     81    fprintf (stderr, "falied to map %s secfilt photcodes to output photcodes table %s\n", input2,
     82      codesFilename);
     83    exit (1);
     84  }
     85
     86  // save the output photcodes in the output catdir
     87  SetPhotcodeTable(outputPhotcodes);
    3188  sprintf (filename, "%s/Photcodes.dat", output);
    3289  if (!check_file_access (filename, TRUE, TRUE, VERBOSE)) {
    … …  
    64121    if (VERBOSE) fprintf (stderr, "output: %s\n", outsky[0].regions[i].name);
    65122
     123    if (outputPhotcodes) {
     124        SetPhotcodeTable(outputPhotcodes);
     125    }
    66126    // load / create output catalog
    67127    LoadCatalog (&outcatalog, &outsky[0].regions[i], outsky[0].filename[i], "w");
    68128
     129    if (input1Photcodes) {
     130        SetPhotcodeTable(input1Photcodes);
     131    }
    69132    // combine only tables at equal or larger depth
    70133     
    … …  
    84147      }
    85148      dvo_update_image_IDs (&IDmap1, &incatalog);
    86       merge_catalogs_new (&outsky[0].regions[i], &outcatalog, &incatalog);
     149      merge_catalogs_new (&outsky[0].regions[i], &outcatalog, &incatalog, secfiltMap1);
    87150      dvo_catalog_unlock (&incatalog);
    88151      dvo_catalog_free (&incatalog);
    89152    }
    90153    SkyListFree (inlist);
     154
     155    if (input2Photcodes) {
     156      SetPhotcodeTable(input2Photcodes);
     157    }
    91158
    92159    // load in all of the tables from input2 for this region
    … …  
    95162      if (VERBOSE) fprintf (stderr, "input 2: %s\n", inlist[0].regions[j][0].name);
    96163
    97       // load input catalog (1)
     164      // load input catalog (2)
    98165      LoadCatalog (&incatalog, inlist[0].regions[j], inlist[0].filename[j], "r");
    99166
    … …  
    105172      }
    106173      dvo_update_image_IDs (&IDmap2, &incatalog);
    107       merge_catalogs_old (&outsky[0].regions[i], &outcatalog, &incatalog, RADIUS);
     174      merge_catalogs_old (&outsky[0].regions[i], &outcatalog, &incatalog, RADIUS, secfiltMap2);
    108175      dvo_catalog_unlock (&incatalog);
    109176      dvo_catalog_free (&incatalog);

  • trunk/Ohana/src/dvomerge/src/dvomergeUpdate.c

    r28354 r28855  
    1010  char filename[256], *input, *output;
    1111  IDmapType IDmap;
     12  PhotCodeData *inputPhotcodes;
     13  PhotCodeData *outputPhotcodes;
     14  int *secfiltMap = NULL;
    1215
    1316  if (strcasecmp (argv[2], "into")) dvomerge_usage();
    … …  
    1619  output = argv[3];
    1720
     21  if (ALTERNATE_PHOTCODE_FILE) {
     22    fprintf (stderr, "cannot specify photcodes when merging into an existing catdir\n");
     23    exit (1);
     24  }
     25
    1826  // since we are merging the input db into the output db, the output defines the photcode
    19   // table & db layout but, this requires the output to exist.  if it does not, instead use the
    20   // input.
    2127  sprintf (filename, "%s/Photcodes.dat", output);
    2228  if (!LoadPhotcodes (filename, NULL, FALSE)) {
    23     sprintf (filename, "%s/Photcodes.dat", input);
    24     if (!LoadPhotcodes (filename, NULL, FALSE)) {
    25       fprintf (stderr, "error loading photcode table: tried %s/Photcodes.dat and %s/Photcodes.dat\n", output, input);
    26       exit (1);
    27     }
    28     if (!check_dir_access (output, VERBOSE)) {
    29       fprintf (stderr, "error creating output database directory %s\n", output);
    30       exit (1);
    31     }   
    32     sprintf (filename, "%s/Photcodes.dat", output);
    33     if (!SavePhotcodesFITS (filename)) {
    34       fprintf (stderr, "error saving photcode table in %s/Photcodes.dat\n", output);
    35       exit (1);
    36     }   
     29    fprintf (stderr, "error loading photcode table: %s/Photcodes.dat\n", output);
     30    exit (1);
     31  }
     32  outputPhotcodes = GetPhotcodeTable();
     33
     34  SetPhotcodeTable(NULL);
     35  sprintf (filename, "%s/Photcodes.dat", input);
     36  if (!LoadPhotcodes (filename, NULL, FALSE)) {
     37    fprintf (stderr, "error reading input database directory %s\n", input);
     38    exit (1);
     39  }     
     40  inputPhotcodes = GetPhotcodeTable();
     41  secfiltMap = GetSecFiltMap(outputPhotcodes, inputPhotcodes);
     42  if (!secfiltMap) {
     43    fprintf (stderr, "failed to map input secfilt photcodes to output photcodes table\n");
     44    exit (1);
    3745  }
    3846
    … …  
    6371    if (VERBOSE) fprintf (stderr, "output: %s\n", outsky[0].regions[i].name);
    6472
     73    SetPhotcodeTable(outputPhotcodes);
    6574    // load / create output catalog (if catalog does not exist, it will be created)
    6675    LoadCatalog (&outcatalog, &outsky[0].regions[i], outsky[0].filename[i], "w");
    … …  
    7786
    7887      // load input catalog
     88      SetPhotcodeTable(inputPhotcodes);
    7989      LoadCatalog (&incatalog, inlist[0].regions[j], inlist[0].filename[j], "r");
    8090
    … …  
    8696      }
    8797      dvo_update_image_IDs (&IDmap, &incatalog);
    88       merge_catalogs_old (&outsky[0].regions[i], &outcatalog, &incatalog, RADIUS);
     98      merge_catalogs_old (&outsky[0].regions[i], &outcatalog, &incatalog, RADIUS, secfiltMap);
    8999      dvo_catalog_unlock (&incatalog);
    90100      dvo_catalog_free (&incatalog);

  • trunk/Ohana/src/dvomerge/src/merge_catalogs_new.c

    r28241 r28855  
    77// merge the input data into the output catalog
    88// input entries always define new objects
    9 int merge_catalogs_new (SkyRegion *region, Catalog *output, Catalog *input) {
     9
     10#define notyet 1
     11
     12int merge_catalogs_new (SkyRegion *region, Catalog *output, Catalog *input, int *secfiltMap) {
    1013 
    1114  off_t i, j, offset;
    12   off_t NAVERAGE, NMEASURE, Naverage, Nmeasure, Nsecfilt, Nm;
     15  off_t NAVERAGE, NMEASURE, Naverage, Nmeasure, NsecfiltIn, NsecfiltOut, Nm;
    1316
    1417  Naverage = output[0].Naverage;
    1518  Nmeasure = output[0].Nmeasure;
    16   Nsecfilt = output[0].Nsecfilt;
    17   assert (output[0].Nsecfilt == input[0].Nsecfilt);
     19
     20  NsecfiltOut = output[0].Nsecfilt;
     21  NsecfiltIn  = input[0].Nsecfilt;
    1822
    1923  NAVERAGE = Naverage + 100;
    … …  
    2125
    2226  REALLOCATE (output[0].average, Average, NAVERAGE);
    23   REALLOCATE (output[0].secfilt, SecFilt, NAVERAGE*Nsecfilt);
     27  REALLOCATE (output[0].secfilt, SecFilt, NAVERAGE*NsecfiltOut);
    2428  REALLOCATE (output[0].measure, Measure, NMEASURE);
    2529
    … …  
    3135      output[0].average[Naverage].measureOffset = Nmeasure;
    3236
    33       // XXX input[0].Nsecfilt == output[0].Nsecfilt ???
    34       for (j = 0; j < Nsecfilt; j++) {
    35           output[0].secfilt[Nsecfilt*Naverage+j] = input[0].secfilt[Nsecfilt*i+j];
     37      if (secfiltMap) {
     38        // the output secfilt photcodes may be in a different place in the output table
     39        assert (output[0].Nsecfilt >= input[0].Nsecfilt);
     40        for (j = 0; j < NsecfiltIn; j++) {
     41          output[0].secfilt[NsecfiltOut*Naverage + secfiltMap[j]] = input[0].secfilt[NsecfiltIn*i+j];
     42        }
     43      } else {
     44        assert (output[0].Nsecfilt == input[0].Nsecfilt);
     45        for (j = 0; j < NsecfiltIn; j++) {
     46            output[0].secfilt[NsecfiltOut*Naverage+j] = input[0].secfilt[NsecfiltIn*i+j];
     47        }
    3648      }
    3749
    … …  
    5567          NAVERAGE += 50;
    5668          REALLOCATE (output[0].average, Average, NAVERAGE);
    57           REALLOCATE (output[0].secfilt, SecFilt, NAVERAGE*Nsecfilt);
     69          REALLOCATE (output[0].secfilt, SecFilt, NAVERAGE*NsecfiltOut);
    5870      }
    5971  }
    6072  REALLOCATE (output[0].average, Average, MAX (Naverage, 1));
    6173  REALLOCATE (output[0].measure, Measure, MAX (Nmeasure, 1));
    62   REALLOCATE (output[0].secfilt, SecFilt, Nsecfilt*MAX (Naverage, 1));
     74  REALLOCATE (output[0].secfilt, SecFilt, NsecfiltOut*MAX (Naverage, 1));
    6375  output[0].Naverage = Naverage;
    6476  output[0].Nmeasure = Nmeasure;
    65   output[0].Nsecf_mem = Naverage * Nsecfilt;
     77  output[0].Nsecf_mem = Naverage * NsecfiltOut;
    6678
    6779  // If we are using dvomergeCreate to split an existing catalog, then the max objID in

  • trunk/Ohana/src/dvomerge/src/merge_catalogs_old.c

    r28354 r28855  
    1616// merge the input data into the output catalog
    1717// input entries are matched by position to the new objects
    18 int merge_catalogs_old (SkyRegion *region, Catalog *output, Catalog *input, double RADIUS) {
     18int merge_catalogs_old (SkyRegion *region, Catalog *output, Catalog *input, double RADIUS, int *secfiltMap) {
    1919
    2020  off_t i, j, k, Nin, offset, J, Jmin, status, Nstars;
    … …  
    2424  off_t *N1, *N2,  *next_meas;
    2525  off_t Nave, NAVE, Nmeas, NMEAS, Nmatch;
    26   int Nsecfilt;
     26  int NsecfiltIn;
     27  int NsecfiltOut;
    2728  unsigned int objID, catID;
    2829  Coords tcoords;
    … …  
    3132  // gettimeofday (&start, (void *) NULL);
    3233
    33   Nsecfilt = output[0].Nsecfilt;
     34  NsecfiltOut = output[0].Nsecfilt;
     35  NsecfiltIn  = input[0].Nsecfilt;
     36  assert (secfiltMap || (NsecfiltOut == NsecfiltIn));
    3437
    3538  /** allocate local arrays (stars) **/
    … …  
    228231    // update the average properties to reflect the incoming entries:
    229232    // if the original value is NAN but the input value is not, accept the input:
    230     for (k = 0; k < Nsecfilt; k++) {
    231       if ( isfinite(output[0].secfilt[n*Nsecfilt+k].M)) continue;
    232       if (!isfinite( input[0].secfilt[N*Nsecfilt+k].M)) continue;
    233       output[0].secfilt[n*Nsecfilt+k] = input[0].secfilt[N*Nsecfilt+k];
     233    if (!secfiltMap) {
     234      // secfilt tables of input and output are the same
     235      for (k = 0; k < NsecfiltIn; k++) {
     236        if ( isfinite(output[0].secfilt[n*NsecfiltIn+k].M)) continue;
     237        if (!isfinite( input[0].secfilt[N*NsecfiltIn+k].M)) continue;
     238        output[0].secfilt[n*NsecfiltIn+k] = input[0].secfilt[N*NsecfiltIn+k];
     239      }
     240    } else {
     241      for (k = 0; k < NsecfiltIn; k++) {
     242        // index for this entry in output's secfilt list
     243        int outputIndex = n * NsecfiltOut + secfiltMap[k];
     244       
     245        if ( isfinite(output[0].secfilt[outputIndex].M)) continue;
     246        if (!isfinite( input[0].secfilt[N*NsecfiltIn+k].M)) continue;
     247        output[0].secfilt[outputIndex] = input[0].secfilt[N*NsecfiltIn+k];
     248      }
    234249    }
    235250
    … …  
    255270      NAVE = Nave + 1000;
    256271      REALLOCATE (output[0].average, Average, NAVE);
    257       REALLOCATE (output[0].secfilt, SecFilt, NAVE*output[0].Nsecfilt);
     272      REALLOCATE (output[0].secfilt, SecFilt, NAVE*NsecfiltOut);
    258273    }
    259274
    … …  
    302317    objID ++;
    303318
    304     for (j = 0; j < Nsecfilt; j++) {
    305       if (isfinite(input[0].secfilt[N*Nsecfilt+j].M)) {
    306         output[0].secfilt[Nave*Nsecfilt+j].M     = input[0].secfilt[N*Nsecfilt+j].M;
    307         output[0].secfilt[Nave*Nsecfilt+j].dM    = input[0].secfilt[N*Nsecfilt+j].dM;
    308         output[0].secfilt[Nave*Nsecfilt+j].Xm    = input[0].secfilt[N*Nsecfilt+j].Xm;
    309         output[0].secfilt[Nave*Nsecfilt+j].M_20  = input[0].secfilt[N*Nsecfilt+j].M_20;
    310         output[0].secfilt[Nave*Nsecfilt+j].M_80  = input[0].secfilt[N*Nsecfilt+j].M_80;
    311         output[0].secfilt[Nave*Nsecfilt+j].Ncode = input[0].secfilt[N*Nsecfilt+j].Ncode;
    312         output[0].secfilt[Nave*Nsecfilt+j].Nused = input[0].secfilt[N*Nsecfilt+j].Nused;
     319    for (j = 0; j < NsecfiltIn; j++) {
     320      int outputIndex;
     321      if (secfiltMap) {
     322        outputIndex = (Nave * NsecfiltOut) + secfiltMap[j];
    313323      } else {
    314         output[0].secfilt[Nave*Nsecfilt+j].M     = NAN;
    315         output[0].secfilt[Nave*Nsecfilt+j].dM    = NAN;
    316         output[0].secfilt[Nave*Nsecfilt+j].Xm    = NAN_S_SHORT;
    317         output[0].secfilt[Nave*Nsecfilt+j].M_20  = NAN_S_SHORT;
    318         output[0].secfilt[Nave*Nsecfilt+j].M_80  = NAN_S_SHORT;
    319         output[0].secfilt[Nave*Nsecfilt+j].Ncode = 0;
    320         output[0].secfilt[Nave*Nsecfilt+j].Nused = 0;
     324        outputIndex = (Nave * NsecfiltIn) + j;
     325      }
     326   
     327      if (isfinite(input[0].secfilt[N*NsecfiltIn+j].M)) {
     328        output[0].secfilt[outputIndex].M     = input[0].secfilt[N*NsecfiltIn+j].M;
     329        output[0].secfilt[outputIndex].dM    = input[0].secfilt[N*NsecfiltIn+j].dM;
     330        output[0].secfilt[outputIndex].Xm    = input[0].secfilt[N*NsecfiltIn+j].Xm;
     331        output[0].secfilt[outputIndex].M_20  = input[0].secfilt[N*NsecfiltIn+j].M_20;
     332        output[0].secfilt[outputIndex].M_80  = input[0].secfilt[N*NsecfiltIn+j].M_80;
     333        output[0].secfilt[outputIndex].Ncode = input[0].secfilt[N*NsecfiltIn+j].Ncode;
     334        output[0].secfilt[outputIndex].Nused = input[0].secfilt[N*NsecfiltIn+j].Nused;
     335      } else {
     336        output[0].secfilt[outputIndex].M     = NAN;
     337        output[0].secfilt[outputIndex].dM    = NAN;
     338        output[0].secfilt[outputIndex].Xm    = NAN_S_SHORT;
     339        output[0].secfilt[outputIndex].M_20  = NAN_S_SHORT;
     340        output[0].secfilt[outputIndex].M_80  = NAN_S_SHORT;
     341        output[0].secfilt[outputIndex].Ncode = 0;
     342        output[0].secfilt[outputIndex].Nused = 0;
    321343      }
    322344    }
    … …  
    376398  output[0].Naverage = Nave;
    377399  output[0].Nmeasure = Nmeas;
    378   output[0].Nsecf_mem = Nave*Nsecfilt;
     400  output[0].Nsecf_mem = Nave*NsecfiltOut;
    379401  if (VERBOSE) fprintf (stderr, "Nstars, Nave, Nmeas: "OFF_T_FMT" "OFF_T_FMT" "OFF_T_FMT", ("OFF_T_FMT" matches)\n",  Nstars,  Nave,  Nmeas,  Nmatch);
    380402

  • trunk/Ohana/src/libdvo/include/dvo.h

    r27817 r28855  
    318318
    319319PhotCodeData *GetPhotcodeTable (void);
     320void SetPhotcodeTable (PhotCodeData *);
     321
     322int *GetSecFiltMap(PhotCodeData *ouput, PhotCodeData* input);
    320323
    321324int LoadPhotcodes (char *catdir_file, char *master_file, int readwrite);

  • trunk/Ohana/src/libdvo/src/dvo_photcode_ops.c

    r24974 r28855  
    2121PhotCodeData *GetPhotcodeTable () {
    2222
    23   /* allocate space to photcode table, free existing data */
    2423  if (photcodes == NULL) {
     24    /* allocate space to photcode table */
    2525    ALLOCATE (photcodes, PhotCodeData, 1);
    2626    photcodes[0].code = NULL;
    2727  }
    2828  return photcodes;
     29}
     30
     31// set the photcode table. This is used to switch between serveral previously
     32// allocated tables
     33void SetPhotcodeTable (PhotCodeData *new) {
     34  photcodes = new; 
    2935}
    3036
    … …  
    531537}
    532538
     539// Create a map between the secfilt values from one photcode table to another
     540int *GetSecFiltMap(PhotCodeData *output, PhotCodeData *input) {
     541  int i;
     542  int NsecfiltIn  = input[0].Nsecfilt;
     543  int *map;
     544  ALLOCATE(map, int, NsecfiltIn);
     545
     546  // loop over entries
     547  for (i = 0; i < NsecfiltIn; i++) {
     548    int code  = input[0].codeNsec[i];
     549    int entry = output[0].hashcode[code];
     550    if (entry == -1) {
     551      // entry is missing
     552      free(map);
     553      return(FALSE);
     554    }
     555    map[i] = entry;
     556  }
     557
     558  return map;
     559}
     560
    533561/* photcode table should have the following format:
    534562
Note: See TracChangeset for help on using the changeset viewer.