IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 28674


Ignore:
Timestamp:
Jul 15, 2010, 6:27:55 AM (16 years ago)
Author:
eugene
Message:

better checks for valid WCS; fix order in dvomerge loop; update test for dvomerge

Location:
branches/eam_branches/ipp-20100621/Ohana/src
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ipp-20100621/Ohana/src/addstar/test/dvomerge.dvo

    r28368 r28674  
    9595    subset T = t if (id == imageID[$i])
    9696    set dT = T - time[$i]
    97     vstat -q dT
     97    vstat dT
    9898    tapOK {abs($MEAN)  < 0.00001} "time for measure ID $i (MEAN)"
    9999    tapOK {abs($SIGMA) < 0.00001} "time for measure ID $i (SIGMA)"

  • branches/eam_branches/ipp-20100621/Ohana/src/dvomerge/src/dvomergeUpdate.c

    r28354 r28674  
    66  off_t i, j, Ns, Ne;
    77  SkyTable *outsky, *insky;
    8   SkyList *inlist;
     8  SkyList *outlist;
    99  Catalog incatalog, outcatalog;
    1010  char filename[256], *input, *output;
    … …  
    5858  depth = insky[0].regions[Ns].depth;
    5959 
    60   // loop over the populatable output regions
     60  // loop over the populated input regions
    6161  for (i = 0; i < outsky[0].Nregions; i++) {
    62     if (!outsky[0].regions[i].table) continue;
    63     if (VERBOSE) fprintf (stderr, "output: %s\n", outsky[0].regions[i].name);
     62    if (!insky[0].regions[i].table) continue;
     63    if (VERBOSE) fprintf (stderr, "input: %s\n", insky[0].regions[i].name);
    6464
    6565    // load / create output catalog (if catalog does not exist, it will be created)
    66     LoadCatalog (&outcatalog, &outsky[0].regions[i], outsky[0].filename[i], "w");
     66    LoadCatalog (&incatalog, &insky[0].regions[i], insky[0].filename[i], "r");
     67    // skip empty input catalogs
     68    if (!incatalog.Naves_disk) {
     69        dvo_catalog_unlock (&incatalog);
     70        dvo_catalog_free (&incatalog);
     71        continue;
     72    }
    6773
    6874    // combine only tables at equal or larger depth
    … …  
    7278    // compare to a slightly reduced footprint
    7379    float dPos = 2.0/3600.0;
    74     inlist = SkyListByBounds (insky, depth, outsky[0].regions[i].Rmin + dPos, outsky[0].regions[i].Rmax - dPos, outsky[0].regions[i].Dmin + dPos, outsky[0].regions[i].Dmax - dPos);
    75     for (j = 0; j < inlist[0].Nregions; j++) {
    76       if (VERBOSE) fprintf (stderr, "input : %s\n", inlist[0].regions[j][0].name);
     80    outlist = SkyListByBounds (outsky, depth, insky[0].regions[i].Rmin + dPos, insky[0].regions[i].Rmax - dPos, insky[0].regions[i].Dmin + dPos, insky[0].regions[i].Dmax - dPos);
     81    for (j = 0; j < outlist[0].Nregions; j++) {
     82      if (VERBOSE) fprintf (stderr, "output : %s\n", outlist[0].regions[j][0].name);
    7783
    7884      // load input catalog
    79       LoadCatalog (&incatalog, inlist[0].regions[j], inlist[0].filename[j], "r");
     85      LoadCatalog (&outcatalog, outlist[0].regions[j], outlist[0].filename[j], "w");
    8086
    81       // skip empty input catalogs
    82       if (!incatalog.Naves_disk) {
    83         dvo_catalog_unlock (&incatalog);
    84         dvo_catalog_free (&incatalog);
    85         continue;
    86       }
    8787      dvo_update_image_IDs (&IDmap, &incatalog);
    88       merge_catalogs_old (&outsky[0].regions[i], &outcatalog, &incatalog, RADIUS);
    89       dvo_catalog_unlock (&incatalog);
    90       dvo_catalog_free (&incatalog);
     88      merge_catalogs_old (&outsky[0].regions[j], &outcatalog, &incatalog, RADIUS);
    9189
    92       fprintf (stderr, "merged %s into %s\n", outsky[0].regions[i].name, inlist[0].regions[j][0].name);
     90      outcatalog.catflags = LOAD_AVES | LOAD_MEAS | LOAD_MISS | LOAD_SECF;
     91      dvo_catalog_save (&outcatalog, VERBOSE);
     92      dvo_catalog_unlock (&outcatalog);
     93      dvo_catalog_free (&outcatalog);
     94
     95      fprintf (stderr, "merged %s into %s\n", insky[0].regions[i].name, outlist[0].regions[j][0].name);
    9396    }
    94     SkyListFree (inlist);
     97    SkyListFree (outlist);
    9598
    96     outcatalog.catflags = LOAD_AVES | LOAD_MEAS | LOAD_MISS | LOAD_SECF;
    97     dvo_catalog_save (&outcatalog, VERBOSE);
    98     dvo_catalog_unlock (&outcatalog);
    99     dvo_catalog_free (&outcatalog);
     99    dvo_catalog_unlock (&incatalog);
     100    dvo_catalog_free (&incatalog);
    100101  }
    101102

  • branches/eam_branches/ipp-20100621/Ohana/src/libdvo/src/coordops.c

    r27941 r28674  
    421421}
    422422
     423enum {COORD_TYPE_NONE, COORD_TYPE_PC, COORD_TYPE_ROT, COORD_TYPE_CD, COORD_TYPE_LIN};
     424
    423425int GetCoords (Coords *coords, Header *header) {
    424426 
    … …  
    427429  double equinox;
    428430  char *ctype;
     431  int mode;
    429432 
    430433  rotate = 0.0;
    … …  
    436439  strcpy (coords[0].ctype, "NONE");
    437440 
    438   status = FALSE;
    439   if (gfits_scan (header, "CTYPE2", "%s", 1, coords[0].ctype)) {
    440     status  = gfits_scan (header, "CRVAL1", "%lf", 1, &coords[0].crval1);
    441     status &= gfits_scan (header, "CRPIX1", "%f", 1, &coords[0].crpix1);
    442     status &= gfits_scan (header, "CRVAL2", "%lf", 1, &coords[0].crval2); 
    443     status &= gfits_scan (header, "CRPIX2", "%f", 1, &coords[0].crpix2);
    444 
    445     if (gfits_scan (header, "CDELT1", "%f", 1, &coords[0].cdelt1)) {
    446       status &= gfits_scan (header, "CDELT2", "%f", 1, &coords[0].cdelt2);
    447       if (gfits_scan (header, "CROTA2", "%lf", 1, &rotate)) {
    448         Lambda = coords[0].cdelt2 / coords[0].cdelt1;
    449         coords[0].pc1_1 =  cos(rotate*RAD_DEG);
    450         coords[0].pc1_2 = -sin(rotate*RAD_DEG) * Lambda;
    451         coords[0].pc2_1 =  sin(rotate*RAD_DEG) / Lambda;
    452         coords[0].pc2_2 =  cos(rotate*RAD_DEG);
    453       }
    454       if (gfits_scan (header, "PC001001", "%f", 1, &coords[0].pc1_1)) {
    455         status &= gfits_scan (header, "PC001002", "%f", 1, &coords[0].pc1_2);
    456         status &= gfits_scan (header, "PC002001", "%f", 1, &coords[0].pc2_1);
    457         status &= gfits_scan (header, "PC002002", "%f", 1, &coords[0].pc2_2);
    458       }
     441  mode = COORD_TYPE_NONE;
     442  {   
     443    int haveCTYPE, haveCDELT, haveCROTA, haveCDij, havePCij, haveRAo;
     444    float tmp;
     445    char stmp[80];
     446   
     447    // there are a few different representations for scale and rotation.  choose an appropriate
     448    // set: (CDELTi + CROTAi), (CDELTi + PCij), (CDij),
     449
     450    haveCTYPE = gfits_scan (header, "CTYPE2",   "%s", 1, stmp);
     451    haveCDELT = gfits_scan (header, "CDELT1",   "%f", 1, &tmp);
     452    haveCROTA = gfits_scan (header, "CROTA1",   "%f", 1, &tmp);
     453    haveCDij  = gfits_scan (header, "CD1_1",    "%f", 1, &tmp);
     454    havePCij  = gfits_scan (header, "PC001001", "%f", 1, &tmp);
     455    haveRAo   = gfits_scan (header, "RA_O",     "%f", 1, &tmp);
     456   
     457    if (haveCTYPE && havePCij  && haveCDELT) { mode = COORD_TYPE_PC;   goto gotit; }
     458    if (haveCTYPE && haveCROTA && haveCDELT) { mode = COORD_TYPE_ROT;  goto gotit; }
     459    if (haveCTYPE && haveCDij)               { mode = COORD_TYPE_CD;   goto gotit; }
     460    if (haveRAo)                             { mode = COORD_TYPE_LIN;  goto gotit; }
     461    // fprintf (stderr, "no valid WCS keywords\n");
     462    return (FALSE);
     463  }
     464 
     465gotit:
     466
     467  status = TRUE;
     468  switch (mode) {
     469    case COORD_TYPE_PC:
     470      status &= gfits_scan (header, "CTYPE2",   "%s",  1, coords[0].ctype);
     471      status &= gfits_scan (header, "CRVAL1",   "%lf", 1, &coords[0].crval1);
     472      status &= gfits_scan (header, "CRPIX1",   "%f",  1, &coords[0].crpix1);
     473      status &= gfits_scan (header, "CRVAL2",   "%lf", 1, &coords[0].crval2); 
     474      status &= gfits_scan (header, "CRPIX2",   "%f",  1, &coords[0].crpix2);
     475
     476      status &= gfits_scan (header, "CDELT1",   "%f",  1, &coords[0].cdelt1);
     477      status &= gfits_scan (header, "CDELT2",   "%f",  1, &coords[0].cdelt2);
     478      status &= gfits_scan (header, "PC001001", "%f",  1, &coords[0].pc1_1);
     479      status &= gfits_scan (header, "PC001002", "%f",  1, &coords[0].pc1_2);
     480      status &= gfits_scan (header, "PC002001", "%f",  1, &coords[0].pc2_1);
     481      status &= gfits_scan (header, "PC002002", "%f",  1, &coords[0].pc2_2);
    459482
    460483      /* set NPLYTERM based on header.  if NPLYTERM is missing, it should have a
    … …  
    490513          break;
    491514      }
    492     } else {
    493       if (gfits_scan (header, "CD1_1", "%f", 1, &coords[0].pc1_1)) {
    494         status &= gfits_scan (header, "CD1_2", "%f", 1, &coords[0].pc1_2);
    495         status &= gfits_scan (header, "CD2_1", "%f", 1, &coords[0].pc2_1);
    496         status &= gfits_scan (header, "CD2_2", "%f", 1, &coords[0].pc2_2);
    497         /* renormalize */
    498         scale = hypot (coords[0].pc1_1, coords[0].pc1_2);
    499         coords[0].cdelt1 = coords[0].cdelt2 = scale;
    500         coords[0].pc1_1 /= scale;
    501         coords[0].pc1_2 /= scale;
    502         coords[0].pc2_1 /= scale;
    503         coords[0].pc2_2 /= scale;
    504       } else {
    505         status = FALSE;
    506       }
    507     }
    508   } else {
    509     /* some of my thesis data uses this simple linear model - convert on read? */
    510     if (gfits_scan (header, "RA_O", "%lf", 1, &coords[0].crval1)) {
    511       status  = gfits_scan (header, "RA_X", "%f", 1, &coords[0].pc1_1);
     515      break;
     516
     517    case COORD_TYPE_ROT:
     518      status &= gfits_scan (header, "CTYPE2", "%s",  1, coords[0].ctype);
     519      status &= gfits_scan (header, "CRVAL1", "%lf", 1, &coords[0].crval1);
     520      status &= gfits_scan (header, "CRPIX1", "%f",  1, &coords[0].crpix1);
     521      status &= gfits_scan (header, "CRVAL2", "%lf", 1, &coords[0].crval2); 
     522      status &= gfits_scan (header, "CRPIX2", "%f",  1, &coords[0].crpix2);
     523
     524      status &= gfits_scan (header, "CDELT1", "%f", 1, &coords[0].cdelt1);
     525      status &= gfits_scan (header, "CDELT2", "%f", 1, &coords[0].cdelt2);
     526
     527      status &= gfits_scan (header, "CROTA2", "%lf", 1, &rotate);
     528      Lambda = coords[0].cdelt2 / coords[0].cdelt1;
     529      coords[0].pc1_1 =  cos(rotate*RAD_DEG);
     530      coords[0].pc1_2 = -sin(rotate*RAD_DEG) * Lambda;
     531      coords[0].pc2_1 =  sin(rotate*RAD_DEG) / Lambda;
     532      coords[0].pc2_2 =  cos(rotate*RAD_DEG);
     533      break;
     534
     535    case COORD_TYPE_CD:
     536      status &= gfits_scan (header, "CTYPE2", "%s",  1, coords[0].ctype);
     537      status &= gfits_scan (header, "CRVAL1", "%lf", 1, &coords[0].crval1);
     538      status &= gfits_scan (header, "CRPIX1", "%f",  1, &coords[0].crpix1);
     539      status &= gfits_scan (header, "CRVAL2", "%lf", 1, &coords[0].crval2); 
     540      status &= gfits_scan (header, "CRPIX2", "%f",  1, &coords[0].crpix2);
     541
     542      status &= gfits_scan (header, "CD1_1", "%f", 1, &coords[0].pc1_1);
     543      status &= gfits_scan (header, "CD1_2", "%f", 1, &coords[0].pc1_2);
     544      status &= gfits_scan (header, "CD2_1", "%f", 1, &coords[0].pc2_1);
     545      status &= gfits_scan (header, "CD2_2", "%f", 1, &coords[0].pc2_2);
     546      /* renormalize */
     547      scale = hypot (coords[0].pc1_1, coords[0].pc1_2);
     548      coords[0].cdelt1 = coords[0].cdelt2 = scale;
     549      coords[0].pc1_1 /= scale;
     550      coords[0].pc1_2 /= scale;
     551      coords[0].pc2_1 /= scale;
     552      coords[0].pc2_2 /= scale;
     553      break;
     554
     555    case COORD_TYPE_LIN:
     556      /* some of my thesis data uses this simple linear model - convert on read? */
     557      status &= gfits_scan (header, "RA_O", "%lf", 1, &coords[0].crval1);
     558      status &= gfits_scan (header, "RA_X", "%f", 1, &coords[0].pc1_1);
    512559      status &= gfits_scan (header, "RA_Y", "%f", 1, &coords[0].pc1_2);
    513560      status &= gfits_scan (header, "DEC_O", "%lf", 1, &coords[0].crval2); 
    … …  
    517564      coords[0].cdelt1 = coords[0].cdelt2 = 1.0;
    518565      strcpy (coords[0].ctype, "GENE");
    519     }
    520   }
     566      break;
     567  }
     568
    521569  if (status) {
    522570    if (!gfits_scan (header, "EQUINOX", "%lf", 1, &equinox)) {
    … …  
    529577    }
    530578  }
     579
    531580  if (!status) {
    532     fprintf (stderr, "error getting all elements for coordinate mode %s\n", coords[0].ctype);
     581    // fprintf (stderr, "error getting all elements for coordinate mode %s\n", coords[0].ctype);
    533582    coords[0].crval1 = coords[0].crpix1 = coords[0].cdelt1 = 0.0;
    534583    coords[0].crval2 = coords[0].crpix2 = coords[0].cdelt2 = 0.0;
Note: See TracChangeset for help on using the changeset viewer.