IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 39471 for trunk


Ignore:
Timestamp:
Mar 16, 2016, 4:05:21 PM (10 years ago)
Author:
eugene
Message:

add code to fix the warp image IDs

Location:
trunk/Ohana/src/relastro
Files:
2 added
7 edited

Legend:

Unmodified
Added
Removed

  • trunk/Ohana/src/relastro/Makefile

    r39457 r39471  
    2525
    2626RELASTRO = \
     27$(SRC)/RepairWarps.$(ARCH).o \
     28$(SRC)/RepairWarpMeasures.$(ARCH).o \
     29$(SRC)/WarpImageMaps.$(ARCH).o \
     30$(SRC)/myIndex.$(ARCH).o \
    2731$(SRC)/ConfigInit.$(ARCH).o          \
    2832$(SRC)/FitChip.$(ARCH).o             \
    … …  
    105109
    106110RELASTRO_CLIENT = \
     111$(SRC)/RepairWarps.$(ARCH).o \
     112$(SRC)/RepairWarpMeasures.$(ARCH).o \
     113$(SRC)/WarpImageMaps.$(ARCH).o \
     114$(SRC)/myIndex.$(ARCH).o \
    107115$(SRC)/ConfigInit.$(ARCH).o          \
    108116$(SRC)/FitSimple.$(ARCH).o           \

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

    r39466 r39471  
    1010# define IDX_T int
    1111
     12# ifndef MAX_INT
     13# define MAX_INT 2147483647
     14# endif
     15
    1216# define LOGRTIME(MSG,...) {                            \
    1317    gettimeofday (&stopTimer, (void *) NULL);           \
    … …  
    2529typedef enum {FIT_NONE, FIT_AVERAGE, FIT_PM_ONLY, FIT_PAR_ONLY, FIT_PM_AND_PAR} FitMode;
    2630
    27 typedef enum {OP_NONE, OP_IMAGES, OP_HIGH_SPEED, OP_MERGE_SOURCE, OP_UPDATE_OBJECTS, OP_UPDATE_OFFSETS, OP_LOAD_OBJECTS, OP_HPM, OP_PARALLEL_REGIONS, OP_PARALLEL_IMAGES} RelastroOp;
     31typedef enum {OP_NONE, OP_IMAGES, OP_HIGH_SPEED, OP_MERGE_SOURCE, OP_UPDATE_OBJECTS, OP_UPDATE_OFFSETS, OP_LOAD_OBJECTS, OP_HPM, OP_PARALLEL_REGIONS, OP_PARALLEL_IMAGES, OP_REPAIR_WARPS} RelastroOp;
    2832
    2933typedef enum {TARGET_NONE, TARGET_SIMPLE, TARGET_CHIPS, TARGET_MOSAICS} FitTarget;
    … …  
    4549  MARK_BIG_OFFSET    = 0x0010,
    4650} MeasurementMask;
     51
     52typedef struct {
     53  int minID;
     54  int maxID;
     55  int *index;
     56  int Nindex;
     57  int NINDEX;
     58} myIndexType;
    4759
    4860typedef struct {
    … …  
    787799double dpsi_cauchy (double x);
    788800void my_memdump (char *message);
     801
     802int RepairWarps (SkyList *skylist, int hostID, char *hostpath);
     803int RepairWarpMeasures (Catalog *catalog);
     804
     805int FindWarpGroups (void);
     806void FreeWarpGroups (void);
     807int GetWarpID (Image *image, int obstime, unsigned short photcode, double Rave, double Dave, float X, float Y);
     808
     809myIndexType *myIndexInit ();
     810int myIndexFree (myIndexType *myIndex);
     811int myIndexUpdateLimits (myIndexType *myIndex, int value);
     812int myIndexSetRange (myIndexType *myIndex);
     813int myIndexSetEntry (myIndexType *myIndex, int value, int entry);
     814int myIndexGetEntry (myIndexType *myIndex, int value);

  • trunk/Ohana/src/relastro/src/RepairWarpMeasures.c

    r39470 r39471  
    22
    33// XXX need to load_images first and generate some lookup tables
    4 int RepairWarpMeasures (Catalog *catalog, int Ncatalog) {
     4int RepairWarpMeasures (Catalog *catalog) {
    55
    6   off_t i, j, Nimage;
    7   Image *image;
    8   double DPOS_MAX_ASEC;
     6  off_t Nimage;
     7  Image *image = getimages (&Nimage, NULL);
    98
    10   image = getimages (&Nimage, NULL);
     9  myAssert (!catalog->Nmeasure || catalog->measure, "programming error");
    1110
    12   for (i = 0; i < Ncatalog; i++) {
    13     myAssert (!catalog[i].Nmeasure || catalog[i].measureT, "programming error");
     11  for (off_t j = 0; j < catalog->Nmeasure; j++) {
     12    Measure *measure = &catalog->measure[j];
    1413
    15     for (j = 0; j < catalog[i].Nmeasure; j++) {
    16       Measure *measure = &catalog[i].measure[j];
    17 
    18       // we are only going to repair warp detections
    19       if (!isGPC1warp(measure->photcode))
    20         // XXX fix extID
    21         continue;
     14    // we are only going to repair warp detections
     15    if (!isGPC1warp(measure->photcode))  {
     16      // XXX fix extID
     17      continue;
    2218    }
    2319
    24       // get the associated average value:
    25       int averef = measure->averef;
    26       Average *average = &catalog[i].average[averef];
     20    // get the associated average value:
     21    int averef = measure->averef;
     22    Average *average = &catalog->average[averef];
    2723
    28       // double-check for consistency
    29       myAssert (average->objID == measure->objID, "mismatch");
     24    // double-check for consistency
     25    myAssert (average->objID == measure->objID, "mismatch");
    3026
    31       // check if this detection is far from it correct location
    32       double Rwrp = measure->R;
    33       double Dwrp = measure->D;
     27    // check if this detection is far from its correct location
     28    double Rwrp = measure->R;
     29    double Dwrp = measure->D;
    3430
    35       double Rave = average->R;
    36       double Dave = average->D;
     31    double Rave = average->R;
     32    double Dave = average->D;
    3733     
    38       float csdec = cos(Dave * RAD_DEG);
     34    float csdec = cos(Dave * RAD_DEG);
    3935
    40       // find the ra,dec displacement in arcsec:
    41       double dR = 3600.0*(Rwrp - Rave)*csdec;
    42       double dD = 3600.0*(Dwrp - Dave);
     36    // find the ra,dec displacement in arcsec:
     37    double dR = 3600.0*(Rwrp - Rave)*csdec;
     38    double dD = 3600.0*(Dwrp - Dave);
    4339
    44       // skip detections which are within a small distance of the expected location
    45       double dPos = hypot(dR,dD);
    46       if (dPos < 2.0) continue;
     40    // skip detections which are within a small distance of the expected location
     41    double dPos = hypot(dR,dD);
     42    if (dPos < 2.0) continue;
    4743
    48       // find the corrected warp ID:
    49       int warpID = getWarpGroup (measure->t, measure->photcode, Rstk, Dstk);
     44    // warp coordinates to confirm warp
     45    double X = measure->Xccd;
     46    double Y = measure->Yccd;
    5047
    51       // assert on coords being TAN?
    52       Coords *imcoords = &image[warpID].coords;
     48    // find the corrected warp ID:
     49    int warpID = GetWarpID (image, measure->t, measure->photcode, Rave, Dave, X, Y);
     50    myAssert (warpID >= 0, "error finding match");
    5351
    54       double X = measure->Xccd;
    55       double Y = measure->Yccd;
     52    // assert on coords being TAN?
     53    Coords *imcoords = &image[warpID].coords;
    5654
    57       double Rnew, Dnew;
    58       XY_to_RD (&Rnew, &Dnew, X, Y, imcoords);
     55    double Rnew, Dnew;
     56    XY_to_RD (&Rnew, &Dnew, X, Y, imcoords);
    5957
    60       Rnew = ohana_normalize_angle_to_midpoint(Rnew, 180.0);
    61       Dnew = ohana_normalize_angle_to_midpoint(Dnew,   0.0);
     58    Rnew = ohana_normalize_angle_to_midpoint(Rnew, 180.0);
     59    Dnew = ohana_normalize_angle_to_midpoint(Dnew,   0.0);
    6260
    63       // find the ra,dec displacement in arcsec:
    64       double dRnew = 3600.0*(Rnew - Rave)*csdec;
    65       double dDnew = 3600.0*(Dnew - Dave);
     61    // find the ra,dec displacement in arcsec:
     62    double dRnew = 3600.0*(Rnew - Rave)*csdec;
     63    double dDnew = 3600.0*(Dnew - Dave);
    6664
    67       // skip detections which are within a small distance of the expected location
    68       double dPos = hypot(dRnew,dDnew);
    69       if (dPos > 2.0) {
    70         // complain. abort?
    71         fprintf (stderr, "measurement still far from average location\n");
    72         abort();
    73       }
     65    // skip detections which are within a small distance of the expected location
     66    dPos = hypot(dRnew,dDnew);
     67    if (dPos > 2.0) {
     68      // complain. abort?
     69      fprintf (stderr, "measurement still far from average location\n");
     70      abort();
     71    }
    7472
    75       measure->R = Rnew;
    76       measure->D = Dnew;
    77     }
     73    measure->R = Rnew;
     74    measure->D = Dnew;
     75
     76    // XXX fix extID
    7877  }
    7978  return (TRUE);
    8079}
    81 
    82 

  • trunk/Ohana/src/relastro/src/WarpImageMaps.c

    r39470 r39471  
    1 
    2 int GetWarpID () {
    3 
     1# include "relastro.h"
     2
     3// we have a list of warp images, each with obstime & photcode.  I need to make a list of the obstime groups
     4
     5# define D_NOBSTIMES 10000
     6
     7typedef struct {
     8  int obstime;
     9  unsigned short photcode;
     10  int Nwarps;
     11  int NWARPS;
     12  off_t *warps;
     13  double *Rmin;
     14  double *Rmax;
     15  double *Dmin;
     16  double *Dmax;
     17} WarpGroup;
     18
     19WarpGroup *warpgroup = NULL;
     20int Nwarpgroup = 0;
     21int NWARPGROUP = 0;
     22
     23myIndexType *warpObstimeIndex = NULL;
     24
     25int FindWarpGroups (void) {
     26
     27  off_t Nimage;
     28  Image *image = getimages (&Nimage, NULL);
     29
     30  int *obstimes = NULL;
     31  int Nobstimes = 0;
     32  int NOBSTIMES = D_NOBSTIMES;
     33  ALLOCATE (obstimes, int, NOBSTIMES);
     34
     35  // first generate a list of obstime values
     36  for (off_t i = 0; i < Nimage; i++) {
     37    if (!isGPC1warp(image[i].photcode)) continue;
     38
     39    obstimes[Nobstimes] = image[i].tzero;
     40    Nobstimes ++;
     41
     42    CHECK_REALLOCATE (obstimes, int, NOBSTIMES, Nobstimes, D_NOBSTIMES);
     43  }
     44
     45  // sort the list:
     46  isort (obstimes, Nobstimes);
     47
     48  // find the unique obstimes
     49  int *uniqtimes = NULL;
     50  int *uniqcount = NULL;
     51  int Nuniqtimes = 0;
     52  ALLOCATE (uniqtimes, int, Nobstimes);
     53  ALLOCATE (uniqcount, int, Nobstimes);
     54
     55  // generate a uniq set of obstimes
     56  for (int i = 0; i < Nobstimes; Nuniqtimes ++) {
     57    uniqtimes[Nuniqtimes] = obstimes[i];
     58    int Ndup = 0;
     59    int lastValue = uniqtimes[Nuniqtimes];
     60    while ((i < Nobstimes) && (obstimes[i] == lastValue)) {
     61      i++;
     62      Ndup ++;
     63      uniqcount[Nuniqtimes] = Ndup;
     64    }
     65  }
     66  REALLOCATE (uniqtimes, int, Nuniqtimes);
     67  REALLOCATE (uniqcount, int, Nuniqtimes);
     68 
     69  // now I need to assign all warp images to a warpgroup
     70  // I have 2 options to go from warp obstime to warpgroup:
     71  // a) use a bisection lookup [may be slow]
     72  // b) use an index on obstime [require ~5*3.14*1e7*4 bytes ~ 600MB for index]
     73
     74  // generate an index on obstime
     75  warpObstimeIndex = myIndexInit ();
     76  warpObstimeIndex->minID = obstimes[0];
     77  warpObstimeIndex->maxID = obstimes[Nobstimes-1];
     78
     79  myIndexSetRange (warpObstimeIndex);
     80
     81  // generate the warp groups (create the index as we go)
     82  Nwarpgroup = Nuniqtimes;
     83  ALLOCATE (warpgroup, WarpGroup, Nwarpgroup);
     84 
     85  for (int i = 0; i < Nuniqtimes; i++) {
     86    warpgroup[i].obstime = uniqtimes[i];
     87    warpgroup[i].photcode = 0; // not yet set
     88    warpgroup[i].Nwarps = 0;
     89    warpgroup[i].NWARPS = 100;
     90    ALLOCATE (warpgroup[i].warps, off_t, warpgroup[i].NWARPS);
     91    myIndexSetEntry (warpObstimeIndex, uniqtimes[i], i);
     92  }
     93
     94  // assign all warps to one of the warpgroups
     95  for (off_t i = 0; i < Nimage; i++) {
     96    if (!isGPC1warp(image[i].photcode)) continue;
     97
     98    int seq = myIndexGetEntry (warpObstimeIndex, image[i].tzero);
     99    myAssert (warpgroup[seq].obstime == image[i].tzero, "oops");
     100
     101    if (!warpgroup[seq].photcode) {
     102      warpgroup[seq].photcode = image[i].photcode;
     103    } else {
     104      myAssert (warpgroup[seq].photcode == image[i].photcode, "oops");
     105    }
     106
     107    int N = warpgroup[seq].Nwarps;
     108    warpgroup[seq].warps[N] = i;
     109    warpgroup[seq].Nwarps ++;
     110   
     111    CHECK_REALLOCATE (warpgroup[seq].warps, off_t, warpgroup[seq].NWARPS, warpgroup[seq].Nwarps, 100);
     112  }
     113
     114  // we now have the warps assigned to groups.  now we need to generate the grid data to find the warp assignments
     115  for (int i = 0; i < Nwarpgroup; i++) {
     116
     117    myAssert (warpgroup[i].Nwarps == uniqcount[i], "failure");
     118    ALLOCATE (warpgroup[i].Rmin, double, warpgroup[i].Nwarps);
     119    ALLOCATE (warpgroup[i].Rmax, double, warpgroup[i].Nwarps);
     120    ALLOCATE (warpgroup[i].Dmin, double, warpgroup[i].Nwarps);
     121    ALLOCATE (warpgroup[i].Dmax, double, warpgroup[i].Nwarps);
     122
     123    for (int j = 0; j < warpgroup[i].Nwarps; j++) {
     124
     125      off_t N = warpgroup[i].warps[j];
     126     
     127      int Nx = image[N].NX;
     128      int Ny = image[N].NY;
     129
     130      double R, D;
     131      XY_to_RD (&R, &D, 0.0, 0.0, &image[N].coords);
     132      warpgroup[i].Rmin[j] = R;
     133      warpgroup[i].Rmax[j] = R;
     134      warpgroup[i].Dmin[j] = D;
     135      warpgroup[i].Dmax[j] = D;
     136
     137      // XXX need to worry about 0,360 boundary
     138      XY_to_RD (&R, &D, Nx, 0.0, &image[N].coords);
     139      warpgroup[i].Rmin[j] = MIN (warpgroup[i].Rmin[j], R);
     140      warpgroup[i].Rmax[j] = MAX (warpgroup[i].Rmax[j], R);
     141      warpgroup[i].Dmin[j] = MIN (warpgroup[i].Dmin[j], D);
     142      warpgroup[i].Dmax[j] = MAX (warpgroup[i].Dmax[j], D);
     143
     144      // XXX need to worry about 0,360 boundary
     145      XY_to_RD (&R, &D, 0.0, Ny, &image[N].coords);
     146      warpgroup[i].Rmin[j] = MIN (warpgroup[i].Rmin[j], R);
     147      warpgroup[i].Rmax[j] = MAX (warpgroup[i].Rmax[j], R);
     148      warpgroup[i].Dmin[j] = MIN (warpgroup[i].Dmin[j], D);
     149      warpgroup[i].Dmax[j] = MAX (warpgroup[i].Dmax[j], D);
     150
     151      // XXX need to worry about 0,360 boundary
     152      XY_to_RD (&R, &D, Nx, Ny, &image[N].coords);
     153      warpgroup[i].Rmin[j] = MIN (warpgroup[i].Rmin[j], R);
     154      warpgroup[i].Rmax[j] = MAX (warpgroup[i].Rmax[j], R);
     155      warpgroup[i].Dmin[j] = MIN (warpgroup[i].Dmin[j], D);
     156      warpgroup[i].Dmax[j] = MAX (warpgroup[i].Dmax[j], D);
     157      // bump Rmin,Rmax,Dmin,Dmax by some padding?
     158    }
     159  }
     160
     161  free (uniqtimes);
     162  free (uniqcount);
     163  free (obstimes);
     164  return TRUE;
    4165}
     166
     167int GetWarpID (Image *image, int obstime, unsigned short photcode, double Rave, double Dave, float X, float Y) {
     168
     169  int seq = myIndexGetEntry (warpObstimeIndex, obstime);
     170  myAssert (seq > -1, "ooops");
     171  myAssert (warpgroup[seq].photcode == photcode, "oops");
     172 
     173  // we now have the warp group, but which is the correct warp?
     174  for (int i = 0; i < warpgroup[seq].Nwarps; i++) {
     175    if (Rave < warpgroup[seq].Rmin[i]) continue;
     176    if (Rave > warpgroup[seq].Rmax[i]) continue;
     177    if (Dave < warpgroup[seq].Dmin[i]) continue;
     178    if (Dave > warpgroup[seq].Dmax[i]) continue;
     179
     180    // this is a possible image: check the coordinates:
     181    int N = warpgroup[seq].warps[i];
     182
     183    double Rwrp, Dwrp;
     184    XY_to_RD (&Rwrp, &Dwrp, X, Y, &image[N].coords);
     185
     186    float csdec = cos(Dave * RAD_DEG);
     187   
     188    // find the ra,dec displacement in arcsec:
     189    double dR = 3600.0*(Rwrp - Rave)*csdec;
     190    double dD = 3600.0*(Dwrp - Dave);
     191   
     192    // skip detections which are within a small distance of the expected location
     193    double dPos = hypot(dR,dD);
     194    if (dPos < 2.0) return N;
     195
     196    // check for multiple possible matches??
     197  }
     198  return -1;
     199}
     200
     201void FreeWarpGroups (void) {
     202
     203  for (int i = 0; i < Nwarpgroup; i++) {
     204    free (warpgroup[i].warps);
     205    free (warpgroup[i].Rmin);
     206    free (warpgroup[i].Rmax);
     207    free (warpgroup[i].Dmin);
     208    free (warpgroup[i].Dmax);
     209  }   
     210
     211  free (warpgroup);
     212  myIndexFree (warpObstimeIndex);
     213}

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

    r39466 r39471  
    5151    remove_argument (N, &argc, argv);
    5252    RELASTRO_OP = OP_UPDATE_OFFSETS;
     53  }
     54
     55  if ((N = get_argument (argc, argv, "-repair-warps"))) {
     56    remove_argument (N, &argc, argv);
     57    RELASTRO_OP = OP_REPAIR_WARPS;
    5358  }
    5459
    … …  
    666671    remove_argument (N, &argc, argv);
    667672    RELASTRO_OP = OP_UPDATE_OFFSETS;
     673  }
     674
     675  if ((N = get_argument (argc, argv, "-repair-warps"))) {
     676    remove_argument (N, &argc, argv);
     677    RELASTRO_OP = OP_REPAIR_WARPS;
    668678  }
    669679

  • trunk/Ohana/src/relastro/src/relastro.c

    r39457 r39471  
    6969      exit (0);
    7070
     71    case OP_REPAIR_WARPS: {
     72      FITS_DB db;
     73      if (!PARALLEL) {
     74     
     75        set_db (&db);
     76        gfits_db_init (&db);
     77
     78        /* lock and load the image db table */
     79        int status = dvo_image_lock (&db, ImageCat, 60.0, LCK_SOFT);
     80        if (!status) Shutdown ("ERROR: failure to lock image catalog %s", db.filename);
     81        if (db.dbstate == LCK_EMPTY) Shutdown ("ERROR: No images in catalog %s (1)", db.filename);
     82        if (!dvo_image_load (&db, VERBOSE, FALSE)) Shutdown ("can't read image catalog %s", db.filename);
     83        // the raw FITS data is freed by dvo_image_load, leaving only the Image structure data
     84
     85        /* load regions and images based on specified sky patch (default depth) */
     86        load_images (&db, skylist, FALSE, USE_ALL_IMAGES);
     87      }
     88
     89      // iterate over catalogs to make detection coordinates consistant
     90      RepairWarps (skylist, 0, NULL);
     91
     92      if (!PARALLEL) {
     93        gfits_db_free (&db);
     94        freeImages();
     95        freeMosaics ();
     96      }
     97
     98      relastro_free (sky, skylist);
     99      exit (0);
     100    }
     101
    71102    case OP_HIGH_SPEED:
    72103      /* high-speed is a 2pt cross-correlation process for linking moving objects (high PM) */

  • trunk/Ohana/src/relastro/src/relastro_client.c

    r39457 r39471  
    100100    }
    101101
     102      // XXX loading the images is fairly costly -- see if we can do an image subset?
     103    case OP_REPAIR_WARPS: {
     104      FITS_DB db;
     105     
     106      set_db (&db);
     107      gfits_db_init (&db);
     108
     109      /* lock and load the image db table */
     110      int status = dvo_image_lock (&db, ImageCat, 60.0, LCK_SOFT);
     111      if (!status) Shutdown ("ERROR: failure to lock image catalog %s", db.filename);
     112      if (db.dbstate == LCK_EMPTY) Shutdown ("ERROR: No images in catalog %s (1)", db.filename);
     113      if (!dvo_image_load (&db, VERBOSE, FALSE)) Shutdown ("can't read image catalog %s", db.filename);
     114      // the raw FITS data is freed by dvo_image_load, leaving only the Image structure data
     115
     116      /* load regions and images based on specified sky patch (default depth) */
     117      load_images (&db, skylist, FALSE, USE_ALL_IMAGES);
     118
     119      RepairWarps (skylist, HOST_ID, HOSTDIR);
     120     
     121      gfits_db_free (&db);
     122      freeImages();
     123      freeMosaics ();
     124
     125      relastro_client_free (sky, skylist);
     126      break;
     127    }
     128
    102129    case OP_HIGH_SPEED: {
    103130      // USAGE: relastro_client -high-speed
Note: See TracChangeset for help on using the changeset viewer.