IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 39474


Ignore:
Timestamp:
Mar 23, 2016, 12:26:55 PM (10 years ago)
Author:
eugene
Message:

finish code to repair warps

Location:
trunk/Ohana/src/relastro
Files:
1 added
8 edited

Legend:

Unmodified
Added
Removed

  • trunk/Ohana/src/relastro/Makefile

    r39471 r39474  
    2828$(SRC)/RepairWarpMeasures.$(ARCH).o \
    2929$(SRC)/WarpImageMaps.$(ARCH).o \
     30$(SRC)/psps_ids.$(ARCH).o \
    3031$(SRC)/myIndex.$(ARCH).o \
    3132$(SRC)/ConfigInit.$(ARCH).o          \
    … …  
    112113$(SRC)/RepairWarpMeasures.$(ARCH).o \
    113114$(SRC)/WarpImageMaps.$(ARCH).o \
     115$(SRC)/psps_ids.$(ARCH).o \
    114116$(SRC)/myIndex.$(ARCH).o \
    115117$(SRC)/ConfigInit.$(ARCH).o          \

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

    r39471 r39474  
    805805int FindWarpGroups (void);
    806806void FreeWarpGroups (void);
    807 int GetWarpID (Image *image, int obstime, unsigned short photcode, double Rave, double Dave, float X, float Y);
     807int GetWarpSeq (Image *image, int obstime, unsigned short photcode, double Rave, double Dave, float X, float Y);
    808808
    809809myIndexType *myIndexInit ();
    … …  
    813813int myIndexSetEntry (myIndexType *myIndex, int value, int entry);
    814814int myIndexGetEntry (myIndexType *myIndex, int value);
     815
     816uint64_t CreatePSPSObjectID(double ra, double dec);
     817uint64_t CreatePSPSStackDetectionID(int sourceID, int imageID, int detID);
     818uint64_t CreatePSPSDetectionID(double tobs, int ccdid, int detID);

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

    r39457 r39474  
    134134# endif
    135135
     136  free (image);
    136137  free_astrom_table();
    137138}

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

    r39471 r39474  
    99  myAssert (!catalog->Nmeasure || catalog->measure, "programming error");
    1010
     11  int NfixChipID = 0, NfixStackID = 0, NfixWarpID = 0, NfixWarpImageID = 0, NmissWarp = 0, NmissStack = 0, NbadWarp = 0;
     12
     13  int onePercent = catalog->Nmeasure / 100;
     14
    1115  for (off_t j = 0; j < catalog->Nmeasure; j++) {
     16    if (j % onePercent == 0) fprintf (stderr, ".");
     17
    1218    Measure *measure = &catalog->measure[j];
    1319
    14     // we are only going to repair warp detections
    15     if (!isGPC1warp(measure->photcode))  {
    16       // XXX fix extID
     20    // repair extID for non-warps:
     21    if (isGPC1chip(measure->photcode))  {
     22      double mjd = ohana_sec_to_mjd (measure->t);
     23      int ccdnum = measure->photcode % 100;
     24      uint64_t extID = CreatePSPSDetectionID (mjd, ccdnum, measure->detID);
     25      if (extID != measure->extID) {
     26        measure->extID = extID;
     27        NfixChipID ++;
     28      }
    1729      continue;
    1830    }
     31
     32    // repair extID for non-warps:
     33    if (isGPC1stack(measure->photcode))  {
     34      int im = getImageByID (measure->imageID);
     35      if (im < 0) {
     36        if (NmissStack < 10) fprintf (stderr, "missing stack exposure: %f %f : %d %d\n", measure->R, measure->D, measure->imageID, measure->photcode);
     37        NmissStack ++;
     38        continue;
     39      }
     40      uint64_t extID = CreatePSPSStackDetectionID (35, image[im].externID, measure->detID);
     41      if (extID != measure->extID) {
     42        measure->extID = extID;
     43        NfixStackID ++;
     44      }
     45      continue;
     46    }
     47
     48    // we are only going to repair warp detections
     49    if (!isGPC1warp(measure->photcode)) continue;
    1950
    2051    // get the associated average value:
    … …  
    2354
    2455    // double-check for consistency
    25     myAssert (average->objID == measure->objID, "mismatch");
     56    myAssert (average->objID == measure->objID, "objID mismatch: %f %f : %d %d\n", average->R, average->D, average->objID, measure->objID);
     57
     58    // warp coordinates to confirm warp
     59    double X = measure->Xccd;
     60    double Y = measure->Yccd;
    2661
    2762    // check if this detection is far from its correct location
    … …  
    3974
    4075    // skip detections which are within a small distance of the expected location
     76    // NOTE: this actually works surprisingly well near the pole
    4177    double dPos = hypot(dR,dD);
    42     if (dPos < 2.0) continue;
    43 
    44     // warp coordinates to confirm warp
    45     double X = measure->Xccd;
    46     double Y = measure->Yccd;
     78    if (dPos < 2.0) {
     79      int im = getImageByID (measure->imageID);
     80      uint64_t extID = CreatePSPSStackDetectionID (34, image[im].externID, measure->detID);
     81      if (extID != measure->extID) {
     82        measure->extID = extID;
     83        NfixWarpID ++;
     84      }
     85      continue;
     86    }
    4787
    4888    // 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");
     89    int warpSeq = GetWarpSeq (image, measure->t, measure->photcode, Rave, Dave, X, Y);
     90    if (warpSeq < 0) {
     91      if (NmissWarp < 10) fprintf (stderr, "missing warp exposure: %f %f : %d %d\n", Rave, Dave, measure->t, measure->photcode);
     92      NmissWarp ++;
     93      continue;
     94    }
    5195
    5296    // assert on coords being TAN?
    53     Coords *imcoords = &image[warpID].coords;
     97    Coords *imcoords = &image[warpSeq].coords;
    5498
    5599    double Rnew, Dnew;
    … …  
    63107    double dDnew = 3600.0*(Dnew - Dave);
    64108
    65     // skip detections which are within a small distance of the expected location
     109    // detections should now be repaired; fail if not
    66110    dPos = hypot(dRnew,dDnew);
    67     if (dPos > 2.0) {
    68       // complain. abort?
    69       fprintf (stderr, "measurement still far from average location\n");
    70       abort();
     111    if (dPos > 5.0) {
     112      fprintf (stderr, "measurement still far from average location: %f %f vs %f %f : %d %d\n", Rave, Dave, Rnew, Dnew, measure->t, measure->photcode);
     113      NbadWarp ++;
    71114    }
    72115
     116    measure->imageID = image[warpSeq].imageID;
    73117    measure->R = Rnew;
    74118    measure->D = Dnew;
    75119
    76     // XXX fix extID
     120    uint64_t extID = CreatePSPSStackDetectionID (34, image[warpSeq].externID, measure->detID);
     121    if (extID != measure->extID) {
     122      measure->extID = extID;
     123      NfixWarpID ++;
     124    }
     125    NfixWarpImageID ++;
    77126  }
     127
     128  fprintf (stderr, "\n");
     129  fprintf (stderr, "repaired %s : warp image ID %d : ext ID chip %d stack %d warp %d : missed warp %d stack %d : bad warp: %d\n", catalog->filename, NfixWarpImageID, NfixChipID, NfixStackID, NfixWarpID, NmissWarp, NmissStack, NbadWarp);
     130
    78131  return (TRUE);
    79132}

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

    r39471 r39474  
    109109               UserPatch.Rmin, UserPatch.Rmax, UserPatch.Dmin, UserPatch.Dmax);
    110110
    111     if (VERBOSE)       strextend (&command, "-v");
    112     if (VERBOSE2)      strextend (&command, "-vv");
    113     if (UPDATE)        strextend (&command, "-update");
     111    if (VERBOSE)         strextend (&command, "-v");
     112    if (VERBOSE2)        strextend (&command, "-vv");
     113    if (UPDATE)          strextend (&command, "-update");
     114    if (USE_BASIC_CHECK) strextend (&command, "-basic-image-search");
    114115
    115116    fprintf (stderr, "command: %s\n", command);

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

    r39471 r39474  
    1515  double *Dmin;
    1616  double *Dmax;
     17  int    *onBoundary;
    1718} WarpGroup;
    1819
    … …  
    120121    ALLOCATE (warpgroup[i].Dmin, double, warpgroup[i].Nwarps);
    121122    ALLOCATE (warpgroup[i].Dmax, double, warpgroup[i].Nwarps);
     123    ALLOCATE (warpgroup[i].onBoundary, int, warpgroup[i].Nwarps);
    122124
    123125    for (int j = 0; j < warpgroup[i].Nwarps; j++) {
    … …  
    128130      int Ny = image[N].NY;
    129131
     132      warpgroup[i].onBoundary[j] = FALSE;
     133
    130134      double R, D;
    131135      XY_to_RD (&R, &D, 0.0, 0.0, &image[N].coords);
     136      R = ohana_normalize_angle_to_midpoint(R, 180.0);
    132137      warpgroup[i].Rmin[j] = R;
    133138      warpgroup[i].Rmax[j] = R;
    … …  
    137142      // XXX need to worry about 0,360 boundary
    138143      XY_to_RD (&R, &D, Nx, 0.0, &image[N].coords);
     144      R = ohana_normalize_angle_to_midpoint(R, 180.0);
    139145      warpgroup[i].Rmin[j] = MIN (warpgroup[i].Rmin[j], R);
    140146      warpgroup[i].Rmax[j] = MAX (warpgroup[i].Rmax[j], R);
    … …  
    144150      // XXX need to worry about 0,360 boundary
    145151      XY_to_RD (&R, &D, 0.0, Ny, &image[N].coords);
     152      R = ohana_normalize_angle_to_midpoint(R, 180.0);
    146153      warpgroup[i].Rmin[j] = MIN (warpgroup[i].Rmin[j], R);
    147154      warpgroup[i].Rmax[j] = MAX (warpgroup[i].Rmax[j], R);
    … …  
    151158      // XXX need to worry about 0,360 boundary
    152159      XY_to_RD (&R, &D, Nx, Ny, &image[N].coords);
     160      R = ohana_normalize_angle_to_midpoint(R, 180.0);
    153161      warpgroup[i].Rmin[j] = MIN (warpgroup[i].Rmin[j], R);
    154162      warpgroup[i].Rmax[j] = MAX (warpgroup[i].Rmax[j], R);
    155163      warpgroup[i].Dmin[j] = MIN (warpgroup[i].Dmin[j], D);
    156164      warpgroup[i].Dmax[j] = MAX (warpgroup[i].Dmax[j], D);
    157       // bump Rmin,Rmax,Dmin,Dmax by some padding?
     165
     166      // bump Rmin,Rmax,Dmin,Dmax 10 arcsec worth of padding
     167      double dR = 10.0/3600.0 / cos (RAD_DEG*warpgroup[i].Dmin[j]);
     168      double dD = 10.0/3600.0;
     169
     170      // at north pole, force test of all nearby skycells
     171      if (warpgroup[i].Dmax[j] > 89.8) {
     172        warpgroup[i].Rmin[j] =   0.0;
     173        warpgroup[i].Rmax[j] = 360.0;
     174        warpgroup[i].Dmax[j] =  90.0;
     175        warpgroup[i].Dmin[j] -=  dD;
     176        continue;
     177      }
     178
     179      if (warpgroup[i].Rmax[j] - warpgroup[i].Rmin[j] > 270.0) {
     180        // Rmin and Rmax are in the range 0 - 360.  For images at the 0,360 boundary,
     181        // "Rmax" is the lower edge, and "Rmin" is the upper edge.  we need to flip them
     182        // and then break the 0-360 range:
     183        double tmp = warpgroup[i].Rmin[j];
     184        warpgroup[i].Rmin[j] = warpgroup[i].Rmax[j] - 360.0;
     185        warpgroup[i].Rmax[j] = tmp;
     186        warpgroup[i].onBoundary[j] = TRUE;
     187      }
     188
     189      warpgroup[i].Rmin[j] -=  dR;
     190      warpgroup[i].Rmax[j] +=  dR;
     191      warpgroup[i].Dmin[j] -=  dD;
     192      warpgroup[i].Dmax[j] +=  dD;
    158193    }
    159194  }
    … …  
    165200}
    166201
    167 int GetWarpID (Image *image, int obstime, unsigned short photcode, double Rave, double Dave, float X, float Y) {
     202int GetWarpSeq (Image *image, int obstime, unsigned short photcode, double Rave, double Dave, float X, float Y) {
    168203
    169204  int seq = myIndexGetEntry (warpObstimeIndex, obstime);
    … …  
    171206  myAssert (warpgroup[seq].photcode == photcode, "oops");
    172207 
     208  double dPosMin = NAN;
     209  int nPosMin = -1;
     210
    173211  // we now have the warp group, but which is the correct warp?
    174212  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;
     213    if (warpgroup[seq].onBoundary[i]) {
     214      int inRange1 = (Rave >= warpgroup[seq].Rmin[i]) && (Rave <= warpgroup[seq].Rmax[i]);
     215      int inRange2 = (Rave >= warpgroup[seq].Rmin[i] + 360.0) && (Rave <= warpgroup[seq].Rmax[i] + 360.0);
     216      if (!inRange1 && !inRange2) continue;
     217    } else {
     218      if (Rave < warpgroup[seq].Rmin[i]) continue;
     219      if (Rave > warpgroup[seq].Rmax[i]) continue;
     220    }
    177221    if (Dave < warpgroup[seq].Dmin[i]) continue;
    178222    if (Dave > warpgroup[seq].Dmax[i]) continue;
    … …  
    181225    int N = warpgroup[seq].warps[i];
    182226
    183     double Rwrp, Dwrp;
    184     XY_to_RD (&Rwrp, &Dwrp, X, Y, &image[N].coords);
    185 
    186     float csdec = cos(Dave * RAD_DEG);
     227    // project Rave,Dave to image pixels
     228    double Xtst, Ytst;
     229    RD_to_XY (&Xtst, &Ytst, Rave, Dave, &image[N].coords);
     230
     231    // find the pixel offset
     232    double dX = (Xtst - X);
     233    double dY = (Ytst - Y);
    187234   
    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;
     235    // skip detections which are within a small distance from the expected location
     236    double dPos = hypot(dX,dY);
     237    if (dPos < 20.0) return N; // 20.0 pixels = 5.0 arcsec
    195238
    196239    // check for multiple possible matches??
    197   }
    198   return -1;
     240
     241    // record the min pos so we can see how things fail
     242    if (isnan(dPosMin)) {
     243      dPosMin = dPos;
     244      nPosMin = i;
     245    } else {
     246      if (dPos < dPosMin) {
     247        dPosMin = dPos;
     248        nPosMin = i;
     249      }
     250    }
     251  }
     252 
     253  // if we did not find any matched, give up
     254  if (nPosMin < 0) return -1;
     255
     256  // return closest image
     257  return warpgroup[seq].warps[nPosMin];
    199258}
    200259
    … …  
    207266    free (warpgroup[i].Dmin);
    208267    free (warpgroup[i].Dmax);
     268    free (warpgroup[i].onBoundary);
    209269  }   
    210270

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

    r39457 r39474  
    4444   
    4545  /* unlock, if we can (else, unlocked below) and free, if we can */
    46   int unlockImages = !UPDATE || (RELASTRO_OP == OP_UPDATE_OFFSETS);
     46  int unlockImages = !UPDATE || (RELASTRO_OP == OP_UPDATE_OFFSETS) || (RELASTRO_OP == OP_REPAIR_WARPS);
    4747  if (unlockImages) {
    4848    if (subset != image) {

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

    r37807 r39474  
    110110  for (i = 0; i < Ntimage; i++) {
    111111     
    112     if (FALSE && !strncmp(timage[i].name, "o5479g0238o", 10)) {
     112    if (FALSE && !strncmp(timage[i].name, "o6406g0242o", 10)) {
    113113      fprintf (stderr, "test image 1\n");
    114114    }
    115115    if (FALSE && !strncmp(timage[i].name, "o6227g0311o", 10)) {
    116116      fprintf (stderr, "test image 2\n");
     117    }
     118
     119    if (FALSE && (i >= 42819857) && (i < 42819972)) {
     120      fprintf (stderr, "test image %s\n", timage[i].name);
    117121    }
    118122
    … …  
    173177      DmaxImage = MAX(DmaxImage, Di[j]);
    174178    }
    175     if (RmaxImage - RminImage > 180.0) {
     179    if (FALSE && (RmaxImage - RminImage > 180.0)) {
    176180        double tmp = RminImage;
    177181        RmaxImage = RminImage;
Note: See TracChangeset for help on using the changeset viewer.